qrcm (version 3.1)

iqrL: 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)

Value

An object of class “iqrL”, a list containing the following items:

theta, phi

estimates of \(\theta\) and \(\phi\).

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 \(\theta\) and \(\phi\), respectively. Note that mf.theta is sorted by increasing id and, within each id, by increasing values of the response variable y, while mf.phi is sorted by increasing id.

s.theta, s.phi

the used ‘s.theta’ and ‘s.phi’ matrices.

fit

a data.frame with the following variables:

• id the cluster identifier.

• y the response variable.

• alpha the estimated individual effects.

• y_alpha = y - alpha[id], the estimated responses purged of the individual effects.

• v estimates of \(V_i\).

• u estimates of \(U_{it}\).

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.

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 fx and fz.

tol

convergence criterion for numerical optimization.

maxit

maximum number of iterations. If missing, a default is computed.

Author

Paolo Frumento

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.

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
  # \donttest{
  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)
  # }