[R-meta] Differences between metafor / robumeta / lme4
James Pustejovsky
jepu@to @end|ng |rom gm@||@com
Sat Jul 11 15:23:49 CEST 2020
Bernard,
Could you clarify what is the goal of your analysis? What question(s) are you trying to answer? And in particular, are you aiming to draw a random effects inference (to a broader population of studies) or a fixed effects inference (limited to the set of studies in the analysis)?
Generally, I would recommend using an approach that captures the correlations between the effect size estimates (i.e., the correlation between the estimated correlation coefficients). As far as I can tell, none of the approaches you listed really does that. The robumeta and metafor models make arbitrary assumptions about the correlation, then use RVE to protect against this. But given that you have the raw data, you should be able to obtain all the information you would need to understand the correlations—either through the old Steiger 1980 formulas or (probably better?) by bootstrapping the set of ES estimates from each study.
James
> On Jul 11, 2020, at 6:15 AM, Bernard Fernou <bernard.fernou using gmail.com> wrote:
>
> Hi everyone,
> I hope this email finds you well.
>
> I am trying to model the relationship between depression and BMI over
> several studies conducted in our research lab. We have access to all
> individual patient data.
>
> We have 7 seven studies in which BMI is always assessed using a
> digital scale and a meter but in which we used various measures to assess
> depression (self reported questionnaire / interviews / medical file). Each
> study included either 1, 2 or the 3 measures of depression. The effect
> size used is correlation.
>
> We have analyzed our data using 3 methods (2 two-stage approaches: robu
> from robumeta / 3 level meta analysis with robust SE using metafor; and 1
> one-stage approach: a mixed model using lme4).
>
> A complete reproducible example is provided below.
>
> Results from metafor and lme4 agreed (p values were not significant and
> very similar: p = .13, .14). However, results from robumeta are quite
> different (the p value is barely significant p = .03). I tried to use
> several rho values in robumeta but it did not change anything. Note that
> heterogeneity/inconsistency (assessed using tau² / I²) are very low.
>
> Unfortunately, as this meta analysis did not followed any systematic
> review, we made the mistake to not pre register it and to directly analyze
> the data. Therefore, we have not any proof of what analysis what thought as
> primary. Even if the effet sizes are similar across 3 methods, it is
> probable that readers will be sensitive to p-values. The choice of the
> statistical analysis is thus critical.
>
> I know that one-stage and two stages approaches could differ because of the
> model specification (Burke et al., 2017). Thus; I suspect that we missed a
> specification differing between robumeta on the one side and both
> metafor+lme4 on the other.
>
> I do not know if this mailing list is the right place to ask for help but
> if someone could give me a feedback on the code used and on the best
> approach to analyze our data, it would be greatly appreciated!
>
> Best,
>
> Bernard
>
>
> set.seed(1234)
> BMI<-rnorm(1000, 23, 4)
> Depression<-rnorm(1000, 10, 2)+0.1*BMI
> cor.test(~BMI+Depression)
> Study<-rep(1:6, times=c(100,100,100,300,200,200))
> Outcome<-c(rep(c("Questionnaire", "Interview"), times=c(100,200)),#first 3
> studies
> rep(c("Questionnaire", "Interview", "MedicalFile"),
> each=100),#4th study
> rep(c("Questionnaire", "Interview"), each=100, times=2))#5&6th
> studies
>
> df<-data.frame(cbind(BMI,Depression,Study,Outcome))
>
> #see the structure of the data
> df %>%
> group_by(Study ,Outcome) %>%
> summarise(N=n())
>
> corFUNC=function(x){
> list(
> cor.test(~as.numeric(x$BMI)+as.numeric(x$Depression))$estimate,
> (cor.test(~as.numeric(x$BMI)+as.numeric(x$Depression))$parameter+2))
> }
>
> df.cor<-split(df, list(df$Study, df$Outcome), drop=TRUE)
> lapply(df.cor,corFUNC)
>
> df.meta.1<-as.data.frame(do.call(rbind, lapply(df.cor,corFUNC)))
>
> df.meta<-df.meta.1 %>%
> dplyr::mutate(Study.Outcome=row.names(df.meta.1)) %>%
> tidyr::separate(Study.Outcome, into=c("Study", "Outcome")) %>%
> dplyr::rename("Cor"=V1,"N"=V2) %>%
> dplyr::arrange(Study) %>%
> dplyr::mutate(
> yi=yi<-escalc(measure="ZCOR", ri=as.numeric(Cor), ni=as.numeric(N),
> data=df.meta)$yi,
> vi=yi<-escalc(measure="ZCOR", ri=as.numeric(Cor), ni=as.numeric(N),
> data=df.meta)$vi)
>
> Approach1=robumeta::robu(
> formula = yi ~ 1,
> v=vi,
> studynum = Study,
> small = TRUE,
> rho=0.3,
> data = df.meta)
>
> df$BMI<-as.numeric(df$BMI); df$Depression<-as.numeric(df$Depression)
>
> Approach2=lme4::lmer(BMI~Depression+(1+Depression|Study)+(1+Depression|Outcome:Study),
> data=df)
>
> Approach3=rma.mv(yi, vi, random = ~ 1 |Study/Outcome, data=df.meta)
> robust(Approach3, cluster = df.meta$Study)
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
More information about the R-sig-meta-analysis
mailing list