Load the sem package. This example also requires the DiagrammeR package for sem::pathDiagram().

if (!require("DiagrammeR")) install.packages("DiagrammeR")
library(sem)

Data

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
  ")

Model A from Joreskog & Sorbom (1984)

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
")

Fit the model using 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 =  85

Examine 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.75154

Model B

It 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 =  88

Compare models

Compare 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 ' ' 1

Path diagrams

A 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)