[R-meta] Violation in non-independece of errors (head to head studies and mutlilevel meta-analysis)?
Gerta Ruecker
ruecker at imbi.uni-freiburg.de
Wed Apr 25 13:27:01 CEST 2018
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_language=en
More information about the R-sig-meta-analysis
mailing list