Mike Cheung December 24, 2022
library(metaSEM)
## Sample data from Becker (1992)
Becker92$data[1:2]## $`Berry (1957)`
## Math Spatial Verbal
## Math 1.00 0.46 0.31
## Spatial 0.46 1.00 0.19
## Verbal 0.31 0.19 1.00
##
## $`Rosenberg (1981)`
## Math Spatial Verbal
## Math 1.00 0.46 0.55
## Spatial 0.46 1.00 0.32
## Verbal 0.55 0.32 1.00
## Proposed regression model
model1 <- "## Regression paths
Math ~ Spatial2Math*Spatial + Verbal2Math*Verbal
Spatial ~~ CorSpatialVerbal*Verbal
## Fix the variances of Spatial and Verbal at 1
Spatial ~~ 1*Spatial
Verbal ~~ 1*Verbal
## Label the error variance of Math
Math ~~ ErrorVarMath*Math"
## Display the model
plot(model1)
## Convert it to RAM
RAM1 <- lavaan2RAM(model1, obs.variables=c("Math", "Spatial", "Verbal"))
RAM1## $A
## Math Spatial Verbal
## Math "0" "0.1*Spatial2Math" "0.1*Verbal2Math"
## Spatial "0" "0" "0"
## Verbal "0" "0" "0"
##
## $S
## Math Spatial Verbal
## Math "0.5*ErrorVarMath" "0" "0"
## Spatial "0" "1" "0*CorSpatialVerbal"
## Verbal "0" "0*CorSpatialVerbal" "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 1 TSSEM with random-effects model
rand1 <- tssem1(Becker92$data, Becker92$n, method="REM")
summary(rand1)##
## 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.3804588 0.0421940 0.2977601 0.4631575 9.0169 < 2.2e-16 ***
## Intercept2 0.2982765 0.0930979 0.1158080 0.4807451 3.2039 0.001356 **
## Intercept3 0.1570785 0.0503387 0.0584165 0.2557405 3.1204 0.001806 **
## Tau2_1_1 0.0020084 0.0053837 -0.0085434 0.0125602 0.3731 0.709109
## Tau2_2_2 0.0418303 0.0314264 -0.0197644 0.1034250 1.3311 0.183171
## Tau2_3_3 0.0041962 0.0099467 -0.0152989 0.0236914 0.4219 0.673118
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Q statistic on the homogeneity of effect sizes: 46.95764
## Degrees of freedom of the Q statistic: 15
## P value of the Q statistic: 3.739907e-05
##
## Heterogeneity indices (based on the estimated Tau2):
## Estimate
## Intercept1: I2 (Q statistic) 0.1928
## Intercept2: I2 (Q statistic) 0.8118
## Intercept3: I2 (Q statistic) 0.2790
##
## Number of studies (or clusters): 6
## Number of observed statistics: 18
## Number of estimated parameters: 6
## Degrees of freedom: 12
## -2 log likelihood: -18.20738
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Stage 2 TSSEM
## Two equivalent versions to calculate the R^2 and its 95% LBCI
rand2a <- tssem2(rand1, RAM=RAM1, intervals.type="LB",
mx.algebras=list(R2a=mxAlgebra(Spatial2Math^2+Verbal2Math^2
+2*CorSpatialVerbal*Spatial2Math*Verbal2Math,
name="R2a"),
## R^2 = 1 - error variance of Math
R2b=mxAlgebra(One-Smatrix[1,1], name="R2b"),
## Define a 1x1 matrix of 1
One=mxMatrix("Iden", ncol=1, nrow=1, name="One")))
summary(rand2a)##
## 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,
## mxModel.Args = mxModel.Args, subset.variables = subset.variables,
## 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.342045 NA 0.252776 0.429839 NA NA
## Verbal2Math 0.244549 NA 0.058172 0.430551 NA NA
## CorSpatialVerbal 0.157079 NA 0.058007 0.256012 NA NA
##
## mxAlgebras objects (and their 95% likelihood-based CIs):
## lbound Estimate ubound
## R2a[1,1] 0.1177919 0.2030773 0.3351317
## R2b[1,1] 0.1177919 0.2030773 0.3351317
##
## Goodness-of-fit indices:
## Value
## Sample size 538.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 90.489
## 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.)
## Display the model with the parameter estimates
plot(rand2a, color="green")
## An easier approach is to define the R2 in the model
model2 <- "## Regression paths
Math ~ Spatial2Math*Spatial + Verbal2Math*Verbal
Spatial ~~ CorSpatialVerbal*Verbal
## Fix the variances of Spatial and Verbal at 1
Spatial ~~ 1*Spatial
Verbal ~~ 1*Verbal
## Label the error variance of Math
Math ~~ ErrorVarMath*Math
## Define an R2
R2 := Spatial2Math^2+Verbal2Math^2+2*CorSpatialVerbal*Spatial2Math*Verbal2Math"
## Convert it to RAM
RAM2 <- lavaan2RAM(model2, obs.variables=c("Math", "Spatial", "Verbal"))
RAM2## $A
## Math Spatial Verbal
## Math "0" "0.1*Spatial2Math" "0.1*Verbal2Math"
## Spatial "0" "0" "0"
## Verbal "0" "0" "0"
##
## $S
## Math Spatial Verbal
## Math "0.5*ErrorVarMath" "0" "0"
## Spatial "0" "1" "0*CorSpatialVerbal"
## Verbal "0" "0*CorSpatialVerbal" "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
##
## $mxalgebras
## $mxalgebras$R2
## mxAlgebra 'R2'
## $formula: Spatial2Math^2 + Verbal2Math^2 + 2 * CorSpatialVerbal * Spatial2Math * Verbal2Math
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
rand2b <- tssem2(rand1, RAM=RAM2, intervals.type="LB")
summary(rand2b)##
## 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,
## mxModel.Args = mxModel.Args, subset.variables = subset.variables,
## 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.342045 NA 0.252776 0.429839 NA NA
## Verbal2Math 0.244549 NA 0.058172 0.430551 NA NA
## CorSpatialVerbal 0.157079 NA 0.058007 0.256012 NA NA
##
## mxAlgebras objects (and their 95% likelihood-based CIs):
## lbound Estimate ubound
## R2[1,1] 0.1177919 0.2030773 0.3351317
##
## Goodness-of-fit indices:
## Value
## Sample size 538.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 90.489
## 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 data to the right format
df <- Cor2DataFrame(Becker92)
osmasem.fit <- osmasem(RAM=RAM2, data=df, intervals.type="LB")
summary(osmasem.fit)## Summary of osmasem
##
## free parameters:
## name matrix row col Estimate Std.Error A z value
## 1 Spatial2Math A0 Math Spatial 0.3420455 0.04489421 7.618922
## 2 Verbal2Math A0 Math Verbal 0.2445485 0.09485628 2.578095
## 3 CorSpatialVerbal S0 Verbal Spatial 0.1570785 0.05033867 3.120434
## 4 Tau1_1 vecTau1 1 1 -6.2104200 2.68058695 -2.316813
## 5 Tau1_2 vecTau1 2 1 -3.1741347 0.75128419 -4.224945
## 6 Tau1_3 vecTau1 3 1 -5.4735694 2.37038631 -2.309147
## Pr(>|z|)
## 1 2.553513e-14
## 2 9.934665e-03
## 3 1.805845e-03
## 4 2.051389e-02
## 5 2.389987e-05
## 6 2.093545e-02
##
## confidence intervals:
## lbound estimate ubound note
## osmasem.Amatrix[1,1] NA 0.000000000 NA !!!
## osmasem.Amatrix[2,1] NA 0.000000000 NA !!!
## osmasem.Amatrix[3,1] NA 0.000000000 NA !!!
## osmasem.Amatrix[1,2] 0.253361009 0.342045488 0.44654265
## osmasem.Amatrix[2,2] NA 0.000000000 NA !!!
## osmasem.Amatrix[3,2] NA 0.000000000 NA !!!
## osmasem.Amatrix[1,3] 0.023298797 0.244548520 0.46816790
## osmasem.Amatrix[2,3] 0.000000000 0.000000000 NA !!!
## osmasem.Amatrix[3,3] NA 0.000000000 NA !!!
## osmasem.Smatrix[1,1] 0.636400030 0.796922700 0.88533443
## osmasem.Smatrix[2,1] NA 0.000000000 NA !!!
## osmasem.Smatrix[3,1] NA 0.000000000 NA !!!
## osmasem.Smatrix[1,2] NA 0.000000000 NA !!!
## osmasem.Smatrix[2,2] NA 1.000000000 NA !!!
## osmasem.Smatrix[3,2] 0.034120353 0.157078520 0.27569326
## osmasem.Smatrix[1,3] NA 0.000000000 NA !!!
## osmasem.Smatrix[2,3] 0.034120353 0.157078520 0.27569326
## osmasem.Smatrix[3,3] NA 1.000000000 NA !!!
## osmasem.Tau2[1,1] NA 0.002008394 0.03105435 !!!
## osmasem.Tau2[2,1] NA 0.000000000 NA !!!
## osmasem.Tau2[3,1] NA 0.000000000 NA !!!
## osmasem.Tau2[1,2] NA 0.000000000 NA !!!
## osmasem.Tau2[2,2] 0.009750432 0.041830285 0.20669836
## osmasem.Tau2[3,2] NA 0.000000000 NA !!!
## osmasem.Tau2[1,3] NA 0.000000000 NA !!!
## osmasem.Tau2[2,3] NA 0.000000000 NA !!!
## osmasem.Tau2[3,3] NA 0.004196227 0.05583177 !!!
## osmasem.R2[1,1] 0.114665573 0.203077300 0.36359997
## To investigate missing CIs, run summary() again, with verbose=T, to see CI details.
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 6 12 -18.20738
## Saturated: 9 9 NA
## Independence: 6 12 NA
## Number of observations/statistics: 538/18
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: -42.20738 -6.207377 -6.049185
## BIC: -93.66168 19.519774 0.473715
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2022-12-24 10:31:22
## Wall clock time: 0.2338519 secs
## optimizer: SLSQP
## OpenMx version number: 2.20.7
## Need help? See help(mxSummary)
plot(osmasem.fit, color="yellow")
sessionInfo()## R version 4.2.2 Patched (2022-11-10 r83330)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.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] metaSEM_1.2.6.2 OpenMx_2.20.7
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-160 RColorBrewer_1.1-3 mi_1.1
## [4] tools_4.2.2 backports_1.4.1 utf8_1.2.2
## [7] R6_2.5.1 rpart_4.1.19 Hmisc_4.7-2
## [10] colorspace_2.0-3 nnet_7.3-18 gridExtra_2.3
## [13] mnormt_2.1.1 compiler_4.2.2 qgraph_1.9.3
## [16] fdrtool_1.2.17 cli_3.5.0 htmlTable_2.4.1
## [19] scales_1.2.1 checkmate_2.1.0 mvtnorm_1.1-3
## [22] psych_2.2.9 pbapply_1.6-0 sem_3.1-15
## [25] stringr_1.5.0 digest_0.6.31 pbivnorm_0.6.0
## [28] foreign_0.8-82 minqa_1.2.5 rmarkdown_2.19
## [31] base64enc_0.1-3 jpeg_0.1-10 pkgconfig_2.0.3
## [34] htmltools_0.5.4 lme4_1.1-31 lisrelToR_0.1.5
## [37] highr_0.9 fastmap_1.1.0 htmlwidgets_1.6.0
## [40] rlang_1.0.6 rstudioapi_0.14 gtools_3.9.4
## [43] zip_2.2.2 magrittr_2.0.3 Formula_1.2-4
## [46] interp_1.1-3 Matrix_1.5-1 Rcpp_1.0.9
## [49] munsell_0.5.0 fansi_1.0.3 abind_1.4-5
## [52] rockchalk_1.8.157 lifecycle_1.0.3 stringi_1.7.8
## [55] yaml_2.3.6 carData_3.0-5 MASS_7.3-58
## [58] plyr_1.8.8 lavaan_0.6-12 grid_4.2.2
## [61] parallel_4.2.2 deldir_1.0-6 lattice_0.20-45
## [64] semPlot_1.1.6 kutils_1.70 splines_4.2.2
## [67] knitr_1.41 pillar_1.8.1 igraph_1.3.5
## [70] boot_1.3-28 corpcor_1.6.10 reshape2_1.4.4
## [73] stats4_4.2.2 XML_3.99-0.9 glue_1.6.2
## [76] evaluate_0.19 latticeExtra_0.6-30 data.table_1.14.6
## [79] RcppParallel_5.1.5 png_0.1-8 vctrs_0.5.1
## [82] nloptr_2.0.3 gtable_0.3.1 ggplot2_3.4.0
## [85] xfun_0.35 openxlsx_4.2.5.1 xtable_1.8-4
## [88] coda_0.19-4 survival_3.4-0 glasso_1.11
## [91] tibble_3.1.8 arm_1.13-1 ellipse_0.4.3
## [94] cluster_2.1.4