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.roe...@med.uni-goettingen.de]
> Sent: Wednesday, January 23, 2013 11:10
> To: r-help@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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
