workshop_part2.Rmd

---
title: "Workshop Part II"
author: "Fabio Votta"
date: "3 9 2018"
output: html_document
---

## Start

In diesem Tutorial beschäftigen uns mit den grundlegenden Komponenten von `lavaan`: die Modellsyntax, die Modellfunktionen (cfa, sem und growth) und die wichtigsten Extraktionsfunktionen (summary, parameterEstimates, fitmeasures).

## Modellsyntax

Im Zentrum des lavaan-Pakets steht die "Modellsyntax". Die Modellsyntax ist eine Beschreibung des zu schätzenden Modells.

In der R-Umgebung hat eine Regressionsformel die folgende Form:

```{r}
"y ~ x1 + x2 + x3 + x4"
```

Wollen wir mehrere Modelle spezifizieren, so können wir das auf folgende Art und Weise tun:

```{r}
"y1 ~ f1 + f2 + x1 + x2
f1 ~ f2 + f3
f2 ~ f3 + x1 + x2"
```

Wenn wir latente Variablen in einer der Regressionsformeln haben, müssen wir sie definieren, indem wir ihre (manifesten oder latenten) Indikatoren auflisten. Wir tun dies, indem wir den speziellen Operator "=~" verwenden.

```{r}
"f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6
f3 =~ y7 + y8 + y9 + y10"
```


Darüber hinaus werden Varianzen und Kovarianzen mit einem 'Double Tilde'-Operator angegeben, zum Beispiel so:

```{r}
"y1 ~~ y1 # Varianz
y1 ~~ y2 # Kovarianz
f1 ~~ f2 # Kovarianz"
```


Intercepts für beobachtete und latente Variablen sind einfache Regressionsformeln mit nur einem Schnittpunkt (explizit mit der Nummer '1' bezeichnet):


```{r}
"y1 ~ 1
f1 ~ 1"
```

Unter Verwendung dieser vier Formulierungstypen kann eine große Vielfalt latenter Variablenmodelle beschrieben werden. Die aktuelle Menge der Formel-Typen ist in der folgenden Tabelle zusammengefasst.

| formula type               | operator | mnemonic           |
|----------------------------|----------|--------------------|
| latent variable definition | =~       | is measured by     |
| regression                 | ~        | is regressed on    |
| (residual) (co)variance    | ~~       | is correlated with |
| intercept                  | ~ 1      | intercept          |


Eine vollständige Lavan-Modell-Syntax ist einfach eine Kombination dieser Formel-Typen, die zwischen einfachen Anführungszeichen stehen. Beispielsweise:

```{r}
myModel <- '# Regressionen
y1 + y2 ~ f1 + f2 + x1 + x2
f1 ~ f2 + f3
f2 ~ f3 + x1 + x2

# latente Variablendefinitionen
f1 = ~ y1 + y2 + y3
f2 = ~ y4 + y5 + y6
f3 = ~ y7 + y8 + y9 + y10

# Varianzen und Kovarianzen
y1 ~~ y1
y1 ~~ y2
f1 ~~ f2

# Intercepts
y1 ~ 1
f1 ~ 1
'
```

Der Codeteil oben erzeugt ein Modell-Syntax-Objekt namens myModel.

## CFA

Lasst uns zunächst einmal das Package `lavaan` einladen.

(Und das `tidyverse` zu haben ist immer gut)

```{r}
pacman::p_load(tidyverse, lavaan)
```


Wir beginnen mit einem einfachen Beispiel der konfirmatorischen Faktorenanalyse mit der Funktion `cfa()`. Das lavaan-Paket enthält einen integrierten Datensatz namens `HolzingerSwineford1939`. Rufen Sie die Hilfeseite für dieses Dataset auf:

```{r}
?HolzingerSwineford1939
```


Die Daten bestehen aus Mentalfähigkeitstest-Scores von Kindern der siebten und achten Klasse aus zwei verschiedenen Schulen (Pasteur und Grant-White). In unserer Version des Datasets sind nur 9 der ursprünglichen 26 Tests enthalten. Ein CFA-Modell, das häufig für diese 9 Variablen vorgeschlagen wird, besteht aus drei latenten Variablen (oder Faktoren) mit jeweils drei Indikatoren:

+ ein visueller Faktor, der durch 3 Variablen gemessen wird: x1, x2 und x3
+ ein Textfaktor, der durch 3 Variablen gemessen wird: x4, x5 und x6
+ ein Geschwindigkeitsfaktor, der durch 3 Variablen gemessen wird: x7, x8 und x9

Ein 3-Faktor-CFA-Beispiel Die entsprechende lavaan-Syntax zur Angabe dieses Modells lautet wie folgt:

```{r}
HS.model <- "visuell =~ x1 + x2 + x3
             textuell =~ x4 + x5 + x6
             Geschwindigkeit =~ x7 + x8 + x9"
```

In diesem Beispiel enthält die Modellsyntax nur drei "latente Variablendefinitionen". Jede Formel hat folgendes Format:

`latente Variable =~ Indikator1 + Indikator2 + Indikator3`

Wir können das Modell nun wie folgt berechnen:

```{r}
fit1 <- cfa(HS.model, data = HolzingerSwineford1939)
```


Die Funktion `cfa()` ist eine dedizierte Funktion zur Anpassung konfirmatorischer Faktorenanalyse-Modelle. Das erste Argument ist das benutzerdefinierte Modell. Das zweite Argument ist der Datensatz, der die beobachteten Variablen enthält. Sobald das Modell angepasst wurde, bietet die Funktion `summary()` eine schöne Zusammenfassung des angepassten Modells:


```{r}
summary(fit1, fit.measures = TRUE)
```

## SEM

Ein volles Strukturgleichungsmodell kann sehr ähnlich und sehr einfach geschätzt werden. Hier ein Beispiel:

```{r}
HS.model2 <- "visuell =~ x1 + x2 + x3
             textuell =~ x4 + x5 + x6
             Geschwindigkeit =~ x7 + x8 + x9
             Geschwindigkeit ~ visuell"
```

Wir benutzen nun einfach die `sem()` Funktion um die Regression zu schätzen.

```{r}
fit2 <- sem(HS.model2, data = HolzingerSwineford1939)
```

Lassen wir uns doch einfach mal die `summary()` ausgeben:

```{r}
summary(fit2, fit.measures = TRUE, standardized = T, rsquare = T)

```

Wir können uns auch direkt die `fitmeasures()` ausgeben lassen:

```{r}
fitmeasures(fit2, c("chisq", "df", "rmsea", "cfi", "srmr"))
```

Oder die Koeffizienten/Ladungen:

```{r}
parameterEstimates(fit2)
```

### Modiciationindices

Modiciationindices werden für jeden nicht freien (oder nichtredundanten) Parameter angezeigt. Die Modiciationindices werden um die Parameter der erwarteten Parameteränderung (`epc`) ergänzt. Die letzten drei Spalten enthalten die standardisierten EPC-Werte:

`sepc.lv`: nur Standardisierung der latenten Variablen; 
`sepc.all`: Standardisierung aller Variablen; 
`sepc.nox`: Standardisierung aller nicht-exogenen manifesten Variablen

```{r}
modificationindices(fit2, sort. = T)
```
### Moderatoren

Moderatoren zwischen manifesten Variablen zu erstellen ist ziemlich einfach. Wir können einfach eine Variable mal die andere Variable nehmen und diese als neue Variable abspeichern. Das sieht ungefähr so aus:

```{r}
HolzingerSwineford1939 %>% 
  mutate(x1_x2 = x1 * x2)
```

Wollen wir allerdings eine Interaktion zwischen zwei latenten Variablen berechnen, müssen wir etwas komplexer vorgehen. Dazu brauchen wir die `indProd` Funktion aus dem `semTools` package.

```{r}
pacman::p_load(semTools)

HolzingerSwineford1939 <- indProd(HolzingerSwineford1939, 
              var1 = c("x1", "x2", "x3"), 
              var2 = c("x4", "x5", "x6"))

HS.model3 <- "visuell =~ x1 + x2 + x3
             textuell =~ x4 + x5 + x6
             Geschwindigkeit =~ x7 + x8 + x9

             # interaction latent factor 
             vis.text =~ x1.x4 + x2.x5 + x3.x6

             Geschwindigkeit ~ visuell + textuell + vis.text"

fit3 <- sem(HS.model3, data = HolzingerSwineford1939)

summary(fit3)
```

Mehr zu Moderatoren:

+ https://rpubs.com/mkearney/103098
+ https://www.rdocumentation.org/packages/semTools/versions/0.5-0/topics/indProd

### Nested Models

Wenn wir ein sogenanntes *nested* Modell rechnen dann können wir mit Hilfe der `anova()` Funktion den Modellfit testen.

```{r}
anova(fit2, fit3)
```

**Ein signifikanter Chi Square Test zeigt, dass `fit3` besser fitted als `fit2`.

```{r}
fitmeasures(fit2, c("chisq", "df", "rmsea", "cfi", "srmr"))
fitmeasures(fit3, c("chisq", "df", "rmsea", "cfi", "srmr"))
```

Wenn wir die Fitmeasures vergleichen kommen wir zur etwa gleichen Konklusion wie zuvor.

### Mediation

Das Beispiel veranschaulicht die Verwendung des Operators `:=` in der Lavan-Modell-Syntax. Dieser Operator definiert neue Parameter, die eine beliebige Funktion der ursprünglichen Modellparameter sind. Die Funktion muss jedoch in Bezug auf die Parameterbezeichnungen angegeben werden, die explizit in der Modellsyntax erwähnt werden. Standardmäßig werden die Standardfehler für diese definierten Parameter mit der so genannten Delta-Methode berechnet.

```{r}
HS.model4 <- "
             # Faktor
             Geschwindigkeit =~ x7 + x8 + x9

             # direkter Effekt
             Geschwindigkeit ~ c*x1

             # mediator
             x2 ~ a*x1
             Geschwindigkeit ~ b*x2

            # indirekter Effect (a*b)
             ab := a*b

            # total effect
             total := c + (a*b)"

fit4 <- sem(HS.model4, data = HolzingerSwineford1939)
summary(fit4)
```



## Growth

### Data

```{r}
lower <- '
 0.63 
 0.50 0.60 
 0.48 0.48 0.58
 0.47 0.48 0.51 0.67 
-0.06 -0.07 -0.08 -0.09 0.24
 1.87 1.80 1.64 1.57 −0.04 64.04 
 1.16 1.11 1.04 0.93 −0.41 42.59 43.23 '

cov_dat <- getCov(lower,  names = c("Time1", "Time2", "Time3", "Time4",
                                    "State", "Poverty", "StPov"))

cov_dat

sample_means <- c(5.17, 5.32, 5.40, 5.52, 0.62, -1.67, -1.08)



```

### Schritt für Schritt

#### 1. Latenter Intercept

1. Modell für einen mittleren latenten Intercept und
eingeschränkte Residualvarianzen.

```{r}
crime_model1 <- 
  '# intercept
   i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4
   i ~~ 0*i
   # Residualvarianzen
   Time1 ~~ r*Time1
   Time2 ~~ r*Time2
   Time3 ~~ r*Time3
   Time4 ~~ r*Time4'

crime_fit1 <- growth(crime_model1, 
                     sample.cov = cov_dat,
                     sample.mean = sample_means,
                     sample.nobs = 952)

summary(crime_fit1)
```

#### 2. Var. Latenter Intercept

2. Modell für einen mittleren latenten Intercept, der variieren darf und Residualvarianzen beschränkt.

```{r}
crime_model2 <- 
  '# intercept
   i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4
   # Residualvarianzen
   Time1 ~~ r*Time1
   Time2 ~~ r*Time2
   Time3 ~~ r*Time3
   Time4 ~~ r*Time4'

crime_fit2 <- growth(crime_model2, 
                     sample.cov = cov_dat,
                     sample.mean = sample_means,
                     sample.nobs = 952)

summary(crime_fit2)
```

##### Modellvergleich

```{r}
anova(crime_fit1, crime_fit2)
```

Fitmeasures Vergleich

 + Little Helperfunction

```{r}
fit_m <- function(x) {
  fitmeasures(x, c("chisq", "df", "rmsea", "cfi", "srmr", "AIC"))
}

fit_compare <- function(...) {
  cbind(...) %>% 
  data.frame() %>% 
  rownames_to_column("fitmeasure") 
}
```


```{r}
fit_compare(fit_m(crime_fit1), 
            fit_m(crime_fit2))
```

#### 3. Latenter Slope

Modell mit einem mittleren latenten Intercept, der variieren darf, latente Steigung und eingeschränkten Varianzen

```{r}
crime_model3 <- 
  '# intercept
   i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4
   # steigung (slope)
   s =~ 0*Time1 + 1*Time2 + 2*Time3 + 3*Time4
   s ~ 0*1
   s ~~ 0*i
   # Residualvarianzen
   Time1 ~~ r*Time1
   Time2 ~~ r*Time2
   Time3 ~~ r*Time3
   Time4 ~~ r*Time4'

crime_fit3 <- growth(crime_model3, 
                     sample.cov = cov_dat,
                     sample.mean = sample_means,
                     sample.nobs = 952)

summary(crime_fit3)
```

##### Modellvergleich

```{r}
anova(crime_fit1, crime_fit2, crime_fit3)
```


```{r}
fit_compare(fit_m(crime_fit1), 
            fit_m(crime_fit2), 
            fit_m(crime_fit3))
```


#### 4. Var. Latenter Slope

Modell mit einem mittleren latenten Intercept, der variieren darf, variiender latenter Steigung und eingeschränkten Varianzen

```{r}
crime_model4 <- 
  '# intercept
   i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4
   # steigung (slope)
   s =~ 0*Time1 + 1*Time2 + 2*Time3 + 3*Time4
   # Residualvarianzen
   Time1 ~~ r*Time1
   Time2 ~~ r*Time2
   Time3 ~~ r*Time3
   Time4 ~~ r*Time4'

crime_fit4 <- growth(crime_model4, 
                     sample.cov = cov_dat,
                     sample.mean = sample_means,
                     sample.nobs = 952)

summary(crime_fit4)
```

##### Modellvergleich

```{r}
anova(crime_fit1, crime_fit2, crime_fit3, crime_fit4)
```


```{r}
fit_compare(fit_m(crime_fit1), 
            fit_m(crime_fit2), 
            fit_m(crime_fit3), 
            fit_m(crime_fit4))
```


#### 5. Uneingeschränkt

```{r}
crime_model5 <- 
  '# intercept
   i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4
   # steigung (slope)
   s =~ 0*Time1 + 1*Time2 + 2*Time3 + 3*Time4'

crime_fit5 <- growth(crime_model5, 
                     sample.cov = cov_dat,
                     sample.mean = sample_means,
                     sample.nobs = 952)

summary(crime_fit5)
```


##### Modellvergleich

```{r}
anova(crime_fit1, crime_fit2, 
      crime_fit3, crime_fit4, 
      crime_fit5)
```


```{r}
fit_compare(fit_m(crime_fit1), 
            fit_m(crime_fit2), 
            fit_m(crime_fit3), 
            fit_m(crime_fit4), 
            fit_m(crime_fit5))
```


#### 6. Unabhängige Variablen

```{r}
crime_model6 <- 
  '# intercept
   i =~ 1*Time1 + 1*Time2 + 1*Time3 + 1*Time4
   # steigung (slope)
   s =~ 0*Time1 + 1*Time2 + 2*Time3 + 3*Time4
   # regression
   i + s ~ State + Poverty + StPov'

crime_fit6 <- growth(crime_model6, 
                     sample.cov = cov_dat,
                     sample.mean = sample_means,
                     sample.nobs = 952)

summary(crime_fit6)
```

##### Modellvergleich

```{r}
anova(crime_fit1, crime_fit2, 
      crime_fit3, crime_fit4, 
      crime_fit5, crime_fit6)
```


```{r}
fit_compare(fit_m(crime_fit1), 
            fit_m(crime_fit2), 
            fit_m(crime_fit3), 
            fit_m(crime_fit4), 
            fit_m(crime_fit5), 
            fit_m(crime_fit6))
```