Load the sem package. This example also requires the DiagrammeR package for sem::pathDiagram().
if (!require("DiagrammeR")) install.packages("DiagrammeR")
library(sem)Two-wave longitudinal data on measures of Anomie and Powerless in 1967 and 1971 from Wheaton et al. (1997). These are considered measures of latent factors for Alienation, and a main question is whether an Alienation construct remains the same over time.
Education and Socioeconomic Index (SEI) are considered measures of a latent factor, SES, which may predict Alientation at both time points.
S.wheaton <- readMoments(
names=c('Anomia67','Powerless67',
'Anomia71','Powerless71',
'Education','SEI'),
text="
11.834
6.947 9.364
6.819 5.091 12.532
4.783 5.028 7.495 9.986
-3.839 -3.889 -3.841 -3.625 9.610
-21.899 -18.831 -21.748 -18.775 35.522 450.288
")This model posits three latent variables of two variables each, and with a reference solution setting the coefficient of the first indicator to 1. Alienation67 and Alienation71 are each predicted by the latent SES variable, but Alienation71 is considered to depend on Alienation67 as well.
wh.model.A <- specifyEquations(text="
Anomia67 = 1*Alienation67
Powerless67 = lamy1*Alienation67
Anomia71 = 1*Alienation71
Powerless71 = lamy2*Alienation71
Education = 1*SES
SEI = lamx*SES
Alienation67 = gam1*SES
Alienation71 = gam2*SES + beta*Alienation67
V(Anomia67) = the1
V(Anomia71) = the2
V(Powerless67) = the3
V(Powerless71) = the4
V(SES) = phi
")sem::sem()sem.wh.A <- sem(wh.model.A, S.wheaton, N=932)Examine fit indices. This model doesn’t fit very well. Why not?
summary(sem.wh.A, fit.indices = c("RMSEA", "NNFI", "CFI"))##
## Model Chisquare = 71.46973 Df = 6 Pr(>Chisq) = 2.041707e-13
## RMSEA index = 0.1082604 90% CI: (0.08658466, 0.1314454)
## Tucker-Lewis NNFI = 0.922665
## Bentler CFI = 0.969066
##
## Normalized Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.2578633 -0.2117727 -0.0000127 -0.0153411 0.2443546 1.3307183
##
## R-square for Endogenous Variables
## Alienation67 Anomia67 Powerless67 Alienation71 Anomia71
## 0.3212 0.6607 0.6592 0.5763 0.7047
## Powerless71 Education SEI
## 0.6370 0.6936 0.4204
##
## Parameter Estimates
## Estimate Std Error z value Pr(>|z|)
## lamy1 0.8885364 0.04150033 21.410348 1.070118e-101
## lamy2 0.8487223 0.03995708 21.240851 4.005674e-100
## lamx 5.3289571 0.42976842 12.399601 2.626270e-35
## gam1 -0.6138170 0.05645164 -10.873326 1.544659e-27
## gam2 -0.1741787 0.05391357 -3.230702 1.234864e-03
## beta 0.7047276 0.05353754 13.163242 1.428235e-39
## the1 4.0155181 0.34315626 11.701719 1.248976e-31
## the2 3.7010811 0.37341979 9.911315 3.717211e-23
## the3 3.1913382 0.27145244 11.756528 6.536739e-32
## the4 3.6248251 0.29208013 12.410379 2.295635e-35
## phi 6.6658511 0.64105526 10.398247 2.525357e-25
## V[Alienation67] 5.3069765 0.47260170 11.229279 2.928570e-29
## V[Alienation71] 3.7412397 0.38756392 9.653220 4.763613e-22
## V[Education] 2.9441577 0.49980006 5.890671 3.846307e-09
## V[SEI] 260.9929854 18.24177314 14.307435 1.966326e-46
##
## lamy1 Powerless67 <--- Alienation67
## lamy2 Powerless71 <--- Alienation71
## lamx SEI <--- SES
## gam1 Alienation67 <--- SES
## gam2 Alienation71 <--- SES
## beta Alienation71 <--- Alienation67
## the1 Anomia67 <--> Anomia67
## the2 Anomia71 <--> Anomia71
## the3 Powerless67 <--> Powerless67
## the4 Powerless71 <--> Powerless71
## phi SES <--> SES
## V[Alienation67] Alienation67 <--> Alienation67
## V[Alienation71] Alienation71 <--> Alienation71
## V[Education] Education <--> Education
## V[SEI] SEI <--> SEI
##
## Iterations = 85Examine modification indices (A: regression coef.; P: covariances)
print(modIndices(sem.wh.A), n.largest=3)##
## 3 largest modification indices, A matrix (regression coefficients):
## Anomia71<-Anomia67 Anomia67<-Anomia71 Powerless71<-Anomia67
## 58.72444 51.91211 46.15612
##
## 3 largest modification indices, P matrix (variances/covariances):
## Anomia71<->Anomia67 Powerless71<->Anomia67 Anomia71<->Powerless67
## 63.70627 49.82948 49.75154It makes sense that the error terms for Anomia and/or Powerlessness could be correlated across years. The largest MI is the covariance for Anomia71<->Anomia67 – set it free.
In the sem package, the update() function makes it easy to add (or remove) parameters. <-> specifies a covariance. You could do the same using copy/paste of wh.model.A and adding a line C(Anomia67, Anomia71) = the13.
wh.model.B <- update(wh.model.A, text="
add, Anomia67 <-> Anomia71, the13"
)Fit model B: This model fits very well, in absolute terms
sem.wh.B <- sem(wh.model.B, S.wheaton, N=932)
summary(sem.wh.B, fit.indices = c("RMSEA", "NNFI", "CFI"))##
## Model Chisquare = 6.330708 Df = 5 Pr(>Chisq) = 0.2753563
## RMSEA index = 0.01690758 90% CI: (NA, 0.05090479)
## Tucker-Lewis NNFI = 0.9981137
## Bentler CFI = 0.9993712
##
## Normalized Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.9566754 -0.1342772 -0.0426406 -0.0917168 -0.0000002 0.5456945
##
## R-square for Endogenous Variables
## Alienation67 Anomia67 Powerless67 Alienation71 Anomia71
## 0.3065 0.5725 0.7635 0.5009 0.6168
## Powerless71 Education SEI
## 0.7313 0.7160 0.4073
##
## Parameter Estimates
## Estimate Std Error z value Pr(>|z|)
## lamy1 1.0265323 0.05313420 19.319614 3.673908e-83
## lamy2 0.9709171 0.04941846 19.646852 6.151147e-86
## lamx 5.1627509 0.42135535 12.252724 1.624460e-34
## gam1 -0.5497509 0.05341119 -10.292805 7.593158e-25
## gam2 -0.2114319 0.04926346 -4.291861 1.771817e-05
## beta 0.6173408 0.04970423 12.420288 2.028311e-35
## the1 5.0652661 0.37106898 13.650470 2.005802e-42
## the2 4.8120224 0.39524465 12.174795 4.234586e-34
## the3 2.2145647 0.31763205 6.972107 3.122290e-12
## the4 2.6830185 0.32988524 8.133188 4.181466e-16
## phi 6.8804462 0.65795344 10.457345 1.356019e-25
## V[Alienation67] 4.7051856 0.43313877 10.862998 1.729752e-27
## V[Alienation71] 3.8663366 0.34348374 11.256244 2.157612e-29
## V[Education] 2.7295609 0.51625560 5.287228 1.241840e-07
## V[SEI] 266.8968034 18.19343380 14.669952 1.004158e-48
## the13 1.8873901 0.24000157 7.864074 3.718375e-15
##
## lamy1 Powerless67 <--- Alienation67
## lamy2 Powerless71 <--- Alienation71
## lamx SEI <--- SES
## gam1 Alienation67 <--- SES
## gam2 Alienation71 <--- SES
## beta Alienation71 <--- Alienation67
## the1 Anomia67 <--> Anomia67
## the2 Anomia71 <--> Anomia71
## the3 Powerless67 <--> Powerless67
## the4 Powerless71 <--> Powerless71
## phi SES <--> SES
## V[Alienation67] Alienation67 <--> Alienation67
## V[Alienation71] Alienation71 <--> Alienation71
## V[Education] Education <--> Education
## V[SEI] SEI <--> SEI
## the13 Anomia71 <--> Anomia67
##
## Iterations = 88Compare the two models by a likelihood-ratio test, using the sem::anova() function.
anova(sem.wh.A, sem.wh.B)## LR Test for Difference Between Models
##
## Model Df Model Chisq Df LR Chisq Pr(>Chisq)
## sem.wh.A 6 71.470
## sem.wh.B 5 6.331 1 65.139 6.98e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A simple version, just showing the parameters in the model. For simplicity, this omits the covariance between Anomia in 67 and 71.
For some unknown reason, the graphs do not appear in this output.
pathDiagram(sem.wh.B,
same.rank=c("Alienation67, Alienation71"),
min.rank=c("Education", "SEI"))A version showing the standardized fitted coefficients, with edge color and thickness indicating the sign and magnitude of the coefficient.
pathDiagram(sem.wh.B,
same.rank=c("Alienation67, Alienation71"),
min.rank=c("Education", "SEI"),
edge.labels = "values",
edge.colors = c("blue", "red"),
node.colors = c("pink", "lightblue1"),
edge.weight="proportional",
standardize=TRUE)