LMMstar (version 1.1.0)

residuals: Extract The Residuals From a Linear Mixed Model

Description

Extract or compute the residuals of a linear mixed model.

Usage

# S3 method for lmm
residuals(
  object,
  type = "response",
  variable = NULL,
  at = NULL,
  newdata = NULL,
  p = NULL,
  format = "long",
  keep.data = FALSE,
  fitted.ci = FALSE,
  simplify = TRUE,
  var,
  ...
)
# S3 method for clmm
residuals(object, ...)
# S3 method for mlmm
residuals(object, simplify = TRUE, ...)

Value

lmm: a vector or a data.frame when format="long" (one line per observation, one column per type of residual), a matrix when format="wide" (one line per cluster, one column per timepoint).

Arguments

object

a lmm object.

type

[character] type of residual to output such as raw residuals ("response"), normalized residuals ("normalized", partial residuals ("partial"). Can also be "all" to output all except partial residuals. See detail section.

variable

[character vector] name of the variable relative to which the partial residuals should be computed.

at

[data.frame] values for the covariates at which to evaluate the partial residuals.

newdata

[data.frame] dataset relative to which the residuals should be computed. Only relevant if differs from the dataset used to fit the model.

p

[numeric vector] value of the model coefficients at which to evaluate the residuals. Only relevant if differs from the fitted values.

format

[character] Should the residuals be output in a matrix format with clusters in row and timepoints in columns ("wide"), or in a data.frame/vector with as many rows as observations ("long")

keep.data

[logical] Should the dataset relative to which the residuals are evaluated be output along side the residual values? Only possible in the long format.

fitted.ci

[logical] Should the confidence intervals relative to the fitted values be added to the output. Only relevant when argument keep.data=TRUE.

simplify

[logical] Simplify the data format (vector instead of data.frame) and column names (no mention of the time variable) when possible. Otherwise, information about the call and reference values used for partial residuals be added as an attribute.

var

[Deprecated] Not used anymore, replaced by argument variable.

...

Not used. For compatibility with the generic method.

Details

The argument type defines how the residuals are computed:

• "response": raw residual, i.e. observed outcome minus fitted value \(\varepsilon_{ij} = Y_{ij} - X_{ij} \hat{\beta}\).

• "pearson": each raw residual is divided by its modeled standard deviation \(\varepsilon_{ij} = \frac{Y_{ij} - X_{ij} \hat{\beta}}{\sqrt{\hat{\omega}_{ij}}}\).

• "studentized": same as "pearson" but excluding the contribution of the cluster in the modeled standard deviation \(\varepsilon_{ij} = \frac{Y_{ij} - X_{ij} \hat{\beta}}{\sqrt{\hat{\omega}_{ij}-\hat{q}_{ij}}}\).

• "normalized": raw residuals are multiplied, within clusters, by the inverse of the (upper) Cholesky factor of the modeled residual variance covariance matrix \(\varepsilon_{ij} = ( Y_{i} - X_{i} \hat{\beta} )\hat{C}^{-1}\).

• "normalized2": raw residuals are multiplied, within clusters, by the inverse of the modeled residual variance covariance matrix \(\varepsilon_{ij} = ( Y_{i} - X_{i} \hat{\beta} )\hat{Omega}^{-1}\).

• "scaled": scaled residuals (see PROC MIXED in SAS). Numerically identical to "normalized" but computed by sequentially scaling and centering the residuals, to make them conditionally independent of previous residuals from the same cluster at previous repetitions.

• "partial": partial residuals (\(\gamma E + \hat{\varepsilon}\)). A reference level can be also be specified via the attribute "reference" to change the absolute level of the partial residuals. "partial-center": partial residuals with centered continuous covariates (\(\gamma E + \hat{\varepsilon}\) where \(E\) has been centered, i.e., has 0-mean)

where

• \(X=(E,W)\) the design matrix. For partial residuals, it is split according to the variable(s) in argument variable (\(E\)) and the rest (\(W\)).

• \(Y\) the outcome

• \(\hat{\beta}=(\hat{\gamma},\hat{\delta})\) the estimated mean coefficients relative to \(X=(E,W)\)

• \(\hat{\Omega}\) the modeled variance-covariance of the residuals and \(\hat{\omega}\) its diagonal elements

• \(\hat{C}\) the upper Cholesky factor of \(\hat{\Omega}\), i.e. upper triangular matrix satisfying \(\hat{C}^{t} \hat{C} = \hat{\Omega}\)

• \(\hat{Q}_i= X_i (X^{t}\hat{\Omega}X)^{-1}X_i^{t}\) a cluster specific correction factor, approximating the contribution of cluster i to \(\hat{\Omega}\). Its diagonal elements are denoted \(\hat{q}_i\).

• \(\hat{D}_i\) the upper Cholesky factor of \(\hat{\Omega}-\hat{Q}_i\)

Setting argument fitted.ci to TRUE, simplify to FALSE, format to "long" returns an attribute "grad" containing the first order partial derivatives of the residuals with respect to the model parameters.

Examples

#### simulate data in the long format ####
set.seed(10)
dL <- sampleRem(100, n.times = 3, format = "long")

#### Linear Model ####
e.lm <- lmm(Y ~ visit + X1 + X2 + X6, data = dL)

## partial residuals
pRes <- residuals(e.lm, type = "partial", variable = "X6")
range(residuals(e.lm) + dL$X6 * coef(e.lm)["X6"] - pRes)
e.reslm <- residuals(e.lm, type = "partial", variable = "X6", keep.data = TRUE, simplify = FALSE)
plot(e.reslm)

## partial residuals with specific reference
residuals(e.lm, type = "partial", variable = "X1",
          at = data.frame(visit=factor(2,1:3),X2=0,X6=3))

## partial residuals with centered covariates
residuals(e.lm, type = "partial-center", variable = "X1")

#### Linear Mixed Model ####
eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5 + X6,
               repetition = ~visit|id, structure = "UN", data = dL)

## residuals
e.resL <- residuals(eUN.lmm, type = "normalized",
                    keep.data = TRUE, simplify = FALSE)
plot(e.resL, type = "qqplot")
plot(e.resL, type = "scatterplot", labeller = ggplot2::label_both)
e.resW <- residuals(eUN.lmm, format = "wide", type = "normalized",
                    simplify = FALSE)
plot(e.resW, type = "correlation")

## residuals and predicted values
residuals(eUN.lmm, type = "all")
residuals(eUN.lmm, type = "all", keep.data = TRUE)

## partial residuals
residuals(eUN.lmm, type = "partial", variable = c("(Intercept)","visit","X6"))