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)

# library(dynamic)
# dynamiccon <- cfaHB(outcon, plot = TRUE)
# readr::write_rds(dynamiccon, "lecture11_dynamiccon.rds")
dynamiccon <- readr::read_rds("lecture11_dynamiccon.rds")

psych::describe(datcon$Age)
hist(datcon$Age, breaks=15)

mcon2 <- '
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
achi ~ Age
caut ~ Age
duti ~ Age
orde ~ Age
disc ~ Age
effi ~ Age
'
outcon2 <- cfa(mcon2, data=datcon)
summary(outcon2, fit=TRUE, std=TRUE)
# Std.all means all variables including exogenous variables are standardized
anova(outcon2, outcon) # This is wrong!!!

mcon2a <- '
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
achi ~ 0*Age
caut ~ 0*Age
duti ~ 0*Age
orde ~ 0*Age
disc ~ 0*Age
effi ~ 0*Age
'
outcon2a <- cfa(mcon2a, data=datcon)
anova(outcon2, outcon2a)

library(nonnest2)
vuongtest(outcon2, outcon2a, nested=TRUE)

table(datcon$Gender)

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

mcon3 <- '
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
achi ~ male + ugender + agec
caut ~ male + ugender + agec
duti ~ male + ugender + agec
orde ~ male + ugender + agec
disc ~ male + ugender + agec
effi ~ male + ugender + agec
'
outcon3 <- cfa(mcon3, data=datcon)
summary(outcon3, fit=TRUE, std=TRUE)

mcon3a <- '
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
achi ~ 0*male + 0*ugender + agec
caut ~ 0*male + 0*ugender + agec
duti ~ 0*male + 0*ugender + agec
orde ~ 0*male + 0*ugender + agec
disc ~ 0*male + 0*ugender + agec
effi ~ 0*male + 0*ugender + agec
'
outcon3a <- cfa(mcon3a, data=datcon)
anova(outcon3, outcon3a)

mcon3b <- '
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
achi ~ male + ugender + 0*agec
caut ~ male + ugender + 0*agec
duti ~ male + ugender + 0*agec
orde ~ male + ugender + 0*agec
disc ~ male + ugender + 0*agec
effi ~ male + ugender + 0*agec
'
outcon3b <- cfa(mcon3b, data=datcon)
anova(outcon3, outcon3b)

library(nonnest2)
vuongtest(outcon3, outcon3a, nested=TRUE)
vuongtest(outcon3, outcon3b, nested=TRUE)


datcon$male_agec <- datcon$male * datcon$agec
datcon$ugen_agec <- datcon$ugender * datcon$agec

mcon4 <- '
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
achi ~ male + ugender + agec + male_agec + ugen_agec
caut ~ male + ugender + agec + male_agec + ugen_agec
duti ~ male + ugender + agec + male_agec + ugen_agec
orde ~ male + ugender + agec + male_agec + ugen_agec
disc ~ male + ugender + agec + male_agec + ugen_agec
effi ~ male + ugender + agec + male_agec + ugen_agec
'
outcon4 <- cfa(mcon4, data=datcon)
summary(outcon4, fit=TRUE, std=TRUE)

mcon4a <- '
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
achi ~ male + ugender + agec + 0*male_agec + 0*ugen_agec
caut ~ male + ugender + agec + 0*male_agec + 0*ugen_agec
duti ~ male + ugender + agec + 0*male_agec + 0*ugen_agec
orde ~ male + ugender + agec + 0*male_agec + 0*ugen_agec
disc ~ male + ugender + agec + 0*male_agec + 0*ugen_agec
effi ~ male + ugender + agec + 0*male_agec + 0*ugen_agec
'
outcon4a <- cfa(mcon4a, data=datcon)
anova(outcon4, outcon4a)
vuongtest(outcon4, outcon4a, nested=TRUE)

mconprune1 <- '
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
achi ~ male + ugender + agec
caut ~ male + ugender + agec
duti ~ male + ugender + agec + male_agec + ugen_agec
orde ~ male + ugender + agec
disc ~ male + ugender + agec
effi ~ male + ugender + agec
'
outconprune1 <- cfa(mconprune1, data=datcon)
summary(outconprune1, fit=TRUE, std=TRUE)



mconprune2 <- '
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
duti ~ male + ugender + agec + male_agec + ugen_agec
orde ~ male + ugender
disc ~ agec
'
outconprune2 <- cfa(mconprune2, data=datcon)
summary(outconprune2, 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
burn =~ o7 + o8
o1 ~ achi + caut + duti + orde + disc + effi
o4 ~ achi + caut + duti + orde + disc + effi
sat ~ achi + caut + duti + orde + disc + effi
burn ~ achi + caut + duti + orde + disc + effi
o1 ~~ o4 + sat + burn
o4 ~~ sat + burn
sat ~~ burn
'
outsem <- sem(msem, data=datcon)
summary(outsem, std=TRUE, fit=TRUE)

mcfa2 <- '
sat =~ o5 + o6
burn =~ o7 + o8
sat ~~ o1 + o4 + burn
burn ~~ o1 + o4
o1 ~~ o4
'
outcfa2 <- cfa(mcfa2, data=datcon)
summary(outcfa2, std=TRUE, fit=TRUE)


msem2 <- '
sat =~ o5 + o6
burn =~ o7 + o8
sat ~ o1 + o4
burn ~ sat + o1 + o4
o1 ~~ o4
'
outsem2 <- sem(msem2, data=datcon)
summary(outsem2, std=TRUE, fit=TRUE)

msem3 <- '
sat =~ o5 + o6
burn =~ o7 + o8
sat ~ o1 + o4
burn ~ sat
o1 ~~ o4
'
outsem3 <- sem(msem3, data=datcon)
library(nonnest2)
vuongtest(outsem2, outsem3, nested = TRUE)

msem4 <- '
sat =~ o5 + o6
burn =~ o7 + o8
sat ~ o1 + o4
burn ~ sat + o1
o1 ~~ o4
'
outsem4 <- sem(msem4, data=datcon)
vuongtest(outsem2, outsem4, nested = TRUE)
vuongtest(outsem3, outsem4, nested = TRUE)
summary(outsem4)

msem5 <- '
sat =~ o5 + o6
burn =~ o7 + o8
sat ~ a1*o1 + a2*o4
burn ~ b*sat + c1*o1 + c2*o4
ind1 := a1*b
ind2 := a2*b
tot1 := c1 + ind1
tot2 := c2 + ind2
'
outsem5 <- sem(msem5, data=datcon, se="boot")
summary(outsem5, std=TRUE, fit=TRUE)
parameterEstimates(outsem5)




dat1 <- read.table("lecture11ex1.csv", sep=",", header=TRUE)

mcfa <- '
f =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
'
outcfa <- cfa(mcfa, data=dat1)
summary(outcfa, std=TRUE)

msem <- '
f =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
f ~ x1 + x2
'
outsem <- sem(msem, data=dat1)
summary(outsem, std=TRUE)

mcfa2 <- '
f =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
f ~~ x1 + x2
x1 ~~ x2
'
outcfa2 <- cfa(mcfa2, data=dat1)
summary(outcfa2, std=TRUE)

mcfa3 <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6 + y7 + y8
'
outcfa3 <- cfa(mcfa3, data=dat1)
summary(outcfa3, std=TRUE)

mcfa4 <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6 + y7 + y8
f1 ~~ f2 + x1 + x2
f2 ~~ x1 + x2 
x1 ~~ x2
'
outcfa4 <- cfa(mcfa4, data=dat1)
summary(outcfa4, std=TRUE)



msam5 <- '
f =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
f ~ x1 + x2
'
outsam5 <- sam(msam5, data=dat1)
summary(outsam5, std=TRUE)
inspect(outsam5, "est")

msem5 <- '
f =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
f ~ x1 + x2
'
outsem5 <- sem(msem5, data=dat1)
inspect(outsem5, "est")
summary(outsem5, std=TRUE)

mcfa5 <- '
f =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
'
outcfa5 <- cfa(mcfa5, data=dat1)
inspect(outcfa5, "est")

library(QuantPsyc)
dat1$sumall <- with(dat1, y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)
summary(lm("sumall ~ x1 + x2", data=dat1))
lm.beta(lm("sumall ~ x1 + x2", data=dat1))

outcfa5a <- cfa(mcfa5, data=dat1, bounds = "wide")

outsam5a <- sam(msam5, data=dat1, mm.args = list(bounds = "none"))


msam2 <- '
f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6 + y7 + y8
f1 ~~ f2 + x1 + x2
f2 ~~ x1 + x2 
x1 ~~ x2
'
outsam2 <- sam(msam2, data=dat1)
summary(outsam2, std=TRUE)

mcfa21 <- '
f1 =~ y1 + y2 + y3
'
mcfa22 <- '
f2 =~ y4 + y5 + y6 + y7 + y8
'

msem2b <- '
sat =~ o5 + o6
leave =~ o7 + o8
sat ~ o1 + o4
leave ~ sat
'
outsem2b <- sem(msem2b, data=datcon)
summary(outsem2b, std=TRUE, fit=TRUE)
anova(outsem2, outsem2b)














################# Not used yet

datcon$edyear <- NA
datcon[datcon$Education == 1, "edyear"] <- 3
datcon[datcon$Education == 2, "edyear"] <- 9
datcon[datcon$Education == 3, "edyear"] <- 12
datcon[datcon$Education == 4, "edyear"] <- 14
datcon[datcon$Education == 5, "edyear"] <- 16
datcon[datcon$Education == 6, "edyear"] <- 18
datcon[datcon$Education == 7, "edyear"] <- 21

psych::describe(datcon$edyear)

table(datcon$Job)

datcon$jobstudent <- datcon$Job == 2
datcon$jobprivate <- datcon$Job == 3
datcon$jobowner <- datcon$Job == 4
datcon$jobgov <- datcon$Job == 5
datcon$jobfreelancer <- datcon$Job == 6

datcon$salarynum <- NA
datcon[!is.na(datcon$Salary) & datcon$Salary == 1, "salarynum"] <- 5
datcon[!is.na(datcon$Salary) & datcon$Salary == 2, "salarynum"] <- 12.5
datcon[!is.na(datcon$Salary) & datcon$Salary == 3, "salarynum"] <- 17.5
datcon[!is.na(datcon$Salary) & datcon$Salary == 4, "salarynum"] <- 22.5
datcon[!is.na(datcon$Salary) & datcon$Salary == 5, "salarynum"] <- 27.5
datcon[!is.na(datcon$Salary) & datcon$Salary == 6, "salarynum"] <- 35
datcon[!is.na(datcon$Salary) & datcon$Salary == 7, "salarynum"] <- 45
datcon[!is.na(datcon$Salary) & datcon$Salary == 8, "salarynum"] <- 62.5
datcon[!is.na(datcon$Salary) & datcon$Salary == 9, "salarynum"] <- 87.5
datcon[!is.na(datcon$Salary) & datcon$Salary == 10, "salarynum"] <- 112.5

psych::describe(datcon$salarynum)

mmimic <- '
f =~ o5 + o6 + o7 + o8
f ~ male + ugender + agec + edyear 
f ~ jobstudent + jobprivate + jobowner + jobgov + jobfreelancer  + salarynum
'
outmimic <- sem(mmimic, data=datcon)
summary(outmimic, std=TRUE, fit=TRUE)

mmimic2 <- '
f =~ o5 + o6 + o7 + o8
'
outmimic2 <- sem(mmimic2, data=datcon)
summary(outmimic2, std=TRUE, fit=TRUE)