---
title: "Quantile-based estimators with Qest"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quantile-based estimators with Qest}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4)
```

```{r setup}
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.

```{r user}
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)
```

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:

```{r compare}
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)))
```

### Censored data

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

```{r weibull}
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)
```

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

```{r plot}
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.

```{r families}
x <- runif(n)
yg <- rgamma(n, shape = 2, scale = exp(1 + x))
coef(Qest(yg ~ x, Q = Qgamma))
yp <- rpois(n, exp(1 + x)) + runif(n)
coef(Qest(yp ~ x, Q = Qpois))
yu <- runif(n, 0, 1 + 2 * x)
coef(Qest(yu ~ x, Q = Qunif(min = FALSE)))
```

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

```{r predict}
fg <- Qest(yg ~ x, Q = Qgamma)
predict(fg, newdata = data.frame(x = c(0, 1)), p = c(0.1, 0.5, 0.9))
```

## 4. Robustness

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

```{r robust}
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)))
```

## 5. Proportional hazards models

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

```{r cox}
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)
cbind(coxph = coef(coxph(Surv(y, d) ~ x, data = dat)), Qcoxph = coef(fc))
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