R help - Jan 2013 - 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

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