qgam: quantile non-parametric additive models

This R package offers methods for fitting additive quantile regression models based on splines, using the methods described in Fasiolo et al., 2017.

The main fitting functions are:

• qgam() fits an additive quantile regression model to a single quantile. Very similar to mgcv::gam(). It returns an object of class qgam, which inherits from mgcv::gamObject.
• mqgam() fits the same additive quantile regression model to several quantiles. It is more efficient that calling qgam() several times, especially in terms of memory usage.
• tuneLearn() useful for tuning the learning rate of the Gibbs posterior. It evaluates a calibration loss function on a grid of values provided by the user.
• tuneLearnFast() similar to tuneLearn(), but here the learning rate is selected by minimizing the calibration loss, using Brent method.

1 A first example: smoothing the motorcycle dataset

Let's start with a simple example. Here we are fitting a regression model with an adaptive spline basis to quantile 0.8 of the motorcycle dataset.

library(qgam); library(MASS)
if( suppressWarnings(require(RhpcBLASctl)) ){ blas_set_num_threads(1) } # Optional

fit <- qgam(accel~s(times, k=20, bs="ad"), 
            data = mcycle, 
            qu = 0.8)

## Estimating learning rate. Each dot corresponds to a loss evaluation. 
## qu = 0.8............done

# Plot the fit
xSeq <- data.frame(cbind("accel" = rep(0, 1e3), "times" = seq(2, 58, length.out = 1e3)))
pred <- predict(fit, newdata = xSeq, se=TRUE)
plot(mcycle$times, mcycle$accel, xlab = "Times", ylab = "Acceleration", ylim = c(-150, 80))
lines(xSeq$times, pred$fit, lwd = 1)
lines(xSeq$times, pred$fit + 2*pred$se.fit, lwd = 1, col = 2)
lines(xSeq$times, pred$fit - 2*pred$se.fit, lwd = 1, col = 2)   

qgam automatically calls tuneLearnFast to select the learning rate. The results of the calibrations are stored in fit$calibr. We can check whether the optimization succeded as follows:

check(fit$calibr, 2)

The plot suggest that the calibration criterion has a single minimum, and that the optimizer has converged to its neighbourhood. Alternatively, we could have selected the learning rate by evaluating the loss function on a grid.

set.seed(6436)
cal <- tuneLearn(accel~s(times, k=20, bs="ad"), 
                 data = mcycle, 
                 qu = 0.8,
                 lsig = seq(1, 3, length.out = 20), 
                 control = list("progress" = "none")) #<- sequence of values for learning rate
                 
check(cal)