[R-meta] Violation in non-independece of errors (head to head studies and mutlilevel meta-analysis)?
Natan Gosmann
natan.gosmann at gmail.com
Thu Apr 26 03:15:38 CEST 2018
Dear Gerta,
Thank you very much for answering.
Our main purpose in this project is to perform multilevel multiple
meta-regression models. Would be possible to perform multilevel models
adjusting for multiple moderators and also network meta-analysis?
The idea of the example data was more to explain our current data
structure, not so much its content, thats why I didn't include sample size
and covariates in the example.
Here it goes the example with sample size and two covariates as an example:
study = c("1", "1", "2","2", "3", "3", "3","3")
trt = c( "1", "2", "1","1", "1", "1", "2","2" )
drug= c( "flx", "esc", "cit","cit", "dul", " dul", "flx","flx" )
outcome = c( "HAM", "HAM", "SDS","HAM", "MADRS", "HAM", "MADRS","HAM" )
pcb_n= c( "252", "252", "127","127", "184", "184", "184","184" )
pcb_m1= c( "25", "25", "13","24", "31", "20", "31","20" )
pcb_sd1 = c( "3", "3", "2","2.5", "6", "5", "6","5" )
pcb_m2 = c( "19", "19", "9","20", "24", "17", "24","17" )
drug_n= c( "249", "255", "130","130", "180", "180", "187","187" )
drug_m1= c( "26", "24", "12","25", "29", "21", "30","20" )
drug_sd1 = c( "3", "2.5", "1.5","2", "6", "4.5", "5.5","4" )
drug_m2 = c( "12", "13", "3","14", "17", "14", "16","13" )
sampling= c( "outpatients", "outpatients", "inpatients","inpatients",
"outpatients", "outpatients", "outpatients","outpatients" )
publication_year= c( "1995", "1995", "2013","2013", "2016", "2016",
"2016","2016" )
data = data.frame( study, trt, drug, outcome, pcb_n, pcb_m1, pcb_sd1,
pcb_m2, drug_n, drug_m1, drug_sd1, drug_m2, samping, publication_year)
Thank you very much for your help.
Best,
Natan
2018-04-25 8:27 GMT-03:00 Gerta Ruecker <ruecker at imbi.uni-freiburg.de>:
> Dear Natan,
>
> notwithstanding your problem with network meta-analysis (which I do not
> understand):
>
> I do not see any standard errors or sample sizes in your example data (and
> also no covariates). Could you provide the sample sizes?
>
> Best,
>
> Gerta
>
>
>
> Am 24.04.2018 um 22:17 schrieb Natan Gosmann:
>
>> Hello all,
>>
>> Here goes a better example of our data structure (maybe the previous one
>> was not so clear because of the text format):
>>
>> study = c("1", "1", "2","2", "3", "3", "3","3")
>> trt = c( "1", "2", "1","1", "1", "1", "2","2" )
>> drug= c( "flx", "esc", "cit","cit", "dul", " dul", "flx","flx" )
>> outcome = c( "HAM", "HAM", "SDS","HAM", "MADRS", "HAM", "MADRS","HAM" )
>> pcb_m1= c( "25", "25", "13","24", "31", "20", "31","20" )
>> pcb_sd1 = c( "3", "3", "2","2.5", "6", "5", "6","5" )
>> pcb_m2 = c( "19", "19", "9","20", "24", "17", "24","17" )
>> drug_m1= c( "26", "24", "12","25", "29", "21", "30","20" )
>> drug_sd1 = c( "3", "2.5", "1.5","2", "6", "4.5", "5.5","4" )
>> drug_m2 = c( "12", "13", "3","14", "17", "14", "16","13" )
>>
>> data = data.frame( study, trt, drug, outcome, pcb_m1, pcb_sd1, pcb_m2,
>> drug_m1, drug_sd1, drug_m2)
>>
>>
>> Any advice or suggestions on the questions described in my previous post
>> would be greatly appreciated and we would be really thankful if someone
>> could clarify these points for us.
>>
>> Best,
>> Natan
>>
>> 2018-04-03 21:25 GMT-03:00 Natan Gosmann <natan.gosmann at gmail.com>:
>>
>> Dear Wolfgang,
>>>
>>>
>>>
>>> Thank you very much for your response and all comments. Here it goes an
>>> example of what we are doing.
>>>
>>>
>>>
>>> Our full data frame is composed of 1100 rows e 700 columns, but an
>>> example
>>> of the data structure could be described like this (basically each study
>>> has one or more trt group, identified by drug, compared to placebo
>>> considering multiple outcomes, identified by measurement scales):
>>>
>>>
>>>
>>> study /trt/drug/ outcome/pcb _m1/pcb _sd1/pcb_m2/drug _m1/drug
>>> _sd1/drug_m2
>>> ------------------------------------------------------------
>>> -------------------------------------------------------
>>> 1 1 flx HAM 25 3
>>> 19 26 3 12
>>> 1 2 esc HAM 25 3
>>> 19 24 2.5 13
>>> 2 1 cit SDS 13 2
>>> 9 12 1.5 3
>>> 2 1 cit HAM 24 2.5
>>> 20 25 2 14
>>> 3 1 dul MADRS 31 6
>>> 24
>>> 29 6 17
>>> 3 1 dul HAM 20 5
>>> 17 21 4.5 14
>>> 3 2 flx MADRS 31 6
>>> 24
>>> 30 5.5 16
>>> 3 2 flx HAM 20 5
>>> 17 20 4 13
>>>
>>>
>>>
>>>
>>>
>>> (1) We are calculating Differences in Standardized Mean Change from Pre
>>> and Post Means and initial SD according to the Gleser & Olkin 2009 Study;
>>>
>>>
>>>
>>> # assuming a 0.25 correlation between pre and post means (ri)
>>>
>>> meta_pcb <- escalc(measure="SMCR", m1i= pcb_m2, m2i= pcb _m1, sd1i= pcb
>>> _sd1, ni=n_pcb, ri=r1,data=mydata2)
>>>
>>>
>>>
>>> meta_drug<- escalc(measure="SMCR", m1i= drug _m2, m2i= drug _m1,
>>> sd1i=drug_sd1, ni=n_drug, ri=r1,data=mydata2)
>>>
>>>
>>>
>>> meta <- data.frame(yi = meta_drug$yi - meta_pcb$yi, vi = meta_drug$vi +
>>> meta_pcb$vi)
>>>
>>>
>>>
>>> (2) Considering that we also want to include studies with more than one
>>> treatment group, we constructed a full V matrix as you suggested;
>>>
>>>
>>>
>>> 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(meta, meta$study), calc.v))
>>>
>>>
>>>
>>>
>>>
>>> (3) So far we were specifying as random variables: Study ID and
>>> Measurement Scale (outcome) as exemplified;
>>>
>>>
>>>
>>> meta_out1<-rma.mv(yi=yi, V=V, data=meta,
>>>
>>> random=list(~1|outcome, ~1|study),
>>>
>>> slab=paste(author, outcome, sep=", "),
>>>
>>> mods = ~relevel(factor(drug), ref="flx"))
>>>
>>>
>>>
>>>
>>>
>>> However, we have doubts if our current analysis is being performed
>>> correctly. Considering our current data structure (as exemplified above),
>>> isn’t problematic to construct a full V matrix to compute the covariance
>>> for various effect sizes of different treatment groups of the same study,
>>> since we also have different rows for different outcomes of the same
>>> treatment group? Should we also include trt as a random variable?
>>>
>>>
>>>
>>> Any advice or suggestions on that would be greatly appreciated.
>>>
>>>
>>>
>>> Best,
>>>
>>>
>>>
>>> Natan
>>>
>>>
>>> 2018-03-31 12:52 GMT-03:00 Viechtbauer Wolfgang (SP) <
>>> wolfgang.viechtbauer at maastrichtuniversity.nl>:
>>>
>>> Convergence problems are difficult to anticipate. But in general, they
>>>> are more likely to appear when fitting complex models and/or when the
>>>> dataset is small. It is also possible that one is trying to fit an
>>>> overparameterized model, that is, a model where certain parameters are
>>>> not
>>>> identifiable. The issue of identifiability is complex, but some articles
>>>> that deal with this are:
>>>>
>>>> Kreutz, C., Raue, A., Kaschek, D., & Timmer, J. (2013). Profile
>>>> likelihood in systems biology. The FEBS Journal, 280(11), 2564-2571.
>>>>
>>>> Lavielle, M., & Aarons, L. (2016). What do we mean by identifiability in
>>>> mixed effects models? Journal of Pharmacokinetics and Pharmacodynamics,
>>>> 43(1), 111-122.
>>>>
>>>> Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M.,
>>>> Klingmuller, U., & Timmer, J. (2009). Structural and practical
>>>> identifiability analysis of partially observed dynamical models by
>>>> exploiting the profile likelihood. Bioinformatics, 25(15), 1923-1929.
>>>>
>>>> Raue, A., Kreutz, C., Maiwald, T., Klingmuller, U., & Timmer, J. (2011).
>>>> Addressing parameter identifiability by model-based experimentation. IET
>>>> Systems Biolology, 5(2), 120-130.
>>>>
>>>> Wang, W. (2013). Identifiability of linear mixed effects models.
>>>> Electronic Journal of Statistics, 7, 244-263.
>>>>
>>>> One way of assessing parameter identifiability is to examine/plot
>>>> profile
>>>> likelihoods. This is what the profile() function is for. When fitting
>>>> complex models, I would always recommend to profile all
>>>> variance/correlation components.
>>>>
>>>> Even if all components are identifiable, it may be difficult to find the
>>>> ML/REML estimates. Complex models require optimization over a large
>>>> number
>>>> of parameters and this is not a trivial task. rma.mv() uses nlminb() by
>>>> default, but that is not always the best option. One can try many other
>>>> optimizers (using the control arguments 'optimizer' and 'optmethod').
>>>> See
>>>> the 'Note' section under help(rma.mv).
>>>>
>>>> As for power, there are these two articles:
>>>>
>>>> Hedges, L. V., & Pigott, T. D. (2001). The power of statistical tests in
>>>> meta-analysis. Psychological Methods, 6(3), 203-217.
>>>> Hedges, L. V., & Pigott, T. D. (2004). The power of statistical tests
>>>> for
>>>> moderators in meta-analysis. Psychological Methods, 9(4), 426-445.
>>>>
>>>> But they do not deal with complex models (just standard
>>>> random/mixed-effects models). Indeed, for complex models, one would
>>>> need to
>>>> take a simulation approach.
>>>>
>>>> Best,
>>>> Wolfgang
>>>>
>>>> -----Original Message-----
>>>> From: Emily Finne [mailto:emily.finne at uni-bielefeld.de]
>>>> Sent: Wednesday, 28 March, 2018 10:35
>>>> To: Viechtbauer Wolfgang (SP); r-sig-meta-analysis at r-project.org
>>>> Subject: Re: [R-meta] Violation in non-independece of errors (head to
>>>> head studies and mutlilevel meta-analysis)?
>>>>
>>>> I really wished, everything WOULD be so obvious to me... !
>>>>
>>>> For this analysis the results turned out to be nearly unchanged when
>>>> including the crossed random effects, although you guessed right that
>>>> convergence problems could emerge. These were related to those
>>>> parameters
>>>> estimating the covariaces between effects within studies.
>>>>
>>>> I wonder how one can anticipate such problems in advance or rather
>>>> determine how complex a model can be with given data to have enough
>>>> power
>>>> to test (fixed) moderator effects of interest and to make sure that
>>>> confidence intervals are reliable.
>>>> Is there something like a formal power analysis for meta-analysis or
>>>> meta-regression? I am aware that this is complex and I think in mixed
>>>> effects models, in general, one would use simulations.
>>>>
>>>> Any advice on literature I could read to get prepared for further
>>>> projects?
>>>>
>>>> Best,
>>>> Emily
>>>>
>>>> Am 25.03.2018 um 12:16 schrieb Viechtbauer Wolfgang (SP):
>>>> See comments below.
>>>>
>>>> Best,
>>>> Wolfgang
>>>>
>>>> -----Original Message-----
>>>> From: Emily Finne [mailto:emily.finne at uni-bielefeld.de]
>>>> Sent: Saturday, 24 March, 2018 21:41
>>>> To: Viechtbauer Wolfgang (SP); r-sig-meta-analysis at r-project.org
>>>> Subject: Re: [R-meta] Violation in non-independece of errors (head to
>>>> head studies and mutlilevel meta-analysis)?
>>>>
>>>> Dear Wolfgang,
>>>>
>>>> oh yes, many people were sick during the last weeks, here too. I hope
>>>> you're feeling better by now.
>>>>
>>>> Thanks - not 100% but good enough to catch up on work.
>>>>
>>>> Yes, that is exactly the data structure I have. I completely missed out
>>>> to think of the problem as crossed random effects! Thank you!
>>>>
>>>> I am quite sure that I constructed the var-cov-matrix V right. I used
>>>> the
>>>> formulas by Gleser & Olkin and James Pustejovsky. I double-checked the
>>>> resulting matrix. Additionally, I used robust estimation, since most
>>>> correlations between outcomes were only a best guess.
>>>>
>>>> Great!
>>>>
>>>> Only to make sure, that I understand correctly the point about the
>>>> random
>>>> effects: I code two different treatment groups within one study with
>>>> different numbers starting with 1 (for example) and than use the code
>>>> you
>>>> provided for the crossed random effects. But the numbers given to
>>>> different treatments are arbitrary and don't mean that the group with
>>>> 'treatment = 1' always got the same treatment. It is only to code that
>>>> treatment 1 and 2 within one study are different (say medication A and B
>>>> each compared against a placebo control), not that 1 and 2 always means
>>>> the same thing in different studies (it could also stand for medication
>>>> B
>>>> and C vs. control in another study). Am I right?
>>>>
>>>> Correct!
>>>>
>>>> On the other hand, if treatment 1 always stood for medication A and
>>>> treatment 2 always stood for medication B (across all studies), then it
>>>> would make sense to distiniguish the two and, for example, allow a
>>>> different tau^2 for treatment 1 vs treatment 2. I assumed that for
>>>> 'outcome', this is actually the the case (so, for example, outcome 1
>>>> always
>>>> stands for measure X and outcome 2 always stands for measure Y).
>>>>
>>>> One thing I forgot to mention: At the moment, the 'inner' terms in a '~
>>>> inner | outer' formula under random must be a character/factor variable.
>>>> So, if 'outcome' and 'trt' are numeric, then use:
>>>>
>>>> random = list(~ factor(outcome) | study, ~ factor(trt) | study),
>>>> struct=c("UN","CS")
>>>>
>>>> Or you could code 'outcome' and 'trt' as character variables to begin
>>>> with.
>>>>
>>>> The treatments we look at in our analysis, in fact, are all different in
>>>> some aspects although pursuing the same goal. We use characteristics of
>>>> the treatments as moderators then and hope to explain differences in
>>>> effect sizes.
>>>>
>>>> Again, spot on. Using this example, it also illustrates why it would
>>>> make
>>>> sense to include a fixed effect for 'outcome' (since outcome 1 and 2 are
>>>> uniquely defined), while it would not make sense to include a fixed
>>>> effect
>>>> for 'trt'. For the latter, as you say, we can use characteristics of the
>>>> treatments as moderators.
>>>>
>>>> Sorry if this is all obvious to you, but having this written down here
>>>> is
>>>> useful for future reference. Also, to come back to the Konstantopoulos
>>>> (2011) and Berkey et al. (1998) examples, this also explains why we
>>>> would
>>>> use:
>>>>
>>>> rma.mv(yi, vi, random = ~ factor(study) | district, struct="CS",
>>>> data=dat)
>>>>
>>>> in the Konstantopoulos example (https://tinyurl.com/ybpzn5ra) (so no
>>>> fixed effect for 'study' and struct="CS" -- which is actually the
>>>> default,
>>>> but added here for clarity)
>>>>
>>>> and
>>>>
>>>> rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial,
>>>> struct="UN", data=dat)
>>>>
>>>> in the Berkey example (https://tinyurl.com/y9yv366v) (so with fixed
>>>> effect for 'outcome' and struct="UN").
>>>>
>>>> In the first case, the coding of 'study' within 'district' is arbitary.
>>>> In the second case, the coding of 'outcome' within 'trial' is
>>>> meaningful.
>>>>
>>>> Again, thank you so much for you detailed help.
>>>>
>>>> I will try, if the model with the crossed effects converge. Otherwise, I
>>>> would stick to the old model (only random = ~ outcome | study,
>>>> struct="UN") and discuss this as a limitation.
>>>>
>>>> Best,
>>>> Emily
>>>> _______________________________________________
>>>> R-sig-meta-analysis mailing list
>>>> R-sig-meta-analysis at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>>>>
>>>>
>>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-meta-analysis mailing list
>> R-sig-meta-analysis at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>>
>
> --
>
> Dr. rer. nat. Gerta Rücker, Dipl.-Math.
>
> Institute of Medical Biometry and Statistics,
> Faculty of Medicine and Medical Center - University of Freiburg
>
> Stefan-Meier-Str. 26, D-79104 Freiburg, Germany
>
> Phone: +49/761/203-6673
> Fax: +49/761/203-6680
> Mail: ruecker at imbi.uni-freiburg.de
> Homepage: https://portal.uni-freiburg.de/imbi/persons/ruecker?set_lang
> uage=en
>
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list