Power for mixed-effects models

library(lattice)
library(lme4)

Reanalysis

Application context: Depression and type of diagnosis

• Reisby et al. (1977) studied the effect of Imipramin on 66 inpatients treated for depression
• Patients’ depression was classified into endogenous and non-endogenous
• Depressive symptoms were measured with the Hamilton depression rating scale
• Depression was measured weekly for 6 time points

Data: reisby.txt

dat      <- read.table("../data/reisby.txt", header = TRUE)
dat$id   <- factor(dat$id)
dat$diag <- factor(dat$diag, levels = c("nonen", "endog"))
dat      <- na.omit(dat)     # drop missing values
head(dat, n = 13)

##     id hamd week  diag endweek
## 1  101   26    0 nonen       0
## 2  101   22    1 nonen       0
## 3  101   18    2 nonen       0
## 4  101    7    3 nonen       0
## 5  101    4    4 nonen       0
## 6  101    3    5 nonen       0
## 7  103   33    0 nonen       0
## 8  103   24    1 nonen       0
## 9  103   15    2 nonen       0
## 10 103   24    3 nonen       0
## 11 103   15    4 nonen       0
## 12 103   13    5 nonen       0
## 13 104   29    0 endog       0

xyplot(hamd ~ week | id, data = dat, type = c("g", "r", "p"),
  pch = 16, strip = strip.custom(bg = "gray96"),
  par.strip.text = list(cex = .8), layout = c(11, 6),
  ylab = "HDRS score", xlab = "Time (week)")

Random-intercept model

\[ \begin{aligned} Y_{ij} &= \beta_0 + \beta_1 \, \mathtt{week}_{ij} + \upsilon_{0i} + \varepsilon_{ij} \\ \upsilon_{0i} &\sim N(0, \sigma^2_{\upsilon_0}) \text{ i.i.d.} \\ \mathbf{\varepsilon}_i &\sim N(0, \, \sigma^2) \text{ i.i.d.} \\ i &= 1, \ldots, I, \quad j = 1, \ldots n_i \end{aligned} \]

m1 <- lmer(hamd ~ week + (1 | id), data = dat, REML = FALSE)
print(m1)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: hamd ~ week + (1 | id)
##    Data: dat
##       AIC       BIC    logLik -2*log(L)  df.resid 
##  2293.189  2308.896 -1142.594  2285.189       371 
## Random effects:
##  Groups   Name        Std.Dev.
##  id       (Intercept) 4.019   
##  Residual             4.363   
## Number of obs: 375, groups:  id, 66
## Fixed Effects:
## (Intercept)         week  
##      23.552       -2.376

Random-slope model

\[ \begin{aligned} Y_{ij} &= \beta_0 + \beta_1 \, \mathtt{week}_{ij} + \upsilon_{0i} + \upsilon_{1i}\, \mathtt{week}_{ij} + \varepsilon_{ij} \\ \begin{pmatrix} \upsilon_{0i}\\ \upsilon_{1i} \end{pmatrix} &\sim N \left(\begin{pmatrix} 0\\ 0 \end{pmatrix}, \, \mathbf{\Sigma}_\upsilon = \begin{pmatrix} \sigma^2_{\upsilon_0} & \sigma_{\upsilon_0 \upsilon_1} \\ \sigma_{\upsilon_0 \upsilon_1} & \sigma^2_{\upsilon_1} \\ \end{pmatrix} \right) \text{ i.i.d.} \\ \mathbf{\varepsilon}_i &\sim N(\mathbf{0}, \, \sigma^2 \mathbf{I}_{n_i}) \text{ i.i.d.} \\ i &= 1, \ldots, I, \quad j = 1, \ldots n_i \end{aligned} \]

m2 <- lmer(hamd ~ week + (week | id), data = dat, REML = FALSE)
print(m2)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: hamd ~ week + (week | id)
##    Data: dat
##       AIC       BIC    logLik -2*log(L)  df.resid 
##  2231.037  2254.599 -1109.519  2219.037       369 
## Random effects:
##  Groups   Name        Std.Dev. Corr  
##  id       (Intercept) 3.554          
##           week        1.442    -0.28 
##  Residual             3.495          
## Number of obs: 375, groups:  id, 66
## Fixed Effects:
## (Intercept)         week  
##      23.577       -2.377

Partial pooling

indiv <- unlist(
  sapply(unique(dat$id),
         function(i) predict(lm(hamd ~ week, dat[dat$id == i, ])))
)
xyplot(hamd + predict(m2, re.form = ~ 0) + predict(m2) + indiv ~ week | id,
  data = dat, type = c("p", "l", "l", "l"),
  pch = 16, strip = strip.custom(bg = "gray96"), grid = TRUE,
  distribute.type = TRUE, par.strip.text = list(cex = .8),
  layout = c(11, 6), ylab = "HDRS score", xlab = "Time (week)")

By-group random-slope model

m3 <- lmer(hamd ~ week * diag + (week | id), data = dat, REML = FALSE)
print(m3)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: hamd ~ week * diag + (week | id)
##    Data: dat
##       AIC       BIC    logLik -2*log(L)  df.resid 
##  2230.929  2262.345 -1107.465  2214.929       367 
## Random effects:
##  Groups   Name        Std.Dev. Corr  
##  id       (Intercept) 3.412          
##           week        1.441    -0.29 
##  Residual             3.495          
## Number of obs: 375, groups:  id, 66
## Fixed Effects:
##    (Intercept)            week       diagendog  week:diagendog  
##       22.47626        -2.36569         1.98802        -0.02706

m3a <- lmer(hamd ~ week + diag + (week | id), data = dat, REML = FALSE)
anova(m3a, m3)

## Data: dat
## Models:
## m3a: hamd ~ week + diag + (week | id)
## m3: hamd ~ week * diag + (week | id)
##     npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
## m3a    7 2228.9 2256.4 -1107.5    2214.9                     
## m3     8 2230.9 2262.3 -1107.5    2214.9 0.0042  1     0.9486

Means and predicted HDRS score by group

dat2 <- aggregate(hamd ~ week + diag, dat, mean)
dat2$m3 <- predict(m3, dat2, re.form = ~ 0)

plot(m3 ~ week, dat2[dat2$diag == "endog", ], type = "l",
     ylim = c(0, 28), xlab = "Week", ylab = "HDRS score")
lines(m3 ~ week, dat2[dat2$diag == "nonen", ], lty = 2)
points(hamd ~ week, dat2[dat2$diag == "endog", ], pch = 16)
points(hamd ~ week, dat2[dat2$diag == "nonen", ], pch = 21, bg = "white")
legend("topright", c("Endogenous", "Non endogenous"),
       lty = 1:2, pch = c(16, 21), bty = "n")

Gory details

fixef(m3)

##    (Intercept)           week      diagendog week:diagendog 
##    22.47626332    -2.36568746     1.98802087    -0.02705576

getME(m3, "theta")

##      id.(Intercept) id.week.(Intercept)             id.week 
##           0.9760823          -0.1175194           0.3951992

t(chol(VarCorr(m3)$id))[lower.tri(diag(2), diag = TRUE)] / sigma(m3)

## [1]  0.9760823 -0.1175194  0.3951992

Power simulation: A new treatment for depression

Setup

## Study design and sample sizes
n_week <- 6
n_subj <- 80
n <- n_week * n_subj
dat <- data.frame(
     id = rep(seq_len(n_subj), each = n_week),
   week = rep(0:(n_week - 1), times = n_subj),
  treat = rep(0:1, each = n/2)
)

## Fixed effects and variance components
beta <- c("(Intercept)" = 23, week = -0.5, treat = 0, "week:treat" = -1)
si <- 3.5
ss <- 1.5
r  <- -0.3
se <- 3.5                # residual sd

## Covariance matrix of random effects
Su <- matrix(c(si^2, r * si * ss, r * si * ss, ss^2), nrow = 2)

## Means from structural part of the model
X <- model.matrix( ~ week * treat, dat)
means <- as.vector(X %*% beta)

Power

pval <- replicate(200, {
  # Data generation
  e <- rnorm(n, mean = 0, sd = se)
  u <- MASS::mvrnorm(n_subj, mu = c(0, 0), Sigma = Su)
  y <- means + u[dat$id, 1] + u[dat$id, 2] * dat$week + e

  # Fitting models to test H0
  m1 <- lmer(y ~ week + treat + (week | id), data = dat, REML = FALSE)
  m2 <- lmer(y ~ week * treat + (week | id), data = dat, REML = FALSE)
  anova(m1, m2)$"Pr(>Chisq)"[2]
})
mean(pval < 0.05)

## [1] 0.745

Parameter recovery

par <- replicate(200, {
  means <- model.matrix( ~ week * treat, dat) %*% beta
  e <- rnorm(n, mean = 0, sd = se)
  u <- MASS::mvrnorm(n_subj, mu = c(0, 0), Sigma = Su)
  y <- means + u[dat$id, 1] + u[dat$id, 2] * dat$week + e

  m1 <- lmer(y ~ week * treat + (1 + week | id), data = dat, REML = FALSE)
  list(
       fixef = fixef(m1),
       theta = getME(m1, "theta"),
       sigma = sigma(m1)
  )
  }, simplify = FALSE
)

rowMeans(sapply(par, function(x) x$fixef))

## (Intercept)        week       treat  week:treat 
##  22.9384614  -0.4527344   0.1297528  -1.0627984

rowMeans(sapply(par, function(x) x$theta))

##      id.(Intercept) id.week.(Intercept)             id.week 
##           0.9814121          -0.1229056           0.3978173

mean(sapply(par, function(x) x$sigma))

## [1] 3.500463

beta

## (Intercept)        week       treat  week:treat 
##        23.0        -0.5         0.0        -1.0

Lt <- chol(Su)
t(Lt)[lower.tri(Lt, diag = TRUE)] / se

## [1]  1.0000000 -0.1285714  0.4088311

se

## [1] 3.5

Reference

Reisby, N., L. F. Gram, P. Bech, A. Nagy, G. O. Petersen, J. Ortmann, I. Ibsen, et al. 1977. “Imipramine: Clinical Effects and Pharmacokinetic Variability.” Psychopharmacology 54: 263–72. https://doi.org/10.1007/BF00426574.