[R-meta] Two ways to calculate subgroup and overall average effect sizes - R-SIG-meta-analysis

Wed, Aug 16, 2017 9:34 AM #

Dear Naike,

Take a look at:

http://www.metafor-project.org/doku.php/tips:comp_two_independent_estimates

and especially the last section (Meta-Regression with All Studies but Different Amounts of (Residual) Heterogeneity). To get an overall estimate, but allowing for different amounts of heterogeneity across the levels of some factor (e.g., 'alloc' here), you would use:

rma.mv(yi, vi, random = ~ alloc | trial, struct="DIAG", data=dat)

Adapt to your data accordingly.

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    

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Naike Wang
Sent: Wednesday, August 16, 2017 18:29
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] Two ways to calculate subgroup and overall average effect sizes

ATTACHMENT(S) REMOVED: ATT00001.txt | do not assume same variance.jpg | assume same variance.jpg | no subgrouping.jpg | two options you can select.jpg | data.csv 

Hi all,
I noticed that in Comprehensive Meta Analysis (CMA), there are two options you can select to calculate subgroup average effect sizes. The first approach is assume a common among-study variance component across subgroups and the second one is not assume a common among-study variance component. When you select different approaches, the subgroup average effect sizes, the overall average effect size and their confidence intervals will change accordingly. This has something to do with the tau-squared calculated by different approaches. This video here explains the reason behind it in detail if you fast forward to 22:40.

I have attached my dataset here as an example. I want to do a meta-analysis of proportions so the data doesn't have a control group. It includes a moderator called studytype that divides the studies into two subgroups: birth cohort and others.

Using R, the subgroup average effect sizes I calculated are the same as obtained in CMA. My question is that how to use R to get the same overall average effect size as calculated in CMA using different approaches?

My code:
dat=read.csv("data.csv",header=T,sep=",")
es=escalc(xi=cases,ni=total,measure="PLO",data=dat)
overall.es=rma(yi,vi,data=es,method="DL",weighted=TRUE)
pred=predict(overall.es,transf=transf.ilogit,digits=5)
print(pred,digits=5)

The overall average effect size:
?pred ? ci.lb ? ci.ub ? cr.lb ? cr.ub
?0.00042 0.00032 0.00057 0.00013 0.00135?
The overall average effect size obtained above is derived when you don't do subgrouping. In CMA, if you do subgrouping, the overall effect size will change according to the approach you select. If you assume a common among-study variance, the overall average will become 0.000430 with a 95% CI [0.000307, 0.000602]; if you do not assume that, the overall average will become 0.000453 with a 95% CI [0.000335, 0.000611], which are the values I try to obtain in R (I don't know how to do this, which is why I'm asking this question).

If assuming a common among-study variance component across subgroups:
studytype=rma(yi,vi,data=ies.logit,mods=~studytype,method="DL")
pred.studytype=predict(studytype,transf=transf.ilogit,digits=5)
print(pred.studytype[c(1,17)])
then the subgroup average effect sizes are:
? ? ? pred ? ci.lb ? ci.ub ? cr.lb ? cr.ub
1 ?0.00034 0.00020 0.00061 0.00008 0.00142 #Birth cohort
17 0.00049 0.00032 0.00074 0.00012 0.00191 #Others

If not assuming a common among-study variance component across subgroups:
birthcohort=rma(yi,vi,data=ies.logit,subset=studytype=="Birth cohort",method="DL")
others=rma(yi,vi,data=ies.logit,subset=studytype=="Others",method="DL")
pred.birthcohort=predict(birthcohort,transf=transf.ilogit,digits=5)
pred.others=predict(others,transf=transf.ilogit,digits=5)
print(pred.birthcohort,digits=5)
print(pred.others,digits=5)
then the subgroup average effect sizes are:
?pred ? ci.lb ? ci.ub ? cr.lb ? cr.ub
?0.00035 0.00016 0.00078 0.00005 0.00274
?pred ? ci.lb ? ci.ub ? cr.lb ? cr.ub
?0.00047 0.00034 0.00065 0.00017 0.00132

The screen shots I attached here shows the same results I obtained above.

Please let me know if my question makes sense. It's kind of difficult to explain this very clearly in an email.

Thank you for your time!

Cheers,

Naike

Nathan Leon Pace, MD, MStat

Wed, Aug 23, 2017 3:44 PM #

Hi Wolfgang, 

I have a k = 29 SMD meta analysis.

The moderator is a three level factor.

painearly1surgery.rma <- rma(yi = yi, vi = vi, mods = ~ surgery,
                            data = painearly.df, test = 'knha', digits = 3)


painearly2surgery.rma.mv <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                         random = ~ surgery | study, struct = 'DIAG',
                         data = painearly.df, test = 't', digits = 3)

There is a nearly 10 fold variation in the individual tau^2s.

Variance Components: 

outer factor: study   (nlvls = 29)
inner factor: surgery (nlvls = 3)

           estim   sqrt  k.lvl  fixed  level
tau^2.1    0.082  0.287      8     no   open
tau^2.2    0.759  0.871     10     no    lap
tau^2.3    0.086  0.294     11     no  other


The average tau^2 is:

tau^2 (estimated amount of residual heterogeneity):     0.312 (SE = 0.110)
tau (square root of estimated tau^2 value):             0.559



The omnibus test of moderators  is not rejected in either model.

Test of Moderators (coefficient(s) 2:3):  (average tau^2)
F(df1 = 2, df2 = 26) = 2.082, p-val = 0.145

Test of Moderators (coefficient(s) 2:3):  (individual tau^2)
F(df1 = 2, df2 = 26) = 2.643, p-val = 0.090


Is there a meaningful statistical comparison of the individual tau^2s?



Are there other ways to compare the model fits (AIC or BIC)?

The anova function won?t mix rma.uni and rma.mv objects.


Nathan

Viechtbauer Wolfgang (STAT)

Thu, Aug 24, 2017 12:51 AM #

You can do a LRT to compare the models. But you have to fit the 'standard' model also with rma.mv(). See here:

http://www.metafor-project.org/doku.php/tips:rma.uni_vs_rma.mv

So, in your case:

painearly.df$id <- 1:nrow(painearly.df)

painearly1surgery.rma <- rma.mv(yi = yi, V = vi, mods = ~ surgery, random = ~ 1 | id,
                            data = painearly.df, test = 't', digits = 3)

painearly2surgery.rma.mv <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                         random = ~ surgery | study, struct = 'DIAG',
                         data = painearly.df, test = 't', digits = 3)

anova(painearly1surgery.rma, painearly2surgery.rma.mv)

So, that will give you a LRT of H0: tau^2.1 = tau^2.2 = tau^2.3.

Another way to think about this is that this is a comparison between an 'ID' and a 'DIAG' structure.

res0 <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
               random = ~ surgery | study, struct = 'ID',
               data = painearly.df, test = 't', digits = 3)

res1 <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
               random = ~ surgery | study, struct = 'DIAG',
               data = painearly.df, test = 't', digits = 3)

anova(res0, res1)

But that will give you the same results (assuming that there are no studies where the same level of 'surgery' occurs more than once).

Best,
Wolfgang

-----Original Message-----
From: Nathan Pace [mailto:n.l.pace at utah.edu] 
Sent: Thursday, August 24, 2017 00:44
To: Viechtbauer Wolfgang (SP); r-sig-meta-analysis at r-project.org
Subject: Re: [R-meta] Two ways to calculate subgroup and overall average effect sizes

Hi Wolfgang, 

I have a k = 29 SMD meta analysis.

The moderator is a three level factor.

painearly1surgery.rma <- rma(yi = yi, vi = vi, mods = ~ surgery,
                            data = painearly.df, test = 'knha', digits = 3)

painearly2surgery.rma.mv <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                         random = ~ surgery | study, struct = 'DIAG',
                         data = painearly.df, test = 't', digits = 3)

There is a nearly 10 fold variation in the individual tau^2s.

Variance Components: 

outer factor: study   (nlvls = 29)
inner factor: surgery (nlvls = 3)

           estim   sqrt  k.lvl  fixed  level
tau^2.1    0.082  0.287      8     no   open
tau^2.2    0.759  0.871     10     no    lap
tau^2.3    0.086  0.294     11     no  other

The average tau^2 is:

tau^2 (estimated amount of residual heterogeneity):     0.312 (SE = 0.110)
tau (square root of estimated tau^2 value):             0.559

The omnibus test of moderators  is not rejected in either model.

Test of Moderators (coefficient(s) 2:3):  (average tau^2)
F(df1 = 2, df2 = 26) = 2.082, p-val = 0.145

Test of Moderators (coefficient(s) 2:3):  (individual tau^2)
F(df1 = 2, df2 = 26) = 2.643, p-val = 0.090

Is there a meaningful statistical comparison of the individual tau^2s?

Are there other ways to compare the model fits (AIC or BIC)?

The anova function won?t mix rma.uni and rma.mv objects.

Nathan

Nathan Leon Pace, MD, MStat

Thu, Aug 24, 2017 4:17 AM #

Thanks for the link.

Using struct = ?ID? works very nicely.

Best,

Nathan

On 2408//2017, 1:51 AM, "Viechtbauer Wolfgang (SP)" <wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:

You can do a LRT to compare the models. But you have to fit the 'standard' model also with rma.mv(). See here:
    
    http://www.metafor-project.org/doku.php/tips:rma.uni_vs_rma.mv
    
    So, in your case:
    
    painearly.df$id <- 1:nrow(painearly.df)
    
    painearly1surgery.rma <- rma.mv(yi = yi, V = vi, mods = ~ surgery, random = ~ 1 | id,
                                data = painearly.df, test = 't', digits = 3)
    
    painearly2surgery.rma.mv <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                             random = ~ surgery | study, struct = 'DIAG',
                             data = painearly.df, test = 't', digits = 3)
    
    anova(painearly1surgery.rma, painearly2surgery.rma.mv)
    
    So, that will give you a LRT of H0: tau^2.1 = tau^2.2 = tau^2.3.
    
    Another way to think about this is that this is a comparison between an 'ID' and a 'DIAG' structure.
    
    res0 <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                   random = ~ surgery | study, struct = 'ID',
                   data = painearly.df, test = 't', digits = 3)
    
    res1 <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                   random = ~ surgery | study, struct = 'DIAG',
                   data = painearly.df, test = 't', digits = 3)
    
    anova(res0, res1)
    
    But that will give you the same results (assuming that there are no studies where the same level of 'surgery' occurs more than once).
    
    Best,
    Wolfgang
    
    -----Original Message-----
    From: Nathan Pace [mailto:n.l.pace at utah.edu] 
    Sent: Thursday, August 24, 2017 00:44
    To: Viechtbauer Wolfgang (SP); r-sig-meta-analysis at r-project.org
    Subject: Re: [R-meta] Two ways to calculate subgroup and overall average effect sizes
    
    Hi Wolfgang, 
    
    I have a k = 29 SMD meta analysis.
    
    The moderator is a three level factor.
    
    painearly1surgery.rma <- rma(yi = yi, vi = vi, mods = ~ surgery,
                                data = painearly.df, test = 'knha', digits = 3)
    
    painearly2surgery.rma.mv <- rma.mv(yi = yi, V = vi, mods = ~ surgery,
                             random = ~ surgery | study, struct = 'DIAG',
                             data = painearly.df, test = 't', digits = 3)
    
    There is a nearly 10 fold variation in the individual tau^2s.
    
    Variance Components: 
    
    outer factor: study   (nlvls = 29)
    inner factor: surgery (nlvls = 3)
    
               estim   sqrt  k.lvl  fixed  level
    tau^2.1    0.082  0.287      8     no   open
    tau^2.2    0.759  0.871     10     no    lap
    tau^2.3    0.086  0.294     11     no  other
    
    The average tau^2 is:
    
    tau^2 (estimated amount of residual heterogeneity):     0.312 (SE = 0.110)
    tau (square root of estimated tau^2 value):             0.559
    
    The omnibus test of moderators  is not rejected in either model.
    
    Test of Moderators (coefficient(s) 2:3):  (average tau^2)
    F(df1 = 2, df2 = 26) = 2.082, p-val = 0.145
    
    Test of Moderators (coefficient(s) 2:3):  (individual tau^2)
    F(df1 = 2, df2 = 26) = 2.643, p-val = 0.090
    
    Is there a meaningful statistical comparison of the individual tau^2s?
    
    Are there other ways to compare the model fits (AIC or BIC)?
    
    The anova function won?t mix rma.uni and rma.mv objects.
    
    Nathan