This is an example of my dataset
> test SITE YEAR CO2 RING BLOCK NPP 1 Europe 1997 amb 1 1 930.0000 2 Europe 1997 amb 5 3 1214.0000 3 Europe 1997 amb 6 2 990.0000 4 Europe 1997 elev 2 1 1185.0000 5 Europe 1997 elev 3 3 1558.0000 6 Europe 1997 elev 4 2 1223.0000 7 Europe 1998 amb 1 1 830.0000 8 Europe 1998 amb 5 3 1078.0000 9 Europe 1998 amb 6 2 859.0000 10 Europe 1998 elev 2 1 1194.0000 11 Europe 1998 elev 3 3 1329.0000 12 Europe 1998 elev 4 2 1074.0000 13 Europe 1999 amb 1 1 1097.0000 14 Europe 1999 amb 5 3 1129.0000 15 Europe 1999 amb 6 2 1082.0000 16 Europe 1999 elev 2 1 1378.0000 17 Europe 1999 elev 3 3 1495.0000 18 Europe 1999 elev 4 2 1386.0000 19 Europe 2000 amb 1 1 1165.0000 20 Europe 2000 amb 5 3 1229.0000 21 Europe 2000 amb 6 2 1153.0000 22 Europe 2000 elev 2 1 1445.0000 23 Europe 2000 elev 3 3 1632.0000 24 Europe 2000 elev 4 2 1467.0000 25 Europe 2001 amb 1 1 964.0000 26 Europe 2001 amb 5 3 1139.0000 27 Europe 2001 amb 6 2 928.0000 28 Europe 2001 elev 2 1 1196.0000 29 Europe 2001 elev 3 3 1456.0000 30 Europe 2001 elev 4 2 1210.0000 31 Europe 2002 amb 1 1 577.0000 32 Europe 2002 amb 5 3 766.0000 33 Europe 2002 amb 6 2 607.0000 34 Europe 2002 elev 2 1 766.0000 35 Europe 2002 elev 3 3 1223.0000 36 Europe 2002 elev 4 2 749.0000 37 Europe 2003 amb 1 1 964.0000 38 Europe 2003 amb 5 3 937.0000 39 Europe 2003 amb 6 2 862.0000 40 Europe 2003 elev 2 1 1091.0000 41 Europe 2003 elev 3 3 1402.0000 42 Europe 2003 elev 4 2 1267.0000 43 Europe 2004 amb 1 1 960.0000 44 Europe 2004 amb 5 3 1141.0000 45 Europe 2004 amb 6 2 976.0000 46 Europe 2004 elev 2 1 1170.0000 47 Europe 2004 elev 3 3 1613.0000 48 Europe 2004 elev 4 2 1343.0000 49 US 1998 amb 3 1 786.3420 50 US 1998 amb 4 1 845.4510 51 US 1998 amb 5 1 779.8370 52 US 1998 elev 1 1 1002.4910 53 US 1998 elev 2 1 1011.1800 54 US 1999 amb 3 1 777.1900 55 US 1999 amb 4 1 941.4630 56 US 1999 amb 5 1 775.3520 57 US 1999 elev 1 1 927.0810 58 US 1999 elev 2 1 977.8410 59 US 2000 amb 3 1 900.7640 60 US 2000 amb 4 1 1024.2440 61 US 2000 amb 5 1 943.6750 62 US 2000 elev 1 1 1119.3830 63 US 2000 elev 2 1 1108.8620 64 US 2001 amb 3 1 855.8670 65 US 2001 amb 4 1 1047.4500 66 US 2001 amb 5 1 989.7610 67 US 2001 elev 1 1 1271.2010 68 US 2001 elev 2 1 1029.1540 69 US 2002 amb 3 1 829.4200 70 US 2002 amb 4 1 946.2660 71 US 2002 amb 5 1 847.3990 72 US 2002 elev 1 1 1133.0680 73 US 2002 elev 2 1 1029.9070 74 US 2003 amb 3 1 926.6870 75 US 2003 amb 4 1 1082.4480 76 US 2003 amb 5 1 900.4960 77 US 2003 elev 1 1 1031.3780 78 US 2003 elev 2 1 1110.4880 79 US 2004 amb 3 1 761.4971 80 US 2004 amb 4 1 866.0739 81 US 2004 amb 5 1 754.0631 82 US 2004 elev 1 1 874.9087 83 US 2004 elev 2 1 977.7755 84 US 2005 amb 3 1 695.4085 85 US 2005 amb 4 1 880.3115 86 US 2005 amb 5 1 635.9845 87 US 2005 elev 1 1 807.2665 88 US 2005 elev 2 1 839.8101 89 US 2006 amb 3 1 652.3512 90 US 2006 amb 4 1 645.5581 91 US 2006 amb 5 1 542.4291 92 US 2006 elev 1 1 622.3796 93 US 2006 elev 2 1 695.6319 94 US 2007 amb 3 1 619.3959 95 US 2007 amb 4 1 669.2196 96 US 2007 amb 5 1 574.7596 97 US 2007 elev 1 1 616.6445 98 US 2007 elev 2 1 704.6839 99 US 2008 amb 3 1 532.6612 100 US 2008 amb 4 1 567.9665 101 US 2008 amb 5 1 424.7423 102 US 2008 elev 1 1 535.3117 103 US 2008 elev 2 1 573.9919 104 Asia 2000 amb 2 1 1183.6667 105 Asia 2000 amb 3 2 1291.0000 106 Asia 2000 amb 6 3 1224.6667 107 Asia 2000 elev 1 1 1652.6667 108 Asia 2000 elev 4 2 1453.6667 109 Asia 2000 elev 5 3 1326.3333 110 Asia 2001 amb 2 1 1650.0000 111 Asia 2001 amb 3 2 1646.3333 112 Asia 2001 amb 6 3 1712.3333 113 Asia 2001 elev 1 1 2133.6667 114 Asia 2001 elev 4 2 1961.3333 115 Asia 2001 elev 5 3 2200.3333 116 Australia 2001 amb 2 1 236.0000 117 Australia 2001 amb 5 2 217.0000 118 Australia 2001 amb 8 3 235.0000 119 Australia 2001 elev 29 1 611.0000 120 Australia 2001 elev 32 2 319.0000 121 Australia 2001 elev 35 3 409.0000 122 Australia 2002 amb 11 1 322.0000 123 Australia 2002 amb 14 2 270.0000 124 Australia 2002 amb 17 3 324.0000 125 Australia 2002 elev 38 1 538.0000 126 Australia 2002 elev 41 2 436.0000 127 Australia 2002 elev 44 3 579.0000 128 Australia 2003 amb 19 1 409.0000 129 Australia 2003 amb 20 1 339.0000 130 Australia 2003 amb 22 2 432.0000 131 Australia 2003 amb 23 2 291.0000 132 Australia 2003 amb 25 3 431.0000 133 Australia 2003 amb 26 3 314.0000 134 Australia 2003 elev 46 1 544.0000 135 Australia 2003 elev 47 1 545.0000 136 Australia 2003 elev 49 2 434.0000 137 Australia 2003 elev 50 2 433.0000 138 Australia 2003 elev 52 3 646.0000 139 Australia 2003 elev 53 3 507.0000 I want to perform a meta-analysis with this dataset, to estimate an overall effect response, based on the annual NPP (Net Primary Production) measurements in each SITE, which are independent. The effect would be calculated as the ratio for each treatment (elev/amb). The problem is that I would like to consider BLOCK as a random effect, so the effect ratio is calculated as the ratio for each site and within the same BLOCK.
I can rearrange the previous dataset aggregating by BLOCK:
test2<- summarySE(test,measurevar="NPP", groupvars=c("SITE","CO2","BLOCK”)) test3<- reshape(test2[,1:6],v.names=c("NPP","N","sd"),timevar= "CO2",idvar=c("SITE","BLOCK"), direction="wide") > test3 SITE BLOCK NPP.amb N.amb sd.amb NPP.elev N.elev sd.elev 1 Europe 1 935.8750 8 177.54632 1178.1250 8 203.29250 2 Europe 2 932.1250 8 165.83592 1214.8750 8 223.27397 3 Europe 3 1079.1250 8 155.28907 1463.5000 8 141.64039 7 US 1 788.5616 33 164.14384 909.1108 22 207.28113 9 Asia 1 1416.8333 2 329.74746 1893.1667 2 340.11836 10 Asia 2 1468.6667 2 251.25861 1707.5000 2 358.97454 11 Asia 3 1468.5000 2 344.83241 1763.3333 2 618.01133 15 Australia 1 326.5000 4 71.11727 559.5000 4 34.47221 16 Australia 2 302.5000 4 91.77690 405.5000 4 57.68015 17 Australia 3 326.0000 4 80.52743 535.2500 4 101.51642 Please, could you help me design the meta-analysis including the random effect term?
Best Answer
I will tentatively suggest one approach to analyzing these data, without having a full understanding of how the data were collected or the relevance of the YEAR and RING variables in the original dataset. So, I will just assume that the aggregation you have done is reasonable, yielding a final dataset with multiple levels of BLOCK within the various levels of SITE, and the goal is to account for the multilevel structure of these data.
Before we get to the model, we need to compute the observed outcomes and corresponding sampling variances. Based on what you wrote and code you had shown earlier, you are working with log-transformed ratios of means as the outcome measure. So, we would use:
library(metafor) dat <- escalc(measure="ROM", m1i=NPP.amb, sd1i=sd.amb, n1i=N.amb, m2i=NPP.elev, sd2i=sd.elev, n2i=N.elev, data=test3) dat[,c(1,2,9,10)] This yields:
SITE BLOCK yi vi 1 Asia 1 -0.2898 0.0432 2 Asia 2 -0.1507 0.0367 3 Asia 3 -0.1830 0.0890 4 Australia 1 -0.5386 0.0128 5 Australia 2 -0.2930 0.0281 6 Australia 3 -0.4958 0.0242 7 Europe 1 -0.2302 0.0082 8 Europe 2 -0.2649 0.0082 9 Europe 3 -0.3047 0.0038 10 US 1 -0.1423 0.0037 So yi is the vector with the log-transformed ratios of means and vi is the vector with the corresponding sampling variances.
Now we can proceed with the analysis. In essence, you can then use a three-level meta-analytic model along the lines described by Konstantopoulos (2011) and others. We therefore want to add random effects at the SITE level and for the various levels of BLOCK within the levels of SITE. This would mean a random effects structure such as:
res <- rma.mv(yi, vi, random = ~ 1 | SITE/BLOCK, data=dat) res This yields:
Multivariate Meta-Analysis Model (k = 10; method: REML) Variance Components: estim sqrt nlvls fixed factor estim sqrt nlvls fixed factor sigma^2.1 0.0146 0.1207 4 no SITE sigma^2.2 0.0000 0.0000 10 no SITE/BLOCK Test for Heterogeneity: Q(df = 9) = 13.1783, p-val = 0.1547 Model Results: estimate se zval pval ci.lb ci.ub -0.2752 0.0716 -3.8453 0.0001 -0.4154 -0.1349 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 So, what we find is that all of the heterogeneity appears to be between the SITE levels and not among the BLOCK within SITE levels. While the Q-test even suggests that there is no significant amount of heterogeneity in these data to begin with ($p = .15$), this result should be treated with caution as this is a small dataset. So, I would probably stick to these results.
Whenever working with more complex models, it's always a good idea to check the profile likelihood plots of the variance components in the model. You can easily obtain these plots with:
profile(res) 
Not surprisingly, the second plot peaks at zero. More importantly, both plots do have a clear peak, so this gives some reassurance that we are not fitting an overparameterized model.
You can find further discussion of this (and another example of the model) on the metafor package website:
http://www.metafor-project.org/doku.php/analyses:konstantopoulos2011
Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61-76.
Similar Posts:
- Solved – two-way multiple comparisons in R
- Solved – Meta-analysis in R with standard errors instead of standard deviations {metafor}
- Solved – Meta-analysis in R with standard errors instead of standard deviations {metafor}
- Solved – Minimum number of repeated measures and levels per nested random effect
- Solved – Split dataset by categorical variable or use as a dumthe/factor variable