estimate.lmm: Delta Method for a Linear Mixed Model

In bozenne/repeated: Repeated Measurement Models for Discrete Times

Delta Method for a Linear Mixed Model

Description

Estimate standard errors, confidence intervals, and p-values for a smooth transformation of parameters from a linear mixed model.

Usage

## S3 method for class 'lmm'
estimate(
  x,
  f,
  df = !is.null(x$df) & !robust,
  robust = FALSE,
  type.information = NULL,
  level = 0.95,
  method.numDeriv = NULL,
  average = FALSE,
  transform.sigma = NULL,
  transform.k = NULL,
  transform.rho = NULL,
  ...
)

Arguments

x	a lmm object.
f	[function] function taking as input p, the linear mixed model parameters, which outputs the parameter(s) of interest. Can accept extra-arguments.
df	[logical] Should degree-of-freedom, computed using Satterthwaite approximation, for the parameter of interest be output. Can also be a numeric vector providing providing the degrees-of-freedom relative to each estimate.
robust	[logical] Should robust standard errors (aka sandwich estimator) be output instead of the model-based standard errors. Can also be 2 compute the degrees-of-freedom w.r.t. robust standard errors instead of w.r.t. model-based standard errors.
type.information	[character] Should the expected information be used (i.e. minus the expected second derivative) or the observed inforamtion (i.e. minus the second derivative).
level	[numeric,0-1] the confidence level of the confidence intervals.
method.numDeriv	[character] method used to approximate the gradient: either "simple" or "Richardson". Passed to numDeriv::jacobian.
average	[logical] is the estimand the average output of argument f? Otherwise consider each individual output of argument f.
transform.sigma	[character] Transformation used on the variance coefficient for the reference level. One of "none", "log", "square", "logsquare" - see details.
transform.k	[character] Transformation used on the variance coefficients relative to the other levels. One of "none", "log", "square", "logsquare", "sd", "logsd", "var", "logvar" - see details.
transform.rho	[character] Transformation used on the correlation coefficients. One of "none", "atanh", "cov" - see details.
...	extra arguments passed to f.

Details

Based a first order delta method to evaluate the variance of the estimate. The derivative of the transformation is evaluated using numerical differentiation (numDeriv::jacobian).

Argument robust: the Satterhwaite approximation for the degrees-of-freedom of robust standard errors is often unreliable. This is why the default is to use the degrees-of-freedom of the modeled based standard error instead.

Examples

if(require(lava) && require(nlme)){

#### Random effect ####
set.seed(10)
dL <- sampleRem(1e2, n.times = 3, format = "long")
e.lmm1 <- lmm(Y ~ X1+X2+X3 + (1|id), repetition = ~visit|id, data = dL)
nlme::ranef(e.lmm1, se = TRUE)
e.ranef <- estimate(e.lmm1, f  = function(object, p){nlme::ranef(object, p = p)})
e.ranef

if(require(ggplot2)){
df.gg <- cbind(index = 1:NROW(e.ranef), e.ranef)
gg.ranef <- ggplot(df.gg, aes(x = index, y=estimate, ymin=lower, ymax = upper))
gg.ranef + geom_point() + geom_errorbar() + ylab("estimated random effect") + xlab("id")
}

#### ANCOVA via mixed model ####
set.seed(10)
d <- sampleRem(1e2, n.time = 2)
e.ANCOVA1 <- lm(Y2~Y1+X1, data = d)

dL2 <- reshape(d, direction = "long", idvar = c("id","X1"), 
               timevar = "time", times = c("1","2"), varying = c("Y1","Y2"), 
               v.names = "Y")

## estimated treatment effect (no baseline constraint)
e.lmm <- lmm(Y ~ time + time:X1, data = dL2, repetition = ~time|id)

e.delta <- estimate(e.lmm, function(p){
     c(Y1 = p["rho(1,2)"]*p["k.2"],
      X1 = p["time2:X1"]-p["k.2"]*p["rho(1,2)"]*p["time1:X1"])
}) ## same estimate and similar standard errors. 
e.delta ## Degrees-of-freedom are a bit off though
cbind(summary(e.ANCOVA1)$coef, df = df.residual(e.ANCOVA1))

## estimated treatment effect (baseline constraint)
dL2$time2 <- as.numeric(dL2$time=="2")
e.lmmC <- lmm(Y ~ time2 + time2:X1, data = dL2, repetition = ~time|id)
e.deltaC <- estimate(e.lmmC, function(p){
       c(Y1 = p["rho(1,2)"]*p["k.2"],
         X1 = p["time2:X1"])
})
e.deltaC ## Degrees-of-freedom are a bit more accurate
}

bozenne/repeated documentation built on July 16, 2025, 11:16 p.m.