Conducting Three-level Meta-analyses using the metaSEM Package

Table of Contents

1 Introduction

This file illustrates how to conduct three-level meta-analyses using the metaSEM and OpenMx packages available in the R environment. The metaSEM package was written to simplify the procedures to conduct meta-analysis. Most readers may only need to use the metaSEM package to conduct the analysis. The next section shows how to conduct two- and three-level meta-analyses with the meta() and meta3() functions. The third section demonstrates more complicated three-level meta-analyses using a dataset with more predictors. The final section shows how to implement three-level meta-analyses as structural equation models using the OpenMx package. It provides detailed steps on how three-level meta-analyses can be formulated as structural equation models.

This file also demonstrates the advantages of using the SEM approach to conduct three-level meta-analyses. These include flexibility on imposing constraints for model comparisons and construction of likelihood-based confidence interval (LBCI). I also demonstrate how to conduct three-level meta-analysis with restricted (or residual) maximum likelihood (REML) using the reml3() function and handling missing covariates with full information maximum likelihood (FIML) using the meta3X() function. Readers may refer to Cheung (2012a; 2012b) for the design and implementation of the metaSEM package and Cheung (2013) for the theory and issues on how to formulate three-level meta-analyses as structural equation models.

Two datasets from published meta-analyses were used in the illustrations. The first dataset was based on Cooper et al. (2003) and Konstantopoulos (2011). Konstantopoulos (2011) selected part of the dataset to illustrate how to conduct three-level meta-analysis. The second dataset was reported by Bornmann et al. (2007) and Marsh et al. (2009). They conducted a three-level meta-analysis on gender effects in peer reviews of grant proposals.

2 Comparisons between Two- and Three-Level Models with Cooper et al.'s (2003) Dataset

As an illustration, I first conduct the tradition (two-level) meta-analysis using the meta() function. Then I conduct a three-level meta-analysis using the meta3() function. We may compare the similarities and differences between these two sets of results.

2.1 Inspecting the data

Before running the analyses, we need to load the metaSEM library. The datasets are stored in the library. It is always a good idea to inspect the data before the analyses. We may display the first few cases of the dataset by using the head() command.

#### Cooper et al. (2003)

library("metaSEM")
head(Cooper03)
Loading required package: OpenMx
  District Study     y     v Year
1       11     1 -0.18 0.118 1976
2       11     2 -0.22 0.118 1976
3       11     3  0.23 0.144 1976
4       11     4 -0.30 0.144 1976
5       12     5  0.13 0.014 1989
6       12     6 -0.26 0.014 1989

2.2 Two-level meta-analysis

Similar to other R packages, we may use summary() to extract the results after running the analyses. I first conduct a random-effects meta-analysis and then a fixed- and mixed-effects meta-analyses.

  1. Random-effects model

    The Q statistic on testing the homogeneity of effect sizes was 578.86, df = 55, p < .001. The estimated heterogeneity \(\tau^2\) (labeled Tau2_1_1 in the output) and \(I^2\) were 0.0866 and 0.9459, respectively. This indicates that the between-study effect explains about 95% of the total variation. The average population effect (labeled Intercept1 in the output; and its 95% Wald CI) was 0.1280 (0.0428, 0.2132).

    #### Two-level meta-analysis
    
    ## Random-effects model  
    summary( meta(y=y, v=v, data=Cooper03) )
    
    Running Meta analysis with ML 
    
    Call:
    meta(y = y, v = v, data = Cooper03)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
    Intercept1 0.128003  0.043472 0.042799 0.213207  2.9445  0.003235 ** 
    Tau2_1_1   0.086537  0.019485 0.048346 0.124728  4.4411 8.949e-06 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 578.864
    Degrees of freedom of the Q statistic: 55
    P value of the Q statistic: 0
    Heterogeneity indices (based on the estimated Tau2):
                                 Estimate
    Intercept1: I2 (Q statistic)   0.9459
    
    Number of studies (or clusters): 56
    Number of observed statistics: 56
    Number of estimated parameters: 2
    Degrees of freedom: 54
    -2 log likelihood: 33.2919 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    
  2. Fixed-effects model

    A fixed-effects meta-analysis can be conducted by fixing the heterogeneity of the random effects at 0 with the RE.constraints argument (random-effects constraints). The estimated common effect (and its 95% Wald CI) was 0.0464 (0.0284, 0.0644).

    ## Fixed-effects model
    summary( meta(y=y, v=v, data=Cooper03, RE.constraints=0) )
    
    Running Meta analysis with ML 
    
    Call:
    meta(y = y, v = v, data = Cooper03, RE.constraints = 0)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                Estimate Std.Error    lbound    ubound z value Pr(>|z|)    
    Intercept1 0.0464072 0.0091897 0.0283957 0.0644186  5.0499 4.42e-07 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 578.864
    Degrees of freedom of the Q statistic: 55
    P value of the Q statistic: 0
    Heterogeneity indices (based on the estimated Tau2):
                                 Estimate
    Intercept1: I2 (Q statistic)        0
    
    Number of studies (or clusters): 56
    Number of observed statistics: 56
    Number of estimated parameters: 1
    Degrees of freedom: 55
    -2 log likelihood: 434.2075 
    OpenMx status1: 1 ("0" and "1": considered fine; other values indicate problems)
    
  3. Mixed-effects model

    Year was used as a covariate. It is easier to interpret the intercept by centering the Year with scale(Year, scale=FALSE). The scale=FALSE argument means that it is centered, but not standardized. The estimated regression coefficient (labeled Slope1_1 in the output; and its 95% Wald CI) was 0.0051 (-0.0033, 0.0136) which is not significant at \(\alpha=.05\). The \(R^2\) is 0.0164.

    ## Mixed-effects model
    summary( meta(y=y, v=v, x=scale(Year, scale=FALSE), data=Cooper03) )
    
    Running Meta analysis with ML 
    
    Call:
    meta(y = y, v = v, x = scale(Year, scale = FALSE), data = Cooper03)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                 Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
    Intercept1  0.1259125  0.0432028  0.0412367  0.2105884  2.9145  0.003563 ** 
    Slope1_1    0.0051307  0.0043248 -0.0033457  0.0136071  1.1864  0.235483    
    Tau2_1_1    0.0851153  0.0190462  0.0477856  0.1224451  4.4689 7.862e-06 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 578.864
    Degrees of freedom of the Q statistic: 55
    P value of the Q statistic: 0
    Explained variances (R2):
                               y1
    Tau2 (no predictor)    0.0865
    Tau2 (with predictors) 0.0851
    R2                     0.0164
    
    Number of studies (or clusters): 56
    Number of observed statistics: 56
    Number of estimated parameters: 3
    Degrees of freedom: 53
    -2 log likelihood: 31.88635 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    

2.3 Three-level meta-analysis

  1. Random-effects model

    The Q statistic on testing the homogeneity of effect sizes was the same as that under the two-level meta-analysis. The estimated heterogeneity at level 2 \(\tau^2_{(2)}\) (labeled Tau2_2 in the output) and at level 3 \(\tau^2_{(3)}\) (labeled Tau2_3 in the output) were 0.0329 and 0.0577, respectively. The level 2 \(I^2_{(2)}\) (labeled I2_2 in the output) and the level 3 \(I^2_{(3)}\) (labeled I2_3 in the output) were 0.3440 and 0.6043, respectively. Schools (level 2) and districts (level 3) explain about 34% and 60% of the total variation, respectively. The average population effect (and its 95% Wald CI) was 0.1845 (0.0266, 0.3423).

    #### Three-level meta-analysis
    
    ## Random-effects model
    summary( meta3(y=y, v=v, cluster=District, data=Cooper03) )
    
    Running Meta analysis with ML 
    
    Call:
    meta3(y = y, v = v, cluster = District, data = Cooper03)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
    Intercept  0.1844554  0.0805411  0.0265977  0.3423130  2.2902 0.022010 * 
    Tau2_2     0.0328648  0.0111397  0.0110314  0.0546982  2.9502 0.003175 **
    Tau2_3     0.0577384  0.0307423 -0.0025154  0.1179921  1.8781 0.060362 . 
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 578.864
    Degrees of freedom of the Q statistic: 55
    P value of the Q statistic: 0
    Heterogeneity indices (based on the estimated Tau2):
                                  Estimate
    I2_2 (Typical v: Q statistic)   0.3440
    I2_3 (Typical v: Q statistic)   0.6043
    
    Number of studies (or clusters): 11
    Number of observed statistics: 56
    Number of estimated parameters: 3
    Degrees of freedom: 53
    -2 log likelihood: 16.78987 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    
  2. Mixed-effects model

    Year was used as a covariate. The estimated regression coefficient (labeled Slope_1 in the output; and its 95% Wald CI) was 0.0051 (-0.0116, 0.0218) which is not significant at \(\alpha=.05\). The estimated level 2 \(R^2_{(2)}\) and level 3 \(R^2_{(3)}\) were 0.0000 and 0.0221, respectively.

    ## Mixed-effects model
    summary( meta3(y=y, v=v, cluster=District, x=scale(Year, scale=FALSE), 
                   data=Cooper03) )
    
     Running Meta analysis with ML 
    
    Call:
    meta3(y = y, v = v, cluster = District, x = scale(Year, scale = FALSE), 
        data = Cooper03)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
    Intercept  0.1780268  0.0805219  0.0202067  0.3358469  2.2109 0.027042 * 
    Slope_1    0.0050737  0.0085266 -0.0116382  0.0217856  0.5950 0.551814   
    Tau2_2     0.0329390  0.0111620  0.0110618  0.0548162  2.9510 0.003168 **
    Tau2_3     0.0564628  0.0300330 -0.0024007  0.1153264  1.8800 0.060104 . 
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 578.864
    Degrees of freedom of the Q statistic: 55
    P value of the Q statistic: 0
    Explained variances (R2):
                            Level 2 Level 3
    Tau2 (no predictor)    0.032865  0.0577
    Tau2 (with predictors) 0.032939  0.0565
    R2                     0.000000  0.0221
    
    Number of studies (or clusters): 11
    Number of observed statistics: 56
    Number of estimated parameters: 4
    Degrees of freedom: 52
    -2 log likelihood: 16.43629 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    

2.4 Model comparisons

Many research hypotheses involve model comparisons among nested models. anova(), a generic function to comparing nested models, may be used to conduct a likelihood ratio test which is also known as a chi-square difference test.

  1. Testing \(H_0: \tau^2_{(3)} = 0\)
    • Based on the data structure, it is clear that a 3-level meta-analysis is preferred to a traditional 2-level meta-analysis. It is still of interest to test whether the 3-level model is statistically better than the 2-level model by testing \(H_0: \tau^2_{(3)}=0\). Since the models with \(\tau^2_{(3)}\) being freely estimated and with \(\tau^2_{(3)}=0\) are nested, we may compare them by the use of a likelihood ratio test.
    • It should be noted, however, that \(H_0: \tau^2_{(3)}=0\) is tested on the boundary. The likelihood ratio test is not distributed as a chi-square variate with 1 df. A simple strategy to correct this bias is to reject the null hypothesis when the observed p value is larger than .10 for \(\alpha=.05\).
    • The likelihood-ratio test was 16.5020 (df =1), p < .001. This clearly demonstrates that the three-level model is statistically better than the two-level model.
    ## Model comparisons
    
    model2 <- meta(y=y, v=v, data=Cooper03, model.name="2 level model", silent=TRUE) 
    #### An equivalent model by fixing tau2 at level 3=0 in meta3()
    ## model2 <- meta3(y=y, v=v, cluster=District, data=Cooper03, 
    ##                 model.name="2 level model", RE3.constraints=0) 
    model3 <- meta3(y=y, v=v, cluster=District, data=Cooper03, 
                    model.name="3 level model", silent=TRUE) 
    anova(model3, model2)
    
               base    comparison ep minus2LL df       AIC   diffLL diffdf
    1 3 level model          <NA>  3 16.78987 53 -89.21013       NA     NA
    2 3 level model 2 level model  2 33.29190 54 -74.70810 16.50203      1
                 p
    1           NA
    2 4.859793e-05
    
  2. Testing \(H_0: \tau^2_{(2)} = \tau^2_{(3)}\)
    • From the results of the 3-level random-effects meta-analysis, it appears the level 3 heterogeneity is much larger than that at level 2.
    • We may test the null hypothesis \(H_0: \tau^2_{(2)} = \tau^2_{(3)}\) by the use of a likelihood-ratio test.
    • We may impose an equality constraint on \(\tau^2_{(2)} = \tau^2_{(3)}\) by using the same label in meta3(). For example, Eq_tau2 is used as the label in RE2.constraints and RE3.constraints meaning that both the level 2 and level 3 random effects heterogeneity variances are constrained equally. The value of 0.1 was used as the starting value in the constraints.
    • The likelihood-ratio test was 0.6871 (df = 1), p = 0.4072. This indicates that there is not enough evidence to reject \(H_0: \tau^2_2=\tau^2_3\).
    ## Testing \tau^2_2 = \tau^2_3
    modelEqTau2 <- meta3(y=y, v=v, cluster=District, data=Cooper03, 
                         model.name="Equal tau2 at both levels",
                         RE2.constraints="0.1*Eq_tau2", RE3.constraints="0.1*Eq_tau2") 
    anova(model3, modelEqTau2)
    
     Running Equal tau2 at both levels
               base                comparison ep minus2LL df       AIC    diffLL
    1 3 level model                      <NA>  3 16.78987 53 -89.21013        NA
    2 3 level model Equal tau2 at both levels  2 17.47697 54 -90.52303 0.6870959
      diffdf         p
    1     NA        NA
    2      1 0.4071539
    

2.5 Likelihood-based confidence interval

  • A Wald CI is constructed by \(\hat{\theta} \pm 1.96 SE\) where \(\hat{\theta}\) and \(SE\) are the parameter estimate and its estimated standard error.
  • A LBCI can be constructed by the use of the likelihood ratio statistic (e.g., Cheung, 2009; Neal & Miller, 1997).
  • It is well known that the performance of Wald CI on variance components is very poor. For example, the 95% Wald CI on \(\hat{\tau}^2_{(3)}\) in the three-level random-effects meta-analysis was (-0.0025, 0.1180). The lower bound falls outside 0.
  • A LBCI on the heterogeneity variance is preferred. Since \(I^2_{(2)}\) and \(I^2_{(3)}\) are functions of \(\tau^2_{(2)}\) and \(\tau^2_{(3)}\), LBCI on these indices may also be requested and used to indicate the precision of these indices.
  • LBCI may be requested by specifying LB in the intervals.type argument.
  • The 95% LBCI on \(\hat{\tau}^2_{(3)}\) is (0.0198, 0.1763) that stay inside the meaningful boundaries. Regarding the \(I^2\), the 95% LBCIs on \(I^2_{(2)}\) and \(I^2_{(3)}\) were (0.1274, 0.6573) and (0.2794, 0.8454), respectively.
## Likelihood-based CI
summary( meta3(y=y, v=v, cluster=District, data=Cooper03, 
               I2=c("I2q", "ICC"), intervals.type="LB") )
 Running Meta analysis with ML 

Call:
meta3(y = y, v = v, cluster = District, data = Cooper03, intervals.type = "LB", 
    I2 = c("I2q", "ICC"))

95% confidence intervals: Likelihood-based statistic
Coefficients:
          Estimate Std.Error   lbound   ubound z value Pr(>|z|)
Intercept 0.184455        NA 0.011668 0.358535      NA       NA
Tau2_2    0.032865        NA 0.016330 0.063126      NA       NA
Tau2_3    0.057738        NA 0.019805 0.176328      NA       NA

Q statistic on homogeneity of effect sizes: 578.864
Degrees of freedom of the Q statistic: 55
P value of the Q statistic: 0
Heterogeneity indices (I2) and their 95% likelihood-based CIs:
                               lbound Estimate ubound
I2_2 (Typical v: Q statistic) 0.12736  0.34396 0.6573
ICC_2 (tau^2/(tau^2+tau^3))   0.13102  0.36273 0.7015
I2_3 (Typical v: Q statistic) 0.27937  0.60429 0.8454
ICC_3 (tau^3/(tau^2+tau^3))   0.29847  0.63727 0.8690

Number of studies (or clusters): 11
Number of observed statistics: 56
Number of estimated parameters: 3
Degrees of freedom: 53
-2 log likelihood: 16.78987 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  • A LBCI may also be requested in mixed-effects meta-analysis.
summary( meta3(y=y, v=v, cluster=District, x=scale(Year, scale=FALSE), 
               data=Cooper03, intervals.type="LB") )
 Running Meta analysis with ML 

Call:
meta3(y = y, v = v, cluster = District, x = scale(Year, scale = FALSE), 
    data = Cooper03, intervals.type = "LB")

95% confidence intervals: Likelihood-based statistic
Coefficients:
            Estimate Std.Error     lbound     ubound z value Pr(>|z|)
Intercept  0.1780268        NA  0.0052929  0.3515622      NA       NA
Slope_1    0.0050737        NA -0.0128356  0.0237979      NA       NA
Tau2_2     0.0329390        NA  0.0163732  0.0632779      NA       NA
Tau2_3     0.0564628        NA  0.0192314  0.1720022      NA       NA

Q statistic on homogeneity of effect sizes: 578.864
Degrees of freedom of the Q statistic: 55
P value of the Q statistic: 0
Explained variances (R2):
                        Level 2 Level 3
Tau2 (no predictor)    0.032865  0.0577
Tau2 (with predictors) 0.032939  0.0565
R2                     0.000000  0.0221

Number of studies (or clusters): 11
Number of observed statistics: 56
Number of estimated parameters: 4
Degrees of freedom: 52
-2 log likelihood: 16.43629 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)

2.6 Restricted maximum likelihood estimation

  • REML may also be used in three-level meta-analysis. The parameter estimates for \(\tau^2_{(2)}\) and \(\tau^2_{(3)}\) were 0.0327 and 0.0651, respectively.
## REML
summary( reml1 <- reml3(y=y, v=v, cluster=District, data=Cooper03) )
Running Variance component with REML 

Call:
reml3(y = y, v = v, cluster = District, data = Cooper03)

95% confidence intervals: z statistic approximation
Coefficients:
         Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Tau2_2  0.0327365  0.0110922  0.0109963  0.0544768  2.9513 0.003164 **
Tau2_3  0.0650619  0.0355102 -0.0045368  0.1346607  1.8322 0.066921 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of studies (or clusters): 56
Number of observed statistics: 55
Number of estimated parameters: 2
Degrees of freedom: 53
-2 log likelihood: -81.14044 
OpenMx status: 0 ("0" and "1": considered fine; other values indicate problems)
  • We may impose an equality constraint on \(\tau^2_{(2)}\) and \(\tau^2_{(3)}\) and test whether this constraint is statistically significant. The estimated value for \(\tau^2_{(2)}=\tau^2_{(3)}\) was 0.0404. When this model is compared against the unconstrained model, the test statistic was 1.0033 (df = 1), p = .3165, which is not significant.
summary( reml0 <- reml3(y=y, v=v, cluster=District, data=Cooper03,
                        RE.equal=TRUE, model.name="Equal Tau2") )
anova(reml1, reml0)
 Running Equal Tau2 

Call:
reml3(y = y, v = v, cluster = District, data = Cooper03, RE.equal = TRUE, 
    model.name = "Equal Tau2")

95% confidence intervals: z statistic approximation
Coefficients:
     Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
Tau2 0.040418  0.010290 0.020249 0.060587  3.9277 8.576e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of studies (or clusters): 56
Number of observed statistics: 55
Number of estimated parameters: 1
Degrees of freedom: 54
-2 log likelihood: -80.1371 
OpenMx status: 0 ("0" and "1": considered fine; other values indicate problems)
                          base comparison ep  minus2LL df AIC   diffLL diffdf
1 Variance component with REML       <NA>  2 -81.14044 53  NA       NA     NA
2 Variance component with REML Equal Tau2  1 -80.13710 54  NA 1.003336      1
          p
1        NA
2 0.3165046
  • We may also estimate the residual heterogeneity after controlling for the covariate. The estimated residual heterogeneity for \(\tau^2_{(2)}\) and \(\tau^2_{(3)}\) were 0.0327 and 0.0723, respectively.
summary( reml3(y=y, v=v, cluster=District, x=scale(Year, scale=FALSE),
               data=Cooper03) )
 Running Variance component with REML 

Call:
reml3(y = y, v = v, cluster = District, x = scale(Year, scale = FALSE), 
    data = Cooper03)

95% confidence intervals: z statistic approximation
Coefficients:
         Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Tau2_2  0.0326502  0.0110529  0.0109870  0.0543134  2.9540 0.003137 **
Tau2_3  0.0722656  0.0405349 -0.0071813  0.1517125  1.7828 0.074619 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of studies (or clusters): 56
Number of observed statistics: 54
Number of estimated parameters: 2
Degrees of freedom: 52
-2 log likelihood: -72.09405 
OpenMx status: 0 ("0" and "1": considered fine; other values indicate problems)

3 More Complex 3-Level Meta-Analyses with Bornmann et al.'s (2007) Dataset

This section replicates the findings in Table 3 of Marsh et al. (2009). Several additional analyses on model comparisons were conducted. Missing data were artificially introduced to illustrate how missing data might be handled with FIML.

3.1 Inspecting the data

The effect size and its sampling variance are logOR (log of the odds ratio) and v, respectively. Cluster is the variable representing the cluster effect, whereas the potential covariates are Year (year of publication), Type (Grants vs. Fellowship), Discipline (Physical sciences, Life sciences/biology, Social sciences/humanities and Multidisciplinary) and Country (United States, Canada, Australia, United Kingdom and Europe).

#### Bornmann et al. (2007)

library("metaSEM")
head(Bornmann07)
  Id                       Study Cluster    logOR          v Year       Type
1  1 Ackers (2000a; Marie Curie)       1 -0.40108 0.01391692 1996 Fellowship
2  2 Ackers (2000b; Marie Curie)       1 -0.05727 0.03428793 1996 Fellowship
3  3 Ackers (2000c; Marie Curie)       1 -0.29852 0.03391122 1996 Fellowship
4  4 Ackers (2000d; Marie Curie)       1  0.36094 0.03404025 1996 Fellowship
5  5 Ackers (2000e; Marie Curie)       1 -0.33336 0.01282103 1996 Fellowship
6  6 Ackers (2000f; Marie Curie)       1 -0.07173 0.01361189 1996 Fellowship
                  Discipline Country
1          Physical sciences  Europe
2          Physical sciences  Europe
3          Physical sciences  Europe
4          Physical sciences  Europe
5 Social sciences/humanities  Europe
6          Physical sciences  Europe

3.2 Model 0: Intercept

The Q statistic was 221.2809 (df = 65), p < .001. The estimated average effect (and its 95% Wald CI) was -0.1008 (-0.1794, -0.0221). The \(\hat{\tau}^2_{(2)}\) and \(\hat{\tau}^3_{(3)}\) were 0.0038 and 0.0141, respectively. The \(I^2_{(2)}\) and \(I^2_{(3)}\) were 0.1568 and 0.5839, respectively.

## Model 0: Intercept  
summary( Model0 <- meta3(y=logOR, v=v, cluster=Cluster, data=Bornmann07, 
                         model.name="3 level model") )
 Running 3 level model 

Call:
meta3(y = logOR, v = v, cluster = Cluster, data = Bornmann07, 
    model.name = "3 level model")

95% confidence intervals: z statistic approximation
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)  
Intercept -0.1007784  0.0401327 -0.1794371 -0.0221198 -2.5111  0.01203 *
Tau2_2     0.0037965  0.0027210 -0.0015367  0.0091297  1.3952  0.16295  
Tau2_3     0.0141352  0.0091445 -0.0037877  0.0320580  1.5458  0.12216  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Heterogeneity indices (based on the estimated Tau2):
                              Estimate
I2_2 (Typical v: Q statistic)   0.1568
I2_3 (Typical v: Q statistic)   0.5839

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 3
Degrees of freedom: 63
-2 log likelihood: 25.80256 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Testing \(H_0: \tau^2_3 = 0\)

    We may test whether the three-level model is necessary by testing \(H_0: \tau^2_{(3)} = 0\). The likelihood ratio statistic was 10.2202 (df = 1), p < .01. Thus, the three-level model is statistically better than the two-level model.

    ## Testing tau^2_3 = 0
    Model0a <- meta3(logOR, v, cluster=Cluster, data=Bornmann07, 
                     RE3.constraints=0, model.name="2 level model")
    anova(Model0, Model0a)
    
     Running 2 level model
               base    comparison ep minus2LL df        AIC   diffLL diffdf
    1 3 level model          <NA>  3 25.80256 63 -100.19744       NA     NA
    2 3 level model 2 level model  2 36.02279 64  -91.97721 10.22024      1
                p
    1          NA
    2 0.001389081
    
  2. Testing \(H_0: \tau^2_2 = \tau^2_3\)

    The likelihood-ratio statistic in testing \(H_0: \tau^2_{(2)} = \tau^2_{(3)}\) was 1.3591 (df = 1), p = 0.2437. Thus, there is no evidence to reject the null hypothesis.

    ## Testing tau^2_2 = tau^2_3
    Model0b <- meta3(logOR, v, cluster=Cluster, data=Bornmann07, 
                     RE2.constraints="0.1*Eq_tau2", RE3.constraints="0.1*Eq_tau2", 
                     model.name="tau2_2 equals tau2_3")
    anova(Model0, Model0b)
    
     Running tau2_2 equals tau2_3
               base           comparison ep minus2LL df       AIC   diffLL diffdf
    1 3 level model                 <NA>  3 25.80256 63 -100.1974       NA     NA
    2 3 level model tau2_2 equals tau2_3  2 27.16166 64 -100.8383 1.359103      1
             p
    1       NA
    2 0.243693
    

3.3 Model 1: Type as a covariate

  • Conventionally, one level (e.g., Grants) is used as the reference group. The estimated intercept (labeled Intercept in the output) represents the estimated effect size for Grants and the regression coefficient (labeled Slope_1 in the output) is the difference between Fellowship and Grants.
  • The estimated slope on Type (and its 95% Wald CI) was -0.1956 (-0.3018, -0.0894) which is statistically significant at \(\alpha=.05\). This is the difference between Fellowship and Grants. The \(R^2_{(2)}\) and \(R^2_{(3)}\) were 0.0693 and 0.7943, respectively.
## Model 1: Type as a covariate  
## Convert characters into a dummy variable
## Type2=0 (Grants); Type2=1 (Fellowship)    
Type2 <- ifelse(Bornmann07$Type=="Fellowship", yes=1, no=0)
summary( Model1 <- meta3(logOR, v, x=Type2, cluster=Cluster, data=Bornmann07))
Running Meta analysis with ML 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = Type2, data = Bornmann07)

95% confidence intervals: z statistic approximation
Coefficients:
            Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
Intercept -0.0066071  0.0371125 -0.0793462  0.0661320 -0.1780 0.8587001    
Slope_1   -0.1955869  0.0541649 -0.3017483 -0.0894256 -3.6110 0.0003051 ***
Tau2_2     0.0035335  0.0024306 -0.0012303  0.0082974  1.4538 0.1460058    
Tau2_3     0.0029079  0.0031183 -0.0032039  0.0090197  0.9325 0.3510704    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0035335  0.0029
R2                     0.0692595  0.7943

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 4
Degrees of freedom: 62
-2 log likelihood: 17.62569 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Alternative model: Grants and Fellowship as indicator variables

    If we want to estimate the effects for both Grants and Fellowship, we may create two indicator variables to represent them. Since we cannot estimate the intercept and these coefficients at the same time, we need to fix the intercept at 0 by specifying the intercept.constraints=0 argument in meta3(). We are now able to include both Grants and Fellowship in the analysis. The estimated effects (and their 95% CIs) for Grants and Fellowship were -0.0066 (-0.0793, 0.0661) and -0.2022 (-0.2805, -0.1239), respectively.

    ## Alternative model: Grants and Fellowship as indicators  
    ## Indicator variables
    Grants <- ifelse(Bornmann07$Type=="Grants", yes=1, no=0)
    Fellowship <- ifelse(Bornmann07$Type=="Fellowship", yes=1, no=0)
    
    summary(Model1b <- meta3(logOR, v, x=cbind(Grants, Fellowship), 
                             cluster=Cluster, data=Bornmann07,
                             intercept.constraints=0, model.name="Model 1"))
    
     Running Model 1 
    
    Call:
    meta3(y = logOR, v = v, cluster = Cluster, x = cbind(Grants, 
        Fellowship), data = Bornmann07, intercept.constraints = 0, 
        model.name = "Model 1")
    
    95% confidence intervals: z statistic approximation
    Coefficients:
              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
    Slope_1 -0.0066071  0.0371125 -0.0793462  0.0661320 -0.1780    0.8587    
    Slope_2 -0.2021940  0.0399473 -0.2804893 -0.1238987 -5.0615 4.159e-07 ***
    Tau2_2   0.0035335  0.0024306 -0.0012303  0.0082974  1.4538    0.1460    
    Tau2_3   0.0029079  0.0031183 -0.0032039  0.0090197  0.9325    0.3511    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 221.2809
    Degrees of freedom of the Q statistic: 65
    P value of the Q statistic: 0
    Explained variances (R2):
                             Level 2 Level 3
    Tau2 (no predictor)    0.0037965  0.0141
    Tau2 (with predictors) 0.0035335  0.0029
    R2                     0.0692595  0.7943
    
    Number of studies (or clusters): 21
    Number of observed statistics: 66
    Number of estimated parameters: 4
    Degrees of freedom: 62
    -2 log likelihood: 17.62569 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    

3.4 Model 2: Year and Year^2 as covariates

  • When there are several covariates, we may combine them with the cbind() command. For example, cbind(Year, Year^2) includes both Year and its squared as covariates. In the output, Slope_1 and Slope_2 refer to the regression coefficients for Year and Year^2, respectively. To increase the numerical stability, the covariates are usually centered before creating the quadratic terms. Since Marsh et al. (2009) standardized Year in their analysis, I follow this practice here.
  • The estimated regression coefficients (and their 95% CIs) for Year and Year^2 were -0.0010 (-0.0473, 0.0454) and -0.0118 (-0.0247, 0.0012), respectively. The \(R^2_{(2)}\) and \(R^2_{(3)}\) were 0.2430 and 0.0000, respectively.
## Model 2: Year and Year^2 as covariates
summary( Model2 <- meta3(logOR, v, x=cbind(scale(Year), scale(Year)^2), 
                         cluster=Cluster, data=Bornmann07,
                         model.name="Model 2") )
 Running Model 2 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(scale(Year), 
    scale(Year)^2), data = Bornmann07, model.name = "Model 2")

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate   Std.Error      lbound      ubound z value Pr(>|z|)  
Intercept -0.08627312  0.04125581 -0.16713302 -0.00541322 -2.0912  0.03651 *
Slope_1   -0.00095287  0.02365224 -0.04731040  0.04540466 -0.0403  0.96786  
Slope_2   -0.01176840  0.00659995 -0.02470407  0.00116727 -1.7831  0.07457 .
Tau2_2     0.00287389  0.00206817 -0.00117965  0.00692744  1.3896  0.16466  
Tau2_3     0.01479446  0.00926095 -0.00335666  0.03294558  1.5975  0.11015  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0028739  0.0148
R2                     0.2430134  0.0000

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 5
Degrees of freedom: 61
-2 log likelihood: 22.3836 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Testing \(H_0: \beta_{Year} = \beta_{Year^2}=0\)

    The test statistic was 3.4190 (df = 2), p = 0.1810. Thus, there is no evidence supporting that Year has a quadratic effect on the effect size.

    ## Testing beta_{Year} = beta_{Year^2}=0
    anova(Model2, Model0)
    
         base    comparison ep minus2LL df       AIC   diffLL diffdf         p
    1 Model 2          <NA>  5 22.38360 61  -99.6164       NA     NA        NA
    2 Model 2 3 level model  3 25.80256 63 -100.1974 3.418955      2 0.1809603
    

3.5 Model 3: Discipline as a covariate

  • There are four categories in Discipline. multidisciplinary is used as the reference group in the analysis.
  • The estimated regression coefficients (and their 95% Wald CIs) for DisciplinePhy, DisciplineLife and DisciplineSoc were -0.0091 (-0.2041, 0.1859), -0.1262 (-0.2804, 0.0280) and -0.2370 (-0.4746, 0.0007), respectively. The \(R^2_2\) and \(R^2_3\) were 0.0000 and 0.4975, respectively.
## Model 3: Discipline as a covariate
## Create dummy variables using multidisciplinary as the reference group
DisciplinePhy <- ifelse(Bornmann07$Discipline=="Physical sciences", yes=1, no=0)
DisciplineLife <- ifelse(Bornmann07$Discipline=="Life sciences/biology", yes=1, no=0)
DisciplineSoc <- ifelse(Bornmann07$Discipline=="Social sciences/humanities", yes=1, no=0)
summary( Model3 <- meta3(logOR, v, x=cbind(DisciplinePhy, DisciplineLife, DisciplineSoc), 
                         cluster=Cluster, data=Bornmann07,
                         model.name="Model 3") )
 Running Model 3 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(DisciplinePhy, 
    DisciplineLife, DisciplineSoc), data = Bornmann07, model.name = "Model 3")

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate   Std.Error      lbound      ubound z value Pr(>|z|)  
Intercept -0.01474783  0.06389945 -0.13998846  0.11049280 -0.2308  0.81747  
Slope_1   -0.00913064  0.09949198 -0.20413133  0.18587005 -0.0918  0.92688  
Slope_2   -0.12617957  0.07866274 -0.28035572  0.02799657 -1.6041  0.10870  
Slope_3   -0.23695698  0.12123091 -0.47456520  0.00065125 -1.9546  0.05063 .
Tau2_2     0.00390942  0.00283948 -0.00165587  0.00947471  1.3768  0.16857  
Tau2_3     0.00710338  0.00643210 -0.00550331  0.01971006  1.1044  0.26944  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0039094  0.0071
R2                     0.0000000  0.4975

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 6
Degrees of freedom: 60
-2 log likelihood: 20.07571 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Testing whether Discipline is significant

    The test statistic was 5.7268 (df = 3), p = 0.1257 which is not significant. Therefore, there is no evidence supporting that Discipline explains the variation of the effect sizes.

    ## Testing whether Discipline is significant
    anova(Model3, Model0)
    
         base    comparison ep minus2LL df        AIC   diffLL diffdf         p
    1 Model 3          <NA>  6 20.07571 60  -99.92429       NA     NA        NA
    2 Model 3 3 level model  3 25.80256 63 -100.19744 5.726842      3 0.1256832
    

3.6 Model 4: Country as a covariate

  • There are five categories in Country. United States is used as the reference group in the analysis.
  • The estimated regression coefficients (and their 95% Wald CIs) for CountryAus, CountryCan and CountryEur CountryUK are -0.0240 (-0.2405, 0.1924), -0.1341 (-0.3674, 0.0993), -0.2211 (-0.3660, -0.0762) and 0.0537 (-0.1413, 0.2487), respectively. The \(R^2_2\) and \(R^2_3\) were 0.1209 and 0.6606, respectively.
## Model 4: Country as a covariate
## Create dummy variables using the United States as the reference group
CountryAus <- ifelse(Bornmann07$Country=="Australia", yes=1, no=0)
CountryCan <- ifelse(Bornmann07$Country=="Canada", yes=1, no=0)
CountryEur <- ifelse(Bornmann07$Country=="Europe", yes=1, no=0)
CountryUK <- ifelse(Bornmann07$Country=="United Kingdom", yes=1, no=0)

summary( Model4 <- meta3(logOR, v, x=cbind(CountryAus, CountryCan, CountryEur, 
                         CountryUK), cluster=Cluster, data=Bornmann07,
                         model.name="Model 4") )
 Running Model 4 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(CountryAus, 
    CountryCan, CountryEur, CountryUK), data = Bornmann07, model.name = "Model 4")

95% confidence intervals: z statistic approximation
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Intercept  0.0025681  0.0597768 -0.1145923  0.1197285  0.0430 0.965732   
Slope_1   -0.0240109  0.1104328 -0.2404552  0.1924333 -0.2174 0.827876   
Slope_2   -0.1340800  0.1190667 -0.3674465  0.0992866 -1.1261 0.260127   
Slope_3   -0.2210801  0.0739174 -0.3659556 -0.0762046 -2.9909 0.002782 **
Slope_4    0.0537251  0.0994803 -0.1412527  0.2487030  0.5401 0.589157   
Tau2_2     0.0033376  0.0023492 -0.0012667  0.0079420  1.4208 0.155383   
Tau2_3     0.0047979  0.0044818 -0.0039862  0.0135820  1.0705 0.284379   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0033376  0.0048
R2                     0.1208598  0.6606

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 7
Degrees of freedom: 59
-2 log likelihood: 14.18259 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Testing whether Discipline is significant

    The test statistic was 11.6200 (df = 4), p = 0.0204 which is statistically significant.

    ## Testing whether Discipline is significant
    anova(Model4, Model0)
    
         base    comparison ep minus2LL df       AIC   diffLL diffdf          p
    1 Model 4          <NA>  7 14.18259 59 -103.8174       NA     NA         NA
    2 Model 4 3 level model  3 25.80256 63 -100.1974 11.61996      4 0.02041284
    

3.7 Model 5: Type and Discipline as covariates

The \(R^2_{(2)}\) and \(R^2_{(3)}\) were 0.3925 and 1.0000, respectively. The \(\hat{\tau}^2_{(3)}\) was near 0 after controlling for the covariates.

## Model 5: Type and Discipline as covariates
summary( Model5 <- meta3(logOR, v, x=cbind(Type2, DisciplinePhy, DisciplineLife, 
                         DisciplineSoc), cluster=Cluster, data=Bornmann07,
                         model.name="Model 5") )
 Running Model 5 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(Type2, DisciplinePhy, 
    DisciplineLife, DisciplineSoc), data = Bornmann07, model.name = "Model 5")

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
Intercept  6.7036e-02  1.8553e-02  3.0672e-02  1.0340e-01  3.6132 0.0003025 ***
Slope_1   -1.9004e-01  4.0234e-02 -2.6890e-01 -1.1118e-01 -4.7233  2.32e-06 ***
Slope_2    1.9511e-02  6.5942e-02 -1.0973e-01  1.4875e-01  0.2959 0.7673209    
Slope_3   -1.2779e-01  3.5914e-02 -1.9818e-01 -5.7400e-02 -3.5582 0.0003734 ***
Slope_4   -2.3950e-01  9.4054e-02 -4.2384e-01 -5.5154e-02 -2.5464 0.0108849 *  
Tau2_2     2.3062e-03  1.4270e-03 -4.9059e-04  5.1030e-03  1.6162 0.1060586    
Tau2_3     1.0000e-10          NA          NA          NA      NA        NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0023062  0.0000
R2                     0.3925434  1.0000

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 7
Degrees of freedom: 59
-2 log likelihood: 4.66727 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Testing whether Discipline is significant after controlling for Type

    The test statistic was 12.9584 (df = 3), p = 0.0047 which is significant. Therefore, Discipline is still significant after controlling for Type.

    ## Testing whether Discipline is significant after controlling for Type
    anova(Model5, Model1)
    
         base            comparison ep minus2LL df       AIC   diffLL diffdf
    1 Model 5                  <NA>  7  4.66727 59 -113.3327       NA     NA
    2 Model 5 Meta analysis with ML  4 17.62569 62 -106.3743 12.95842      3
                p
    1          NA
    2 0.004727388
    

3.8 Model 6: Type and Country as covariates

The \(R^2_{(2)}\) and \(R^2_{(3)}\) were 0.3948 and 1.0000, respectively. The \(\hat{\tau}^2_{(3)}\) was near 0 after controlling for the covariates.

## Model 6: Type and Country as covariates
summary( Model6 <- meta3(logOR, v, x=cbind(Type2, CountryAus, CountryCan, 
                         CountryEur, CountryUK), cluster=Cluster, data=Bornmann07,
                         model.name="Model 6") )
 Running Model 6 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(Type2, CountryAus, 
    CountryCan, CountryEur, CountryUK), data = Bornmann07, model.name = "Model 6")

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
Intercept  6.7507e-02  1.8933e-02  3.0399e-02  1.0461e-01  3.5656 0.0003631 ***
Slope_1   -1.5167e-01  4.1113e-02 -2.3225e-01 -7.1092e-02 -3.6892 0.0002250 ***
Slope_2   -6.9580e-02  8.5164e-02 -2.3650e-01  9.7339e-02 -0.8170 0.4139267    
Slope_3   -1.4231e-01  7.5204e-02 -2.8970e-01  5.0879e-03 -1.8923 0.0584497 .  
Slope_4   -1.6116e-01  4.0203e-02 -2.3995e-01 -8.2361e-02 -4.0086 6.108e-05 ***
Slope_5    9.0419e-03  7.0074e-02 -1.2830e-01  1.4639e-01  0.1290 0.8973315    
Tau2_2     2.2976e-03  1.4407e-03 -5.2618e-04  5.1213e-03  1.5947 0.1107693    
Tau2_3     1.0000e-10          NA          NA          NA      NA        NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0022976  0.0000
R2                     0.3948192  1.0000

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 8
Degrees of freedom: 58
-2 log likelihood: 5.076592 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
  1. Testing whether Country is significant after controlling for Type

    The test statistic was 12.5491 (df = 4), p = 0.0137. Thus, Country is significant after controlling for Type.

    ## Testing whether Country is significant after controlling for Type
    anova(Model6, Model1)
    
         base            comparison ep  minus2LL df       AIC  diffLL diffdf
    1 Model 6                  <NA>  8  5.076592 58 -110.9234      NA     NA
    2 Model 6 Meta analysis with ML  4 17.625692 62 -106.3743 12.5491      4
               p
    1         NA
    2 0.01370262
    

3.9 Model 7: Discipline and Country as covariates

The \(R^2_{(2)}\) and \(R^2_{(3)}\) were 0.1397 and 0.7126, respectively.

## Model 7: Discipline and Country as covariates
summary( meta3(logOR, v, x=cbind(DisciplinePhy, DisciplineLife, DisciplineSoc,
                         CountryAus, CountryCan, CountryEur, CountryUK), 
                         cluster=Cluster, data=Bornmann07,
                         model.name="Model 7") )
 Running Model 7 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(DisciplinePhy, 
    DisciplineLife, DisciplineSoc, CountryAus, CountryCan, CountryEur, 
    CountryUK), data = Bornmann07, model.name = "Model 7")

95% confidence intervals: z statistic approximation
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)  
Intercept  0.0029305  0.0576743 -0.1101090  0.1159700  0.0508  0.95948  
Slope_1    0.1742169  0.1702554 -0.1594774  0.5079113  1.0233  0.30618  
Slope_2    0.0826806  0.1599802 -0.2308748  0.3962359  0.5168  0.60528  
Slope_3   -0.0462265  0.1715773 -0.3825118  0.2900589 -0.2694  0.78761  
Slope_4   -0.0486321  0.1306918 -0.3047834  0.2075191 -0.3721  0.70981  
Slope_5   -0.2169132  0.1915703 -0.5923841  0.1585576 -1.1323  0.25751  
Slope_6   -0.3036578  0.1670720 -0.6311129  0.0237974 -1.8175  0.06914 .
Slope_7   -0.0605272  0.1809419 -0.4151668  0.2941124 -0.3345  0.73799  
Tau2_2     0.0032661  0.0022784 -0.0011994  0.0077317  1.4335  0.15171  
Tau2_3     0.0040618  0.0038459 -0.0034759  0.0115996  1.0562  0.29090  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0032661  0.0041
R2                     0.1396973  0.7126

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 10
Degrees of freedom: 56
-2 log likelihood: 10.31105 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)

3.10 Model 8: Type, Discipline and Country as covariates

The \(R^2_{(2)}\) and \(R^2_{(3)}\) were 0.4466 and 1.0000, respectively. The \(\hat{\tau}^2_{(3)}\) was near 0 after controlling for the covariates.

## Model 8: Type, Discipline and Country as covariates
Model8 <- meta3(logOR, v, x=cbind(Type2, DisciplinePhy, DisciplineLife, DisciplineSoc,
                           CountryAus, CountryCan, CountryEur, CountryUK), 
                           cluster=Cluster, data=Bornmann07,
                           model.name="Model 8") 
## There was an estimation error. The model was rerun again.
summary(rerun(Model8))
 Running Model 8
Running Model 8 

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(Type2, DisciplinePhy, 
    DisciplineLife, DisciplineSoc, CountryAus, CountryCan, CountryEur, 
    CountryUK), data = Bornmann07, model.name = "Model 8")

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
Intercept  6.8563e-02  1.8630e-02  3.2049e-02  1.0508e-01  3.6802  0.000233 ***
Slope_1   -1.6885e-01  4.1545e-02 -2.5028e-01 -8.7425e-02 -4.0643 4.818e-05 ***
Slope_2    2.5329e-01  1.5814e-01 -5.6670e-02  5.6325e-01  1.6016  0.109240    
Slope_3    1.2689e-01  1.4774e-01 -1.6268e-01  4.1646e-01  0.8589  0.390411    
Slope_4   -8.3549e-03  1.5796e-01 -3.1795e-01  3.0124e-01 -0.0529  0.957818    
Slope_5   -1.1530e-01  1.1147e-01 -3.3377e-01  1.0317e-01 -1.0344  0.300948    
Slope_6   -2.6412e-01  1.6402e-01 -5.8559e-01  5.7344e-02 -1.6103  0.107324    
Slope_7   -2.9029e-01  1.4859e-01 -5.8152e-01  9.5221e-04 -1.9536  0.050754 .  
Slope_8   -1.5975e-01  1.6285e-01 -4.7893e-01  1.5943e-01 -0.9810  0.326609    
Tau2_2     2.1010e-03  1.2925e-03 -4.3226e-04  4.6342e-03  1.6255  0.104051    
Tau2_3     1.0000e-10          NA          NA          NA      NA        NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on homogeneity of effect sizes: 221.2809
Degrees of freedom of the Q statistic: 65
P value of the Q statistic: 0
Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0037965  0.0141
Tau2 (with predictors) 0.0021010  0.0000
R2                     0.4466073  1.0000

Number of studies (or clusters): 21
Number of observed statistics: 66
Number of estimated parameters: 11
Degrees of freedom: 55
-2 log likelihood: -1.645211 
OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)

3.11 Handling missing covariates with FIML

When there are missing data in the covariates, data with missing values are excluded before the analysis in meta3(). The missing covariates can be handled by the use of FIML in meta3X(). We illustrate two examples on how to analyze data with missing covariates with missing completely at random (MCAR) and missing at random (MAR) data.

  1. MCAR

    About 25% of the level-2 covariate Type was introduced by the MCAR mechanism.

    #### Handling missing covariates with FIML
    
    ## MCAR
    ## Set seed for replication
    set.seed(1000000)
    
    ## Copy Bornmann07 to my.df
    my.df <- Bornmann07
    ## "Fellowship": 1; "Grant": 0
    my.df$Type_MCAR <- ifelse(Bornmann07$Type=="Fellowship", yes=1, no=0)
    
    ## Create 17 out of 66 missingness with MCAR
    my.df$Type_MCAR[sample(1:66, 17)] <- NA
    
    summary(meta3(y=logOR, v=v, cluster=Cluster, x=Type_MCAR, data=my.df))
    
    Running Meta analysis with ML 
    
    Call:
    meta3(y = logOR, v = v, cluster = Cluster, x = Type_MCAR, data = my.df)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                 Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
    Intercept -0.00484542  0.03934429 -0.08195881  0.07226797 -0.1232    0.9020    
    Slope_1   -0.21090081  0.05346221 -0.31568483 -0.10611680 -3.9449 7.985e-05 ***
    Tau2_2     0.00446788  0.00549282 -0.00629784  0.01523361  0.8134    0.4160    
    Tau2_3     0.00092884  0.00336491 -0.00566625  0.00752394  0.2760    0.7825    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 151.643
    Degrees of freedom of the Q statistic: 48
    P value of the Q statistic: 1.115552e-12
    Explained variances (R2):
                             Level 2 Level 3
    Tau2 (no predictor)    0.0042664  0.0145
    Tau2 (with predictors) 0.0044679  0.0009
    R2                     0.0000000  0.9361
    
    Number of studies (or clusters): 20
    Number of observed statistics: 49
    Number of estimated parameters: 4
    Degrees of freedom: 45
    -2 log likelihood: 13.13954 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    

    There is no need to specify whether the covariates are level 2 or level 3 in meta3() because the covariates are treated as a design matrix. When meta3X() is used, users need to specify whether the covariates are at level 2 (x2) or level 3 (x3).

    summary(meta3X(y=logOR, v=v, cluster=Cluster, x2=Type_MCAR, data=my.df))
    
    Running Meta analysis with ML 
    
    Call:
    meta3X(y = logOR, v = v, cluster = Cluster, x2 = Type_MCAR, data = my.df)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
    Intercept -0.0106318  0.0397685 -0.0885766  0.0673131 -0.2673 0.789207   
    SlopeX2_1 -0.1753249  0.0582619 -0.2895161 -0.0611336 -3.0093 0.002619 **
    Tau2_2     0.0030338  0.0026839 -0.0022266  0.0082941  1.1304 0.258324   
    Tau2_3     0.0036839  0.0042817 -0.0047082  0.0120759  0.8604 0.389586   
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Explained variances (R2):
                             Level 2 Level 3
    Tau2 (no predictor)    0.0037965  0.0141
    Tau2 (with predictors) 0.0030338  0.0037
    R2                     0.2009070  0.7394
    
    Number of studies (or clusters): 21
    Number of observed statistics: 115
    Number of estimated parameters: 7
    Degrees of freedom: 108
    -2 log likelihood: 49.76372 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    
  2. MAR

    For the case for missing covariates with MAR, the missingness in Type depends on the values of Year. Type is missing when Year is smaller than 1996.

    ## MAR
    Type_MAR <- ifelse(Bornmann07$Type=="Fellowship", yes=1, no=0)
    
    ## Create 27 out of 66 missingness with MAR for cases Year<1996
    index_MAR <- ifelse(Bornmann07$Year<1996, yes=TRUE, no=FALSE)
    Type_MAR[index_MAR] <- NA
    
    summary(meta3(logOR, v, x=Type_MAR, cluster=Cluster, data=Bornmann07))
    
    Running Meta analysis with ML 
    
    Call:
    meta3(y = logOR, v = v, cluster = Cluster, x = Type_MAR, data = Bornmann07)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                 Estimate   Std.Error      lbound      ubound z value Pr(>|z|)   
    Intercept -0.01587052  0.03952546 -0.09333901  0.06159796 -0.4015 0.688032   
    Slope_1   -0.17573647  0.06328327 -0.29976940 -0.05170355 -2.7770 0.005487 **
    Tau2_2     0.00259266  0.00171596 -0.00077056  0.00595588  1.5109 0.130811   
    Tau2_3     0.00278384  0.00267150 -0.00245221  0.00801989  1.0421 0.297388   
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Q statistic on homogeneity of effect sizes: 151.11
    Degrees of freedom of the Q statistic: 38
    P value of the Q statistic: 1.998401e-15
    Explained variances (R2):
                             Level 2 Level 3
    Tau2 (no predictor)    0.0029593  0.0097
    Tau2 (with predictors) 0.0025927  0.0028
    R2                     0.1238926  0.7121
    
    Number of studies (or clusters): 12
    Number of observed statistics: 39
    Number of estimated parameters: 4
    Degrees of freedom: 35
    -2 log likelihood: -24.19956 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    

    It is possible to include level 2 (av2) and level 3 (av3) auxiliary variables. Auxiliary variables are those that predict the missing values or are correlated with the variables that contain missing values. The inclusion of auxiliary variables can improve the efficiency of the estimation and the parameter estimates.

    ## Include auxiliary variable
    summary(meta3X(y=logOR, v=v, cluster=Cluster, x2=Type_MAR, av2=Year, data=my.df))
    
    Running Meta analysis with ML 
    
    Call:
    meta3X(y = logOR, v = v, cluster = Cluster, x2 = Type_MAR, av2 = Year, 
        data = my.df)
    
    95% confidence intervals: z statistic approximation
    Coefficients:
                Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
    Intercept -0.0264057  0.0572068 -0.1385291  0.0857176 -0.4616 0.644380   
    SlopeX2_1 -0.2003999  0.0691101 -0.3358532 -0.0649466 -2.8997 0.003735 **
    Tau2_2     0.0029970  0.0022371 -0.0013877  0.0073817  1.3397 0.180359   
    Tau2_3     0.0030212  0.0032463 -0.0033415  0.0093839  0.9306 0.352036   
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Explained variances (R2):
                             Level 2 Level 3
    Tau2 (no predictor)    0.0049237  0.0088
    Tau2 (with predictors) 0.0029970  0.0030
    R2                     0.3913242  0.6571
    
    Number of studies (or clusters): 21
    Number of observed statistics: 171
    Number of estimated parameters: 14
    Degrees of freedom: 157
    -2 log likelihood: 377.3479 
    OpenMx status1: 0 ("0" and "1": considered fine; other values indicate problems)
    

4 Implementing Three-Level Meta-Analyses as Structural Equation Models in OpenMx

This section illustrates how to formulate three-level meta-analyses as structural equation models using the OpenMx package. The steps outline how to create the model-implied mean vector and the model-implied covariance matrix to fit the three-level meta-analyses. y is the effect size (standardized mean difference on the modified school calendars) and v is its sampling variance. District is the variable for the cluster effect, whereas Year is the year of publication.

4.1 Preparing data

  • Data in a three-level meta-analysis are usually stored in the long format, e.g., Cooper03 in this example, whereas the SEM approach uses the wide format.
  • Suppose the maximum number of effect sizes in the level-2 unit is \(k\) (\(k=11\) in this example). Each cluster is represented by one row with \(k=11\) variables representing the outcome effect size, say y1 to y11 in this example. The incomplete data are represented by NA (missing value).
  • Similarly, \(k=11\) variables are required to represent the known sampling variances, say v1 to v11 in this example.
  • If the covariates are at level 2, \(k=11\) variables are also required to represent each of them. For example, Year is a level-2 covariate, Year1 to Year11 are required to represent it.
  • Several extra steps are required to handle missing values. Missing values (represented by NA in R) are not allowed in v1 to v11 as they are definition variables. The missing data are converted into some arbitrary values, say 1e10 in this example. The actual value does not matter because the missing values will be removed before the analysis. It is because missing values in y1 to y11 (and v1 to v11) will be filtered out automatically by the use of FIML.
#### Steps in Analyzing Three-level Meta-analysis in OpenMx

#### Preparing data
## Load the library
library("OpenMx")

## Get the dataset from the metaSEM library
data(Cooper03, package="metaSEM")

## Make a copy of the original data
my.long <- Cooper03

## Show the first few cases in my.long
head(my.long)
  District Study     y     v Year
1       11     1 -0.18 0.118 1976
2       11     2 -0.22 0.118 1976
3       11     3  0.23 0.144 1976
4       11     4 -0.30 0.144 1976
5       12     5  0.13 0.014 1989
6       12     6 -0.26 0.014 1989
## Center the Year to increase numerical stability
my.long$Year <- scale(my.long$Year, scale=FALSE)

## maximum no. of effect sizes in level-2
k <- 11

## Create a variable called "time" to store: 1, 2, 3, ... k
my.long$time <- c(unlist(sapply(split(my.long$y, my.long$District), 
                                function(x) 1:length(x))))

## Convert long format to wide format by "District"
my.wide <- reshape(my.long, timevar="time", idvar=c("District"), 
                   sep="_", direction="wide")

## NA in v is due to NA in y in wide format
## Replace NA with 1e10 in "v"
temp <- my.wide[, paste("v_", 1:k, sep="")]
temp[is.na(temp)] <- 1e10
my.wide[, paste("v_", 1:k, sep="")] <- temp

## Replace NA with 0 in "Year"
temp <- my.wide[, paste("Year_", 1:k, sep="")]
temp[is.na(temp)] <- 0
my.wide[, paste("Year_", 1:k, sep="")] <- temp

## Show the first few cases in my.wide
head(my.wide)
   District Study_1   y_1   v_1      Year_1 Study_2   y_2   v_2      Year_2
1        11       1 -0.18 0.118 -13.5535714       2 -0.22 0.118 -13.5535714
5        12       5  0.13 0.014  -0.5535714       6 -0.26 0.014  -0.5535714
9        18       9  0.45 0.023   4.4464286      10  0.38 0.043   4.4464286
12       27      12  0.16 0.020 -13.5535714      13  0.65 0.004 -13.5535714
16       56      16  0.08 0.019   7.4464286      17  0.04 0.007   7.4464286
20       58      20 -0.18 0.020 -13.5535714      21  0.00 0.018 -13.5535714
   Study_3  y_3   v_3      Year_3 Study_4   y_4      v_4      Year_4 Study_5
1        3 0.23 0.144 -13.5535714       4 -0.30 1.44e-01 -13.5535714      NA
5        7 0.19 0.015  -0.5535714       8  0.32 2.40e-02  -0.5535714      NA
9       11 0.29 0.012   4.4464286      NA    NA 1.00e+10   0.0000000      NA
12      14 0.36 0.004 -13.5535714      15  0.60 7.00e-03 -13.5535714      NA
16      18 0.19 0.005   7.4464286      19 -0.06 4.00e-03   7.4464286      NA
20      22 0.00 0.019 -13.5535714      23 -0.28 2.20e-02 -13.5535714      24
     y_5   v_5    Year_5 Study_6  y_6     v_6    Year_6 Study_7  y_7   v_7
1     NA 1e+10   0.00000      NA   NA 1.0e+10   0.00000      NA   NA 1e+10
5     NA 1e+10   0.00000      NA   NA 1.0e+10   0.00000      NA   NA 1e+10
9     NA 1e+10   0.00000      NA   NA 1.0e+10   0.00000      NA   NA 1e+10
12    NA 1e+10   0.00000      NA   NA 1.0e+10   0.00000      NA   NA 1e+10
16    NA 1e+10   0.00000      NA   NA 1.0e+10   0.00000      NA   NA 1e+10
20 -0.04 2e-02 -13.55357      25 -0.3 2.1e-02 -13.55357      26 0.07 6e-03
      Year_7 Study_8 y_8   v_8    Year_8 Study_9  y_9   v_9    Year_9 Study_10
1    0.00000      NA  NA 1e+10   0.00000      NA   NA 1e+10   0.00000       NA
5    0.00000      NA  NA 1e+10   0.00000      NA   NA 1e+10   0.00000       NA
9    0.00000      NA  NA 1e+10   0.00000      NA   NA 1e+10   0.00000       NA
12   0.00000      NA  NA 1e+10   0.00000      NA   NA 1e+10   0.00000       NA
16   0.00000      NA  NA 1e+10   0.00000      NA   NA 1e+10   0.00000       NA
20 -13.55357      27   0 7e-03 -13.55357      28 0.05 7e-03 -13.55357       29
    y_10  v_10   Year_10 Study_11  y_11  v_11   Year_11
1     NA 1e+10   0.00000       NA    NA 1e+10   0.00000
5     NA 1e+10   0.00000       NA    NA 1e+10   0.00000
9     NA 1e+10   0.00000       NA    NA 1e+10   0.00000
12    NA 1e+10   0.00000       NA    NA 1e+10   0.00000
16    NA 1e+10   0.00000       NA    NA 1e+10   0.00000
20 -0.08 7e-03 -13.55357       30 -0.09 7e-03 -13.55357

4.2 Random-effects model

  • To implement a three-level meta-analysis as a structural equation model, we need to specify both the model-implied mean vector \(\mu(\theta)\), say expMean, and the model-implied covariance matrix \(\Sigma(\theta)\), say expCov.
  • When there is no covariate, the expected mean is a \(k \times 1\) vector with all elements of beta0 (the intercept), i.e., \( \mu(\theta) = \left[ \begin{array}{c} 1 \\ \vdots \\ 1 \end{array} \right]\beta_0 \). Since OpenMx expects a row vector rather than a column vector in the model-implied means, we need to transpose the expMean in the analysis.
  • Tau2 (\(T^2_{(2)}\)) and Tau3 (\(T^2_{(3)}\)) are the level 2 and level 3 matrices of heterogeneity, respectively. Tau2 is a diagonal matrix with elements of \(\tau^2_{(2)}\), whereas Tau3 is a full matrix with elements of \(\tau^2_{(3)}\). V is a diagonal matrix of the known sampling variances \(v_{ij}\).
  • The model-implied covariance matrix is \( \Sigma(\theta) = T^2_{(3)} + T^2_{(2)} + V \).
  • All of these matrices are stored into a model called random.model.
#### Random-effects model  
## Intercept
Beta0 <- mxMatrix("Full", ncol=1, nrow=1, free=TRUE, labels="beta0", 
                  name="Beta0")
## 1 by k row vector of ones
Ones <- mxMatrix("Unit", nrow=k, ncol=1, name="Ones")

## Model implied mean vector 
## OpenMx expects a row vector rather than a column vector.
expMean <- mxAlgebra( t(Ones %*% Beta0), name="expMean")

## Tau2_2
Tau2 <- mxMatrix("Symm", ncol=1, nrow=1, values=0.01, free=TRUE, labels="tau2_2", 
                 name="Tau2")
Tau3 <- mxMatrix("Symm", ncol=1, nrow=1, values=0.01, free=TRUE, labels="tau2_3", 
                 name="Tau3")

## k by k identity matrix
Iden <- mxMatrix("Iden", nrow=k, ncol=k, name="Iden")

## Conditional sampling variances
## data.v_1, data.v_2, ... data.v_k represent values for definition variables
V <- mxMatrix("Diag", nrow=k, ncol=k, free=FALSE, 
              labels=paste("data.v", 1:k, sep="_"), name="V")

## Model implied covariance matrix
expCov <- mxAlgebra( Ones%*% Tau3 %*% t(Ones) + Iden %x% Tau2 + V, name="expCov")

## Model stores everthing together
random.model <- mxModel(model="Random effects model", 
                        mxData(observed=my.wide, type="raw"), 
                        Iden, Ones, Beta0, Tau2, Tau3, V, expMean, expCov,
                        mxFIMLObjective("expCov","expMean", 
                        dimnames=paste("y", 1:k, sep="_")))
  • We perform a random-effects three-level meta-analysis by running the model with the mxRun() command. The parameter estimates (and their SEs) for \(\beta_0\), \(\tau^2_{(2)}\) and \(\tau^2_{(3)}\) were 0.1845 (0.0805), 0.0329 (0.0111) and 0.0577 (0.0307), respectively.
summary( mxRun(random.model) )
Running Random effects model 
data:
$`Random effects model.data`
    District        Study_1           y_1                v_1         
 Min.   : 11.0   Min.   : 1.00   Min.   :-0.52000   Min.   :0.00100  
 1st Qu.: 22.5   1st Qu.:10.50   1st Qu.:-0.12500   1st Qu.:0.01450  
 Median : 58.0   Median :20.00   Median : 0.12000   Median :0.02000  
 Mean   :107.5   Mean   :24.64   Mean   : 0.07182   Mean   :0.03255  
 3rd Qu.: 88.5   3rd Qu.:38.00   3rd Qu.: 0.23000   3rd Qu.:0.02700  
 Max.   :644.0   Max.   :53.00   Max.   : 0.50000   Max.   :0.11800  
                                                                     
     Year_1           Study_2           y_2               v_2         
 Min.   :-13.554   Min.   : 2.00   Min.   :-0.2600   Min.   :0.00100  
 1st Qu.: -7.054   1st Qu.:11.50   1st Qu.:-0.0250   1st Qu.:0.00900  
 Median :  5.446   Median :21.00   Median : 0.3800   Median :0.01400  
 Mean   :  1.083   Mean   :25.64   Mean   : 0.3173   Mean   :0.03091  
 3rd Qu.:  7.446   3rd Qu.:39.00   3rd Qu.: 0.6550   3rd Qu.:0.03700  
 Max.   : 10.446   Max.   :54.00   Max.   : 0.9800   Max.   :0.11800  
                                                                      
     Year_2           Study_3           y_3               v_3         
 Min.   :-13.554   Min.   : 3.00   Min.   :-0.0300   Min.   :0.00100  
 1st Qu.: -7.054   1st Qu.:12.50   1st Qu.: 0.0200   1st Qu.:0.00750  
 Median :  5.446   Median :22.00   Median : 0.1900   Median :0.01200  
 Mean   :  1.083   Mean   :26.64   Mean   : 0.2409   Mean   :0.02882  
 3rd Qu.:  7.446   3rd Qu.:40.00   3rd Qu.: 0.2600   3rd Qu.:0.02450  
 Max.   : 10.446   Max.   :55.00   Max.   : 1.1900   Max.   :0.14400  
                                                                      
     Year_3            Study_4           y_4                v_4           
 Min.   :-13.5536   Min.   : 4.00   Min.   :-0.30000   Min.   :0.000e+00  
 1st Qu.: -7.0536   1st Qu.:15.00   1st Qu.:-0.06000   1st Qu.:0.000e+00  
 Median :  4.4464   Median :23.00   Median : 0.00000   Median :0.000e+00  
 Mean   :  0.9919   Mean   :28.67   Mean   : 0.05778   Mean   :1.818e+09  
 3rd Qu.:  7.4464   3rd Qu.:45.00   3rd Qu.: 0.27000   3rd Qu.:0.000e+00  
 Max.   : 10.4464   Max.   :56.00   Max.   : 0.60000   Max.   :1.000e+10  
                    NA's   :2       NA's   :2                             
     Year_4             Study_5          y_5              v_5           
 Min.   :-13.55357   Min.   :24.0   Min.   :-0.340   Min.   :0.000e+00  
 1st Qu.: -7.05357   1st Qu.:34.5   1st Qu.:-0.115   1st Qu.:0.000e+00  
 Median :  0.00000   Median :42.0   Median :-0.035   Median :1.000e+10  
 Mean   : -0.08929   Mean   :40.0   Mean   :-0.090   Mean   :6.364e+09  
 3rd Qu.:  7.44643   3rd Qu.:47.5   3rd Qu.:-0.010   3rd Qu.:1.000e+10  
 Max.   : 10.44643   Max.   :52.0   Max.   : 0.050   Max.   :1.000e+10  
                     NA's   :7      NA's   :7                           
     Year_5           Study_6        y_6                v_6           
 Min.   :-13.554   Min.   :25   Min.   :-0.30000   Min.   :0.000e+00  
 1st Qu.:  0.000   1st Qu.:32   1st Qu.:-0.15000   1st Qu.:5.000e+09  
 Median :  0.000   Median :39   Median : 0.00000   Median :1.000e+10  
 Mean   :  1.344   Mean   :37   Mean   :-0.07667   Mean   :7.273e+09  
 3rd Qu.:  3.723   3rd Qu.:43   3rd Qu.: 0.03500   3rd Qu.:1.000e+10  
 Max.   : 10.446   Max.   :47   Max.   : 0.07000   Max.   :1.000e+10  
                   NA's   :8    NA's   :8                             
     Year_6            Study_7          y_7             v_7           
 Min.   :-13.5536   Min.   :26.0   Min.   :0.010   Min.   :0.000e+00  
 1st Qu.:  0.0000   1st Qu.:29.5   1st Qu.:0.025   1st Qu.:1.000e+10  
 Median :  0.0000   Median :33.0   Median :0.040   Median :1.000e+10  
 Mean   :  0.3945   Mean   :33.0   Mean   :0.040   Mean   :8.182e+09  
 3rd Qu.:  0.0000   3rd Qu.:36.5   3rd Qu.:0.055   3rd Qu.:1.000e+10  
 Max.   : 10.4464   Max.   :40.0   Max.   :0.070   Max.   :1.000e+10  
                    NA's   :9      NA's   :9                          
     Year_7            Study_8          y_8              v_8           
 Min.   :-13.5536   Min.   :27.0   Min.   :-0.100   Min.   :0.000e+00  
 1st Qu.:  0.0000   1st Qu.:30.5   1st Qu.:-0.075   1st Qu.:1.000e+10  
 Median :  0.0000   Median :34.0   Median :-0.050   Median :1.000e+10  
 Mean   : -0.5552   Mean   :34.0   Mean   :-0.050   Mean   :8.182e+09  
 3rd Qu.:  0.0000   3rd Qu.:37.5   3rd Qu.:-0.025   3rd Qu.:1.000e+10  
 Max.   :  7.4464   Max.   :41.0   Max.   : 0.000   Max.   :1.000e+10  
                    NA's   :9      NA's   :9                           
     Year_8            Study_9        y_9            v_9           
 Min.   :-13.5536   Min.   :28   Min.   :0.05   Min.   :0.000e+00  
 1st Qu.:  0.0000   1st Qu.:28   1st Qu.:0.05   1st Qu.:1.000e+10  
 Median :  0.0000   Median :28   Median :0.05   Median :1.000e+10  
 Mean   : -0.5552   Mean   :28   Mean   :0.05   Mean   :9.091e+09  
 3rd Qu.:  0.0000   3rd Qu.:28   3rd Qu.:0.05   3rd Qu.:1.000e+10  
 Max.   :  7.4464   Max.   :28   Max.   :0.05   Max.   :1.000e+10  
                    NA's   :10   NA's   :10                        
     Year_9           Study_10       y_10            v_10          
 Min.   :-13.554   Min.   :29   Min.   :-0.08   Min.   :0.000e+00  
 1st Qu.:  0.000   1st Qu.:29   1st Qu.:-0.08   1st Qu.:1.000e+10  
 Median :  0.000   Median :29   Median :-0.08   Median :1.000e+10  
 Mean   : -1.232   Mean   :29   Mean   :-0.08   Mean   :9.091e+09  
 3rd Qu.:  0.000   3rd Qu.:29   3rd Qu.:-0.08   3rd Qu.:1.000e+10  
 Max.   :  0.000   Max.   :29   Max.   :-0.08   Max.   :1.000e+10  
                   NA's   :10   NA's   :10                         
    Year_10           Study_11       y_11            v_11          
 Min.   :-13.554   Min.   :30   Min.   :-0.09   Min.   :0.000e+00  
 1st Qu.:  0.000   1st Qu.:30   1st Qu.:-0.09   1st Qu.:1.000e+10  
 Median :  0.000   Median :30   Median :-0.09   Median :1.000e+10  
 Mean   : -1.232   Mean   :30   Mean   :-0.09   Mean   :9.091e+09  
 3rd Qu.:  0.000   3rd Qu.:30   3rd Qu.:-0.09   3rd Qu.:1.000e+10  
 Max.   :  0.000   Max.   :30   Max.   :-0.09   Max.   :1.000e+10  
                   NA's   :10   NA's   :10                         
    Year_11       
 Min.   :-13.554  
 1st Qu.:  0.000  
 Median :  0.000  
 Mean   : -1.232  
 3rd Qu.:  0.000  
 Max.   :  0.000  
                  

free parameters:
    name matrix row col   Estimate  Std.Error lbound ubound
1  beta0  Beta0   1   1 0.18445537 0.08054111              
2 tau2_2   Tau2   1   1 0.03286479 0.01113968              
3 tau2_3   Tau3   1   1 0.05773836 0.03074229              

observed statistics:  56 
estimated parameters:  3 
degrees of freedom:  53 
-2 log likelihood:  16.78987 
saturated -2 log likelihood:  NA 
number of observations:  11 
chi-square:  NA 
p:  NA 
Information Criteria: 
     df Penalty Parameters Penalty Sample-Size Adjusted
AIC:  -89.21013           22.78987                   NA
BIC: -110.29858           23.98356             14.95056
CFI: NA 
TLI: NA 
RMSEA:  NA 
timestamp: 2013-11-12 22:55:46 
frontend time: 0.07711315 secs 
backend time: 0.01225734 secs 
independent submodels time: 3.385544e-05 secs 
wall clock time: 0.08940434 secs 
cpu time: 0.08940434 secs 
openmx version number: 1.3.2-2301

4.3 Mixed-effects model

  • We may extend a random-effects model to a mixed-effects model by including a covariate (Year in this example).
  • beta1 is the regression coefficient, whereas X stores the value of Year via definition variables.
  • The conditional model-implied mean vector is \( \mu(\theta|Year_{ij}) = \left[ \begin{array}{c} 1 \\ \vdots \\ 1 \end{array} \right]\beta_0 + \left[ \begin{array}{c} Year_{1j} \\ \vdots \\ Year_{kj} \end{array} \right]\beta_1 \).
  • The conditional model-implied covariance matrix is the same as that in the random-effects model, i.e., \( \Sigma(\theta|Year_{ij}) = T^2_{(3)} + T^2_{(2)} + V \).
#### Mixed-effects model

## Design matrix via definition variable
X <- mxMatrix("Full", nrow=k, ncol=1, free=FALSE, 
              labels=paste("data.Year_", 1:k, sep=""), name="X")

## Regression coefficient
Beta1 <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=0,
                  labels="beta1", name="Beta1")

## Model implied mean vector
expMean <- mxAlgebra( t(Ones%*%Beta0 + X%*%Beta1), name="expMean")

mixed.model <- mxModel(model="Mixed effects model", 
                       mxData(observed=my.wide, type="raw"), 
                       Iden, Ones, Beta0, Beta1, Tau2, Tau3, V, expMean, expCov, 
                       X, mxFIMLObjective("expCov","expMean", 
                       dimnames=paste("y", 1:k, sep="_")))
  • The parameter estimates (and their SEs) for \(\beta_0\), \(\beta_1\), \(\tau^2_2\) and \(\tau^2_3\) were 0.1780 (0.0805), 0.0051 (0.0085), 0.0329 (0.0112) and 0.0565 (0.0300), respectively.
summary ( mxRun(mixed.model) )
Running Mixed effects model 
data:
$`Mixed effects model.data`
    District        Study_1           y_1                v_1         
 Min.   : 11.0   Min.   : 1.00   Min.   :-0.52000   Min.   :0.00100  
 1st Qu.: 22.5   1st Qu.:10.50   1st Qu.:-0.12500   1st Qu.:0.01450  
 Median : 58.0   Median :20.00   Median : 0.12000   Median :0.02000  
 Mean   :107.5   Mean   :24.64   Mean   : 0.07182   Mean   :0.03255  
 3rd Qu.: 88.5   3rd Qu.:38.00   3rd Qu.: 0.23000   3rd Qu.:0.02700  
 Max.   :644.0   Max.   :53.00   Max.   : 0.50000   Max.   :0.11800  
                                                                     
     Year_1           Study_2           y_2               v_2         
 Min.   :-13.554   Min.   : 2.00   Min.   :-0.2600   Min.   :0.00100  
 1st Qu.: -7.054   1st Qu.:11.50   1st Qu.:-0.0250   1st Qu.:0.00900  
 Median :  5.446   Median :21.00   Median : 0.3800   Median :0.01400  
 Mean   :  1.083   Mean   :25.64   Mean   : 0.3173   Mean   :0.03091  
 3rd Qu.:  7.446   3rd Qu.:39.00   3rd Qu.: 0.6550   3rd Qu.:0.03700  
 Max.   : 10.446   Max.   :54.00   Max.   : 0.9800   Max.   :0.11800  
                                                                      
     Year_2           Study_3           y_3               v_3         
 Min.   :-13.554   Min.   : 3.00   Min.   :-0.0300   Min.   :0.00100  
 1st Qu.: -7.054   1st Qu.:12.50   1st Qu.: 0.0200   1st Qu.:0.00750  
 Median :  5.446   Median :22.00   Median : 0.1900   Median :0.01200  
 Mean   :  1.083   Mean   :26.64   Mean   : 0.2409   Mean   :0.02882  
 3rd Qu.:  7.446   3rd Qu.:40.00   3rd Qu.: 0.2600   3rd Qu.:0.02450  
 Max.   : 10.446   Max.   :55.00   Max.   : 1.1900   Max.   :0.14400  
                                                                      
     Year_3            Study_4           y_4                v_4           
 Min.   :-13.5536   Min.   : 4.00   Min.   :-0.30000   Min.   :0.000e+00  
 1st Qu.: -7.0536   1st Qu.:15.00   1st Qu.:-0.06000   1st Qu.:0.000e+00  
 Median :  4.4464   Median :23.00   Median : 0.00000   Median :0.000e+00  
 Mean   :  0.9919   Mean   :28.67   Mean   : 0.05778   Mean   :1.818e+09  
 3rd Qu.:  7.4464   3rd Qu.:45.00   3rd Qu.: 0.27000   3rd Qu.:0.000e+00  
 Max.   : 10.4464   Max.   :56.00   Max.   : 0.60000   Max.   :1.000e+10  
                    NA's   :2       NA's   :2                             
     Year_4             Study_5          y_5              v_5           
 Min.   :-13.55357   Min.   :24.0   Min.   :-0.340   Min.   :0.000e+00  
 1st Qu.: -7.05357   1st Qu.:34.5   1st Qu.:-0.115   1st Qu.:0.000e+00  
 Median :  0.00000   Median :42.0   Median :-0.035   Median :1.000e+10  
 Mean   : -0.08929   Mean   :40.0   Mean   :-0.090   Mean   :6.364e+09  
 3rd Qu.:  7.44643   3rd Qu.:47.5   3rd Qu.:-0.010   3rd Qu.:1.000e+10  
 Max.   : 10.44643   Max.   :52.0   Max.   : 0.050   Max.   :1.000e+10  
                     NA's   :7      NA's   :7                           
     Year_5           Study_6        y_6                v_6           
 Min.   :-13.554   Min.   :25   Min.   :-0.30000   Min.   :0.000e+00  
 1st Qu.:  0.000   1st Qu.:32   1st Qu.:-0.15000   1st Qu.:5.000e+09  
 Median :  0.000   Median :39   Median : 0.00000   Median :1.000e+10  
 Mean   :  1.344   Mean   :37   Mean   :-0.07667   Mean   :7.273e+09  
 3rd Qu.:  3.723   3rd Qu.:43   3rd Qu.: 0.03500   3rd Qu.:1.000e+10  
 Max.   : 10.446   Max.   :47   Max.   : 0.07000   Max.   :1.000e+10  
                   NA's   :8    NA's   :8                             
     Year_6            Study_7          y_7             v_7           
 Min.   :-13.5536   Min.   :26.0   Min.   :0.010   Min.   :0.000e+00  
 1st Qu.:  0.0000   1st Qu.:29.5   1st Qu.:0.025   1st Qu.:1.000e+10  
 Median :  0.0000   Median :33.0   Median :0.040   Median :1.000e+10  
 Mean   :  0.3945   Mean   :33.0   Mean   :0.040   Mean   :8.182e+09  
 3rd Qu.:  0.0000   3rd Qu.:36.5   3rd Qu.:0.055   3rd Qu.:1.000e+10  
 Max.   : 10.4464   Max.   :40.0   Max.   :0.070   Max.   :1.000e+10  
                    NA's   :9      NA's   :9                          
     Year_7            Study_8          y_8              v_8           
 Min.   :-13.5536   Min.   :27.0   Min.   :-0.100   Min.   :0.000e+00  
 1st Qu.:  0.0000   1st Qu.:30.5   1st Qu.:-0.075   1st Qu.:1.000e+10  
 Median :  0.0000   Median :34.0   Median :-0.050   Median :1.000e+10  
 Mean   : -0.5552   Mean   :34.0   Mean   :-0.050   Mean   :8.182e+09  
 3rd Qu.:  0.0000   3rd Qu.:37.5   3rd Qu.:-0.025   3rd Qu.:1.000e+10  
 Max.   :  7.4464   Max.   :41.0   Max.   : 0.000   Max.   :1.000e+10  
                    NA's   :9      NA's   :9                           
     Year_8            Study_9        y_9            v_9           
 Min.   :-13.5536   Min.   :28   Min.   :0.05   Min.   :0.000e+00  
 1st Qu.:  0.0000   1st Qu.:28   1st Qu.:0.05   1st Qu.:1.000e+10  
 Median :  0.0000   Median :28   Median :0.05   Median :1.000e+10  
 Mean   : -0.5552   Mean   :28   Mean   :0.05   Mean   :9.091e+09  
 3rd Qu.:  0.0000   3rd Qu.:28   3rd Qu.:0.05   3rd Qu.:1.000e+10  
 Max.   :  7.4464   Max.   :28   Max.   :0.05   Max.   :1.000e+10  
                    NA's   :10   NA's   :10                        
     Year_9           Study_10       y_10            v_10          
 Min.   :-13.554   Min.   :29   Min.   :-0.08   Min.   :0.000e+00  
 1st Qu.:  0.000   1st Qu.:29   1st Qu.:-0.08   1st Qu.:1.000e+10  
 Median :  0.000   Median :29   Median :-0.08   Median :1.000e+10  
 Mean   : -1.232   Mean   :29   Mean   :-0.08   Mean   :9.091e+09  
 3rd Qu.:  0.000   3rd Qu.:29   3rd Qu.:-0.08   3rd Qu.:1.000e+10  
 Max.   :  0.000   Max.   :29   Max.   :-0.08   Max.   :1.000e+10  
                   NA's   :10   NA's   :10                         
    Year_10           Study_11       y_11            v_11          
 Min.   :-13.554   Min.   :30   Min.   :-0.09   Min.   :0.000e+00  
 1st Qu.:  0.000   1st Qu.:30   1st Qu.:-0.09   1st Qu.:1.000e+10  
 Median :  0.000   Median :30   Median :-0.09   Median :1.000e+10  
 Mean   : -1.232   Mean   :30   Mean   :-0.09   Mean   :9.091e+09  
 3rd Qu.:  0.000   3rd Qu.:30   3rd Qu.:-0.09   3rd Qu.:1.000e+10  
 Max.   :  0.000   Max.   :30   Max.   :-0.09   Max.   :1.000e+10  
                   NA's   :10   NA's   :10                         
    Year_11       
 Min.   :-13.554  
 1st Qu.:  0.000  
 Median :  0.000  
 Mean   : -1.232  
 3rd Qu.:  0.000  
 Max.   :  0.000  
                  

free parameters:
    name matrix row col    Estimate   Std.Error lbound ubound
1  beta0  Beta0   1   1 0.178026735 0.080521934              
2  beta1  Beta1   1   1 0.005073714 0.008526627              
3 tau2_2   Tau2   1   1 0.032939012 0.011162043              
4 tau2_3   Tau3   1   1 0.056462856 0.030032981              

observed statistics:  56 
estimated parameters:  4 
degrees of freedom:  52 
-2 log likelihood:  16.43629 
saturated -2 log likelihood:  NA 
number of observations:  11 
chi-square:  NA 
p:  NA 
Information Criteria: 
     df Penalty Parameters Penalty Sample-Size Adjusted
AIC:  -87.56371           24.43629                   NA
BIC: -108.25427           26.02787             13.98387
CFI: NA 
TLI: NA 
RMSEA:  NA 
timestamp: 2013-11-12 22:55:47 
frontend time: 0.1024585 secs 
backend time: 0.02352667 secs 
independent submodels time: 3.767014e-05 secs 
wall clock time: 0.1260228 secs 
cpu time: 0.1260228 secs 
openmx version number: 1.3.2-2301

References

Bornmann, L., Mutz, R., & Daniel, H.-D. (2007). Gender differences in grant peer review: A meta-analysis. Journal of Informetrics, 1(3), 226–238.

Cheung, M.W.-L. (2009). Constructing approximate confidence intervals for parameters with structural equation models. Structural Equation Modeling, 16(2), 267-294.

Cheung, M.W.-L. (2012a). metaSEM: Meta-analysis using structural equation modeling. Retrieved from http://courses.nus.edu.sg/course/psycwlm/Internet/metaSEM/

Cheung, M.W.-L. (2012b). metaSEM: An R package for meta-analysis using strutural equation modeling. Manuscript submitted for publication. Retrieved from http://courses.nus.edu.sg/course/psycwlm/internet/metaSEM/metasem.pdf

Cheung, M.W.-L. (2013, July 8). Modeling dependent effect sizes with three-level meta-analyses: A structural equation modeling approach. Psychological Methods. Advance online publication.

Cooper, H., Valentine, J.C., Charlton, K., & Melson, A. (2003). The effects of modified school calendars on student achievement and on school and community attitudes. Review of Educational Research, 73(1), 1 –52.

Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61–76.

Marsh, H.W., Bornmann, L., Mutz, R., Daniel, H.-D., & O’Mara, A. (2009). Gender effects in the peer reviews of grant proposals: A comprehensive meta-analysis comparing traditional and multilevel approaches. Review of Educational Research, 79(3), 1290–1326.

Neale, M.C., & Miller, M.B. (1997). The use of likelihood-based confidence intervals in genetic models. Behavior Genetics, 27(2), 113–120.

Author: Mike W.-L. Cheung

Created: 2013-11-12 Tue 22:55

Emacs 24.3.1 (Org mode 8.2.1)

Validate