Basics of Using MIIVsem to Estimate
Structural Equation Models
Overview
MIIVsem is an R package to estimate and test Structural Equation Models (SEMs) using Model Implied Instrumental Variables (MIIVs). These notes are oriented toward using the package not to explaining the method. See references at end for the basis for the method.
To use MIIVsem requires R and the lavaan package. The first part of the instructions explain how to load and use lavaan to generate its estimates. This is followed by the commands for loading and using MIIVsem. In these notes only some basic commands of MIIVsem for simultaneous equation models are shown.
Data Management.
Load the dataset.
data <- read.csv("~/Desktop/EducData.txt", header=FALSE)Add the variable names
names(data) <- c( "prestg80", "age", "educ", "female", "black",
"asian", "hispanic", "othrace", "napaeduc", "paeducM",
"namaeduc", "maeducM", "paNOeduc", "maNOeduc", "napapr80",
"papr80M", "namapr80", "mapr80M", "lnsibs" )Preview the data
You can preview the data using the head() function from base R.
head(data)## pr80 ag ed fm bl as hs ot np pM nm mM pN mN np80 p80M nm80 m8 ln
## 1 65 31 16 0 0 0 1 0 0 8.0 0 3 0 0 0 31.0 0 35.0 1.1
## 2 32 71 8 1 0 0 1 0 1 11.5 0 6 0 0 1 43.9 1 41.5 1.8
## 3 42 40 6 0 1 0 0 0 0 8.0 0 8 0 0 0 45.0 1 41.5 2.5
## 4 51 46 16 1 1 0 0 0 0 12.0 0 16 0 0 0 45.0 0 66.0 1.1
## 5 66 80 15 1 1 0 0 0 0 9.0 0 12 0 0 0 40.0 0 23.0 1.1
## 6 25 31 14 1 1 0 0 0 0 8.0 0 8 0 0 0 29.0 0 47.0 1.9Model 1
Model 1 Path Diagram
Code for producing this path diagram can be found in Appendix A.
Syntax for Model 1
model1 <- 'educ ~ paeducM + maeducM + paNOeduc + maNOeduc + napaeduc + namaeduc +
papr80M + mapr80M + napapr80 + namapr80 + age + lnsibs +
female + black + asian + hispanic + othrace
prestg80 ~ educ 'Additional Syntax Examples
=~ can be read as “measured by”
~ can be read as “regressed on”
~~ can be read as “covariance with”
== can be read as “is set equal to”
Maximum Likelihood (ML) Analysis using lavaan
Load lavaan into the R workspace.
library("lavaan")Estimate Model 1 using lavaan
First, we save the estimated fit for model1 as model1_fit.
model1_fit <- sem(model = model1, data = data, meanstructure = TRUE)Parameter Estimates
We can save the parameter estimates of model1_fit using the parameterEstimates() function from lavaan.
model1_params <- parameterEstimates(model1_fit)
write.csv(model1_params, "~/Desktop/model1_params.csv")Summary Statistics
Most of the time, the summary function with additional argument fit.measures = TRUE will be sufficient.
summary(model1_fit, fit.measures = TRUE) ## lavaan (0.5-18) converged normally after 115 iterations
##
## Number of observations 3921
##
## Estimator ML
## Minimum Function Test Statistic 161.856
## Degrees of freedom 17
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 2687.274
## Degrees of freedom 35
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.945
## Tucker-Lewis Index (TLI) 0.888
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -83181.002
## Loglikelihood unrestricted model (H1) -83100.074
##
## Number of free parameters 22
## Akaike (AIC) 166406.003
## Bayesian (BIC) 166544.033
## Sample-size adjusted Bayesian (BIC) 166474.127
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.047
## 90 Percent Confidence Interval 0.040 0.053
## P-value RMSEA <= 0.05 0.791
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.013
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|)
## Regressions:
## educ ~
## paeducM 0.128 0.016 7.939 0.000
## maeducM 0.165 0.018 9.157 0.000
## paNOeduc -1.183 0.397 -2.982 0.003
## maNOeduc 0.360 0.389 0.925 0.355
## napaeduc -2.555 0.913 -2.799 0.005
## namaeduc -1.179 0.223 -5.289 0.000
## papr80M 0.023 0.004 5.750 0.000
## mapr80M 0.008 0.004 2.105 0.035
## napapr80 1.891 0.916 2.066 0.039
## namapr80 -0.004 0.099 -0.042 0.967
## age 0.012 0.003 4.412 0.000
## lnsibs -0.660 0.072 -9.118 0.000
## female 0.075 0.081 0.931 0.352
## black -0.087 0.129 -0.675 0.499
## asian 1.481 0.265 5.587 0.000
## hispanic -0.847 0.152 -5.574 0.000
## othrace -0.807 0.355 -2.272 0.023
## prestg80 ~
## educ 2.445 0.065 37.808 0.000
##
## Intercepts:
## educ 9.568 0.321 29.762 0.000
## prestg80 10.912 0.917 11.901 0.000
##
## Variances:
## educ 6.250 0.141
## prestg80 143.024 3.230MIIV-2SLS Analysis using MIIVsem
Load MIIVsem.
library("MIIVsem")We can obtain MIIV-2SLS estimates for Model 1 using the miive() function.
miive(model = model1, data = data)##
## MIIVsem results
##
## DV EV Estimate StdErr Z P(|Z|) Sargan df P(Chi)
## educ Int 9.5681 0.3215 29.762 0.00
## paeducM 0.1278 0.0161 7.939 0.00
## maeducM 0.1649 0.0180 9.157 0.00
## paNOeduc -1.1834 0.3968 -2.982 0.00
## maNOeduc 0.3598 0.3889 0.925 0.35
## napaeduc -2.5550 0.9129 -2.799 0.01
## namaeduc -1.1785 0.2228 -5.289 0.00
## papr80M 0.0234 0.0041 5.750 0.00
## mapr80M 0.0082 0.0039 2.105 0.04
## napapr80 1.8913 0.9155 2.066 0.04
## namapr80 -0.0041 0.0985 -0.042 0.97
## age 0.0122 0.0028 4.412 0.00
## lnsibs -0.6603 0.0724 -9.118 0.00
## female 0.0751 0.0806 0.931 0.35
## black -0.0874 0.1294 -0.675 0.50
## asian 1.4812 0.2651 5.587 0.00
## hispanic -0.8470 0.1520 -5.574 0.00
## othrace -0.8065 0.3550 -2.272 0.02
## prestg80 Int 10.9116 0.9169 11.901 0.00 159 17 0
## educ 2.4450 0.0647 37.808 0.00Similarly, we can save the MIIV-2SLS parameter estimates as a data.frame or .csv file.
model1_miiv <- miive(model = model1, data = data)$dat
write.csv(model1_miiv, "~/Desktop/model1_miiv.csv")Model 1: Maximum Likelihood and MIIV-2SLS Estimates
It is convenient to put the maximum likelihood estimates from lavaan side-by-side with the MIIV-2SLS estimates of MIIVsem. Appendix B to these slides gives the commands that were used to create this table.
## DV EV lavaan Model 1 MIIVsem Model 1
## 1 educ age 0.012 (0.003) 0.012 (0.003)
## 2 educ asian 1.481 (0.265) 1.481 (0.265)
## 3 educ black -0.087 (0.129) -0.087 (0.129)
## 4 educ female 0.075 (0.081) 0.075 (0.081)
## 5 educ hispanic -0.847 (0.152) -0.847 (0.152)
## 6 educ lnsibs -0.66 (0.072) -0.66 (0.072)
## 7 educ maeducM 0.165 (0.018) 0.165 (0.018)
## 8 educ maNOeduc 0.36 (0.389) 0.36 (0.389)
## 9 educ mapr80M 0.008 (0.004) 0.008 (0.004)
## 10 educ namaeduc -1.179 (0.223) -1.179 (0.223)
## 11 educ namapr80 -0.004 (0.099) -0.004 (0.099)
## 12 educ napaeduc -2.555 (0.913) -2.555 (0.913)
## 13 educ napapr80 1.891 (0.916) 1.891 (0.916)
## 14 educ othrace -0.807 (0.355) -0.807 (0.355)
## 15 educ paeducM 0.128 (0.016) 0.128 (0.016)
## 16 educ paNOeduc -1.183 (0.397) -1.183 (0.397)
## 17 educ papr80M 0.023 (0.004) 0.023 (0.004)
## 18 prestg80 educ 2.445 (0.065) 2.445 (0.065)## Likelihood Chi Square Sargan Chi Square
## Model 1 161.86, df = 17, p = 0 ---------------------
## Model 1: educ equation --------------------- 0 (exactly identified)
## Model 1: prestg80 equation --------------------- 158.56, df = 17, p = 0Model 2
We can also respecify Model 1 with a correlation between the errors of educ and prestg80.
Model 2 Path Diagram
Code for producing this path diagram can be found in Appendix A.
Syntax for Model 2
model2 <- 'educ ~ paeducM + maeducM + paNOeduc + maNOeduc + napaeduc + namaeduc +
papr80M + mapr80M + napapr80 + namapr80 + age + lnsibs +
female + black + asian + hispanic + othrace
prestg80 ~ educ
prestg80 ~~ educ 'Model 2: Maximum Likelihood and MIIV-2SLS Estimates
We can again compare the maximum likelihood estimates from lavaan with the MIIV-2SLS estimates of MIIVsem. Appendix B to these slides gives the commands used to generate this table.
## DV EV lavaan Model 2 MIIVsem Model 2
## 1 educ age 0.015 (0.003) 0.012 (0.003)
## 2 educ asian 1.527 (0.263) 1.481 (0.265)
## 3 educ black -0.158 (0.128) -0.087 (0.129)
## 4 educ female 0.079 (0.08) 0.075 (0.081)
## 5 educ hispanic -0.817 (0.151) -0.847 (0.152)
## 6 educ lnsibs -0.66 (0.072) -0.66 (0.072)
## 7 educ maeducM 0.163 (0.018) 0.165 (0.018)
## 8 educ maNOeduc 0.349 (0.385) 0.36 (0.389)
## 9 educ mapr80M 0.009 (0.004) 0.008 (0.004)
## 10 educ namaeduc -1.166 (0.221) -1.179 (0.223)
## 11 educ namapr80 -0.003 (0.098) -0.004 (0.099)
## 12 educ napaeduc -2.385 (0.905) -2.555 (0.913)
## 13 educ napapr80 1.714 (0.907) 1.891 (0.916)
## 14 educ othrace -0.743 (0.352) -0.807 (0.355)
## 15 educ paeducM 0.127 (0.016) 0.128 (0.016)
## 16 educ paNOeduc -1.057 (0.393) -1.183 (0.397)
## 17 educ papr80M 0.025 (0.004) 0.023 (0.004)
## 18 prestg80 educ 2.915 (0.122) 2.855 (0.122)## Likelihood Chi Square Sargan Chi Square
## Model 2 143.58, df = 16, p = 0 ---------------------
## Model 2: educ equation --------------------- 0 (exactly identified)
## Model 2: prestg80 equation --------------------- 141.19, df = 16, p = 0Additional Features
MIIVsem includess the Industrialization-Democracy dataset from Bollen (1989). Users can view the dataset by typing bollen1989a into the R console. We use this example to demonstrate some additional capabilities of MIIVsem.
head(bollen1989a)## y1 y2 y3 y4 y5 y6 y7 y8 x1 x2 x3
## 1 2.50 0.000 3.333 0.000 1.250 0.000 3.726 3.333 4.443 3.638 2.558
## 2 1.25 0.000 3.333 0.000 6.250 1.100 6.667 0.737 5.384 5.063 3.568
## 3 7.50 8.800 10.000 9.200 8.750 8.094 10.000 8.212 5.961 6.256 5.224
## 4 8.90 8.800 10.000 9.200 8.908 8.128 10.000 4.615 6.286 7.568 6.267
## 5 10.00 3.333 10.000 6.667 7.500 3.333 10.000 6.667 5.864 6.819 4.574
## 6 7.50 3.333 6.667 6.667 6.250 1.100 6.667 0.368 5.533 5.136 3.892Industrialization-Democracy Example
Path Diagram
Model Specification
model <- ' ### "=~" can be read as "measured by"
Eta1 =~ y1 + y2 + y3 + y4
Xi1 =~ x1 + x2 + x3
Eta2 =~ y5 + y6 + y7 + y8
### "~" can be read as "regressed on"
Eta1 ~ Xi1
Eta2 ~ Xi1
Eta2 ~ Eta1
### "~~" can be read as "covariance with"
y1 ~~ y5
y2 ~~ y4
y2 ~~ y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'Estimation with Equality Constraints
Cross-equation equality constraints can be placed on the coefficients using labels in the model synatx.
model_con <-
' ### "=~" can be read as "measured by"
Eta1 =~ y1 + label1*y2 + label2*y3 + label3*y4
Xi1 =~ x1 + x2 + x3
Eta2 =~ y5 + label1*y6 + label2*y7 + label3*y8
### "~" can be read as "regressed on"
Eta1 ~ Xi1
Eta2 ~ Xi1
Eta2 ~ Eta1
### "~~" can be read as "covariance with"
y1 ~~ y5
y2 ~~ y4
y2 ~~ y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'miive(model = model_con, data = bollen1989a)##
## MIIVsem results
##
## DV EV Estimate StdErr Z P(|Z|) Sargan df P(Chi)
## Eta1 Int -0.91 2.17 -0.42 0.68 0.5 1 0.48
## Xi1 1.26 0.43 2.96 0.00
## Eta2 Int -4.50 1.42 -3.16 0.00 0.8 3 0.85
## Eta1 0.72 0.10 7.14 0.00
## Xi1 1.12 0.31 3.60 0.00
## y2 Int -1.73 0.76 -2.28 0.02 8.6 5 0.13
## Eta1 1.10 0.12 9.01 0.00
## y3 Int 0.70 0.63 1.12 0.26 6.2 6 0.41
## Eta1 1.07 0.10 10.49 0.00
## y4 Int -2.14 0.64 -3.37 0.00 4.3 5 0.51
## Eta1 1.21 0.10 11.69 0.00
## y6 Int -2.65 0.71 -3.73 0.00 8.6 5 0.13
## Eta2 1.10 0.12 9.01 0.00
## y7 Int 0.69 0.60 1.15 0.25 10.7 6 0.10
## Eta2 1.07 0.10 10.49 0.00
## y8 Int -2.15 0.62 -3.49 0.00 2.8 5 0.73
## Eta2 1.21 0.10 11.69 0.00
## x2 Int -5.71 0.65 -8.73 0.00 8.3 8 0.40
## Xi1 2.08 0.13 16.17 0.00
## x3 Int -5.29 0.76 -6.99 0.00 8.7 8 0.36
## Xi1 1.75 0.15 11.78 0.00
##
## Tests of the Linear Restrictions
##
## Restrictions:
## y2 = y6 1
## y3 = y7 1
## y4 = y8 1
##
## Wald Test: Asymptotic Chi-squared
## df Chi-sq P(Chi)
## 3 1.2 0.76Model Implied Instrumental Variable (MIIV) Search
With the model specified, the miivs() function can be used to identify MIIVs. Additionally, by setting the miivs.out argument to TRUE you can print the list of MIIVs as an object to be used in future analyses.
miivs(model = model, miivs.out = TRUE)## Model Equation Information
##
## LHS RHS Composite Disturbance MIIVs
## y1 x1 e.x1, e.y1, e.Eta1 x2, x3
## y5 y1, x1 e.y1, e.x1, e.y5, e.Eta2 y2, y3, y4, x2, x3
## y2 y1 e.y1, e.y2 y3, y7, y8, x2, x3, x1
## y3 y1 e.y1, e.y3 y2, y4, y6, y8, x2, x3, x1
## y4 y1 e.y1, e.y4 y3, y6, y7, x2, x3, x1
## y6 y5 e.y5, e.y6 y3, y4, y7, x2, x3, x1
## y7 y5 e.y5, e.y7 y2, y4, y6, y8, x2, x3, x1
## y8 y5 e.y5, e.y8 y2, y3, y7, x2, x3, x1
## x2 x1 e.x2, e.x1 y1, y5, y2, y3, y4, y6, y7, y8, x3
## x3 x1 e.x3, e.x1 y1, y5, y2, y3, y4, y6, y7, y8, x2
##
##
## instruments <- '
## y1 ~ x2 + x3
## y5 ~ y2 + y3 + y4 + x2 + x3
## y2 ~ y3 + y7 + y8 + x2 + x3 + x1
## y3 ~ y2 + y4 + y6 + y8 + x2 + x3 + x1
## y4 ~ y3 + y6 + y7 + x2 + x3 + x1
## y6 ~ y3 + y4 + y7 + x2 + x3 + x1
## y7 ~ y2 + y4 + y6 + y8 + x2 + x3 + x1
## y8 ~ y2 + y3 + y7 + x2 + x3 + x1
## x2 ~ y1 + y5 + y2 + y3 + y4 + y6 + y7 + y8 + x3
## x3 ~ y1 + y5 + y2 + y3 + y4 + y6 + y7 + y8 + x2
## 'Bootstrap Standard Errors
MIIVsem provides two options for obtaining bootstrap standard errors, with accompanying significant tests and p-values: (1) the pairs and (2) residual bootstrap.
miive(model = model, data = bollen1989a, bootstrap = "pairs")##
## MIIVsem results
##
## DV EV Estimate BtSD P(t) Sargan df P(Chi)
## Eta1 Int -0.91 1.970 0.01 0.5 1 0.48
## Xi1 1.26 0.396 0.02
## Eta2 Int -4.50 1.356 0.02 0.8 3 0.85
## Eta1 0.72 0.096 0.00
## Xi1 1.12 0.290 0.01
## y2 Int -1.97 0.663 0.04 8.4 5 0.14
## Eta1 1.14 0.128 0.01
## y3 Int 1.27 0.877 0.04 5.9 6 0.44
## Eta1 0.97 0.136 0.00
## y4 Int -2.16 0.702 0.06 4.3 5 0.51
## Eta1 1.21 0.118 0.00
## y6 Int -2.42 0.749 0.04 8.7 5 0.12
## Eta2 1.05 0.142 0.00
## y7 Int 0.14 0.876 0.06 9.5 6 0.15
## Eta2 1.18 0.137 0.00
## y8 Int -2.14 0.813 0.09 2.8 5 0.73
## Eta2 1.20 0.146 0.00
## x2 Int -5.71 0.695 0.08 8.3 8 0.40
## Xi1 2.08 0.132 0.00
## x3 Int -5.29 0.726 0.03 8.7 8 0.36
## Xi1 1.75 0.139 0.00References
Bollen, K. A. 1996. An Alternative 2SLS Estimator for Latent Variable Models. Psychometrika, 61, 109-121.
Bollen, K. A. 2001. Two-stage Least Squares and Latent Variable Models: Simultaneous Estimation and Robustness to Misspecifications. In R. Cudeck, S. Du Toit, and D. Sorbom (Eds.), Structural Equation Modeling: Present and Future, A Festschrift in Honor of Karl Joreskog (pp. 119-138). Lincoln, IL: Scientific Software.
Bollen, K. A. and D. J. Bauer. 2004. Automating the Selection of Model-Implied Instrumental Variables. Sociological Methods and Research, 32, 425-52.
Appendix A
Code for producing Model 1 path diagram.
library(DiagrammeR)
grViz("
digraph SEM {
graph [layout = neato,
overlap = FALSE,
outputorder = nodesfirst]
node [shape = rectangle, fixedsize = TRUE,
height = 1.6, fontsize = 24, color = 'blue',
style ='filled', fontcolor = 'white']
z1 [pos = '0, 5!', label = 'age', shape = square]
z2 [pos = '2, 5!', label = 'female', shape = square]
z3 [pos = '4, 5!', label = 'black', shape = square]
z4 [pos = '6, 5!', label = 'asian', shape = square]
z5 [pos = '8, 5!', label = 'hispanic', shape = square]
z6 [pos = '10, 5!', label = 'othrace', shape = square]
z7 [pos = '12, 5!', label = 'napaeduc', shape = square]
z8 [pos = '14, 5!', label = 'paeducM', shape = square]
z9 [pos = '16, 5!', label = 'namaeduc', shape = square]
z10 [pos = '1, -5!', label = 'maeducM', shape = square]
z11 [pos = '3, -5!', label = 'paNOeduc', shape = square]
z12 [pos = '5,-5!', label = 'maNOeduc', shape = square]
z13 [pos = '7,-5!', label = 'napapr80', shape = square]
z14 [pos = '9,-5!', label = 'papr80M', shape = square]
z15 [pos = '11,-5!', label = 'namapr80', shape = square]
z16 [pos = '13,-5!', label = 'mapr80M', shape = square]
z17 [pos = '15,-5!', label = 'lnsibs', shape = square]
y1 [pos = '8, 0!', label = 'educ', shape = square]
y2 [pos = '12,0!', label = 'prestg80', shape = square]
z1->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z2->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z3->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z4->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z5->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z6->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z7->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z8->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z9->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z10->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z11->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z12->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z13->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z14->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z15->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z16->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z17->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
y1 ->y2 [label = '', headport = 'w',penwidth = 1.5, arrowsize = 2]
node [shape = circle, fixedsize = TRUE,
height = .8, fontsize = 12, color = 'blue',
style ='filled', fontcolor = 'white']
e1 [pos = '10, -1!', label = 'e1', shape = circle]
e2 [pos = '14, -1!', label = 'e2', shape = circle]
e1->y1 [label = '', tailport = 'nw',headport = 'e', penwidth = 1.5, arrowsize = 2]
e2->y2 [label = '', tailport = 'nw',headport = 'e', penwidth = 1.5, arrowsize = 2]
}
")Code for producing Model 2 path diagram.
library(DiagrammeR)
grViz("
digraph SEM {
graph [layout = neato,
overlap = FALSE,
outputorder = nodesfirst]
node [shape = rectangle, fixedsize = TRUE,
height = 1.6, fontsize = 24, color = 'blue',
style ='filled', fontcolor = 'white']
z1 [pos = '0, 5!', label = 'age', shape = square]
z2 [pos = '2, 5!', label = 'female', shape = square]
z3 [pos = '4, 5!', label = 'black', shape = square]
z4 [pos = '6, 5!', label = 'asian', shape = square]
z5 [pos = '8, 5!', label = 'hispanic', shape = square]
z6 [pos = '10, 5!', label = 'othrace', shape = square]
z7 [pos = '12, 5!', label = 'napaeduc', shape = square]
z8 [pos = '14, 5!', label = 'paeducM', shape = square]
z9 [pos = '16, 5!', label = 'namaeduc', shape = square]
z10 [pos = '1, -5!', label = 'maeducM', shape = square]
z11 [pos = '3, -5!', label = 'paNOeduc', shape = square]
z12 [pos = '5,-5!', label = 'maNOeduc', shape = square]
z13 [pos = '7,-5!', label = 'napapr80', shape = square]
z14 [pos = '9,-5!', label = 'papr80M', shape = square]
z15 [pos = '11,-5!', label = 'namapr80', shape = square]
z16 [pos = '13,-5!', label = 'mapr80M', shape = square]
z17 [pos = '15,-5!', label = 'lnsibs', shape = square]
y1 [pos = '8, 0!', label = 'educ', shape = square]
y2 [pos = '12,0!', label = 'prestg80', shape = square]
z1->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z2->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z3->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z4->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z5->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z6->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z7->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z8->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z9->y1 [label = '', tailport = 's',headport = 'n',penwidth = 1.5]
z10->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z11->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z12->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z13->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z14->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z15->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z16->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
z17->y1 [label = '', tailport = 'n',headport = 's',penwidth = 1.5]
y1 ->y2 [label = '', headport = 'w',penwidth = 1.5, arrowsize = 2]
node [shape = circle, fixedsize = TRUE,
height = .8, fontsize = 12, color = 'blue',
style ='filled', fontcolor = 'white']
e1 [pos = '10, -1!', label = 'e1', shape = circle]
e2 [pos = '14, -1!', label = 'e2', shape = circle]
e1->y1 [label = '', tailport = 'nw',headport = 'e', penwidth = 1.5, arrowsize = 2]
e2->y2 [label = '', tailport = 'nw',headport = 'e', penwidth = 1.5, arrowsize = 2]
e1->e2 [style = 'invis', tailport = 'se', headport = 'sw',
dir = both, penwidth = 1.5, arrowsize = 2]
e1->e2 [tailport = 'se', headport = 'sw', dir = both, penwidth = 1.5, arrowsize = 2]
}
")Appendix B
Code for producing Model 1 comparison table.
## lavaan
model1_fit <- sem(model = model1, data = data, meanstructure = TRUE)
model1_params <- parameterEstimates(model1_fit)
lava1 <- model1_params[model1_params$op == "~",]
lava1 <- lava1 [,c("lhs", "rhs", "est", "se")]
lava1[,3:4] <- round(lava1[,3:4], 3)
lava1$L <- paste0(lava1[,3], " (", lava1[,4], ") ")
lava1 <- lava1[,c("lhs", "rhs", "L")]
colnames(lava1) <- c("DV", "EV", "L")
## MIIVsem
model1_miiv <- miive(model = model1, data = data)$dat
miiv1 <- model1_miiv
miiv1 <- miiv1[,c("DV", "EV", "Estimate", "StdErr")]
miiv1[,3:4] <- round(miiv1[,3:4], 3)
miiv1$M <- paste0(miiv1[,3], " (", miiv1[,4], ") ")
miiv1 <- miiv1[,c("DV", "EV", "M")]
miiv1 <- miiv1[miiv1 $EV != "Int",]
## merge lava1 and miiv1
mod1 <- merge(lava1, miiv1, by = c("DV", "EV"), all.x = TRUE, all.y = TRUE)
colnames(mod1) <- c("DV", "EV", "lavaan Model 1", "MIIVsem Model 1")
## extract lavaan fit estimates
chi <- round(fitMeasures(model1_fit)["chisq"], 2)
chidf <- fitMeasures(model1_fit)["df"]
chip <- round(fitMeasures(model1_fit)["pvalue"], 2)
## extract MIIVsem fit estimates
sarg <- round(model1_miiv[19, "Sargan"],2)
sargdf <- model1_miiv[19, "df"]
sargp <- round(model1_miiv[19, "P(Chi)"], 2)
## create comparison table for display
likchi_lavaan <- paste0(chi, ", df = ", chidf, ", p = ", chip)
eq1sarg_lavaan <- "---------------------"
eq2sarg_lavaan <- "---------------------"
likchi_miiv <- "---------------------"
eq1sarg_miiv <- paste0("0 (exactly identified)")
eq2sarg_miiv <- paste0(sarg, ", df = ", sargdf, ", p = ", sargp)
col1 <- rbind(likchi_lavaan, eq1sarg_lavaan, eq2sarg_lavaan)
col2 <- rbind(likchi_miiv, eq1sarg_miiv, eq2sarg_miiv)
fit0 <- cbind(col1, col2)
rownames(fit0) <- c("Model 1", "Model 1: educ equation", "Model 1: prestg80 equation")
colnames(fit0) <- c("Likelihood Chi Square", "Sargan Chi Square")
print(mod1)
print(fit0, quote = FALSE, right = FALSE, row.names = TRUE)Code for producing Model 2 comparison table.
## lavaan
model2_fit <- sem(model = model2, data = data, meanstructure = TRUE)
model2_params <- parameterEstimates(model2_fit)
lava2 <- model2_params[model2_params$op == "~",]
lava2 <- lava2[,c("lhs", "rhs", "est", "se")]
lava2[,3:4] <- round(lava2[,3:4], 3)
lava2$L <- paste0(lava2[,3], " (", lava2[,4], ") ")
lava2 <- lava2[,c("lhs", "rhs", "L")]
colnames(lava2) <- c("DV", "EV", "L")
## MIIVsem
model2_miiv <- miive(model = model2, data = data)$dat
miiv2 <- model2_miiv
miiv2 <- miiv2[,c("DV", "EV", "Estimate", "StdErr")]
miiv2[,3:4] <- round(miiv2[,3:4], 3)
miiv2$M <- paste0(miiv2[,3], " (", miiv2[,4], ") ")
miiv2 <- miiv2[,c("DV", "EV", "M")]
miiv2 <- miiv2[miiv2$EV != "Int",]
## merge lava2 and miiv2
mod2 <- merge(lava2, miiv2, by = c("DV", "EV"), all.x = TRUE, all.y = TRUE)
colnames(mod2) <- c("DV", "EV", "lavaan Model 2", "MIIVsem Model 2")
## extract lavaan fit estimates
chi <- round(fitMeasures(model2_fit)["chisq"], 2)
chidf <- fitMeasures(model2_fit)["df"]
chip <- round(fitMeasures(model2_fit)["pvalue"], 2)
## extract MIIVsem fit estimates
sarg <- round(model2_miiv[19, "Sargan"],2)
sargdf <- model2_miiv[19, "df"]
sargp <- round(model2_miiv[19, "P(Chi)"], 2)
## create comparison table for display
likchi_lavaan <- paste0(chi, ", df = ", chidf, ", p = ", chip)
eq1sarg_lavaan <- "---------------------"
eq2sarg_lavaan <- "---------------------"
likchi_miiv <- "---------------------"
eq1sarg_miiv <- paste0("0 (exactly identified)")
eq2sarg_miiv <- paste0(sarg, ", df = ", sargdf, ", p = ", sargp)
col1 <- rbind(likchi_lavaan, eq1sarg_lavaan, eq2sarg_lavaan)
col2 <- rbind(likchi_miiv, eq1sarg_miiv, eq2sarg_miiv)
fit0 <- cbind(col1, col2)
rownames(fit0) <- c("Model 2", "Model 2: educ equation", "Model 2: prestg80 equation")
colnames(fit0) <- c("Likelihood Chi Square", "Sargan Chi Square")
print(mod2)
print(fit0, quote = FALSE, right = FALSE, row.names = TRUE)