Quantile-based estimators with Qest

library(Qest)
library(survival)

1. The Q-estimator

A parametric model can be described by its quantile function \(Q(\tau \mid x, \theta)\). The Q-estimator of Sottile and Frumento (2022) minimizes the integrated loss of quantile regression, \[ L(\theta) = \sum_{i=1}^n \int_0^1 w(\tau)\,\rho_\tau\{y_i - Q(\tau \mid x_i, \theta)\}\,d\tau, \qquad \rho_\tau(u) = u\{\tau - I(u < 0)\}, \] where \(w(\tau) \ge 0\) is a weighting function. The estimator is consistent and asymptotically normal for any correctly specified model, and it is usually more robust than maximum likelihood; truncated weights such as \(w(\tau) = I(\delta \le \tau \le 1 - \delta)\) (function wtrunc()) increase its robustness further. Censored and truncated data are handled by replacing the check function with its conditional expectation.

2. User-defined quantile functions

Q(theta, tau, data) must return \(Q(\tau \mid x_i, \theta)\) for each row of data, where tau is either a matrix (one row for each observation) or a vector. Here is a normal linear model.

Qlinmod <- function(theta, tau, data) {
  X <- model.matrix(~ x1 + x2, data = data)
  qnorm(tau, c(X %*% theta[-1]), exp(theta[1]))
}
set.seed(1)
n <- 200
dat <- data.frame(x1 = rnorm(n), x2 = runif(n, 0, 3))
dat$y <- 4 + dat$x1 + 2 * dat$x2 + rnorm(n)

m1 <- Qest(y ~ x1 + x2, Q = Qlinmod, start = rep(NA, 4), data = dat)
summary(m1)
#> 
#> Call: Qest(formula = y ~ x1 + x2, Q = Qlinmod, start = rep(NA, 4), 
#>     data = dat)
#> 
#> n. of observations: 200 
#> n. of parameters: 4 
#> n. of iterations: 4 
#> 
#> Coefficients:
#>        Estimate Std. Error z value Pr(>|z|)    
#> Theta1  0.06302    0.05689   1.108    0.268    
#> Theta2  3.90198    0.16266  23.988   <2e-16 ***
#> Theta3  1.08346    0.08052  13.456   <2e-16 ***
#> Theta4  2.06345    0.09412  21.923   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Minimized loss function: 60.43505

The length of start defines the number of parameters; NA values are replaced by data-driven starting points (Algorithm 1 of the paper). The derivatives of Q with respect to \(\theta\) are computed by central differences, and all integrals in \(\tau\) on a grid of ntau points.

The same model is available as a family (Qnorm) and, for complete data, as Qlm(), both with closed-form estimating equations. The user-defined version agrees with them:

m2 <- Qest(y ~ x1 + x2, Q = Qnorm, data = dat)
m3 <- Qlm(y ~ x1 + x2, data = dat)
cbind(user = coef(m1), Qnorm = coef(m2), Qlm = c(log(m3$sigma), coef(m3)))
#>              user      Qnorm        Qlm
#> Theta1 0.06301937 0.06312662 0.06312662
#> Theta2 3.90197624 3.90196715 3.90196715
#> Theta3 1.08346264 1.08344991 1.08344991
#> Theta4 2.06345128 2.06346288 2.06346288

Censored data

With a Surv() response, right-censored and left-truncated data are allowed. A Weibull proportional hazards model:

QWeibPH <- function(theta, tau, data) {
  shape <- exp(theta[1])
  X <- model.matrix(~ x, data = data)
  qweibull(tau, shape = shape, scale = exp(-c(X %*% theta[-1]) / shape))
}
dat <- data.frame(x = runif(n, 0, 3))
t <- QWeibPH(c(log(2), -0.5, 1), runif(n), dat)
cens <- runif(n, 0, 2)
dat$y <- pmin(t, cens)
dat$d <- as.numeric(t <= cens)
m4 <- Qest(Surv(y, d) ~ x, Q = QWeibPH, start = rep(NA, 3), data = dat)
summary(m4)
#> 
#> Call: Qest(formula = Surv(y, d) ~ x, Q = QWeibPH, start = rep(NA, 3), 
#>     data = dat)
#> 
#> n. of observations: 200 
#> n. of parameters: 3 
#> n. of iterations: 2 
#> 
#> Coefficients:
#>        Estimate Std. Error z value Pr(>|z|)    
#> Theta1  0.65603    0.07749   8.467  < 2e-16 ***
#> Theta2 -0.49349    0.17256  -2.860  0.00424 ** 
#> Theta3  1.03640    0.12404   8.356  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Loss (not the function being minimized): 11.18698

If the model is correct, \(F(y_i \mid x_i)\) is uniformly distributed: plot() compares its (Kaplan-Meier) distribution with the uniform one.

plot(m4)

3. Families

Qnorm, Qgamma, Qpois (a continuous version of the Poisson distribution, to be used with jittered counts \(y + U\)) and Qunif have analytical estimating equations and are much faster.

x <- runif(n)
yg <- rgamma(n, shape = 2, scale = exp(1 + x))
coef(Qest(yg ~ x, Q = Qgamma))
#>  log(shape) (Intercept)           x 
#>   0.5144618   1.2627170   0.8061452
yp <- rpois(n, exp(1 + x)) + runif(n)
coef(Qest(yp ~ x, Q = Qpois))
#> (Intercept)           x 
#>   0.9540547   1.0419677
yu <- runif(n, 0, 1 + 2 * x)
coef(Qest(yu ~ x, Q = Qunif(min = FALSE)))
#> max: (Intercept)           max: x 
#>        0.9470187        2.0461569

predict() returns fitted quantiles or the fitted distribution function:

fg <- Qest(yg ~ x, Q = Qgamma)
predict(fg, newdata = data.frame(x = c(0, 1)), p = c(0.1, 0.5, 0.9))
#>       p0.1      p0.5     p0.9
#> 1 1.308812  4.785588 12.00040
#> 2 2.930770 10.716174 26.87201

4. Robustness

We contaminate 5% of a normal sample and compare maximum likelihood (least squares) with Q-estimators with and without truncated weights.

set.seed(2)
x <- runif(n)
y <- 1 + 2 * x + rnorm(n)
out <- sample(n, 10)
y[out] <- y[out] + 15
rbind(lm = coef(lm(y ~ x)),
      Qlm = coef(Qlm(y ~ x)),
      Qlm.trunc = coef(Qlm(y ~ x, wtau = wtrunc, delta.left = 0.1, delta.right = 0.1)))
#>           (Intercept)        x
#> lm           1.604800 2.475819
#> Qlm          1.162206 1.999379
#> Qlm.trunc    1.157532 1.968077

5. Proportional hazards models

Qcoxph() fits a proportional hazards model with a piecewise linear baseline cumulative hazard, using closed-form estimating equations.

set.seed(3)
dat <- data.frame(x = runif(n, 0, 3))
t <- rweibull(n, shape = 2, scale = (1 / exp(2 + 2 * dat$x))^(1 / 2))
cens <- runif(n, 0, 1)
dat$y <- pmin(t, cens); dat$d <- as.numeric(t <= cens)
fc <- Qcoxph(Surv(y, d) ~ x, data = dat)
summary(fc)
#> Call:
#> Qcoxph(formula = Surv(y, d) ~ x, data = dat)
#> 
#>   n = 200, number of events = 185
#> 
#>     coef exp(coef) se(coef)     z Pr(>|z|)    
#> x 1.7122    5.5413   0.1503 11.39   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>   exp(coef) exp(-coef) lower .95 upper .95
#> x     5.541     0.1805     4.127      7.44
#> 
#> Concordance = 0.799 (se = 0.015)
#> Objective function = 2.59 (not the function being minimized)
cbind(coxph = coef(coxph(Surv(y, d) ~ x, data = dat)), Qcoxph = coef(fc))
#>     coxph   Qcoxph
#> x 1.92666 1.712227
sf <- survfit(fc, newdata = data.frame(x = c(0.5, 1.5)))
plot(sf, col = 1:2, conf.int = TRUE, xlab = "time", ylab = "survival")

6. Numerical details

  • The grid of quantile orders has ntau = 199 points (denser in the tails) in \([5 \cdot 10^{-4}, 1 - 5 \cdot 10^{-4}]\); the tails are integrated with rectangles.
  • The estimating equations are solved by Newton-Raphson with backtracking, with Levenberg-Marquardt steps when needed; the Jacobian is the exact derivative of the discretized equations. Standard errors use the sandwich formula \(J^{-1}\Omega J^{-\top}\), where \(J\) is not symmetric with censored or truncated data.
  • Qest.control() sets the tolerance, the maximum number of iterations and the step of the numerical derivatives.

References

Sottile, G. and Frumento, P. (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis, 171, 107471. doi:10.1016/j.csda.2022.107471