---
title: "L5 PFA SFA lavaan"
author: "Stephanie Lane & Kathleen M. Gates"
output:
  pdf_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE)
library(lavaan)
```

## Generating Data

```{r, echo = TRUE}
set.seed(12345678)
library(mvtnorm)
# setting up matrices

npad      <- 100 #Occasions to throw out to wash away the effects of initial condition
time      <- 200 + npad 
ne        <- 2 #Number of latent variables
ny        <- 6 #Number of manifest variables
psi       <- matrix(c(2.77, 2.47,  # Residual variance-covariance matrix
                      2.47, 8.40),  
                    ncol = ne, byrow = T)
lambda    <- matrix(c(1, 0,        # Lambda matrix containing contemporaneous relations among  
                      2, 0,        # observed variables and 2 latent variables. 
                      1, 0, 
                      0, 1, 
                      0, 2, 
                      0, 1), 
                    ncol = ne, byrow = TRUE)
theta     <- diag(.5, ncol = ny, nrow = ny)  # Measurement error variances. 
epsilon   <- rmvnorm(time, mean = c(0, 0, 0, 0, 0, 0), sigma = theta)  # Measurement errors.
beta      <- matrix(c(0.5, 0,    # Lagged directed relations among variables. 
                      0.4, 0.5), 
                    ncol = ne, byrow = TRUE)
zeta      <- mvtnorm::rmvnorm(time+npad, mean = c(0, 0), sigma = psi)   # Latent variable residuals.
etaC      <- matrix(0, nrow = ne, ncol = time + npad) # Set up matrix for contemporaneous variables. 
etaL      <- matrix(0, nrow = ne, ncol = time + npad + 1) # Set up matrix for lagged variables. 

etaL[,1]  <- c(0,0)
# generate factors
for (i in 1:(time+npad)){
  etaC[ ,i]   <- beta %*% etaL[ ,i] + zeta[i, ]
  etaL[ ,i+1] <- etaC[,i]
}
etaC <- etaC[,(npad+1):(npad+time)]
etaL <- etaL[,(npad+1):(npad+time)]
eta <- t(etaC)

# generate observed series
y   <- matrix(0, nrow = time, ncol = ny)
for (p in 1:nrow(y)){
  y[p, ] <- lambda %*% eta[p, ] + epsilon[p, ]
}

y <- y[101:time,]
```

\newpage 

##Process Factor Analysis: lavaan
Prepare data
```{r, echo = TRUE, comment = ''}
# create block toeplitz data structure
y1   <- y[1:(nrow(y)-1), ]
y0   <- y[2:nrow(y), ]
y_bt <- data.frame(y1, y0)
names(y_bt) <- c(paste0("y", seq(1:6), "L"), paste0("y", seq(1:6), "C"))

```
\newpage
CM5 (no cross-lag)

```{r, echo = TRUE}
model5 <- '
eta1C =~ eq1*y1C + eq2*y2C + eq3*y3C 
eta2C =~ eq4*y4C + eq5*y5C + eq6*y6C
eta1L =~ eq1*y1L + eq2*y2L + eq3*y3L
eta2L =~ eq4*y4L + eq5*y5L + eq6*y6L

eta1C ~~ eta2C
eta1L ~~ eta2L

eta1C ~ eta1L
eta2C ~ eta2L
'


fit  <- sem(model5, y_bt)
# summary(fit) # See full results from SEM.
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)
```

\newpage

CM 6 (incorrect cross-lag)

```{r, echo= TRUE}
model6 <- '
eta1C =~ eq1*y1C + eq2*y2C + eq3*y3C 
eta2C =~ eq4*y4C + eq5*y5C + eq6*y6C
eta1L =~ eq1*y1L + eq2*y2L + eq3*y3L
eta2L =~ eq4*y4L + eq5*y5L + eq6*y6L

eta1C ~~ eta2C
eta1L ~~ eta2L

eta1C ~ eta1L
eta2C ~ eta2L
eta1C ~ eta2L
'


fit  <- sem(model6, y_bt)
# summary(fit) # See full results from SEM.
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)
```
\newpage

Correct model (CM7)
```{r, echo = TRUE, comment = ''}
# create block toeplitz data structure
yL   <- y[1:(nrow(y)-1), ]
yC   <- y[2:nrow(y), ]
y_bt <- data.frame(yL, yC)
names(y_bt) <- c(paste0("y", seq(1:6), "L"), paste0("y", seq(1:6), "C"))
library(lavaan)
model <- '
eta1C =~ eq1*y1C + eq2*y2C + eq3*y3C 
eta2C =~ eq4*y4C + eq5*y5C + eq6*y6C

eta1C ~~ eta2C

'


fit  <- sem(model, y_bt)
# summary(fit) # See full results from SEM.
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)
```
\newpage
##SFA

# Lag one models
Prepare data

```{r, echo = TRUE, comment = ''}
# create block toeplitz data structure
y1   <- y[1:(nrow(y)-1), ]
y0   <- y[2:nrow(y), ]
y_bt <- data.frame(y1, y0)
names(y_bt) <- c(paste0("y", seq(1:6), "_1"), paste0("y", seq(1:6), "_0"))

```
CM1: 1 LV, 1 lag
```{r}
library(lavaan)
model1 <- '

eta11 =~ 1*y1_1 + eq1*y2_1 + eq2*y3_1  + eq3*y4_1 + eq4*y5_1 + eq5*y6_1 + y1_0 + y2_0 + y3_0  + y4_0 + y5_0 + y6_0

eta10 =~ 1*y1_0 + eq1*y2_0 + eq2*y3_0 + eq3*y4_0 + eq4*y5_0 + eq5*y6_0

eta11 ~~  0*eta10 

'

fit  <- sem(model1, y_bt)
# summary(fit)
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)

```
\newpage
CM 3: 2 LVs, lag 1
```{r, echo = TRUE}
model3 <- '

eta11 =~ 1*y1_1 + eq1*y2_1 + eq2*y3_1  +eq3*y1_0 + eq4*y2_0 + eq5*y3_0
eta10 =~ 1*y1_0 + eq1*y2_0 + eq2*y3_0

eta21 =~ 1*y4_1 + eq6*y5_1 + eq7*y6_1 + eq8*y4_0 + eq9*y5_0 + eq10*y6_0
eta20 =~ 1*y4_0 + eq6*y5_0 + eq7*y6_0


eta11 ~~ eta21 + 0*eta10 + 0*eta20
eta10 ~~ eta20 + 0*eta21 

eta21 ~~ 0*eta20
'

fit  <- sem(model3, y_bt)
# summary(fit)
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)

```

\newpage
# 2 lags
Prep data
```{r, echo = TRUE, comment = ''}
# create block toeplitz data structure
y2   <- y[1:(nrow(y)-2), ]
y1   <- y[2:(nrow(y)-1), ]
y0   <- y[3:nrow(y), ]
y_bt <- data.frame(y2, y1, y0)
names(y_bt) <- c(paste0("y", seq(1:6), "_2"), paste0("y", seq(1:6), "_1"), paste0("y", seq(1:6), "_0"))
```
CM 2: 1 LV, 2 lags
```{r}
model <- '

eta12 =~ 1*y1_2 + eq1*y2_2 + eq2*y3_2 + eq3*y4_2 + eq4*y5_2 +
eq5*y6_2 + eq6*y1_1 + eq7*y2_1 + eq8*y3_1  + y1_0 + y2_0 + y3_0 + eq10*y4_1 + eq11*y5_1 + eq12*y6_1 + y4_0 + y5_0 + y6_0

eta11 =~ 1*y1_1 + eq1*y2_1 + eq2*y3_1 + eq3*y4_1 + eq4*y5_1 +
eq5*y6_1 + eq6*y1_0 + eq7*y2_0 + eq8*y3_0 +  eq10*y4_0 + eq11*y5_0 + eq12*y6_0 

eta10 =~ 1*y1_0 + eq1*y2_0 + eq2*y3_0 + eq3*y4_0 + eq4*y5_0 + eq5*y6_0

eta12 ~~ 0*eta11 + 0*eta10 
eta11 ~~ 0*eta10 


'

fit  <- sem(model, y_bt)
# summary(fit)
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)

```

\newpage
CM4: 2 LVs, 2 lags

```{r, echo = TRUE, comment = ''}
# create block toeplitz data structure

model4 <- '
eta12 =~ 1*y1_2 + eq1*y2_2 + eq2*y3_2 + eq3*y1_1 + eq4*y2_1 + eq5*y3_1 + y1_0 + y2_0 + y3_0 
eta11 =~ 1*y1_1 + eq1*y2_1 + eq2*y3_1 + eq3*y1_0 + eq4*y2_0 + eq5*y3_0
eta10 =~ 1*y1_0 + eq1*y2_0 + eq2*y3_0  


eta22 =~ 1*y4_2 + eq6*y5_2 + eq7*y6_2 + eq8*y4_1 + eq9*y5_1 + eq10*y6_1 + y4_0 + y5_0 + y6_0  
eta21 =~ 1*y4_1 + eq6*y5_1 + eq7*y6_1 + eq8*y4_0 + eq9*y5_0 + eq10*y6_0
eta20 =~ 1*y4_0 + eq6*y5_0 + eq7*y6_0  


eta12 ~~ eta22 + 0*eta11 + 0*eta10 + 0*eta21 + 0*eta20 
eta11 ~~ eta21 + 0*eta10 + 0*eta20 + 0*eta22 
eta20 ~~ 0*eta22
eta21 ~~ 0*eta20 + eta22
eta10 ~~ eta20 + 0*eta21 + 0*eta22
'

fit  <- sem(model4, y_bt)
summary(fit)
lava <- parameterEstimates(fit)[ ,c(1, 2, 3, 5, 6,8)]
lava <- lava[lava$op %in% c("~", "=~"), ]
lava <- lava[,-2]
lava[,3:4] <- round(lava[,3:4], 3)
rownames(lava) <- NULL
knitr::kable(lava)
round(fitMeasures(fit)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "bic")], digits = 3)

```