---
title: "Untitled"
author: "Kmg"
date: "1/7/2019"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#install.packages("nFactors")
#install.packages("lavaan")

```


```{r, echo = FALSE}
library(nFactors)
library(lavaan)

setwd("~/Dropbox/MIIVsem/Structural Search/R code/borkenau")

all <- list.files("~/Dropbox/MIIVsem/Structural Search/R code/borkenau", full.names = TRUE)

lowvar <- NULL
for (i in 1:length(all)){
  data <- read.csv(all[i], header = F)
  if (min(apply(data, 2, sd)) <.5) {
    lowvar[i] <- i
  }
}

# see if linear detrending needed
time <- seq(1:90)
data <- read.csv(all[15], header = F)
coeffs <- NULL
for (p in 1: 30){
  fita<- lm(data[,p] ~ time)  #run a simple regression 
  coeffs[p] <- fita$coefficients[2]
}
max(coeffs) # pretty low, not a problem

### lavaan ####

model1 <- '
FAC1 =~ NA*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC1 ~~ 1*FAC1'

model2 <- '
FAC1 =~ NA*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC2 =~ 0*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC1 ~~ 1*FAC1 + 0*FAC2
FAC2 ~~ 1*FAC2 
'

model3 <- '
FAC1 =~ NA*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC2 =~ 0*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC1 ~~ 1*FAC1 + 0*FAC2 + 0*FAC3
FAC2 ~~ 1*FAC2 + 0*FAC3
FAC3 =~ 0*V1 + 0*V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 
FAC3 ~~ 1*FAC3
'
model4 <- '
FAC1 =~ NA*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC2 =~ 0*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC1 ~~ 1*FAC1 + 0*FAC2 + 0*FAC3
FAC2 ~~ 1*FAC2 + 0*FAC3
FAC3 =~ 0*V1 + 0*V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 
FAC3 ~~ 1*FAC3

FAC4 =~ 0*V1 + 0*V2 + 0*V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 

FAC4 ~~ 1*FAC4 + 0*FAC1 + 0*FAC2 + 0*FAC3 
'

model5 <- '
FAC1 =~ NA*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC2 =~ NA*V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC1 ~~ 1*FAC1 + 0*FAC2 + 0*FAC3
FAC2 ~~ 1*FAC2 + 0*FAC3
FAC3 =~ NA*V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 
FAC3 ~~ 1*FAC3

FAC4 =~ 0*V1 + 0*V2 + 0*V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 

FAC4 ~~ 1*FAC4 + 0*FAC1 + 0*FAC2 + 0*FAC3 

FAC5 =~ 0*V1 + 0*V2 + 0*V3 + 0*V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 
FAC4 ~~ 0*FAC5
FAC5 ~~ 0*FAC1
FAC5 ~~ 0*FAC2
FAC5 ~~ 0*FAC3
FAC5 ~~ 1*FAC5'

#for (p in 1:length(all)){
p <- 6
bic_ind <- matrix(,5,1)
data <- read.csv(all[p], header = F)
fit1 <- lavaan::sem(model1, as.data.frame(data)) 
fit2 <- lavaan::sem(model2, as.data.frame(data)) 
fit3 <- lavaan::sem(model3, as.data.frame(data)) 
fit4 <- lavaan::sem(model4, as.data.frame(data)) 
fit5 <- lavaan::sem(model5, as.data.frame(data)) 

bic_ind[1] <- round(lavaan::fitMeasures(fit1)[c("bic")], digits = 2)
bic_ind[2] <- round(lavaan::fitMeasures(fit2)[c("bic")], digits = 2)
bic_ind[3] <- round(lavaan::fitMeasures(fit3)[c("bic")], digits = 2)
bic_ind[4] <- round(lavaan::fitMeasures(fit4)[c("bic")], digits = 2)
bic_ind[5] <- round(lavaan::fitMeasures(fit5)[c("bic")], digits = 2)
which(bic_ind == min(bic_ind)) #3


# lavaan::parameterEstimates(fit)
indices <-rbind(round(lavaan::fitMeasures(fit1)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr","bic")], digits = 2), round(lavaan::fitMeasures(fit2)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr","bic")], digits = 2), round(lavaan::fitMeasures(fit3)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr","bic")], digits = 2), round(lavaan::fitMeasures(fit4)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr","bic")], digits = 2), round(lavaan::fitMeasures(fit5)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr","bic")], digits = 2))
#}

knitr::kable(indices)


```
## Inspect results.

```{r}

summary(fit3)
```

## The loadings closest to zero for each factor are: 
```{r}
which(inspect(fit3, what = "std")$lambda[,1] == max(inspect(fit3,what="std")$lambda[,1]))

which(inspect(fit3, what = "std")$lambda[,2] == max(inspect(fit3,what="std")$lambda[,2]))

which(inspect(fit3, what = "std")$lambda[,3] == max(inspect(fit3,what="std")$lambda[,3]))
```

## Remove variable 18 from factors 1 and 3, remove variable 10 from factors 1 and 2. 
```{r}

model3r <- '
FAC1 =~ NA*V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC2 =~  NA*V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30
FAC3 =~  NA*V3 + V4 + V5 + V6 + V7 + V8 + V9 +V10 + V11 +V12 + V13 + V14 + V15 + V16 + V17 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 
FAC1 ~~ 1*FAC1
FAC2 ~~ 1*FAC2 + FAC1
FAC3 ~~ 1*FAC3 + FAC2 + FAC3
'
fit3r <- lavaan::sem(model3r, as.data.frame(data), std.lv=T, std.ov=TRUE) 
summary(fit3r)

```
## Fit final model by removing variables from factor equations if their loadings are below 0.50 to zero. 

```{r}
FA1 <- paste("FA1 =~", paste(names(which(abs(inspect(fit3r,what="std")$lambda[,1]) >= 0.3)), collapse = "+"))
FA2 <- paste("FA2 =~", paste(names(which(abs(inspect(fit3r,what="std")$lambda[,2]) >= 0.3)), collapse = "+"))
FA3 <- paste("FA3 =~", paste(names(which(abs(inspect(fit3r,what="std")$lambda[,3]) >= 0.3)), collapse = "+"))

model <- paste(FA1, FA2, FA3, 'FA1~~FA1 + FA2 + FA3', 'FA2~~FA2 + FA3', 'FA3~~FA3',sep = "\n")

fitfinal <- lavaan::sem(model, as.data.frame(data)) 

round(lavaan::fitMeasures(fitfinal)[c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr","bic")], digits = 2)

lambdas <- inspect(fitfinal, what = "std")$lambda
lambdas[which(lambdas == 0)] = NA
lambdas
```