[R-meta] inner-outer-grouping-factor-structure for multiple treatment studies

Fri Sep 1 16:41:22 CEST 2017

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 <http://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 <http://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 <http://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

On Fri, Sep 1, 2017 at 3:19 AM, Michael Dewey <lists at dewey.myzen.co.uk>
wrote:

> Dear Vivien
>
> Your email got rather scrambled I fear. Can you try re-sending it after
> setting your mailer to send plain text instead of HTML to see whether that
> fixes the problem?
>
> Michael
>
> On 01/09/2017 01:48, Vivien Bonnesoeur wrote:
>
>> Dear Wolfgang
>> first of all, thank you very much for your quick and kind answer. I'm new
>> in this community and I find your work very useful. I'm sure I will manage
>> to perform this analysis with your help but for the moment I still have a
>> few difficulties.
>>
>> 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 is :
>> matrix <http://stat.ethz.ch/R-manual/R-devel/library/base/html/matr
>> ix.html>(
>> 1/x$n2i[1] + outer
>> <http://stat.ethz.ch/R-manual/R-devel/library/base/html/outer.html>(x$yi,
>> x$yi, "*")/(2*x$Ni[1]), nrow
>> <http://stat.ethz.ch/R-manual/R-devel/library/base/html/nrow.html>=nrow
>> <http://stat.ethz.ch/R-manual/R-devel/library/base/html/nrow.html>(x),
>> ncol
>> <http://stat.ethz.ch/R-manual/R-devel/library/base/html/ncol.html>=nrow
>> <http://stat.ethz.ch/R-manual/R-devel/library/base/html/nrow.html>(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 I used the
>> metamodel3 <- rma.mv(yi, ], mods = ~ Land_use_change - 1, random = ~
>> factor(trial) | article, rho=1/2, data=ma.grass) to fit the model and I
>> get
>> the error message :
>>
>>
>>
>>
>>
>>
>>
>> *Error in .ll.rma.mv <http://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 messages:1: In rma.mv <http://rma.mv>(yi, V[k,
>> k], mods = ~Luchange - 1, random = ~factor(trial) |  :  'V' appears to be
>> not positive definite.2: In rma.mv <http://rma.mv>(yi, V[k, k], mods =
>> ~Luchange - 1, random = ~factor(trial) |  :  Redundant predictors dropped
>> from the model.*
>>
>>
>> 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
>>
>>
>> On Thu, Aug 31, 2017 at 2:18 PM, Wolfgang Viechtbauer <
>> wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
>>
>> Dear Vivien,
>>>
>>> This can be thought of as a network meta-analysis, so indeed, the
>>> dat.hasselblad1998 and dat.senn2013 datasets and the corresponding code
>>> in
>>> the example section of the corresponding help files are quite relevant.
>>>
>>> The dataset as given is in the format for a 'contrast-based' analysis.
>>> And
>>> when you have studies where the same group is being reused for computing
>>> the outcomes, then this results in dependency in the sampling errors of
>>> those outcomes. You may want to study carefully:
>>>
>>> http://www.metafor-project.org/doku.php/analyses:gleser2009
>>>
>>> You are dealing with multiple treatment studies. The code given there
>>> illustrates the computation of the covariance of outcomes that share a
>>> common group. I do not cover measure="ROM" there, but the covariance is:
>>>
>>> sd^2/(n*mean^2)
>>>
>>> from the group whose data is being re-used. So, if you want to compute
>>> this programmatically, you will have to do some coding to make that work
>>> (and it's a bit more complex in your case, since sometimes the first and
>>> sometimes the second group appears to be the 'reference' group). But the
>>> idea is that you need to build the entire V matrix.
>>>
>>> And indeed, that will capture the higher dependency in such studies.
>>>
>>> Then you need to go back to help(dat.hasselblad1998) and
>>> help(dat.senn2013) to see how the modeling is done there using a
>>> contrast-based model. It is shown there, so I rather not repeat all of
>>> that
>>> here.
>>>
>>> A final note:
>>>
>>> escalc() does not have a 'method' argument, so the 'method = "REML"' part
>>> of your call to escalc() doesn't do anything.
>>>
>>> Best,
>>> Wolfgang
>>>
>>> --
>>> Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
>>> Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
>>> Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com
>>>
>>> On 08/31/2017 05:07 PM, Vivien Bonnesoeur wrote:
>>> Hi,
>>> I would need help about the way to model mutliple treatment studies.
>>> I have been posting a question on
>>> https://stats.stackexchange.com/questions/300425/metafor-rma
>>> -mv-function-inner-outer-grouping-factor-structure
>>> and was told to write to this list :
>>>
>>> I would like to know how to use the inner/outer grouping factor structure
>>> when my data set present some multiple treatment studies (e.g., when
>>> multiple treatment groups are compared with a common control/reference
>>> group, such that the data from the control/reference group is used
>>> multiple
>>> times to compute the effect sizes or outcomes)? I have been practicing
>>> the
>>> metafor examples "dat.hasselblad1998" and "dat.senn2013" where such
>>> multiple treatment studies exist but I could not understood exactly how
>>> the
>>> reference group was designated.
>>>
>>> I'm performing a meta-analysis on the impact of reforestation/forest
>>> cover
>>> change on the soil organic matter content. Here is the dataset I'm using
>>> (sep=";") :
>>>
>>>      1.
>>>    article;Land_use_change;forest_N;forest_mean;forest_sd;ref_
>>> N;ref_mean;ref_sd;inner_group_struct
>>>      2. Hes;Forestation – native forest;5;2,5;0,6;5;8,5;3,1;1
>>>      3. Hes;Forestation – native forest;5;9,1;1,2;5;22,4;3,2;2
>>>      4. Hes;Forestation – native forest;5;9,9;1,2;5;13,0;2,7;3
>>>      5. Hes;Afforestation – ungrazed grassland;5;0,6;0,2;5;0,6;0,0;4
>>>      6. Hes;Afforestation – ungrazed grassland;5;1,6;0,4;5;2,8;1,2;5
>>>      7. Henr;Afforestation – ungrazed grassland;3;6,0;0,5;3;7,0;0,4;6
>>>      8. Henr;Afforestation – grazed grassland;3;6,0;0,5;3;5,9;1,0;6
>>>      9. Henr;Afforestation – grazed grassland;3;5,9;3,2;3;5,6;4,8;7
>>>      10. Konin;Regrowth – grazed grassland;8;6,7;2,0;8;4,2;1,1;8
>>>      11. Konin;Regrowth – grazed grassland;8;9,0;3,3;8;6,4;1,7;9
>>>      12. Konin;Regrowth – grazed grassland;8;8,0;0,8;8;8,6;3,0;10
>>>      13. Konin;Regrowth – grazed grassland;8;6,3;1,6;8;5,3;0,9;11
>>>      14. Konin;Regrowth – grazed grassland;8;6,8;1,5;8;5,3;2,3;12
>>>      15. Dub;Afforestation – native forest;6;7,6;1,0;6;10,6;1,1;13
>>>      16. Dub;Afforestation – grazed grassland;6;7,6;1,0;6;13,5;0,8;13
>>>      17. Chaco;Afforestation – native forest;6;36,6;12,1;6;40,7;17,9;14
>>>      18. Chaco;Afforestation – grazed grassland;6;40,5;7,7;6;40,3;7,7;15
>>>      19. Rhoade;Forestation – native forest;10;9,2;4,1;10;11,3;2,5;16
>>>      20. Rhoade;Forestation – native forest;10;12,9;4,1;10;11,3;3,5;16
>>>      21. Rhoade;Regrowth – grazed grassland;10;9,2;4,1;10;9,3;2,7;16
>>>      22. Rhoade;Regrowth – grazed grassland;10;12,9;4,1;10;9,3;2,7;16
>>>      23. Schlatte;Forestation – native forest;10;14,3;3,1;12;12,6;2,9;17
>>>      24. Farle;Afforestation – ungrazed grassland;10;5,7;1,6;30;7,2;1,
>>> 9;18
>>>      25. Farle;Afforestation – ungrazed grassland;30;6,0;0,9;30;7,2;1,
>>> 9;18
>>>      26. Farle;Afforestation – ungrazed grassland;30;4,7;1,1;30;7,2;1,
>>> 9;18
>>>      27. Nosett;Afforestation – ungrazed grassland;5;8,6;6,9;3;8,8;7,0;
>>> 19
>>>      28. Nosett;Afforestation – grazed grassland;5;8,6;6,9;5;9,0;7,0;19
>>>      29. ManN;Afforestation – grazed grassland;24;7,5;0,8;24;5,7;1,0;20
>>>      30. Breme;Afforestation – grazed grassland;60;45,1;3,4;66;39,7;
>>> 5,4;21
>>>
>>> article refers to the study where the data are extracted, Land-use change
>>> is the fixed effect I'm interested in. In each row, I'm comparing a
>>> forestation situation with the referent situation and I've used escalc to
>>> compute the effect size= log ratio of the means
>>>
>>> ma.grass=escalc(m1 = forest_mean, m2 = ref_mean,
>>>               sd1 = forest_sd, sd2 = ref_sd,
>>>               n1 = forest_N, n2 = ref_N,
>>>               method = "REML", measure = "ROM",slab=article,
>>>               )
>>>
>>> The articles "Dub", "Rhoade", "Farle" and "Nosett" are multiple treatment
>>> studies. For example in "Farle", the the soil matter content of 3
>>> different
>>> pine plantations are compared to a single control. In the same time, an
>>> article like "Hes" gives 5 contrasts from 5 paired-site measurement.
>>> Indeed
>>> those 5 contrasts are not independent (same methodology, they come from
>>> the
>>> same region, etc...) but I believe the correlations between the contrast
>>> are less strong than between contrasts calculated from the same reference
>>> group as in the above multiple treatment studies. How can I model
>>> correctly
>>> these *a prior* different type of correlations within studies? My guess
>>> is
>>>
>>> I need to use article as the outer factor but I don't know what to write
>>> in
>>> the inner factor. I have been trying many different inner group factors
>>> but
>>> in any cases the results were satisfying. For example, like in
>>> "dat.hasselblad1998" and "dat.senn2013", I used :
>>>
>>> metamodel1=rma.mv
>>> (yi,vi,data=ma.grass,mods=~Land_use_change-1,random=~factor(
>>> id)|article,method="REML")
>>>
>>> but the results were not satisfying to me since the Land use change
>>> "Afforestation - ungrazed grassland" should be significantly <0
>>>
>>> I then created a "inner_group_struct"
>>>
>>> metamodel3=rma.mv(yi,vi,data=ma.grass,mods=~Land_use_change-
>>> 1,random=~factor(inner_group_struct)|article,struct="UN",method="REML")
>>>
>>> but it is still not satysfing : the Land.use_change "Forestation - native
>>> forest" is estimated to be >0 while almost all the raw Effect Size from
>>> this treatment apart 1 are negative.
>>>
>>> any help will be much appreciated for modeling the inner group structure.
>>> Kind regards
>>>
>>> --
>>> Vivien BONNESOEUR
>>> Docteur en biologie forestière
>>>
>>>
>>
>>
>>
> --
> Michael
> http://www.dewey.myzen.co.uk/home.html
>

-- 
Vivien BONNESOEUR
Docteur en biologie forestière

	[[alternative HTML version deleted]]