# The metaSEM package: Examples

#### Mike W.-L. Cheung

#### 29 November, 2020

# Introduction

- The `metaSEM` package provides functions for 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 a fixed- and random-effects meta-analytic structural equation modeling (MASEM) on correlation/covariance matrices.

- The `metaSEM` package 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(3), 182–202. https://doi.org/10.1037/a0013163
  - Cheung, M. W.-L. (2009). Constructing approximate confidence intervals for parameters with structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 16(2), 267–294. https://doi.org/10.1080/10705510902751291
  - Cheung, M. W.-L. (2010). Fixed-effects meta-analyses as multiple-group structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 17(3), 481–509. https://doi.org/10.1080/10705511.2010.489367
  - Cheung, M. W.-L. (2013). Implementing restricted maximum likelihood estimation in structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 20(1), 157–167. https://doi.org/10.1080/10705511.2013.742404
  - Cheung, M. W.-L. (2013). Multivariate meta-analysis as structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 20(3), 429–454. https://doi.org/10.1080/10705511.2013.797827
  - Cheung, M. W.-L. (2014). Modeling dependent effect sizes with three-level meta-analyses: A structural equation modeling approach. Psychological Methods, 19(2), 211–229. https://doi.org/10.1037/a0032968
  - Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29–40. https://doi.org/10.3758/s13428-013-0361-y
  - Cheung, M. W.-L. (2015). metaSEM: an R package for meta-analysis using structural equation modeling. Frontiers in Psychology, 5(1521). https://doi.org/10.3389/fpsyg.2014.01521
  - Cheung, M. W.-L. (2015). Meta-Analysis: A Structural Equation Modeling Approach. Wiley.
  - Cheung, M.W.-L. (2018). Computing multivariate effect sizes and their sampling covariance matrices with structural equation modeling: Theory, examples, and computer simulations. Frontiers in Psychology, 9(1387). https://doi.org/10.3389/fpsyg.2018.01387
  - Cheung, M. W.-L. (2018). Issues in solving the problem of effect size heterogeneity in meta-analytic structural equation modeling: A commentary and simulation study on Yu, Downes, Carter, and O’Boyle (2016). Journal of Applied Psychology, 103(7), 787–803. https://doi.org/10.1037/apl0000284
  - Cheung, M. W.-L. (2019). Some reflections on combining meta-analysis and structural equation modeling. Research Synthesis Methods, 10(1), 15-22. https://doi.org/10.1002/jrsm.1321
  - Cheung, M. W.-L., & Chan, W. (2004). Testing dependent correlation coefficients via structural equation modeling. Organizational Research Methods, 7(2), 206–223. https://doi.org/10.1177/1094428104264024
  - Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40–64. https://doi.org/10.1037/1082-989X.10.1.40
  - Cheung, M. W.-L., & Chan, W. (2005). Classifying correlation matrices into relatively homogeneous subgroups: A cluster analytic approach. Educational and Psychological Measurement, 65(6), 954–979. https://doi.org/10.1177/0013164404273946
  - Cheung, M. W.-L., & Chan, W. (2009). A two-stage approach to synthesizing covariance matrices in meta-analytic structural equation modeling. Structural Equation Modeling: A Multidisciplinary Journal, 16(1), 28–53. https://doi.org/10.1080/10705510802561295
  - Cheung, M. W.-L., & Cheung, S. F. (2016). Random-effects models for meta-analytic structural equation modeling: Review, issues, and illustrations. Research Synthesis Methods, 7(2), 140–155. https://doi.org/10.1002/jrsm.1166
  - Cheung, M. W.-L., & Hafdahl, A. R. (2016). Special issue on meta-analytic structural equation modeling: Introduction from the guest editors. Research Synthesis Methods, 7(2), 112–120. https://doi.org/10.1002/jrsm.1212
  - Cheung, M. W.-L., & Hong, R. Y. (2017). Applications of meta-analytic structural equation modelling in health psychology: examples, issues, and recommendations. Health Psychology Review, 11(3), 265–279. https://doi.org/10.1080/17437199.2017.1343678
  - Jak, S., & Cheung, M. W.-L. (2018). Testing moderator hypotheses in meta-analytic structural equation modeling using subgroup analysis. Behavior Research Methods, 50(4), 1359–1373. https://doi.org/10.3758/s13428-018-1046-3
  - Jak, S., & Cheung, M. W.-L. (2019). Meta-analytic structural equation modeling with moderating effects on SEM parameters.Psychological Methods. https://doi.org/10.1037/met0000245

# 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  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 may indicate problems.)

## Handling missing covariates with FIML

    ## Sample data from Tenenbaum and Leaper (2002, Table 1).
    Tenenbaum02 <- Tenenbaum02[, c("r", "v", "Offspring_age", "Year_pub")]

    ## Set seed for reproducibility
    set.seed(1234567)

    ## Let's drop 40% in Offspring_age
    missing_per <- 0.4

    ## MCAR
    index <- round(nrow(Tenenbaum02)*missing_per)
    index <- rep(c(TRUE, FALSE), times=c(index, nrow(Tenenbaum02)-index))
    index <- sample(index)
    my.MCAR <- Tenenbaum02
    my.MCAR[index, "Offspring_age"] <- NA
    my.MCAR$Offspring_age <- scale(my.MCAR$Offspring_age, scale=FALSE)
    my.MCAR$Year_pub <- scale(my.MCAR$Year_pub, scale=FALSE)

    my.MCAR

           r            v Offspring_age     Year_pub
    1   0.12 0.0029084053            NA  -4.91666667
    2  -0.08 0.0099721309   -48.3448276   3.08333333
    3  -0.05 0.0103646484   -72.3448276  -7.91666667
    4  -0.08 0.0094022949   -69.3448276  -8.91666667
    5   0.15 0.0013089127   209.6551724  -1.91666667
    6   0.12 0.0404753067   -60.3448276   0.08333333
    7   0.17 0.0052101393            NA   0.08333333
    8   0.34 0.0181898456   107.6551724  -3.91666667
    9  -0.01 0.0020238867            NA   5.08333333
    10  0.33 0.0048124801            NA  10.08333333
    11  0.40 0.0147000000   -93.3448276   1.08333333
    12  0.30 0.0138016667   -90.3448276   4.08333333
    13  0.07 0.0162331805   -57.3448276  13.08333333
    14 -0.02 0.0099920016            NA  -2.91666667
    15  0.19 0.0092910321    -5.3448276  -2.91666667
    16  0.15 0.0006964331    23.6551724  11.08333333
    17  0.19 0.0032148900   -18.3448276   7.08333333
    18  0.02 0.0057097152   -30.3448276  -7.91666667
    19  0.47 0.0085492508            NA  -6.91666667
    20  0.19 0.0046455161     5.6551724  11.08333333
    21  0.33 0.0124071752            NA  -7.91666667
    22  0.06 0.0157589359   -36.3448276  -7.91666667
    23  0.05 0.0033166875            NA   2.08333333
    24  0.24 0.0024877248            NA   5.08333333
    25  0.14 0.0123228738            NA -10.91666667
    26 -0.02 0.0285485760            NA   5.08333333
    27  0.28 0.0066877682            NA  11.08333333
    28  0.14 0.0034825513    35.6551724   9.08333333
    29  0.27 0.0245575546   -42.3448276  -3.91666667
    30  0.38 0.0140779108   -12.3448276 -10.91666667
    31  0.52 0.0212926464            NA  -5.91666667
    32  0.66 0.0127418944    23.6551724  -5.91666667
    33  0.36 0.0303038464    23.6551724  -5.91666667
    34  0.21 0.0042698356    11.6551724  10.08333333
    35  0.19 0.0042231964    11.6551724  10.08333333
    36  0.14 0.0184843108            NA -12.91666667
    37  0.36 0.0102377859    95.6551724  -7.91666667
    38  0.09 0.0006314927    23.6551724  11.08333333
    39  0.16 0.0053946327            NA  13.08333333
    40 -0.07 0.0319427100            NA  -0.91666667
    41  0.36 0.0140295585    89.6551724  -8.91666667
    42  0.36 0.0008270700    95.6551724  -4.91666667
    43  0.22 0.0238300674    -0.3448276   8.08333333
    44  0.04 0.0066899501   -72.3448276   7.08333333
    45  0.06 0.0146001906            NA  -8.91666667
    46  0.41 0.0150447307            NA  -3.91666667
    47  0.04 0.0140394727            NA  -3.91666667
    48  0.23 0.0067954425   -48.3448276   2.08333333

    ## y: effect size
    ## v: sampling variance
    ## x: covariate
    ## av: auxiliary variable
    fit <- metaFIML(y=r, v=v, x=Offspring_age, av=Year_pub, data=my.MCAR)
    summary(fit)

    Call:
    metaFIML(y = r, v = v, x = Offspring_age, av = Year_pub, data = my.MCAR)

    95% confidence intervals: z statistic approximation (robust=FALSE)
    Coefficients:
                  Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)
    Tau2_1_1    1.2136e-02  4.5317e-03  3.2543e-03  2.1018e-02  2.6781 0.0074043
    CovX1_X1    4.6116e+03  1.2071e+03  2.2457e+03  6.9775e+03  3.8203 0.0001333
    CovX2_X1   -6.0799e+01  9.3447e+01 -2.4395e+02  1.2235e+02 -0.6506 0.5152919
    CovX2_X2    5.7118e+01  1.1659e+01  3.4267e+01  7.9969e+01  4.8990 9.632e-07
    CovX2_Y1   -1.3300e-01  1.6016e-01 -4.4691e-01  1.8090e-01 -0.8305 0.4062771
    Slope1_1    6.8319e-04  3.5911e-04 -2.0647e-05  1.3870e-03  1.9025 0.0571100
    Intercept1  1.8420e-01  2.1882e-02  1.4131e-01  2.2709e-01  8.4179 < 2.2e-16
    MeanX1      7.7303e-02  1.2286e+01 -2.4004e+01  2.4158e+01  0.0063 0.9949800
    MeanX2      1.3438e-07  1.0907e+00 -2.1378e+00  2.1378e+00  0.0000 0.9999999

    Tau2_1_1   **
    CovX1_X1   ***
    CovX2_X1
    CovX2_X2   ***
    CovX2_Y1
    Slope1_1   .
    Intercept1 ***
    MeanX1
    MeanX2
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Q statistic on the homogeneity of effect sizes: 179.0075
    Degrees of freedom of the Q statistic: 47
    P value of the Q statistic: 0

    Explained variances (R2):
                               y1
    Tau2 (no predictor)    0.0142
    Tau2 (with predictors) 0.0121
    R2                     0.1474

    Number of studies (or clusters): 48
    Number of observed statistics: 125
    Number of estimated parameters: 9
    Degrees of freedom: 116
    -2 log likelihood: 613.2553
    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
    ---
    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.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 may indicate problems.)

# Three-level meta-analysis

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

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

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

## 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 ***
    ---
    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    diffLL
    1 3 level model                      <NA>  3 16.78987 53 -89.21013        NA
    2 3 level model Equal tau2 at both levels  2 17.47697 54 -90.52303 0.6870959
      diffdf         p
    1     NA        NA
    2      1 0.4071539

- 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.158399      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.34963  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.0602799      NA       NA
    Tau2_3     0.0564628        NA  0.0192365  0.1716333      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 diffdf
    1 Variance component with REML       <NA>  2 -81.14044 -2  NA       NA     NA
    2 Variance component with REML Equal Tau2  1 -80.13710 -1  NA 1.003336      1
              p
    1        NA
    2 0.3165046

- We may also estimate the residual heterogeneity after controlling for the covariate. The estimated residual heterogeneity for \(\tau^2_{(2)}\) and \(\tau^2_{(3)}\) were 0.0327 and 0.0723, respectively.

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

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

## 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 diffdf
    1 3 level model                 <NA>  3 25.80256 63 -100.1974       NA     NA
    2 3 level model tau2_2 equals tau2_3  2 27.16166 64 -100.8383 1.359103      1
             p
    1       NA
    2 0.243693

## 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
    ---
    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.00095286  0.02365224 -0.04731040  0.04540467 -0.0403  0.96786
    Slope_2   -0.01176840  0.00659995 -0.02470407  0.00116727 -1.7831  0.07457 .
    Tau2_2     0.00287389  0.00206817 -0.00117965  0.00692744  1.3896  0.16466
    Tau2_3     0.01479446  0.00926095 -0.00335666  0.03294558  1.5975  0.11015
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Q statistic on 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.2430133  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.13998846  0.11049279 -0.2308  0.81747
    Slope_1   -0.00913064  0.09949198 -0.20413135  0.18587006 -0.0918  0.92688
    Slope_2   -0.12617957  0.07866274 -0.28035571  0.02799657 -1.6041  0.10870
    Slope_3   -0.23695698  0.12123092 -0.47456521  0.00065126 -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
    ---
    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         p
    1 Model 3          <NA>  6 20.07571 60  -99.92429       NA     NA        NA
    2 Model 3 3 level model  3 25.80256 63 -100.19744 5.726842      3 0.1256832

## 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.0992866 -1.1261 0.260127
    Slope_3   -0.2210801  0.0739174 -0.3659556 -0.0762046 -2.9909 0.002782 **
    Slope_4    0.0537251  0.0994803 -0.1412527  0.2487030  0.5401 0.589157
    Tau2_2     0.0033376  0.0023492 -0.0012667  0.0079420  1.4208 0.155383
    Tau2_3     0.0047979  0.0044818 -0.0039862  0.0135820  1.0705 0.284379
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Q statistic on 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          p
    1 Model 4          <NA>  7 14.18259 59 -103.8174       NA     NA         NA
    2 Model 4 3 level model  3 25.80256 63 -100.1974 11.61996      4 0.02041284

## 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  Pr(>|z|)
    Intercept  6.7039e-02  1.8555e-02  3.0673e-02  1.0341e-01  3.6131 0.0003026 ***
    Slope_1   -1.9004e-01  4.0234e-02 -2.6889e-01 -1.1118e-01 -4.7233 2.321e-06 ***
    Slope_2    1.9491e-02  6.5942e-02 -1.0975e-01  1.4874e-01  0.2956 0.7675469
    Slope_3   -1.2779e-01  3.5915e-02 -1.9818e-01 -5.7401e-02 -3.5582 0.0003734 ***
    Slope_4   -2.3955e-01  9.4054e-02 -4.2389e-01 -5.5207e-02 -2.5469 0.0108674 *
    Tau2_2     2.3061e-03  1.4269e-03 -4.9055e-04  5.1028e-03  1.6162 0.1060566
    Tau2_3     1.6024e-08          NA          NA          NA      NA        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.0023061  0.0000
    R2                     0.3925670  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.667327
    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.

- 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.667327 59 -113.3327       NA     NA
    2 Model 5 Meta analysis with ML  4 17.625692 62 -106.3743 12.95836      3
                p
    1          NA
    2 0.004727515

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

    Q statistic on 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.2308748  0.3962360  0.5168  0.60528
    Slope_3   -0.0462265  0.1715773 -0.3825119  0.2900589 -0.2694  0.78761
    Slope_4   -0.0486321  0.1306918 -0.3047834  0.2075191 -0.3721  0.70981
    Slope_5   -0.2169132  0.1915703 -0.5923841  0.1585576 -1.1323  0.25751
    Slope_6   -0.3036578  0.1670720 -0.6311130  0.0237974 -1.8175  0.06914 .
    Slope_7   -0.0605272  0.1809419 -0.4151668  0.2941124 -0.3345  0.73799
    Tau2_2     0.0032661  0.0022784 -0.0011994  0.0077317  1.4335  0.15171
    Tau2_3     0.0040618  0.0038459 -0.0034759  0.0115996  1.0562  0.29090
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Q statistic on 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  Pr(>|z|)
    Intercept  6.8563e-02  1.8630e-02  3.2049e-02  1.0508e-01  3.6802  0.000233 ***
    Slope_1   -1.6885e-01  4.1545e-02 -2.5028e-01 -8.7424e-02 -4.0643 4.818e-05 ***
    Slope_2    2.5329e-01  1.5814e-01 -5.6670e-02  5.6325e-01  1.6016  0.109240
    Slope_3    1.2689e-01  1.4774e-01 -1.6268e-01  4.1646e-01  0.8589  0.390411
    Slope_4   -8.3548e-03  1.5796e-01 -3.1795e-01  3.0124e-01 -0.0529  0.957818
    Slope_5   -1.1530e-01  1.1147e-01 -3.3377e-01  1.0317e-01 -1.0344  0.300948
    Slope_6   -2.6412e-01  1.6402e-01 -5.8559e-01  5.7344e-02 -1.6103  0.107324
    Slope_7   -2.9029e-01  1.4859e-01 -5.8152e-01  9.5211e-04 -1.9536  0.050754 .
    Slope_8   -1.5975e-01  1.6285e-01 -4.7893e-01  1.5943e-01 -0.9810  0.326609
    Tau2_2     2.1010e-03  1.2925e-03 -4.3226e-04  4.6342e-03  1.6255  0.104051
    Tau2_3     1.0000e-10          NA          NA          NA      NA        NA
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

## 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
    ---
    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.9461940
    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
    ---
    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.5592669  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.17573647  0.06328327 -0.29976940 -0.05170355 -2.7770 0.005487 **
    Tau2_2     0.00259266  0.00171596 -0.00077056  0.00595588  1.5109 0.130811
    Tau2_3     0.00278384  0.00267150 -0.00245221  0.00801989  1.0421 0.297388
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Q statistic on 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.1238925  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.0572059 -0.1385272  0.0857157 -0.4616 0.644374
    SlopeX2_1 -0.2003999  0.0691096 -0.3358521 -0.0649476 -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.1760 < 2.2e-16 ***
    S[1,3] 0.390375  0.012880 30.3076 < 2.2e-16 ***
    S[1,4] 0.103572  0.015048  6.8830 5.862e-12 ***
    S[1,5] 0.092286  0.015047  6.1330 8.621e-10 ***
    S[2,3] 0.416070  0.012519 33.2344 < 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. An alternative and easier way to specify the model is to use `lavann`’s syntax.

    model1 <- "## Factor loadings
               Alpha =~ A+C+ES
               Beta =~ E+I
               ## Factor correlation
               Alpha ~~ Beta"

    plot(model1)

    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

- 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, RAM=RAM1, 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
    EonBeta        0.78141        NA 0.71872 0.85498      NA       NA
    ESonAlpha      0.71920        NA 0.68876 0.75032      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.7228
    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, RAM=RAM1, 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 ***
    EonBeta       0.967136  0.058751 0.851986 1.082287  16.462 < 2.2e-16 ***
    ESonAlpha     0.732173  0.018929 0.695074 0.769273  38.680 < 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.3313
    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.3081   <2e-16 ***
    ConAlpha       0.911705  0.019864  0.872772  0.950638 45.8970   <2e-16 ***
    EonBeta        0.152564  0.159116 -0.159297  0.464426  0.9588   0.3376
    ESonAlpha      0.677435  0.025863  0.626743  0.728126 26.1928   <2e-16 ***
    IonBeta        3.283814  3.362956 -3.307458  9.875086  0.9765   0.3288
    AlphawithBeta  0.117258  0.125412 -0.128546  0.363061  0.9350   0.3498
    ---
    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.2324
    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  Pr(>|z|)
    Intercept1   0.38971913  0.05429255  0.28330768  0.49613059  7.1781 7.068e-13
    Intercept2   0.43265886  0.04145137  0.35141568  0.51390205 10.4377 < 2.2e-16
    Intercept3   0.04945639  0.06071079 -0.06953457  0.16844736  0.8146   0.41529
    Intercept4   0.09603714  0.04478711  0.00825602  0.18381826  2.1443   0.03201
    Intercept5   0.42724247  0.03911620  0.35057612  0.50390882 10.9224 < 2.2e-16
    Intercept6   0.11929323  0.04106204  0.03881311  0.19977335  2.9052   0.00367
    Intercept7   0.19292431  0.04757962  0.09966997  0.28617866  4.0548 5.018e-05
    Intercept8   0.22690171  0.03240893  0.16338138  0.29042204  7.0012 2.538e-12
    Intercept9   0.18105573  0.04258856  0.09758369  0.26452777  4.2513 2.126e-05
    Intercept10  0.43614971  0.03205960  0.37331405  0.49898536 13.6043 < 2.2e-16
    Tau2_1_1     0.03648371  0.01513054  0.00682839  0.06613903  2.4113   0.01590
    Tau2_2_2     0.01963098  0.00859789  0.00277942  0.03648254  2.2832   0.02242
    Tau2_3_3     0.04571438  0.01952285  0.00745030  0.08397845  2.3416   0.01920
    Tau2_4_4     0.02236121  0.00995083  0.00285794  0.04186447  2.2472   0.02463
    Tau2_5_5     0.01729072  0.00796404  0.00168149  0.03289995  2.1711   0.02992
    Tau2_6_6     0.01815482  0.00895897  0.00059557  0.03571408  2.0264   0.04272
    Tau2_7_7     0.02604882  0.01130266  0.00389601  0.04820162  2.3047   0.02119
    Tau2_8_8     0.00988761  0.00513713 -0.00018097  0.01995619  1.9247   0.05426
    Tau2_9_9     0.01988244  0.00895053  0.00233973  0.03742515  2.2214   0.02633
    Tau2_10_10   0.01049221  0.00501577  0.00066147  0.02032295  2.0918   0.03645

    Intercept1  ***
    Intercept2  ***
    Intercept3
    Intercept4  *
    Intercept5  ***
    Intercept6  **
    Intercept7  ***
    Intercept8  ***
    Intercept9  ***
    Intercept10 ***
    Tau2_1_1    *
    Tau2_2_2    *
    Tau2_3_3    *
    Tau2_4_4    *
    Tau2_5_5    *
    Tau2_6_6    *
    Tau2_7_7    *
    Tau2_8_8    .
    Tau2_9_9    *
    Tau2_10_10  *
    ---
    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.38971913  0.43265886  0.04945639  0.09603714  0.42724247  0.11929323
     Intercept7  Intercept8  Intercept9 Intercept10
     0.19292431  0.22690171  0.18105573  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   Tau2_7_7
    0.03648371 0.01963098 0.04571438 0.02236121 0.01729072 0.01815482 0.02604882
      Tau2_8_8   Tau2_9_9 Tau2_10_10
    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, RAM=RAM1, 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
    EonBeta        0.67996        NA 0.54690 0.81751      NA       NA
    ESonAlpha      0.76045        NA 0.64855 0.89661      NA       NA
    IonBeta        0.64184        NA 0.50459 0.79834      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")

## Becker and Schram (1994)

### Inspect the data

- 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 (not required for the analysis)
    head(Becker94$data)

    $`Becker (1978) Females`
             Math Spatial Verbal
    Math     1.00    0.47  -0.21
    Spatial  0.47    1.00  -0.15
    Verbal  -0.21   -0.15   1.00

    $`Becker (1978) Males`
            Math Spatial Verbal
    Math    1.00    0.28   0.19
    Spatial 0.28    1.00   0.18
    Verbal  0.19    0.18   1.00

    $`Berry (1957) Females`
            Math Spatial Verbal
    Math    1.00    0.48   0.41
    Spatial 0.48    1.00   0.26
    Verbal  0.41    0.26   1.00

    $`Berry (1957) Males`
            Math Spatial Verbal
    Math    1.00    0.37   0.40
    Spatial 0.37    1.00   0.27
    Verbal  0.40    0.27   1.00

    $`Rosenberg (1981) Females`
            Math Spatial Verbal
    Math    1.00    0.42   0.48
    Spatial 0.42    1.00   0.23
    Verbal  0.48    0.23   1.00

    $`Rosenberg (1981) Males`
            Math Spatial Verbal
    Math    1.00    0.41   0.74
    Spatial 0.41    1.00   0.44
    Verbal  0.74    0.44   1.00

    head(Becker94$n)

    [1]  74 153  48  55  51  18

## Fixed-effects TSSEM

### Stage 1 analysis

- The test statistic of the Stage 1 analysis based on the TSSEM approach was \(\chi^2(27, N=538) = 62.50, p < .001\), CFI=0.7943, RMSEA=0.1565 and SRMR=0.2011. Based on the test statistic and the goodness-of-fit indices, the hypothesis of homogeneity of correlation matrices was rejected.

    #### Fixed-effects model
    ## Stage 1 analysis
    fixed1 <- tssem1(Becker94$data, Becker94$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.379961  0.037123 10.2351 < 2.2e-16 ***
    S[1,3] 0.334562  0.039947  8.3751 < 2.2e-16 ***
    S[2,3] 0.176461  0.042334  4.1683 3.069e-05 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Goodness-of-fit indices:
                                         Value
    Sample size                       538.0000
    Chi-square of target model         63.6553
    DF of target model                 27.0000
    p value of target model             0.0001
    Chi-square of independence model  207.7894
    DF of independence model           30.0000
    RMSEA                               0.1590
    RMSEA lower 95% CI                  0.1096
    RMSEA upper 95% CI                  0.2117
    SRMR                                0.1586
    TLI                                 0.7709
    CFI                                 0.7938
    AIC                                 9.6553
    BIC                              -106.1169
    OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
    Other values may indicate problems.)

    coef(fixed1)

                 Math   Spatial   Verbal
    Math    1.0000000 0.3799605 0.334562
    Spatial 0.3799605 1.0000000 0.176461
    Verbal  0.3345620 0.1764610 1.000000

### Stage 2 analysis

We may specify the model via the RAM formulation. If all variables are observed, there is no need to specify the F matrix. Since the df of the regression model is 0, the proposed model always fits the data perfectly.

    #### Prepare models for stage 2 analysis
    model2 <- "## Math is modeled by Spatial and Verbal
               Math ~ Spatial2Math*Spatial + Verbal2Math*Verbal
               ## Variances of predictors are fixed at 1
               Spatial ~~ 1*Spatial
               Verbal ~~ 1*Verbal
               ## Correlation between the predictors
               Spatial ~~ SpatialCorVerbal*Verbal
               ## Error variance
               Math ~~ ErrorVarMath*Math"
    plot(model2)

    RAM2 <- lavaan2RAM(model2)
    RAM2

    $A
            Math Spatial          Verbal
    Math    "0"  "0*Spatial2Math" "0*Verbal2Math"
    Spatial "0"  "0"              "0"
    Verbal  "0"  "0"              "0"

    $S
            Math             Spatial              Verbal
    Math    "0*ErrorVarMath" "0"                  "0"
    Spatial "0"              "1"                  "0*SpatialCorVerbal"
    Verbal  "0"              "0*SpatialCorVerbal" "1"

    $F
            Math Spatial Verbal
    Math       1       0      0
    Spatial    0       1      0
    Verbal     0       0      1

    $M
      Math Spatial Verbal
    1    0       0      0

    ## Stage 2 analysis
    fixed2 <- tssem2(fixed1, RAM=RAM2, intervals="LB")
    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|)
    Spatial2Math     0.331238        NA 0.258216 0.404145      NA       NA
    Verbal2Math      0.276111        NA 0.199501 0.352776      NA       NA
    SpatialCorVerbal 0.176461        NA 0.093488 0.259434      NA       NA

    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           160.47
    DF of independence model                     3.00
    RMSEA                                        0.00
    RMSEA lower 95% CI                           0.00
    RMSEA upper 95% CI                           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.)

    plot(fixed2)

## Fixed-effects TSSEM with cluster

### Stage 1 analysis

- The goodness-of-fit indices of the Stage 1 analysis for the female and male participants were \(\chi^2(12, N=235) = 42.41, p < .001\), CFI=0.7116, RMSEA=0.2327 and SRMR=0.2339, and \(\chi^2(12, N=303) = 16.13, p = .1852\), CFI=0.9385, RMSEA=0.0755 and SRMR=0.1508, respectively.

    #### Fixed-effects model with cluster
    ## Stage 1 analysis
    cluster1 <- tssem1(Becker94$data, Becker94$n, method="FEM", cluster=Becker94$gender)
    summary(cluster1)

    $Females

    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.455896  0.051993  8.7685 < 2.2e-16 ***
    S[1,3] 0.341583  0.061943  5.5144 3.499e-08 ***
    S[2,3] 0.171931  0.064731  2.6561  0.007905 **
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Goodness-of-fit indices:
                                        Value
    Sample size                      235.0000
    Chi-square of target model        43.1898
    DF of target model                12.0000
    p value of target model            0.0000
    Chi-square of independence model 123.4399
    DF of independence model          15.0000
    RMSEA                              0.2357
    RMSEA lower 95% CI                 0.1637
    RMSEA upper 95% CI                 0.3161
    SRMR                               0.2141
    TLI                                0.6405
    CFI                                0.7124
    AIC                               19.1898
    BIC                              -22.3252
    OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
    Other values may indicate problems.)

    $Males

    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.318051  0.051698  6.1521 7.646e-10 ***
    S[1,3] 0.328286  0.052226  6.2858 3.261e-10 ***
    S[2,3] 0.179549  0.055944  3.2094   0.00133 **
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Goodness-of-fit indices:
                                        Value
    Sample size                      303.0000
    Chi-square of target model        16.4819
    DF of target model                12.0000
    p value of target model            0.1701
    Chi-square of independence model  84.3496
    DF of independence model          15.0000
    RMSEA                              0.0786
    RMSEA lower 95% CI                 0.0000
    RMSEA upper 95% CI                 0.1643
    SRMR                               0.1025
    TLI                                0.9192
    CFI                                0.9354
    AIC                               -7.5181
    BIC                              -52.0829
    OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
    Other values may indicate problems.)

    coef(cluster1)

    $Females
                 Math   Spatial    Verbal
    Math    1.0000000 0.4558958 0.3415826
    Spatial 0.4558958 1.0000000 0.1719309
    Verbal  0.3415826 0.1719309 1.0000000

    $Males
                 Math   Spatial    Verbal
    Math    1.0000000 0.3180507 0.3282856
    Spatial 0.3180507 1.0000000 0.1795489
    Verbal  0.3282856 0.1795489 1.0000000

### Stage 2 analysis

    ## Stage 2 analysis
    cluster2 <- tssem2(cluster1, RAM=RAM2, intervals="LB")
    summary(cluster2)

    $Females

    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|)
    Spatial2Math     0.409265        NA 0.305292 0.512536      NA       NA
    Verbal2Math      0.271217        NA 0.156125 0.386561      NA       NA
    SpatialCorVerbal 0.171931        NA 0.045061 0.298800      NA       NA

    Goodness-of-fit indices:
                                                Value
    Sample size                                235.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           105.57
    DF of independence model                     3.00
    RMSEA                                        0.00
    RMSEA lower 95% CI                           0.00
    RMSEA upper 95% CI                           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.)

    $Males

    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|)
    Spatial2Math     0.267739        NA 0.166906 0.368546      NA       NA
    Verbal2Math      0.280213        NA 0.178023 0.382428      NA       NA
    SpatialCorVerbal 0.179549        NA 0.069892 0.289197      NA       NA

    Goodness-of-fit indices:
                                                 Value
    Sample size                                303.000
    Chi-square of target model                   0.000
    DF of target model                           0.000
    p value of target model                      0.000
    Number of constraints imposed on "Smatrix"   0.000
    DF manually adjusted                         0.000
    Chi-square of independence model            68.564
    DF of independence model                     3.000
    RMSEA                                        0.000
    RMSEA lower 95% CI                           0.000
    RMSEA upper 95% CI                           0.000
    SRMR                                         0.000
    TLI                                           -Inf
    CFI                                          1.000
    AIC                                          0.000
    BIC                                          0.000
    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.

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

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

## Random-effects TSSEM

### Stage 1 analysis

- The \(I^2\) indices for the correlations between spatial and math, verbal and math, and spatial and verbal are .00, .81 and .23, respectively.

    #### Random-effects model
    ## Stage 1 analysis: A diagonal matrix for random effects
    random1 <- tssem1(Becker94$data, Becker94$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  Pr(>|z|)
    Intercept1  0.3777491  0.0395030  0.3003246  0.4551736  9.5625 < 2.2e-16 ***
    Intercept2  0.3807843  0.0784956  0.2269357  0.5346328  4.8510 1.228e-06 ***
    Intercept3  0.1704927  0.0513545  0.0698398  0.2711457  3.3199 0.0009004 ***
    Tau2_1_1    0.0005038  0.0042009 -0.0077298  0.0087374  0.1199 0.9045411
    Tau2_2_2    0.0416264  0.0257388 -0.0088206  0.0920734  1.6173 0.1058209
    Tau2_3_3    0.0067540  0.0102792 -0.0133928  0.0269008  0.6571 0.5111468
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Q statistic on the homogeneity of effect sizes: 61.02636
    Degrees of freedom of the Q statistic: 27
    P value of the Q statistic: 0.0001932118

    Heterogeneity indices (based on the estimated Tau2):
                                 Estimate
    Intercept1: I2 (Q statistic)   0.0337
    Intercept2: I2 (Q statistic)   0.7224
    Intercept3: I2 (Q statistic)   0.2676

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

    coef(random1, select="fixed")

    Intercept1 Intercept2 Intercept3
     0.3777491  0.3807843  0.1704927

    coef(random1, select="random")

        Tau2_1_1     Tau2_2_2     Tau2_3_3
    0.0005038018 0.0416263981 0.0067539586

### Stage 2 analysis

    ## Stage 2 analysis
    random2 <- tssem2(random1, RAM=RAM2, 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|)
    Spatial2Math      0.32219        NA 0.23725 0.40432      NA       NA
    Verbal2Math       0.32585        NA 0.16874 0.48257      NA       NA
    SpatialCorVerbal  0.17049        NA 0.06984 0.27115      NA       NA

    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           110.81
    DF of independence model                     3.00
    RMSEA                                        0.00
    RMSEA lower 95% CI                           0.00
    RMSEA upper 95% CI                           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.)

- We may also plot the models with the labels or the parameter estimates.

    ## Plot the model with labels
    plot(random2, whatLabels="path", color="red")

    ## Plot the parameter estimates
    plot(random2, color="green")

# One-stage meta-analytic structural equation modeling (OSMASEM)

## Nohe15

- Data: Nohe, C., Meier, L. L., Sonntag, K., & Michel, A. (2015). The chicken or the egg? A meta-analysis of panel studies of the relationship between work–family conflict and strain. Journal of Applied Psychology, 100(2), 522–536. https://doi.org/10.1037/a0038012

### Data preparation

    ## Convert correlation matrices into a dataframe
    my.df <- Cor2DataFrame(Nohe15A1$data, Nohe15A1$n)

    ## Add the standardized Lag (moderator) to the data
    my.df$data <- data.frame(my.df$data, Lag=scale(Nohe15A1$Lag), check.names=FALSE)
    head(my.df$data)

                                           S1_W1 W2_W1 S2_W1 W2_S1 S2_S1 S2_W2
    Britt...Dawson..2005.                   0.29  0.58  0.22  0.24  0.57  0.27
    Demerouti.et.al...2004.                 0.53  0.57  0.41  0.41  0.68  0.54
    Ford..2010.                             0.35  0.75  0.32  0.26  0.74  0.30
    Hammer.et.al...2005...female.subsample  0.32  0.57  0.22  0.30  0.43  0.30
    Hammer.et.al...2005...male.subsample    0.19  0.54  0.17  0.21  0.60  0.30
    Innstrand.et.al...2008.                 0.42  0.63  0.31  0.30  0.62  0.44
                                           C(S1_W1 S1_W1) C(W2_W1 S1_W1)
    Britt...Dawson..2005.                    0.0014224768   2.033620e-04
    Demerouti.et.al...2004.                  0.0020763976   2.968524e-04
    Ford..2010.                              0.0021207071   3.031847e-04
    Hammer.et.al...2005...female.subsample   0.0029726186   4.249795e-04
    Hammer.et.al...2005...male.subsample     0.0029726152   4.249774e-04
    Innstrand.et.al...2008.                  0.0003112269   4.449439e-05
                                           C(S2_W1 S1_W1) C(W2_S1 S1_W1)
    Britt...Dawson..2005.                    0.0008791213   0.0008736583
    Demerouti.et.al...2004.                  0.0012832625   0.0012752882
    Ford..2010.                              0.0013106441   0.0013025001
    Hammer.et.al...2005...female.subsample   0.0018371479   0.0018257305
    Hammer.et.al...2005...male.subsample     0.0018371437   0.0018257263
    Innstrand.et.al...2008.                  0.0001923453   0.0001911500
                                           C(S2_S1 S1_W1) C(S2_W2 S1_W1)
    Britt...Dawson..2005.                    2.047329e-04   0.0004890867
    Demerouti.et.al...2004.                  2.988530e-04   0.0007139274
    Ford..2010.                              3.052294e-04   0.0007291591
    Hammer.et.al...2005...female.subsample   4.278443e-04   0.0010220746
    Hammer.et.al...2005...male.subsample     4.278412e-04   0.0010220712
    Innstrand.et.al...2008.                  4.479435e-05   0.0001070088
                                           C(W2_W1 W2_W1) C(S2_W1 W2_W1)
    Britt...Dawson..2005.                    0.0007981101   3.750373e-04
    Demerouti.et.al...2004.                  0.0011650046   5.474480e-04
    Ford..2010.                              0.0011898662   5.591278e-04
    Hammer.et.al...2005...female.subsample   0.0016678476   7.837379e-04
    Hammer.et.al...2005...male.subsample     0.0016678468   7.837358e-04
    Innstrand.et.al...2008.                  0.0001746202   8.205568e-05
                                           C(W2_S1 W2_W1) C(S2_S1 W2_W1)
    Britt...Dawson..2005.                    3.671002e-04   1.042801e-04
    Demerouti.et.al...2004.                  5.358621e-04   1.522206e-04
    Ford..2010.                              5.472951e-04   1.554681e-04
    Hammer.et.al...2005...female.subsample   7.671507e-04   2.179215e-04
    Hammer.et.al...2005...male.subsample     7.671487e-04   2.179201e-04
    Innstrand.et.al...2008.                  8.031902e-05   2.281597e-05
                                           C(S2_W2 W2_W1) C(S2_W1 S2_W1)
    Britt...Dawson..2005.                    2.035363e-04   0.0016496782
    Demerouti.et.al...2004.                  2.971074e-04   0.0024080460
    Ford..2010.                              3.034447e-04   0.0024594317
    Hammer.et.al...2005...female.subsample   4.253442e-04   0.0034474152
    Hammer.et.al...2005...male.subsample     4.253427e-04   0.0034474108
    Innstrand.et.al...2008.                  4.453258e-05   0.0003609371
                                           C(W2_S1 S2_W1) C(S2_S1 S2_W1)
    Britt...Dawson..2005.                    0.0005769062   0.0003592754
    Demerouti.et.al...2004.                  0.0008421211   0.0005244397
    Ford..2010.                              0.0008600869   0.0005356297
    Hammer.et.al...2005...female.subsample   0.0012055975   0.0007507994
    Hammer.et.al...2005...male.subsample     0.0012055924   0.0007507962
    Innstrand.et.al...2008.                  0.0001262232   0.0000786071
                                           C(S2_W2 S2_W1) C(W2_S1 W2_S1)
    Britt...Dawson..2005.                    0.0008607757   0.0016601375
    Demerouti.et.al...2004.                  0.0012564833   0.0024233139
    Ford..2010.                              0.0012832933   0.0024750262
    Hammer.et.al...2005...female.subsample   0.0017988102   0.0034692713
    Hammer.et.al...2005...male.subsample     0.0017988067   0.0034692668
    Innstrand.et.al...2008.                  0.0001883315   0.0003632253
                                           C(S2_S1 W2_S1) C(S2_W2 W2_S1)
    Britt...Dawson..2005.                    3.727395e-04   0.0008738993
    Demerouti.et.al...2004.                  5.440932e-04   0.0012756405
    Ford..2010.                              5.557032e-04   0.0013028593
    Hammer.et.al...2005...female.subsample   7.789350e-04   0.0018262350
    Hammer.et.al...2005...male.subsample     7.789319e-04   0.0018262313
    Innstrand.et.al...2008.                  8.155286e-05   0.0001912028
                                           C(S2_S1 S2_S1) C(S2_W2 S2_S1)
    Britt...Dawson..2005.                    0.0007805631   1.951848e-04
    Demerouti.et.al...2004.                  0.0011393911   2.849157e-04
    Ford..2010.                              0.0011637067   2.909942e-04
    Hammer.et.al...2005...female.subsample   0.0016311787   4.078912e-04
    Hammer.et.al...2005...male.subsample     0.0016311776   4.078885e-04
    Innstrand.et.al...2008.                  0.0001707811   4.270531e-05
                                           C(S2_W2 S2_W2)        Lag
    Britt...Dawson..2005.                    0.0013981125 -0.6794521
    Demerouti.et.al...2004.                  0.0020408333 -0.7711151
    Ford..2010.                              0.0020843835 -0.8016694
    Hammer.et.al...2005...female.subsample   0.0029217041 -0.1294740
    Hammer.et.al...2005...male.subsample     0.0029217016 -0.1294740
    Innstrand.et.al...2008.                  0.0003058963  0.6038301

    ## Display the pairwise no. of studies
    pattern.na(my.df, show.na=FALSE, type="osmasem")

          S1_W1 W2_W1 S2_W1 W2_S1 S2_S1 S2_W2
    S1_W1    32    32    32    32    32    32
    W2_W1    32    32    32    32    32    32
    S2_W1    32    32    32    32    32    32
    W2_S1    32    32    32    32    32    32
    S2_S1    32    32    32    32    32    32
    S2_W2    32    32    32    32    32    32

    ## Proposed model
    model3 <- 'W2 ~ w2w*W1 + s2w*S1
               S2 ~ w2s*W1 + s2s*S1
               W1 ~~ w1WITHs1*S1
               W2 ~~ w2WITHs2*S2
               W1 ~~ 1*W1
               S1 ~~ 1*S1
               W2 ~~ Errw2*W2
               S2 ~~ Errs2*S2'

    ## Plot the model
    plot(model3, col="yellow", layout="spring")

    ## Convert the lavaan syntax into the RAM specification
    RAM3 <- lavaan2RAM(model3, obs.variables=c("W1", "S1", "W2", "S2"))
    RAM3

    $A
       W1      S1      W2  S2
    W1 "0"     "0"     "0" "0"
    S1 "0"     "0"     "0" "0"
    W2 "0*w2w" "0*s2w" "0" "0"
    S2 "0*w2s" "0*s2s" "0" "0"

    $S
       W1           S1           W2           S2
    W1 "1"          "0*w1WITHs1" "0"          "0"
    S1 "0*w1WITHs1" "1"          "0"          "0"
    W2 "0"          "0"          "0*Errw2"    "0*w2WITHs2"
    S2 "0"          "0"          "0*w2WITHs2" "0*Errs2"

    $F
       W1 S1 W2 S2
    W1  1  0  0  0
    S1  0  1  0  0
    W2  0  0  1  0
    S2  0  0  0  1

    $M
      W1 S1 W2 S2
    1  0  0  0  0

### Fitting a just-identified model (df=0)

    # Fit the OSMASEM
    fit0 <- osmasem(model.name="No moderator", RAM=RAM3, data=my.df)
    summary(fit0)

    Summary of No moderator

    free parameters:
           name  matrix row col    Estimate  Std.Error A    z value     Pr(>|z|)
    1       w2w      A0  W2  W1  0.57247128 0.02226456    25.712224 0.000000e+00
    2       w2s      A0  S2  W1  0.08023683 0.02484212     3.229870 1.238464e-03
    3       s2w      A0  W2  S1  0.08584126 0.02479588     3.461916 5.363452e-04
    4       s2s      A0  S2  S1  0.58612398 0.02078999    28.192599 0.000000e+00
    5  w1WITHs1      S0  S1  W1  0.38045214 0.02256155    16.862856 0.000000e+00
    6  w2WITHs2      S0  S2  W2  0.16888457 0.02523190     6.693296 2.181988e-11
    7    Tau1_1 vecTau1   1   1 -2.15335935 0.14358424   -14.997185 0.000000e+00
    8    Tau1_2 vecTau1   2   1 -2.36882282 0.14419317   -16.428121 0.000000e+00
    9    Tau1_3 vecTau1   3   1 -2.47296602 0.15796636   -15.655017 0.000000e+00
    10   Tau1_4 vecTau1   4   1 -2.47676250 0.15669553   -15.806210 0.000000e+00
    11   Tau1_5 vecTau1   5   1 -2.46245722 0.14519715   -16.959405 0.000000e+00
    12   Tau1_6 vecTau1   6   1 -2.19983872 0.14187319   -15.505669 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             12                    180             -300.1702
       Saturated:             27                    165                    NA
    Independence:             12                    180                    NA
    Number of observations/statistics: 12906/192

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -660.1702              -276.1702                -276.1460
    BIC:     -2003.9507              -186.5848                -224.7196
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:38
    Wall clock time: 0.2403653 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Plot the fitted model
    plot(fit0, layout="spring", color="yellow")

    ## Get the estimated coefficients
    coef(fit0)

           w2w        w2s        s2w        s2s   w1WITHs1   w2WITHs2
    0.57247128 0.08023683 0.08584126 0.58612398 0.38045214 0.16888457

    ## A matrix
    mxEval(Amatrix, fit0$mx.fit)

               [,1]       [,2] [,3] [,4]
    [1,] 0.00000000 0.00000000    0    0
    [2,] 0.00000000 0.00000000    0    0
    [3,] 0.57247128 0.08584126    0    0
    [4,] 0.08023683 0.58612398    0    0

    ## S matrix
    mxEval(Smatrix, fit0$mx.fit)

              [,1]      [,2]      [,3]      [,4]
    [1,] 1.0000000 0.3804521 0.0000000 0.0000000
    [2,] 0.3804521 1.0000000 0.0000000 0.0000000
    [3,] 0.0000000 0.0000000 0.6275158 0.1688846
    [4,] 0.0000000 0.0000000 0.1688846 0.6142364

    ## Extract the heterogeneity variance-covariance matrix
    VarCorr(fit0)

              Tau2_1      Tau2_2      Tau2_3      Tau2_4      Tau2_5    Tau2_6
    Tau2_1 0.0134777 0.000000000 0.000000000 0.000000000 0.000000000 0.0000000
    Tau2_2 0.0000000 0.008759244 0.000000000 0.000000000 0.000000000 0.0000000
    Tau2_3 0.0000000 0.000000000 0.007112283 0.000000000 0.000000000 0.0000000
    Tau2_4 0.0000000 0.000000000 0.000000000 0.007058484 0.000000000 0.0000000
    Tau2_5 0.0000000 0.000000000 0.000000000 0.000000000 0.007263348 0.0000000
    Tau2_6 0.0000000 0.000000000 0.000000000 0.000000000 0.000000000 0.0122813

### Using `Lag` as a moderator on the A matrix

    A1 <- create.modMatrix(RAM=RAM3, output="A", mod="Lag")
    A1

       W1           S1           W2  S2
    W1 "0"          "0"          "0" "0"
    S1 "0"          "0"          "0" "0"
    W2 "0*data.Lag" "0*data.Lag" "0" "0"
    S2 "0*data.Lag" "0*data.Lag" "0" "0"

    fit1 <- osmasem(model.name="Lag as moderator on Ax", RAM=RAM3, Ax=A1, data=my.df)
    summary(fit1)

    Summary of Lag as moderator on Ax

    free parameters:
           name  matrix row col     Estimate  Std.Error A     z value     Pr(>|z|)
    1       w2w      A0  W2  W1  0.573039754 0.01839765    31.1474450 0.000000e+00
    2       w2s      A0  S2  W1  0.079844957 0.02419487     3.3000773 9.665820e-04
    3       s2w      A0  W2  S1  0.085391029 0.02393611     3.5674563 3.604635e-04
    4       s2s      A0  S2  S1  0.586234240 0.01962560    29.8708936 0.000000e+00
    5  w1WITHs1      S0  S1  W1  0.381183684 0.02282276    16.7019094 0.000000e+00
    6  w2WITHs2      S0  S2  W2  0.166974807 0.02500437     6.6778237 2.425171e-11
    7     w2w_1      A1  W2  W1 -0.062015581 0.01850573    -3.3511558 8.047503e-04
    8     w2s_1      A1  S2  W1 -0.025933696 0.02096723    -1.2368680 2.161361e-01
    9     s2w_1      A1  W2  S1 -0.002382794 0.02055896    -0.1159005 9.077314e-01
    10    s2s_1      A1  S2  S1 -0.027809752 0.01974171    -1.4086798 1.589299e-01
    11   Tau1_1 vecTau1   1   1 -2.138191073 0.14360107   -14.8897993 0.000000e+00
    12   Tau1_2 vecTau1   2   1 -2.630518388 0.16155260   -16.2827362 0.000000e+00
    13   Tau1_3 vecTau1   3   1 -2.524194455 0.16007545   -15.7687793 0.000000e+00
    14   Tau1_4 vecTau1   4   1 -2.519909061 0.15983967   -15.7652298 0.000000e+00
    15   Tau1_5 vecTau1   5   1 -2.537475873 0.14712279   -17.2473344 0.000000e+00
    16   Tau1_6 vecTau1   6   1 -2.198863896 0.14144102   -15.5461544 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             16                    176             -323.6921
       Saturated:             27                    165                    NA
    Independence:             12                    180                    NA
    Number of observations/statistics: 12906/192

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -675.6921              -291.6921                -291.6499
    BIC:     -1989.6109              -172.2450                -223.0914
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:38
    Wall clock time: 0.398262 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Get the estimated coefficients
    coef(fit1)

             w2w          w2s          s2w          s2s     w1WITHs1     w2WITHs2
     0.573039754  0.079844957  0.085391029  0.586234240  0.381183684  0.166974807
           w2w_1        w2s_1        s2w_1        s2s_1
    -0.062015581 -0.025933696 -0.002382794 -0.027809752

    ## Extract the residual heterogeneity variance-covariance matrix
    VarCorr(fit1)

               Tau2_1      Tau2_2      Tau2_3      Tau2_4      Tau2_5     Tau2_6
    Tau2_1 0.01389283 0.000000000 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_2 0.00000000 0.005189921 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_3 0.00000000 0.000000000 0.006419668 0.000000000 0.000000000 0.00000000
    Tau2_4 0.00000000 0.000000000 0.000000000 0.006474926 0.000000000 0.00000000
    Tau2_5 0.00000000 0.000000000 0.000000000 0.000000000 0.006251388 0.00000000
    Tau2_6 0.00000000 0.000000000 0.000000000 0.000000000 0.000000000 0.01230527

    ## Calculate the R2
    osmasemR2(fit1, fit0)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.013477701 0.008759244 0.007112283 0.007058484 0.007263348 0.012281301

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.013892834 0.005189921 0.006419668 0.006474926 0.006251388 0.012305268

    $R2
      Tau2_1_1   Tau2_2_2   Tau2_3_3   Tau2_4_4   Tau2_5_5   Tau2_6_6
    0.00000000 0.40749214 0.09738292 0.08267469 0.13932412 0.00000000

    ## Compare the models with and without the moderator
    anova(fit1, fit0)

                        base   comparison ep  minus2LL  df       AIC   diffLL
    1 Lag as moderator on Ax         <NA> 16 -323.6921 176 -675.6921       NA
    2 Lag as moderator on Ax No moderator 12 -300.1702 180 -660.1702 23.52198
      diffdf            p
    1     NA           NA
    2      4 9.957542e-05

### Using `Lag` as a moderator on the S matrix

    S1 <- create.modMatrix(RAM=RAM3, output="S", mod="Lag")
    S1

       W1           S1           W2           S2
    W1 "0"          "0*data.Lag" "0"          "0"
    S1 "0*data.Lag" "0"          "0"          "0"
    W2 "0"          "0"          "0"          "0*data.Lag"
    S2 "0"          "0"          "0*data.Lag" "0"

    fit2 <- osmasem(model.name="Lag as moderator on Sx", RAM=RAM3, Sx=S1, data=my.df)
    summary(fit2)

    Summary of Lag as moderator on Sx

    free parameters:
             name  matrix row col     Estimate  Std.Error A     z value
    1         w2w      A0  W2  W1  0.569234670 0.02189962    25.9928959
    2         w2s      A0  S2  W1  0.085032454 0.02395174     3.5501572
    3         s2w      A0  W2  S1  0.092336298 0.02405803     3.8380655
    4         s2s      A0  S2  S1  0.584608176 0.02049270    28.5276318
    5    w1WITHs1      S0  S1  W1  0.376986112 0.02230556    16.9009918
    6    w2WITHs2      S0  S2  W2  0.165484765 0.02477816     6.6786535
    7  w1WITHs1_1      S1  S1  W1 -0.053146775 0.01794364    -2.9618720
    8  w2WITHs2_1      S1  S2  W2 -0.008169624 0.02172729    -0.3760074
    9      Tau1_1 vecTau1   1   1 -2.169523120 0.14503078   -14.9590532
    10     Tau1_2 vecTau1   2   1 -2.376926159 0.14470313   -16.4262249
    11     Tau1_3 vecTau1   3   1 -2.524376081 0.15992215   -15.7850309
    12     Tau1_4 vecTau1   4   1 -2.516292980 0.15999135   -15.7276818
    13     Tau1_5 vecTau1   5   1 -2.464565302 0.14547984   -16.9409410
    14     Tau1_6 vecTau1   6   1 -2.210555114 0.14170819   -15.5993466
           Pr(>|z|)
    1  0.000000e+00
    2  3.850012e-04
    3  1.240074e-04
    4  0.000000e+00
    5  0.000000e+00
    6  2.411471e-11
    7  3.057749e-03
    8  7.069114e-01
    9  0.000000e+00
    10 0.000000e+00
    11 0.000000e+00
    12 0.000000e+00
    13 0.000000e+00
    14 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             14                    178             -309.5339
       Saturated:             27                    165                    NA
    Independence:             12                    180                    NA
    Number of observations/statistics: 12906/192

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -665.5339              -281.5339                -281.5013
    BIC:     -1994.3835              -177.0176                -221.5082
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:39
    Wall clock time: 0.3515871 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Get the estimated coefficients
    coef(fit2)

             w2w          w2s          s2w          s2s     w1WITHs1     w2WITHs2
     0.569234670  0.085032454  0.092336298  0.584608176  0.376986112  0.165484765
      w1WITHs1_1   w2WITHs2_1
    -0.053146775 -0.008169624

    ## Extract the residual heterogeneity variance-covariance matrix
    VarCorr(fit2)

               Tau2_1     Tau2_2      Tau2_3      Tau2_4      Tau2_5     Tau2_6
    Tau2_1 0.01304897 0.00000000 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_2 0.00000000 0.00861843 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_3 0.00000000 0.00000000 0.006417336 0.000000000 0.000000000 0.00000000
    Tau2_4 0.00000000 0.00000000 0.000000000 0.006521923 0.000000000 0.00000000
    Tau2_5 0.00000000 0.00000000 0.000000000 0.000000000 0.007232789 0.00000000
    Tau2_6 0.00000000 0.00000000 0.000000000 0.000000000 0.000000000 0.01202088

    ## Calculate the R2
    osmasemR2(fit2, fit0)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.013477701 0.008759244 0.007112283 0.007058484 0.007263348 0.012281301

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.013048968 0.008618430 0.006417336 0.006521923 0.007232789 0.012020879

    $R2
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.031810588 0.016076056 0.097710736 0.076016399 0.004207279 0.021204745

    ## Compare the models with and without the moderator
    anova(fit2, fit0)

                        base   comparison ep  minus2LL  df       AIC  diffLL diffdf
    1 Lag as moderator on Sx         <NA> 14 -309.5339 178 -665.5339      NA     NA
    2 Lag as moderator on Sx No moderator 12 -300.1702 180 -660.1702 9.36368      2
                p
    1          NA
    2 0.009261957

### Using `Lag` as a moderator on both the A and S matrices

    fit3 <- osmasem(model.name="Lag as moderator on Ax and Sx",
                    RAM=RAM3, Ax=A1, Sx=S1, data=my.df)
    summary(fit3)

    Summary of Lag as moderator on Ax and Sx

    free parameters:
             name  matrix row col     Estimate  Std.Error A     z value
    1         w2w      A0  W2  W1  0.574425646 0.01843827    31.1539875
    2         w2s      A0  S2  W1  0.079662705 0.02387730     3.3363362
    3         s2w      A0  W2  S1  0.083087383 0.02363556     3.5153550
    4         s2s      A0  S2  S1  0.585858041 0.01964429    29.8233211
    5    w1WITHs1      S0  S1  W1  0.381044061 0.02197099    17.3430504
    6    w2WITHs2      S0  S2  W2  0.167582968 0.02478852     6.7605063
    7       w2w_1      A1  W2  W1 -0.063025034 0.01656924    -3.8037377
    8       w2s_1      A1  S2  W1 -0.013043444 0.02072888    -0.6292402
    9       s2w_1      A1  W2  S1  0.006865330 0.02021927     0.3395440
    10      s2s_1      A1  S2  S1 -0.031632455 0.01771820    -1.7853084
    11 w1WITHs1_1      S1  S1  W1 -0.038859418 0.02182853    -1.7802124
    12 w2WITHs2_1      S1  S2  W2  0.007050356 0.02335410     0.3018895
    13     Tau1_1 vecTau1   1   1 -2.184627864 0.14427007   -15.1426274
    14     Tau1_2 vecTau1   2   1 -2.633170334 0.16172683   -16.2815925
    15     Tau1_3 vecTau1   3   1 -2.531970054 0.16055055   -15.7705472
    16     Tau1_4 vecTau1   4   1 -2.534620769 0.16068834   -15.7735198
    17     Tau1_5 vecTau1   5   1 -2.539156303 0.14724186   -17.2447994
    18     Tau1_6 vecTau1   6   1 -2.205765683 0.14170877   -15.5654843
           Pr(>|z|)
    1  0.000000e+00
    2  8.489044e-04
    3  4.391665e-04
    4  0.000000e+00
    5  0.000000e+00
    6  1.375100e-11
    7  1.425291e-04
    8  5.291918e-01
    9  7.342000e-01
    10 7.421130e-02
    11 7.504120e-02
    12 7.627363e-01
    13 0.000000e+00
    14 0.000000e+00
    15 0.000000e+00
    16 0.000000e+00
    17 0.000000e+00
    18 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             18                    174             -327.0718
       Saturated:             27                    165                    NA
    Independence:             12                    180                    NA
    Number of observations/statistics: 12906/192

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -675.0718              -291.0718                -291.0187
    BIC:     -1974.0597              -156.6938                -213.8959
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:39
    Wall clock time: 0.4757283 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    coef(fit3)

             w2w          w2s          s2w          s2s     w1WITHs1     w2WITHs2
     0.574425646  0.079662705  0.083087383  0.585858041  0.381044061  0.167582968
           w2w_1        w2s_1        s2w_1        s2s_1   w1WITHs1_1   w2WITHs2_1
    -0.063025034 -0.013043444  0.006865330 -0.031632455 -0.038859418  0.007050356

    VarCorr(fit3)

               Tau2_1      Tau2_2      Tau2_3      Tau2_4      Tau2_5     Tau2_6
    Tau2_1 0.01266066 0.000000000 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_2 0.00000000 0.005162467 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_3 0.00000000 0.000000000 0.006320607 0.000000000 0.000000000 0.00000000
    Tau2_4 0.00000000 0.000000000 0.000000000 0.006287187 0.000000000 0.00000000
    Tau2_5 0.00000000 0.000000000 0.000000000 0.000000000 0.006230413 0.00000000
    Tau2_6 0.00000000 0.000000000 0.000000000 0.000000000 0.000000000 0.01213658

    osmasemR2(fit3, fit0)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.013477701 0.008759244 0.007112283 0.007058484 0.007263348 0.012281301

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.012660660 0.005162467 0.006320607 0.006287187 0.006230413 0.012136579

    $R2
      Tau2_1_1   Tau2_2_2   Tau2_3_3   Tau2_4_4   Tau2_5_5   Tau2_6_6
    0.06062172 0.41062642 0.11131111 0.10927232 0.14221188 0.01178395

    anova(fit3, fit0)

                               base   comparison ep  minus2LL  df       AIC
    1 Lag as moderator on Ax and Sx         <NA> 18 -327.0718 174 -675.0718
    2 Lag as moderator on Ax and Sx No moderator 12 -300.1702 180 -660.1702
        diffLL diffdf            p
    1       NA     NA           NA
    2 26.90165      6 0.0001510812

## Roorda11

- Data: Roorda, D. L., Koomen, H. M. Y., Spilt, J. L., & Oort, F. J. (2011). The influence of affective teacher-student relationships on students’ school engagement and achievement a meta-analytic approach. Review of Educational Research, 81(4), 493-529.

### Data preparation

    my.df <- Cor2DataFrame(Roorda11$data, Roorda11$n)

    ## Add centered SES as a moderator (standardizing SES helps to improve the stability of the results)
    my.df$data <- data.frame(my.df$data, SES=scale(Roorda11$SES), check.names=FALSE)
    head(my.df$data)

      neg_pos enga_pos achiev_pos enga_neg achiev_neg achiev_enga
    1   -0.54       NA       0.18       NA      -0.29          NA
    2      NA     0.64       0.29       NA         NA        0.23
    3      NA     0.29         NA       NA         NA          NA
    4      NA     0.29         NA       NA         NA          NA
    5      NA       NA         NA       NA      -0.24          NA
    6      NA     0.06      -0.09       NA         NA        0.20
      C(neg_pos neg_pos) C(enga_pos neg_pos) C(achiev_pos neg_pos)
    1        0.000598950       -0.0001515536         -9.072864e-05
    2        0.001837528       -0.0004649533         -2.783479e-04
    3        0.006379060       -0.0016141045         -9.662965e-04
    4        0.011888252       -0.0030081066         -1.800830e-03
    5        0.004383376       -0.0011091328         -6.639929e-04
    6        0.008436836       -0.0021347947         -1.278014e-03
      C(enga_neg neg_pos) C(achiev_neg neg_pos) C(achiev_enga neg_pos)
    1        0.0001507644          3.883486e-05          -2.214629e-05
    2        0.0004625319          1.191421e-04          -6.794345e-05
    3        0.0016056989          4.136080e-04          -2.358654e-04
    4        0.0029924421          7.708164e-04          -4.395731e-04
    5        0.0011033565          2.842102e-04          -1.620764e-04
    6        0.0021236752          5.470317e-04          -3.119591e-04
      C(enga_pos enga_pos) C(achiev_pos enga_pos) C(enga_neg enga_pos)
    1         0.0006397572           0.0001647858        -0.0001910375
    2         0.0019627200           0.0005055490        -0.0005860855
    3         0.0068136688           0.0017550373        -0.0020346178
    4         0.0126982034           0.0032707579        -0.0037917943
    5         0.0046820179           0.0012059757        -0.0013980893
    6         0.0090116374           0.0023211804        -0.0026909598
      C(achiev_neg enga_pos) C(achiev_enga enga_pos) C(achiev_pos achiev_pos)
    1          -5.134717e-05            3.127032e-05             0.0007526508
    2          -1.575275e-04            9.593538e-05             0.0023090669
    3          -5.468632e-04            3.330453e-04             0.0080160312
    4          -1.019159e-03            6.206760e-04             0.0149389747
    5          -3.757785e-04            2.288513e-04             0.0055082245
    6          -7.232771e-04            4.404744e-04             0.0106018629
      C(enga_neg achiev_pos) C(achiev_neg achiev_pos) C(achiev_enga achiev_pos)
    1          -6.395578e-05            -0.0002463391              0.0001998713
    2          -1.962098e-04            -0.0007557443              0.0006131868
    3          -6.811476e-04            -0.0026236014              0.0021287066
    4          -1.269423e-03            -0.0048894512              0.0039671428
    5          -4.680534e-04            -0.0018028119              0.0014627429
    6          -9.008851e-04            -0.0034699401              0.0028153918
      C(enga_neg enga_neg) C(achiev_neg enga_neg) C(achiev_enga enga_neg)
    1         0.0006389864           0.0001558886           -7.403054e-05
    2         0.0019603547           0.0004782522           -2.271188e-04
    3         0.0068054518           0.0016602714           -7.884458e-04
    4         0.0126829005           0.0030941499           -1.469389e-03
    5         0.0046763738           0.0011408581           -5.417835e-04
    6         0.0090007789           0.0021958510           -1.042796e-03
      C(achiev_neg achiev_neg) C(achiev_enga achiev_neg) C(achiev_enga achiev_enga)
    1             0.0007298084             -0.0001921435               0.0006716402
    2             0.0022389876             -0.0005894775               0.0020605339
    3             0.0077727451             -0.0020463951               0.0071532324
    4             0.0144855852             -0.0038137535               0.0133310385
    5             0.0053410506             -0.0014061840               0.0049153508
    6             0.0102800992             -0.0027065365               0.0094607378
             SES
    1  0.4464188
    2  0.7145383
    3  0.8821129
    4 -0.8941783
    5 -0.9947231
    6 -0.5925439

    ## Display the pairwise no. of studies
    pattern.na(my.df, show.na=FALSE, type="osmasem")

                neg_pos enga_pos achiev_pos enga_neg achiev_neg achiev_enga
    neg_pos          15       30         28       17         20          24
    enga_pos         30       21         37       23         35          24
    achiev_pos       28       37         26       30         31          29
    enga_neg         17       23         30        8         21          16
    achiev_neg       20       35         31       21         18          26
    achiev_enga      24       24         29       16         26          13

    ## Proposed model
    model4 <- 'enga ~ b31*neg + b32*pos
               achiev ~ b41*neg + b42*pos + b43*enga
               pos ~~ p21*neg
               pos ~~ 1*pos
               neg ~~ 1*neg
               enga ~~ p33*enga
               achiev ~~ p44*achiev'

    plot(model4, color="yellow")

    ## Convert the lavaan model to the RAM specification
    RAM4 <- lavaan2RAM(model4, obs.variables=c("neg", "pos", "enga", "achiev"))
    RAM4

    $A
           neg     pos     enga    achiev
    neg    "0"     "0"     "0"     "0"
    pos    "0"     "0"     "0"     "0"
    enga   "0*b31" "0*b32" "0"     "0"
    achiev "0*b41" "0*b42" "0*b43" "0"

    $S
           neg     pos     enga    achiev
    neg    "1"     "0*p21" "0"     "0"
    pos    "0*p21" "1"     "0"     "0"
    enga   "0"     "0"     "0*p33" "0"
    achiev "0"     "0"     "0"     "0*p44"

    $F
           neg pos enga achiev
    neg      1   0    0      0
    pos      0   1    0      0
    enga     0   0    1      0
    achiev   0   0    0      1

    $M
      neg pos enga achiev
    1   0   0    0      0

### Fitting a just-identified model (df=0)

    mx.fit0a <- osmasem(model.name="Just identified model", RAM=RAM4, data= my.df)
    summary(mx.fit0a)

    Summary of Just identified model

    free parameters:
         name  matrix    row  col    Estimate  Std.Error A    z value     Pr(>|z|)
    1     b31      A0   enga  neg  0.25513471 0.04084654     6.246177 4.206211e-10
    2     b41      A0 achiev  neg  0.04088993 0.02704882     1.511708 1.306081e-01
    3     b32      A0   enga  pos -0.24326592 0.04795947    -5.072323 3.929891e-07
    4     b42      A0 achiev  pos -0.09342504 0.03205930    -2.914132 3.566786e-03
    5     b43      A0 achiev enga  0.23653110 0.04498219     5.258328 1.453714e-07
    6     p21      S0    pos  neg -0.23871423 0.04044945    -5.901545 3.601131e-09
    7  Tau1_1 vecTau1      1    1 -1.97354078 0.21976079    -8.980404 0.000000e+00
    8  Tau1_2 vecTau1      2    1 -1.87778194 0.18152775   -10.344325 0.000000e+00
    9  Tau1_3 vecTau1      3    1 -2.55758768 0.22104151   -11.570622 0.000000e+00
    10 Tau1_4 vecTau1      4    1 -2.21995381 0.31525223    -7.041834 1.897371e-12
    11 Tau1_5 vecTau1      5    1 -2.50426684 0.26585163    -9.419791 0.000000e+00
    12 Tau1_6 vecTau1      6    1 -2.11922097 0.23882306    -8.873603 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             12                     89             -127.1991
       Saturated:             27                     74                    NA
    Independence:             12                     89                    NA
    Number of observations/statistics: 29438/101

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -305.1991            -103.199109               -103.18851
    BIC:     -1043.0128              -3.718609                -41.85444
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:41
    Wall clock time: 0.2561717 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Plot the fitted model
    plot(mx.fit0a, layout="spring", color="yellow")

    coef(mx.fit0a)

            b31         b41         b32         b42         b43         p21
     0.25513471  0.04088993 -0.24326592 -0.09342504  0.23653110 -0.23871423

    ## Extract the variance component
    VarCorr(mx.fit0a)

               Tau2_1     Tau2_2      Tau2_3     Tau2_4      Tau2_5     Tau2_6
    Tau2_1 0.01931098 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000
    Tau2_2 0.00000000 0.02338726 0.000000000 0.00000000 0.000000000 0.00000000
    Tau2_3 0.00000000 0.00000000 0.006004925 0.00000000 0.000000000 0.00000000
    Tau2_4 0.00000000 0.00000000 0.000000000 0.01179703 0.000000000 0.00000000
    Tau2_5 0.00000000 0.00000000 0.000000000 0.00000000 0.006680692 0.00000000
    Tau2_6 0.00000000 0.00000000 0.000000000 0.00000000 0.000000000 0.01443006

### Fitting an over-identified model with b41=b42=0 (df=2))

    ## Proposed model
    model5 <- 'enga ~ b31*neg + b32*pos
               achiev ~ 0*neg + 0*pos + b43*enga
               pos ~~ p21*neg
               pos ~~ 1*pos
               neg ~~ 1*neg
               enga ~~ p33*enga
               achiev ~~ p44*achiev'

    plot(model5, layout="tree", color="yellow")

    RAM5 <- lavaan2RAM(model5, obs.variables=c("neg", "pos", "enga", "achiev"))
    RAM5

    $A
           neg     pos     enga    achiev
    neg    "0"     "0"     "0"     "0"
    pos    "0"     "0"     "0"     "0"
    enga   "0*b31" "0*b32" "0"     "0"
    achiev "0"     "0"     "0*b43" "0"

    $S
           neg     pos     enga    achiev
    neg    "1"     "0*p21" "0"     "0"
    pos    "0*p21" "1"     "0"     "0"
    enga   "0"     "0"     "0*p33" "0"
    achiev "0"     "0"     "0"     "0*p44"

    $F
           neg pos enga achiev
    neg      1   0    0      0
    pos      0   1    0      0
    enga     0   0    1      0
    achiev   0   0    0      1

    $M
      neg pos enga achiev
    1   0   0    0      0

    mx.fit0b <- osmasem(model.name="Over identified model", RAM=RAM5, data= my.df)

    ## Get the summary and test statistics
    summary(mx.fit0b, fitIndices=TRUE)

    Summary of Over identified model

    free parameters:
         name  matrix    row  col   Estimate  Std.Error A    z value     Pr(>|z|)
    1     b31      A0   enga  neg  0.2649761 0.03780500     7.009023 2.399858e-12
    2     b32      A0   enga  pos -0.2964373 0.05474967    -5.414413 6.149005e-08
    3     b43      A0 achiev enga  0.3500633 0.04052992     8.637157 0.000000e+00
    4     p21      S0    pos  neg -0.2386957 0.04057830    -5.882348 4.044864e-09
    5  Tau1_1 vecTau1      1    1 -1.9694734 0.21977208    -8.961436 0.000000e+00
    6  Tau1_2 vecTau1      2    1 -1.8706614 0.18373396   -10.181359 0.000000e+00
    7  Tau1_3 vecTau1      3    1 -2.5635350 0.22701101   -11.292558 0.000000e+00
    8  Tau1_4 vecTau1      4    1 -2.0996828 0.34898874    -6.016477 1.782536e-09
    9  Tau1_5 vecTau1      5    1 -2.3799825 0.27470083    -8.663907 0.000000e+00
    10 Tau1_6 vecTau1      6    1 -1.9542599 0.26333231    -7.421269 1.159073e-13

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             10                     91            -117.55799
       Saturated:             12                     89            -127.19911
    Independence:              6                     95               8.93733
    Number of observations/statistics: 29438/101

    chi-square:  χ² ( df=2 ) = 9.641117,  p = 0.008062282
    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:       -299.558              -97.55799                -97.55052
    BIC:      -1053.952              -14.65758                -46.43743
    CFI: 0.9412838
    TLI: 0.8238514   (also known as NNFI)
    RMSEA:  0.01139224  [95% CI (0.003446383, 0.02034633)]
    Prob(RMSEA <= 0.05): 1
    timestamp: 2020-11-29 10:21:41
    Wall clock time: 0.2394009 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Get the SRMR
    osmasemSRMR(mx.fit0b)

    [1] 0.04402732

    ## Get the coefficients
    coef(mx.fit0b)

           b31        b32        b43        p21
     0.2649761 -0.2964373  0.3500633 -0.2386957

    ## Extract the variance component
    VarCorr(mx.fit0b)

               Tau2_1    Tau2_2      Tau2_3     Tau2_4     Tau2_5     Tau2_6
    Tau2_1 0.01946871 0.0000000 0.000000000 0.00000000 0.00000000 0.00000000
    Tau2_2 0.00000000 0.0237227 0.000000000 0.00000000 0.00000000 0.00000000
    Tau2_3 0.00000000 0.0000000 0.005933922 0.00000000 0.00000000 0.00000000
    Tau2_4 0.00000000 0.0000000 0.000000000 0.01500509 0.00000000 0.00000000
    Tau2_5 0.00000000 0.0000000 0.000000000 0.00000000 0.00856591 0.00000000
    Tau2_6 0.00000000 0.0000000 0.000000000 0.00000000 0.00000000 0.02007019

### Using `SES` as a moderator on the A matrix on the just-identified model

    ## data.SES as a moderator
    A1 <- create.modMatrix(RAM=RAM4, output="A", mod="SES")
    A1

           neg          pos          enga         achiev
    neg    "0"          "0"          "0"          "0"
    pos    "0"          "0"          "0"          "0"
    enga   "0*data.SES" "0*data.SES" "0"          "0"
    achiev "0*data.SES" "0*data.SES" "0*data.SES" "0"

    mx.fit1 <- osmasem(model.name="Ax as moderator", RAM=RAM4, Ax=A1, data= my.df)
    summary(mx.fit1)

    Summary of Ax as moderator

    free parameters:
         name  matrix    row  col     Estimate  Std.Error A      z value
    1     b31      A0   enga  neg  0.247627711 0.04297917     5.76157441
    2     b41      A0 achiev  neg  0.039945846 0.02754367     1.45027309
    3     b32      A0   enga  pos -0.225527642 0.04086165    -5.51929859
    4     b42      A0 achiev  pos -0.098288930 0.03138366    -3.13185076
    5     b43      A0 achiev enga  0.244497002 0.04242094     5.76359242
    6     p21      S0    pos  neg -0.239388682 0.04053847    -5.90522249
    7   b31_1      A1   enga  neg  0.019830422 0.04592343     0.43181489
    8   b41_1      A1 achiev  neg -0.005822283 0.02831309    -0.20563923
    9   b32_1      A1   enga  pos -0.058565588 0.04046453    -1.44733129
    10  b42_1      A1 achiev  pos  0.001098254 0.03183874     0.03449425
    11  b43_1      A1 achiev enga -0.050966282 0.04162396    -1.22444567
    12 Tau1_1 vecTau1      1    1 -1.970678728 0.21993273    -8.96037041
    13 Tau1_2 vecTau1      2    1 -1.889914072 0.18061610   -10.46370740
    14 Tau1_3 vecTau1      3    1 -2.557254935 0.21908249   -11.67256669
    15 Tau1_4 vecTau1      4    1 -2.506901597 0.42195793    -5.94111739
    16 Tau1_5 vecTau1      5    1 -2.511236292 0.26946316    -9.31940481
    17 Tau1_6 vecTau1      6    1 -2.192215522 0.24436292    -8.97114631
           Pr(>|z|)
    1  8.333291e-09
    2  1.469824e-01
    3  3.403554e-08
    4  1.737082e-03
    5  8.234213e-09
    6  3.521711e-09
    7  6.658760e-01
    8  8.370727e-01
    9  1.478042e-01
    10 9.724830e-01
    11 2.207841e-01
    12 0.000000e+00
    13 0.000000e+00
    14 0.000000e+00
    15 2.830858e-09
    16 0.000000e+00
    17 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             17                     84             -131.9935
       Saturated:             27                     74                    NA
    Independence:             12                     89                    NA
    Number of observations/statistics: 29438/101

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -299.9935              -97.99352                -97.97272
    BIC:      -996.3570               42.93718                -11.08858
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:43
    Wall clock time: 1.364671 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Get the coefficients
    coef(mx.fit1)

             b31          b41          b32          b42          b43          p21
     0.247627711  0.039945846 -0.225527642 -0.098288930  0.244497002 -0.239388682
           b31_1        b41_1        b32_1        b42_1        b43_1
     0.019830422 -0.005822283 -0.058565588  0.001098254 -0.050966282

    ## Get the R2
    osmasemR2(mx.fit1, mx.fit0a)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.019310977 0.023387259 0.006004925 0.011797028 0.006680692 0.014430057

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.019421833 0.022826614 0.006008922 0.006645581 0.006588217 0.012469981

    $R2
      Tau2_1_1   Tau2_2_2   Tau2_3_3   Tau2_4_4   Tau2_5_5   Tau2_6_6
    0.00000000 0.02397225 0.00000000 0.43667332 0.01384221 0.13583288

    ## Extract the variance component
    VarCorr(mx.fit1)

               Tau2_1     Tau2_2      Tau2_3      Tau2_4      Tau2_5     Tau2_6
    Tau2_1 0.01942183 0.00000000 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_2 0.00000000 0.02282661 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_3 0.00000000 0.00000000 0.006008922 0.000000000 0.000000000 0.00000000
    Tau2_4 0.00000000 0.00000000 0.000000000 0.006645581 0.000000000 0.00000000
    Tau2_5 0.00000000 0.00000000 0.000000000 0.000000000 0.006588217 0.00000000
    Tau2_6 0.00000000 0.00000000 0.000000000 0.000000000 0.000000000 0.01246998

    ## Comparing with the model without the moderator
    anova(mx.fit1, mx.fit0a)

                 base            comparison ep  minus2LL df       AIC   diffLL
    1 Ax as moderator                  <NA> 17 -131.9935 84 -299.9935       NA
    2 Ax as moderator Just identified model 12 -127.1991 89 -305.1991 4.794415
      diffdf         p
    1     NA        NA
    2      5 0.4414819

### Using `SES` as a moderator on the S matrix

    S1 <- create.modMatrix(RAM=RAM4, output="S", mod="SES")
    S1

           neg          pos          enga achiev
    neg    "0"          "0*data.SES" "0"  "0"
    pos    "0*data.SES" "0"          "0"  "0"
    enga   "0"          "0"          "0"  "0"
    achiev "0"          "0"          "0"  "0"

    mx.fit2 <- osmasem(model.name="Sx as moderator", RAM=RAM4, Sx=S1, data= my.df)
    summary(mx.fit2)

    Summary of Sx as moderator

    free parameters:
         name  matrix    row  col    Estimate  Std.Error A    z value     Pr(>|z|)
    1     b31      A0   enga  neg  0.25341804 0.04086685     6.201066 5.608183e-10
    2     b41      A0 achiev  neg  0.04264173 0.02700073     1.579281 1.142717e-01
    3     b32      A0   enga  pos -0.24065291 0.04559656    -5.277874 1.306913e-07
    4     b42      A0 achiev  pos -0.09173710 0.03147241    -2.914842 3.558691e-03
    5     b43      A0 achiev enga  0.23723933 0.04503766     5.267577 1.382363e-07
    6     p21      S0    pos  neg -0.23700732 0.03856954    -6.144936 7.999579e-10
    7   p21_1      S1    pos  neg -0.06195800 0.04074913    -1.520474 1.283919e-01
    8  Tau1_1 vecTau1      1    1 -2.03486183 0.22008855    -9.245651 0.000000e+00
    9  Tau1_2 vecTau1      2    1 -1.88594578 0.18130177   -10.402247 0.000000e+00
    10 Tau1_3 vecTau1      3    1 -2.54500464 0.22122876   -11.503950 0.000000e+00
    11 Tau1_4 vecTau1      4    1 -2.30296748 0.33650537    -6.843776 7.713163e-12
    12 Tau1_5 vecTau1      5    1 -2.51530485 0.27003026    -9.314900 0.000000e+00
    13 Tau1_6 vecTau1      6    1 -2.11483013 0.23847689    -8.868071 0.000000e+00

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             13                     88             -129.4354
       Saturated:             27                     74                    NA
    Independence:             12                     89                    NA
    Number of observations/statistics: 29438/101

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -305.4354            -103.435442               -103.42307
    BIC:     -1034.9591               4.335099                -36.97872
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:44
    Wall clock time: 0.8297784 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Get the coefficients
    coef(mx.fit2)

            b31         b41         b32         b42         b43         p21
     0.25341804  0.04264173 -0.24065291 -0.09173710  0.23723933 -0.23700732
          p21_1
    -0.06195800

    ## R2
    osmasemR2(mx.fit2, mx.fit0a)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.019310977 0.023387259 0.006004925 0.011797028 0.006680692 0.014430057

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.017082108 0.023008500 0.006157963 0.009992355 0.006534825 0.014557335

    $R2
      Tau2_1_1   Tau2_2_2   Tau2_3_3   Tau2_4_4   Tau2_5_5   Tau2_6_6
    0.11541979 0.01619509 0.00000000 0.15297692 0.02183413 0.00000000

    ## Extract the variance component
    VarCorr(mx.fit2)

               Tau2_1    Tau2_2      Tau2_3      Tau2_4      Tau2_5     Tau2_6
    Tau2_1 0.01708211 0.0000000 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_2 0.00000000 0.0230085 0.000000000 0.000000000 0.000000000 0.00000000
    Tau2_3 0.00000000 0.0000000 0.006157963 0.000000000 0.000000000 0.00000000
    Tau2_4 0.00000000 0.0000000 0.000000000 0.009992355 0.000000000 0.00000000
    Tau2_5 0.00000000 0.0000000 0.000000000 0.000000000 0.006534825 0.00000000
    Tau2_6 0.00000000 0.0000000 0.000000000 0.000000000 0.000000000 0.01455734

    ## Comparing with the model without the moderator
    anova(mx.fit2, mx.fit0a)

                 base            comparison ep  minus2LL df       AIC   diffLL
    1 Sx as moderator                  <NA> 13 -129.4354 88 -305.4354       NA
    2 Sx as moderator Just identified model 12 -127.1991 89 -305.1991 2.236333
      diffdf         p
    1     NA        NA
    2      1 0.1348003

## Scalco17

- Data: Scalco, A., Noventa, S., Sartori, R., & Ceschi, A. (2017). Predicting organic food consumption: A meta-analytic structural equation model based on the theory of planned behavior. Appetite, 112, 235-248.

### Data preparation

    my.df <- Cor2DataFrame(Scalco17$data, Scalco17$n, acov = "weighted")

    ## Add the centered Age as a moderator
    my.df$data <- data.frame(my.df$data, Age=scale(Scalco17$Age,  scale=FALSE), check.names=FALSE)
    head(my.df$data)

                                          SN_ATT PBC_ATT BI_ATT BEH_ATT PBC_SN
    Al-Swidi et al., 2014                  0.562    0.18  0.798      NA  0.314
    Arvola et al., 2008 (Study A from IT)  0.690    0.44  0.730      NA  0.460
    Arvola et al., 2008 (Study A from UK)  0.520    0.22  0.600      NA  0.280
    Arvola et al., 2008 (Study A from FI)  0.570    0.40  0.670      NA  0.340
    Arvola et al., 2008 (Study B from IT)  0.760    0.35  0.710      NA  0.360
    Arvola et al., 2008 (Study B from UK)  0.460    0.03  0.550      NA  0.150
                                          BI_SN BEH_SN BI_PBC BEH_PBC BEH_BI
    Al-Swidi et al., 2014                 0.696     NA  0.216      NA     NA
    Arvola et al., 2008 (Study A from IT) 0.620     NA  0.410      NA     NA
    Arvola et al., 2008 (Study A from UK) 0.560     NA  0.310      NA     NA
    Arvola et al., 2008 (Study A from FI) 0.550     NA  0.360      NA     NA
    Arvola et al., 2008 (Study B from IT) 0.640     NA  0.240      NA     NA
    Arvola et al., 2008 (Study B from UK) 0.580     NA  0.100      NA     NA
                                          C(SN_ATT SN_ATT) C(PBC_ATT SN_ATT)
    Al-Swidi et al., 2014                      0.003943223      0.0007678278
    Arvola et al., 2008 (Study A from IT)      0.003591841      0.0006994029
    Arvola et al., 2008 (Study A from UK)      0.002687242      0.0005232731
    Arvola et al., 2008 (Study A from FI)      0.003627768      0.0007064044
    Arvola et al., 2008 (Study B from IT)      0.003591841      0.0006993999
    Arvola et al., 2008 (Study B from UK)      0.002687225      0.0005232525
                                          C(BI_ATT SN_ATT) C(BEH_ATT SN_ATT)
    Al-Swidi et al., 2014                     0.0011950591      0.0012021339
    Arvola et al., 2008 (Study A from IT)     0.0010885646      0.0010950075
    Arvola et al., 2008 (Study A from UK)     0.0008144157      0.0008192398
    Arvola et al., 2008 (Study A from FI)     0.0010994556      0.0011059662
    Arvola et al., 2008 (Study B from IT)     0.0010885637      0.0010950059
    Arvola et al., 2008 (Study B from UK)     0.0008144047      0.0008192225
                                          C(PBC_SN SN_ATT) C(BI_SN SN_ATT)
    Al-Swidi et al., 2014                     0.0010869714     0.001853515
    Arvola et al., 2008 (Study A from IT)     0.0009901025     0.001688345
    Arvola et al., 2008 (Study A from UK)     0.0007407627     0.001263144
    Arvola et al., 2008 (Study A from FI)     0.0010000130     0.001705237
    Arvola et al., 2008 (Study B from IT)     0.0009901030     0.001688346
    Arvola et al., 2008 (Study B from UK)     0.0007407392     0.001263129
                                          C(BEH_SN SN_ATT) C(BI_PBC SN_ATT)
    Al-Swidi et al., 2014                      0.001499225     0.0005440999
    Arvola et al., 2008 (Study A from IT)      0.001365624     0.0004956093
    Arvola et al., 2008 (Study A from UK)      0.001021705     0.0003708079
    Arvola et al., 2008 (Study A from FI)      0.001379292     0.0005005723
    Arvola et al., 2008 (Study B from IT)      0.001365626     0.0004956081
    Arvola et al., 2008 (Study B from UK)      0.001021684     0.0003707851
                                          C(BEH_PBC SN_ATT) C(BEH_BI SN_ATT)
    Al-Swidi et al., 2014                      0.0004085998     0.0006118210
    Arvola et al., 2008 (Study A from IT)      0.0003721824     0.0005572991
    Arvola et al., 2008 (Study A from UK)      0.0002784627     0.0004169530
    Arvola et al., 2008 (Study A from FI)      0.0003759111     0.0005628791
    Arvola et al., 2008 (Study B from IT)      0.0003721815     0.0005572983
    Arvola et al., 2008 (Study B from UK)      0.0002784457     0.0004169373
                                          C(PBC_ATT PBC_ATT) C(BI_ATT PBC_ATT)
    Al-Swidi et al., 2014                        0.004547980      0.0007764244
    Arvola et al., 2008 (Study A from IT)        0.004142713      0.0007072360
    Arvola et al., 2008 (Study A from UK)        0.003099374      0.0005291247
    Arvola et al., 2008 (Study A from FI)        0.004184142      0.0007143110
    Arvola et al., 2008 (Study B from IT)        0.004142704      0.0007072310
    Arvola et al., 2008 (Study B from UK)        0.003099359      0.0005291132
                                          C(BEH_ATT PBC_ATT) C(PBC_SN PBC_ATT)
    Al-Swidi et al., 2014                       0.0013663034       0.001663768
    Arvola et al., 2008 (Study A from IT)       0.0012445510       0.001515507
    Arvola et al., 2008 (Study A from UK)       0.0009311184       0.001133840
    Arvola et al., 2008 (Study A from FI)       0.0012570002       0.001530667
    Arvola et al., 2008 (Study B from IT)       0.0012445421       0.001515501
    Arvola et al., 2008 (Study B from UK)       0.0009311011       0.001133818
                                          C(BI_SN PBC_ATT) C(BEH_SN PBC_ATT)
    Al-Swidi et al., 2014                     0.0004362854      0.0005789466
    Arvola et al., 2008 (Study A from IT)     0.0003974052      0.0005273539
    Arvola et al., 2008 (Study A from UK)     0.0002973308      0.0003945554
    Arvola et al., 2008 (Study A from FI)     0.0004013863      0.0005326355
    Arvola et al., 2008 (Study B from IT)     0.0003974024      0.0005273502
    Arvola et al., 2008 (Study B from UK)     0.0002973135      0.0003945311
                                          C(BI_PBC PBC_ATT) C(BEH_PBC PBC_ATT)
    Al-Swidi et al., 2014                       0.002592417        0.001688969
    Arvola et al., 2008 (Study A from IT)       0.002361408        0.001538463
    Arvola et al., 2008 (Study A from UK)       0.001766696        0.001151009
    Arvola et al., 2008 (Study A from FI)       0.002385024        0.001553850
    Arvola et al., 2008 (Study B from IT)       0.002361399        0.001538456
    Arvola et al., 2008 (Study B from UK)       0.001766679        0.001150996
                                          C(BEH_BI PBC_ATT) C(BI_ATT BI_ATT)
    Al-Swidi et al., 2014                      0.0006653377      0.002028710
    Arvola et al., 2008 (Study A from IT)      0.0006060497      0.001847932
    Arvola et al., 2008 (Study A from UK)      0.0004534243      0.001382531
    Arvola et al., 2008 (Study A from FI)      0.0006121135      0.001866414
    Arvola et al., 2008 (Study B from IT)      0.0006060431      0.001847930
    Arvola et al., 2008 (Study B from UK)      0.0004534083      0.001382524
                                          C(BEH_ATT BI_ATT) C(PBC_SN BI_ATT)
    Al-Swidi et al., 2014                      0.0012950270     0.0004215199
    Arvola et al., 2008 (Study A from IT)      0.0011796260     0.0003839524
    Arvola et al., 2008 (Study A from UK)      0.0008825405     0.0002872638
    Arvola et al., 2008 (Study A from FI)      0.0011914262     0.0003877963
    Arvola et al., 2008 (Study B from IT)      0.0011796222     0.0003839500
    Arvola et al., 2008 (Study B from UK)      0.0008825309     0.0002872490
                                          C(BI_SN BI_ATT) C(BEH_SN BI_ATT)
    Al-Swidi et al., 2014                    0.0005366953     0.0005252386
    Arvola et al., 2008 (Study A from IT)    0.0004888678     0.0004784312
    Arvola et al., 2008 (Study A from UK)    0.0003657523     0.0003579463
    Arvola et al., 2008 (Study A from FI)    0.0004937613     0.0004832217
    Arvola et al., 2008 (Study B from IT)    0.0004888671     0.0004784302
    Arvola et al., 2008 (Study B from UK)    0.0003657430     0.0003579326
                                          C(BI_PBC BI_ATT) C(BEH_PBC BI_ATT)
    Al-Swidi et al., 2014                     0.0005461550      0.0003593272
    Arvola et al., 2008 (Study A from IT)     0.0004974844      0.0003273035
    Arvola et al., 2008 (Study A from UK)     0.0003722009      0.0002448788
    Arvola et al., 2008 (Study A from FI)     0.0005024612      0.0003305790
    Arvola et al., 2008 (Study B from IT)     0.0004974805      0.0003273006
    Arvola et al., 2008 (Study B from UK)     0.0003721880      0.0002448692
                                          C(BEH_BI BI_ATT) C(BEH_ATT BEH_ATT)
    Al-Swidi et al., 2014                     0.0005717437        0.003394126
    Arvola et al., 2008 (Study A from IT)     0.0005207946        0.003091676
    Arvola et al., 2008 (Study A from UK)     0.0003896368        0.002313040
    Arvola et al., 2008 (Study A from FI)     0.0005260058        0.003122599
    Arvola et al., 2008 (Study B from IT)     0.0005207921        0.003091670
    Arvola et al., 2008 (Study B from UK)     0.0003896277        0.002313025
                                          C(PBC_SN BEH_ATT) C(BI_SN BEH_ATT)
    Al-Swidi et al., 2014                      0.0006189004     0.0006147049
    Arvola et al., 2008 (Study A from IT)      0.0005637424     0.0005599251
    Arvola et al., 2008 (Study A from UK)      0.0004217775     0.0004189166
    Arvola et al., 2008 (Study A from FI)      0.0005693872     0.0005655319
    Arvola et al., 2008 (Study B from IT)      0.0005637377     0.0005599233
    Arvola et al., 2008 (Study B from UK)      0.0004217558     0.0004189021
                                          C(BEH_SN BEH_ATT) C(BI_PBC BEH_ATT)
    Al-Swidi et al., 2014                      0.0010606127      0.0007797350
    Arvola et al., 2008 (Study A from IT)      0.0009660980      0.0007102493
    Arvola et al., 2008 (Study A from UK)      0.0007227971      0.0005313840
    Arvola et al., 2008 (Study A from FI)      0.0009757684      0.0007173553
    Arvola et al., 2008 (Study B from IT)      0.0009660956      0.0007102417
    Arvola et al., 2008 (Study B from UK)      0.0007227765      0.0005313645
                                          C(BEH_PBC BEH_ATT) C(BEH_BI BEH_ATT)
    Al-Swidi et al., 2014                       0.0007167352      0.0014083867
    Arvola et al., 2008 (Study A from IT)       0.0006528621      0.0012828855
    Arvola et al., 2008 (Study A from UK)       0.0004884476      0.0009597956
    Arvola et al., 2008 (Study A from FI)       0.0006593946      0.0012957187
    Arvola et al., 2008 (Study B from IT)       0.0006528562      0.0012828807
    Arvola et al., 2008 (Study B from UK)       0.0004884331      0.0009597820
                                          C(PBC_SN PBC_SN) C(BI_SN PBC_SN)
    Al-Swidi et al., 2014                      0.004844332    0.0010602621
    Arvola et al., 2008 (Study A from IT)      0.004412653    0.0009657759
    Arvola et al., 2008 (Study A from UK)      0.003301335    0.0007225586
    Arvola et al., 2008 (Study A from FI)      0.004456786    0.0009754417
    Arvola et al., 2008 (Study B from IT)      0.004412651    0.0009657763
    Arvola et al., 2008 (Study B from UK)      0.003301314    0.0007225395
                                          C(BEH_SN PBC_SN) C(BI_PBC PBC_SN)
    Al-Swidi et al., 2014                      0.001559438      0.002233172
    Arvola et al., 2008 (Study A from IT)      0.001420471      0.002034172
    Arvola et al., 2008 (Study A from UK)      0.001062741      0.001521879
    Arvola et al., 2008 (Study A from FI)      0.001434685      0.002054518
    Arvola et al., 2008 (Study B from IT)      0.001420471      0.002034167
    Arvola et al., 2008 (Study B from UK)      0.001062716      0.001521857
                                          C(BEH_PBC PBC_SN) C(BEH_BI PBC_SN)
    Al-Swidi et al., 2014                       0.001572052     0.0006097933
    Arvola et al., 2008 (Study A from IT)       0.001431961     0.0005554507
    Arvola et al., 2008 (Study A from UK)       0.001071334     0.0004155717
    Arvola et al., 2008 (Study A from FI)       0.001446286     0.0005610109
    Arvola et al., 2008 (Study B from IT)       0.001431958     0.0005554473
    Arvola et al., 2008 (Study B from UK)       0.001071317     0.0004155527
                                          C(BI_SN BI_SN) C(BEH_SN BI_SN)
    Al-Swidi et al., 2014                    0.002860045     0.001725473
    Arvola et al., 2008 (Study A from IT)    0.002605187     0.001571714
    Arvola et al., 2008 (Study A from UK)    0.001949075     0.001175887
    Arvola et al., 2008 (Study A from FI)    0.002631245     0.001587440
    Arvola et al., 2008 (Study B from IT)    0.002605187     0.001571715
    Arvola et al., 2008 (Study B from UK)    0.001949063     0.001175870
                                          C(BI_PBC BI_SN) C(BEH_PBC BI_SN)
    Al-Swidi et al., 2014                    0.0005120235     0.0003559346
    Arvola et al., 2008 (Study A from IT)    0.0004663928     0.0003242119
    Arvola et al., 2008 (Study A from UK)    0.0003489458     0.0002425705
    Arvola et al., 2008 (Study A from FI)    0.0004710632     0.0003274601
    Arvola et al., 2008 (Study B from IT)    0.0004663917     0.0003242109
    Arvola et al., 2008 (Study B from UK)    0.0003489269     0.0002425566
                                          C(BEH_BI BI_SN) C(BEH_SN BEH_SN)
    Al-Swidi et al., 2014                    0.0006331380      0.003763989
    Arvola et al., 2008 (Study A from IT)    0.0005767177      0.003428580
    Arvola et al., 2008 (Study A from UK)    0.0004314786      0.002565100
    Arvola et al., 2008 (Study A from FI)    0.0005824908      0.003462875
    Arvola et al., 2008 (Study B from IT)    0.0005767166      0.003428580
    Arvola et al., 2008 (Study B from UK)    0.0004314659      0.002565079
                                          C(BI_PBC BEH_SN) C(BEH_PBC BEH_SN)
    Al-Swidi et al., 2014                     0.0006727061      0.0005788168
    Arvola et al., 2008 (Study A from IT)     0.0006127572      0.0005272330
    Arvola et al., 2008 (Study A from UK)     0.0004584516      0.0003944635
    Arvola et al., 2008 (Study A from FI)     0.0006188919      0.0005325122
    Arvola et al., 2008 (Study B from IT)     0.0006127554      0.0005272313
    Arvola et al., 2008 (Study B from UK)     0.0004584256      0.0003944440
                                          C(BEH_BI BEH_SN) C(BI_PBC BI_PBC)
    Al-Swidi et al., 2014                     0.0011957508      0.004234764
    Arvola et al., 2008 (Study A from IT)     0.0010891967      0.003857407
    Arvola et al., 2008 (Study A from UK)     0.0008148909      0.002885925
    Arvola et al., 2008 (Study A from FI)     0.0011000953      0.003895983
    Arvola et al., 2008 (Study B from IT)     0.0010891950      0.003857399
    Arvola et al., 2008 (Study B from UK)     0.0008148736      0.002885907
                                          C(BEH_PBC BI_PBC) C(BEH_BI BI_PBC)
    Al-Swidi et al., 2014                       0.002227324     0.0009451229
    Arvola et al., 2008 (Study A from IT)       0.002028846     0.0008609023
    Arvola et al., 2008 (Study A from UK)       0.001517889     0.0006440932
    Arvola et al., 2008 (Study A from FI)       0.002049137     0.0008695145
    Arvola et al., 2008 (Study B from IT)       0.002028840     0.0008608970
    Arvola et al., 2008 (Study B from UK)       0.001517875     0.0006440759
                                          C(BEH_PBC BEH_PBC) C(BEH_BI BEH_PBC)
    Al-Swidi et al., 2014                        0.003778506      0.0006339266
    Arvola et al., 2008 (Study A from IT)        0.003441803      0.0005774350
    Arvola et al., 2008 (Study A from UK)        0.002574988      0.0004320159
    Arvola et al., 2008 (Study A from FI)        0.003476224      0.0005832128
    Arvola et al., 2008 (Study B from IT)        0.003441799      0.0005774309
    Arvola et al., 2008 (Study B from UK)        0.002574977      0.0004320028
                                          C(BEH_BI BEH_BI)       Age
    Al-Swidi et al., 2014                      0.002176137 -2.398421
    Arvola et al., 2008 (Study A from IT)      0.001982224  2.871579
    Arvola et al., 2008 (Study A from UK)      0.001483003  2.871579
    Arvola et al., 2008 (Study A from FI)      0.002002049  2.871579
    Arvola et al., 2008 (Study B from IT)      0.001982220  2.871579
    Arvola et al., 2008 (Study B from UK)      0.001482992  2.871579

    ## Display the pairwise no. of studies
    pattern.na(my.df, show.na=FALSE, type="osmasem")

            SN_ATT PBC_ATT BI_ATT BEH_ATT PBC_SN BI_SN BEH_SN BI_PBC BEH_PBC BEH_BI
    SN_ATT      22      22     23      22     22    23     22     23      22     22
    PBC_ATT     22      22     23      22     22    23     22     23      22     22
    BI_ATT      23      23     23      23     23    23     23     23      23     23
    BEH_ATT     22      22     23       6     22    23      6     23       6      6
    PBC_SN      22      22     23      22     22    23     22     23      22     22
    BI_SN       23      23     23      23     23    23     23     23      23     23
    BEH_SN      22      22     23       6     22    23      6     23       6      6
    BI_PBC      23      23     23      23     23    23     23     23      23     23
    BEH_PBC     22      22     23       6     22    23      6     23       6      6
    BEH_BI      22      22     23       6     22    23      6     23       6      6

    ## Proposed model
    model6 <- 'BI ~ b41*ATT + b42*SN + b43*PBC
               BEH ~ b54*BI + b53*PBC
               ATT ~~ p21*SN
               ATT ~~ p31*PBC
               SN ~~ p32*PBC
               ATT ~~ 1*ATT
               SN ~~ 1*SN
               PBC ~~ 1*PBC
               BI ~~ p44*BI
               BEH ~~ p44*BEH'

    plot(model6, color="yellow")

    RAM6 <- lavaan2RAM(model6, obs.variables=c("ATT", "SN", "PBC", "BI","BEH"))
    RAM6

    $A
        ATT     SN      PBC     BI      BEH
    ATT "0"     "0"     "0"     "0"     "0"
    SN  "0"     "0"     "0"     "0"     "0"
    PBC "0"     "0"     "0"     "0"     "0"
    BI  "0*b41" "0*b42" "0*b43" "0"     "0"
    BEH "0"     "0"     "0*b53" "0*b54" "0"

    $S
        ATT     SN      PBC     BI      BEH
    ATT "1"     "0*p21" "0*p31" "0"     "0"
    SN  "0*p21" "1"     "0*p32" "0"     "0"
    PBC "0*p31" "0*p32" "1"     "0"     "0"
    BI  "0"     "0"     "0"     "0*p44" "0"
    BEH "0"     "0"     "0"     "0"     "0*p44"

    $F
        ATT SN PBC BI BEH
    ATT   1  0   0  0   0
    SN    0  1   0  0   0
    PBC   0  0   1  0   0
    BI    0  0   0  1   0
    BEH   0  0   0  0   1

    $M
      ATT SN PBC BI BEH
    1   0  0   0  0   0

## Fitting an over-identified model

    mx.fit0 <- osmasem(model.name="No moderator", RAM=RAM6, data=my.df)

    ## Get the chi-square statistic on the proposed model
    summary(mx.fit0, fitIndices=TRUE)

    Summary of No moderator

    free parameters:
          name  matrix row col   Estimate  Std.Error A    z value     Pr(>|z|)
    1      b41      A0  BI ATT  0.4548365 0.03958523    11.490056 0.000000e+00
    2      b42      A0  BI  SN  0.2781751 0.04707429     5.909279 3.436090e-09
    3      b43      A0  BI PBC  0.1369545 0.03334843     4.106776 4.012195e-05
    4      b53      A0 BEH PBC  0.1305297 0.07190047     1.815423 6.945898e-02
    5      b54      A0 BEH  BI  0.5825793 0.07659635     7.605836 2.819966e-14
    6      p21      S0  SN ATT  0.4127349 0.04003478    10.309409 0.000000e+00
    7      p31      S0 PBC ATT  0.2549528 0.03031158     8.411071 0.000000e+00
    8      p32      S0 PBC  SN  0.2367971 0.02308961    10.255571 0.000000e+00
    9   Tau1_1 vecTau1   1   1 -1.7087420 0.16320416   -10.469966 0.000000e+00
    10  Tau1_2 vecTau1   2   1 -2.0237771 0.17270815   -11.717901 0.000000e+00
    11  Tau1_3 vecTau1   3   1 -1.9815611 0.15824811   -12.521863 0.000000e+00
    12  Tau1_4 vecTau1   4   1 -2.1785702 0.32607797    -6.681133 2.371037e-11
    13  Tau1_5 vecTau1   5   1 -2.3625783 0.20417531   -11.571322 0.000000e+00
    14  Tau1_6 vecTau1   6   1 -1.8009593 0.15775994   -11.415822 0.000000e+00
    15  Tau1_7 vecTau1   7   1 -2.0852314 0.32604091    -6.395613 1.599048e-10
    16  Tau1_8 vecTau1   8   1 -2.1423355 0.17392978   -12.317244 0.000000e+00
    17  Tau1_9 vecTau1   9   1 -1.9887403 0.29545685    -6.731068 1.684208e-11
    18 Tau1_10 vecTau1  10   1 -1.7197174 0.35066166    -4.904207 9.380549e-07

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             18                    141             -153.2011
       Saturated:             20                    139             -156.9808
    Independence:             10                    149              159.8486
    Number of observations/statistics: 11349/159

    chi-square:  χ² ( df=2 ) = 3.779663,  p = 0.1510972
    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -435.2011             -117.20110               -117.14073
    BIC:     -1469.7019               14.86283                -42.33897
    CFI: 0.9941998
    TLI: 0.9709991   (also known as NNFI)
    RMSEA:  0.008854726  [95% CI (0, 0.02466774)]
    Prob(RMSEA <= 0.05): 1
    timestamp: 2020-11-29 10:21:45
    Wall clock time: 0.5672746 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    ## Get the SRMR
    osmasemSRMR(mx.fit0)

    [1] 0.03592283

    ## Extract the variance component
    diag(VarCorr(mx.fit0))

         Tau2_1      Tau2_2      Tau2_3      Tau2_4      Tau2_5      Tau2_6
    0.032794846 0.017465040 0.019003689 0.012814982 0.008869325 0.027271351
         Tau2_7      Tau2_8      Tau2_9     Tau2_10
    0.015445112 0.013778153 0.018732777 0.032082815

### Using `Age` as a moderator on the A matrix

    ## An index for data without missing value in Age
    indexAge <- !is.na(Scalco17$Age)
    indexAge

     [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE
    [13]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE

    A1 <- create.modMatrix(RAM=RAM6, output="A", mod="Age")
    A1

        ATT          SN           PBC          BI           BEH
    ATT "0"          "0"          "0"          "0"          "0"
    SN  "0"          "0"          "0"          "0"          "0"
    PBC "0"          "0"          "0"          "0"          "0"
    BI  "0*data.Age" "0*data.Age" "0*data.Age" "0"          "0"
    BEH "0"          "0"          "0*data.Age" "0*data.Age" "0"

    ## Refit the model without any moderator with indexAge
    mx.fit0.Age <- osmasem(model.name="No moderator", RAM=RAM6, data=my.df, subset.rows=indexAge)

    mx.fit1 <- osmasem(model.name="Ax as moderator", RAM=RAM6, Ax=A1, data=my.df, subset.rows=indexAge)
    summary(mx.fit1)

    Summary of Ax as moderator

    free parameters:
          name  matrix row col     Estimate   Std.Error A     z value     Pr(>|z|)
    1      b41      A0  BI ATT  0.494316458 0.037016479    13.3539569 0.000000e+00
    2      b42      A0  BI  SN  0.270283661 0.045625149     5.9240061 3.141919e-09
    3      b43      A0  BI PBC  0.106333949 0.035297714     3.0124883 2.591154e-03
    4      b53      A0 BEH PBC  0.231819524 0.214809189     1.0791881 2.805039e-01
    5      b54      A0 BEH  BI  0.505361776 0.128009169     3.9478561 7.885418e-05
    6      p21      S0  SN ATT  0.436652318 0.045342884     9.6300076 0.000000e+00
    7      p31      S0 PBC ATT  0.254545016 0.033163651     7.6754219 1.643130e-14
    8      p32      S0 PBC  SN  0.245076491 0.027056722     9.0578781 0.000000e+00
    9    b41_1      A1  BI ATT -0.004274833 0.005017516    -0.8519821 3.942240e-01
    10   b42_1      A1  BI  SN  0.018113821 0.005343634     3.3897948 6.994495e-04
    11   b43_1      A1  BI PBC  0.001718695 0.004055060     0.4238397 6.716827e-01
    12   b53_1      A1 BEH PBC -0.013887873 0.037502815    -0.3703155 7.111474e-01
    13   b54_1      A1 BEH  BI  0.018597143 0.021048819     0.8835243 3.769530e-01
    14  Tau1_1 vecTau1   1   1 -1.683369084 0.178878361    -9.4106916 0.000000e+00
    15  Tau1_2 vecTau1   2   1 -2.034762249 0.191728002   -10.6127547 0.000000e+00
    16  Tau1_3 vecTau1   3   1 -2.168585817 0.175989934   -12.3222151 0.000000e+00
    17  Tau1_4 vecTau1   4   1 -2.526114259 0.414673842    -6.0918100 1.116411e-09
    18  Tau1_5 vecTau1   5   1 -2.292340087 0.218273573   -10.5021422 0.000000e+00
    19  Tau1_6 vecTau1   6   1 -2.043288463 0.180149895   -11.3421574 0.000000e+00
    20  Tau1_7 vecTau1   7   1 -2.518005904 0.464483984    -5.4210823 5.923928e-08
    21  Tau1_8 vecTau1   8   1 -2.220179181 0.195574302   -11.3521008 0.000000e+00
    22  Tau1_9 vecTau1   9   1 -1.869188022 0.365488980    -5.1142117 3.150538e-07
    23 Tau1_10 vecTau1  10   1 -2.653198679 0.479619579    -5.5318815 3.168139e-08

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             23                    104             -149.3285
       Saturated:             65                     62                    NA
    Independence:             20                    107                    NA
    Number of observations/statistics: 9507/127

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -357.3285             -103.32854               -103.21212
    BIC:     -1101.9460               61.34649                -11.74391
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:49
    Wall clock time: 1.607328 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    coef(mx.fit1)

             b41          b42          b43          b53          b54          p21
     0.494316458  0.270283661  0.106333949  0.231819524  0.505361776  0.436652318
             p31          p32        b41_1        b42_1        b43_1        b53_1
     0.254545016  0.245076491 -0.004274833  0.018113821  0.001718695 -0.013887873
           b54_1
     0.018597143

    osmasemR2(mx.fit1, mx.fit0.Age)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.034362514 0.016891805 0.013896527 0.006192528 0.010061669 0.032521082
       Tau2_7_7    Tau2_8_8    Tau2_9_9  Tau2_10_10
    0.007863517 0.013393079 0.023999247 0.009817024

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.034501994 0.017085511 0.013073452 0.006395066 0.010207014 0.016796631
       Tau2_7_7    Tau2_8_8    Tau2_9_9  Tau2_10_10
    0.006499618 0.011791712 0.023792710 0.004959763

    $R2
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.000000000 0.000000000 0.059228826 0.000000000 0.000000000 0.483515598
       Tau2_7_7    Tau2_8_8    Tau2_9_9  Tau2_10_10
    0.173446376 0.119566728 0.008605956 0.494779411

    ## Extract the variance component
    diag(VarCorr(mx.fit1))

         Tau2_1      Tau2_2      Tau2_3      Tau2_4      Tau2_5      Tau2_6
    0.034501994 0.017085511 0.013073452 0.006395066 0.010207014 0.016796631
         Tau2_7      Tau2_8      Tau2_9     Tau2_10
    0.006499618 0.011791712 0.023792710 0.004959763

    anova(mx.fit1, mx.fit0.Age)

                 base   comparison ep  minus2LL  df       AIC   diffLL diffdf
    1 Ax as moderator         <NA> 23 -149.3285 104 -357.3285       NA     NA
    2 Ax as moderator No moderator 18 -131.9769 109 -349.9769 17.35165      5
                p
    1          NA
    2 0.003878977

### Using `Age` as a moderator on the S matrix

    S1 <- create.modMatrix(RAM=RAM6, output="S", mod="Age")
    S1

        ATT          SN           PBC          BI  BEH
    ATT "0"          "0*data.Age" "0*data.Age" "0" "0"
    SN  "0*data.Age" "0"          "0*data.Age" "0" "0"
    PBC "0*data.Age" "0*data.Age" "0"          "0" "0"
    BI  "0"          "0"          "0"          "0" "0"
    BEH "0"          "0"          "0"          "0" "0"

    mx.fit2 <- osmasem(model.name="Sx as moderator", RAM=RAM6, Sx=S1, data=my.df, subset.rows=indexAge)
    summary(mx.fit2)

    Summary of Sx as moderator

    free parameters:
          name  matrix row col     Estimate   Std.Error A    z value     Pr(>|z|)
    1      b41      A0  BI ATT  0.502034463 0.036526112    13.744536 0.000000e+00
    2      b42      A0  BI  SN  0.271247591 0.044979799     6.030431 1.635226e-09
    3      b43      A0  BI PBC  0.112593416 0.034449120     3.268397 1.081584e-03
    4      b53      A0 BEH PBC  0.159775431 0.094407255     1.692406 9.056854e-02
    5      b54      A0 BEH  BI  0.608563313 0.058046451    10.484074 0.000000e+00
    6      p21      S0  SN ATT  0.419470520 0.038367664    10.932918 0.000000e+00
    7      p31      S0 PBC ATT  0.242344119 0.029628741     8.179359 2.220446e-16
    8      p32      S0 PBC  SN  0.236407032 0.024730115     9.559480 0.000000e+00
    9    p21_1      S1  SN ATT  0.020822499 0.004735239     4.397349 1.095812e-05
    10   p31_1      S1 PBC ATT  0.009453332 0.003791840     2.493072 1.266430e-02
    11   p32_1      S1 PBC  SN  0.008224392 0.003490346     2.356326 1.845674e-02
    12  Tau1_1 vecTau1   1   1 -1.869313976 0.185085301   -10.099743 0.000000e+00
    13  Tau1_2 vecTau1   2   1 -2.179173619 0.197283295   -11.045910 0.000000e+00
    14  Tau1_3 vecTau1   3   1 -2.146493293 0.176182666   -12.183340 0.000000e+00
    15  Tau1_4 vecTau1   4   1 -2.594880011 0.407063687    -6.374629 1.834064e-10
    16  Tau1_5 vecTau1   5   1 -2.422205209 0.224718485   -10.778843 0.000000e+00
    17  Tau1_6 vecTau1   6   1 -1.997519337 0.183052377   -10.912283 0.000000e+00
    18  Tau1_7 vecTau1   7   1 -2.515028564 0.467182848    -5.383392 7.309497e-08
    19  Tau1_8 vecTau1   8   1 -2.198277019 0.194076569   -11.326854 0.000000e+00
    20  Tau1_9 vecTau1   9   1 -1.871291067 0.366273650    -5.108997 3.238731e-07
    21 Tau1_10 vecTau1  10   1 -2.411961800 0.418193854    -5.767569 8.042329e-09

    Model Statistics:
                   |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
           Model:             21                    106             -161.4281
       Saturated:             65                     62                    NA
    Independence:             20                    107                    NA
    Number of observations/statistics: 9507/127

    Information Criteria:
          |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
    AIC:      -373.4281             -119.42814                -119.3307
    BIC:     -1132.3652               30.92731                 -35.8074
    To get additional fit indices, see help(mxRefModels)
    timestamp: 2020-11-29 10:21:50
    Wall clock time: 1.225012 secs
    optimizer:  SLSQP
    OpenMx version number: 2.18.1
    Need help?  See help(mxSummary)

    coef(mx.fit2)

            b41         b42         b43         b53         b54         p21
    0.502034463 0.271247591 0.112593416 0.159775431 0.608563313 0.419470520
            p31         p32       p21_1       p31_1       p32_1
    0.242344119 0.236407032 0.020822499 0.009453332 0.008224392

    osmasemR2(mx.fit2, mx.fit0.Age)

    $Tau2.0
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.034362514 0.016891805 0.013896527 0.006192528 0.010061669 0.032521082
       Tau2_7_7    Tau2_8_8    Tau2_9_9  Tau2_10_10
    0.007863517 0.013393079 0.023999247 0.009817024

    $Tau2.1
       Tau2_1_1    Tau2_2_2    Tau2_3_3    Tau2_4_4    Tau2_5_5    Tau2_6_6
    0.023786717 0.012799525 0.013664055 0.005573344 0.007872257 0.018406735
       Tau2_7_7    Tau2_8_8    Tau2_9_9  Tau2_10_10
    0.006538437 0.012319720 0.023692846 0.008035198

    $R2
      Tau2_1_1   Tau2_2_2   Tau2_3_3   Tau2_4_4   Tau2_5_5   Tau2_6_6   Tau2_7_7
    0.30777134 0.24226424 0.01672879 0.09998889 0.21759923 0.43400608 0.16850983
      Tau2_8_8   Tau2_9_9 Tau2_10_10
    0.08014278 0.01276709 0.18150368

    ## Extract the variance component
    diag(VarCorr(mx.fit2))

         Tau2_1      Tau2_2      Tau2_3      Tau2_4      Tau2_5      Tau2_6
    0.023786717 0.012799525 0.013664055 0.005573344 0.007872257 0.018406735
         Tau2_7      Tau2_8      Tau2_9     Tau2_10
    0.006538437 0.012319720 0.023692846 0.008035198

    anova(mx.fit2, mx.fit0.Age)

                 base   comparison ep  minus2LL  df       AIC   diffLL diffdf
    1 Sx as moderator         <NA> 21 -161.4281 106 -373.4281       NA     NA
    2 Sx as moderator No moderator 18 -131.9769 109 -349.9769 29.45126      3
                 p
    1           NA
    2 1.800082e-06

# Installation and help

- First of all, you need R to run it. You can install the `metaSEM` package by running the following command inside an `R` session:

    install.packages("metaSEM")

- The developmental version can be installed from GitHub by running:

    library(devtools)
    install_github("mikewlcheung/metasem")

## Potential issues with the OpenMx available at CRAN

- The OpenMx available at CRAN includes the open source `SLSQP` and `CSOLNP` optimizers.
- If the `SLSQP` optimizer (default in the `metaSEM` package) does not work well for you, e.g., there are many error codes, you may try to rerun the analysis. For example,

    random1 <- meta(y=yi, v=vi, data=Hox02)
    random1 <- rerun(random1)
    summary(random1)

- If you still prefer the non-free `NPSOL` optimizer, you may install it from the OpenMx website and call it by running:

    mxOption(NULL, "Default optimizer", "NPSOL")

## Help

- If you need help, please send your questions to OpenMx discussion forum, a discussion forum for the `metaSEM` package in `OpenMx`. Please include information on the R session. It will be helpful if you can include a reproducible example. You may save a copy of your data, say `my.df`, and attach the content of `myData.R` in the post by using

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

The files are based on the following versions of `R` and `R` packages:

    sessionInfo()

    R version 4.0.3 (2020-10-10)
    Platform: x86_64-pc-linux-gnu (64-bit)
    Running under: Ubuntu 20.04.1 LTS

    Matrix products: default
    BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
    LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

    locale:
     [1] LC_CTYPE=en_SG.UTF-8       LC_NUMERIC=C
     [3] LC_TIME=en_SG.UTF-8        LC_COLLATE=en_SG.UTF-8
     [5] LC_MONETARY=en_SG.UTF-8    LC_MESSAGES=en_SG.UTF-8
     [7] LC_PAPER=en_SG.UTF-8       LC_NAME=C
     [9] LC_ADDRESS=C               LC_TELEPHONE=C
    [11] LC_MEASUREMENT=en_SG.UTF-8 LC_IDENTIFICATION=C

    attached base packages:
    [1] stats     graphics  grDevices utils     datasets  methods   base

    other attached packages:
    [1] semPlot_1.1.2 metafor_2.4-0 Matrix_1.2-18 metaSEM_1.2.5 OpenMx_2.18.1
    [6] knitr_1.30    rmarkdown_2.5

    loaded via a namespace (and not attached):
     [1] nlme_3.1-149        RColorBrewer_1.1-2  mi_1.0
     [4] tools_4.0.3         backports_1.1.10    R6_2.5.0
     [7] d3Network_0.5.2.1   rpart_4.1-15        Hmisc_4.4-1
    [10] colorspace_1.4-1    nnet_7.3-14         tidyselect_1.1.0
    [13] gridExtra_2.3       mnormt_2.0.2        compiler_4.0.3
    [16] fdrtool_1.2.15      qgraph_1.6.5        htmlTable_2.1.0
    [19] regsem_1.6.2        scales_1.1.1        checkmate_2.0.0
    [22] psych_2.0.9         mvtnorm_1.1-1       pbapply_1.4-3
    [25] sem_3.1-11          stringr_1.4.0       digest_0.6.27
    [28] pbivnorm_0.6.0      foreign_0.8-80      minqa_1.2.4
    [31] base64enc_0.1-3     jpeg_0.1-8.1        pkgconfig_2.0.3
    [34] htmltools_0.5.0     lme4_1.1-23         lisrelToR_0.1.4
    [37] htmlwidgets_1.5.2   rlang_0.4.8         huge_1.3.4.1
    [40] rstudioapi_0.11     generics_0.0.2      gtools_3.8.2
    [43] dplyr_1.0.2         zip_2.1.1           magrittr_1.5
    [46] Formula_1.2-3       Rcpp_1.0.5          munsell_0.5.0
    [49] abind_1.4-5         rockchalk_1.8.144   lifecycle_0.2.0
    [52] whisker_0.4         stringi_1.5.3       yaml_2.2.1
    [55] carData_3.0-4       MASS_7.3-53         plyr_1.8.6
    [58] matrixcalc_1.0-3    lavaan_0.6-7        grid_4.0.3
    [61] parallel_4.0.3      crayon_1.3.4        lattice_0.20-41
    [64] kutils_1.70         splines_4.0.3       tmvnsim_1.0-2
    [67] pillar_1.4.6        igraph_1.2.6        rjson_0.2.20
    [70] boot_1.3-25         corpcor_1.6.9       BDgraph_2.63
    [73] reshape2_1.4.4      stats4_4.0.3        XML_3.99-0.5
    [76] glue_1.4.2          evaluate_0.14       latticeExtra_0.6-29
    [79] data.table_1.13.2   png_0.1-7           vctrs_0.3.4
    [82] nloptr_1.2.2.2      gtable_0.3.0        purrr_0.3.4
    [85] ggplot2_3.3.2       xfun_0.19           openxlsx_4.2.2
    [88] xtable_1.8-4        coda_0.19-4         Rsolnp_1.16
    [91] glasso_1.11         survival_3.2-7      truncnorm_1.0-8
    [94] tibble_3.0.4        arm_1.11-2          ellipse_0.4.2
    [97] cluster_2.1.0       statmod_1.4.34      ellipsis_0.3.1