Selecting variables
Mike Cheung December 25, 2022
- Sometimes, we want to select some of the variables in fitting structural equation models. The following code illustrates the procedures.
library(metaSEM)
## Sample data from Digman (1997)
Digman97$data[1:2]## $`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
## Stage 1 TSSEM with random-effects model
rand1 <- tssem1(Digman97$data, Digman97$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.38971910 0.05429256 0.28330764 0.49613055 7.1781 7.068e-13
## Intercept2 0.43265879 0.04145136 0.35141562 0.51390197 10.4377 < 2.2e-16
## Intercept3 0.04945635 0.06071079 -0.06953461 0.16844731 0.8146 0.41529
## Intercept4 0.09603708 0.04478711 0.00825595 0.18381820 2.1443 0.03201
## Intercept5 0.42724243 0.03911620 0.35057609 0.50390878 10.9224 < 2.2e-16
## Intercept6 0.11929318 0.04106203 0.03881309 0.19977328 2.9052 0.00367
## Intercept7 0.19292425 0.04757962 0.09966991 0.28617859 4.0548 5.018e-05
## Intercept8 0.22690164 0.03240892 0.16338132 0.29042196 7.0012 2.538e-12
## Intercept9 0.18105567 0.04258855 0.09758364 0.26452770 4.2513 2.126e-05
## Intercept10 0.43614968 0.03205960 0.37331401 0.49898535 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.02236122 0.00995083 0.00285794 0.04186449 2.2472 0.02463
## Tau2_5_5 0.01729072 0.00796404 0.00168149 0.03289995 2.1711 0.02992
## Tau2_6_6 0.01815481 0.00895896 0.00059557 0.03571405 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.00018098 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.01049222 0.00501578 0.00066148 0.02032296 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.334
## 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.)
## Stage 2 TSSEM
## Suppose we only want to fit a one-factor model with 3 variables for illustration.
model <- "## One-factor model with 3 indicators
Alpha=~A+C+ES"
plot(model) 
RAM <- lavaan2RAM(model, obs.variables=c("A","C","ES"))
RAM## $A
## A C ES Alpha
## A "0" "0" "0" "0.1*AONAlpha"
## C "0" "0" "0" "0.1*CONAlpha"
## ES "0" "0" "0" "0.1*ESONAlpha"
## Alpha "0" "0" "0" "0"
##
## $S
## A C ES Alpha
## A "0.5*AWITHA" "0" "0" "0"
## C "0" "0.5*CWITHC" "0" "0"
## ES "0" "0" "0.5*ESWITHES" "0"
## Alpha "0" "0" "0" "1"
##
## $F
## A C ES Alpha
## A 1 0 0 0
## C 0 1 0 0
## ES 0 0 1 0
##
## $M
## A C ES Alpha
## 1 0 0 0 0
rand2 <- tssem2(rand1, RAM=RAM, subset.variables=c("A","C","ES"))
summary(rand2)##
## 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: z statistic approximation
## Coefficients:
## Estimate Std.Error lbound ubound z value Pr(>|z|)
## AONAlpha 0.628220 0.059860 0.510896 0.745543 10.495 < 2.2e-16 ***
## CONAlpha 0.620355 0.058846 0.505020 0.735690 10.542 < 2.2e-16 ***
## ESONAlpha 0.688706 0.064815 0.561670 0.815742 10.626 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Goodness-of-fit indices:
## Value
## Sample size 4496.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 257.46
## 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.)
## Display the model with the parameter estimates
plot(rand2, color="green")
## Convert the data to the right format
df <- Cor2DataFrame(Digman97)
osmasem.fit <- osmasem(RAM=RAM, data=df, subset.variables=c("A","C","ES"))
summary(osmasem.fit)## Summary of osmasem
##
## free parameters:
## name matrix row col Estimate Std.Error A z value Pr(>|z|)
## 1 AONAlpha A0 A Alpha 0.6273017 0.05952844 10.537847 0.000000e+00
## 2 CONAlpha A0 C Alpha 0.6208633 0.05863687 10.588275 0.000000e+00
## 3 ESONAlpha A0 ES Alpha 0.6887714 0.06453852 10.672253 0.000000e+00
## 4 Tau1_1 vecTau1 1 1 -3.3087018 0.41451430 -7.982118 1.332268e-15
## 5 Tau1_2 vecTau1 2 1 -3.9119737 0.43715848 -8.948640 0.000000e+00
## 6 Tau1_3 vecTau1 3 1 -4.1377834 0.46195970 -8.957022 0.000000e+00
##
## To obtain confidence intervals re-run with intervals=TRUE
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 6 36 -32.37633
## Saturated: 9 33 NA
## Independence: 6 36 NA
## Number of observations/statistics: 4496/42
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: -104.3763 -20.37633 -20.3576200
## BIC: -335.1703 18.08933 -0.9763266
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2022-12-25 11:59:22
## Wall clock time: 0.08998942 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