[R-meta] How to conduct a meta-analysis on multiple-treatment studies with a repeated measure designs?
Gram, Gil (IITA)
G@Gr@m @ending from cgi@r@org
Thu May 24 10:52:38 CEST 2018
Dear all,
I found the discussion of Koenraad, Gabriele, and Wolfgang very interesting because it applies to my case as well.
However, it didnt answer all of my questions. Im also trying to perform a multivariate mixed effect model with repeated measures over time with rma.mv. Hopefully you can shed some light on it...
Here is a simplified exerpt of my data:
dat1 = data.table(site = c('Chinyika', 'Chinyika', 'Chinyika', 'Chinyika', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2',
'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2', 'Dodoma2'),
time = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2),
treatment = c('Control', 'MR', 'OR', 'ORMR', 'Control', 'MR', 'MR', 'OR', 'OR', 'ORMR', 'ORMR', 'ORMR', 'ORMR', 'Control', 'MR', 'MR', 'OR', 'OR', 'ORMR', 'ORMR', 'ORMR', 'ORMR'),
idGroup = c(34,34,34,34,270,270,270,270,270,270,270,270,270,271,271,271,271,271,271,271,271,271),
yield = c(0.75,3.41,1.66,3.61,0.85,1.1,1.5,1.6,2.25,1.8,2.1,2.3,2.5,0.6,1.25,1.7,0.95,1.25,1.4,1.7,1.9,2),
sd = c(0.51,2.04,0.98,2.25,0.15,0.15,0.15,0.18,0.10,0.10,0.13,0.18,0.30,0.05,0.10,0.10,0.10,0.08,0.10,0.05,0.10,0.10),
nRep = c(6,6,6,6,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3),
orgRate = c(0,0,NA,NA,0,0,0,65,130,65,130,65,130,0,0,0,65,130,65,130,65,130),
minRate = c(0,NA,0,NA,0,40,80,0,0,40,40,80,80,0,40,80,0,0,40,40,80,80))
I have studies with one or more sites, over 1 or more years, and with essentially 4 fertiliser treatments. But in practice more treatments are possible if
different fertiliser rates are used (as shown in the example data set, by the organic Rate and mineral Rate columns).
According to https://rdrr.io/cran/metafor/man/escalc.html (Measures for Quantitative Variables) I use the escalc to calculate the sampling variances (vi)
and use measure = 'MN' because I want the raw mean of each group to be the observed outcome for the meta-analysis
dat = escalc(measure="MN", mi=yield, sdi=sd, ni=nRep, data=dat1)
Then following the example of Gleser & Olkin (2009) and Gabriele/Wolfgang, I should compute a var-covar matrix in order to account for the multiple (correlated) treatments.
How should I adapt the function below for a varying number of treatments (minimum 4, maximum 9 for dat1)? is that even possible?
calc.v <- function(x) {
v <- matrix(1/x$n2i[1] + outer(x$yi, x$yi, "*")/(2*x$Ni[1]), nrow=nrow(x), ncol=nrow(x))
diag(v) <- x$vi
v
}
V <- bldiag(lapply(split(dat, dat$idGroup), calc.v))
V
Thus, the model I think I need is the following:
rma.mv(yi, vi, mods = ~ treatment, random = ~ treatment | site/time, struct="UN", data=dat)
but with vi that should be replaced with V for multiple (correlated) treatments
struct = "UN" because random effects are expected to have different variances for each outcome and might be correlated.
random = ~ treatment | site/time, because the effect of treatment is likely to be depeninding on site and will differ over time (random slope),
and to account for correlations within sites, but also whithin time within sites.
However, I understood from Wolfgang's reply that adding random effects wont be sufficient, but that I should account for covariances
in the sampling errors due to the repeated measurements as wel. But how to do that? integrate a new V matrix?
Weighting: so far we have used sample errors in our model, that I think rma.mv is using as weights (?), but how can we tell the model to also
account for number of years or study errors in the weighting?
lastly, the output interpretation: I'm interested in 2 things, the (fixed) effect of the treatments, and the variance for each treatment.
When I run the model above without time in the random structure, I get a variance for each treatment:
Multivariate Meta-Analysis Model (k = 22; method: REML)
Variance Components:
outer factor: site (nlvls = 2)
inner factor: treatment (nlvls = 4)
estim sqrt k.lvl fixed level
tau^2.1 0.0070 0.0840 3 no Control
tau^2.2 1.7843 1.3358 5 no MR
tau^2.3 0.0183 0.1352 5 no OR
tau^2.4 1.4943 1.2224 9 no ORMR
rho.Cntr rho.MR rho.OR rho.ORMR Cntr MR OR ORMR
Control 1 1.0000 1.0000 1.0000 - no no no
MR 1.0000 1 1.0000 1.0000 2 - no no
OR 1.0000 1.0000 1 1.0000 2 2 - no
ORMR 1.0000 1.0000 1.0000 1 2 2 2 -
Test for Residual Heterogeneity:
QE(df = 18) = 482.3498, p-val < .0001
Test of Moderators (coefficient(s) 2:4):
QM(df = 3) = 532.2245, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.6816 0.0677 10.0681 <.0001 0.5489 0.8143 ***
treatmentMR 1.6370 0.9284 1.7633 0.0778 -0.1825 3.4566 .
treatmentOR 0.8680 0.0551 15.7403 <.0001 0.7599 0.9761 ***
treatmentORMR 1.9295 0.8442 2.2855 0.0223 0.2748 3.5841 *
I think this can be interpreted as the treatment Control has the lowest yield variability (tau^2.1), and OR the highest (tau^2.2).
Meaning that the yield under the MR treatment is (not significantly, p 0.05) higher, but less reliabel than the control (?)
when I run the model with time in the random structure, I don't get a variance for each treatment. So how do I assess then yield variability?
Multivariate Meta-Analysis Model (k = 22; method: REML)
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0763 0.2763 4 no treatment
sigma^2.2 0.0000 0.0001 8 no treatment/site
sigma^2.3 0.1288 0.3589 12 no treatment/site/time
Test for Residual Heterogeneity:
QE(df = 18) = 482.3498, p-val < .0001
Test of Moderators (coefficient(s) 2:4):
QM(df = 3) = 6.4643, p-val = 0.0911
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.7297 0.3525 2.0700 0.0385 0.0388 1.4206 *
treatmentMR 0.8082 0.5113 1.5805 0.1140 -0.1940 1.8104
treatmentOR 0.8916 0.5039 1.7694 0.0768 -0.0960 1.8792 .
treatmentORMR 1.2450 0.5113 2.4348 0.0149 0.2428 2.2472 *
Any help is greatly appreciated! thanks…
> On 19May, 18, at 13:00, r-sig-meta-analysis-request at r-project.org wrote:
>
> Send R-sig-meta-analysis mailing list submissions to
> r-sig-meta-analysis at r-project.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
> or, via email, send a message with subject or body 'help' to
> r-sig-meta-analysis-request at r-project.org
>
> You can reach the person managing the list at
> r-sig-meta-analysis-owner at r-project.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-meta-analysis digest..."
>
>
> Today's Topics:
>
> 1. Re: How to conduct a meta-analysis on multiple-treatment
> studies with a repeated measure designs? (Viechtbauer, Wolfgang (SP))
> 2. When to remove the intercept and other questions (Rafael Rios)
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Fri, 18 May 2018 14:02:16 +0000
> From: "Viechtbauer, Wolfgang (SP)"
> <wolfgang.viechtbauer at maastrichtuniversity.nl>
> To: Koenraad van Meerbeek <koenraad.vanmeerbeek at bios.au.dk>,
> "r-sig-meta-analysis at r-project.org"
> <r-sig-meta-analysis at r-project.org>
> Subject: Re: [R-meta] How to conduct a meta-analysis on
> multiple-treatment studies with a repeated measure designs?
> Message-ID: <6ed0c395aa0041eb976c0df020c16e00 at UM-MAIL3214.unimaas.nl>
> Content-Type: text/plain; charset="utf-8"
>
> In that case, you ideally would want to account for covariances in the sampling errors due to the repeated measurements as well -- adding random effects is not sufficient then. If computing the covariances is not possible, then you could browse through the archives to see suggestions on how to deal with this issue (this is pretty much the #1 issue that comes up on this mailing list).
>
> P.S.: Koenraad, please actually sign up for the mailing list. Every time you try to post, one of the admins has to manually approve your post before it goes through.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: Koenraad van Meerbeek [mailto:koenraad.vanmeerbeek at bios.au.dk]
> Sent: Friday, 18 May, 2018 15:55
> To: Viechtbauer, Wolfgang (SP); r-sig-meta-analysis at r-project.org
> Cc: Gabriele Midolo
> Subject: RE: [R-meta] How to conduct a meta-analysis on multiple-treatment studies with a repeated measure designs?
>
> Dear all,
>
> The units of analysis stay the same across the years, so we´ll definitely want to add this to the model. We already succeeded in doing the analysis according to Gabriele´s suggestions. I will have a look at those new examples you mentioned.
>
> Thanks for all the help.
>
> Koenraad
>
> -----Original Message-----
> From: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl>
> Sent: vrijdag 18 mei 2018 15:47
> To: r-sig-meta-analysis at r-project.org
> Cc: Koenraad van Meerbeek <koenraad.vanmeerbeek at bios.au.dk>; Gabriele Midolo <gabriele.midolo at gmail.com>
> Subject: RE: [R-meta] How to conduct a meta-analysis on multiple-treatment studies with a repeated measure designs?
>
> Hi Koenraad,
>
> With this type of data structure, you ideally would not just account for the covariances in the sampling errors due to repeated use of the control group information, but also due to measurements taken over multiple years. The former induces covariance for sure. Whether the latter is also an issue for your data is hard to say without knowing more about what you are studying. For example, if I measure the same patients multiple times, then these measurements will tend to be correlated, which in turn induces correlation between effect size measures computed with these measurements. But it sounds like you are in ecology, in which case you are probably taking rather different types of measurements. If the units of analysis are actually different across the years, then there might not be covariance in the sampling errors over time after all. However, one would still want to allow for covariance in the underlying true effects. This can be accomplished by adding appropriate random effects to the model. As Gabriele mentioned, for example something like 'random = ~ Study / Year / ID' (where ID is just a unique value for each row in the dataset). One could also do fancier things allowing for autocorrelation in the random effects over time, since one commonly sees decaying correlation as the time difference increases. In rma.mv(), there are the "AR", "HAR", and "CAR" structures for this purpose. Examples based on this can be found under help(dat.fine1993) and help(dat.ishak2007). There is also:
>
> Musekiwa, A., Manda, S. O., Mwambi, H. G., & Chen, D. G. (2016). Meta-analysis of effect sizes reported at multiple time points using general linear mixed model. PLOS ONE, 11(10), e0164898.
>
> which walks you through these types of models.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Koenraad van Meerbeek
> Sent: Tuesday, 15 May, 2018 20:37
> To: Gabriele Midolo; r-sig-meta-analysis at r-project.org
> Subject: Re: [R-meta] How to conduct a meta-analysis on multiple-treatment studies with a repeated measure designs?
>
> Hi,
>
> The data is just a simplified version of the actual dataset. I included it to show the structure of the data. We also thought of only using the last sampling year, but then you´ll lose a lot of information. So I was looking for a way to include all data. So if I only have to include it in the random structure of the model, as you suggest, then I don´t see any further problems.
>
> Thanks a lot!
>
> Koenraad
>
> From: Gabriele Midolo <gabriele.midolo at gmail.com>
> Sent: dinsdag 15 mei 2018 18:42
> To: Koenraad van Meerbeek <koenraad.vanmeerbeek at bios.au.dk>; r-sig-meta-analysis at r-project.org
> Subject: Re: [R-meta] How to conduct a meta-analysis on multiple-treatment studies with a repeated measure designs?
>
> Hi Koenraad,
> I've put r-sig-meta-analysis back in the discussion too.
> Data are fine to me but I don't see the number of replicates and sd values, so I can't calculate the variance-covariance matrix myself. However, assuming you have these values, you should be able to easly calculate it with the code I've sent you?
> Note that I still think it could be worthy to account for replication due to different years as I wrote in the previous email (with random effect accounting for different year levels nested within each study). Something I saw quite often in ecological meta-analysis is to drop all "previous" observation from experiments and compare treatment and control of the last year only for each experiment reported by the study to avoid replication, but this can lead to important loss of information, I think. Of course, this also depend on the context of your study?
> Maybe some of the other people will have better inputs on this.
> With my best,
> Gabriele
>
> On 15 May 2018 at 13:34, Koenraad van Meerbeek <koenraad.vanmeerbeek at bios.au.dk<mailto:koenraad.vanmeerbeek at bios.au.dk>> wrote:
> Hi Gabriele,
>
> Thanks for the reply. I will look at the paper you mentioned.
>
> So in response to Michael´s mail, should I post my example data table like this:
> structure(list(Study = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), Treatment = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 4L, 5L, 1L, 4L, 5L), .Label = c("Control", "TR1", "TR2", "TRa", "TRb"), class = "factor"), Year = c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L), Species_div = c(1.35, 0.78, 0.23, 1.3, 0.65, 0.2, 1.74, 1.34, 1.12, 1.69, 1.21, 0.98), Magnitude = c(0, 0.75, 1.5, 0, 0.75, 1.5, 0, 0.5, 1, 0, 0.5, 1), Duration = c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L)), class = "data.frame", row.names = c(NA,
> -12L))
>
> This is the output of the dput I get.
>
> Best,
>
> Koenraad
>
> From: Gabriele Midolo <gabriele.midolo at gmail.com<mailto:gabriele.midolo at gmail.com>>
> Sent: dinsdag 15 mei 2018 13:25
> To: Koenraad van Meerbeek <koenraad.vanmeerbeek at bios.au.dk<mailto:koenraad.vanmeerbeek at bios.au.dk>>
> Cc: r-sig-meta-analysis at r-project.org<mailto:r-sig-meta-analysis at r-project.org>
> Subject: Re: [R-meta] How to conduct a meta-analysis on multiple-treatment studies with a repeated measure designs?
>
> Hi Koenraad,
> I agree with what Michael put above.
> Not sure this is what you are looking for, but I posted a similar question some time ago and got a code from Wolfgang to build variance-covariance matrix based on Lajeunesse (2011) Ecology, 92(11), pp. 2049–2055:
>
> calc.v <- function(x) {
> v <- matrix(x$SD_C[1]^2 / (x$N_C[1] * x$X_C[1]^2), nrow=nrow(x), ncol=nrow(x))
> diag(v) <- x$vi
> v
> }
>
> V <- bldiag(lapply(split(dat, dat$common_ID), calc.v)) V Where 'common_ID' is the column that codes groups of effect sizes that share the mean ( 'X_C' ), standard deviation ( 'SD_C' ),and n ('N_C') of a control group.
> When calling rma.mv<http://rma.mv>(), V is what you should then give as the second argument instead of vi.
> You should also be able to deal with repeated measure over time by adding a nested element ( e.g. "sampling year") to the multi-level structure of your model in rma.mv<http://rma.mv> (something like "random = ~ Study / year / ID "), but I could be wrong here...
>
> Cheers,
> Gabriele
More information about the R-sig-meta-analysis
mailing list