# Strukturgleichungsmodell
#
# Beispiel aus Fox, J. (2006). Structural equation modeling with the sem
# package in R. Structural Equation Modeling, 13, 465-486.
#
# ==============================================================================

# Laden des Pakets "sem" mittels library() (muss zuvor installiert werden):
library(sem)

# ------------------------------------------------------------------------------

# Kovarianzmatrix der Variablen (die Funktion sem() erlaubt es, Matrizen in nur
# in der unteren (oder oberen) Dreiecksform zu spezifizieren):

S.wheaton = matrix(c( 11.834,      0,      0,      0,      0,      0,
                       6.947,  9.364,      0,      0,      0,      0,
                       6.819,  5.091, 12.532,      0,      0,      0,
                       4.783,  5.028,  7.495,  9.986,      0,      0,
                      -3.839, -3.889, -3.841, -3.625,  9.610,      0,
                     -21.899,-18.831,-21.748,-18.775, 35.522,450.288),
                   6, 6, byrow=T)
rownames(S.wheaton)=colnames(S.wheaton)=c('Anomia67','Powerless67','Anomia71',
                                          'Powerless71','Education','SEI')

# Modellspezifikation in RAM-Notation (Ladungen werden hier mit L (Loading),
# Messfehler (bzw. Residuen) der beobachteten Variablen mit E (Error),
# Zusammenhänge zwischen latenten Variablen mit C (Contingency) und Residuen
# der latenten Variablen mit D (Disturbance) bezeichnet, wobei die Nummern
# sich auf die Nummer der Variable bzw. des Faktors beziehen):

mod.wheaton.1 = specify.model()
  Anomia67     <-  Alienation67,   NA,  1
  Powerless67  <-  Alienation67, L021, NA
  Anomia71     <-  Alienation71,   NA,  1
  Powerless71  <-  Alienation71, L042, NA
  Education    <-  SES         ,   NA,  1
  SEI          <-  SES         , L063, NA
  Alienation67 <-  SES         ,  C13, NA
  Alienation71 <-  Alienation67,  C21, NA
  Alienation71 <-  SES         ,  C23, NA
  Anomia67     <-> Anomia67    ,   E1, NA
  Powerless67  <-> Powerless67 ,   E2, NA
  Anomia71     <-> Anomia71    ,   E3, NA
  Powerless71  <-> Powerless71 ,   E4, NA
  Education    <-> Education   ,   E5, NA
  SEI          <-> SEI         ,   E6, NA
  Alienation67 <-> Alienation67,   D1, NA
  Alienation71 <-> Alienation71,   D2, NA
  SES          <-> SES         ,   D3, NA

# Schätzen der Modellparameter:
sem.wheaton.1 = sem(mod.wheaton.1, S.wheaton, N=932)
summary(sem.wheaton.1)
std.coef(sem.wheaton.1)
mod.indices(sem.wheaton.1)

# Die Modifikationsindizes weisen auf eine hohe Korrelation der Messfehler von
# Anomia und Powerless über die Zeit hin, was theoretisch Sinn macht. Eine ent-
# sprechende Modifikation des Modells liefert:
mod.wheaton.2 = specify.model()
  Anomia67     <-  Alienation67,   NA,  1
  Powerless67  <-  Alienation67, L021, NA
  Anomia71     <-  Alienation71,   NA,  1
  Powerless71  <-  Alienation71, L042, NA
  Education    <-  SES         ,   NA,  1
  SEI          <-  SES         , L063, NA
  Alienation67 <-  SES         ,  C13, NA
  Alienation71 <-  Alienation67,  C21, NA
  Alienation71 <-  SES         ,  C23, NA
  Anomia67     <-> Anomia67    ,   E1, NA
  Powerless67  <-> Powerless67 ,   E2, NA
  Anomia71     <-> Anomia71    ,   E3, NA
  Powerless71  <-> Powerless71 ,   E4, NA
  Education    <-> Education   ,   E5, NA
  SEI          <-> SEI         ,   E6, NA
  Anomia67     <-> Anomia71    , E1E3, NA
  Powerless67  <-> Powerless71 , E2E4, NA
  Alienation67 <-> Alienation67,   D1, NA
  Alienation71 <-> Alienation71,   D2, NA
  SES          <-> SES         ,   D3, NA

sem.wheaton.2 = sem(mod.wheaton.2, S.wheaton, N=932)
summary(sem.wheaton.2)
std.coef(sem.wheaton.2)
mod.indices(sem.wheaton.2)

# Auch wenn die Kovarianz zwischen den Messfehlern von Powerless nicht signi-
# fikant ist, sollte sie hier aus theoretischen Gründen beibehalten werden. Wie
# der Comparative-Fit-Index (CFI) and der Root-Mean-Square-Error-of-Approximat-
# ion (RMSEA), passt das Gesamtmodell sehr gut zu den Daten.