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