library(simstandard)

genmodel <- "
f1 =~ 0.6*x1 + 0.6*x2 + 0.6*x3 + 0.6*x4
f2 =~ 0.6*x5 + 0.6*x6 + 0.6*x7 + 0.6*x8
f3 =~ 0.6*x8 + 0.6*x9 + 0.6*x10 + 0.6*x11
f2 ~ 0.3*f1 
f3 ~ 0.1*f1 + -0.2*f2
"
popdat <- sim_standardized(genmodel, n = 1000000)
popdat <- popdat[,paste0("x", 1:11)]
model1 <- "
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7 + x8
f3 =~ x8 + x9 + x10 + x11
f2 ~ a*f1
f3 ~ c*f1 + b*f2
ind := a*b
"
out <- sem(model1, popdat)
simout <- sim(nRep=20, model=model1, n=2000, rawData=popdat, lavaanfun="cfa")
summary(simout)
summaryParam(simout)
summaryParam(simout, std=TRUE)
summaryParam(simout, detail=TRUE)



library(simsem)

path <- matrix(0, 3, 3)
path[2,1] <- 0.3
path[3,1] <- 0.1
path[3,2] <- -0.2
resFacVar <- findFactorResidualVar(path, diag(3), rep(1, 3))
resFacVar
facCov <- findFactorTotalCov(path, diag(3), resFacVar)
loading <- matrix(0, 11, 3)
loading[1:4, 1] <- 0.6
loading[5:8, 2] <- 0.6
loading[8:11, 3] <- 0.6
findIndResidualVar(loading, facCov, rep(1, 11))

mod <- '
f1 =~ 0.6*x1 + 0.6*x2 + 0.6*x3 + 0.6*x4
f2 =~ 0.6*x5 + 0.6*x6 + 0.6*x7 + 0.6*x8
f3 =~ 0.6*x8 + 0.6*x9 + 0.6*x10 + 0.6*x11
f1 ~~ 1*f1
f2 ~~ 0.910*f2
f3 ~~ 0.962*f3
f2 ~ 0.3*f1 
f3 ~ 0.1*f1 + -0.2*f2


'



library(lavaan)
library(psych)
describe(HolzingerSwineford1939) 


library(semTools)
mardiaSkew(HolzingerSwineford1939[,paste0("x", 1:9)])
mardiaKurtosis(HolzingerSwineford1939[,paste0("x", 1:9)])

mhs <- ' 
visual  =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 
'
ouths <- cfa(mhs, data = HolzingerSwineford1939, estimator="mlm")
summary(ouths, fit.measures = TRUE)

mhs2 <- ' 
visual  =~ 1*x1 + 1*x2 + 1*x3
textual =~ 1*x4 + 1*x5 + 1*x6
speed   =~ 1*x7 + 1*x8 + 1*x9 
'
ouths2 <- cfa(mhs2, data = HolzingerSwineford1939, estimator="mlm")
anova(ouths, ouths2)


datsurvey <- read.table("lecture14survey.csv", sep=",", header=TRUE)
head(datsurvey)
table(datsurvey$fpc1)
with(datsurvey, table(group, fpc2))

mdisagg <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6
'
outdisagg <- cfa(mdisagg, data=datsurvey)
summary(outdisagg, fit=TRUE, std=TRUE)

# Dummy coding
datsurvey$groupfac <- factor(datsurvey$group) 
datsurvey$groupfac <- relevel(datsurvey$groupfac, ref="4")
dummies <- model.matrix(lm(id ~ groupfac, data=datsurvey))
dummies2 <- cbind(dummies, group=datsurvey$group)
head(dummies2)
tail(dummies2)
colnames(dummies) <- c("intcept", "school1vs4", "school2vs4", "school3vs4", "school5vs4")
datsurvey <- data.frame(datsurvey, dummies)


mfixed <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6
y1 ~ school1vs4 + school2vs4 + school3vs4 + school5vs4
y2 ~ school1vs4 + school2vs4 + school3vs4 + school5vs4
y3 ~ school1vs4 + school2vs4 + school3vs4 + school5vs4
y4 ~ school1vs4 + school2vs4 + school3vs4 + school5vs4
y5 ~ school1vs4 + school2vs4 + school3vs4 + school5vs4
y6 ~ school1vs4 + school2vs4 + school3vs4 + school5vs4
'
outfixed <- cfa(mfixed, data=datsurvey)
summary(outfixed, fit=TRUE, std=TRUE)


dummieseff <- model.matrix(lm(id ~ groupfac, data=datsurvey, 
                           contrasts = list(groupfac = contr.sum)))
dummieseff2 <- cbind(dummieseff, datsurvey$group)
head(dummieseff2)
tail(dummieseff2)
colnames(dummieseff) <- c("intcepteff", "school1eff", "school2eff", "school3eff", "school4eff")
datsurvey <- data.frame(datsurvey, dummieseff)

mfixedeff <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6
y1 ~ school1eff + school2eff + school3eff + school4eff
y2 ~ school1eff + school2eff + school3eff + school4eff
y3 ~ school1eff + school2eff + school3eff + school4eff
y4 ~ school1eff + school2eff + school3eff + school4eff
y5 ~ school1eff + school2eff + school3eff + school4eff
y6 ~ school1eff + school2eff + school3eff + school4eff
'
outfixedeff <- cfa(mfixedeff, data=datsurvey)
summary(outfixedeff, fit=TRUE, std=TRUE)

library(survey)
designdat <- svydesign(id=~group+id, data=datsurvey, fpc=~fpc1+fpc2)
designdat
datsurvey$w <- weights(designdat)
mweights <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6
'
outweights <- cfa(mweights, data=datsurvey, sampling.weights="w")
summary(outweights, fit=TRUE, std=TRUE)

datsurveyg <- aggregate(cbind(y1, y2, y3, y4, y5, y6) ~ group, 
                        data=datsurvey, FUN=mean)
datsurveyg$f1 <- with(datsurveyg, y1 + y2 + y3)
datsurveyg$f2 <- with(datsurveyg, y4 + y5 + y6)
cor(datsurveyg$f1, datsurveyg$f2)

mmsem <- '
level: 1
fw1 =~ y1 + y2 + y3
fw2 =~ y4 + y5 + y6
level: 2
fb1 =~ y1 + y2 + y3
fb2 =~ y4 + y5 + y6
'
outmsem <- sem(mmsem, data = datsurvey, cluster = "group")
summary(outmsem)



datcon <- read.table("lecture11consci.csv", sep=",", header=TRUE, na.strings="999")
mcon <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55'

datcon2 <- datcon[,paste0("c", c(1, 7, 13, 19, 25, 31, 37, 43, 49, 55))]
mh <- mahalanobis(datcon2, center=apply(datcon2, 2, mean), cov=cov(datcon2))
mh <- data.frame(id=datcon$id, mh=mh)
mh <- mh[order(mh$mh, decreasing=TRUE),]
mh$p <- pchisq(mh$mh, df=10, lower.tail=FALSE)
head(mh)

library(faoutlier)
set.seed(123321)
mhrobust <- robustMD(datcon2)
mhrobust

library(lavaan)

outcongof <- GOF(datcon2, mcon)
outcongof
outcongcd <- gCD(datcon2, mcon)
outcongcd

mcon <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55'
outcon <- cfa(mcon, data=datcon2)
outcon408 <- cfa(mcon, data=datcon[-408,])
outcon282 <- cfa(mcon, data=datcon[-282,])
outcon455 <- cfa(mcon, data=datcon[-455,])
fitmeasures(outcon)["chisq"]
fitmeasures(outcon408)["chisq"]
fitmeasures(outcon282)["chisq"]
fitmeasures(outcon455)["chisq"]


set.seed(123321)
datcon <- read.table("lecture11consci.csv", sep=",", header=TRUE, na.strings="999")
datconcopy <- datcon
n <- nrow(datconcopy)
pmissing <- 1/(1 + exp(-(-1 - 1*scale(datconcopy$Age))))
datconcopy[runif(n) < pmissing, "c49"] <- NA

library(mice)
md.pattern(datconcopy[,setdiff(colnames(datconcopy), "Job")])

library(lavaan)
mcon <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55'
outcon <- cfa(mcon, data=datconcopy)
summary(outcon, fit=TRUE, std=TRUE)


outconfiml <- cfa(mcon, data=datconcopy, missing="ml")
summary(outconfiml, fit=TRUE, std=TRUE)


datconcopy$dgender2 <- datconcopy$Gender == 2
datconcopy$dgender3 <- datconcopy$Gender == 3
mconaux <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
dgender2 ~~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
dgender3 ~~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
Age ~~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
Education ~~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
Salary ~~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
dgender2 ~~ dgender3 + Age + Education + Salary
dgender3 ~~ Age + Education + Salary
Age ~~ Education + Salary
Education ~~ Salary
'
outconaux <- cfa(mconaux, data=datconcopy, missing="ml")
summary(outconaux, fit=TRUE, std=TRUE)

mconfixed <- '
achi =~ c1 + c7 + c13 + c19 + c25 + c31 + c37 + c43 + c49 + c55
c1 ~ dgender2 + dgender3 + Age + Education + Salary
c7 ~ dgender2 + dgender3 + Age + Education + Salary
c13 ~ dgender2 + dgender3 + Age + Education + Salary
c19 ~ dgender2 + dgender3 + Age + Education + Salary
c25 ~ dgender2 + dgender3 + Age + Education + Salary
c31 ~ dgender2 + dgender3 + Age + Education + Salary
c37 ~ dgender2 + dgender3 + Age + Education + Salary
c43 ~ dgender2 + dgender3 + Age + Education + Salary
c49 ~ dgender2 + dgender3 + Age + Education + Salary
c55 ~ dgender2 + dgender3 + Age + Education + Salary
'
outconfixed <- cfa(mconfixed, data=datconcopy, missing="ml", fixed.x=FALSE)
summary(outconfixed, fit=TRUE, std=TRUE)


msem <- '
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
sat =~ o5 + o6
sat ~ achi + caut + duti + orde + disc + effi
'
outsem <- sem(msem, data=datcon)
inspect(outsem, "cov.lv")
covlv <- inspect(outsem, "cov.lv")
library(psych)
lmCor(achi ~ caut + duti + orde + disc + effi, data=cov2cor(covlv), n.obs=475)$R2
lmCor(caut ~ achi + duti + orde + disc + effi, data=cov2cor(covlv), n.obs=475)$R2
lmCor(duti ~ achi + caut + orde + disc + effi, data=cov2cor(covlv), n.obs=475)$R2
lmCor(orde ~ achi + caut + duti + disc + effi, data=cov2cor(covlv), n.obs=475)$R2
lmCor(disc ~ achi + caut + duti + orde + effi, data=cov2cor(covlv), n.obs=475)$R2
lmCor(effi ~ achi + caut + duti + orde + disc, data=cov2cor(covlv), n.obs=475)$R2