Fit multiple smooth additive quantile regression models

This function fits a smooth additive regression model to several quantiles.

Usage

mqgam(
  form,
  data,
  qu,
  discrete = FALSE,
  lsig = NULL,
  err = NULL,
  multicore = !is.null(cluster),
  cluster = NULL,
  ncores = detectCores() - 1,
  paropts = list(),
  control = list(),
  argGam = NULL
)

Arguments

form: A GAM formula, or a list of formulae. See ?mgcv::gam details.
data: A data frame or list containing the model response variable and covariates required by the formula. By default the variables are taken from environment(formula): typically the environment from which gam is called.
qu: A vectors of quantiles of interest. Each entry should be in (0, 1).
discrete: If TRUE then covariate discretisation is used for faster model fitting. See mgcv::bam for details.
lsig: The value of the log learning rate used to create the Gibbs posterior. By defauls lsig=NULL and this parameter is estimated by posterior calibration described in Fasiolo et al. (2017). Obviously, the function is much faster if the user provides a value.
err: An upper bound on the error of the estimated quantile curve. Should be in (0, 1). If it is a vector, it should be of the same length of qu. Since qgam v1.3 it is selected automatically, using the methods of Fasiolo et al. (2017). The old default was err=0.05.
multicore: If TRUE the calibration will happen in parallel.
cluster: An object of class c("SOCKcluster", "cluster"). This allowes the user to pass her own cluster, which will be used if multicore == TRUE. The user has to remember to stop the cluster.
ncores: Number of cores used. Relevant if multicore == TRUE.
paropts: a list of additional options passed into the foreach function when parallel computation is enabled. This is important if (for example) your code relies on external data or packages: use the .export and .packages arguments to supply them so that all cluster nodes have the correct environment set up for computing.
control: A list of control parameters. The only one relevant here is link, which is the link function used (see ?elf and ?elflss for defaults). All other control parameters are used by tuneLearnFast. See ?tuneLearnFast for details.
argGam: A list of parameters to be passed to mgcv::gam. This list can potentially include all the arguments listed in ?gam, with the exception of formula, family and data.

Value

A list with entries:

fit = a gamObject, one for each entry of qu. Notice that the slots model and smooth of each object has been removed to save memory. See ?gamObject.
model = the model slot of the gamObjects in the fit slot. This is the same for every fit, hence only one copy is stored.
smooth = the smooth slot of the gamObjects in the fit slot. This is the same for every fit, hence only one copy is stored.
calibr = a list which is the output of an internal call to tuneLearnFast, which is used for calibrating the learning rate. See ?tuneLearnFast for details.

References

Fasiolo, M., Wood, S.N., Zaffran, M., Nedellec, R. and Goude, Y., 2020. Fast calibrated additive quantile regression. Journal of the American Statistical Association (to appear). doi:10.1080/01621459.2020.1725521 .

Fasiolo, M., Wood, S.N., Zaffran, M., Nedellec, R. and Goude, Y., 2021. qgam: Bayesian Nonparametric Quantile Regression Modeling in R. Journal of Statistical Software, 100(9), 1-31, doi:10.18637/jss.v100.i09 .

Author

Matteo Fasiolo <matteo.fasiolo@gmail.com>.

Examples


#####
# Multivariate Gaussian example
####
library(qgam)
set.seed(2)
dat <- gamSim(1, n=300, dist="normal", scale=2)
#> Gu & Wahba 4 term additive model

fit <- mqgam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat, qu = c(0.2, 0.8))
#> Estimating learning rate. Each dot corresponds to a loss evaluation. 
#> qu = 0.2...........done 
#> qu = 0.8...............done 

invisible( qdo(fit, 0.2, plot, pages = 1) )


#####
# Univariate "car" example
####
library(qgam); library(MASS)

# Fit for quantile 0.8 using the best sigma
quSeq <- c(0.2, 0.4, 0.6, 0.8)
set.seed(6436)
fit <- mqgam(accel~s(times, k=20, bs="ad"), data = mcycle, qu = quSeq)
#> Estimating learning rate. Each dot corresponds to a loss evaluation. 
#> qu = 0.4............done 
#> qu = 0.6..........done 
#> qu = 0.2..............done 
#> qu = 0.8.........done 

# Plot the fit
xSeq <- data.frame(cbind("accel" = rep(0, 1e3), "times" = seq(2, 58, length.out = 1e3)))
plot(mcycle$times, mcycle$accel, xlab = "Times", ylab = "Acceleration", ylim = c(-150, 80))
for(iq in quSeq){
  pred <- qdo(fit, iq, predict, newdata = xSeq)
  lines(xSeq$times, pred, col = 2)
}