Penalized and nonlinear quantile regression coefficients modeling with qrcmNP

library(qrcmNP)

1. Quantile regression coefficients modeling

In quantile regression the conditional quantile function of a response \(y\) given covariates \(x\) is \(Q(p \mid x) = x^\top\beta(p)\), \(0 < p < 1\). Standard quantile regression estimates one quantile at a time; quantile regression coefficients modeling (QRCM, Frumento and Bottai, 2016) instead describes the coefficients as parametric functions of \(p\), \[\beta(p \mid \theta) = \theta\, b(p),\] where \(b(p) = (b_1(p), \dots, b_k(p))^\top\) are known functions and \(\theta\) is a \(q \times k\) matrix of parameters. All quantiles are estimated simultaneously by minimizing the integrated loss \[L(\theta) = \int_0^1 \sum_{i=1}^n \rho_p\{y_i - x_i^\top\theta\, b(p)\}\,dp,\] which is a smooth convex function of \(\theta\) whose gradient and Hessian are available in closed form. The package ‘qrcm’ implements this estimator (qrcm::iqr()); ‘qrcmNP’ adds penalized estimation for covariate selection (piqr()) and nonlinear parametric models (niqr()).

2. Penalized QRCM: `piqr()`

Sottile et al. (2020) proposed to select covariates by minimizing \[L(\theta) + \lambda\, \mathrm{pen}(\theta)\] with standardized \(x\) and \(b(p)\). Since \(\beta_j(p \mid \theta) \equiv 0\) if and only if \(\theta_{j1} = \dots = \theta_{jk} = 0\), selection uses the information on all quantiles at once. Two penalties are available:

penalty = "lasso" (default): \(\lambda\sum_{j,h} w_j|\theta_{jh}|\), the penalty of Sottile et al. (2020, eq. 3.4);
penalty = "group": \(\lambda\sum_j w_j\sqrt{k_j}\,\|\theta_{j\cdot}\|_2\), a group lasso that keeps or removes each covariate as a whole.

The weights \(w_j\) are set by penalty.factor (default 1; 0 leaves a covariate unpenalized). The intercept is never penalized.

A simulated example

We simulate data from the model \(\beta_0(p) = 1 + \Phi^{-1}(p)\), \(\beta_1(p) = 1 + 2p\), \(\beta_2(p) = 1\), with three irrelevant covariates.

set.seed(1234)
n <- 400
d <- data.frame(x1 = rexp(n), x2 = runif(n, 0, 5), x3 = rnorm(n),
                x4 = runif(n), x5 = rbinom(n, 1, 0.5))
u <- runif(n)
d$y <- 1 + qnorm(u) + (1 + 2 * u) * d$x1 + d$x2

The coefficient functions are modelled by shifted Legendre polynomials of degree 2 (slp(p, 2), plus a constant). piqr() computes the whole path of solutions, starting from \(\lambda_{max}\) (only the intercept in the model):

fit <- piqr(y ~ x1 + x2 + x3 + x4 + x5, formula.p = ~ slp(p, 2), data = d,
            nlambda = 40)
fit
#> 
#> Call:
#> piqr(formula = y ~ x1 + x2 + x3 + x4 + x5, formula.p = ~slp(p, 
#>     2), data = d, nlambda = 40)
#> 
#> Penalty: lasso   n. of lambda values: 40
#> 
#>           Df  Npar  ObjFunction     Lambda
#> lambda1    1     3        334.6  80.490000
#> lambda2    3     6        259.0  63.560000
#> lambda3    3     6        225.7  50.190000
#> lambda4    3     7        210.7  39.630000
#> lambda5    3     7        196.1  31.300000
#> lambda6    3     7        188.5  24.710000
#> lambda7    3     8        183.6  19.520000
#> lambda8    3     8        180.1  15.410000
#> lambda9    4    10        177.8  12.170000
#> lambda10   4    10        176.3   9.609000
#> lambda11   4    10        175.4   7.588000
#> lambda12   4     9        174.9   5.992000
#> lambda13   5    10        174.5   4.731000
#> lambda14   5    11        174.3   3.736000
#> lambda15   5    11        174.1   2.950000
#> lambda16   5    11        174.1   2.330000
#> lambda17   5    11        174.0   1.840000
#> lambda18   5    11        174.0   1.453000
#> lambda19   5    11        173.9   1.147000
#> lambda20   6    13        173.9   0.905800
#> lambda21   6    13        173.8   0.715300
#> lambda22   6    15        173.8   0.564800
#> lambda23   6    14        173.8   0.446000
#> lambda24   6    15        173.7   0.352200
#> lambda25   6    17        173.7   0.278100
#> lambda26   6    17        173.7   0.219600
#> lambda27   6    17        173.7   0.173400
#> lambda28   6    17        173.7   0.136900
#> lambda29   6    17        173.7   0.108100
#> lambda30   6    17        173.7   0.085390
#> lambda31   6    17        173.7   0.067430
#> lambda32   6    17        173.7   0.053240
#> lambda33   6    17        173.6   0.042040
#> lambda34   6    17        173.6   0.033200
#> lambda35   6    17        173.6   0.026220
#> lambda36   6    18        173.6   0.020700
#> lambda37   6    18        173.6   0.016350
#> lambda38   6    18        173.6   0.012910
#> lambda39   6    18        173.6   0.010190
#> lambda40   6    18        173.6   0.008049

oldpar <- par(mfrow = c(1, 2))
plot(fit, xvar = "lambda")
plot(fit, xvar = "grad")

par(oldpar)

Selecting \(\lambda\)

gof.piqr() computes the AIC and BIC of Sottile et al. (2020, Section 3.1), \(\log L(\hat\theta_\lambda) + c_n\,df(\lambda)/n\), with \(c_n = 2\) (AIC) or \(c_n = \log n\) (BIC):

best <- gof.piqr(fit, method = "BIC")

best
#> 
#> Selection of lambda by BIC 
#> 
#> Selected lambda: 5.992  (position 12 )
#> Degrees of freedom: 9 
#> 
#> Coefficients:
#>              (Intercept)  slp1      slp2    
#> (Intercept)  -0.36617      1.65368   0.10827
#> x1            1.01840      0.83515  -0.13759
#> x2            0.94377           NA        NA
#> x3                 NA           NA        NA
#> x4                 NA           NA        NA
#> x5                 NA      0.09813  -0.10019

Cross-validation of the integrated loss is also available:

cv <- cv_piqr(y ~ x1 + x2 + x3 + x4 + x5, formula.p = ~ slp(p, 2), data = d,
              nfolds = 5, nlambda = 40)
cv
#> 
#> Cross-validated integrated loss (5 folds)
#> 
#>     Lambda Index Measure      SE Df
#> min  4.731    13  0.4451 0.01506  5
#> 1se 15.410     8  0.4557 0.01121  3
plot(cv)

Group lasso

fitg <- piqr(y ~ x1 + x2 + x3 + x4 + x5, formula.p = ~ slp(p, 2), data = d,
             nlambda = 40, penalty = "group")
bestg <- gof.piqr(fitg, plot = FALSE)
bestg
#> 
#> Selection of lambda by BIC 
#> 
#> Selected lambda: 14.37  (position 8 )
#> Degrees of freedom: 9 
#> 
#> Coefficients:
#>              (Intercept)  slp1     slp2   
#> (Intercept)   1.1980       1.5147   0.9307
#> x1            0.5802       0.7119  -0.4972
#> x2            0.3155       0.1943  -0.2793
#> x3                NA           NA       NA
#> x4                NA           NA       NA
#> x5                NA           NA       NA

Inference after selection

As suggested by Sottile et al. (2020, Section 3.2), the penalized estimator is used to select the covariates, and the selected model is re-estimated without penalty. refit_piqr() does this with qrcm::iqr(), which provides standard errors, Wald tests and goodness-of-fit diagnostics (the standard errors do not account for the selection step):

m <- refit_piqr(fit, pos.lambda = best$posMinLambda)
summary(m)
#> 
#> Call: qrcm::iqr(formula = y ~ x1 + x2 + x3 + x4 + x5, formula.p = ~slp(p, 
#>     2), data = d, s = structure(c(1, 1, 1, 0, 0, 1, 1, 1, 1, 
#> 0, 0, 1, 1, 1, 1, 0, 0, 1), dim = c(6L, 3L)))
#> 
#> converged: TRUE 
#> n. of iterations: 10 
#> n. of free parameters: 12 
#> n. of observations: 400 
#> 
#> ###################### 
#> ###################### 
#> 
#> Coefficients:
#>              (Intercept)  slp1      slp2    
#> (Intercept)  -0.99415      1.84750  -0.01593
#> x1            1.21977      0.80738  -0.08405
#> x2            1.12352     -0.08471   0.03807
#> x3                 NA           NA        NA
#> x4                 NA           NA        NA
#> x5            0.08917      0.12613  -0.12671
#> 
#> Standard errors:
#>              (Intercept)  slp1     slp2   
#> (Intercept)  0.38322      0.21992  0.20119
#> x1           0.17633      0.09931  0.14888
#> x2           0.10656      0.06213  0.05731
#> x3                NA           NA       NA
#> x4                NA           NA       NA
#> x5           0.33200      0.19311  0.16572
#> 
#> ###################### 
#> ###################### 
#> 
#> Wald test for x:
#>             chi-square df P(> chi)    
#> (Intercept)    126.241  3   <2e-16 ***
#> x1             890.682  3   <2e-16 ***
#> x2             432.139  3   <2e-16 ***
#> x3                  NA  0       NA    
#> x4                  NA  0       NA    
#> x5               6.509  3   0.0893 .  
#> 
#> 
#> Wald test for b(p):
#>             chi-square df P(> chi)    
#> (Intercept)    616.893  4   <2e-16 ***
#> slp1           784.854  4   <2e-16 ***
#> slp2             1.582  4    0.812    
#> 
#> Minimized loss function: 173.9093
summary(m, p = c(0.1, 0.5, 0.9))
#> 
#> Call: qrcm::iqr(formula = y ~ x1 + x2 + x3 + x4 + x5, formula.p = ~slp(p, 
#>     2), data = d, s = structure(c(1, 1, 1, 0, 0, 1, 1, 1, 1, 
#> 0, 0, 1, 1, 1, 1, 0, 0, 1), dim = c(6L, 3L)))
#> 
#> p = 0.1 
#> 
#> Coefficients:
#>             Estimate  std.err z value p(>|z|))    
#> (Intercept) -0.61605  0.28968  -2.127   0.0334 *  
#> x1           1.42664  0.12511  11.403   <2e-16 ***
#> x2           1.08602  0.08075  13.449   <2e-16 ***
#> x3                NA       NA      NA       NA    
#> x4                NA       NA      NA       NA    
#> x5           0.18281  0.25353   0.721   0.4709    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> 
#> 
#> p = 0.5 
#> 
#> Coefficients:
#>             Estimate std.err z value p(>|z|))    
#> (Intercept)  0.87724 0.22157   3.959 7.52e-05 ***
#> x1           2.15323 0.16310  13.202  < 2e-16 ***
#> x2           0.98170 0.05992  16.385  < 2e-16 ***
#> x3                NA      NA      NA       NA    
#> x4                NA      NA      NA       NA    
#> x5           0.40536 0.17269   2.347   0.0189 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> 
#> 
#> p = 0.9 
#> 
#> Coefficients:
#>             Estimate std.err z value p(>|z|))    
#> (Intercept)  2.33995 0.21604  10.831   <2e-16 ***
#> x1           2.71845 0.09903  27.451   <2e-16 ***
#> x2           0.95048 0.06296  15.097   <2e-16 ***
#> x3                NA      NA      NA       NA    
#> x4                NA      NA      NA       NA    
#> x5           0.38462 0.19612   1.961   0.0499 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

oldpar <- par(mfrow = c(1, 3))
plot(m, which = 1:3, ask = FALSE)

par(oldpar)

Censored data

Right-censored and left-truncated responses are specified with survival::Surv(), as in qrcm::iqr():

cens <- rexp(n, 0.05)
d$time <- pmin(d$y + 5, cens)
d$status <- as.numeric(d$y + 5 <= cens)
fitc <- piqr(survival::Surv(time, status) ~ x1 + x2 + x3 + x4 + x5,
             formula.p = ~ slp(p, 2), data = d, nlambda = 20)
gof.piqr(fitc, plot = FALSE)
#> 
#> Selection of lambda by BIC 
#> 
#> Selected lambda: 2.723  (position 7 )
#> Degrees of freedom: 8 
#> 
#> Coefficients:
#>              (Intercept)  slp1       slp2     
#> (Intercept)   4.223977     1.798755   0.079512
#> x1            1.069666     0.848069         NA
#> x2            1.016954           NA  -0.002393
#> x3                  NA           NA         NA
#> x4                  NA           NA         NA
#> x5                  NA           NA  -0.189861

Algorithm

For each value of \(\lambda\), piqr() uses a proximal Newton method: the integrated loss is replaced by its second-order expansion (its gradient and Hessian are available in closed form), and the resulting penalized quadratic problem is solved by coordinate descent (soft-thresholding for the lasso, group soft-thresholding for the group lasso) followed by an exact “polishing” step on the selected coefficients; a backtracking line search guarantees a decrease of the penalized objective. Along the path the algorithm uses warm starts, the sequential strong rules of Tibshirani et al. (2012) and a final check of the Karush-Kuhn-Tucker conditions. The conditional CDF of each observation is found by bisection on a grid of 1023 values of \(p\), denser in the tails, on which \(b(p)\) and its integrals are pre-computed (as in ‘qrcm’). The algorithm of ‘qrcmNP’ 0.2.x is available with method = "legacy".

3. Nonlinear parametric quantile models: `niqr()`

Bottai and Cilluffo (2020) considered parametric models \(Q(p \mid x, \theta)\) of any functional form, estimated by minimizing \[L_n(\theta) = \frac1n\sum_{i=1}^n\int_0^1 \rho_p\{y_i - Q(p \mid x_i, \theta)\}\,dp.\] In niqr(), the model is specified by two functions:

fun(theta, p) returns the matrix of coefficient functions \(\beta(p \mid \theta)\) (one row for each value of p);
fun2(beta, X) (optional) combines them with the covariates, row by row; by default \(Q(p \mid x) = x^\top\beta(p \mid \theta)\).

Nonlinear in the parameters

Following Bottai and Cilluffo (2020, Section 6), the intercept can be modelled as a transformation of a normal quantile function, \(z = \Phi^{-1}(p)\), with polynomial effects of the covariates. Here the data have heavier tails than the normal distribution:

set.seed(1)
n <- 400
x <- runif(n)
z <- qnorm(runif(n))
y <- 2 + 0.5 * z + 0.1 * z^3 + (1 + pnorm(z)) * x
dn <- data.frame(y, x)

# Model 1: normal intercept
m1 <- niqr(y ~ x, data = dn, x0 = c(2, 0, 1, 0),
           fun = function(theta, p) cbind(theta[1] + exp(theta[2]) * qnorm(p),
                                          theta[3] + theta[4] * p))
# Model 2: a cubic term accounts for the kurtosis (exp() keeps it positive,
# so that the quantile function is increasing in p)
m2 <- niqr(y ~ x, data = dn, x0 = c(2, 0, -2, 1, 0),
           fun = function(theta, p) {
             z <- qnorm(p)
             cbind(theta[1] + exp(theta[2]) * z + exp(theta[3]) * z^3, theta[4] + theta[5] * p)
           })
gof.niqr(m1, m2)
#>    npar      loss      AIC      BIC   AIC.log   BIC.log
#> m1    4 0.2537794 109.5117 125.4776 -1.361290 -1.341333
#> m2    5 0.2528031 111.1212 131.0785 -1.362644 -1.337698
summary(m2)
#> 
#> Call: niqr(formula = y ~ x, data = dn, fun = function(theta, p) {
#>     z <- qnorm(p)
#>     cbind(theta[1] + exp(theta[2]) * z + exp(theta[3]) * z^3, 
#>         theta[4] + theta[5] * p)
#> }, x0 = c(2, 0, -2, 1, 0))
#> 
#> n. of observations: 400    n. of iterations: 4 
#> Standard errors: sandwich 
#> 
#>        Estimate Std. Error z value Pr(>|z|)    
#> theta0  2.03442    0.07919  25.690  < 2e-16 ***
#> theta1 -0.59066    0.16849  -3.506 0.000456 ***
#> theta2 -2.22999    0.21187 -10.525  < 2e-16 ***
#> theta3  1.23175    0.26375   4.670 3.01e-06 ***
#> theta4  0.55099    0.42081   1.309 0.190417    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Minimized loss function: 0.2528031 
#> AIC = 111.121   BIC = 131.079   AIC* = -1.36264   BIC* = -1.3377

The log-criteria \(\mathrm{AIC}^* = \log L_n + k/n\) and \(\mathrm{BIC}^* = \log L_n + k\log(n)/(2n)\) (Bottai and Cilluffo, 2020, Section 3) are reported as AIC.log and BIC.log; here AIC* prefers the second model, while BIC* slightly prefers the simpler one. The goodness-of-fit test below is more conclusive. The estimated coefficient functions are plotted with pointwise 95% confidence bands:

oldpar <- par(mfrow = c(1, 3))
plot(m2, which = 1:3, ask = FALSE)

par(oldpar)

Under a correctly specified model the fitted values of the conditional CDF are uniformly distributed; test.fit() computes the Kolmogorov-Smirnov and Cramer-von Mises statistics, with bootstrap p-values. The normal model is rejected, the model with heavy tails is not:

test.fit(m1, R = 30)
#>                     statistic p-value
#> Kolmogorov-Smirnov 0.05454231       0
#> Cramer-Von Mises   0.20379778       0
test.fit(m2, R = 30)
#> Warning: 1 bootstrap replications failed
#>                     statistic   p-value
#> Kolmogorov-Smirnov 0.02994790 0.3103448
#> Cramer-Von Mises   0.03403249 0.5862069

Nonlinear in the covariates: a dose-response model

The Hill quantile function (Bottai and Cilluffo, 2020, eq. 2) \[Q(p \mid \text{dose}) = \beta_0(p) + \frac{\beta_1(p) - \beta_0(p)}{1 + 10^{\beta_3(p) - \log_{10}\text{dose}}}\] is nonlinear in both the covariate and the parameters. Here the minimal response \(\beta_0(p)\) is normal, and the maximal response and the (log) dose at half the maximal effect are constant:

set.seed(2)
dose <- rep(10^seq(-2, 2, length.out = 20), each = 15)
fun.h <- function(theta, p) cbind(theta[1] + exp(theta[2]) * qnorm(p), theta[3], theta[4])
fun2.h <- function(beta, X) beta[, 1] + (beta[, 2] - beta[, 1]) / (1 + 10^(beta[, 3] - log10(X[, 2])))
yh <- fun2.h(fun.h(c(1, log(0.5), 10, 0), runif(length(dose))), cbind(1, dose))
mh <- niqr(yh ~ dose, fun = fun.h, fun2 = fun2.h, x0 = c(0, 0, 5, 0.5))
summary(mh)
#> 
#> Call: niqr(formula = yh ~ dose, fun = fun.h, fun2 = fun2.h, x0 = c(0, 
#>     0, 5, 0.5))
#> 
#> n. of observations: 300    n. of iterations: 16 
#> Standard errors: sandwich 
#> 
#>         Estimate Std. Error  z value Pr(>|z|)    
#> theta0  1.003353   0.056302   17.821   <2e-16 ***
#> theta1 -0.632827   0.054445  -11.623   <2e-16 ***
#> theta2 10.001469   0.002756 3629.444   <2e-16 ***
#> theta3  0.003041   0.006897    0.441    0.659    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Minimized loss function: 0.073441 
#> AIC = 30.0323   BIC = 44.8474   AIC* = -2.59794   BIC* = -2.57325

dd <- data.frame(dose = 10^seq(-2, 2, length.out = 100))
qf <- predict(mh, type = "QF", p = c(0.05, 0.5, 0.95), newdata = dd, se = FALSE)
plot(dose, yh, log = "x", pch = 20, col = "gray", xlab = "dose", ylab = "response")
matlines(dd$dose, qf, lty = c(2, 1, 2), col = 1)

Inference

The covariance matrix of \(\hat\theta\) is estimated by the sandwich formula \(H^{-1}\Omega H^{-1}/n\) (Bottai and Cilluffo, 2020, Theorem 3), where \(H\) is the Hessian of the expected loss, \(H = n^{-1}\sum_i \partial_\theta Q(p_i^*)\,\partial_\theta Q(p_i^*)^\top/Q'(p_i^*)\), \(p_i^* = F(y_i \mid x_i, \hat\theta)\). With repeated measures, the argument cluster gives cluster-robust standard errors; the (cluster) bootstrap used in the paper is available with control = niqr.control(se = "boot").

References

Bottai, M. and Cilluffo, G. (2020). Nonlinear parametric quantile models. Statistical Methods in Medical Research, 29(12), 3757-3769.
Frumento, P. and Bottai, M. (2016). Parametric modeling of quantile regression coefficient functions. Biometrics, 72(1), 74-84.
Frumento, P. and Bottai, M. (2017). Parametric modeling of quantile regression coefficient functions with censored and truncated data. Biometrics, 73(4), 1179-1188.
Sottile, G., Frumento, P., Chiodi, M. and Bottai, M. (2020). A penalized approach to covariate selection through quantile regression coefficient models. Statistical Modelling, 20(4), 369-385.
Tibshirani, R., Bien, J., Friedman, J., Hastie, T., Simon, N., Taylor, J. and Tibshirani, R. J. (2012). Strong rules for discarding predictors in lasso-type problems. Journal of the Royal Statistical Society: Series B, 74(2), 245-266.