| iqrL {qrcm} | R Documentation |
Quantile Regression Coefficients Modeling with Longitudinal Data
Description
This function implements Frumento et al's (2021) method for quantile regression coefficients modeling with longitudinal data.
Usage
iqrL(fx, fu = ~ slp(u,3), fz = ~ 1, fv = ~ -1 + I(qnorm(v)),
id, weights, s.theta, s.phi, data, tol = 1e-5, maxit)
Arguments
fx, fu, fz, fv |
formulas that describe the model (see ‘Details’). |
id |
a vector of cluster identifiers. |
weights |
an optional vector of weights to be used in the fitting process. |
s.theta, s.phi |
optional 0/1 matrices that permit excluding some model coefficients. |
data |
an optional data frame, list or environment containing the variables in |
tol |
convergence criterion for numerical optimization. |
maxit |
maximum number of iterations. If missing, a default is computed. |
Details
New users are recommended to read Frumento and Bottai's (2016) paper for details
on notation and modeling, and to have some familiarity with the
iqr command, of which iqrL is a natural expansion.
The following data-generating process is assumed:
Y_{it} = x_{it}\beta(U_{it}) + z_i\gamma(V_i)
where x_{it} are level-1 covariates, z_i are level-2 covariates,
and (U_{it}, V_i) are independent U(0,1) random variables.
This model implies that \alpha_i = z_i\gamma(V_i) are cluster-level
effects with quantile function z_i\gamma(v), while x_{it}\beta(u)
is the quantile function of Y_{it} - \alpha_i.
Both \beta(u) and \gamma(v)
are modeled parametrically, using a linear combination of known “basis”
functions b(u) and c(v) such that
\beta(u) = \beta(u | \theta) = \theta b(u),
\gamma(u) = \gamma(u | \phi) = \phi c(v),
where \theta and \phi are matrices of model parameters.
Model specification is implemented as follows.
-
fx is a two-sided formula of the form y ~ x.
-
fu is a one-sided formula that describes
b(u). -
fz is a one-sided formula of the form ~ z.
-
fv is a one-sided formula that describes
c(v).
By default, fu = ~ slp(u,3), a shifted Legendre's polynomial (see slp),
and the distribution of \alpha_i is assumed to be Normal (fv = ~ -1 + I(qnorm(v)))
and to not depend on covariates (fz = ~ 1).
Restrictions on \theta and \phi are imposed by setting to zero the corresponding elements
of s.theta and s.phi.
Value
An object of class “iqrL”, a list containing the following items:
theta, phi |
estimates of |
obj.function |
the value of the minimized loss function, and, separately, the level-1 and the level-2 loss. The number of model parameters (excluding the individual effects) is returned as an attribute. |
call |
the matched call. |
converged |
logical. The convergence status. |
n.it |
the number of iterations. |
covar.theta, covar.phi |
the estimated covariance matrices. |
mf.theta, mf.phi |
the model frames used to fit |
s.theta, s.phi |
the used ‘s.theta’ and ‘s.phi’ matrices. |
fit |
a data.frame with the following variables:
Observations are sorted by increasing id and, within each id, by increasing y. |
Use summary.iqrL, plot.iqrL, and predict.iqrL
for summary information, plotting, and predictions from the fitted model. The function
test.fit.iqrL can be used for goodness-of-fit assessment.
The generic accessory functions coefficients, formula, terms, model.matrix, vcov are available to extract information from the fitted model.
Author(s)
Paolo Frumento paolo.frumento@unipi.it
References
Frumento, P., and Bottai, M. (2016). Parametric modeling of quantile regression coefficient functions. Biometrics, 72 (1), 74-84.
Frumento, P., Bottai, M., and Fernandez-Val, I. (2021). Parametric modeling of quantile regression coefficient functions with longitudinal data. Journal of the American Statistical Association, 116 (534), 783-797.
See Also
summary.iqrL, plot.iqrL, predict.iqrL,
for summary, plotting, and prediction, and test.fit.iqrL for goodness-of-fit assessment.
plf and slp to define b(u) or c(v)
to be piecewise linear functions and shifted Legendre polynomials, respectively.
Examples
##### Also see ?iqr for a tutorial on modeling
##### Using simulated data in all examples
##### Example 1
n <- 1000 # n. of observations
n.id <- 100 # n. of clusters
id <- rep(1:n.id, each = n/n.id) # cluster id
x1 <- runif(n) # a level-1 covariate
z1 <- rbinom(n,1,0.5)[id] # a level-2 covariate
V <- runif(n.id) # V_i
U <- runif(n) # U_it
alpha <- (0.5 + z1)*qnorm(V) # or alpha = rnorm(n.id, 0, 0.5 + z1)
y_alpha <- qexp(U) + 3*x1 # or y_alpha = 3*x1 + rexp(n)
y <- y_alpha + alpha[id] # observed outcome
mydata <- data.frame(id = id, y = y, x1 = x1, z1 = z1[id])
# true quantile function: beta0(u) + beta1(u)*x1 + gamma0(v) + gamma1(v)*z1
# beta0(u) = qexp(u)
# beta1(u) = 3
# gamma0(v) = 0.5*qnorm(v)
# gamma1(v) = qnorm(v)
##### Example 1 (cont.) fitting the model
model1 <- iqrL(fx = y ~ x1, fu = ~ I(qexp(u)), fz = ~ z1, fv = ~ -1 + I(qnorm(v)),
id = id, data = mydata)
summary(model1) # theta, phi
summary(model1, level = 1, p = c(0.1,0.9)) # beta
summary(model1, level = 2, p = c(0.1,0.9)) # gamma
par(mfrow = c(2,2)); plot(model1, ask = FALSE)
##### Example 1 (cont.) - excluding coefficients
s.theta <- rbind(0:1,1:0) # beta0(u) has no intercept, and beta1(u) does not depend on u.
model2 <- iqrL(fx = y ~ x1, fu = ~ I(qexp(u)), fz = ~ z1, fv = ~ -1 + I(qnorm(v)),
id = id, s.theta = s.theta, data = mydata)
summary(model2)
test.fit(model2) # testing goodness-of-fit
##### Example 1 (cont.) - flexible modeling using slp for lev. 1, asymm. logistic for lev. 2
model3 <- iqrL(fx = y ~ x1, fu = ~ slp(u,3),
fz = ~ z1, fv = ~ -1 + I(log(2*v)) + I(-log(2*(1 - v))),
id = id, data = mydata)
par(mfrow = c(2,2)); plot(model3, ask = FALSE)
##### Example 2 - revisiting the classical linear random-effects model
n <- 1000 # n. of observations
n.id <- 100 # n. of clusters
id <- rep(1:n.id, each = n/n.id) # id
x1 <- runif(n,0,5)
E <- rnorm(n) # level-1 error
W <- rnorm(n.id, 0, 0.5) # level-2 error
y <- 2 + 3*x1 + E + W[id] # linear random-intercept model
s.theta <- rbind(1, 1:0)
linmod <- iqrL(fx = y ~ x1, fu = ~ I(qnorm(u)), id = id, s.theta = s.theta)
summary(linmod)