library(lavaan)
datcon <- read.table("lecture11consci.csv", sep=",", header=TRUE, na.strings="999")
mcon <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
caut =~ c2 + c8 + c14 + c20 + c26 + c32 + c38 + c44 + c50 + c56
duti =~ c3 + c9 + c15 + c21 + c27 + c33 + c39 + c45 + c51 + c57
orde =~ c4 + c10 + c16 + c22 + c28 + c34 + c40 + c46 + c52 + c58
disc =~ c5 + c11 + c17 + c23 + c29 + c35 + c41 + c47 + c53 + c59
effi =~ c6 + c12 + c18 + c24 + c30 + c36 + c42 + c48 + c54 + c60'
outcon <- cfa(mcon, data=datcon)
summary(outcon, fit=TRUE, std=TRUE)
inspect(outcon, "cov.lv")
inspect(outcon, "cor.lv")

mconh <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
caut =~ c2 + c8 + c14 + c20 + c26 + c32 + c38 + c44 + c50 + c56
duti =~ c3 + c9 + c15 + c21 + c27 + c33 + c39 + c45 + c51 + c57
orde =~ c4 + c10 + c16 + c22 + c28 + c34 + c40 + c46 + c52 + c58
disc =~ c5 + c11 + c17 + c23 + c29 + c35 + c41 + c47 + c53 + c59
effi =~ c6 + c12 + c18 + c24 + c30 + c36 + c42 + c48 + c54 + c60
con =~ achi + caut + duti + orde + disc + effi'
outconh <- cfa(mconh, data=datcon)
summary(outconh, fit=TRUE, std=TRUE)

library(nonnest2)
vuongtest(outcon, outconh, nested=TRUE)

library(readxl)
detailscon <- read_xlsx("lecture11conscidetails.xlsx")
detailscon <- data.frame(detailscon[,c("VariableName", "Direction")])
detailscon <- detailscon[detailscon$VariableName %in% paste0("c", 1:60),]

for(i in 1:nrow(detailscon)) {
  if(detailscon[i, "Direction"] == "-"){
    temp <- datcon[,detailscon[i, "VariableName"]]
    temp <- 6 - temp
    datcon[,detailscon[i, "VariableName"]] <- temp
  }
}
outcon <- cfa(mcon, data=datcon)
outconh <- cfa(mconh, data=datcon)

library(semTools)
compRelSEM(outconh, higher="con")


mconb <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
caut =~ c2 + c8 + c14 + c20 + c26 + c32 + c38 + c44 + c50 + c56
duti =~ c3 + c9 + c15 + c21 + c27 + c33 + c39 + c45 + c51 + c57
orde =~ c4 + c10 + c16 + c22 + c28 + c34 + c40 + c46 + c52 + c58
disc =~ c5 + c11 + c17 + c23 + c29 + c35 + c41 + c47 + c53 + c59
effi =~ c6 + c12 + c18 + c24 + c30 + c36 + c42 + c48 + c54 + c60
con =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55 
+ c2 + c8 + c14 + c20 + c26 + c32 + c38 + c44 + c50 + c56 
+ c3 + c9 + c15 + c21 + c27 + c33 + c39 + c45 + c51 + c57
+ c4 + c10 + c16 + c22 + c28 + c34 + c40 + c46 + c52 + c58
+ c5 + c11 + c17 + c23 + c29 + c35 + c41 + c47 + c53 + c59
+ c6 + c12 + c18 + c24 + c30 + c36 + c42 + c48 + c54 + c60
con ~~ 0*achi + 0*caut + 0*duti + 0*orde + 0*disc + 0*effi
achi ~~ 0*caut + 0*duti + 0*orde + 0*disc + 0*effi
caut ~~ 0*duti + 0*orde + 0*disc + 0*effi
duti ~~ 0*orde + 0*disc + 0*effi
orde ~~ 0*disc + 0*effi
disc ~~ 0*effi
'
outconb <- cfa(mconb, data=datcon, std.lv=TRUE)
summary(outconb, fit=TRUE, std=TRUE)

library(nonnest2)
vuongtest(outcon, outconb)
vuongtest(outconb, outconh, nested=TRUE)

compRelSEM(outconb)

###################################

x <- seq(-2.5, 2.5, by=0.01)
y <- dnorm(x)
plot(x, y, type="l", xlab="ystar", ylab="Density")

lines(c(0.25, 0.25), c(0, dnorm(0.25)))

x <- seq(-2.5, 2.5, by=0.01)
y <- dnorm(x)
plot(x, y, type="l", xlab="ystar", ylab="Density")

lines(c(-1.25, -1.25), c(0, dnorm(-1.25)))
lines(c(-0.15, -0.15), c(0, dnorm(-0.15)))
lines(c(0.9, 0.9), c(0, dnorm(0.9)))

set.seed(123321)
f <- rnorm(10, 0, 1)
e <- rnorm(10, 0, 0.6)
ystar <- 0.8*f + e
ycut <- cut(ystar, c(-Inf, 0.25, Inf))
y <- as.numeric(ycut) - 1
data.frame(f, e, ystar, ycut, y)


set.seed(113113)
f <- rnorm(10, 0, 1)
e <- rnorm(10, 0, 0.6)
ystar <- 0.8*f + e
ycut <- cut(ystar, c(-Inf, -1.25, -0.15, 0.90, Inf))
y <- as.numeric(ycut) - 1
data.frame(f, e, ystar, ycut, y)


tau <- 0.25
lambda <- 0.8
f <- seq(from=-3, to=3, by=0.01)
prob1 <- pnorm(-tau + lambda*f)
prob0 <- pnorm(tau - lambda*f)
plot(f, prob1, type="l", col="blue", xlab="F", ylab="Probability")
lines(f, prob0, col="red")
legend(-3, 0.8, c("c=0","c=1"), col=c("red", "blue"), lty=1)

tau <- c(-1.25, -0.15, 0.90)
lambda <- 0.8
f <- seq(from=-3, to=3, by=0.01)
prob1 <- pnorm(tau[1] - lambda*f)
prob2 <- pnorm(tau[2] - lambda*f) - pnorm(tau[1] - lambda*f)
prob3 <- pnorm(tau[3] - lambda*f) - pnorm(tau[2] - lambda*f)
prob4 <- 1 - pnorm(tau[3] - lambda*f)
plot(f, prob1, type="l", col="red", xlab="F", ylab="Probability")
lines(f, prob2, col="blue")
lines(f, prob3, col="purple")
lines(f, prob4, col="darkgreen")
legend(-1.2, 0.9, c("c=1","c=2", "c=3", "c=4"), col=c("red", "blue", "purple", "darkgreen"), lty=1)

library(lavaan)
datdicho <- read.table("lecture12ex1.csv", sep=",", header=TRUE)

mdicho <- '
f1 =~ y1 + y2 + y3 + y4
f2 =~ y5 + y6 + y7 + y8
'
outdicho <- cfa(mdicho, data=datdicho, ordered=TRUE)
summary(outdicho, std=TRUE, fit=TRUE)
modindices(outdicho)

outdicho2 <- cfa(mdicho, data=datdicho, ordered=TRUE, estimator="WLSMV")
outdicho3 <- cfa(mdicho, data=datdicho, ordered=TRUE, estimator="DWLS",
                 se="standard", test="scaled.shifted")
outdicho3

outdichotheta <- cfa(mdicho, data=datdicho, ordered=TRUE, 
                     parameterization="theta")
summary(outdichotheta, fit=TRUE, std=TRUE)

library(lavaan)
datlikert <- read.table("lecture12ex2.csv", sep=",", header=TRUE)
head(datlikert)


mlikert <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6
f3 =~ y7 + y8 + y9
f4 =~ y10 + y11 + y12
'
outlikert <- cfa(mlikert, data=datlikert, ordered=TRUE)
summary(outlikert, fit=TRUE, std=TRUE)


mlikert2 <- '
f1 =~ y1 + lambda2*y2 + lambda3*y3
f2 =~ y4 + lambda2*y5 + lambda3*y6
f3 =~ y7 + lambda2*y8 + lambda3*y9
f4 =~ y10 + lambda2*y11 + lambda3*y12
'
outlikert2 <- cfa(mlikert2, data=datlikert, ordered=TRUE)
summary(outlikert2, fit=TRUE, std=TRUE)
anova(outlikert, outlikert2)

mlikert3 <- '
f1 =~ y1 + lambda2*y2 + lambda3*y3
f2 =~ y4 + lambda2*y5 + lambda3*y6
f3 =~ y7 + lambda2*y8 + lambda3*y9
f4 =~ y10 + lambda2*y11 + lambda3*y12
f2 ~ f1
f3 ~ f2
f4 ~ f3
'
outlikert3 <- cfa(mlikert3, data=datlikert, ordered=TRUE)
summary(outlikert3, fit=TRUE, std=TRUE)
anova(outlikert2, outlikert3)

mlikert4 <- '
f1 =~ y1 + lambda2*y2 + lambda3*y3
f2 =~ y4 + lambda2*y5 + lambda3*y6
f3 =~ y7 + lambda2*y8 + lambda3*y9
f4 =~ y10 + lambda2*y11 + lambda3*y12
f2 ~ beta*f1
f3 ~ beta*f2
f4 ~ beta*f3
f2 ~~ psi*f2
f3 ~~ psi*f3
f4 ~~ psi*f4
'
outlikert4 <- cfa(mlikert4, data=datlikert, ordered=TRUE)
summary(outlikert4, fit=TRUE, std=TRUE)
anova(outlikert3, outlikert4)

set.seed(121211)
n <- 10
x1 <- rnorm(n, 0, 1)
x2 <- 0.545*x1 + rnorm(n, 0, sqrt(0.338))
x3 <- 0.545*x2 + rnorm(n, 0, sqrt(0.338))
x4 <- 0.545*x3 + rnorm(n, 0, sqrt(0.338))
X <- cbind(x1, x2, x3, x4)

plot(1:4, seq(min(X), max(X), length.out=4), type="n", xaxt='n', xlab="Time", ylab="f")
axis(side=1, at=1:4)
for(i in 1:n) {
  lines(1:4, X[i,], col=rainbow(n)[i])
}

library(lavaan)
library(psych)
datcon <- read.table("lecture11consci.csv", sep=",", header=TRUE, na.strings="999")
datcon2 <- datcon[,paste0("c", 1:60)]
fout7 <- fa(datcon2, 7, rotate="bifactor")
fout7
loadings(fout7)

negcon <- c(20, 26, 32:35, 38:42, 44:60)

for(i in negcon) {
  datcon[,paste0("c", i)] <- 6 - datcon[,paste0("c", i)]
}

datcon$achi <- with(datcon, c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55)
datcon$caut <- with(datcon, c2 + c8 + c14 + c20 + c26 + c32 + c38 + c44 + c50 + c56)
datcon$duti <- with(datcon, c3 + c9 + c15 + c21 + c27 + c33 + c39 + c45 + c51 + c57)
datcon$orde <- with(datcon, c4 + c10 + c16 + c22 + c28 + c34 + c40 + c46 + c52 + c58)
datcon$disc <- with(datcon, c5 + c11 + c17 + c23 + c29 + c35 + c41 + c47 + c53 + c59)
datcon$effi <- with(datcon, c6 + c12 + c18 + c24 + c30 + c36 + c42 + c48 + c54 + c60)

mcon <- '
con =~ achi + caut + duti + orde + disc + effi
'
outcon <- cfa(mcon, data=datcon)
summary(outcon, fit=TRUE, std=TRUE)
modindices(outcon)

mcon2 <- '
con =~ achi + caut + duti + orde + disc + effi
achi ~~ effi
'
outcon2 <- cfa(mcon2, data=datcon)
summary(outcon2, fit=TRUE, std=TRUE)
modindices(outcon2)

mcon3 <- '
con =~ achi + caut + duti + orde + disc + effi
achi ~~ effi
caut ~~ orde
'
outcon3 <- cfa(mcon3, data=datcon)
summary(outcon3, fit=TRUE, std=TRUE)
modindices(outcon3)

mcon4 <- '
con =~ achi + caut + duti + orde + disc + effi
achi ~~ effi
caut ~~ orde + duti
'
outcon4 <- cfa(mcon4, data=datcon)
summary(outcon4, fit=TRUE, std=TRUE)

datcon$male <- datcon$Gender == 1
datcon$ugender <- datcon$Gender == 3
datcon$agec <- datcon$Age - 25

mcon5 <- '
con =~ achi + caut + duti + orde + disc + effi
achi ~~ effi
caut ~~ orde + duti
con ~ male + ugender + agec
'
outcon5 <- sem(mcon5, data=datcon)
summary(outcon5, fit=TRUE, std=TRUE, rsquare=TRUE)


datagg <- read.table("lecture12aggression.csv", sep=",", header=TRUE, na.strings="99")
datagg <- datagg[datagg$educationtype == 1,]
head(datagg)
# A : Mindfulness: balancing

library(lavaan)
library(psych)
outmind1 <- fa(datagg[,paste0("A", 1:21)], 1)
outmind2 <- fa(datagg[,paste0("A", 1:21)], 2)
outmind3 <- fa(datagg[,paste0("A", 1:21)], 3)
outmind4 <- fa(datagg[,paste0("A", 1:21)], 4)
outmind5 <- fa(datagg[,paste0("A", 1:21)], 5)
rbind(do.call(c, outmind1[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outmind2[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outmind3[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outmind4[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outmind5[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]))
fa.parallel(datagg[,paste0("A", 1:21)])
outmind4 <- fa(datagg[,paste0("A", 1:21)], 4, rotate="bifactor")
loadings(outmind4)

outmind4a <- fa(datagg[,paste0("A", 1:21)], 4, rotate="quartimin")
loadings(outmind4a)
outmind4a$Phi

datagg$AP1 <- apply(datagg[,paste0("A", 13:21)], 1, mean)
datagg$AP2 <- apply(datagg[,paste0("A", 1:3)], 1, mean)
datagg$AP3 <- apply(datagg[,paste0("A", 4:7)], 1, mean)
datagg$AP4 <- apply(datagg[,paste0("A", 8:12)], 1, mean)

outfriend1 <- fa(datagg[,paste0("B", 1:25)], 1)
outfriend2 <- fa(datagg[,paste0("B", 1:25)], 2)
outfriend3 <- fa(datagg[,paste0("B", 1:25)], 3)
outfriend4 <- fa(datagg[,paste0("B", 1:25)], 4)
outfriend5 <- fa(datagg[,paste0("B", 1:25)], 5)
rbind(do.call(c, outfriend1[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriend2[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriend3[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriend4[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriend5[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]))
fa.parallel(datagg[,paste0("B", 1:25)])
outfriend4 <- fa(datagg[,paste0("B", 1:25)], 4, rotate="bifactor")
loadings(outfriend4)

friendB <- c(4, 5, 10, 11, 18, 22, 23)
friendA <- setdiff(1:25, friendB)

outfriendA1 <- fa(datagg[,paste0("B", friendA)], 1)
outfriendA2 <- fa(datagg[,paste0("B", friendA)], 2)
outfriendA3 <- fa(datagg[,paste0("B", friendA)], 3)
outfriendA4 <- fa(datagg[,paste0("B", friendA)], 4)
outfriendA5 <- fa(datagg[,paste0("B", friendA)], 5)
rbind(do.call(c, outfriendA1[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendA2[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendA3[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendA4[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendA5[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]))
fa.parallel(datagg[,paste0("B", friendA)])
outfriendA3 <- fa(datagg[,paste0("B", friendA)], 3, rotate="bifactor")
loadings(outfriendA3)
outfriendA3a <- fa(datagg[,paste0("B", friendA)], 3, rotate="quartimin")
loadings(outfriendA3a)
outfriendA3a$Phi
outfriendA2a <- fa(datagg[,paste0("B", friendA)], 2, rotate="quartimin")
loadings(outfriendA2a)
outfriendA2a$Phi
loadings(outfriendA1)
sort(loadings(outfriendA1)[,1])

datagg$BPA1 <- apply(datagg[,paste0("B", c(1, 24, 25, 20, 12, 17))], 1, mean)
datagg$BPA2 <- apply(datagg[,paste0("B", c(9, 3, 6, 15, 14, 21))], 1, mean)
datagg$BPA3 <- apply(datagg[,paste0("B", c(2, 7, 13, 8, 16, 19))], 1, mean)

# Two factor is not that great so keep three factors

outfriendB1 <- fa(datagg[,paste0("B", friendB)], 1)
outfriendB2 <- fa(datagg[,paste0("B", friendB)], 2)
outfriendB3 <- fa(datagg[,paste0("B", friendB)], 3)
outfriendB4 <- fa(datagg[,paste0("B", friendB)], 4)
rbind(do.call(c, outfriendB1[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendB2[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendB3[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outfriendB4[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]))
fa.parallel(datagg[,paste0("B", friendB)])
outfriendB2 <- fa(datagg[,paste0("B", friendB)], 2, rotate="quartimin")
loadings(outfriendB2)
outfriendB2$Phi
sort(loadings(outfriendB1)[,1])

datagg$BPB1 <- apply(datagg[,paste0("B", c(4, 5, 11))], 1, mean)
datagg$BPB2 <- apply(datagg[,paste0("B", c(10, 18))], 1, mean)
datagg$BPB3 <- apply(datagg[,paste0("B", c(22, 23))], 1, mean)



# C : School Connectedness
outschool1 <- fa(datagg[,paste0("C", 1:5)], 1)
summary(outschool1)
loadings(outschool1)

datagg$CP1 <- apply(datagg[,paste0("C", c(1, 3))], 1, mean)
datagg$CP2 <- apply(datagg[,paste0("C", c(2, 5))], 1, mean)
datagg$CP3 <- datagg[,"C4"]

# G : Aggression


outagg1 <- fa(datagg[,paste0("G", 1:40)], 1)
outagg2 <- fa(datagg[,paste0("G", 1:40)], 2)
outagg3 <- fa(datagg[,paste0("G", 1:40)], 3)
outagg4 <- fa(datagg[,paste0("G", 1:40)], 4)
outagg5 <- fa(datagg[,paste0("G", 1:40)], 5)
rbind(do.call(c, outagg1[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outagg2[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outagg3[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outagg4[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]),
      do.call(c, outagg5[c("chi", "dof", "PVAL", "RMSEA", "TLI", "BIC")]))
fa.parallel(datagg[,paste0("G", 1:40)])
outagg3 <- fa(datagg[,paste0("G", 1:40)], 3, rotate="bifactor")
loadings(outagg3)


outagg3a <- fa(datagg[,paste0("G", 1:40)], 3, rotate="quartimin")
loadings(outagg3a)
outagg2 <- fa(datagg[,paste0("G", 1:40)], 2, rotate="quartimin")
loadings(outagg2)

datagg$GP1 <- apply(datagg[,paste0("G", c(3, 8, 11, 14, 16, 20, 25, 30, 36, 37))], 1, mean)
datagg$GP2 <- apply(datagg[,paste0("G", c(1, 5, 12, 18, 21, 24, 27, 28, 33, 35))], 1, mean)
datagg$GP3 <- apply(datagg[,paste0("G", c(4, 7, 10, 15, 17, 22, 31, 34, 38, 40))], 1, mean)
datagg$GP4 <- apply(datagg[,paste0("G", c(2, 6, 9, 13, 19, 23, 26, 29, 32, 39))], 1, mean)

scriptagg <- '
mind =~ AP1 + AP2 + AP3 + AP4
friendsupport =~ BPA1 + BPA2 + BPA3
friendalienation =~ BPB1 + BPB2 + BPB3
school =~ CP1 + CP2 + CP3
agg =~ GP1 + GP2 + GP3 + GP4
'
outagg <- cfa(scriptagg, data=datagg)
summary(outagg, fit=TRUE, std=TRUE)

scriptagg2 <- '
mind =~ AP1 + AP2 + AP3 + AP4
friendsupport =~ BPA1 + BPA2 + BPA3
friendalienation =~ BPB1 + BPB2 + BPB3
school =~ CP1 + CP2 + CP3
agg =~ GP1 + GP2 + GP3 + GP4
agg ~ mind + friendsupport + friendalienation + school
'
outagg2 <- cfa(scriptagg2, data=datagg)
summary(outagg2, std=TRUE, rsquare=TRUE)


scriptc <- '
cc =~ C1 + C2 + C3 
'
outc <- cfa(scriptc, data=datagg, ordered=TRUE)
summary(outc, std=TRUE, rsquare=TRUE)
scriptc2 <- '
cc =~ C1 + C2
cc2 =~  C3 + C4
'
outc2 <- cfa(scriptc2, data=datagg, ordered=TRUE)
summary(outc2, std=TRUE, rsquare=TRUE)

library(lavaan)
datdicho <- read.table("lecture12ex1.csv", sep=",", header=TRUE)

mdicho <- '
f1 =~ y1 + y2 + y3 + y4
f2 =~ y5 + y6 + y7 + y8
'
outdicho <- cfa(mdicho, data=datdicho, ordered=TRUE)
summary(outdicho, std=TRUE, fit=TRUE)

mdicho <- '
f1 =~ y1 + y2 
f2 =~ y5 + y6 + y7 + y8
'
outdicho <- cfa(mdicho, data=datdicho, ordered=TRUE)
summary(outdicho, std=TRUE, fit=TRUE)