iqrL: Quantile Regression Coefficients Modeling with Longitudinal...
In qrcm: Quantile Regression Coefficients Modeling

Description Usage Arguments Details Value Author(s) References See Also Examples

This function implements Frumento et al's (2021) method for quantile regression coefficients modeling with longitudinal data.

1 2	iqrL(fx, fu = ~ slp(u,3), fz = ~ 1, fv = ~ -1 + I(qnorm(v)), id, weights, s.theta, s.phi, data, tol = 1e-5, maxit)

`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.

New users are recommended to read Frumento and Bottai's (2018) 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]*β(U[it]) + z[i]*γ(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 α[i] = z[i]γ(V[i]) are cluster-level effects with quantile function z[i]γ(v), while x[it]β(u) is the quantile function of Y[it] - α[i].

Both β(u) and γ(v) are modeled parametrically, using a linear combination of known “basis” functions b(u) and c(v) such that

β(u) = β(u | θ) = θ b(u),

γ(u) = γ(u | φ) = φ c(v),

where θ and φ 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 α[i] is assumed to be Normal (fv = ~ -1 + I(qnorm(v))) and to not depend on covariates (fz = ~ 1).

Restrictions on θ and φ are imposed by setting to zero the corresponding elements of s.theta and s.phi.

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

`theta, phi`	estimates of θ and φ.
`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 θ and φ, 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.

Paolo Frumento paolo.frumento@unipi.it

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 [forthcoming].

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.

  ##### 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)