#' --- #' title: "Wheaton et al. Stability of Alienation: sem package" #' author: "Michael Friendly" #' date: "22 May 2019" #' --- #+ include=FALSE knitr::opts_chunk$set(message=FALSE, warning=FALSE, width=90) #' 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")) #' Examine modification indices (**A**: regression coef.; **P**: covariances) print(modIndices(sem.wh.A), n.largest=3) #' ## **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")) #' ## Compare models #' Compare the two models by a likelihood-ratio test, using the `sem::anova()` function. anova(sem.wh.A, sem.wh.B) #' ## 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)