| Title: | Linear Quantile Mixed Models |
|---|---|
| Description: | Functions to fit quantile regression models for hierarchical data (2-level nested designs) as described in Geraci and Bottai (2014, Statistics and Computing) <doi:10.1007/s11222-013-9381-9>. A vignette is given in Geraci (2014, Journal of Statistical Software) <doi:10.18637/jss.v057.i13> and included in the package documents. The packages also provides functions to fit quantile models for independent data and for count responses. |
| Authors: | Marco Geraci |
| Maintainer: | Marco Geraci <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.5.8 |
| Built: | 2024-11-29 08:47:03 UTC |
| Source: | CRAN |
Fit quantile regression models for independent and hierarchical data
| Package: | lqmm |
| Type: | Package |
| Version: | 1.5.8 |
| Date: | 2022-04-05 |
| License: | GPL (>=2) |
| LazyLoad: | yes |
Marco Geraci
Maintainer: Marco Geraci <[email protected]>
Geraci M (2014). Linear quantile mixed models: The lqmm package for Laplace quantile regression. Journal of Statistical Software, 57(13), 1–29. <doi:10.18637/jss.v057.i13>
Geraci M and Bottai M (2007). Quantile regression for longitudinal data using the asymmetric Laplace distribution. Biostatistics 8(1), 140–154. <doi:10.1093/biostatistics/kxj039>
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479. <doi:10.1007/s11222-013-9381-9>.
This function is used to obtain a bootstrap sample of a fitted LQM or LQMM. It is a generic function.
boot(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE) ## S3 method for class 'lqm' boot(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE) ## S3 method for class 'lqmm' boot(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE)boot(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE) ## S3 method for class 'lqm' boot(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE) ## S3 method for class 'lqmm' boot(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE)
object |
an object of |
R |
number of bootstrap replications. |
seed |
optional random number generator seed. |
startQR |
logical flag. If |
An object of class boot.lqm is a data frame with R rows and npars columns containing the bootstrap estimates of theta. If object contains results for multiple quantiles, boot.lqm returns an array of dimension c(R,npars,nt), where nt is the length of tau.
An object of class boot.lqmm is a data frame with R rows and npars columns containing the bootstrap estimates of theta_x, theta_z, and scale. If object contains results for multiple quantiles, boot.lqmm returns an array of dimension c(R,npars,nt), where nt is the length of tau. The elements of theta_z are labelled with reStruct. See function covHandling and the example below on how to derive the variance-covariance matrix of the random effects starting from theta_z.
The following attributes are available:
tau |
index of the quantile(s). |
estimated |
the estimated parameter as given by |
R |
number of bootstrap replications. |
seed |
the random number generator seed used to produce the bootstrap sample. |
npars |
total numer of parameters. |
rdf |
the number of residual degrees of freedom. |
indices |
the bootstrap sample of independent data units. |
Marco Geraci
# boot.lqm set.seed(123) n <- 500 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) fit.lqm <- lqm(y ~ x, data = test, tau = 0.5) fit.boot <- boot(fit.lqm) str(fit.boot) # boot.lqmm data(Orthodont) fit <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = 0.5, data = Orthodont) fit.boot <- boot(fit) str(fit.boot)# boot.lqm set.seed(123) n <- 500 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) fit.lqm <- lqm(y ~ x, data = test, tau = 0.5) fit.boot <- boot(fit.lqm) str(fit.boot) # boot.lqmm data(Orthodont) fit <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = 0.5, data = Orthodont) fit.boot <- boot(fit) str(fit.boot)
coef extracts model coefficients from lqm, lqm.counts objects.
## S3 method for class 'lqm' coef(object, ...)## S3 method for class 'lqm' coef(object, ...)
object |
an |
... |
not used. |
a vector for single quantiles or a matrix for multiple quantiles.
Marco Geraci
coef extracts model coefficients from lqmm objects.
## S3 method for class 'lqmm' coef(object, ...)## S3 method for class 'lqmm' coef(object, ...)
object |
a fitted object of |
... |
not used. |
a vector for single quantiles or a matrix for multiple quantiles.
Marco Geraci
This is an auxiliary function.
covHandling(theta, n, cov_name, quad_type)covHandling(theta, n, cov_name, quad_type)
theta |
unique parameters of the variance-covariance matrix of the random effects as returned by |
n |
dimension of the vector of random effects. |
cov_name |
see argument |
quad_type |
type of quadrature "c("normal","robust")". |
Marco Geraci
Density, distribution function, quantile function and random generation for the asymmetric Laplace distribution.
dal(x, mu = 0, sigma = 1, tau = 0.5, log = FALSE) pal(x, mu = 0, sigma = 1, tau = 0.5) qal(x, mu = 0, sigma = 1, tau = 0.5) ral(n, mu = 0, sigma = 1, tau = 0.5)dal(x, mu = 0, sigma = 1, tau = 0.5, log = FALSE) pal(x, mu = 0, sigma = 1, tau = 0.5) qal(x, mu = 0, sigma = 1, tau = 0.5) ral(n, mu = 0, sigma = 1, tau = 0.5)
x |
vector of quantiles ( |
n |
number of observations. |
mu |
location parameter. |
sigma |
positive scale parameter. |
tau |
skewness parameter (0,1). |
log |
logical; if |
The asymmetric Laplace distribution with parameters (mu, sigma, tau) has density
Marco Geraci
This generic function extracts the fixed and random components of bootstrapped estimates of an lqmm object.
extractBoot(object, which = "fixed") ## S3 method for class 'boot.lqmm' extractBoot(object, which = "fixed")extractBoot(object, which = "fixed") ## S3 method for class 'boot.lqmm' extractBoot(object, which = "fixed")
object |
an object of |
which |
character indicating whether |
The "random" parameters refer to the "raw" parameters of the variance-covariance matrix of the random effects as returned by lqmm.fit.gs and lqmm.fit.df.
a matrix of bootstrapped estimates.
Marco Geraci
boot.lqmm, lqmm.fit.gs, lqmm.fit.df
## Orthodont data data(Orthodont) # Random intercept model fit <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = 0.5, data = Orthodont) fit.boot <- boot(fit) # extract fixed effects B <- extractBoot(fit.boot, which = "fixed") # covariance matrix estimated fixed parameters cov(B)## Orthodont data data(Orthodont) # Random intercept model fit <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = 0.5, data = Orthodont) fit.boot <- boot(fit) # extract fixed effects B <- extractBoot(fit.boot, which = "fixed") # covariance matrix estimated fixed parameters cov(B)
This function calculates nodes and weights for Gaussian quadrature. See help("gauss.quad") from package statmod.
Original version by Gordon Smyth
Gordon Smyth with contributions from Yifang Hu, Peter Dunn and Belinda Phipson. (2011). statmod: Statistical Modeling. R package version 1.4.11. https://CRAN.R-project.org/package=statmod
This function calculates nodes and weights for Gaussian quadrature in terms of probability distributions. See help("gauss.quad.prob") from package statmod.
Original version by Gordon Smyth
Gordon Smyth with contributions from Yifang Hu, Peter Dunn and Belinda Phipson. (2011). statmod: Statistical Modeling. R package version 1.4.11. https://CRAN.R-project.org/package=statmod
This function tests whether all eigenvalues of a symmetric matrix are positive. See help("is.positive.definite") from package corpcor.
Original version by Korbinian Strimmer
Juliane Schaefer, Rainer Opgen-Rhein, Verena Zuber, A. Pedro Duarte Silva and Korbinian Strimmer. (2011). corpcor: Efficient Estimation of Covariance and (Partial) Correlation. R package version 1.6.0. https://CRAN.R-project.org/package=corpcor
The labor data frame has 358 rows and 4 columns of the
change in pain over time for several 83 women in labor.
This data frame contains the following columns:
an ordered factor indicating the subject on which the
measurement was made. The levels are labelled 1
to 83.
a numeric vector of self–reported pain scores on a 100mm line.
a dummy variable with values
1 for subjects who received a pain medication and
0 for subjects who received a placebo.
a numeric vector of times (minutes since randomization) at which pain was measured.
The labor pain data were reported by Davis (1991) and successively analyzed by Jung (1996) and Geraci and Bottai (2007). The data set consists of repeated measurements of self–reported amount of pain on N = 83 women in labor, of which 43 were randomly assigned to a pain medication group and 40 to a placebo group. The response was measured every 30 min on a 100–mm line, where 0 means no pain and 100 means extreme pain. A nearly monotone pattern of missing data was found for the response variable and the maximum number of measurements for each woman was six.
Davis CS (1991). Semi–parametric and non–parametric methods for the analysis of repeated measurements with applications to clinical trials. Statistics in Medicine 10, 1959–80.
Geraci M and Bottai M (2007). Quantile regression for longitudinal data using the asymmetric Laplace distribution. Biostatistics 8(1), 140–154.
Jung S (1996). Quasi–likelihood for median regression models. Journal of the American Statistical Association 91, 251–7.
logLik.lqm extracts the log-likelihood of a fitted LQM.
## S3 method for class 'lqm' logLik(object, ...)## S3 method for class 'lqm' logLik(object, ...)
object |
an object of |
... |
not used. |
Marco Geraci
logLik.lqmm extracts the log-likelihood of a fitted LQMM.
## S3 method for class 'lqmm' logLik(object, ...)## S3 method for class 'lqmm' logLik(object, ...)
object |
an object of |
... |
not used. |
Marco Geraci
lqm is used to fit linear quantile models based on the asymmetric Laplace distribution.
lqm(formula, data, subset, na.action, weights = NULL, tau = 0.5, contrasts = NULL, control = list(), fit = TRUE)lqm(formula, data, subset, na.action, weights = NULL, tau = 0.5, contrasts = NULL, control = list(), fit = TRUE)
formula |
an object of class |
data |
an optional data frame, list or environment (or object coercible by |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
na.action |
a function which indicates what should happen when the data contain NAs. The default is set by the |
weights |
An optional vector of weights to be used in the fitting process. |
tau |
the quantile(s) to be estimated. This must be a number between 0 and 1, otherwise the execution is stopped. If more than one quantile is specified, rounding off to the 4th decimal must give non–duplicated values of |
contrasts |
an optional list. See the contrasts.arg of |
control |
list of control parameters of the fitting process. See |
fit |
logical flag. If |
The function computes an estimate on the tau-th quantile function of the response, conditional on the covariates, as specified by the formula argument. The quantile predictor is assumed to be linear. The function maximizes the (log)likelihood of a Laplace regression which is equivalent to the minimization of the weighted sum of absolute residuals (Koenker and Bassett, 1978). The optimization algorithm is based on the gradient of the Laplace log–likelihood (Bottai, Orsini and Geraci, 2013).
lqm returns an object of class lqm.
The function summary is used to obtain and print a summary of the results.
An object of class lqm is a list containing the following components:
theta |
a vector of coefficients. |
scale |
the scale parameter. |
gradient |
the gradient. |
logLik |
the log–likelihood. |
opt |
details on optimization (see |
call |
the matched call. |
term.labels |
names for theta. |
terms |
the terms object used. |
nobs |
the number of observations. |
edf, dim_theta
|
the length of theta. |
rdf |
the number of residual degrees of freedom. |
tau |
the estimated quantile(s). |
x |
the model matrix. |
y |
the model response. |
weights |
the weights used in the fitting process (a vector of 1's if |
InitialPar |
starting values for theta. |
control |
list of control parameters used for optimization (see |
Updates/FAQ/news are published here https://marcogeraci.wordpress.com/. New versions are usually published here https://github.com/marco-geraci/lqmm/ before going on CRAN.
Marco Geraci
Bottai M, Orsini N, Geraci M (2015). A Gradient Search Maximization Algorithm for the Asymmetric Laplace Likelihood, Journal of Statistical Computation and Simulation, 85(10), 1919-1925.
Chen C (2007). A finite smoothing algorithm for quantile regression. Journal of Computational and Graphical Statistics, 16(1), 136-164.
Koenker R and Bassett G (1978). Regression Quantiles. Econometrica 46(1), 33–50.
summary.lqm, coef.lqm, predict.lqm, residuals.lqm
set.seed(123) n <- 500 p <- 1:3/4 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) fit.lqm <- lqm(y ~ x, data = test, tau = p, control = list(verbose = FALSE, loop_tol_ll = 1e-9), fit = TRUE) fit.lqmset.seed(123) n <- 500 p <- 1:3/4 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) fit.lqm <- lqm(y ~ x, data = test, tau = p, control = list(verbose = FALSE, loop_tol_ll = 1e-9), fit = TRUE) fit.lqm
This function is used to fit a quantile regression model when the response is a count variable.
lqm.counts(formula, data, weights = NULL, offset = NULL, contrasts = NULL, tau = 0.5, M = 50, zeta = 1e-05, B = 0.999, cn = NULL, alpha = 0.05, control = list())lqm.counts(formula, data, weights = NULL, offset = NULL, contrasts = NULL, tau = 0.5, M = 50, zeta = 1e-05, B = 0.999, cn = NULL, alpha = 0.05, control = list())
formula |
an object of class |
data |
an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which lqm is called. |
weights |
an optional vector of weights to be used in the fitting process. |
offset |
an optional offset to be included in the model frame. |
contrasts |
an optional list. See the |
tau |
quantile to be estimated. |
M |
number of dithered samples. |
zeta |
small constant (see References). |
B |
right boundary for uniform random noise U[0,B] to be added to the response variable (see References). |
cn |
small constant to be passed to |
alpha |
significance level. |
control |
list of control parameters of the fitting process. See |
A linear quantile regression model if fitted to the log–transformed response. Additional tranformation functions will be implemented. The notation used here follows closely that of Machado and Santos Silva (2005).
an object of class "lqm.counts" containing the following components
tau |
the estimated quantile. |
theta |
regression quantile (on the log–scale). |
fitted |
predicted quantile (on the response scale). |
tTable |
coefficients, standard errors, etc. |
x |
the model matrix. |
y |
the model response. |
offset |
offset. |
nobs |
the number of observations. |
M |
specified number of dithered samples for standard error estimation. |
Mn |
actual number of dithered samples used for standard error estimation that gave an invertible D matrix (Machado and Santos Silva, 2005). |
term.labels |
names for theta. |
terms |
the terms object used. |
rdf |
the number of residual degrees of freedom. |
InitialPar |
starting values for theta. |
control |
list of control parameters used for optimization (see |
Marco Geraci
Machado JAF and Santos Silva JMC (2005). Quantiles for counts. Journal of the American Statistical Association, 100(472), 1226–1237.
n <- 100 x <- runif(n) test <- data.frame(x = x, y = rpois(n, 2*x)) lqm.counts(y ~ x, data = test, M = 50)n <- 100 x <- runif(n) test <- data.frame(x = x, y = rpois(n, 2*x)) lqm.counts(y ~ x, data = test, M = 50)
This function controls the arguments to be passed to routines written in C for LQM estimation. The optimization algorithm is based on the gradient of the Laplace log–likelihood (Bottai, Orsini and Geraci, 2013).
lqm.fit.gs(theta, x, y, weights, tau, control)lqm.fit.gs(theta, x, y, weights, tau, control)
theta |
starting values for the regression coefficients. |
x |
the model matrix. |
y |
the model response. |
weights |
the weights used in the fitting process. |
tau |
the quantile to be estimated. |
control |
list of control parameters used for optimization (see |
See argument fit in lqm for generating a list of arguments to be called by this function.
An object of class list containing the following components:
theta |
a vector of coefficients. |
scale |
the scale parameter. |
gradient |
the gradient. |
logLik |
the log–likelihood. |
opt |
number of iterations when the estimation algorithm stopped. |
.
Marco Geraci
Bottai M, Orsini N, Geraci M (2014). A Gradient Search Maximization Algorithm for the Asymmetric Laplace Likelihood, Journal of Statistical Computation and Simulation, 85, 1919-1925.
set.seed(123) n <- 500 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) lqm.ls <- lqm(y ~ x, data = test, fit = FALSE) do.call("lqm.fit.gs", lqm.ls)set.seed(123) n <- 500 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) lqm.ls <- lqm(y ~ x, data = test, fit = FALSE) do.call("lqm.fit.gs", lqm.ls)
A list of parameters for controlling the fitting process.
lqmControl(method = "gs1", loop_tol_ll = 1e-5, loop_tol_theta = 1e-3, check_theta = FALSE, loop_step = NULL, beta = 0.5, gamma = 1.25, reset_step = FALSE, loop_max_iter = 1000, smooth = FALSE, omicron = 0.001, verbose = FALSE)lqmControl(method = "gs1", loop_tol_ll = 1e-5, loop_tol_theta = 1e-3, check_theta = FALSE, loop_step = NULL, beta = 0.5, gamma = 1.25, reset_step = FALSE, loop_max_iter = 1000, smooth = FALSE, omicron = 0.001, verbose = FALSE)
method |
character vector that specifies which code to use for carrying out the gradient search algorithm: "gs1" (default) based on C code and "gs2" based on R code. Method "gs3" uses a smoothed loss function. See details. |
loop_tol_ll |
tolerance expressed as relative change of the log-likelihood. |
loop_tol_theta |
tolerance expressed as relative change of the estimates. |
check_theta |
logical flag. If |
loop_step |
step size (default standard deviation of response). |
beta |
decreasing step factor for line search (0,1). |
gamma |
nondecreasing step factor for line search (>= 1). |
reset_step |
logical flag. If |
loop_max_iter |
maximum number of iterations. |
smooth |
logical flag. If |
omicron |
small constant for smoothing the loss function when using |
verbose |
logical flag. |
The methods "gs1" and "gs2" implement the same algorithm (Bottai et al, 2015). The former is based on C code, the latter on R code. While the C code is faster, the R code seems to be more efficient in handling large datasets. For method "gs2", it is possible to replace the classical non-differentiable loss function with a smooth version (Chen, 2007).
a list of control parameters.
Marco Geraci
Bottai M, Orsini N, Geraci M (2015). A Gradient Search Maximization Algorithm for the Asymmetric Laplace Likelihood, Journal of Statistical Computation and Simulation, 85(10), 1919-1925.
Chen C (2007). A finite smoothing algorithm for quantile regression. Journal of Computational and Graphical Statistics, 16(1), 136-164.
lqmm is used to fit linear quantile mixed models based on the asymmetric Laplace distribution.
lqmm(fixed, random, group, covariance = "pdDiag", tau = 0.5, nK = 7, type = "normal", rule = 1, data = sys.frame(sys.parent()), subset, weights, na.action = na.fail, control = list(), contrasts = NULL, fit = TRUE)lqmm(fixed, random, group, covariance = "pdDiag", tau = 0.5, nK = 7, type = "normal", rule = 1, data = sys.frame(sys.parent()), subset, weights, na.action = na.fail, control = list(), contrasts = NULL, fit = TRUE)
fixed |
an object of class |
random |
a one-sided formula of the form |
group |
grouping factor. |
covariance |
variance–covariance matrix of the random effects. Default is |
tau |
the quantile(s) to be estimated. |
nK |
number of quadrature knots. |
type |
type of quadrature "c("normal","robust")" (see details). |
rule |
quadrature rule (see details). |
data |
an optional data frame containing the variables named in
|
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
weights |
an optional vector of weights to be used in the fitting process of the same length as the number of rows of |
na.action |
a function that indicates what should happen when the
data contain |
control |
list of control parameters of the fitting process. See |
contrasts |
not yet implemented. |
fit |
logical flag. If FALSE the function returns a list of arguments to be passed to |
The function computes an estimate on the tau-th quantile function of the response, conditional on the covariates, as specified by the formula argument, and on random effects, as specified by the random argument. The quantile predictor is assumed to be linear. The function maximizes the (log)likelihood of the Laplace regression proposed by Geraci and Bottai (2014). The likelihood is numerically integrated via Gaussian quadrature techniques. The optimization algorithm is based on the gradient of the Laplace log–likelihood (control = list(method = "gs")). An alternative optimization algorithm is based on a Nelder-Mead algorithm (control = list(method = "df")) via optim. The scale parameter is optimized in a refinement step via optimize.
Quadrature approaches include Gauss-Hermite (type = "normal") and Gauss-Laguerre (type = "robust") quadrature. The argument rule takes one of the following: 1 (product rule quadrature), 2 (sparse grid quadrature), 3 (nested quadrature rule - only for type = "normal"), 4 (quadrature rule with the smallest number of nodes between rules 1 or 2). Rules 2 and 3 have not yet been tested extensively.
Different standard types of positive–definite matrices for the random effects can be specified: pdIdent multiple of an identity; pdCompSymm compound symmetry structure (constant diagonal and constant off–diagonal elements); pdDiag diagonal; pdSymm general positive–definite matrix, with no additional structure.
Weights are given to clusters, therefore it is expected that these are constant within cluster. When the weights are specified in the main call, then the first value by group in the vector weights will be replicated for the same length of each group. Alternatively, different weights within the same cluster can be introduced with a direct call to lqmm.fit.gs or lqmm.fit.df.
The lqmm vignette can be accessed by typing help(package = "lqmm") and then following the link 'User guides, package vignettes and other documentation'.
lqmm returns an object of class lqmm.
The function summary is used to obtain and print a summary of the results.
An object of class lqmm is a list containing the following components:
theta |
a vector containing fixed regression coefficients and parameters of the variance–covariance matrix of the random effects. See |
theta_x, theta_z
|
partition of |
scale |
the scale parameter. |
gradient |
the gradient ( |
logLik |
the log–likelihood. |
opt |
details on optimization (see |
call |
the matched call. |
nn |
column names of |
mm |
column names of |
nobs |
the number of observations. |
dim_theta |
the number of columns in |
dim_theta_z |
the length of |
edf |
length of |
rdf |
the number of residual degrees of freedom. |
df |
edf + 1 (scale parameter). |
tau |
the estimated quantile(s). |
mmf |
the model matrix – fixed effects. |
mmr |
the model matrix – random effects. |
y |
the model response. |
revOrder |
original order of observations (now ordered according to |
weights |
the likelihood weights used in the fitting process (a vector of 1's if |
group |
the grouping factor. |
ngroups |
the number of groups. |
QUAD |
quadrature nodes and weights. |
type |
the type of quadrature. |
rule |
quadrature rule. |
InitialPar |
starting values for theta. |
control |
list of control parameters used for optimization (see |
cov_name |
class of variance-covariance matrix for the random effects. |
mfArgs |
arguments for |
Updates/FAQ/news are published here https://marcogeraci.wordpress.com/. New versions are usually published here https://github.com/marco-geraci/lqmm/ before going on CRAN.
Marco Geraci
Genz A, and Keister BD (1996). Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight. Journal of Computational and Applied Mathematics, 71(2), 299–309. <doi:10.1016/0377-0427(95)00232-4>
Geraci M (2014). Linear quantile mixed models: The lqmm package for Laplace quantile regression. Journal of Statistical Software, 57(13), 1–29. <doi:10.18637/jss.v057.i13>
Geraci M and Bottai M (2007). Quantile regression for longitudinal data using the asymmetric Laplace distribution. Biostatistics 8(1), 140–154. <doi:10.1093/biostatistics/kxj039>
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479. <doi:10.1007/s11222-013-9381-9>.
Heiss F, and Winschel V (2008). Likelihood approximation by numerical integration on sparse grids. Journal of Econometrics, 144(1), 62–80. <doi:10.1016/j.jeconom.2007.12.004>
lqm, summary.lqmm, coef.lqmm, VarCorr.lqmm, predict.lqmm, residuals.lqmm
# Test example set.seed(123) M <- 50 n <- 10 test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) fit.lqmm <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, tau = 0.5, nK = 11, type = "normal") fit.lqmm #Call: lqmm(fixed = y ~ x, random = ~1, group = group, tau = 0.5, nK = 11, # type = "normal", data = test) #Quantile 0.5 #Fixed effects: #(Intercept) x # 3.443 9.258 #Covariance matrix of the random effects: #(Intercept) # 3.426 #Residual scale parameter: 0.8697 (standard deviation 2.46) #Log-likelihood: -1178 #Number of observations: 500 #Number of groups: 50 ## Orthodont data data(Orthodont) # Random intercept model fitOi.lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = c(0.1,0.5,0.9), data = Orthodont) coef(fitOi.lqmm) # Random slope model fitOs.lqmm <- lqmm(distance ~ age, random = ~ age, group = Subject, tau = c(0.1,0.5,0.9), cov = "pdDiag", data = Orthodont) # Extract estimates VarCorr(fitOs.lqmm) coef(fitOs.lqmm) ranef(fitOs.lqmm) # AIC AIC(fitOi.lqmm) AIC(fitOs.lqmm)# Test example set.seed(123) M <- 50 n <- 10 test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) fit.lqmm <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, tau = 0.5, nK = 11, type = "normal") fit.lqmm #Call: lqmm(fixed = y ~ x, random = ~1, group = group, tau = 0.5, nK = 11, # type = "normal", data = test) #Quantile 0.5 #Fixed effects: #(Intercept) x # 3.443 9.258 #Covariance matrix of the random effects: #(Intercept) # 3.426 #Residual scale parameter: 0.8697 (standard deviation 2.46) #Log-likelihood: -1178 #Number of observations: 500 #Number of groups: 50 ## Orthodont data data(Orthodont) # Random intercept model fitOi.lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = c(0.1,0.5,0.9), data = Orthodont) coef(fitOi.lqmm) # Random slope model fitOs.lqmm <- lqmm(distance ~ age, random = ~ age, group = Subject, tau = c(0.1,0.5,0.9), cov = "pdDiag", data = Orthodont) # Extract estimates VarCorr(fitOs.lqmm) coef(fitOs.lqmm) ranef(fitOs.lqmm) # AIC AIC(fitOi.lqmm) AIC(fitOs.lqmm)
This function controls the arguments to be passed to optim and optimize for LQMM estimation.
lqmm.fit.df(theta_0, x, y, z, weights, cov_name, V, W, sigma_0, tau, group, control)lqmm.fit.df(theta_0, x, y, z, weights, cov_name, V, W, sigma_0, tau, group, control)
theta_0 |
starting values for the linear predictor. |
x |
the model matrix for fixed effects (see details). |
y |
the model response (see details). |
z |
the model matrix for random effects (see details). |
weights |
the weights used in the fitting process (see details). |
cov_name |
variance–covariance matrix of the random effects. Default is |
V |
nodes of the quadrature. |
W |
weights of the quadrature. |
sigma_0 |
starting value for the scale parameter. |
tau |
the quantile(s) to be estimated. |
group |
the grouping factor (see details). |
control |
list of control parameters used for optimization (see |
In lqmm, see argument fit for generating a list of arguments to be called by this function; see argument covariance for alternative variance–covariance matrices.
NOTE: the data should be ordered by group when passed to lqmm.fit.df (such ordering is performed by lqmm).
An object of class "list" containing the following components:
theta |
a vector of coefficients, including the "raw" variance–covariance parameters (see |
scale |
the scale parameter. |
logLik |
the log–likelihood. |
opt |
number of iterations when the estimation algorithm stopped for lower (theta) and upper (scale) loop. |
.
Marco Geraci
set.seed(123) M <- 50 n <- 10 test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) lqmm.ls <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, fit = FALSE) do.call("lqmm.fit.df", lqmm.ls)set.seed(123) M <- 50 n <- 10 test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) lqmm.ls <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, fit = FALSE) do.call("lqmm.fit.df", lqmm.ls)
This function controls the arguments to be passed to routines written in C for LQMM estimation. The optimization algorithm is based on the gradient of the Laplace log–likelihood (Bottai, Orsini and Geraci, 2014; Geraci and Bottai, 2014).
lqmm.fit.gs(theta_0, x, y, z, weights, cov_name, V, W, sigma_0, tau, group, control)lqmm.fit.gs(theta_0, x, y, z, weights, cov_name, V, W, sigma_0, tau, group, control)
theta_0 |
starting values for the linear predictor. |
x |
the model matrix for fixed effects (see details). |
y |
the model response (see details). |
z |
the model matrix for random effects (see details). |
weights |
the weights used in the fitting process (see details). |
cov_name |
variance–covariance matrix of the random effects. Default is |
V |
nodes of the quadrature. |
W |
weights of the quadrature. |
sigma_0 |
starting value for the scale parameter. |
tau |
the quantile(s) to be estimated. |
group |
the grouping factor (see details). |
control |
list of control parameters used for optimization (see |
In lqmm, see argument fit for generating a list of arguments to be called by this function; see argument covariance for alternative variance–covariance matrices.
NOTE: the data should be ordered by group when passed to lqmm.fit.gs (such ordering is performed by lqmm).
An object of class "list" containing the following components:
theta |
a vector of coefficients, including the "raw" variance–covariance parameters (see |
scale |
the scale parameter. |
gradient |
the gradient. |
logLik |
the log–likelihood. |
opt |
number of iterations when the estimation algorithm stopped for lower (theta) and upper (scale) loop. |
.
Marco Geraci
Bottai M, Orsini N, Geraci M. (2014). A gradient search maximization algorithm for the asymmetric Laplace likelihood, Journal of Statistical Computation and Simulation (in press).
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479.
set.seed(123) M <- 50 n <- 10 test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) lqmm.ls <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, fit = FALSE) do.call("lqmm.fit.gs", lqmm.ls)set.seed(123) M <- 50 n <- 10 test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) lqmm.ls <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, fit = FALSE) do.call("lqmm.fit.gs", lqmm.ls)
A list of parameters for controlling the fitting process.
lqmmControl(method = "gs", LP_tol_ll = 1e-5, LP_tol_theta = 1e-5, check_theta = FALSE, LP_step = NULL, beta = 0.5, gamma = 1, reset_step = FALSE, LP_max_iter = 500, UP_tol = 1e-4, UP_max_iter = 20, startQR = FALSE, verbose = FALSE)lqmmControl(method = "gs", LP_tol_ll = 1e-5, LP_tol_theta = 1e-5, check_theta = FALSE, LP_step = NULL, beta = 0.5, gamma = 1, reset_step = FALSE, LP_max_iter = 500, UP_tol = 1e-4, UP_max_iter = 20, startQR = FALSE, verbose = FALSE)
method |
character vector that specifies the estimation method: "gs" for gradient search (default) and "df" for Nelder-Mead. |
LP_tol_ll |
tolerance expressed as absolute change of the log-likelihood. |
LP_tol_theta |
tolerance expressed as absolute change of |
check_theta |
logical flag. If TRUE the algorithm performs an additional check on the change in the estimates. |
LP_step |
step size (default standard deviation of response). |
beta |
decreasing step factor for line search (0,1). |
gamma |
nondecreasing step factor for line search (>= 1). |
reset_step |
logical flag. If |
LP_max_iter |
maximum number of iterations |
UP_tol |
tolerance expressed as absolute change of the |
UP_max_iter |
maximum number of iterations. |
startQR |
logical flag. If |
verbose |
logical flag. |
LP (lower loop) refers to the estimation of regression coefficients and variance-covariance parameters. UP (upper loop) refers to the estimation of the scale parameter.
a list of control parameters.
Marco Geraci
This function computes the nearest positive definite of a real symmetric matrix.
See help("make.positive.definite") from package corpcor.
Original version by Korbinian Strimmer
Juliane Schaefer, Rainer Opgen-Rhein, Verena Zuber, A. Pedro Duarte Silva and Korbinian Strimmer. (2011). corpcor: Efficient Estimation of Covariance and (Partial) Correlation. R package version 1.6.0. https://CRAN.R-project.org/package=corpcor
Accessory functions.
meanAL(mu, sigma, tau) varAL(sigma, tau) invvarAL(x, tau)meanAL(mu, sigma, tau) varAL(sigma, tau) invvarAL(x, tau)
mu |
location parameter. |
sigma |
scale parameter. |
tau |
skewness parameter. |
x |
numeric value. |
meanAL computes the mean of an asymmetric Laplace with parameters mu, sigma and tau.
varAL computes the variance of an asymmetric Laplace with parameters sigma and tau.
invvarAL computes the scale parameter of an asymmetric Laplace with parameter tau and variance x.
Marco Geraci
Yu K and Zhang J (2005). A three-parameter asymmetric Laplace distribution and its extension. Communications in Statistics-Theory and Methods 34, 1867–1879.
This function estimates the parameters of an asymmetric Laplace distribution for a sample.
mleAL(x)mleAL(x)
x |
a numeric vector. |
an object of class list containing the following components:
m |
location parameter |
sigma |
scale parameter |
tau |
skewness parameter |
r |
number of iterations |
Marco Geraci
Yu K and Zhang J (2005). A three-parameter asymmetric Laplace distribution and its extension. Communications in Statistics-Theory and Methods 34, 1867–1879.
The Orthodont data frame has 108 rows and 4 columns of the
change in an orthdontic measurement over time for several young subjects.
This data frame contains the following columns:
a numeric vector of distances from the pituitary to the pterygomaxillary fissure (mm). These distances are measured on x-ray images of the skull.
a numeric vector of ages of the subject (yr).
an ordered factor indicating the subject on which the
measurement was made. The levels are labelled M01
to M16 for the males and F01 to F13 for
the females. The ordering is by increasing average distance
within sex.
a factor with levels
Male and
Female
Investigators at the University of North Carolina Dental School followed the growth of 27 children (16 males, 11 females) from age 8 until age 14. Every two years they measured the distance between the pituitary and the pterygomaxillary fissure, two points that are easily identified on x-ray exposures of the side of the head.
Pinheiro, J. C. and Bates, D. M. (2000), Mixed-Effects Models in S and S-PLUS, Springer, New York. (Appendix A.17)
Potthoff, R. F. and Roy, S. N. (1964), “A generalized multivariate analysis of variance model useful especially for growth curve problems”, Biometrika, 51, 313–326.
Jose Pinheiro, Douglas Bates, Saikat DebRoy, Deepayan Sarkar and the R Development Core Team (2011). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-100. https://CRAN.R-project.org/package=nlme
This function computes predictions based on fitted linear quantile model.
## S3 method for class 'lqm' predict(object, newdata, interval = FALSE, level = 0.95, na.action = na.pass, ...) ## S3 method for class 'lqm.counts' predict(object, newdata, na.action = na.pass, ...)## S3 method for class 'lqm' predict(object, newdata, interval = FALSE, level = 0.95, na.action = na.pass, ...) ## S3 method for class 'lqm.counts' predict(object, newdata, na.action = na.pass, ...)
object |
an |
newdata |
an optional data frame in which to look for variables with which to predict. If omitted, the fitted values are used. |
interval |
logical flag. If |
level |
confidence level. This argument is for |
na.action |
function determining what should be done with missing values in |
... |
further arguments passed to |
a vector or a matrix or an array of predictions.
Marco Geraci
residuals.lqm, residuals.lqm.counts, lqm, lqm.counts, coef.lqm, boot.lqm
lqmm Object
The predictions at level 0 correspond to predictions based only on the fixed effects estimates. The predictions at level 1 are obtained by adding the best linear predictions of the random effects to the predictions at level 0. See details for interpretation. The function predint will produce 1-alpha confidence intervals based on bootstrap centiles.
## S3 method for class 'lqmm' predict(object, newdata, level = 0, na.action = na.pass, ...) ## S3 method for class 'lqmm' predint(object, level = 0, alpha = 0.05, R = 50, seed = round(runif(1, 1, 10000)))## S3 method for class 'lqmm' predict(object, newdata, level = 0, na.action = na.pass, ...) ## S3 method for class 'lqmm' predint(object, level = 0, alpha = 0.05, R = 50, seed = round(runif(1, 1, 10000)))
object |
an |
newdata |
an optional data frame in which to look for variables with which to predict. If omitted, the fitted values are produced. |
level |
an optional integer vector giving the level of grouping to be used in obtaining the predictions. |
na.action |
function determining what should be done with missing values in |
alpha |
1- |
R |
number of bootstrap replications. |
seed |
optional random number generator seed. |
... |
not used. |
As discussed by Geraci and Bottai (2014), integrating over the random effects will give "weighted averages" of the cluster-specific quantile effects. These may be interpreted strictly as population regression quantiles for the median (tau=0.5) only. Therefore, predictions at the population level (code=0) should be interpreted analogously.
a vector or a matrix of predictions for predict.lqmm. A data frame or a list of data frames for predint.lqmm containing predictions, lower and upper bounds of prediction intervals, and standard errors.
Marco Geraci
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479.
## Orthodont data data(Orthodont) # Random intercept model fitOi.lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = c(0.1,0.5,0.9), data = Orthodont) # Predict (y - Xb) predict(fitOi.lqmm, level = 0) # Predict (y - Xb - Zu) predict(fitOi.lqmm, level = 1) # 95% confidence intervals predint(fitOi.lqmm, level = 0, alpha = 0.05)## Orthodont data data(Orthodont) # Random intercept model fitOi.lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = c(0.1,0.5,0.9), data = Orthodont) # Predict (y - Xb) predict(fitOi.lqmm, level = 0) # Predict (y - Xb - Zu) predict(fitOi.lqmm, level = 1) # 95% confidence intervals predint(fitOi.lqmm, level = 0, alpha = 0.05)
Print an object generated by lqm or lqm.counts.
## S3 method for class 'lqm' print(x, digits = max(6, getOption("digits")), ...)## S3 method for class 'lqm' print(x, digits = max(6, getOption("digits")), ...)
x |
an |
digits |
a non-null value for digits specifies the minimum number of significant digits to be printed in values. |
... |
not used. |
Marco Geraci
lqmm Object
Print an object generated by lqmm.
## S3 method for class 'lqmm' print(x, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'lqmm' print(x, digits = max(3, getOption("digits") - 3), ...)
x |
an |
digits |
a non-null value for digits specifies the minimum number of significant digits to be printed in values. |
... |
not used. |
Marco Geraci
lqm Summary Object
Print summary of an lqm object.
## S3 method for class 'summary.lqm' print(x, ...)## S3 method for class 'summary.lqm' print(x, ...)
x |
a |
... |
not used. |
Marco Geraci
lqmm Summary Object
Print summary of an lqmm object.
## S3 method for class 'summary.lqmm' print(x, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'summary.lqmm' print(x, digits = max(3, getOption("digits") - 3), ...)
x |
a |
digits |
a non-null value for digits specifies the minimum number of significant digits to be printed in values. |
... |
not used. |
Marco Geraci
This function computes random effects for a linear quantile mixed model.
## S3 method for class 'lqmm' ranef(object, ...)## S3 method for class 'lqmm' ranef(object, ...)
object |
an object of |
... |
not used. |
The prediction of the random effects is done via estimated best linear prediction (Geraci and Bottai, 2014). The generic function ranef is imported from the nlme package (Pinheiro et al, 2014).
a data frame or a list of data frames of predicted random effects.
Marco Geraci
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479. doi: 10.1007/s11222-013-9381-9.
Pinheiro J, Bates D, DebRoy S, Sarkar D and R Core Team (2014). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-117, https://CRAN.R-project.org/package=nlme.
This function computes the residuals from a fitted linear quantile model.
## S3 method for class 'lqm' residuals(object, ...)## S3 method for class 'lqm' residuals(object, ...)
object |
an |
... |
not used. |
a vector or matrix of residuals.
Marco Geraci
lqm, lqm.counts, predict.lqm, coef.lqm
lqmm Object
The residuals at level 0 correspond to population residuals (based only on the fixed effects estimates). The residuals at level 1 are obtained by adding the best linear predictions of the random effects to the predictions at level 0 and the subtracting these from the model response.
## S3 method for class 'lqmm' residuals(object, level = 0, ...)## S3 method for class 'lqmm' residuals(object, level = 0, ...)
object |
an |
level |
an optional integer vector giving the level of grouping to be used in obtaining the predictions. Level zero corresponds to the population residuals. |
... |
not used. |
a matrix of residuals.
Marco Geraci
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479. doi: 10.1007/s11222-013-9381-9.
lqmm, predict.lqmm, coef.lqmm, ranef.lqmm,
boot.lqm Object
Summary method for class boot.lqm.
## S3 method for class 'boot.lqm' summary(object, alpha = 0.05, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'boot.lqm' summary(object, alpha = 0.05, digits = max(3, getOption("digits") - 3), ...)
object |
an object of |
alpha |
numeric value for the interval confidence level ( |
digits |
a non-null value for digits specifies the minimum number of significant digits to be printed in values. |
... |
not used. |
Marco Geraci
boot.lqmm Object
This function gives a summary of a botstrapped lqmm object
## S3 method for class 'boot.lqmm' summary(object, alpha = 0.05, digits = max(3, getOption("digits") - 3), ...)## S3 method for class 'boot.lqmm' summary(object, alpha = 0.05, digits = max(3, getOption("digits") - 3), ...)
object |
an object of |
alpha |
numeric value for the interval confidence level ( |
digits |
a non-null value for digits specifies the minimum number of significant digits to be printed in values. |
... |
not used. |
Marco Geraci
Geraci M and Bottai M (2014). Linear quantile mixed models. Statistics and Computing, 24(3), 461–479. doi: 10.1007/s11222-013-9381-9.
lqm ObjectSummary method for class lqm.
## S3 method for class 'lqm' summary(object, method = "boot", alpha = 0.05, covariance = FALSE, ...)## S3 method for class 'lqm' summary(object, method = "boot", alpha = 0.05, covariance = FALSE, ...)
object |
an object of |
method |
specifies the method used to compute standard errors: "boot" for bootstrap (default), "nid" for large sample approximations under nid assumptions. |
alpha |
significance level. |
covariance |
logical flag. If |
... |
see |
print.summary.lqm formats the coefficients, standard errors, etc. and additionally gives ‘significance stars’.
an object of class summary.lqm. The function summary.lqm computes and returns a list of summary statistics of the fitted linear quantile mixed model given in object, using the components (list elements) from its argument, plus
Cov |
the covariance matrix obtained from the bootstrapped estimates (if |
tTable |
a matrix with estimates, standard errors, etc. |
Marco Geraci
The code for the "nid" method has been adapted from the function summary.rq in package quantreg. It depends on the function bandwidth.rq.
Roger Koenker (2016). quantreg: Quantile Regression. R package version 5.29. https://CRAN.R-project.org/package=quantreg
set.seed(12356) n <- 200 p <- 1:3/4 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) fit.lqm <- lqm(y ~ x, data = test, tau = p) summary(fit.lqm, R = 50)set.seed(12356) n <- 200 p <- 1:3/4 test <- data.frame(x = runif(n,0,1)) test$y <- 30 + test$x + rnorm(n) fit.lqm <- lqm(y ~ x, data = test, tau = p) summary(fit.lqm, R = 50)
lqmm ObjectSummary method for class lqmm.
## S3 method for class 'lqmm' summary(object, method = "boot", alpha = 0.05, covariance = FALSE, ...)## S3 method for class 'lqmm' summary(object, method = "boot", alpha = 0.05, covariance = FALSE, ...)
object |
an object of |
method |
specifies the method used to compute standard errors. Currently, only the bootstrap method ("boot") is available. |
alpha |
significance level. |
covariance |
logical flag. If |
... |
see |
print.summary.lqmm formats the coefficients, standard errors, etc. and additionally gives ‘significance stars’.
an object of class summary.lqmm. The function summary.lqmm computes and returns a list of summary statistics of the fitted linear quantile mixed model given in object, using the components (list elements) from its argument, plus
Cov |
the covariance matrix obtained from the bootstrapped estimates (if |
tTable |
a matrix with estimates, standard errors, etc. |
B |
the matrix of all bootstrapped parameters. |
Marco Geraci
data(Orthodont) fitOi.lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = c(0.1,0.5,0.9), data = Orthodont) summary(fitOi.lqmm)data(Orthodont) fitOi.lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject, tau = c(0.1,0.5,0.9), data = Orthodont) summary(fitOi.lqmm)
This function extracts the variance-covariance matrix of the random effects from a fitted lqmm object.
## S3 method for class 'lqmm' VarCorr(x, sigma = NULL, ...)## S3 method for class 'lqmm' VarCorr(x, sigma = NULL, ...)
x |
an object of |
sigma |
not used. |
... |
not used. |
This function returns the variance or the variance-covariance matrix of the random effects. It calls covHandling to manage the output of lqmm.fit.gs or lqmm.fit.df. A post-fitting approximation to the nearest positive (semi)definite matrix (Higham, 2002) is applied if necessary. The generic function VarCorr is imported from the nlme package (Pinheiro et al, 2014).
Marco Geraci
Higham N (2002). Computing the Nearest Correlation Matrix - A Problem from Finance. IMA Journal of Numerical Analysis, 22, 329-343.
Pinheiro J, Bates D, DebRoy S, Sarkar D and R Core Team (2014). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-117, https://CRAN.R-project.org/package=nlme.