4
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ Could you add the code that creates the required dataset without aggregating over BLOCK? If I understand you correctly, this is what you would like to use in the end. $\endgroup$
    – Wolfgang
    Jan 22, 2015 at 13:15
  • $\begingroup$ Sorry, but I am not too sure what you mean. I have modified the dataset to the previous version before aggregating over BLOCK, if that is what you meant. Thanks $\endgroup$
    – fede_luppi
    Jan 22, 2015 at 14:33
  • $\begingroup$ The dataset created by your code has 4 rows, one for each level of SITE. There is no distinction between the various levels of BLOCK anymore. So adding a random effect for BLOCK makes no sense. I thought you meant to analyze a dataset where you have a row for each level of BLOCK within each level of SITE. But maybe I am misunderstanding your question. $\endgroup$
    – Wolfgang
    Jan 23, 2015 at 9:38
  • $\begingroup$ Dear Wolfgang, I think the new edit may be in line with your previous comment. I hope my post is now more understandable. Thanks $\endgroup$
    – fede_luppi
    Jan 23, 2015 at 17:50
  • $\begingroup$ We are getting there. You still need to put the amb and elev values from the same level of SITE and BLOCK in a single row. So, in the end, there should be 10 rows (3 for Europe, Asia, and Australia and 1 for US). $\endgroup$
    – Wolfgang
    Jan 23, 2015 at 18:20

1 Answer 1

Reset to default
6
$\begingroup$

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)

profile likelihood plots

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for such a clear explanation. Indeed, the results are similar to those without including BLOCK as a random effect in a mixed effects model. Besides, in the resulting forest each BLOCK effect within each SITE is shown, which is rather too much information to show. I will probably stick with the model that leaves out BLOCKas a random effect. $\endgroup$
    – fede_luppi
    Jan 26, 2015 at 11:56

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.