inst/scripts/SEM.R
In martscht/PsyMSc1: Übungen für die Seminare im Sommersemester 2020
#### R-Skript zur 4. Sitzung zu PFadanalysen und SEM ####
# R Chunks wurden aus dem HTML exportiert und sind dort alle enthalten
#### Vorbereitung ----
library(lavaan)
library(semPlot) # grafische Darstellung von Pfadanalyse- und Strukturgleichungsmodellen
# Datensatz laden
data("StressAtWork", package = "PsyMSc1")
head(StressAtWork) # gebe die ersten 6 Zeilen aus
names(StressAtWork) # Spalten Namen des Datensatzes
# Skalenmittelwerte bestimmen und dem Datensatz anhängen
StressAtWork$ZDs <- rowMeans(StressAtWork[,paste0("zd",c(1, 2, 6))])
StressAtWork$BOEEs <- rowMeans(StressAtWork[,paste0("bo",c(1, 6, 12, 19))])
StressAtWork$BFs <- rowMeans(StressAtWork[,paste0("bf",1:20)])
##### Pfadanalyse -----
# Model für Skalenmittelwerte (fitte mit sem):
model_paths <- '
BOEEs ~ ZDs
BFs ~ BOEEs + ZDs
'
# So sähe das Modell in lavaan aus (fitte mit lavaan):
model_paths_lavaan <- '
BOEEs ~ ZDs
BFs ~ BOEEs + ZDs
BOEEs ~~ BOEEs
BFs ~~ BFs
'
# Fit und summary mit sem
fit_paths <- sem(model_paths, data = StressAtWork)
summary(fit_paths, rsq = T, fit.measures = T)
## Grafische Veranschaulichung des Modells und der Ergebnisse ---
# Defaults: Modellstruktur:
semPaths(fit_paths)
# Ergebnisse zeichnen:
semPaths(fit_paths, what = "est")
# weitere Ploteinstellungen:
semPaths(object = fit_paths, what = "model", layout = "tree2", rotation = 2,
col = list(man = "skyblue"), edge.label.cex=1, sizeMan = 5)
semPaths(object = fit_paths, what = "est", layout = "tree2", rotation = 2,
col = list(man = "skyblue"), edge.label.cex=1, sizeMan = 5)
## Berechnen und Testen des indirekten Effektes von Zeitdruck über Emotionale Erschöpfung auf Psychosomatische Beschwerden ---
# Pfadmodell mit indirektem und totalem Effekt:
model_paths_IE_TE <- '
BOEEs ~ a*ZDs
BFs ~ b*BOEEs + c*ZDs
# Neue Parameter
IE := a*b
TE := IE + c
'
# Schätzen und Summary des Pfadmodells mit indirektem und totalem Effekt
fit_paths_IE_TE <- sem(model_paths_IE_TE, data = StressAtWork)
summary(fit_paths_IE_TE) # Achtung!
# Leider können wir den angezeigten Standardfehlern und dem zugehörigen p-Wert nicht unbedingt vertrauen,
# da einige Studien gezeigt haben, dass der indirekte Effekt asymptotisch nicht immer normalverteilt bzw.
# symmetrisch verteilt ist, weswegen ein einfach Teilen des Estimates durch SE und der Vergleich mit der
# z-Verteilung (so wie wir dies eigentlich immer tun; also die einfache Parametersignifikanzentscheidung
# in komplexen Modellen) in diesem Fall nicht sinnvoll erscheint. Aus diesem Grund möchten wir dieser
# Problematik entgehen, indem wir die Bootstrapp-Methode verwenden.
## Bootstrapping ---
set.seed(1234)
fit_paths_IE_TE_boot <- sem(model_paths_IE_TE, data = StressAtWork, se = "boot", bootstrap = 1000)
parameterEstimates(fit_paths_IE_TE_boot, ci = TRUE)
parameterEstimates(fit_paths_IE_TE_boot, ci = TRUE)[7:8,] # Zeilen im Output von Interesse
##### Strukturgleichungsmodelle -----
## Das Modell ---
model_sem_IE_TE <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ a*ZD
BFs ~ b*BOEE + c*ZD
# Neue Parameter
IE := a*b
TE := IE + c
'
fit_sem_IE_TE <- sem(model_sem_IE_TE, StressAtWork)
summary(fit_sem_IE_TE, fit.measures = T, rsq = T)
## Berechnen und Testen des indirekten Effektes von Zeitdruck über Emotionale Erschöpfung auf Psychosomatische Beschwerden im SEM ---
set.seed(12345)
fit_sem_IE_TE_boot <- sem(model_sem_IE_TE, data = StressAtWork, se = "boot", bootstrap = 1000)
parameterEstimates(fit_sem_IE_TE_boot, ci = TRUE)
parameterEstimates(fit_sem_IE_TE_boot, ci = TRUE)[21:22,] # Parameter von Interesse
## Grafische Repräsentation ---
semPaths(object = fit_sem_IE_TE, what = "model", layout = "tree2",
rotation = 2, curve = T, col = list(man = "skyblue", lat = "yellow"),
curvePivot = T, edge.label.cex=1.2, sizeMan = 5, sizeLat = 8)
semPaths(object = fit_sem_IE_TE, what = "est", layout = "tree2",
rotation = 2, curve = T, col = list(man = "skyblue", lat = "yellow"),
curvePivot = T, edge.label.cex=1.2, sizeMan = 5, sizeLat = 8)
##### MSA -----
# Modell mit 2 Gruppen: Sex
fit_sem_IE_TE_MSA <- sem(model_sem_IE_TE, data = StressAtWork, group = "sex")
summary(fit_sem_IE_TE_MSA)
# Vollständig benanntes Modell mit 2 indirekten und 2 totalen Effekten:
model_sem_IE_TE_MSA_abc <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ c(a1, a2)*ZD
BFs ~ c(b1, b2)*BOEE + c(c1,c2)*ZD
# Neue Parameter
IE1 := a1*b1
TE1 := IE1 + c1
IE2 := a2*b2
TE2 := IE2 + c2
'
fit_sem_IE_TE_MSA_abc <- sem(model_sem_IE_TE_MSA_abc, StressAtWork, group = "sex")
summary(fit_sem_IE_TE_MSA_abc)
## Invarianztestung über das Geschlecht ---
# entfernen der a,b,c Benennung - dieses Modell wird für die Invarianztestung verwendet!
model_sem <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
'
## konfigurale Invarianz
fit_sem_sex_konfigural <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c(""), group.partial = c("BFs ~ 1", "BFs ~~*BFs"))
summary(fit_sem_sex_konfigural, fit.measures = T)
## metrische Invarianz
fit_sem_sex_metrisch <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
summary(fit_sem_sex_metrisch, fit.measures = T)
# Modellvergleich (metrisch vs. konfigural):
lavTestLRT(fit_sem_sex_metrisch, fit_sem_sex_konfigural) # --> metrisch wird NICHT verworfen
## Skalare Invarianz
fit_sem_sex_skalar <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
summary(fit_sem_sex_skalar, fit.measures = T)
# Modellvergleich (skalar vs. metrisch):
lavTestLRT(fit_sem_sex_skalar, fit_sem_sex_metrisch) # --> skalar wird NICHT verworfen
## Strikte Invarianz
fit_sem_sex_strikt <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
# Modellvergleich (strikt vs. skalar):
lavTestLRT(fit_sem_sex_strikt, fit_sem_sex_skalar) # --> strikt wird NICHT verworfen
## Vollständige Invarianz
fit_sem_sex_voll <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals",
"means", "lv.variances", "lv.covariances", "regressions"))
# Modellvergleich (vollständig vs. strikt):
lavTestLRT(fit_sem_sex_voll, fit_sem_sex_strikt) # --> vollständig WIRD verworfen (keine vollständige Invarianz über das Geschlecht)
# suchen des Effekts im strikt-invarianten Modell:
summary(fit_sem_sex_strikt)
#### Appendix A -----------
# Populationsmodell zur Simulation
pop_model_H0 <- '
# Messmodelle
Xi1 =~ x1 + 0.7*x2 + 0.6*x3
Eta1 =~ y1 + 0.8*y2
Eta2 =~ y3 + 0.9*y4
# Strukturmodell
Eta1 ~ 0.5*Xi1
Eta2 ~ 0.54*Xi1 + 0.4*Eta1
# Fehlerkovarianzen
x1 ~~ 0.4*x2
'
set.seed(123456)
data <- simulateData(model = pop_model_H0, meanstructure = F, sample.nobs = 200)
head(data)
# Modelle zum Schätzen
model_H0 <- '
# Messmodelle
Xi1 =~ x1 + x2 + x3
Eta1 =~ y1 + y2
Eta2 =~ y3 + y4
# Strukturmodell
Eta1 ~ Xi1
Eta2 ~ Xi1 + Eta1
# Fehlerkovarianzen
x1 ~~ x2
'
fit_H0 <- sem(model = model_H0, data = data)
summary(fit_H0)
semPaths(fit_H0, curve = T, curvePivot = T)
fitmeasures(fit_H0, c("chisq", 'df', "pvalue"))
# H1 Modell: keine Fehlerkovarianz
model_H1_kov <- '
# Messmodelle
Xi1 =~ x1 + x2 + x3
Eta1 =~ y1 + y2
Eta2 =~ y3 + y4
# Strukturmodell
Eta1 ~ Xi1
Eta2 ~ Xi1 + Eta1
'
semPaths(sem(model_H1_kov, data))
# H1 Modell: keine direkte Beziehung
model_H1_Struk <- '
# Messmodelle
Xi1 =~ x1 + x2 + x3
Eta1 =~ y1 + y2
Eta2 =~ y3 + y4
# Strukturmodell
Eta1 ~ Xi1
Eta2 ~ Eta1
# Fehlerkovarianzen
x1 ~~ x2
'
semPaths(sem(model_H1_Struk, data))
# Modellvergleich:
fit_H1_kov <- sem(model_H1_kov, data)
fit_H1_Struk <- sem(model_H1_Struk, data)
fitmeasures(fit_H0, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_kov, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_Struk, c("chisq", 'df', "pvalue"))
# Neuer Versuch für N = 1000
set.seed(123456)
data <- simulateData(model = pop_model_H0, meanstructure = F, sample.nobs = 1000)
fit_H0 <- sem(model = model_H0, data = data)
fit_H1_kov <- sem(model_H1_kov, data)
fit_H1_Struk <- sem(model_H1_Struk, data)
fitmeasures(fit_H0, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_kov, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_Struk, c("chisq", 'df', "pvalue"))
## ezCutoffs:
library(ezCutoffs)
ezCutoffs(model = model_sem, data = StressAtWork)
#### Appendix B -----------
# MSA zu Fuß
model_sem <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
'
# Vergleiche händische Modelle mit jenen zuvor!
## Konfigural ---
fit_sem_sex_konfigural <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c(""), group.partial = c("BFs ~ 1", "BFs ~~*BFs"))
fit_sem_sex_konfigural2 <- sem(model_sem, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_konfigural, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_konfigural2, c("chisq", 'df', "pvalue")) # --> gleich
## Metrisch ---
model_sem_metrisch <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD'
fit_sem_sex_metrisch <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
fit_sem_sex_metrisch2 <- sem(model_sem_metrisch, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_metrisch, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_metrisch2, c("chisq", 'df', "pvalue")) # --> gleich
## Skalar ---
model_sem_skalar <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
zd1 ~ c(tau1, tau1)*1
zd2 ~ c(tau2, tau2)*1
zd6 ~ c(tau3, tau3)*1
bo1 ~ c(tau4, tau4)*1
bo6 ~ c(tau5, tau5)*1
bo12 ~ c(tau6, tau6)*1
bo19 ~ c(tau7, tau7)*1
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
BOEE ~ c(0, NA)*1
ZD ~ c(0, NA)*1
'
fit_sem_sex_skalar <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
fit_sem_sex_skalar2 <- sem(model_sem_skalar, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_skalar, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_skalar2, c("chisq", 'df', "pvalue")) # --> gleich
## Strikt ---
model_sem_strikt <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
zd1 ~ c(tau1, tau1)*1
zd2 ~ c(tau2, tau2)*1
zd6 ~ c(tau3, tau3)*1
bo1 ~ c(tau4, tau4)*1
bo6 ~ c(tau5, tau5)*1
bo12 ~ c(tau6, tau6)*1
bo19 ~ c(tau7, tau7)*1
zd1 ~~ c(t1, t1)*zd1
zd2 ~~ c(t2, t2)*zd2
zd6 ~~ c(t3, t3)*zd6
bo1 ~~ c(t4, t4)*bo1
bo6 ~~ c(t5, t5)*bo6
bo12 ~~ c(t6, t6)*bo12
bo19 ~~ c(t7, t7)*bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
BOEE ~ c(0, NA)*1
ZD ~ c(0, NA)*1
'
fit_sem_sex_strikt <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
fit_sem_sex_strikt2 <- sem(model_sem_strikt, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_strikt, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_strikt2, c("chisq", 'df', "pvalue")) # --> gleich
## Vollständig ---
model_sem_voll <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
zd1 ~ c(tau1, tau1)*1
zd2 ~ c(tau2, tau2)*1
zd6 ~ c(tau3, tau3)*1
bo1 ~ c(tau4, tau4)*1
bo6 ~ c(tau5, tau5)*1
bo12 ~ c(tau6, tau6)*1
bo19 ~ c(tau7, tau7)*1
zd1 ~~ c(t1, t1)*zd1
zd2 ~~ c(t2, t2)*zd2
zd6 ~~ c(t3, t3)*zd6
bo1 ~~ c(t4, t4)*bo1
bo6 ~~ c(t5, t5)*bo6
bo12 ~~ c(t6, t6)*bo12
bo19 ~~ c(t7, t7)*bo19
# Strukturmodell
BOEE ~ c(a, a)*ZD
BFs ~ c(b, b)*BOEE + c(c, c)*ZD
BOEE ~ c(0, 0)*1
ZD ~ c(0, 0)*1
BFs ~ c(kappa, kappa)*1
ZD ~~ c(phi, phi)*ZD
BOEE ~~ c(psi1, psi1)*BOEE
BFs ~~ c(psi2, psi2)*BFs
'
fit_sem_sex_voll <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals",
"means", "lv.variances", "lv.covariances", "regressions"))
fit_sem_sex_voll2 <- sem(model_sem_voll, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_voll, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_voll2, c("chisq", 'df', "pvalue")) # --> gleich
#### Appendix C -----------
model_sem_IE_TE <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ a*ZD
BFs ~ b*BOEE + c*ZD
# Neue Parameter
IE := a*b
TE := IE + c
'
model_sem_IE_TE_tau <- '
# Messmodelle
ZD =~ l1*zd1 + l1*zd2 + l1*zd6
BOEE =~ l2*bo1 + l2*bo6 + l2*bo12 + l2*bo19
# Strukturmodell
BOEE ~ a*ZD
BFs ~ b*BOEE + c*ZD
# Neue Parameter
IE := a*b
TE := IE + c
'
fit_sem_IE_TE_tau <- sem(model_sem_IE_TE_tau, StressAtWork)
summary(fit_sem_IE_TE_tau, fit.measures = T, rsq = T)
# Modellvergleich
lavTestLRT(fit_sem_IE_TE_tau, fit_sem_IE_TE) # essentielle Tau-Äquivalenz wird verworfen!
martscht/PsyMSc1 documentation built on July 11, 2020, 3:45 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#### R-Skript zur 4. Sitzung zu PFadanalysen und SEM ####
# R Chunks wurden aus dem HTML exportiert und sind dort alle enthalten
#### Vorbereitung ----
library(lavaan)
library(semPlot) # grafische Darstellung von Pfadanalyse- und Strukturgleichungsmodellen
# Datensatz laden
data("StressAtWork", package = "PsyMSc1")
head(StressAtWork) # gebe die ersten 6 Zeilen aus
names(StressAtWork) # Spalten Namen des Datensatzes
# Skalenmittelwerte bestimmen und dem Datensatz anhängen
StressAtWork$ZDs <- rowMeans(StressAtWork[,paste0("zd",c(1, 2, 6))])
StressAtWork$BOEEs <- rowMeans(StressAtWork[,paste0("bo",c(1, 6, 12, 19))])
StressAtWork$BFs <- rowMeans(StressAtWork[,paste0("bf",1:20)])
##### Pfadanalyse -----
# Model für Skalenmittelwerte (fitte mit sem):
model_paths <- '
BOEEs ~ ZDs
BFs ~ BOEEs + ZDs
'
# So sähe das Modell in lavaan aus (fitte mit lavaan):
model_paths_lavaan <- '
BOEEs ~ ZDs
BFs ~ BOEEs + ZDs
BOEEs ~~ BOEEs
BFs ~~ BFs
'
# Fit und summary mit sem
fit_paths <- sem(model_paths, data = StressAtWork)
summary(fit_paths, rsq = T, fit.measures = T)
## Grafische Veranschaulichung des Modells und der Ergebnisse ---
# Defaults: Modellstruktur:
semPaths(fit_paths)
# Ergebnisse zeichnen:
semPaths(fit_paths, what = "est")
# weitere Ploteinstellungen:
semPaths(object = fit_paths, what = "model", layout = "tree2", rotation = 2,
col = list(man = "skyblue"), edge.label.cex=1, sizeMan = 5)
semPaths(object = fit_paths, what = "est", layout = "tree2", rotation = 2,
col = list(man = "skyblue"), edge.label.cex=1, sizeMan = 5)
## Berechnen und Testen des indirekten Effektes von Zeitdruck über Emotionale Erschöpfung auf Psychosomatische Beschwerden ---
# Pfadmodell mit indirektem und totalem Effekt:
model_paths_IE_TE <- '
BOEEs ~ a*ZDs
BFs ~ b*BOEEs + c*ZDs
# Neue Parameter
IE := a*b
TE := IE + c
'
# Schätzen und Summary des Pfadmodells mit indirektem und totalem Effekt
fit_paths_IE_TE <- sem(model_paths_IE_TE, data = StressAtWork)
summary(fit_paths_IE_TE) # Achtung!
# Leider können wir den angezeigten Standardfehlern und dem zugehörigen p-Wert nicht unbedingt vertrauen,
# da einige Studien gezeigt haben, dass der indirekte Effekt asymptotisch nicht immer normalverteilt bzw.
# symmetrisch verteilt ist, weswegen ein einfach Teilen des Estimates durch SE und der Vergleich mit der
# z-Verteilung (so wie wir dies eigentlich immer tun; also die einfache Parametersignifikanzentscheidung
# in komplexen Modellen) in diesem Fall nicht sinnvoll erscheint. Aus diesem Grund möchten wir dieser
# Problematik entgehen, indem wir die Bootstrapp-Methode verwenden.
## Bootstrapping ---
set.seed(1234)
fit_paths_IE_TE_boot <- sem(model_paths_IE_TE, data = StressAtWork, se = "boot", bootstrap = 1000)
parameterEstimates(fit_paths_IE_TE_boot, ci = TRUE)
parameterEstimates(fit_paths_IE_TE_boot, ci = TRUE)[7:8,] # Zeilen im Output von Interesse
##### Strukturgleichungsmodelle -----
## Das Modell ---
model_sem_IE_TE <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ a*ZD
BFs ~ b*BOEE + c*ZD
# Neue Parameter
IE := a*b
TE := IE + c
'
fit_sem_IE_TE <- sem(model_sem_IE_TE, StressAtWork)
summary(fit_sem_IE_TE, fit.measures = T, rsq = T)
## Berechnen und Testen des indirekten Effektes von Zeitdruck über Emotionale Erschöpfung auf Psychosomatische Beschwerden im SEM ---
set.seed(12345)
fit_sem_IE_TE_boot <- sem(model_sem_IE_TE, data = StressAtWork, se = "boot", bootstrap = 1000)
parameterEstimates(fit_sem_IE_TE_boot, ci = TRUE)
parameterEstimates(fit_sem_IE_TE_boot, ci = TRUE)[21:22,] # Parameter von Interesse
## Grafische Repräsentation ---
semPaths(object = fit_sem_IE_TE, what = "model", layout = "tree2",
rotation = 2, curve = T, col = list(man = "skyblue", lat = "yellow"),
curvePivot = T, edge.label.cex=1.2, sizeMan = 5, sizeLat = 8)
semPaths(object = fit_sem_IE_TE, what = "est", layout = "tree2",
rotation = 2, curve = T, col = list(man = "skyblue", lat = "yellow"),
curvePivot = T, edge.label.cex=1.2, sizeMan = 5, sizeLat = 8)
##### MSA -----
# Modell mit 2 Gruppen: Sex
fit_sem_IE_TE_MSA <- sem(model_sem_IE_TE, data = StressAtWork, group = "sex")
summary(fit_sem_IE_TE_MSA)
# Vollständig benanntes Modell mit 2 indirekten und 2 totalen Effekten:
model_sem_IE_TE_MSA_abc <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ c(a1, a2)*ZD
BFs ~ c(b1, b2)*BOEE + c(c1,c2)*ZD
# Neue Parameter
IE1 := a1*b1
TE1 := IE1 + c1
IE2 := a2*b2
TE2 := IE2 + c2
'
fit_sem_IE_TE_MSA_abc <- sem(model_sem_IE_TE_MSA_abc, StressAtWork, group = "sex")
summary(fit_sem_IE_TE_MSA_abc)
## Invarianztestung über das Geschlecht ---
# entfernen der a,b,c Benennung - dieses Modell wird für die Invarianztestung verwendet!
model_sem <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
'
## konfigurale Invarianz
fit_sem_sex_konfigural <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c(""), group.partial = c("BFs ~ 1", "BFs ~~*BFs"))
summary(fit_sem_sex_konfigural, fit.measures = T)
## metrische Invarianz
fit_sem_sex_metrisch <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
summary(fit_sem_sex_metrisch, fit.measures = T)
# Modellvergleich (metrisch vs. konfigural):
lavTestLRT(fit_sem_sex_metrisch, fit_sem_sex_konfigural) # --> metrisch wird NICHT verworfen
## Skalare Invarianz
fit_sem_sex_skalar <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
summary(fit_sem_sex_skalar, fit.measures = T)
# Modellvergleich (skalar vs. metrisch):
lavTestLRT(fit_sem_sex_skalar, fit_sem_sex_metrisch) # --> skalar wird NICHT verworfen
## Strikte Invarianz
fit_sem_sex_strikt <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
# Modellvergleich (strikt vs. skalar):
lavTestLRT(fit_sem_sex_strikt, fit_sem_sex_skalar) # --> strikt wird NICHT verworfen
## Vollständige Invarianz
fit_sem_sex_voll <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals",
"means", "lv.variances", "lv.covariances", "regressions"))
# Modellvergleich (vollständig vs. strikt):
lavTestLRT(fit_sem_sex_voll, fit_sem_sex_strikt) # --> vollständig WIRD verworfen (keine vollständige Invarianz über das Geschlecht)
# suchen des Effekts im strikt-invarianten Modell:
summary(fit_sem_sex_strikt)
#### Appendix A -----------
# Populationsmodell zur Simulation
pop_model_H0 <- '
# Messmodelle
Xi1 =~ x1 + 0.7*x2 + 0.6*x3
Eta1 =~ y1 + 0.8*y2
Eta2 =~ y3 + 0.9*y4
# Strukturmodell
Eta1 ~ 0.5*Xi1
Eta2 ~ 0.54*Xi1 + 0.4*Eta1
# Fehlerkovarianzen
x1 ~~ 0.4*x2
'
set.seed(123456)
data <- simulateData(model = pop_model_H0, meanstructure = F, sample.nobs = 200)
head(data)
# Modelle zum Schätzen
model_H0 <- '
# Messmodelle
Xi1 =~ x1 + x2 + x3
Eta1 =~ y1 + y2
Eta2 =~ y3 + y4
# Strukturmodell
Eta1 ~ Xi1
Eta2 ~ Xi1 + Eta1
# Fehlerkovarianzen
x1 ~~ x2
'
fit_H0 <- sem(model = model_H0, data = data)
summary(fit_H0)
semPaths(fit_H0, curve = T, curvePivot = T)
fitmeasures(fit_H0, c("chisq", 'df', "pvalue"))
# H1 Modell: keine Fehlerkovarianz
model_H1_kov <- '
# Messmodelle
Xi1 =~ x1 + x2 + x3
Eta1 =~ y1 + y2
Eta2 =~ y3 + y4
# Strukturmodell
Eta1 ~ Xi1
Eta2 ~ Xi1 + Eta1
'
semPaths(sem(model_H1_kov, data))
# H1 Modell: keine direkte Beziehung
model_H1_Struk <- '
# Messmodelle
Xi1 =~ x1 + x2 + x3
Eta1 =~ y1 + y2
Eta2 =~ y3 + y4
# Strukturmodell
Eta1 ~ Xi1
Eta2 ~ Eta1
# Fehlerkovarianzen
x1 ~~ x2
'
semPaths(sem(model_H1_Struk, data))
# Modellvergleich:
fit_H1_kov <- sem(model_H1_kov, data)
fit_H1_Struk <- sem(model_H1_Struk, data)
fitmeasures(fit_H0, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_kov, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_Struk, c("chisq", 'df', "pvalue"))
# Neuer Versuch für N = 1000
set.seed(123456)
data <- simulateData(model = pop_model_H0, meanstructure = F, sample.nobs = 1000)
fit_H0 <- sem(model = model_H0, data = data)
fit_H1_kov <- sem(model_H1_kov, data)
fit_H1_Struk <- sem(model_H1_Struk, data)
fitmeasures(fit_H0, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_kov, c("chisq", 'df', "pvalue"))
fitmeasures(fit_H1_Struk, c("chisq", 'df', "pvalue"))
## ezCutoffs:
library(ezCutoffs)
ezCutoffs(model = model_sem, data = StressAtWork)
#### Appendix B -----------
# MSA zu Fuß
model_sem <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
'
# Vergleiche händische Modelle mit jenen zuvor!
## Konfigural ---
fit_sem_sex_konfigural <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c(""), group.partial = c("BFs ~ 1", "BFs ~~*BFs"))
fit_sem_sex_konfigural2 <- sem(model_sem, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_konfigural, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_konfigural2, c("chisq", 'df', "pvalue")) # --> gleich
## Metrisch ---
model_sem_metrisch <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD'
fit_sem_sex_metrisch <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
fit_sem_sex_metrisch2 <- sem(model_sem_metrisch, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_metrisch, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_metrisch2, c("chisq", 'df', "pvalue")) # --> gleich
## Skalar ---
model_sem_skalar <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
zd1 ~ c(tau1, tau1)*1
zd2 ~ c(tau2, tau2)*1
zd6 ~ c(tau3, tau3)*1
bo1 ~ c(tau4, tau4)*1
bo6 ~ c(tau5, tau5)*1
bo12 ~ c(tau6, tau6)*1
bo19 ~ c(tau7, tau7)*1
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
BOEE ~ c(0, NA)*1
ZD ~ c(0, NA)*1
'
fit_sem_sex_skalar <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
fit_sem_sex_skalar2 <- sem(model_sem_skalar, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_skalar, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_skalar2, c("chisq", 'df', "pvalue")) # --> gleich
## Strikt ---
model_sem_strikt <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
zd1 ~ c(tau1, tau1)*1
zd2 ~ c(tau2, tau2)*1
zd6 ~ c(tau3, tau3)*1
bo1 ~ c(tau4, tau4)*1
bo6 ~ c(tau5, tau5)*1
bo12 ~ c(tau6, tau6)*1
bo19 ~ c(tau7, tau7)*1
zd1 ~~ c(t1, t1)*zd1
zd2 ~~ c(t2, t2)*zd2
zd6 ~~ c(t3, t3)*zd6
bo1 ~~ c(t4, t4)*bo1
bo6 ~~ c(t5, t5)*bo6
bo12 ~~ c(t6, t6)*bo12
bo19 ~~ c(t7, t7)*bo19
# Strukturmodell
BOEE ~ ZD
BFs ~ BOEE + ZD
BOEE ~ c(0, NA)*1
ZD ~ c(0, NA)*1
'
fit_sem_sex_strikt <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals"), group.partial = c("BFs ~ 1", "BFs ~~BFs"))
fit_sem_sex_strikt2 <- sem(model_sem_strikt, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_strikt, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_strikt2, c("chisq", 'df', "pvalue")) # --> gleich
## Vollständig ---
model_sem_voll <- '
# Messmodelle
ZD =~ zd1 + c(l1, l1)*zd2 + c(l2, l2)*zd6
BOEE =~ bo1 + c(l3,l3)*bo6 + c(l4, l4)*bo12 + c(l5, l5)*bo19
zd1 ~ c(tau1, tau1)*1
zd2 ~ c(tau2, tau2)*1
zd6 ~ c(tau3, tau3)*1
bo1 ~ c(tau4, tau4)*1
bo6 ~ c(tau5, tau5)*1
bo12 ~ c(tau6, tau6)*1
bo19 ~ c(tau7, tau7)*1
zd1 ~~ c(t1, t1)*zd1
zd2 ~~ c(t2, t2)*zd2
zd6 ~~ c(t3, t3)*zd6
bo1 ~~ c(t4, t4)*bo1
bo6 ~~ c(t5, t5)*bo6
bo12 ~~ c(t6, t6)*bo12
bo19 ~~ c(t7, t7)*bo19
# Strukturmodell
BOEE ~ c(a, a)*ZD
BFs ~ c(b, b)*BOEE + c(c, c)*ZD
BOEE ~ c(0, 0)*1
ZD ~ c(0, 0)*1
BFs ~ c(kappa, kappa)*1
ZD ~~ c(phi, phi)*ZD
BOEE ~~ c(psi1, psi1)*BOEE
BFs ~~ c(psi2, psi2)*BFs
'
fit_sem_sex_voll <- sem(model_sem, data = StressAtWork, group = "sex",
group.equal = c("loadings", "intercepts", "residuals",
"means", "lv.variances", "lv.covariances", "regressions"))
fit_sem_sex_voll2 <- sem(model_sem_voll, data = StressAtWork, group = "sex")
# chi^2, df, p-Wert
fitmeasures(fit_sem_sex_voll, c("chisq", 'df', "pvalue"))
fitmeasures(fit_sem_sex_voll2, c("chisq", 'df', "pvalue")) # --> gleich
#### Appendix C -----------
model_sem_IE_TE <- '
# Messmodelle
ZD =~ zd1 + zd2 + zd6
BOEE =~ bo1 + bo6 + bo12 + bo19
# Strukturmodell
BOEE ~ a*ZD
BFs ~ b*BOEE + c*ZD
# Neue Parameter
IE := a*b
TE := IE + c
'
model_sem_IE_TE_tau <- '
# Messmodelle
ZD =~ l1*zd1 + l1*zd2 + l1*zd6
BOEE =~ l2*bo1 + l2*bo6 + l2*bo12 + l2*bo19
# Strukturmodell
BOEE ~ a*ZD
BFs ~ b*BOEE + c*ZD
# Neue Parameter
IE := a*b
TE := IE + c
'
fit_sem_IE_TE_tau <- sem(model_sem_IE_TE_tau, StressAtWork)
summary(fit_sem_IE_TE_tau, fit.measures = T, rsq = T)
# Modellvergleich
lavTestLRT(fit_sem_IE_TE_tau, fit_sem_IE_TE) # essentielle Tau-Äquivalenz wird verworfen!
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.