[R] mixed effects meta-regression: nlme vs. metafor

[Viechtbauer Wolfgang (STAT)]

Hello Christian,

First of all, it's good to see that you are well aware of the fact that lme() without lmeControl(sigma=1) will lead to the estimation of the residual variance component, which implies that the sampling variances specified via varFixed() are only assumed to be known up to a proportionality constant -- however, in the usual meta-analytic models, we assume that the sampling variances are exactly known. In fact, trying to disentangle that residual variance component from any random study effects is usually next to impossible. I mention this explicitly one more time, because I have seen some publications using lme() in exactly this way ...

Indeed, lmeControl(sigma=1) is an option only available in (at least some versions of) S-Plus (I know that it is available in version 6). Of course, that's not that helpful unless you happen to have a copy of S-Plus.

In fact, I ended up developing the metafor package (which started out with a function called mima() that is essentially the predecessor of the rma() function) for that very reason -- I needed a function to fit the meta-analytic random- and mixed-effects models.

As you are also aware of, right now, rma() adds a random effect per observation (i.e., observed effect or outcome), while you want a random effect per study (which of course only matters if you have more than one outcome per study -- as in your example). I have a function in the works that will allow you to do just that. It's not in the package yet, but it will be in the future. This essentially relates back to numerous requests I have received for adding functions to the metafor package that will handle multivariate meta-analytic models, dependent outcomes, and things like network meta-analyses. And to my shame, I have said numerous times: Yes, it's in the works, it will be in the package in the future, don't know yet when ("don't hold your breath"). It's probably just as frustrating for me not to find the time to work as much on the package as I would like as it is for those waiting for me to get around to actually adding that functionality to the package.

I actually have picked up quite a bit of steam in terms of working on the package recently. I am very close to releasing an updated version, the package website (http://www.metafor-project.org/ ) has been totally revamped (and is starting to become actually useful), and a bit of grant money is soon trickling in my direction that involves certain developments in the package. The updated version of the package still does not include that aforementioned function, but I may consider putting a pre-alpha version on the website so that the adventurous are able to try it out.

Alternatively, you could try taking a look at MCMCglmm (http://cran.r-project.org/web/packages/MCMCglmm/index.html), which should be able to fit the model that you want. Can't give you any details on how, but if you get stuck, try posting to R-sig-mixed-models and Jarrod Hadfield (the MCMCglmm package author) is very likely to help you further.

Best,
Wolfgang

--   
Wolfgang Viechtbauer, Ph.D., Statistician   
Department of Psychiatry and Psychology   
School for Mental Health and Neuroscience   
Faculty of Health, Medicine, and Life Sciences   
Maastricht University, P.O. Box 616 (VIJV1)   
6200 MD Maastricht, The Netherlands   
+31 (43) 388-4170 | http://www.wvbauer.com   


> -----Original Message-----
> From: Christian Röver [mailto:christian.roever at med.uni-goettingen.de]
> Sent: Wednesday, January 23, 2013 11:10
> To: r-help at r-project.org
> Cc: Wolfgang Viechtbauer
> Subject: mixed effects meta-regression: nlme vs. metafor
> 
> Hi,
> 
> I would like to do a meta-analysis, i.e., a mixed-effects regression,
> but I don't seem to get what I want using both the nlme or metafor
> packages.
> 
> My question: is there indeed no way to do it?
> And if so, is there another package I could use?
> 
> Here are the details:
> 
> In my meta-analysis I'm comparing different studies that report a
> measure at time zero and after a certain followup time. Each reported
> measurement comes with standard error, and each study uses one (or
> several) of a few treatment categories. I want to fit a random effect
> for each study (the study effect) and a treatment-dependent time effect.
> For the moment I use a linear model, i.e., twice the followup time will
> give you twice the effect, etc...
> 
> I get /almost/ what I want using "nlme" via this command:
> 
>   lme01 <- lme(effect ~ treatment + treatment*time - time - 1,
>                random = ~ 1|study,
>                weights = varFixed(~se2),
>                data = dat)
> 
> "effect" is the real-valued measurement, "treatment" is a factor, and
> "time" is the followup time in months. "se2" is the squared standard
> error.
> Problem is: using the "varFixed()" option, "lme()" will fit an
> additional variance parameter scaling the provided standard errors by a
> certain factor to be estimated. According to some discussions on the
> web, you once were able to prevent the fitting of the extra variance
> parameter in some pre-1998 S-plus versions of nlme using a
> "lmeControl(sigma=1)" option, but this does not appear to available any
> more.
> 
> I again get /almost/ what I want using the "metafor" package:
> 
>   rma01 <- rma(yi = effect,
>                vi = se2,
>                mods = ~ treatment + treatment*time - time - 1,
>                data = dat)
> 
> "rma()" will correctly digest the provided standard errors, but this
> time the problem is that "rma()" will always treat each line in the data
> set as a different study, there does not appear to be a way to tell
> "rma()" that several data points belong to the same study, i.e., have a
> common random effect. What I am missing is an equivalent to the "random"
> statement in the "lme()" command above. Adding an option like
> "slab=as.character(dat$study)" only seems to affect the labeling but not
> the actual computation.
> 
> Any ideas?
> 
> Many thanks in advance,
> 
> Christian Roever