1
$\begingroup$

I am doing a meta analysis using R. After running an overall random effect model I wanted to look at subgroups in a mixed-effect model and was using the meta and dmetar package (Doing Meta Analysis in R), using the metagen()function which worked fine. However, I might have issues with 1) non-independend data and 2) also correlated sampling errors in some cases. Thus, I thought about switching to the metafor package to use 1) the rma.mv() function to account for the non-independet data and using a 3-level model and 2) use the robust () function for the cases where I have correlated sampling errors.

In doing so I compared the metagen with the rma model. I could get similar results for the basic (reduced) random effects model. But when I run the mixed-effect subgroup analysis I get very different results. I'm trying to figure out why because using one package indicates there is a significant heterogeneity between the groups whereas the other doesn’t. Does anyone have any experience with these packages to know why?

Then I ran a random effects model with metagen () ...

MA1 <- metagen(TE, seTE, data = mydata, studlab = paste(ID),comb.fixed = FALSE,comb.random = TRUE,
             method.tau = "REML",  hakn = FALSE, sm = "SMD")

and rma ()

MA2 <- rma(TE, var, data = mydata, method="REML")

This yields the same results.

However if I do this:

MA1.taxa <- metagen(TE, seTE, data = mydata, studlab = paste(ID), comb.fixed = FALSE, comb.random = TRUE, method.tau = "REML",  hakn = FALSE, sm = "SMD", byvar = Taxa)

and the same using rma:

MA2.taxa <- rma(TE, var, mod= ~Taxa, data = mydata, method = "REML")

I get very different results. It would be great if someone could help me out here. I have the feeling it might have something to do with some of the arguments which I didn’t specify correctly. But in looking in the packages I couldn’t figure it out.

This is my data frame:

**"ID"**|**"Study"**|**"Organism"**|**"Taxa"**|**"Species"**|**"Ne"**|**"Me"**|**"Se"**|**"Nc"**|**"Mc"**|**"Sc"**|**"TE"**|**"seTE"**|**"var"**
:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:
1|1|1|"A "|"F"|10|259.5|41.1096|10|405.5|15.8114|-4.48970286181781|0.86566827459304|0.74938156163689
2|1|1|"A "|"F"|10|285.9|12.6491|10|452.5|15.8114|-11.1441917779727|1.89336302866627|3.58482355832032
3|1|1|"A "|"F"|10|301.5|12.6491|10|459.1|12.6491|-11.9329315975384|2.02012537466928|4.08090652938268
4|2|2|"B "|"N"|8|40.58|8|8|565.3|8|-62.0123636363636|11.6055591474517|134.689003125
5|2|3|"B "|"Q"|8|15.59|1|8|436.2|5|-110.293157961666|20.6281624673369|425.521086778846
6|2|4|"B "|"R"|8|40.58|8|8|565.3|8|-62.0123636363636|11.6055591474517|134.689003125
7|3|5|"B "|"C"|4|1.887|1.9|4|519.3|26|-24.40757392403|7.05271437893494|49.7407801108356
8|3|6|"B "|"J"|4|18.4|1.5|4|519.3|24|-25.6159441945354|7.39845232761053|54.7370968439256
9|4|7|"B "|"K"|6|0.8174|2.4495|6|296.3|3.1843|-96.0142079773492|21.2399012943102|451.133406992038
10|5|8|"B "|"G"|12|20.2|2.4744|12|30.05|5.8324|-2.12287801050705|0.517088404453606|0.267380418020376
11|6|9|"C"|"A"|4|317.4|44|4|509.2|20|-4.88010696352562|1.5711445799467|2.46849529109589
12|6|10|"C"|"B"|4|333.3|20|4|541.8|16|-10.0108579072198|2.96371122455629|8.78358422256097
13|7|11|"B "|"S"|5|93|8.9443|5|477|11.1803|-34.2583184199475|8.50469827426505|72.329892736287
14|8|12|"B "|"H"|8|175.4|19.799|8|330.4|19.799|-7.4016594042858|1.4714809993973|2.16525633158727
15|9|13|"C"|"E"|6|353|7|6|584|14|-19.2655685142566|4.2992247363139|18.4833333333333
16|10|14|"C"|"O"|5|307.3|8|5|546.798|16|-17.1016587915647|4.28074511636117|18.32477875125
17|11|15|"B "|"D"|12|85.77|9|12|954.8|9|-93.2292720306513|13.9430530966266|194.40872965535
18|11|15|"B "|"D"|12|105.5|9|12|943.5|9|-89.9003831417625|13.4456304955873|180.784979423868
19|11|16|"B "|"I"|12|519.1|9|12|834.2|9|-33.8038314176245|5.06988227434525|25.7037062757202
20|11|16|"B "|"I"|12|504.6|9|12|861.6|9|-38.2988505747126|5.73992676427664|32.9467592592593
21|12|17|"B "|"M"|10|60|12.6491|6|187|31.8434|-5.56823589086092|1.16215466879402|1.35060347419974
22|12|17|"B "|"M"|10|74.6|12.6491|6|227|31.8434|-6.6818830690331|1.35186377624163|1.82753566951429
23|12|17|"B "|"M"|10|92.5|18.9737|12|276|13.8564|-10.7914331134072|1.74445205185438|3.04311296121895
24|13|18|"B "|"P"|6|334.7|29.3939|6|405.3|12.2474|-2.89426691425862|0.861952319419821|0.742961800953209
25|13|19|"B "|"P"|6|397.7|12.2474|6|407.5|39.1918|-0.311566458691052|0.581446713087198|0.338080280159906
26|14|20|"B "|"L"|3|216|36.3731|3|356|39.8372|-2.93620355270232|1.33762126119416|1.78923063839866
27|15|21|"B "|"P"|13|153.4|69.2266|19|163.1|51.8709|-0.159116969627101|0.360514880984436|0.129970979411222
28|16|22|"B "|"M"|5|124.6|2.2361|5|665.6|2.2361|-218.525630021163|54.1029220348491|2927.12617270896
Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
2
$\begingroup$

The meta-regression model fitted via rma() assumes that there is one common $\tau^2$ value within the subgroups, while metagen() with byvar allows $\tau^2$ to differ between subgroups. An illustration:

library(metafor)
library(meta)

dat <- dat.bcg
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat)
dat$sei <- sqrt(dat$vi)

First, we fit the meta-regression model with rma(). I will remove the intercept term, so that the model coefficients are the estimated average effects for each level of the moderator:

rma(yi, vi, data=dat, mods = ~ 0 + alloc)

These are the results:

Mixed-Effects Model (k = 13; tau^2 estimator: REML)

tau^2 (estimated amount of residual heterogeneity):     0.3615 (SE = 0.2111)
tau (square root of estimated tau^2 value):             0.6013
I^2 (residual heterogeneity / unaccounted variability): 88.77%
H^2 (unaccounted variability / sampling variability):   8.91

Test for Residual Heterogeneity:
QE(df = 10) = 132.3676, p-val < .0001

Test of Moderators (coefficients 1:3):
QM(df = 3) = 15.9842, p-val = 0.0011

Model Results:

                 estimate      se     zval    pval    ci.lb    ci.ub 
allocalternate    -0.5180  0.4412  -1.1740  0.2404  -1.3827   0.3468      
allocrandom       -0.9658  0.2672  -3.6138  0.0003  -1.4896  -0.4420  *** 
allocsystematic   -0.4289  0.3449  -1.2434  0.2137  -1.1050   0.2472      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now the same analysis using metagen():

metagen(yi, sei, data=dat, comb.fixed=FALSE, method.tau="REML", byvar=alloc)

These are the results:

Number of studies combined: k = 13

                                         95%-CI     z  p-value
Random effects model -0.7145 [-1.0669; -0.3622] -3.97 < 0.0001

Quantifying heterogeneity:
 tau^2 = 0.3132 [0.1197; 1.1115]; tau = 0.5597 [0.3460; 1.0543];
 I^2 = 92.1% [88.3%; 94.7%]; H = 3.56 [2.93; 4.34]

Quantifying residual heterogeneity:
 I^2 = 92.4% [88.4%; 95.1%]; H = 3.64 [2.94; 4.50]

Test of heterogeneity:
      Q d.f.  p-value
 152.23   12 < 0.0001

Results for subgroups (random effects model):
                     k                     95%-CI  tau^2    tau      Q   I^2
alloc = random       7 -0.9710 [-1.5118; -0.4301] 0.3925 0.6265 110.21 94.6%
alloc = alternate    2 -0.5408 [-1.0927;  0.0111] 0.1326 0.3641   5.56 82.0%
alloc = systematic   4 -0.4242 [-1.1293;  0.2809] 0.4003 0.6327  16.59 81.9%

Test for subgroup differences (random effects model):
                    Q d.f. p-value
Between groups   1.86    2  0.3943

We see minor differences between the Model Results versus the Results for subgroups output.

One can fit a meta-regression model that allows for different $\tau^2$ values across subgroups with:

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

The results are then:

Multivariate Meta-Analysis Model (k = 13; method: REML)

Variance Components:

outer factor: trial (nlvls = 13)
inner factor: alloc (nlvls = 3)

            estim    sqrt  k.lvl  fixed       level 
tau^2.1    0.1326  0.3641      2     no   alternate 
tau^2.2    0.3925  0.6265      7     no      random 
tau^2.3    0.4003  0.6327      4     no  systematic 

Test for Residual Heterogeneity:
QE(df = 10) = 132.3676, p-val < .0001

Test of Moderators (coefficients 1:3):
QM(df = 3) = 17.4587, p-val = 0.0006

Model Results:

                 estimate      se     zval    pval    ci.lb    ci.ub 
allocalternate    -0.5408  0.2816  -1.9204  0.0548  -1.0927   0.0111    . 
allocrandom       -0.9710  0.2760  -3.5185  0.0004  -1.5118  -0.4301  *** 
allocsystematic   -0.4242  0.3597  -1.1792  0.2383  -1.1293   0.2809      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

These results are now identical to the ones from metagen().

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you @Wolfgang, that clarified a lot. I changed my code and got similar results. I have a couple of follow up questions. 1) Is the mods= ~0+alloc the same expression as mods=~alloc-1? 2) If I would want to add a nested structure to the model to account for non-independend data. Could I change the inner|outer structure like this: rma.mv(yi, vi, data=dat, mods = ~ 0 + alloc, random = ~ alloc | study/trial, struct="DIAG") I assume trial is a placeholder for a unique ID. $\endgroup$
    – ipodbielski
    Dec 21, 2020 at 9:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.