Type
as a covariateYear
and Year^2
as covariatesDiscipline
as a covariateCountry
as a covariateType
and Discipline
as covariatesType
and Country
as covariatesDiscipline
and Country
as covariatesType
, Discipline
, and Country
as covariatesThe 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:
## 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 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 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.)
## 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.)
## 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 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.)
## 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.)
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.
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.
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
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.
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.)
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.)
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.)
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.)
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\)
## 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
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.## 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
LB
in the intervals.type
argument.## 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.)
## 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.)
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
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.)
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.
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
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 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 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
Type
as a covariateGrants
) 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
.
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.)
Grants
and Fellowship
as indicator variables
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.)
Year
and Year^2
as covariatescbind()
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.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 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
Discipline
as a covariateDiscipline
. multidisciplinary
is used as the reference group in the analysis.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.)
Discipline
is significant
## 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
Country
as a covariateCountry
. The United States
is used as the reference group in the analysis.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.)
Discipline
is 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
Type
and Discipline
as 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.
Discipline
is significant after controlling for Type
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
Type
and Country
as 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.)
Country
is significant after controlling for Type
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
Discipline
and Country
as covariates## 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.)
Type
, Discipline
, and Country
as 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.)
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.
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.)
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.)
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.)
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.)
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
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
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
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)
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"
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.)
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
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")
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.)
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
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")
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 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
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 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
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.)
## 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 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
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.)
## Plot the model with labels
plot(random2, whatLabels="path", color="red")
## Plot the parameter estimates
plot(random2, color="green")
## 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
# 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-06-13 12:31:02
Wall clock time: 0.247977 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
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
Lag
as a moderator on the A matrixA1 <- matrix(c(0,0,0,0,
0,0,0,0,
"0*data.Lag","0*data.Lag",0,0,
"0*data.Lag","0*data.Lag",0,0),
nrow=4, ncol=4, byrow=TRUE)
A1
[,1] [,2] [,3] [,4]
[1,] "0" "0" "0" "0"
[2,] "0" "0" "0" "0"
[3,] "0*data.Lag" "0*data.Lag" "0" "0"
[4,] "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-06-13 12:31:03
Wall clock time: 0.3948386 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
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
Lag
as a moderator on the S matrixS1 <- matrix(c(0,"0*data.Lag",0,0,
"0*data.Lag",0,0,0,
0,0,0,"0*data.Lag",
0,0,"0*data.Lag",0),
nrow=4, ncol=4, byrow=TRUE)
S1
[,1] [,2] [,3] [,4]
[1,] "0" "0*data.Lag" "0" "0"
[2,] "0*data.Lag" "0" "0" "0"
[3,] "0" "0" "0" "0*data.Lag"
[4,] "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-06-13 12:31:03
Wall clock time: 0.3480837 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
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
Lag
as a moderator on both the A and S matricesfit3 <- 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-06-13 12:31:04
Wall clock time: 0.4624577 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
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
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
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-06-13 12:31:05
Wall clock time: 0.2675793 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
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
## 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="spring", 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-06-13 12:31:05
Wall clock time: 0.2313678 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
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
SES
as a moderator on the A matrix## data.SES as a moderator
A1 <- matrix(c(0,0,0,0,
0,0,0,0,
"0*data.SES","0*data.SES",0,0,
"0*data.SES","0*data.SES","0*data.SES",0),
nrow=4, ncol=4, byrow=TRUE)
A1
[,1] [,2] [,3] [,4]
[1,] "0" "0" "0" "0"
[2,] "0" "0" "0" "0"
[3,] "0*data.SES" "0*data.SES" "0" "0"
[4,] "0*data.SES" "0*data.SES" "0*data.SES" "0"
mx.fit1 <- osmasem(model.name="Ax as moderator", RAM=RAM5, Ax=A1, data= my.df)
summary(mx.fit1)
Summary of Ax as moderator
free parameters:
name matrix row col Estimate Std.Error A
1 b31 A0 enga neg 0.264294952 0.03923351
2 b32 A0 enga pos -0.270675676 0.05551137
3 b43 A0 achiev enga 0.350683664 0.04092856
4 p21 S0 pos neg -0.238994221 0.04060491
5 b31_1 A1 enga neg 0.010983038 0.04750488
6 Ax as moderator.A1[4,1] A1 achiev neg -0.007902378 0.03456298
7 b32_1 A1 enga pos -0.040062980 0.05194494
8 Ax as moderator.A1[4,2] A1 achiev pos 0.007729473 0.04484276
9 b43_1 A1 achiev enga -0.043108942 0.05449184
10 Tau1_1 vecTau1 1 1 -1.968731780 0.21989251
11 Tau1_2 vecTau1 2 1 -1.881526012 0.18354753
12 Tau1_3 vecTau1 3 1 -2.566353837 0.22558656
13 Tau1_4 vecTau1 4 1 -2.314847939 0.45477990
14 Tau1_5 vecTau1 5 1 -2.324124004 0.27482075
15 Tau1_6 vecTau1 6 1 -2.032378061 0.27033768
z value Pr(>|z|)
1 6.7364587 1.622924e-11
2 -4.8760399 1.082367e-06
3 8.5681904 0.000000e+00
4 -5.8858458 3.960235e-09
5 0.2311981 8.171609e-01
6 -0.2286370 8.191510e-01
7 -0.7712586 4.405536e-01
8 0.1723684 8.631479e-01
9 -0.7911082 4.288809e-01
10 -8.9531553 0.000000e+00
11 -10.2508924 0.000000e+00
12 -11.3763597 0.000000e+00
13 -5.0900401 3.579878e-07
14 -8.4568723 0.000000e+00
15 -7.5179239 5.573320e-14
Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 15 86 -120.8623
Saturated: 27 74 NA
Independence: 12 89 NA
Number of observations/statistics: 29438/101
Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: -292.8623 -90.86226 -90.84595
BIC: -1005.8058 33.48836 -14.18143
To get additional fit indices, see help(mxRefModels)
timestamp: 2020-06-13 12:31:07
Wall clock time: 1.270855 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
Need help? See help(mxSummary)
## Get the coefficients
coef(mx.fit1)
b31 b32 b43
0.264294952 -0.270675676 0.350683664
p21 b31_1 Ax as moderator.A1[4,1]
-0.238994221 0.010983038 -0.007902378
b32_1 Ax as moderator.A1[4,2] b43_1
-0.040062980 0.007729473 -0.043108942
## 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.019497607 0.023212786 0.005900562 0.009757726 0.009578368 0.017167175
$R2
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.000000000 0.007460173 0.017379509 0.172865745 0.000000000 0.000000000
## Extract the variance component
VarCorr(mx.fit1)
Tau2_1 Tau2_2 Tau2_3 Tau2_4 Tau2_5 Tau2_6
Tau2_1 0.01949761 0.00000000 0.000000000 0.000000000 0.000000000 0.00000000
Tau2_2 0.00000000 0.02321279 0.000000000 0.000000000 0.000000000 0.00000000
Tau2_3 0.00000000 0.00000000 0.005900562 0.000000000 0.000000000 0.00000000
Tau2_4 0.00000000 0.00000000 0.000000000 0.009757726 0.000000000 0.00000000
Tau2_5 0.00000000 0.00000000 0.000000000 0.000000000 0.009578368 0.00000000
Tau2_6 0.00000000 0.00000000 0.000000000 0.000000000 0.000000000 0.01716718
## Comparing with the model without the moderator
anova(mx.fit1, mx.fit0a)
base comparison ep minus2LL df AIC diffLL
1 Ax as moderator <NA> 15 -120.8623 86 -292.8623 NA
2 Ax as moderator Just identified model 12 -127.1991 89 -305.1991 -6.336844
diffdf p
1 NA NA
2 3 1
SES
as a moderator on the S matrixS1 <- matrix(c(0,"0*data.SES",0,0,
"0*data.SES",0,0,0,
0,0,0,0,
0,0,0,0),
nrow=4, ncol=4, byrow=TRUE)
S1
[,1] [,2] [,3] [,4]
[1,] "0" "0*data.SES" "0" "0"
[2,] "0*data.SES" "0" "0" "0"
[3,] "0" "0" "0" "0"
[4,] "0" "0" "0" "0"
mx.fit2 <- osmasem(model.name="Sx as moderator", RAM=RAM5, 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.26460021 0.03800017 6.963132 3.328005e-12
2 b32 A0 enga pos -0.28569742 0.05296008 -5.394580 6.868394e-08
3 b43 A0 achiev enga 0.35103259 0.04143698 8.471482 0.000000e+00
4 p21 S0 pos neg -0.23709129 0.03862286 -6.138627 8.323799e-10
5 p21_1 S1 pos neg -0.06092496 0.04150602 -1.467858 1.421427e-01
6 Tau1_1 vecTau1 1 1 -2.03327701 0.22020772 -9.233450 0.000000e+00
7 Tau1_2 vecTau1 2 1 -1.87774706 0.18384252 -10.213889 0.000000e+00
8 Tau1_3 vecTau1 3 1 -2.55448455 0.22766473 -11.220379 0.000000e+00
9 Tau1_4 vecTau1 4 1 -2.19263227 0.37117135 -5.907332 3.476933e-09
10 Tau1_5 vecTau1 5 1 -2.36446301 0.27705445 -8.534290 0.000000e+00
11 Tau1_6 vecTau1 6 1 -1.95193319 0.26533855 -7.356387 1.889600e-13
Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 11 90 -119.6397
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.6397 -97.639745 -97.63077
BIC: -1045.7435 -6.449287 -41.40713
To get additional fit indices, see help(mxRefModels)
timestamp: 2020-06-13 12:31:08
Wall clock time: 0.7399418 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
Need help? See help(mxSummary)
## Get the coefficients
coef(mx.fit2)
b31 b32 b43 p21 p21_1
0.26460021 -0.28569742 0.35103259 -0.23709129 -0.06092496
## 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.017136338 0.023388891 0.006042309 0.012459592 0.008835956 0.020163800
$R2
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.1126115 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## Extract the variance component
VarCorr(mx.fit2)
Tau2_1 Tau2_2 Tau2_3 Tau2_4 Tau2_5 Tau2_6
Tau2_1 0.01713634 0.00000000 0.000000000 0.00000000 0.000000000 0.0000000
Tau2_2 0.00000000 0.02338889 0.000000000 0.00000000 0.000000000 0.0000000
Tau2_3 0.00000000 0.00000000 0.006042309 0.00000000 0.000000000 0.0000000
Tau2_4 0.00000000 0.00000000 0.000000000 0.01245959 0.000000000 0.0000000
Tau2_5 0.00000000 0.00000000 0.000000000 0.00000000 0.008835956 0.0000000
Tau2_6 0.00000000 0.00000000 0.000000000 0.00000000 0.000000000 0.0201638
## Comparing with the model without the moderator
anova(mx.fit2, mx.fit0a)
Warning in collectStatistics1(otherStats, ref, other, bootPair): Model 'Sx as
moderator' has more degrees of freedom than Just identified model which means
that the models need to be compared in the opposite order
base comparison ep minus2LL df AIC diffLL
1 Sx as moderator <NA> 11 -119.6397 90 -299.6397 NA
2 Sx as moderator Just identified model 12 -127.1991 89 -305.1991 -7.559364
diffdf p
1 NA NA
2 -1 NA
# select studies without missing value on Age
data <- Scalco17$data[is.na(Scalco17$Age)==FALSE]
n <- Scalco17$n[is.na(Scalco17$Age)==FALSE]
Age <- Scalco17$Age[is.na(Scalco17$Age)==FALSE]
my.df <- Cor2DataFrame(data, n, acov = "weighted")
## Add the centered Age as a moderator
my.df$data <- data.frame(my.df$data, Age=scale(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.003894946 0.0007859557
Arvola et al., 2008 (Study A from IT) 0.003547861 0.0007159400
Arvola et al., 2008 (Study A from UK) 0.002654327 0.0005356288
Arvola et al., 2008 (Study A from FI) 0.003583361 0.0007231242
Arvola et al., 2008 (Study B from IT) 0.003547860 0.0007159322
Arvola et al., 2008 (Study B from UK) 0.002654321 0.0005356212
C(BI_ATT SN_ATT) C(BEH_ATT SN_ATT)
Al-Swidi et al., 2014 0.0011346410 0.0012594666
Arvola et al., 2008 (Study A from IT) 0.0010335433 0.0011472401
Arvola et al., 2008 (Study A from UK) 0.0007732444 0.0008583078
Arvola et al., 2008 (Study A from FI) 0.0010438925 0.0011587346
Arvola et al., 2008 (Study B from IT) 0.0010335432 0.0011472392
Arvola et al., 2008 (Study B from UK) 0.0007732402 0.0008583019
C(PBC_SN SN_ATT) C(BI_SN SN_ATT)
Al-Swidi et al., 2014 0.0010444207 0.001906300
Arvola et al., 2008 (Study A from IT) 0.0009513602 0.001736415
Arvola et al., 2008 (Study A from UK) 0.0007117583 0.001299097
Arvola et al., 2008 (Study A from FI) 0.0009608962 0.001753797
Arvola et al., 2008 (Study B from IT) 0.0009513543 0.001736414
Arvola et al., 2008 (Study B from UK) 0.0007117501 0.001299092
C(BEH_SN SN_ATT) C(BI_PBC SN_ATT)
Al-Swidi et al., 2014 0.001564648 0.0005640044
Arvola et al., 2008 (Study A from IT) 0.001425214 0.0005137615
Arvola et al., 2008 (Study A from UK) 0.001066274 0.0003843701
Arvola et al., 2008 (Study A from FI) 0.001439488 0.0005189262
Arvola et al., 2008 (Study B from IT) 0.001425215 0.0005137553
Arvola et al., 2008 (Study B from UK) 0.001066268 0.0003843636
C(BEH_PBC SN_ATT) C(BEH_BI SN_ATT)
Al-Swidi et al., 2014 0.0003771588 0.0005704365
Arvola et al., 2008 (Study A from IT) 0.0003435592 0.0005196029
Arvola et al., 2008 (Study A from UK) 0.0002570336 0.0003887420
Arvola et al., 2008 (Study A from FI) 0.0003470153 0.0005248166
Arvola et al., 2008 (Study B from IT) 0.0003435572 0.0005196024
Arvola et al., 2008 (Study B from UK) 0.0002570289 0.0003887383
C(PBC_ATT PBC_ATT) C(BI_ATT PBC_ATT)
Al-Swidi et al., 2014 0.004582124 0.0006865668
Arvola et al., 2008 (Study A from IT) 0.004173813 0.0006254033
Arvola et al., 2008 (Study A from UK) 0.003122630 0.0004678942
Arvola et al., 2008 (Study A from FI) 0.004215567 0.0006316700
Arvola et al., 2008 (Study B from IT) 0.004173799 0.0006253985
Arvola et al., 2008 (Study B from UK) 0.003122618 0.0004678877
C(BEH_ATT PBC_ATT) C(PBC_SN PBC_ATT)
Al-Swidi et al., 2014 0.001476461 0.001698940
Arvola et al., 2008 (Study A from IT) 0.001344913 0.001547552
Arvola et al., 2008 (Study A from UK) 0.001006196 0.001157797
Arvola et al., 2008 (Study A from FI) 0.001358382 0.001563044
Arvola et al., 2008 (Study B from IT) 0.001344905 0.001547537
Arvola et al., 2008 (Study B from UK) 0.001006187 0.001157783
C(BI_SN PBC_ATT) C(BEH_SN PBC_ATT)
Al-Swidi et al., 2014 0.0004422697 0.0005909667
Arvola et al., 2008 (Study A from IT) 0.0004028708 0.0005383181
Arvola et al., 2008 (Study A from UK) 0.0003014079 0.0004027429
Arvola et al., 2008 (Study A from FI) 0.0004069198 0.0005437235
Arvola et al., 2008 (Study B from IT) 0.0004028639 0.0005383102
Arvola et al., 2008 (Study B from UK) 0.0003013994 0.0004027328
C(BI_PBC PBC_ATT) C(BEH_PBC PBC_ATT)
Al-Swidi et al., 2014 0.002780122 0.001778309
Arvola et al., 2008 (Study A from IT) 0.002532377 0.001619837
Arvola et al., 2008 (Study A from UK) 0.001894593 0.001211878
Arvola et al., 2008 (Study A from FI) 0.002557719 0.001636049
Arvola et al., 2008 (Study B from IT) 0.002532363 0.001619828
Arvola et al., 2008 (Study B from UK) 0.001894580 0.001211869
C(BEH_BI PBC_ATT) C(BI_ATT BI_ATT)
Al-Swidi et al., 2014 0.0006609423 0.001774417
Arvola et al., 2008 (Study A from IT) 0.0006020510 0.001616325
Arvola et al., 2008 (Study A from UK) 0.0004504248 0.001209251
Arvola et al., 2008 (Study A from FI) 0.0006080873 0.001632496
Arvola et al., 2008 (Study B from IT) 0.0006020450 0.001616324
Arvola et al., 2008 (Study B from UK) 0.0004504182 0.001209248
C(BEH_ATT BI_ATT) C(PBC_SN BI_ATT)
Al-Swidi et al., 2014 0.0013079417 0.0003808364
Arvola et al., 2008 (Study A from IT) 0.0011914075 0.0003469104
Arvola et al., 2008 (Study A from UK) 0.0008913508 0.0002595404
Arvola et al., 2008 (Study A from FI) 0.0012033336 0.0003503919
Arvola et al., 2008 (Study B from IT) 0.0011914069 0.0003469069
Arvola et al., 2008 (Study B from UK) 0.0008913462 0.0002595334
C(BI_SN BI_ATT) C(BEH_SN BI_ATT)
Al-Swidi et al., 2014 0.0004761852 0.0004869782
Arvola et al., 2008 (Study A from IT) 0.0004337567 0.0004435876
Arvola et al., 2008 (Study A from UK) 0.0003245155 0.0003318708
Arvola et al., 2008 (Study A from FI) 0.0004381059 0.0004480367
Arvola et al., 2008 (Study B from IT) 0.0004337568 0.0004435880
Arvola et al., 2008 (Study B from UK) 0.0003245113 0.0003318659
C(BI_PBC BI_ATT) C(BEH_PBC BI_ATT)
Al-Swidi et al., 2014 0.0004943840 0.0003008892
Arvola et al., 2008 (Study A from IT) 0.0004503437 0.0002740850
Arvola et al., 2008 (Study A from UK) 0.0003369241 0.0002050563
Arvola et al., 2008 (Study A from FI) 0.0004548612 0.0002768366
Arvola et al., 2008 (Study B from IT) 0.0004503397 0.0002740834
Arvola et al., 2008 (Study B from UK) 0.0003369180 0.0002050520
C(BEH_BI BI_ATT) C(BEH_ATT BEH_ATT)
Al-Swidi et al., 2014 0.0004709275 0.003023503
Arvola et al., 2008 (Study A from IT) 0.0004289659 0.002754089
Arvola et al., 2008 (Study A from UK) 0.0003209313 0.002060469
Arvola et al., 2008 (Study A from FI) 0.0004332652 0.002781648
Arvola et al., 2008 (Study B from IT) 0.0004289655 0.002754087
Arvola et al., 2008 (Study B from UK) 0.0003209282 0.002060462
C(PBC_SN BEH_ATT) C(BI_SN BEH_ATT)
Al-Swidi et al., 2014 0.0006150517 0.0006149783
Arvola et al., 2008 (Study A from IT) 0.0005602563 0.0005601786
Arvola et al., 2008 (Study A from UK) 0.0004191561 0.0004190992
Arvola et al., 2008 (Study A from FI) 0.0005658778 0.0005657984
Arvola et al., 2008 (Study B from IT) 0.0005602495 0.0005601776
Arvola et al., 2008 (Study B from UK) 0.0004191458 0.0004190929
C(BEH_SN BEH_ATT) C(BI_PBC BEH_ATT)
Al-Swidi et al., 2014 0.0009062449 0.0008205288
Arvola et al., 2008 (Study A from IT) 0.0008254916 0.0007474243
Arvola et al., 2008 (Study A from UK) 0.0006175932 0.0005591861
Arvola et al., 2008 (Study A from FI) 0.0008337666 0.0007549201
Arvola et al., 2008 (Study B from IT) 0.0008254905 0.0007474170
Arvola et al., 2008 (Study B from UK) 0.0006175857 0.0005591773
C(BEH_PBC BEH_ATT) C(BEH_BI BEH_ATT)
Al-Swidi et al., 2014 0.0005476718 0.0011377485
Arvola et al., 2008 (Study A from IT) 0.0004988769 0.0010363578
Arvola et al., 2008 (Study A from UK) 0.0003732353 0.0007753517
Arvola et al., 2008 (Study A from FI) 0.0005038824 0.0010467359
Arvola et al., 2008 (Study B from IT) 0.0004988738 0.0010363563
Arvola et al., 2008 (Study B from UK) 0.0003732291 0.0007753470
C(PBC_SN PBC_SN) C(BI_SN PBC_SN)
Al-Swidi et al., 2014 0.004822479 0.0009802739
Arvola et al., 2008 (Study A from IT) 0.004392744 0.0008929287
Arvola et al., 2008 (Study A from UK) 0.003286423 0.0006680441
Arvola et al., 2008 (Study A from FI) 0.004436686 0.0009018752
Arvola et al., 2008 (Study B from IT) 0.004392733 0.0008929233
Arvola et al., 2008 (Study B from UK) 0.003286411 0.0006680358
C(BEH_SN PBC_SN) C(BI_PBC PBC_SN)
Al-Swidi et al., 2014 0.001673767 0.002320233
Arvola et al., 2008 (Study A from IT) 0.001524628 0.002113483
Arvola et al., 2008 (Study A from UK) 0.001140650 0.001581197
Arvola et al., 2008 (Study A from FI) 0.001539893 0.002134634
Arvola et al., 2008 (Study B from IT) 0.001524622 0.002113469
Arvola et al., 2008 (Study B from UK) 0.001140640 0.001581185
C(BEH_PBC PBC_SN) C(BEH_BI PBC_SN)
Al-Swidi et al., 2014 0.001653831 0.0005719372
Arvola et al., 2008 (Study A from IT) 0.001506456 0.0005209782
Arvola et al., 2008 (Study A from UK) 0.001127051 0.0003897703
Arvola et al., 2008 (Study A from FI) 0.001521532 0.0005262024
Arvola et al., 2008 (Study B from IT) 0.001506447 0.0005209731
Arvola et al., 2008 (Study B from UK) 0.001127042 0.0003897635
C(BI_SN BI_SN) C(BEH_SN BI_SN)
Al-Swidi et al., 2014 0.002738121 0.001836723
Arvola et al., 2008 (Study A from IT) 0.002494127 0.001673052
Arvola et al., 2008 (Study A from UK) 0.001865978 0.001251693
Arvola et al., 2008 (Study A from FI) 0.002519082 0.001689799
Arvola et al., 2008 (Study B from IT) 0.002494126 0.001673051
Arvola et al., 2008 (Study B from UK) 0.001865973 0.001251687
C(BI_PBC BI_SN) C(BEH_PBC BI_SN)
Al-Swidi et al., 2014 0.0005233625 0.0003211938
Arvola et al., 2008 (Study A from IT) 0.0004767413 0.0002925805
Arvola et al., 2008 (Study A from UK) 0.0003566746 0.0002188949
Arvola et al., 2008 (Study A from FI) 0.0004815300 0.0002955224
Arvola et al., 2008 (Study B from IT) 0.0004767357 0.0002925785
Arvola et al., 2008 (Study B from UK) 0.0003566672 0.0002188895
C(BEH_BI BI_SN) C(BEH_SN BEH_SN)
Al-Swidi et al., 2014 0.0005686601 0.003414714
Arvola et al., 2008 (Study A from IT) 0.0005179854 0.003110426
Arvola et al., 2008 (Study A from UK) 0.0003875319 0.002327063
Arvola et al., 2008 (Study A from FI) 0.0005231789 0.003141549
Arvola et al., 2008 (Study B from IT) 0.0005179846 0.003110425
Arvola et al., 2008 (Study B from UK) 0.0003875279 0.002327055
C(BI_PBC BEH_SN) C(BEH_PBC BEH_SN)
Al-Swidi et al., 2014 0.0007000802 0.0004650752
Arvola et al., 2008 (Study A from IT) 0.0006377103 0.0004236437
Arvola et al., 2008 (Study A from UK) 0.0004771038 0.0003169500
Arvola et al., 2008 (Study A from FI) 0.0006441104 0.0004278976
Arvola et al., 2008 (Study B from IT) 0.0006377039 0.0004236413
Arvola et al., 2008 (Study B from UK) 0.0004770948 0.0003169437
C(BEH_BI BEH_SN) C(BI_PBC BI_PBC)
Al-Swidi et al., 2014 0.0009280306 0.004310893
Arvola et al., 2008 (Study A from IT) 0.0008453324 0.003926758
Arvola et al., 2008 (Study A from UK) 0.0006324364 0.002937796
Arvola et al., 2008 (Study A from FI) 0.0008538005 0.003966043
Arvola et al., 2008 (Study B from IT) 0.0008453313 0.003926744
Arvola et al., 2008 (Study B from UK) 0.0006324315 0.002937784
C(BEH_PBC BI_PBC) C(BEH_BI BI_PBC)
Al-Swidi et al., 2014 0.002432731 0.0009285226
Arvola et al., 2008 (Study A from IT) 0.002215950 0.0008457926
Arvola et al., 2008 (Study A from UK) 0.001657859 0.0006327798
Arvola et al., 2008 (Study A from FI) 0.002238123 0.0008542667
Arvola et al., 2008 (Study B from IT) 0.002215941 0.0008457872
Arvola et al., 2008 (Study B from UK) 0.001657850 0.0006327739
C(BEH_PBC BEH_PBC) C(BEH_BI BEH_PBC)
Al-Swidi et al., 2014 0.003390958 0.0004010992
Arvola et al., 2008 (Study A from IT) 0.003088779 0.0003653611
Arvola et al., 2008 (Study A from UK) 0.002310864 0.0002733456
Arvola et al., 2008 (Study A from FI) 0.003119675 0.0003690273
Arvola et al., 2008 (Study B from IT) 0.003088773 0.0003653585
Arvola et al., 2008 (Study B from UK) 0.002310857 0.0002733414
C(BEH_BI BEH_BI) Age
Al-Swidi et al., 2014 0.001556503 -2.398421
Arvola et al., 2008 (Study A from IT) 0.001417802 2.871579
Arvola et al., 2008 (Study A from UK) 0.001060728 2.871579
Arvola et al., 2008 (Study A from FI) 0.001431990 2.871579
Arvola et al., 2008 (Study B from IT) 0.001417801 2.871579
Arvola et al., 2008 (Study B from UK) 0.001060724 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 18 18 19 18 18 19 18 19 18 18
PBC_ATT 18 18 19 18 18 19 18 19 18 18
BI_ATT 19 19 19 19 19 19 19 19 19 19
BEH_ATT 18 18 19 4 18 19 4 19 4 4
PBC_SN 18 18 19 18 18 19 18 19 18 18
BI_SN 19 19 19 19 19 19 19 19 19 19
BEH_SN 18 18 19 4 18 19 4 19 4 4
BI_PBC 19 19 19 19 19 19 19 19 19 19
BEH_PBC 18 18 19 4 18 19 4 19 4 4
BEH_BI 18 18 19 4 18 19 4 19 4 4
## 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
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.4835277 0.04170564 11.593822 0.000000e+00
2 b42 A0 BI SN 0.2962537 0.06020032 4.921132 8.604519e-07
3 b43 A0 BI PBC 0.1027158 0.03693804 2.780760 5.423186e-03
4 b53 A0 BEH PBC 0.1740509 0.09367856 1.857959 6.317484e-02
5 b54 A0 BEH BI 0.6262584 0.06625250 9.452601 0.000000e+00
6 p21 S0 SN ATT 0.4364109 0.04526608 9.641014 0.000000e+00
7 p31 S0 PBC ATT 0.2534395 0.03287804 7.708474 1.265654e-14
8 p32 S0 PBC SN 0.2452304 0.02698188 9.088708 0.000000e+00
9 Tau1_1 vecTau1 1 1 -1.6849473 0.17884383 -9.421334 0.000000e+00
10 Tau1_2 vecTau1 2 1 -2.0449409 0.19199594 -10.650959 0.000000e+00
11 Tau1_3 vecTau1 3 1 -2.1316440 0.17348442 -12.287236 0.000000e+00
12 Tau1_4 vecTau1 4 1 -2.5231920 0.40426060 -6.241499 4.333980e-10
13 Tau1_5 vecTau1 5 1 -2.2968639 0.21839542 -10.516997 0.000000e+00
14 Tau1_6 vecTau1 6 1 -1.7114049 0.17600126 -9.723822 0.000000e+00
15 Tau1_7 vecTau1 7 1 -2.4037343 0.49631419 -4.843171 1.277835e-06
16 Tau1_8 vecTau1 8 1 -2.1618818 0.19340894 -11.177776 0.000000e+00
17 Tau1_9 vecTau1 9 1 -1.8628790 0.36552979 -5.096381 3.462083e-07
18 Tau1_10 vecTau1 10 1 -2.3005518 0.45855526 -5.016957 5.249645e-07
Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 18 109 -131.9660
Saturated: 20 107 -136.5611
Independence: 10 117 135.0936
Number of observations/statistics: 9507/127
chi-square: χ² ( df=2 ) = 4.595126, p = 0.1005035
Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: -349.966 -95.96601 -95.89392
BIC: -1130.382 32.91010 -24.29109
CFI: 0.9900819
TLI: 0.9504093 (also known as NNFI)
RMSEA: 0.01168268 [95% CI (0, 0.02850145)]
Prob(RMSEA <= 0.05): 0.9999995
timestamp: 2020-06-13 12:31:09
Wall clock time: 0.534013 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
Need help? See help(mxSummary)
## Get the SRMR
osmasemSRMR(mx.fit0)
[1] 0.02889252
## Extract the variance component
VarCorr(mx.fit0)
Tau2_1 Tau2_2 Tau2_3 Tau2_4 Tau2_5 Tau2_6
Tau2_1 0.03439326 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_2 0.00000000 0.01674121 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_3 0.00000000 0.00000000 0.01407594 0.000000000 0.00000000 0.00000000
Tau2_4 0.00000000 0.00000000 0.00000000 0.006432551 0.00000000 0.00000000
Tau2_5 0.00000000 0.00000000 0.00000000 0.000000000 0.01011508 0.00000000
Tau2_6 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.03262065
Tau2_7 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_8 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_9 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_10 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_7 Tau2_8 Tau2_9 Tau2_10
Tau2_1 0.000000000 0.00000000 0.00000000 0.00000000
Tau2_2 0.000000000 0.00000000 0.00000000 0.00000000
Tau2_3 0.000000000 0.00000000 0.00000000 0.00000000
Tau2_4 0.000000000 0.00000000 0.00000000 0.00000000
Tau2_5 0.000000000 0.00000000 0.00000000 0.00000000
Tau2_6 0.000000000 0.00000000 0.00000000 0.00000000
Tau2_7 0.008168511 0.00000000 0.00000000 0.00000000
Tau2_8 0.000000000 0.01324992 0.00000000 0.00000000
Tau2_9 0.000000000 0.00000000 0.02409483 0.00000000
Tau2_10 0.000000000 0.00000000 0.00000000 0.01004075
Age
as a moderator on the A matrixA1 <- matrix(c(0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
"0*data.Age","0*data.Age","0*data.Age",0,0,
0,0,"0*data.Age","0*data.Age",0),
nrow=5, ncol=5, byrow=TRUE)
A1
[,1] [,2] [,3] [,4] [,5]
[1,] "0" "0" "0" "0" "0"
[2,] "0" "0" "0" "0" "0"
[3,] "0" "0" "0" "0" "0"
[4,] "0*data.Age" "0*data.Age" "0*data.Age" "0" "0"
[5,] "0" "0" "0*data.Age" "0*data.Age" "0"
mx.fit1 <- osmasem(model.name="Ax as moderator", RAM=RAM6, 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 Pr(>|z|)
1 b41 A0 BI ATT 0.494226472 0.037079425 13.3288603 0.000000e+00
2 b42 A0 BI SN 0.269821508 0.045638256 5.9121783 3.376127e-09
3 b43 A0 BI PBC 0.106890256 0.035049552 3.0496896 2.290780e-03
4 b53 A0 BEH PBC 0.226760082 0.214463585 1.0573361 2.903582e-01
5 b54 A0 BEH BI 0.507277404 0.127673722 3.9732327 7.090373e-05
6 p21 S0 SN ATT 0.436647775 0.045332695 9.6320719 0.000000e+00
7 p31 S0 PBC ATT 0.254398113 0.033049675 7.6974466 1.398881e-14
8 p32 S0 PBC SN 0.244896521 0.027105282 9.0350109 0.000000e+00
9 b41_1 A1 BI ATT -0.004220379 0.005032280 -0.8386614 4.016594e-01
10 b42_1 A1 BI SN 0.018116306 0.005359144 3.3804476 7.236788e-04
11 b43_1 A1 BI PBC 0.001645495 0.004044294 0.4068684 6.841047e-01
12 b53_1 A1 BEH PBC -0.013210850 0.037494829 -0.3523379 7.245848e-01
13 b54_1 A1 BEH BI 0.018242923 0.021048281 0.8667180 3.860965e-01
14 Tau1_1 vecTau1 1 1 -1.683167410 0.178727308 -9.4175167 0.000000e+00
15 Tau1_2 vecTau1 2 1 -2.039207751 0.191967996 -10.6226444 0.000000e+00
16 Tau1_3 vecTau1 3 1 -2.162536461 0.174254181 -12.4102414 0.000000e+00
17 Tau1_4 vecTau1 4 1 -2.507961605 0.407829970 -6.1495275 7.771412e-10
18 Tau1_5 vecTau1 5 1 -2.289117700 0.217645949 -10.5176215 0.000000e+00
19 Tau1_6 vecTau1 6 1 -2.040032123 0.179314150 -11.3768608 0.000000e+00
20 Tau1_7 vecTau1 7 1 -2.499063496 0.453627141 -5.5090696 3.607353e-08
21 Tau1_8 vecTau1 8 1 -2.224783139 0.196211319 -11.3387095 0.000000e+00
22 Tau1_9 vecTau1 9 1 -1.866924527 0.364289844 -5.1248328 2.978020e-07
23 Tau1_10 vecTau1 10 1 -2.623309470 0.461552632 -5.6836627 1.318402e-08
Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 23 104 -149.2628
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.2628 -103.26284 -103.14642
BIC: -1101.8803 61.41218 -11.67822
To get additional fit indices, see help(mxRefModels)
timestamp: 2020-06-13 12:31:12
Wall clock time: 1.593615 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
Need help? See help(mxSummary)
coef(mx.fit1)
b41 b42 b43 b53 b54 p21
0.494226472 0.269821508 0.106890256 0.226760082 0.507277404 0.436647775
p31 p32 b41_1 b42_1 b43_1 b53_1
0.254398113 0.244896521 -0.004220379 0.018116306 0.001645495 -0.013210850
b54_1
0.018242923
osmasemR2(mx.fit1, mx.fit0)
$Tau2.0
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.034393260 0.016741211 0.014075944 0.006432551 0.010115081 0.032620650
Tau2_7_7 Tau2_8_8 Tau2_9_9 Tau2_10_10
0.008168511 0.013249922 0.024094829 0.010040748
$Tau2.1
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.034515913 0.016934277 0.013232585 0.006631507 0.010273008 0.016906379
Tau2_7_7 Tau2_8_8 Tau2_9_9 Tau2_10_10
0.006750579 0.011683633 0.023900664 0.005265291
$R2
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.000000000 0.000000000 0.059914930 0.000000000 0.000000000 0.481727695
Tau2_7_7 Tau2_8_8 Tau2_9_9 Tau2_10_10
0.173585147 0.118211137 0.008058369 0.475607754
## Extract the variance component
VarCorr(mx.fit1)
Tau2_1 Tau2_2 Tau2_3 Tau2_4 Tau2_5 Tau2_6
Tau2_1 0.03451591 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_2 0.00000000 0.01693428 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_3 0.00000000 0.00000000 0.01323259 0.000000000 0.00000000 0.00000000
Tau2_4 0.00000000 0.00000000 0.00000000 0.006631507 0.00000000 0.00000000
Tau2_5 0.00000000 0.00000000 0.00000000 0.000000000 0.01027301 0.00000000
Tau2_6 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.01690638
Tau2_7 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_8 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_9 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_10 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
Tau2_7 Tau2_8 Tau2_9 Tau2_10
Tau2_1 0.000000000 0.00000000 0.00000000 0.000000000
Tau2_2 0.000000000 0.00000000 0.00000000 0.000000000
Tau2_3 0.000000000 0.00000000 0.00000000 0.000000000
Tau2_4 0.000000000 0.00000000 0.00000000 0.000000000
Tau2_5 0.000000000 0.00000000 0.00000000 0.000000000
Tau2_6 0.000000000 0.00000000 0.00000000 0.000000000
Tau2_7 0.006750579 0.00000000 0.00000000 0.000000000
Tau2_8 0.000000000 0.01168363 0.00000000 0.000000000
Tau2_9 0.000000000 0.00000000 0.02390066 0.000000000
Tau2_10 0.000000000 0.00000000 0.00000000 0.005265291
anova(mx.fit1, mx.fit0)
base comparison ep minus2LL df AIC diffLL diffdf
1 Ax as moderator <NA> 23 -149.2628 104 -357.2628 NA NA
2 Ax as moderator No moderator 18 -131.9660 109 -349.9660 17.29683 5
p
1 NA
2 0.003969926
Age
as a moderator on the S matrixS1 <- matrix(c(0,"0*data.Age","0*data.Age",0,0,
"0*data.Age",0,"0*data.Age",0,0,
"0*data.Age","0*data.Age",0,0,0,
0,0,0,0,0,
0,0,0,0,0),
nrow=5, ncol=5, byrow=TRUE)
mx.fit2 <- osmasem(model.name="Sx as moderator", RAM=RAM6, 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 b41 A0 BI ATT 0.502147856 0.036565019 13.733012 0.000000e+00
2 b42 A0 BI SN 0.270704090 0.044927076 6.025411 1.686806e-09
3 b43 A0 BI PBC 0.112908755 0.034262421 3.295411 9.827786e-04
4 b53 A0 BEH PBC 0.158356966 0.094302439 1.679246 9.310417e-02
5 b54 A0 BEH BI 0.608081928 0.058175044 10.452625 0.000000e+00
6 p21 S0 SN ATT 0.419449215 0.038352779 10.936606 0.000000e+00
7 p31 S0 PBC ATT 0.242380404 0.029557721 8.200240 2.220446e-16
8 p32 S0 PBC SN 0.236297207 0.024784303 9.534148 0.000000e+00
9 p21_1 S1 SN ATT 0.020810720 0.004737072 4.393161 1.117141e-05
10 p31_1 S1 PBC ATT 0.009424310 0.003788768 2.487434 1.286685e-02
11 p32_1 S1 PBC SN 0.008231116 0.003499737 2.351924 1.867659e-02
12 Tau1_1 vecTau1 1 1 -1.869040919 0.184899192 -10.108432 0.000000e+00
13 Tau1_2 vecTau1 2 1 -2.182834138 0.197486293 -11.053092 0.000000e+00
14 Tau1_3 vecTau1 3 1 -2.141432241 0.174540422 -12.268976 0.000000e+00
15 Tau1_4 vecTau1 4 1 -2.573067779 0.399609024 -6.438963 1.202924e-10
16 Tau1_5 vecTau1 5 1 -2.418158664 0.224005080 -10.795106 0.000000e+00
17 Tau1_6 vecTau1 6 1 -1.994278343 0.182161413 -10.947864 0.000000e+00
18 Tau1_7 vecTau1 7 1 -2.494802544 0.455839641 -5.472983 4.425229e-08
19 Tau1_8 vecTau1 8 1 -2.201416698 0.194577226 -11.313846 0.000000e+00
20 Tau1_9 vecTau1 9 1 -1.869557637 0.365239431 -5.118718 3.076196e-07
21 Tau1_10 vecTau1 10 1 -2.395655838 0.408378803 -5.866259 4.457373e-09
Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 21 106 -161.3619
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.3619 -119.36187 -119.26445
BIC: -1132.2989 30.99359 -35.74112
To get additional fit indices, see help(mxRefModels)
timestamp: 2020-06-13 12:31:13
Wall clock time: 1.285182 secs
optimizer: SLSQP
OpenMx version number: 2.17.3
Need help? See help(mxSummary)
coef(mx.fit2)
b41 b42 b43 b53 b54 p21
0.502147856 0.270704090 0.112908755 0.158356966 0.608081928 0.419449215
p31 p32 p21_1 p31_1 p32_1
0.242380404 0.236297207 0.020810720 0.009424310 0.008231116
## The model may be misspecified as the R2 of Tau2_7_7 is -0.86.
osmasemR2(mx.fit2, mx.fit0)
$Tau2.0
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.034393260 0.016741211 0.014075944 0.006432551 0.010115081 0.032620650
Tau2_7_7 Tau2_8_8 Tau2_9_9 Tau2_10_10
0.008168511 0.013249922 0.024094829 0.010040748
$Tau2.1
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6
0.023799711 0.012706161 0.013803067 0.005821860 0.007936227 0.018526434
Tau2_7_7 Tau2_8_8 Tau2_9_9 Tau2_10_10
0.006808353 0.012242603 0.023775128 0.008301561
$R2
Tau2_1_1 Tau2_2_2 Tau2_3_3 Tau2_4_4 Tau2_5_5 Tau2_6_6 Tau2_7_7
0.30801234 0.24102498 0.01938609 0.09493769 0.21540648 0.43206422 0.16651243
Tau2_8_8 Tau2_9_9 Tau2_10_10
0.07602453 0.01326842 0.17321290
## Extract the variance component
VarCorr(mx.fit2)
Tau2_1 Tau2_2 Tau2_3 Tau2_4 Tau2_5 Tau2_6
Tau2_1 0.02379971 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000
Tau2_2 0.00000000 0.01270616 0.00000000 0.00000000 0.000000000 0.00000000
Tau2_3 0.00000000 0.00000000 0.01380307 0.00000000 0.000000000 0.00000000
Tau2_4 0.00000000 0.00000000 0.00000000 0.00582186 0.000000000 0.00000000
Tau2_5 0.00000000 0.00000000 0.00000000 0.00000000 0.007936227 0.00000000
Tau2_6 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.01852643
Tau2_7 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000
Tau2_8 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000
Tau2_9 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000
Tau2_10 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000
Tau2_7 Tau2_8 Tau2_9 Tau2_10
Tau2_1 0.000000000 0.0000000 0.00000000 0.000000000
Tau2_2 0.000000000 0.0000000 0.00000000 0.000000000
Tau2_3 0.000000000 0.0000000 0.00000000 0.000000000
Tau2_4 0.000000000 0.0000000 0.00000000 0.000000000
Tau2_5 0.000000000 0.0000000 0.00000000 0.000000000
Tau2_6 0.000000000 0.0000000 0.00000000 0.000000000
Tau2_7 0.006808353 0.0000000 0.00000000 0.000000000
Tau2_8 0.000000000 0.0122426 0.00000000 0.000000000
Tau2_9 0.000000000 0.0000000 0.02377513 0.000000000
Tau2_10 0.000000000 0.0000000 0.00000000 0.008301561
anova(mx.fit2, mx.fit0)
base comparison ep minus2LL df AIC diffLL diffdf
1 Sx as moderator <NA> 21 -161.3619 106 -373.3619 NA NA
2 Sx as moderator No moderator 18 -131.9660 109 -349.9660 29.39586 3
p
1 NA
2 1.849006e-06
metaSEM
package by running the following command inside an R
session:install.packages("metaSEM")
library(devtools)
install_github("mikewlcheung/metasem")
SLSQP
and CSOLNP
optimizers.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)
NPSOL
optimizer, you may install it from the OpenMx website and call it by running:mxOption(NULL, "Default optimizer", "NPSOL")
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 usingsessionInfo()
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.1 (2020-06-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04 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.4 OpenMx_2.17.3
[6] knitr_1.28 rmarkdown_2.2
loaded via a namespace (and not attached):
[1] nlme_3.1-148 RColorBrewer_1.1-2 mi_1.0
[4] tools_4.0.1 backports_1.1.6 R6_2.4.1
[7] d3Network_0.5.2.1 rpart_4.1-15 Hmisc_4.4-0
[10] colorspace_1.4-1 nnet_7.3-14 tidyselect_1.1.0
[13] gridExtra_2.3 mnormt_1.5-7 compiler_4.0.1
[16] fdrtool_1.2.15 qgraph_1.6.5 htmlTable_1.13.3
[19] regsem_1.5.2 scales_1.1.0 checkmate_2.0.0
[22] psych_1.9.12.31 mvtnorm_1.1-0 pbapply_1.4-2
[25] sem_3.1-11 stringr_1.4.0 digest_0.6.25
[28] pbivnorm_0.6.0 foreign_0.8-79 minqa_1.2.4
[31] base64enc_0.1-3 jpeg_0.1-8.1 pkgconfig_2.0.3
[34] htmltools_0.4.0 lme4_1.1-23 lisrelToR_0.1.4
[37] htmlwidgets_1.5.1 rlang_0.4.6 huge_1.3.4.1
[40] rstudioapi_0.11 generics_0.0.2 gtools_3.8.2
[43] acepack_1.4.1 dplyr_1.0.0 zip_2.0.4
[46] magrittr_1.5 Formula_1.2-3 Rcpp_1.0.4.6
[49] munsell_0.5.0 abind_1.4-5 rockchalk_1.8.144
[52] lifecycle_0.2.0 whisker_0.4 stringi_1.4.6
[55] yaml_2.2.1 carData_3.0-3 MASS_7.3-51.6
[58] plyr_1.8.6 matrixcalc_1.0-3 lavaan_0.6-6
[61] grid_4.0.1 parallel_4.0.1 crayon_1.3.4
[64] lattice_0.20-41 kutils_1.70 splines_4.0.1
[67] pillar_1.4.4 igraph_1.2.5 rjson_0.2.20
[70] boot_1.3-25 corpcor_1.6.9 BDgraph_2.62
[73] reshape2_1.4.4 stats4_4.0.1 XML_3.99-0.3
[76] glue_1.4.1 evaluate_0.14 latticeExtra_0.6-29
[79] data.table_1.12.8 png_0.1-7 vctrs_0.3.0
[82] nloptr_1.2.2.1 gtable_0.3.0 purrr_0.3.4
[85] ggplot2_3.3.1 xfun_0.13 openxlsx_4.1.5
[88] xtable_1.8-4 coda_0.19-3 Rsolnp_1.16
[91] glasso_1.11 survival_3.1-12 truncnorm_1.0-8
[94] tibble_3.0.1 arm_1.11-1 ellipse_0.4.2
[97] cluster_2.1.0 statmod_1.4.34 ellipsis_0.3.1