Introduction

Univariate meta-analysis

Random-effects model

## Load the library
library(metaSEM)

## Try to use more than one cores
mxOption(NULL, 'Number of Threads', parallel::detectCores()-2)

## Show the first few studies of the data set
head(Becker83)
  study    di   vi percentage items
1     1 -0.33 0.03         25     2
2     2  0.07 0.03         25     2
3     3 -0.30 0.02         50     2
4     4  0.35 0.02        100    38
5     5  0.69 0.07        100    30
6     6  0.81 0.22        100    45
## Random-effects meta-analysis with ML
summary( meta(y=di, v=vi, data=Becker83) )

Call:
meta(y = di, v = vi, data = Becker83)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate Std.Error    lbound    ubound z value Pr(>|z|)
Intercept1  0.174734  0.113378 -0.047482  0.396950  1.5412   0.1233
Tau2_1_1    0.077376  0.054108 -0.028674  0.183426  1.4300   0.1527

Q statistic on the homogeneity of effect sizes: 30.64949
Degrees of freedom of the Q statistic: 9
P value of the Q statistic: 0.0003399239

Heterogeneity indices (based on the estimated Tau2):
                             Estimate
Intercept1: I2 (Q statistic)   0.6718

Number of studies (or clusters): 10
Number of observed statistics: 10
Number of estimated parameters: 2
Degrees of freedom: 8
-2 log likelihood: 7.928307 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

Fixed-effects model

## Fixed-effects meta-analysis by fixing the heterogeneity variance at 0
summary( meta(y=di, v=vi, data=Becker83, RE.constraints=0) )

Call:
meta(y = di, v = vi, data = Becker83, RE.constraints = 0)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate Std.Error    lbound    ubound z value Pr(>|z|)  
Intercept1  0.100640  0.060510 -0.017957  0.219237  1.6632  0.09627 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 30.64949
Degrees of freedom of the Q statistic: 9
P value of the Q statistic: 0.0003399239

Heterogeneity indices (based on the estimated Tau2):
                             Estimate
Intercept1: I2 (Q statistic)        0

Number of studies (or clusters): 10
Number of observed statistics: 10
Number of estimated parameters: 1
Degrees of freedom: 9
-2 log likelihood: 17.86043 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

Mixed-effects model

## Mixed-effects meta-analysis with "log(items)" as the predictor
summary( meta(y=di, v=vi, x=log(items), data=Becker83) ) 

Call:
meta(y = di, v = vi, x = log(items), data = Becker83)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
              Estimate   Std.Error      lbound      ubound z value
Intercept1 -3.2015e-01  1.0981e-01 -5.3539e-01 -1.0492e-01 -2.9154
Slope1_1    2.1088e-01  4.5084e-02  1.2251e-01  2.9924e-01  4.6774
Tau2_1_1    1.0000e-10  2.0095e-02 -3.9386e-02  3.9386e-02  0.0000
            Pr(>|z|)    
Intercept1  0.003552 ** 
Slope1_1   2.905e-06 ***
Tau2_1_1    1.000000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 30.64949
Degrees of freedom of the Q statistic: 9
P value of the Q statistic: 0.0003399239

Explained variances (R2):
                           y1
Tau2 (no predictor)    0.0774
Tau2 (with predictors) 0.0000
R2                     1.0000

Number of studies (or clusters): 10
Number of observed statistics: 10
Number of estimated parameters: 3
Degrees of freedom: 7
-2 log likelihood: -4.208024 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

Multivariate meta-analysis

Random-effects model

## Show the data set
Berkey98
  trial pub_year no_of_patients   PD    AL var_PD cov_PD_AL var_AL
1     1     1983             14 0.47 -0.32 0.0075    0.0030 0.0077
2     2     1982             15 0.20 -0.60 0.0057    0.0009 0.0008
3     3     1979             78 0.40 -0.12 0.0021    0.0007 0.0014
4     4     1987             89 0.26 -0.31 0.0029    0.0009 0.0015
5     5     1988             16 0.56 -0.39 0.0148    0.0072 0.0304
## Multivariate meta-analysis with a random-effects model
mult1 <- meta(y=cbind(PD, AL), v=cbind(var_PD, cov_PD_AL, var_AL), data=Berkey98) 

summary(mult1)

Call:
meta(y = cbind(PD, AL), v = cbind(var_PD, cov_PD_AL, var_AL), 
    data = Berkey98)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
             Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
Intercept1  0.3448392  0.0536312  0.2397239  0.4499544  6.4298 1.278e-10
Intercept2 -0.3379381  0.0812480 -0.4971812 -0.1786951 -4.1593 3.192e-05
Tau2_1_1    0.0070020  0.0090497 -0.0107351  0.0247391  0.7737    0.4391
Tau2_2_1    0.0094607  0.0099698 -0.0100797  0.0290010  0.9489    0.3427
Tau2_2_2    0.0261445  0.0177409 -0.0086270  0.0609161  1.4737    0.1406
              
Intercept1 ***
Intercept2 ***
Tau2_1_1      
Tau2_2_1      
Tau2_2_2      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 128.2267
Degrees of freedom of the Q statistic: 8
P value of the Q statistic: 0

Heterogeneity indices (based on the estimated Tau2):
                             Estimate
Intercept1: I2 (Q statistic)   0.6021
Intercept2: I2 (Q statistic)   0.9250

Number of studies (or clusters): 5
Number of observed statistics: 10
Number of estimated parameters: 5
Degrees of freedom: 5
-2 log likelihood: -11.68131 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
## Plot the effect sizes
plot(mult1)

## Plot the effect sizes with the forest plots
## Load the library for forest plots 
library("metafor")

## Create extra panels for the forest plots
plot(mult1, diag.panel=TRUE, main="Multivariate meta-analysis",
axis.label=c("PD", "AL"))

## Forest plot for PD
forest( rma(yi=PD, vi=var_PD, data=Berkey98) )
title("Forest plot of PD")

## Forest plot for AL
forest( rma(yi=AL, vi=var_AL, data=Berkey98) )
title("Forest plot of AL")

Fixed-effects model

## Fixed-effects meta-analysis by fixiing the heterogeneity variance component at 
## a 2x2 matrix of 0.
summary( meta(y=cbind(PD, AL), v=cbind(var_PD, cov_PD_AL, var_AL), data=Berkey98,
              RE.constraints=matrix(0, nrow=2, ncol=2)) )

Call:
meta(y = cbind(PD, AL), v = cbind(var_PD, cov_PD_AL, var_AL), 
    data = Berkey98, RE.constraints = matrix(0, nrow = 2, ncol = 2))

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate Std.Error    lbound    ubound z value  Pr(>|z|)    
Intercept1  0.307219  0.028575  0.251212  0.363225  10.751 < 2.2e-16 ***
Intercept2 -0.394377  0.018649 -0.430929 -0.357825 -21.147 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 128.2267
Degrees of freedom of the Q statistic: 8
P value of the Q statistic: 0

Heterogeneity indices (based on the estimated Tau2):
                             Estimate
Intercept1: I2 (Q statistic)        0
Intercept2: I2 (Q statistic)        0

Number of studies (or clusters): 5
Number of observed statistics: 10
Number of estimated parameters: 2
Degrees of freedom: 8
-2 log likelihood: 90.88326 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

Mixed-effects model

## Multivariate meta-analysis with "publication year-1979" as a predictor 
summary( meta(y=cbind(PD, AL), v=cbind(var_PD, cov_PD_AL, var_AL), data=Berkey98,
              x=scale(pub_year, center=1979)) )

Call:
meta(y = cbind(PD, AL), v = cbind(var_PD, cov_PD_AL, var_AL), 
    x = scale(pub_year, center = 1979), data = Berkey98)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
             Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
Intercept1  0.3440001  0.0857659  0.1759020  0.5120982  4.0109 6.048e-05
Intercept2 -0.2918175  0.1312796 -0.5491208 -0.0345141 -2.2229   0.02622
Slope1_1    0.0063540  0.1078235 -0.2049762  0.2176842  0.0589   0.95301
Slope2_1   -0.0705888  0.1620966 -0.3882922  0.2471146 -0.4355   0.66322
Tau2_1_1    0.0080405  0.0101206 -0.0117955  0.0278766  0.7945   0.42692
Tau2_2_1    0.0093413  0.0105515 -0.0113392  0.0300218  0.8853   0.37599
Tau2_2_2    0.0250135  0.0170788 -0.0084603  0.0584873  1.4646   0.14303
              
Intercept1 ***
Intercept2 *  
Slope1_1      
Slope2_1      
Tau2_1_1      
Tau2_2_1      
Tau2_2_2      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 128.2267
Degrees of freedom of the Q statistic: 8
P value of the Q statistic: 0

Explained variances (R2):
                              y1     y2
Tau2 (no predictor)    0.0070020 0.0261
Tau2 (with predictors) 0.0080405 0.0250
R2                     0.0000000 0.0433

Number of studies (or clusters): 5
Number of observed statistics: 10
Number of estimated parameters: 7
Degrees of freedom: 3
-2 log likelihood: -12.00859 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

Three-level meta-analysis

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.

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)
  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

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.

  • 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) )

Call:
meta(y = y, v = v, data = Cooper03)

95% confidence intervals: z statistic approximation (robust=FALSE)
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 the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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) )

Call:
meta(y = y, v = v, data = Cooper03, RE.constraints = 0)

95% confidence intervals: z statistic approximation (robust=FALSE)
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 the 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: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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) )

Call:
meta(y = y, v = v, x = scale(Year, scale = FALSE), data = Cooper03)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
             Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
Intercept1  0.1259126  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
              
Intercept1 ** 
Slope1_1      
Tau2_1_1   ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)

Three-level meta-analysis

  • 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) )

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

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Intercept  0.1844554  0.0805411  0.0265977  0.3423131  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 the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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) )

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

95% confidence intervals: z statistic approximation (robust=FALSE)
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 the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • Model comparisons Many research hypotheses involve model comparisons among nested models. anova(), a generic function for comparing nested models, may be used to conduct a likelihood ratio test which is also known as a chi-square difference test.

  • 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 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
  • 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)
           base                comparison ep minus2LL df       AIC
1 3 level model                      <NA>  3 16.78987 53 -89.21013
2 3 level model Equal tau2 at both levels  2 17.47697 54 -90.52303
     diffLL diffdf         p
1        NA     NA        NA
2 0.6870959      1 0.4071539
  • 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.
    • An 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.
    • An 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.
    • An 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.
## LBCI for random-effects model
summary( meta3(y=y, v=v, cluster=District, data=Cooper03, 
               I2=c("I2q", "ICC"), intervals.type="LB") ) 

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.012025 0.358179      NA       NA
Tau2_2    0.032865        NA 0.016331 0.060489      NA       NA
Tau2_3    0.057738        NA 0.019810 0.172714      NA       NA

Q statistic on the 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.12754  0.34396 0.6565
ICC_2 (tau^2/(tau^2+tau^3))   0.13123  0.36273 0.7005
I2_3 (Typical v: Q statistic) 0.28173  0.60429 0.8451
ICC_3 (tau^3/(tau^2+tau^3))   0.29948  0.63727 0.8688

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" or "1": The optimization is considered fine.
Other values may indicate problems.)
## LBCI for mixed-effects model
summary( meta3(y=y, v=v, cluster=District, x=scale(Year, scale=FALSE), 
               data=Cooper03, intervals.type="LB") ) 

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.0056477  0.3512066      NA       NA
Slope_1    0.0050737        NA -0.0128322  0.0237932      NA       NA
Tau2_2     0.0329390        NA  0.0163748  0.0591198      NA       NA
Tau2_3     0.0564628        NA  0.0192365  0.1429149      NA       NA

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • Restricted maximum likelihood estimation
    • REML may also be used in the 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) )

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 status1: 0 ("0" or "1": The optimization is considered fine.
Other values may 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") )

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 status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
anova(reml1, reml0)
                          base comparison ep  minus2LL df AIC   diffLL
1 Variance component with REML       <NA>  2 -81.14044 -2  NA       NA
2 Variance component with REML Equal Tau2  1 -80.13710 -1  NA 1.003336
  diffdf         p
1     NA        NA
2      1 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) )

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 status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

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.

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

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") )

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

95% confidence intervals: z statistic approximation (robust=FALSE)
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 the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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)
           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
  • 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)
           base           comparison ep minus2LL df       AIC   diffLL
1 3 level model                 <NA>  3 25.80256 63 -100.1974       NA
2 3 level model tau2_2 equals tau2_3  2 27.16166 64 -100.8383 1.359103
  diffdf        p
1     NA       NA
2      1 0.243693

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 of 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)) 

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

95% confidence intervals: z statistic approximation (robust=FALSE)
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
             
Intercept    
Slope_1   ***
Tau2_2       
Tau2_3       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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
Grant <- ifelse(Bornmann07$Type=="Grant", yes=1, no=0)
Fellowship <- ifelse(Bornmann07$Type=="Fellowship", yes=1, no=0)
  
Model1b <- meta3(logOR, v, x=cbind(Grant, Fellowship), 
                 cluster=Cluster, data=Bornmann07,
                 intercept.constraints=0, model.name="Model 1")
summary(Model1b)

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = cbind(Grant, Fellowship), 
    data = Bornmann07, intercept.constraints = 0, model.name = "Model 1")

95% confidence intervals: z statistic approximation (robust=FALSE)
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 the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)

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") ) 

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 (robust=FALSE)
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
           
Intercept *
Slope_1    
Slope_2   .
Tau2_2     
Tau2_3     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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

Model 3: Discipline as a covariate

  • There are four categories of 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") )

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 (robust=FALSE)
Coefficients:
             Estimate   Std.Error      lbound      ubound z value Pr(>|z|)
Intercept -0.01474783  0.06389945 -0.13998845  0.11049278 -0.2308  0.81747
Slope_1   -0.00913064  0.09949198 -0.20413134  0.18587005 -0.0918  0.92688
Slope_2   -0.12617957  0.07866274 -0.28035570  0.02799656 -1.6041  0.10870
Slope_3   -0.23695698  0.12123091 -0.47456520  0.00065125 -1.9546  0.05063
Tau2_2     0.00390942  0.00283949 -0.00165587  0.00947471  1.3768  0.16857
Tau2_3     0.00710338  0.00643210 -0.00550331  0.01971006  1.1044  0.26944
           
Intercept  
Slope_1    
Slope_2    
Slope_3   .
Tau2_2     
Tau2_3     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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
1 Model 3          <NA>  6 20.07571 60  -99.92429       NA     NA
2 Model 3 3 level model  3 25.80256 63 -100.19744 5.726842      3
          p
1        NA
2 0.1256832

Model 4: Country as a covariate

  • There are five categories in Country. The United States is used as the reference group in the analysis.
  • The estimated regression coefficients (and their 95% Wald CIs) for CountryAus, CountryCan, CountryEur, and 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") )  

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 (robust=FALSE)
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.0992865 -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 the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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
1 Model 4          <NA>  7 14.18259 59 -103.8174       NA     NA
2 Model 4 3 level model  3 25.80256 63 -100.1974 11.61996      4
           p
1         NA
2 0.02041284

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
Model5 <- meta3(logOR, v, x=cbind(Type2, DisciplinePhy, DisciplineLife, 
                DisciplineSoc), cluster=Cluster, data=Bornmann07,
                model.name="Model 5")

## Rerun to remove error code
Model5 <- rerun(Model5)
summary(Model5)

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

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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

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
Model6 <- meta3(logOR, v, x=cbind(Type2, CountryAus, CountryCan, 
                CountryEur, CountryUK), cluster=Cluster, data=Bornmann07,
                model.name="Model 6") 

## Rerun to remove error code
Model6 <- rerun(Model6)
summary(Model6)

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 (robust=FALSE)
Coefficients:
             Estimate   Std.Error      lbound      ubound z value
Intercept  6.7507e-02  1.8933e-02  3.0399e-02  1.0461e-01  3.5656
Slope_1   -1.5167e-01  4.1113e-02 -2.3225e-01 -7.1092e-02 -3.6892
Slope_2   -6.9580e-02  8.5164e-02 -2.3650e-01  9.7339e-02 -0.8170
Slope_3   -1.4231e-01  7.5204e-02 -2.8970e-01  5.0878e-03 -1.8923
Slope_4   -1.6116e-01  4.0203e-02 -2.3995e-01 -8.2361e-02 -4.0086
Slope_5    9.0419e-03  7.0074e-02 -1.2830e-01  1.4639e-01  0.1290
Tau2_2     2.2976e-03  1.4407e-03 -5.2618e-04  5.1213e-03  1.5947
Tau2_3     1.0000e-10          NA          NA          NA      NA
           Pr(>|z|)    
Intercept 0.0003631 ***
Slope_1   0.0002250 ***
Slope_2   0.4139266    
Slope_3   0.0584497 .  
Slope_4   6.108e-05 ***
Slope_5   0.8973315    
Tau2_2    0.1107693    
Tau2_3           NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • 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

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") )

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 (robust=FALSE)
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.1594775  0.5079114  1.0233  0.30618  
Slope_2    0.0826806  0.1599802 -0.2308749  0.3962360  0.5168  0.60528  
Slope_3   -0.0462265  0.1715774 -0.3825119  0.2900590 -0.2694  0.78761  
Slope_4   -0.0486321  0.1306918 -0.3047834  0.2075192 -0.3721  0.70981  
Slope_5   -0.2169132  0.1915703 -0.5923841  0.1585577 -1.1323  0.25751  
Slope_6   -0.3036578  0.1670721 -0.6311130  0.0237975 -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 the 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.1396974  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" or "1": The optimization is considered fine.
Other values may indicate problems.)

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))

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

Q statistic on the 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: 6 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
Warning in print.summary.meta(x): OpenMx status1 is neither 0 or 1. You are advised to 'rerun' it again.

Handling missing covariates with FIML

When there are missing data in the covariates, data with missing values are excluded from 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.

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))

Call:
meta3(y = logOR, v = v, cluster = Cluster, x = Type_MCAR, data = my.df)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
Intercept  0.0044909  0.0362672 -0.0665916  0.0755733  0.1238    0.9015
Slope_1   -0.2182446  0.0554287 -0.3268829 -0.1096063 -3.9374 8.237e-05
Tau2_2     0.0014063  0.0021623 -0.0028317  0.0056443  0.6504    0.5155
Tau2_3     0.0031148  0.0035202 -0.0037846  0.0100143  0.8848    0.3762
             
Intercept    
Slope_1   ***
Tau2_2       
Tau2_3       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 154.2762
Degrees of freedom of the Q statistic: 48
P value of the Q statistic: 4.410916e-13

Explained variances (R2):
                         Level 2 Level 3
Tau2 (no predictor)    0.0011603  0.0185
Tau2 (with predictors) 0.0014063  0.0031
R2                     0.0000000  0.8318

Number of studies (or clusters): 20
Number of observed statistics: 49
Number of estimated parameters: 4
Degrees of freedom: 45
-2 log likelihood: 10.56012 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may 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))

Call:
meta3X(y = logOR, v = v, cluster = Cluster, x2 = Type_MCAR, data = my.df)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)
Intercept -0.0024343  0.0360701 -0.0731303  0.0682618 -0.0675 0.9461939
SlopeX2_1 -0.2086677  0.0545138 -0.3155128 -0.1018226 -3.8278 0.0001293
Tau2_2     0.0016732  0.0022114 -0.0026610  0.0060075  0.7567 0.4492584
Tau2_3     0.0035540  0.0035810 -0.0034646  0.0105726  0.9925 0.3209675
             
Intercept    
SlopeX2_1 ***
Tau2_2       
Tau2_3       
---
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.0016732  0.0036
R2                     0.5592670  0.7486

Number of studies (or clusters): 21
Number of observed statistics: 115
Number of estimated parameters: 7
Degrees of freedom: 108
-2 log likelihood: 56.64328 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

MAR

  • For the case of 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)) 

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

95% confidence intervals: z statistic approximation (robust=FALSE)
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.17573648  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
            
Intercept   
Slope_1   **
Tau2_2      
Tau2_3      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the 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" or "1": The optimization is considered fine.
Other values may 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))

Call:
meta3X(y = logOR, v = v, cluster = Cluster, x2 = Type_MAR, av2 = Year, 
    data = my.df)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Intercept -0.0264058  0.0572057 -0.1385269  0.0857154 -0.4616 0.644373   
SlopeX2_1 -0.2003999  0.0691094 -0.3358519 -0.0649479 -2.8997 0.003735 **
Tau2_2     0.0029970  0.0022371 -0.0013877  0.0073817  1.3396 0.180359   
Tau2_3     0.0030212  0.0032463 -0.0033415  0.0093839  0.9307 0.352034   
---
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.3913243  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" or "1": The optimization is considered fine.
Other values may indicate problems.)

Two-stage structural equation modeling (TSSEM)

Inspect the data: Digman (1997)

The correlation matrices and the sample sizes were stored in Digman97$data and Digman97$n, respectively. We may list the first few cases of the data by using the head() command.

#### Load the metaSEM library for TSSEM
library(metaSEM)
  
#### Inspect the data for inspection
head(Digman97$data)
$`Digman 1 (1994)`
       A     C   ES     E    I
A   1.00  0.62 0.41 -0.48 0.00
C   0.62  1.00 0.59 -0.10 0.35
ES  0.41  0.59 1.00  0.27 0.41
E  -0.48 -0.10 0.27  1.00 0.37
I   0.00  0.35 0.41  0.37 1.00

$`Digman 2 (1994)`
       A    C   ES     E     I
A   1.00 0.39 0.53 -0.30 -0.05
C   0.39 1.00 0.59  0.07  0.44
ES  0.53 0.59 1.00  0.09  0.22
E  -0.30 0.07 0.09  1.00  0.45
I  -0.05 0.44 0.22  0.45  1.00

$`Digman 3 (1963c)`
      A     C   ES     E    I
A  1.00  0.65 0.35  0.25 0.14
C  0.65  1.00 0.37 -0.10 0.33
ES 0.35  0.37 1.00  0.24 0.41
E  0.25 -0.10 0.24  1.00 0.41
I  0.14  0.33 0.41  0.41 1.00

$`Digman & Takemoto-Chock (1981b)`
       A     C   ES     E     I
A   1.00  0.65 0.70 -0.26 -0.03
C   0.65  1.00 0.71 -0.16  0.24
ES  0.70  0.71 1.00  0.01  0.11
E  -0.26 -0.16 0.01  1.00  0.66
I  -0.03  0.24 0.11  0.66  1.00

$`Graziano & Ward (1992)`
      A    C   ES    E    I
A  1.00 0.64 0.35 0.29 0.22
C  0.64 1.00 0.27 0.16 0.22
ES 0.35 0.27 1.00 0.32 0.36
E  0.29 0.16 0.32 1.00 0.53
I  0.22 0.22 0.36 0.53 1.00

$`Yik & Bond (1993)`
      A    C   ES    E    I
A  1.00 0.66 0.57 0.35 0.38
C  0.66 1.00 0.45 0.20 0.31
ES 0.57 0.45 1.00 0.49 0.31
E  0.35 0.20 0.49 1.00 0.59
I  0.38 0.31 0.31 0.59 1.00
head(Digman97$n)
[1] 102 149 334 162  91 656

Fixed-effects TSSEM

Stage 1 analysis

To conduct a fixed-effects TSSEM, we may specify method=FEM argument (the default method) in calling the tssem1() function. The results are stored in an object named fixed1. It can be displayed by the summary() command. The \(\chi^2(130, N=4,496) = 1,499.73, p < .001\), CFI=0.6825, RMSEA=0.1812 and SRMR=0.1750. Based on the test statistic and the goodness-of-fit indices, the assumption of homogeneity of correlation matrices was rejected.

## Fixed-effects model: Stage 1 analysis
fixed1 <- tssem1(Cov=Digman97$data, n=Digman97$n, method="FEM")
summary(fixed1)

Call:
tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
    cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
    run = run)

Coefficients:
       Estimate Std.Error z value  Pr(>|z|)    
S[1,2] 0.363278  0.013368 27.1759 < 2.2e-16 ***
S[1,3] 0.390375  0.012880 30.3075 < 2.2e-16 ***
S[1,4] 0.103572  0.015048  6.8830 5.862e-12 ***
S[1,5] 0.092286  0.015047  6.1331 8.621e-10 ***
S[2,3] 0.416070  0.012519 33.2343 < 2.2e-16 ***
S[2,4] 0.135148  0.014776  9.1464 < 2.2e-16 ***
S[2,5] 0.141431  0.014866  9.5135 < 2.2e-16 ***
S[3,4] 0.244459  0.014153 17.2723 < 2.2e-16 ***
S[3,5] 0.138339  0.014834  9.3259 < 2.2e-16 ***
S[4,5] 0.424566  0.012376 34.3070 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Goodness-of-fit indices:
                                     Value
Sample size                      4496.0000
Chi-square of target model       1505.4443
DF of target model                130.0000
p value of target model             0.0000
Chi-square of independence model 4471.4242
DF of independence model          140.0000
RMSEA                               0.1815
RMSEA lower 95% CI                  0.1736
RMSEA upper 95% CI                  0.1901
SRMR                                0.1621
TLI                                 0.6580
CFI                                 0.6824
AIC                              1245.4443
BIC                               412.0217
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

The pooled correlation matrix (the parameter estimates) can be extracted by the use of the coef() command.

coef(fixed1)
            A         C        ES         E          I
A  1.00000000 0.3632782 0.3903748 0.1035716 0.09228557
C  0.36327824 1.0000000 0.4160695 0.1351477 0.14143058
ES 0.39037483 0.4160695 1.0000000 0.2444593 0.13833895
E  0.10357155 0.1351477 0.2444593 1.0000000 0.42456626
I  0.09228557 0.1414306 0.1383390 0.4245663 1.00000000

Stage 2 analysis

As an illustration, I continued to fit the structural model based on the pooled correlation matrix. We may specify the two-factor model with the RAM formulation

## Factor covariance among latent factors
Phi <- matrix(c(1,"0.3*cor","0.3*cor",1), ncol=2, nrow=2)

## Error covariance matrix
Psi <- Diag(c("0.2*e1","0.2*e2","0.2*e3","0.2*e4","0.2*e5"))

## S matrix
S1 <- bdiagMat(list(Psi, Phi))

## This step is not necessary, but it is useful for inspecting the model.
dimnames(S1)[[1]] <- dimnames(S1)[[2]] <- c("A","C","ES","E","I","Alpha","Beta") 

## Display S1
S1
      A        C        ES       E        I        Alpha     Beta     
A     "0.2*e1" "0"      "0"      "0"      "0"      "0"       "0"      
C     "0"      "0.2*e2" "0"      "0"      "0"      "0"       "0"      
ES    "0"      "0"      "0.2*e3" "0"      "0"      "0"       "0"      
E     "0"      "0"      "0"      "0.2*e4" "0"      "0"       "0"      
I     "0"      "0"      "0"      "0"      "0.2*e5" "0"       "0"      
Alpha "0"      "0"      "0"      "0"      "0"      "1"       "0.3*cor"
Beta  "0"      "0"      "0"      "0"      "0"      "0.3*cor" "1"      
## A matrix
Lambda <-
matrix(c(".3*Alpha_A",".3*Alpha_C",".3*Alpha_ES",rep(0,5),".3*Beta_E",".3*Beta_I"),
       ncol=2, nrow=5)
A1 <- rbind( cbind(matrix(0,ncol=5,nrow=5), Lambda),
             matrix(0, ncol=7, nrow=2) )

## This step is not necessary, but it is useful for inspecting the model.
dimnames(A1)[[1]] <- dimnames(A1)[[2]] <- c("A","C","ES","E","I","Alpha","Beta") 

## Display A1
A1
      A   C   ES  E   I   Alpha         Beta       
A     "0" "0" "0" "0" "0" ".3*Alpha_A"  "0"        
C     "0" "0" "0" "0" "0" ".3*Alpha_C"  "0"        
ES    "0" "0" "0" "0" "0" ".3*Alpha_ES" "0"        
E     "0" "0" "0" "0" "0" "0"           ".3*Beta_E"
I     "0" "0" "0" "0" "0" "0"           ".3*Beta_I"
Alpha "0" "0" "0" "0" "0" "0"           "0"        
Beta  "0" "0" "0" "0" "0" "0"           "0"        
## F matrix to select the observed variables
F1 <- create.Fmatrix(c(1,1,1,1,1,0,0), as.mxMatrix=FALSE)

## Display F1
F1
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    0    0    0    0    0
[2,]    0    1    0    0    0    0    0
[3,]    0    0    1    0    0    0    0
[4,]    0    0    0    1    0    0    0
[5,]    0    0    0    0    1    0    0
  • We may also specify the model using lavann’s syntax.
model1 <- "## Factor loadings
           Alpha=~A+C+ES
           Beta=~E+I
           ## Factor correlation
           Alpha~~Beta"

RAM1 <- lavaan2RAM(model1, obs.variables=c("A","C","ES","E","I"),
                   A.notation="on", S.notation="with")
RAM1
$A
      A   C   ES  E   I   Alpha         Beta       
A     "0" "0" "0" "0" "0" "0*AonAlpha"  "0"        
C     "0" "0" "0" "0" "0" "0*ConAlpha"  "0"        
ES    "0" "0" "0" "0" "0" "0*ESonAlpha" "0"        
E     "0" "0" "0" "0" "0" "0"           "0*EonBeta"
I     "0" "0" "0" "0" "0" "0"           "0*IonBeta"
Alpha "0" "0" "0" "0" "0" "0"           "0"        
Beta  "0" "0" "0" "0" "0" "0"           "0"        

$S
      A          C          ES           E          I         
A     "0*AwithA" "0"        "0"          "0"        "0"       
C     "0"        "0*CwithC" "0"          "0"        "0"       
ES    "0"        "0"        "0*ESwithES" "0"        "0"       
E     "0"        "0"        "0"          "0*EwithE" "0"       
I     "0"        "0"        "0"          "0"        "0*IwithI"
Alpha "0"        "0"        "0"          "0"        "0"       
Beta  "0"        "0"        "0"          "0"        "0"       
      Alpha             Beta             
A     "0"               "0"              
C     "0"               "0"              
ES    "0"               "0"              
E     "0"               "0"              
I     "0"               "0"              
Alpha "1"               "0*AlphawithBeta"
Beta  "0*AlphawithBeta" "1"              

$F
   A C ES E I Alpha Beta
A  1 0  0 0 0     0    0
C  0 1  0 0 0     0    0
ES 0 0  1 0 0     0    0
E  0 0  0 1 0     0    0
I  0 0  0 0 1     0    0

$M
  A C ES E I Alpha Beta
1 0 0  0 0 0     0    0
A1 <- RAM1$A
S1 <- RAM1$S
F1 <- RAM1$F
  • Since we are conducting a correlation structure analysis, the error variances are not free parameters. The chi-square test statistic of the Stage 2 analysis was \(\chi^2(4, N=4,496) = 65.45, p < .001\), CFI=0.9802, RMSEA=0.0585 and SRMR=0.0284.
fixed2 <- tssem2(fixed1, Amatrix=A1, Smatrix=S1, Fmatrix=F1, 
                 diag.constraints=FALSE, intervals="LB",
                 model.name="TSSEM2 Digman97")
summary(fixed2)

Call:
wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
    n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
    Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
    intervals.type = intervals.type, mx.algebras = mx.algebras, 
    model.name = model.name, suppressWarnings = suppressWarnings, 
    silent = silent, run = run)

95% confidence intervals: Likelihood-based statistic
Coefficients:
              Estimate Std.Error  lbound  ubound z value Pr(>|z|)
AonAlpha       0.56280        NA 0.53270 0.59302      NA       NA
ConAlpha       0.60522        NA 0.57525 0.63535      NA       NA
ESonAlpha      0.71920        NA 0.68876 0.75032      NA       NA
EonBeta        0.78141        NA 0.71872 0.85497      NA       NA
IonBeta        0.55137        NA 0.50000 0.60270      NA       NA
AlphawithBeta  0.36268        NA 0.31859 0.40646      NA       NA

Goodness-of-fit indices:
                                               Value
Sample size                                4496.0000
Chi-square of target model                   65.4526
DF of target model                            4.0000
p value of target model                       0.0000
Number of constraints imposed on "Smatrix"    0.0000
DF manually adjusted                          0.0000
Chi-square of independence model           3112.6867
DF of independence model                     10.0000
RMSEA                                         0.0585
RMSEA lower 95% CI                            0.0465
RMSEA upper 95% CI                            0.0713
SRMR                                          0.0284
TLI                                           0.9505
CFI                                           0.9802
AIC                                          57.4526
BIC                                          31.8088
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
plot(fixed2)

Fixed-effects TSSEM with cluster

Stage 1 analysis

There are 4 types of sample characteristics in the original cluster. We may group them into either younger or older sample.

#### Display the frequencies of "cluster"
table(Digman97$cluster)

  Adolescents      Children Mature adults  Young adults 
            1             4             6             3 
#### Fixed-effects TSSEM with several clusters
#### Create a variable for different cluster
#### Younger participants: Children and Adolescents
#### Older participants: others
clusters <- ifelse(Digman97$cluster %in% c("Children","Adolescents"),
                   yes="Younger participants", no="Older participants")
  
#### Show the clusters
clusters
 [1] "Younger participants" "Younger participants" "Younger participants"
 [4] "Younger participants" "Younger participants" "Older participants"  
 [7] "Older participants"   "Older participants"   "Older participants"  
[10] "Older participants"   "Older participants"   "Older participants"  
[13] "Older participants"   "Older participants"  
  • We may conduct a fixed-effects model by specifying the cluster=clusters argument. Fixed-effects TSSEM will be conducted according to the labels in the clusters. The goodness-of-fit indices of the Stage 1 analysis for the older and younger participants were \(\chi^2(80, N=3,658) = 823.88, p < .001\), CFI=0.7437, RMSEA=0.1513 and SRMR=0.1528, and \(\chi^2(40, N=838) = 344.18, p < .001\), CFI=0.7845, RMSEA=0.2131 and SRMR=0.1508, respectively.
## Example of Fixed-effects TSSEM with several clusters
cluster1 <- tssem1(Digman97$data, Digman97$n, method="FEM", 
                   cluster=clusters)
summary(cluster1)
$`Older participants`

Call:
tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
    model.name = model.name, suppressWarnings = suppressWarnings)

Coefficients:
       Estimate Std.Error z value  Pr(>|z|)    
S[1,2] 0.297471  0.015436 19.2716 < 2.2e-16 ***
S[1,3] 0.370248  0.014532 25.4774 < 2.2e-16 ***
S[1,4] 0.137694  0.016403  8.3944 < 2.2e-16 ***
S[1,5] 0.098061  0.016724  5.8637 4.528e-09 ***
S[2,3] 0.393692  0.014146 27.8307 < 2.2e-16 ***
S[2,4] 0.183045  0.016055 11.4009 < 2.2e-16 ***
S[2,5] 0.092774  0.016643  5.5743 2.485e-08 ***
S[3,4] 0.260753  0.015554 16.7645 < 2.2e-16 ***
S[3,5] 0.096141  0.016573  5.8009 6.596e-09 ***
S[4,5] 0.411756  0.013900 29.6225 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Goodness-of-fit indices:
                                     Value
Sample size                      3658.0000
Chi-square of target model        825.9826
DF of target model                 80.0000
p value of target model             0.0000
Chi-square of independence model 3000.9731
DF of independence model           90.0000
RMSEA                               0.1515
RMSEA lower 95% CI                  0.1424
RMSEA upper 95% CI                  0.1611
SRMR                                0.1459
TLI                                 0.7117
CFI                                 0.7437
AIC                               665.9826
BIC                               169.6088
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

$`Younger participants`

Call:
tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
    model.name = model.name, suppressWarnings = suppressWarnings)

Coefficients:
        Estimate Std.Error z value  Pr(>|z|)    
S[1,2]  0.604330  0.022125 27.3141 < 2.2e-16 ***
S[1,3]  0.465536  0.027493 16.9327 < 2.2e-16 ***
S[1,4] -0.031381  0.035940 -0.8732   0.38258    
S[1,5]  0.061508  0.034547  1.7804   0.07500 .  
S[2,3]  0.501432  0.026348 19.0311 < 2.2e-16 ***
S[2,4] -0.060678  0.034557 -1.7559   0.07911 .  
S[2,5]  0.320987  0.031064 10.3330 < 2.2e-16 ***
S[3,4]  0.175437  0.033675  5.2097 1.891e-07 ***
S[3,5]  0.305149  0.031586  9.6609 < 2.2e-16 ***
S[4,5]  0.478640  0.026883 17.8045 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Goodness-of-fit indices:
                                     Value
Sample size                       838.0000
Chi-square of target model        346.2810
DF of target model                 40.0000
p value of target model             0.0000
Chi-square of independence model 1470.4511
DF of independence model           50.0000
RMSEA                               0.2139
RMSEA lower 95% CI                  0.1939
RMSEA upper 95% CI                  0.2355
SRMR                                0.1411
TLI                                 0.7305
CFI                                 0.7844
AIC                               266.2810
BIC                                77.0402
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • The pooled correlation matrix for each cluster can be extracted by the use of the coef() command.
coef(cluster1)
$`Older participants`
            A          C         ES         E          I
A  1.00000000 0.29747123 0.37024803 0.1376942 0.09806125
C  0.29747123 1.00000000 0.39369157 0.1830450 0.09277411
ES 0.37024803 0.39369157 1.00000000 0.2607530 0.09614072
E  0.13769424 0.18304500 0.26075304 1.0000000 0.41175622
I  0.09806125 0.09277411 0.09614072 0.4117562 1.00000000

$`Younger participants`
             A          C        ES          E          I
A   1.00000000  0.6043300 0.4655359 -0.0313810 0.06150839
C   0.60433002  1.0000000 0.5014319 -0.0606784 0.32098713
ES  0.46553592  0.5014319 1.0000000  0.1754367 0.30514853
E  -0.03138100 -0.0606784 0.1754367  1.0000000 0.47864004
I   0.06150839  0.3209871 0.3051485  0.4786400 1.00000000

Stage 2 analysis

The goodness-of-fit indices of the Stage 2 analysis for the older and younger participants were \(\chi^2(4, N=3,658) = 21.92, p < .001\), CFI=0.9921, RMSEA=0.0350 and SRMR=0.0160, and \(\chi^2(4, N=838) = 144.87, p < .001\), CFI=0.9427, RMSEA=0.2051 and SRMR=0.1051, respectively.

cluster2 <- tssem2(cluster1, Amatrix=A1, Smatrix=S1, Fmatrix=F1, 
                   diag.constraints=FALSE, intervals.type="z")
#### Please note that the estimates for the younger participants are problematic.
summary(cluster2)
$`Older participants`

Call:
wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
    n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
    Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
    intervals.type = intervals.type, mx.algebras = mx.algebras, 
    model.name = model.name, suppressWarnings = suppressWarnings, 
    silent = silent, run = run)

95% confidence intervals: z statistic approximation
Coefficients:
              Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
AonAlpha      0.512656  0.018206 0.476973 0.548340  28.158 < 2.2e-16 ***
ConAlpha      0.549967  0.017945 0.514795 0.585140  30.647 < 2.2e-16 ***
ESonAlpha     0.732174  0.018929 0.695074 0.769273  38.680 < 2.2e-16 ***
EonBeta       0.967136  0.058751 0.851986 1.082287  16.462 < 2.2e-16 ***
IonBeta       0.430649  0.029634 0.372568 0.488730  14.532 < 2.2e-16 ***
AlphawithBeta 0.349236  0.028118 0.294125 0.404346  12.420 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Goodness-of-fit indices:
                                               Value
Sample size                                3658.0000
Chi-square of target model                   21.9954
DF of target model                            4.0000
p value of target model                       0.0002
Number of constraints imposed on "Smatrix"    0.0000
DF manually adjusted                          0.0000
Chi-square of independence model           2273.3342
DF of independence model                     10.0000
RMSEA                                         0.0351
RMSEA lower 95% CI                            0.0217
RMSEA upper 95% CI                            0.0500
SRMR                                          0.0160
TLI                                           0.9801
CFI                                           0.9920
AIC                                          13.9954
BIC                                         -10.8233
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)

$`Younger participants`

Call:
wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
    n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
    Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
    intervals.type = intervals.type, mx.algebras = mx.algebras, 
    model.name = model.name, suppressWarnings = suppressWarnings, 
    silent = silent, run = run)

95% confidence intervals: z statistic approximation
Coefficients:
               Estimate Std.Error    lbound    ubound z value Pr(>|z|)    
AonAlpha       0.747647  0.023880  0.700842  0.794451 31.3082   <2e-16 ***
ConAlpha       0.911705  0.019864  0.872772  0.950638 45.8972   <2e-16 ***
ESonAlpha      0.677435  0.025862  0.626746  0.728124 26.1940   <2e-16 ***
EonBeta        0.152566  0.159051 -0.159168  0.464300  0.9592   0.3374    
IonBeta        3.283787  3.361518 -3.304667  9.872241  0.9769   0.3286    
AlphawithBeta  0.117259  0.125361 -0.128444  0.362961  0.9354   0.3496    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Goodness-of-fit indices:
                                               Value
Sample size                                 838.0000
Chi-square of target model                  145.6167
DF of target model                            4.0000
p value of target model                       0.0000
Number of constraints imposed on "Smatrix"    0.0000
DF manually adjusted                          0.0000
Chi-square of independence model           2480.2312
DF of independence model                     10.0000
RMSEA                                         0.2057
RMSEA lower 95% CI                            0.1778
RMSEA upper 95% CI                            0.2350
SRMR                                          0.1051
TLI                                           0.8567
CFI                                           0.9427
AIC                                         137.6167
BIC                                         118.6926
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)

We may also plot the parameter estimates of these two groups.

library(semPlot)

## Convert the model to semPlotModel object with 2 plots
my.plots <- lapply(X=cluster2, FUN=meta2semPlot, latNames=c("Alpha","Beta"))

## Setup two plots
layout(t(1:2))
semPaths(my.plots[[1]], whatLabels="est", nCharNodes=10, color="green")
title("Younger") 
semPaths(my.plots[[2]], whatLabels="est", nCharNodes=10, color="yellow")
title("Older")

Random-effects TSSEM

Stage 1 analysis

  • Since there is not enough data to estimate the full variance component of the random effects, I estimate the variance component with a diagonal matrix in the RE.type="Diag" argument. The range of \(I^2\) indices, the percentage of total variance that can be explained by the between study effect, are from .84 to .95.
#### Random-effects TSSEM with random effects on the diagonals
random1 <- tssem1(Digman97$data, Digman97$n, method="REM", RE.type="Diag")
summary(random1)

Call:
meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
    "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
    I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
    silent = silent, run = run)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
               Estimate   Std.Error      lbound      ubound z value
Intercept1   0.38971914  0.05429255  0.28330769  0.49613059  7.1781
Intercept2   0.43265886  0.04145137  0.35141568  0.51390205 10.4377
Intercept3   0.04945639  0.06071079 -0.06953458  0.16844736  0.8146
Intercept4   0.09603713  0.04478711  0.00825601  0.18381825  2.1443
Intercept5   0.42724248  0.03911620  0.35057613  0.50390883 10.9224
Intercept6   0.11929323  0.04106204  0.03881311  0.19977336  2.9052
Intercept7   0.19292431  0.04757962  0.09966997  0.28617866  4.0548
Intercept8   0.22690171  0.03240893  0.16338138  0.29042204  7.0012
Intercept9   0.18105574  0.04258856  0.09758370  0.26452778  4.2513
Intercept10  0.43614971  0.03205960  0.37331405  0.49898536 13.6043
Tau2_1_1     0.03648371  0.01513055  0.00682839  0.06613903  2.4113
Tau2_2_2     0.01963098  0.00859789  0.00277942  0.03648254  2.2832
Tau2_3_3     0.04571438  0.01952285  0.00745030  0.08397846  2.3416
Tau2_4_4     0.02236121  0.00995083  0.00285794  0.04186447  2.2472
Tau2_5_5     0.01729072  0.00796404  0.00168149  0.03289995  2.1711
Tau2_6_6     0.01815482  0.00895897  0.00059557  0.03571408  2.0264
Tau2_7_7     0.02604882  0.01130266  0.00389601  0.04820163  2.3047
Tau2_8_8     0.00988761  0.00513713 -0.00018097  0.01995619  1.9247
Tau2_9_9     0.01988244  0.00895053  0.00233973  0.03742515  2.2214
Tau2_10_10   0.01049221  0.00501577  0.00066147  0.02032295  2.0918
             Pr(>|z|)    
Intercept1  7.068e-13 ***
Intercept2  < 2.2e-16 ***
Intercept3    0.41529    
Intercept4    0.03201 *  
Intercept5  < 2.2e-16 ***
Intercept6    0.00367 ** 
Intercept7  5.018e-05 ***
Intercept8  2.538e-12 ***
Intercept9  2.126e-05 ***
Intercept10 < 2.2e-16 ***
Tau2_1_1      0.01590 *  
Tau2_2_2      0.02242 *  
Tau2_3_3      0.01920 *  
Tau2_4_4      0.02463 *  
Tau2_5_5      0.02992 *  
Tau2_6_6      0.04272 *  
Tau2_7_7      0.02119 *  
Tau2_8_8      0.05426 .  
Tau2_9_9      0.02633 *  
Tau2_10_10    0.03645 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q statistic on the homogeneity of effect sizes: 1220.333
Degrees of freedom of the Q statistic: 130
P value of the Q statistic: 0

Heterogeneity indices (based on the estimated Tau2):
                              Estimate
Intercept1: I2 (Q statistic)    0.9326
Intercept2: I2 (Q statistic)    0.8872
Intercept3: I2 (Q statistic)    0.9325
Intercept4: I2 (Q statistic)    0.8703
Intercept5: I2 (Q statistic)    0.8797
Intercept6: I2 (Q statistic)    0.8480
Intercept7: I2 (Q statistic)    0.8887
Intercept8: I2 (Q statistic)    0.7669
Intercept9: I2 (Q statistic)    0.8590
Intercept10: I2 (Q statistic)   0.8193

Number of studies (or clusters): 14
Number of observed statistics: 140
Number of estimated parameters: 20
Degrees of freedom: 120
-2 log likelihood: -112.4196 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
  • The pooled correlation coefficients (fixed effects) and the variance components (the random effects) can be extracted by the use of the coef() command via the select argument.
## Fixed effects
coef(random1, select="fixed")
 Intercept1  Intercept2  Intercept3  Intercept4  Intercept5  Intercept6 
 0.38971914  0.43265886  0.04945639  0.09603713  0.42724248  0.11929323 
 Intercept7  Intercept8  Intercept9 Intercept10 
 0.19292431  0.22690171  0.18105574  0.43614971 
## Random effects
coef(random1, select="random")
  Tau2_1_1   Tau2_2_2   Tau2_3_3   Tau2_4_4   Tau2_5_5   Tau2_6_6 
0.03648371 0.01963098 0.04571438 0.02236121 0.01729072 0.01815482 
  Tau2_7_7   Tau2_8_8   Tau2_9_9 Tau2_10_10 
0.02604882 0.00988761 0.01988244 0.01049221 

Stage 2 analysis

  • The steps are exactly the same as those in the fixed-effects model. The chi-square test statistic of the Stage 2 analysis was \(\chi^2(4, N=4,496) = 8.51, p < .001\), CFI=0.9911, RMSEA=0.0158 and SRMR=0.0463.
random2 <- tssem2(random1, Amatrix=A1, Smatrix=S1, Fmatrix=F1, 
                  diag.constraints=FALSE, intervals="LB")
summary(random2)

Call:
wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
    Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
    diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
    intervals.type = intervals.type, mx.algebras = mx.algebras, 
    model.name = model.name, suppressWarnings = suppressWarnings, 
    silent = silent, run = run)

95% confidence intervals: Likelihood-based statistic
Coefficients:
              Estimate Std.Error  lbound  ubound z value Pr(>|z|)
AonAlpha       0.56944        NA 0.46903 0.67529      NA       NA
ConAlpha       0.59063        NA 0.48968 0.69723      NA       NA
ESonAlpha      0.76045        NA 0.64855 0.89658      NA       NA
EonBeta        0.67996        NA 0.54690 0.77379      NA       NA
IonBeta        0.64184        NA 0.50459 0.79832      NA       NA
AlphawithBeta  0.37769        NA 0.28683 0.47386      NA       NA

Goodness-of-fit indices:
                                               Value
Sample size                                4496.0000
Chi-square of target model                    7.8204
DF of target model                            4.0000
p value of target model                       0.0984
Number of constraints imposed on "Smatrix"    0.0000
DF manually adjusted                          0.0000
Chi-square of independence model            501.6766
DF of independence model                     10.0000
RMSEA                                         0.0146
RMSEA lower 95% CI                            0.0000
RMSEA upper 95% CI                            0.0297
SRMR                                          0.0436
TLI                                           0.9806
CFI                                           0.9922
AIC                                          -0.1796
BIC                                         -25.8234
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
## Plot the parameter estimates
plot(random2, color="green")