The metaSEM Package

Table of Contents

1 Introduction

The metaSEM package provides functions to conducting univariate and multivariate meta-analysis using a structural equation modeling approach via the OpenMx package. It also implemented the two-stage structural equation modeling (TSSEM) approach to conducting fixed- and random-effects meta-analytic structural equation modeling (MASEM) on correlation/covariance matrices.

It is based on the following papers:

  • Cheung, M.W.-L. (2008). A model for integrating fixed-, random-, and mixed-effects meta-analyses into structural equation modeling. Psychological Methods, 13, 182-202.
  • Cheung, M.W.-L. (2009). Constructing approximate confidence intervals for parameters with structural equation models. Structural Equation Modeling, 16, 267–294.
  • Cheung, M.W.-L. (2010). Fixed-effects meta-analyses as multiple-group structural equation models. Structural Equation Modeling, 17, 481-509.
  • Cheung, M.W.-L. (2013). Implementing restricted maximum likelihood estimation in structural equation models. Structural Equation Modeling, 20, 157-167.
  • Cheung, M.W.-L. (2013). Multivariate meta-analysis as structural equation models. Structural Equation Modeling, 20, 429-454.
  • Cheung, M.W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46, 29-40.
  • Cheung, M.W.-L. (in press). Modeling dependent effect sizes with three-level meta-analyses: A structural equation modeling approach. Psychological Methods.
  • Cheung, M.W.-L., & Chan, W. (2004). Testing dependent correlation coefficients via structural equation modeling. Organizational Research Methods, 7, 206–223.
  • Cheung, M.W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10, 40-64.
  • Cheung, M.W.-L., & Chan, W. (2009). A two-stage approach to synthesizing covariance matrices in meta-analytic structural equation modeling. Structural Equation Modeling, 6, 28-53.

2 Examples

Univariate random-effects model

## Load the library
library(metaSEM)

## 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, model.name="Univariate random effects model") )
Running Univariate random effects model 

Call:
meta(y = di, v = vi, data = Becker83, model.name = "Univariate random effects model")

95% confidence intervals: z statistic approximation
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 indicate problems.)

Univariate mixed-effects model

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

Call:
meta(y = di, v = vi, x = log(items), data = Becker83, model.name = "Univariate mixed effects model")

95% confidence intervals: z statistic approximation
Coefficients:
              Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)
Intercept1 -3.2015e-01  1.0981e-01 -5.3539e-01 -1.0492e-01 -2.9154  0.003552
Slope1_1    2.1088e-01  4.5084e-02  1.2251e-01  2.9924e-01  4.6774 2.905e-06
Tau2_1_1    1.0000e-10  2.0095e-02 -3.9386e-02  3.9386e-02  0.0000  1.000000
              
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: 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 indicate problems.)

Multivariate 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
summary( mult1 <- meta(y=cbind(PD, AL), v=cbind(var_PD, cov_PD_AL, var_AL), 
                       data=Berkey98, model.name="Multivariate random effects model") )
 Running Multivariate random effects model 

Call:
meta(y = cbind(PD, AL), v = cbind(var_PD, cov_PD_AL, var_AL), 
    data = Berkey98, model.name = "Multivariate random effects model")

95% confidence intervals: z statistic approximation
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    
---
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 indicate problems.)
## Plot the effect sizes
plot(mult1)

mult1.png

Multivariate 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), model.name="Multivariate mixed effects model") )
 Running Multivariate mixed effects model 

Call:
meta(y = cbind(PD, AL), v = cbind(var_PD, cov_PD_AL, var_AL), 
    x = scale(pub_year, center = 1979), data = Berkey98, model.name = "Multivariate mixed effects model")

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
Intercept1  0.3440000  0.0857659  0.1759019  0.5120981  4.0109 6.048e-05 ***
Intercept2 -0.2918176  0.1312797 -0.5491210 -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.3882923  0.2471147 -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    
---
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 indicate problems.)

Three-level mixed-effects model

## Show the data set
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
## No predictor
summary( meta3(y=y, v=v, cluster=District, data=Cooper03, model.name="Three level random effects model") )
Running Three level random effects model 

Call:
meta3(y = y, v = v, cluster = District, data = Cooper03, model.name = "Three level random effects model")

95% confidence intervals: z statistic approximation
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Intercept  0.1844554  0.0805411  0.0265977  0.3423130  2.2902 0.022010 * 
Tau2_2     0.0328648  0.0111397  0.0110314  0.0546982  2.9502 0.003175 **
Tau2_3     0.0577384  0.0307423 -0.0025154  0.1179921  1.8781 0.060362 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on 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 indicate problems.)
## Year as a predictor
summary( meta3(y=y, v=v, cluster=District, x=scale(Year, scale=FALSE), data=Cooper03, 
               model.name="Three level mixed effects model") )
 Running Three level mixed effects model 

Call:
meta3(y = y, v = v, cluster = District, x = scale(Year, scale = FALSE), 
    data = Cooper03, model.name = "Three level mixed effects model")

95% confidence intervals: z statistic approximation
Coefficients:
            Estimate  Std.Error     lbound     ubound z value Pr(>|z|)   
Intercept  0.1780268  0.0805219  0.0202067  0.3358469  2.2109 0.027042 * 
Slope_1    0.0050737  0.0085266 -0.0116382  0.0217856  0.5950 0.551814   
Tau2_2     0.0329390  0.0111620  0.0110618  0.0548162  2.9510 0.003168 **
Tau2_3     0.0564628  0.0300330 -0.0024007  0.1153264  1.8800 0.060104 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on 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 indicate problems.)

Fixed-effects MASEM

## First stage analysis
## Replicate Becker's (1992) analysis using 4 studies only
fixed1 <- tssem1(Becker92$data[1:4], Becker92$n[1:4], method="FEM")
summary(fixed1)
Running TSSEM1 Analysis of Correlation Matrix

Call:
tssem1FEM(my.df = my.df, n = n, cor.analysis = cor.analysis, 
    model.name = model.name, cluster = cluster, suppressWarnings = suppressWarnings)

Coefficients:
       Estimate Std.Error z value  Pr(>|z|)    
S[1,2] 0.403620  0.047967  8.4145 < 2.2e-16 ***
S[1,3] 0.446522  0.046116  9.6826 < 2.2e-16 ***
S[2,3] 0.224942  0.054269  4.1449 3.399e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Goodness-of-fit indices:
                                    Value
Sample size                      311.0000
Chi-square of target model         9.3651
DF of target model                 9.0000
p value of target model            0.4043
Chi-square of independence model 132.1225
DF of independence model          12.0000
RMSEA                              0.0229
SRMR                               0.0833
TLI                                0.9959
CFI                                0.9970
AIC                               -8.6349
BIC                              -42.2930
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
## Prepare a regression model using create.mxMatrix()
A1 <- create.mxMatrix(c(0,0,0,"0.2*Spa2Apt",
                        0,0,"0.2*Ver2Apt",0,0), 
                      type="Full", ncol=3, nrow=3, name="A1")
A1
FullMatrix 'A1' 

@labels
     [,1] [,2]      [,3]     
[1,] NA   "Spa2Apt" "Ver2Apt"
[2,] NA   NA        NA       
[3,] NA   NA        NA       

@values
     [,1] [,2] [,3]
[1,]    0  0.2  0.2
[2,]    0  0.0  0.0
[3,]    0  0.0  0.0

@free
      [,1]  [,2]  [,3]
[1,] FALSE  TRUE  TRUE
[2,] FALSE FALSE FALSE
[3,] FALSE FALSE FALSE

@lbound: No lower bounds assigned.

@ubound: No upper bounds assigned.
S1 <- create.mxMatrix(c("0.2*ErrorVarApt",0,0,1,"0.2*CorBetSpaVer",1),
                        type="Symm", name="S1")
S1
SymmMatrix 'S1' 

@labels
     [,1]          [,2]           [,3]          
[1,] "ErrorVarApt" NA             NA            
[2,] NA            NA             "CorBetSpaVer"
[3,] NA            "CorBetSpaVer" NA            

@values
     [,1] [,2] [,3]
[1,]  0.2  0.0  0.0
[2,]  0.0  1.0  0.2
[3,]  0.0  0.2  1.0

@free
      [,1]  [,2]  [,3]
[1,]  TRUE FALSE FALSE
[2,] FALSE FALSE  TRUE
[3,] FALSE  TRUE FALSE

@lbound: No lower bounds assigned.

@ubound: No upper bounds assigned.
## Fixed-effects model: Second stage analysis
## Two equivalent versions to calculate the R2 and its 95% LBCI
fixed2 <- tssem2(fixed1, Amatrix=A1, Smatrix=S1, intervals.type="LB",
       mx.algebras=list(R1=mxAlgebra(Spa2Apt^2+Ver2Apt^2+2*CorBetSpaVer*Spa2Apt*Ver2Apt, name="R1"),
                        R2=mxAlgebra(One-Ematrix[1,1], name="R2"),
                        One=mxMatrix("Iden", ncol=1, nrow=1, name="One")))
summary(fixed2)
 Running TSSEM2 (Fixed Effects Model) Analysis of Correlation Structure

Call:
wls(Cov = tssem1.obj$pooledS, asyCov = tssem1.obj$acovS, n = tssem1.obj$total.n, 
    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)

95% confidence intervals: Likelihood-based statistic
Coefficients:
             Estimate Std.Error  lbound  ubound z value Pr(>|z|)
Amatrix[1,2]  0.31934        NA 0.22585 0.41262      NA       NA
Amatrix[1,3]  0.37469        NA 0.28265 0.46652      NA       NA
Smatrix[2,3]  0.22494        NA 0.11849 0.33141      NA       NA

mxAlgebras objects (and their 95% likelihood-based CIs):
           lbound Estimate    ubound
R1[1,1] 0.2154344 0.296198 0.3876107
R2[1,1] 0.2154344 0.296198 0.3876107

Goodness-of-fit indices:
                                            Value
Sample size                                311.00
Chi-square of target model                   0.00
DF of target model                           0.00
p value of target model                      0.00
Number of constraints imposed on "Smatrix"   0.00
DF manually adjusted                         0.00
Chi-square of independence model           151.14
DF of independence model                     3.00
RMSEA                                        0.00
SRMR                                         0.00
TLI                                          -Inf
CFI                                          1.00
AIC                                          0.00
BIC                                          0.00
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)

Random-effects MASEM

## First stage analysis
## No random effects for off-diagonal elements
random1 <- tssem1(Becker92$data, Becker92$n, method="REM", RE.type="Diag")
summary(random1)
Running TSSEM1 (Random Effects Model) Analysis of Correlation Matrix

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

95% confidence intervals: z statistic approximation
Coefficients:
              Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)
Intercept1  3.6738e-01  3.9364e-02  2.9023e-01  4.4453e-01  9.3328 < 2.2e-16
Intercept2  3.2299e-01  8.5164e-02  1.5608e-01  4.8991e-01  3.7926 0.0001491
Intercept3  1.7591e-01  4.1379e-02  9.4807e-02  2.5701e-01  4.2511 2.127e-05
Tau2_1_1    1.0000e-10  4.0975e-03 -8.0309e-03  8.0309e-03  0.0000 1.0000000
Tau2_2_2    3.4481e-02  2.9090e-02 -2.2535e-02  9.1496e-02  1.1853 0.2359001
Tau2_3_3    1.0000e-10  1.0738e-02 -2.1046e-02  2.1046e-02  0.0000 1.0000000
              
Intercept1 ***
Intercept2 ***
Intercept3 ***
Tau2_1_1      
Tau2_2_2      
Tau2_3_3      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on the homogeneity of effect sizes: 48.11435
Degrees of freedom of the Q statistic: 15
P value of the Q statistic: 2.437274e-05
Heterogeneity indices (based on the estimated Tau2):
                             Estimate
Intercept1: I2 (Q statistic)    0.000
Intercept2: I2 (Q statistic)    0.811
Intercept3: I2 (Q statistic)    0.000

Number of studies (or clusters): 6
Number of observed statistics: 18
Number of estimated parameters: 6
Degrees of freedom: 12
-2 log likelihood: -21.40985 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
## Random-effects model: Second stage analysis
random2 <- tssem2(random1, Amatrix=A1, Smatrix=S1, intervals.type="z")
summary(random2)
Running TSSEM2 (Random Effects Model) Analysis of Correlation Structure

Call:
wls(Cov = pooledS, asyCov = asyCov, n = tssem1.obj$total.n, 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)

95% confidence intervals: z statistic approximation
Coefficients:
             Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
Amatrix[1,2] 0.320476  0.043864 0.234504 0.406448  7.3061 2.749e-13 ***
Amatrix[1,3] 0.266619  0.087133 0.095841 0.437397  3.0599  0.002214 ** 
Smatrix[2,3] 0.175909  0.041379 0.094807 0.257011  4.2511 2.127e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Goodness-of-fit indices:
                                            Value
Sample size                                538.00
Chi-square of target model                   0.00
DF of target model                           0.00
p value of target model                      0.00
Number of constraints imposed on "Smatrix"   0.00
DF manually adjusted                         0.00
Chi-square of independence model           107.23
DF of independence model                     3.00
RMSEA                                        0.00
SRMR                                         0.00
TLI                                          -Inf
CFI                                          1.00
AIC                                          0.00
BIC                                          0.00
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
  1. More examples

3 Installation

The current version is 0.8-4:

First of all, you need R to run it. Since metaSEM uses OpenMx as the workhorse, OpenMx has to be installed. To install OpenMx, run the following command inside an R session:

install.packages('OpenMx', repos='http://openmx.psyc.virginia.edu/packages/')

See http://openmx.psyc.virginia.edu/installing-openmx for the details on how to install OpenMx. Moreover, the metaSEM also depends on the ellipse and the MASS packages that can be installed by the following command inside an R session:

install.packages(c('ellipse','MASS'))

Windows platform

Download the Windows binary of metaSEM. If the file is saved at d:\. Run the following command inside an R session:

install.packages(pkgs="d:/metaSEM_0.8-4.zip", repos=NULL)

Please note that d:\ in Windows is represented by either d:/ or d:\\ in R.

Linux and Mac OS X platform

Download the source package of metaSEM. Run the following command (as Root) inside an R session:

install.packages(pkgs="metaSEM_0.8-4.tar.gz", repos=NULL, type="source")

Alternatively, the metaSEM package may also be installed inside the terminal:

R CMD INSTALL metaSEM_0.8-4.tar.gz

4 Help

  1. OpenMx discussion forum

    There is a discussion forum for the metaSEM package in OpenMx. You may post technical questions related to metaSEM there. Please include information on the session. It will be helpful if you can include a reproducible example. You may save a copy of your data, say my.df, by using

    sessionInfo()
    hdump(c("my.df"), file="myData.R")
    

    Then you may include the content of "myData.R" in the post.

  2. Meta-analysis resources

    There is another discussion forum related to the conceptual issues of meta-analytic structural equation modeling at Meta-analysis resources.

Author: Mike W.-L. Cheung

Created: 2014-03-20 Thu 16:35

Emacs 24.3.50.1 (Org mode 8.2.4)

Validate