Penalized and nonlinear quantile regression coefficients modeling with qrcmNP

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.

Penalized and nonlinear quantile regression coefficients modeling with qrcmNP

1. Quantile regression coefficients modeling

2. Penalized QRCM: piqr()

A simulated example

Selecting \(\lambda\)

Group lasso

Inference after selection

Censored data

Algorithm

3. Nonlinear parametric quantile models: niqr()

Nonlinear in the parameters

Nonlinear in the covariates: a dose-response model

Inference

References

2. Penalized QRCM: `piqr()`

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