[R-meta] inner-outer-grouping-factor-structure for multiple treatment studies
Vivien Bonnesoeur
bonnesoeur.vivien at gmail.com
Tue Sep 5 00:18:07 CEST 2017
thanks for your help
I've been computing the variance-covariance as you indicated :
ma.grass <- escalc(m1 = forest_mean, m2 = ref_mean, sd1 = forest_sd, sd2 =
ref_sd,
n1 = forest_N, n2 = ref_N, measure = "ROM", slab =
article)
V <- diag(ma.grass$vi)
and case by case, i filled the covariance with
V[row,column] <- V[column,row] <- forest_sd[row]^2 / (forest_N[row] *
forest_mean[row]^2)
or
V[row,column] <- V[column,row] <- ref_sd[row]^2 / (ref_N[row] *
ref_mean[row]^2)
I was able to fit a model but I'm not convinced about the result for
several reasons :
-- it does not confirm my intuition for the comparison between
afforestation and ungrazed grassland which did not show a significant
result while most of the contrast where significantly <0
- the weights of the different contrasts also looks strange. For example
Breme's study has only 1 contrast and have a total number of sample = 126,
much more than the other studies but its weight is only 2.81% . On the
contrary, studies with several contrast can present individual contrast
weights higher than 3% (and the total weight of the study higher than 15%)
while the their total number of sample is twice as less as Breme. If I
adjust a fixed-effect model, Breme indeed has a much stronger weight
(41%). A fixed-effect model is likely not relevant in my case but the
difference of weight with the random effect model is so strong that I have
the overall impression that studies with several contrasts are overweighted
in the random effect model.
Actually, in the code you proposed for SMD effect size, the
variance-covariance matrix is calculated from the total number of sample
from the study (called "Ni" in the example), so that each study with
multiple treatment are well weighted. In the code I used for the ROM effect
size, there was no such adjustment to account for the real number of
samples for studies with several contrasts.
Is it possible to do such adjustment to get more reasonable weights?
I hope my questions and remarks are not too irrelevant, it's the first time
I'm doing a meta-analysis.
best regards
Vivien
On Mon, Sep 4, 2017 at 9:25 AM, Viechtbauer Wolfgang (SP) <
wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
> I am not quite sure why you switched to the standardized mean difference
> as the outcome measure if your original plan was to use the (log) ratio of
> means. Switching just because I happen to give some illustrative code for
> the SMD (instead of ROM) on the metafor website is not really a good reason
> for such a decision IMHO.
>
> In the end, this is purely a programming issue now. As I mentioned, your
> dataset is a bit more complex, since in your case, it isn't always the
> first (or second) group whose data is being reused, but sometimes the
> first, sometimes the second. So the code that I give on the metafor website
> isn't going to work for you anyway without some further adjustments.
>
> Again, for ROM, the covariance is simply:
>
> sd^2/(n*mean^2)
>
> from the group whose data is being re-used. escalc() already gives you the
> variances, which go on the diagonal of the V matrix. So, you just need to
> find a way to fill in the off-diagonal elements of V using the formula
> above. You could even get this done with something like:
>
> ma.grass <- escalc(m1 = forest_mean, m2 = ref_mean, sd1 = forest_sd, sd2 =
> ref_sd,
> n1 = forest_N, n2 = ref_N, measure = "ROM", slab =
> article)
>
> V <- diag(ma.grass$vi)
>
> And let's say the first and second row are data from a multi-treatment
> study with 2 treatment groups and 1 control group, where 'ref' corresponds
> to the control group (so ref_mean/sd/N[1] and ref_mean/sd/N[2] are actually
> the same data). Then:
>
> V[1,2] <- V[2,1] <- ref_sd[1]^2 / (ref_N[1] * ref_mean[1]^2)
>
> will fill in the covariance of the first and second row into the V matrix.
>
> Now proceed in the same way for other studies. If the data from the
> 'forest' group are being reused, then use:
>
> V[row,column] <- V[column,row] <- forest_sd[row]^2 / (forest_N[row] *
> forest_mean[row]^2)
>
> This isn't elegant and 'hard coded' but it gets the job done.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: Vivien Bonnesoeur [mailto:bonnesoeur.vivien at gmail.com]
> Sent: Friday, 01 September, 2017 16:41
> To: Michael Dewey
> Cc: Viechtbauer Wolfgang (SP); r-sig-meta-analysis at r-project.org
> Subject: Re: [R-meta] inner-outer-grouping-factor-structure for multiple
> treatment studies
>
> Oh yes, It looked clear when I sent it but it is unreadable now. It should
> be clearer this time
>
> After Wolfgang advices, I've been changing the way I'm computing the
> effect size and I have been using the standardised mean difference instead
> of the log ratio of the mean
> to be able to use the code from "gleser 2009". In this code, the
> covariance in the covariance matrix is calculated as:
>
> matrix(1/x$ref_mean[1] + outer(x$yi, x$yi, "*")/(2*x$Ntot[1]),
> nrow=nrow(x), ncol=nrow(x))
>
> where n2i is the control group reused. I used this code because N was equal
> for the treatment and reference group. There is only the study from Rhoade
> where I'm not sure of the application of this formula because in this case
> the treatment and the reference are reused once.
> I was able to get the variance-covariance V matrix and fit the model:
>
> metamodel3 <-rma.mv (yi, V, mods = ~ Land_use_change - 1, random = ~
> factor(id) | article, rho=1/2, data=ma.grass[k,])
> but it did not work and I get the error message :
>
> Error in .ll.rma.mv(opt.res$par, reml = reml, Y = Y, M = V, A = A, X.fit
> = X, :
> Final variance-covariance matrix not positive definite.
> De plus : Warning message:
> In rma.mv(yi, V[k, k], mods = ~Luchange - 1, random = ~factor(trial) | :
> 'V' appears to be not positive definite.
>
> after some tests restricting the matrix to individual studies, it appears
> that the sub-matrix from the study Farle is not positive definite.
> Could you give a little more help to solve this problem and to confirm
> that the code I used was appropriate?
> I would like not dropping this study which is quite important in the
> meta-analysis.
>
> Thanks in advance
> Vivien
>
--
Vivien BONNESOEUR
Docteur en biologie forestière
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list