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:
The discrepancies arise only with REML estimation (not with ML estimation).
The discrepancies are present whether or not the varFixed specification is used.
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. Wolfgang Viechtbauer, author of the metafor package, identified this problem with lme a few months ago already (see his responses in this thread 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:
Code
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)
f logLik est se tau
1 rma.uni -12.57566 -0.7451778 0.1860279 0.5811816
2 rma.mv -12.57566 -0.7451778 0.1860280 0.5811818
3 lme -13.34043 -0.7471979 0.1916902 0.6030524
Code
bcg_example("REML", constant_var =TRUE)
f logLik est se tau
1 rma.uni -12.96495 -0.7716272 0.1977007 0.5911451
2 rma.mv -12.96495 -0.7716272 0.1977007 0.5911452
3 lme -15.62846 -0.7716272 0.1899448 0.5571060
Code
bcg_example("ML", constant_var =FALSE)
f logLik est se tau
1 rma.uni -13.07276 -0.7419668 0.1779534 0.5499605
2 rma.mv -13.07276 -0.7419669 0.1779534 0.5499608
3 lme -13.07276 -0.7419668 0.1779534 0.5499605
Code
bcg_example("ML", constant_var =TRUE)
f logLik est se tau
1 rma.uni -13.525084 -0.7716272 0.1899447 0.5571059
2 rma.mv -13.525084 -0.7716272 0.1899447 0.5571059
3 lme -2.479133 -0.7716272 0.1899447 0.5571060
Bi-variate random effects model
This example fits a bi-variate random effects model, also to the BCG vaccine data:
---title: Bug in nlme::lme with fixed sigma and REML estimationdate: '2016-11-07'categories:- Rstats- programming- hierarchical models- nlmecode-tools: truedate-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 modelThis 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 modelThis 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 modelThis 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()```