---
title: Bug in nlme::lme with fixed sigma and REML estimation
date: '2016-11-07'
categories:
- Rstats
- programming
- hierarchical models
- nlme
code-tools: true
date-modified: '2024-06-07'
---
About one year ago, the `nlme` package introduced a feature that allowed the user to specify a fixed value for the residual variance in linear mixed effect models fitted with `lme()`. This feature is interesting to me because, when used with the `varFixed()` specification for the residual weights, it allows for estimation of a wide variety of meta-analysis models, including basic random effects models, bivariate models for estimating effects by trial arm, and other sorts of multivariate/multi-level random effects models. However, in kicking the tires on this feature, I noticed that the results that it produces are not quite consistent with the results produced by `metafor`, which is the main package I use for fitting meta-analytic models.
In this post, I document several examples of discrepant estimates between `lme()` and `rma.mv()`, using standard datasets included in the `metafor` package. The main take-aways are:
1. The discrepancies arise only with `REML` estimation (not with `ML` estimation).
2. The discrepancies are present whether or not the `varFixed` specification is used.
3. The discrepancies are mostly small (with minimal impact on the standard errors of the fixed effect estimates), but are larger than I would expect from computational/convergence differences alone.
Another example, based on a different dataset, is documented in [this bug report](https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=16975). Wolfgang Viechtbauer, author of the `metafor` package, identified this problem with `lme` a few months ago already (see his responses in [this thread](https://stat.ethz.ch/pipermail/r-sig-mixed-models/2016q2/024862.html) on the R mixed models mailing list) and noted that the issue was localized to REML estimation. My thanks to Wolfgang for providing feedback on this post.
### Basic random effects model
This example fits a basic random effects model to the BCG vaccine data, available within `metafor`:
```{r, message = FALSE}
library(metafor)
library(nlme)
bcg_example <- function(method = "REML", constant_var = FALSE) {
data(dat.bcg)
dat <- escalc(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
v_bar <- mean(dat$vi)
if (constant_var) dat$vi <- v_bar
# random-effects model using rma.uni()
LOR_uni_fit <- rma(yi, vi, data=dat, method = method)
LOR_uni <- with(LOR_uni_fit,
data.frame(f = "rma.uni",
logLik = logLik(LOR_uni_fit),
est = as.numeric(b),
se = se,
tau = sqrt(tau2)))
# random-effects model using rma.mv()
LOR_mv_fit <- rma.mv(yi, vi, random = ~ 1 | trial, data=dat, method = method)
LOR_mv <- with(LOR_mv_fit,
data.frame(f = "rma.mv",
logLik = logLik(LOR_mv_fit),
est = as.numeric(b),
se = se,
tau = sqrt(sigma2)))
# random-effects model using lme()
if (constant_var) {
LOR_lme_fit <- lme(yi ~ 1, data = dat, method = method,
random = ~ 1 | trial,
control = lmeControl(sigma = sqrt(v_bar)))
tau <- sqrt(as.numeric(coef(LOR_lme_fit$modelStruct$reStruct, unconstrained = FALSE)) * v_bar)
} else {
LOR_lme_fit <- lme(yi ~ 1, data = dat, method = method,
random = ~ 1 | trial,
weights = varFixed(~ vi),
control = lmeControl(sigma = 1))
tau <- sqrt(as.numeric(coef(LOR_lme_fit$modelStruct$reStruct, unconstrained = FALSE)))
}
LOR_lme <- data.frame(f = "lme",
logLik = logLik(LOR_lme_fit),
est = as.numeric(fixef(LOR_lme_fit)),
se = as.numeric(sqrt(vcov(LOR_lme_fit))),
tau = tau)
rbind(LOR_uni, LOR_mv, LOR_lme)
}
bcg_example("REML", constant_var = FALSE)
bcg_example("REML", constant_var = TRUE)
bcg_example("ML", constant_var = FALSE)
bcg_example("ML", constant_var = TRUE)
```
### Bi-variate random effects model
This example fits a bi-variate random effects model, also to the BCG vaccine data:
```{r}
bcg_bivariate <- function(method = "REML", constant_var = FALSE) {
data(dat.bcg)
dat_long <- to.long(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
levels(dat_long$group) <- c("exp", "con")
dat_long$group <- relevel(dat_long$group, ref="con")
dat_long <- escalc(measure="PLO", xi=out1, mi=out2, data=dat_long)
v_bar <- mean(dat_long$vi)
if (constant_var) dat_long$vi <- v_bar
# bivariate random-effects model using rma.mv()
bv_rma_fit <- rma.mv(yi, vi, mods = ~ group,
random = ~ group | study,
struct = "UN", method = method,
data=dat_long)
bv_rma <- with(bv_rma_fit, data.frame(f = "rma.mv",
logLik = logLik(bv_rma_fit),
tau1 = sqrt(tau2[1]),
tau2 = sqrt(tau2[2])))
# bivariate random-effects model using lme()
if (constant_var) {
bv_lme_fit <- lme(yi ~ group, data = dat_long, method = method,
random = ~ group | study,
control = lmeControl(sigma = sqrt(v_bar)))
tau_sq <- colSums(coef(bv_lme_fit$modelStruct$reStruct, unconstrained = FALSE) * matrix(c(1,0,0, 1,2,1), 3, 2)) * v_bar
} else {
bv_lme_fit <- lme(yi ~ group, data = dat_long, method = method,
random = ~ group | study,
weights = varFixed(~ vi),
control = lmeControl(sigma = 1))
tau_sq <- colSums(coef(bv_lme_fit$modelStruct$reStruct, unconstrained = FALSE) * matrix(c(1,0,0, 1,2,1), 3, 2))
}
bv_lme <- data.frame(f = "lme",
logLik = logLik(bv_lme_fit),
tau1 = sqrt(tau_sq[1]),
tau2 = sqrt(tau_sq[2]))
rbind(bv_rma, bv_lme)
}
bcg_bivariate("REML", constant_var = FALSE)
bcg_bivariate("REML", constant_var = TRUE)
bcg_bivariate("ML", constant_var = FALSE)
bcg_bivariate("ML", constant_var = TRUE)
```
### Three-level random-effects model
This example fits a three-level random-effects model to the data from Konstantopoulos (2011):
```{r}
Konstantopoulos <- function(method = "REML", constant_var = FALSE) {
dat <- get(data(dat.konstantopoulos2011))
v_bar <- mean(dat$vi)
if (constant_var) dat$vi <- v_bar
# multilevel random-effects model using rma.mv()
ml_rma_fit <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat, method = method)
ml_rma <- with(ml_rma_fit,
data.frame(f = "rma.mv",
logLik = logLik(ml_rma_fit),
est = as.numeric(b),
se = se,
tau1 = sqrt(sigma2[1]),
tau2 = sqrt(sigma2[2])))
# multilevel random-effects model using lme()
if (constant_var) {
ml_lme_fit <- lme(yi ~ 1, data = dat, method = method,
random = ~ 1 | district / school,
control = lmeControl(sigma = sqrt(v_bar)))
tau <- sqrt(as.numeric(coef(ml_lme_fit$modelStruct$reStruct, unconstrained = FALSE)) * v_bar)
} else {
ml_lme_fit <- lme(yi ~ 1, data = dat, method = method,
random = ~ 1 | district / school,
weights = varFixed(~ vi),
control = lmeControl(sigma = 1))
tau <- sqrt(as.numeric(coef(ml_lme_fit$modelStruct$reStruct, unconstrained = FALSE)))
}
ml_lme <- data.frame(f = "lme",
logLik = logLik(ml_lme_fit),
est = as.numeric(fixef(ml_lme_fit)),
se = as.numeric(sqrt(diag(vcov(ml_lme_fit)))),
tau1 = tau[2],
tau2 = tau[1])
rbind(ml_rma, ml_lme)
}
Konstantopoulos("REML", constant_var = FALSE)
Konstantopoulos("REML", constant_var = TRUE)
Konstantopoulos("ML", constant_var = FALSE)
Konstantopoulos("ML", constant_var = TRUE)
```
### Colophon
```{r}
sessioninfo::session_info()
```