library(lavaan)

dat <- read.table("lecture8ex1.csv", sep=",", header=TRUE)

m1 <- '
f1 =~ x1 + x2 + x3 + x4 + x5
'

out <- cfa(m1, data=dat)
summary(out)
fitmeasures(out)

m2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2 <- cfa(m2, data=dat)
summary(out2)
fitmeasures(out2)

m3 <- '
f1 =~ x1 + x2 + x3 + x4 + x5
f2 =~ x2 + x3 + x4
f1 ~~ 0*f2
'
out3 <- cfa(m3, data=dat)
summary(out3)
fitmeasures(out3)

A <- matrix(c(9, 2, 3,
              4, 5, 6,
              7, 8, 9), 3, 3, byrow=TRUE)
Ainv <- solve(A)
A
Ainv
round(A %*% Ainv, 3)
round(Ainv %*% A, 3)

as.matrix(dat)

X <- as.matrix(dat)
one_n <- as.matrix(rep(1, nrow(X)))
m <- (t(X) %*% one_n) / nrow(X)
m
as.matrix(colMeans(X))

one_n <- matrix(1, nrow(X),1)
dev <- X - one_n %*% t(m)
(t(dev) %*% dev)/(nrow(X)-1)
cov(X)

S <- cov(X)
D <- diag(diag(S))
D2 <- solve(D)^(1/2)
D2 %*% S %*% D2
cor(X)

S1 <- matrix(c(1, 0, 0,
              0, 1, 0,
              0, 0, 1), 3, 3, byrow=TRUE)
det(S1)
S2 <- matrix(c(1, 0.2, 0.3,
               0.2, 1, -0.1,
               0.3, -0.1, 1), 3, 3, byrow=TRUE)
det(S2)
S3 <- matrix(c(1, 0.5, 0.6,
               0.5, 1, 0.7,
               0.6, 0.7, 1), 3, 3, byrow=TRUE)
det(S3)
S4 <- matrix(c(1, 0.95, 0.95,
               0.95, 1, 0.95,
               0.95, 0.95, 1), 3, 3, byrow=TRUE)
det(S4)

S <- matrix(c(1, 0.2, 0.1,
             0.2, 1, 0.3,
             0.1, 0.3, 1), nrow=3)
Sigma1 <- Sigma2 <- Sigma3 <- S
Sigma1[1,2] <- Sigma1[2,1] <- 0
Sigma3[1,2] <- Sigma3[2,1] <- 0.4
c(sum(diag(S %*% solve(Sigma1))), 
  sum(diag(S %*% solve(Sigma2))), 
  sum(diag(S %*% solve(Sigma3))))
c(log(det(Sigma1)) - log(det(S)),
  log(det(Sigma2)) - log(det(S)),
  log(det(Sigma3)) - log(det(S)))
c(sum(diag(S %*% solve(Sigma1))) + log(det(Sigma1)) - log(det(S)),
  sum(diag(S %*% solve(Sigma2))) + log(det(Sigma2)) - log(det(S)),
  sum(diag(S %*% solve(Sigma3))) + log(det(Sigma3)) - log(det(S)))

calcfml <- function(vec, S) {
  lambda <- matrix(0, 5, 2)
  lambda[1,1] <- vec[1]
  lambda[5,1] <- vec[2]
  lambda[2,2] <- vec[3]
  lambda[3,2] <- vec[4]
  lambda[4,2] <- vec[5]
  psi <- diag(2)
  psi[1,2] <- psi[2,1] <- vec[6]
  theta <- diag(vec[7:11])
  Sigma <- lambda %*% psi %*% t(lambda) + theta
  fml <- log(det(Sigma)) + sum(diag(S %*% solve(Sigma))) - log(det(S))
  fml
}

S <- cov(dat)
Sigma <- S
Sigma[1,2] <- Sigma[2,1] <- 0.4
round(S %*% solve(Sigma),2)
diag(S %*% solve(Sigma))
sum(diag(S %*% solve(Sigma)))
log(det(Sigma)) 
log(det(S))
log(det(Sigma)) - log(det(S))

v <- c(0.3, 0.2, 0.5, 0.6, 0.4, 0.7, 0.3, 0.3, 0.2, 0.3, 0.3)
calcfml(v, S)
v2 <- c(0.34144107, 0.25694397, 0.56865379, 0.64570328, 0.36788153, 
        0.83354113, 0.32649055, 0.28933795, 0.16982663, 0.26685149, 0.22520233)
calcfml(v2, S)

x <- optim(v, calcfml, S=S)



m2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2 <- cfa(m2, data=dat, std.lv=TRUE)
summary(out2)
fitmeasures(out2)

calcfml2 <- function(vec, S) {
  lambda <- matrix(0, 5, 2)
  lambda[1,1] <- vec[1]
  lambda[5,1] <- vec[2]
  lambda[2,2] <- vec[3]
  lambda[3,2] <- vec[4]
  lambda[4,2] <- vec[5]
  psi <- diag(2)
  psi[1,2] <- psi[2,1] <- vec[6]
  theta <- diag(vec[7:11])
  Sigma <- lambda %*% psi %*% t(lambda) + theta
  fml <- log(det(Sigma)) + sum(diag(S %*% solve(Sigma))) - log(det(S))
  list(lambda, psi, theta, Sigma, fml)
}
calcfml2(x$par, S)

x <- seq(0, 20, 0.01)
y <- dchisq(x, 4)
plot(x, y, type="l")
crit <- qchisq(0.95, 4)
polygon(c(crit, seq(crit, 20, 0.01), 20), c(0, dchisq(seq(crit, 20, 0.01), 4), 0), col="skyblue")


m2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2 <- cfa(m2, data=dat, std.lv=TRUE)
summary(out2)
fitmeasures(out2)

datneed <- read.table("lecture8manifestneed.csv", sep=",", header=TRUE)
mneed <- '
naff =~ x2 + x6 + x10 + x13
nach =~ x1 + x5 + x9
ndom =~ x4 + x8 + x12 + x15
nauto =~ x3 + x7 + x11 + x14
'
outneed <- cfa(mneed, data=datneed, std.lv=TRUE)
summary(outneed)
fitmeasures(outneed)

m2mv1 <- '
f1 =~ 1*x1 + NA*x5
f2 =~ 1*x2 + NA*x3 + NA*x4
f1 ~~ NA*f1
f2 ~~ NA*f2
f1 ~~ NA*f2
'
out2mv1 <- cfa(m2mv1, data=dat)
summary(out2mv1)

m2ff1 <- '
f1 =~ NA*x1 + NA*x5
f2 =~ NA*x2 + NA*x3 + NA*x4
f1 ~~ 1*f1
f2 ~~ 1*f2
f1 ~~ NA*f2
'
out2ff1 <- cfa(m2ff1, data=dat)
summary(out2ff1)

m2mv2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2mv2 <- cfa(m2mv2, data=dat)
summary(out2mv2)

m2ff2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2ff2 <- cfa(m2ff2, data=dat, std.lv=TRUE)
summary(out2ff2)

m2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2 <- cfa(m2, data=dat, std.lv=TRUE)
summary(out2, std=TRUE)
standardizedSolution(out2)

m2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
out2 <- cfa(m2, data=dat, std.lv=TRUE)
summary(out2, rsquare=TRUE)

msingle0 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3
'
outsingle0 <- cfa(msingle0, data=dat, std.lv=TRUE)
summary(outsingle0)

msingle1 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3
x4 ~~ f1
x4 ~~ f2
'
outsingle1 <- cfa(msingle1, data=dat, std.lv=TRUE)
summary(outsingle1)

msingle2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3
x4 ~~ x1 + x5 + x2 + x3
'
outsingle2 <- cfa(msingle2, data=dat, std.lv=TRUE)
summary(outsingle2)
lavInspect(outsingle2, "theta")
eigen(lavInspect(outsingle2, "theta"))

msingle1a <- '
f1 =~ x1 + x5
f2 =~ x2 + x3
f4 =~ x4
x4 ~~ 0*x4
'
outsingle1a <- cfa(msingle1a, data=dat, std.lv=TRUE)
summary(outsingle1a)

inspect(outsingle1a, "fitted")
inspect(outsingle1a, "sampstat")
residuals(outsingle1a)

lavResiduals(outsingle1a, type = "cor.bollen")
lavResiduals(outsingle1a)

modindices(out2)

m2a <- '
f1 =~ x1 + x5 + x2
f2 =~ x2 + x3 + x4
'
out2a <- cfa(m2a, data=dat, std.lv=TRUE)
summary(out2a)
fitmeasures(out2a)

library(dynamic)
dynamicneed <- dynamic::cfaHB(outneed, plot = TRUE)
dynamicout <- dynamic::cfaOne(out, plot = TRUE)

mtry <- '
f1 =~ x1 + x2 + x6 + x10 
f2 =~ x1 + x5 + x9 + x10
x2 ~~ x5 + x9 + x10
x6 ~~ x5 + x9
x1 ~~  x5
'
outtry <- cfa(mtry, data=datneed, std.lv=TRUE)
outtry
modificationindices(outtry)

mnested <- '
f1 =~ x1 + x2 + x3 + x4 + x5
'
outnested <- cfa(mnested, data=dat, std.lv=TRUE)
outnested
mparent <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
outparent<- cfa(mparent, data=dat, std.lv=TRUE)
outparent
anova(outnested, outparent)
mnested2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
f1 ~~ 1*f2
'
outnested2 <- cfa(mnested2, data=dat, std.lv=TRUE)

mgrandparent <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
x1 ~~ x2
x5 ~~ x3
'
outgrandparent<- cfa(mgrandparent, data=dat, std.lv=TRUE)
outgrandparent
anova(outparent, outgrandparent)
anova(outnested, outgrandparent)

ll <- function(x, vec) {
  mu <- x[1]
  sigma <- x[2]
  -sum(-0.5*((vec-mu)/sigma)^2 - log(sigma*sqrt(2*pi)))
}
optim(c(3.5, 1), ll, vec=c(2,3,4,5))

library(nonnest2)
mnonnested1 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
outnonnested1 <- cfa(mnonnested1, data=dat, std.lv=TRUE)

mnonnested2 <- '
f1 =~ x1 + x4 + x5
f2 =~ x2 + x3
'
outnonnested2 <- cfa(mnonnested2, data=dat, std.lv=TRUE)

fitmeasures(outnonnested1)

fitmeasures(outnonnested1)[c("aic", "bic", "bic2")]
fitmeasures(outnonnested2)[c("aic", "bic", "bic2")]

vuongtest(outnonnested1, outnonnested2)
icci(outnonnested1, outnonnested2)

vuongtest(outnested, outparent, nested=TRUE)


mm1 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
'
outmm1 <- cfa(mm1, data=dat, meanstructure=TRUE)
summary(outmm1)

mm1a <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
x1 ~ NA*1
x2 ~ NA*1
x3 ~ NA*1
x4 ~ NA*1
x5 ~ NA*1
f1 ~ 0*1
f2 ~ 0*1
'
outmm1a <- cfa(mm1a, data=dat)
summary(outmm1a)

mm2 <- '
f1 =~ x1 + x5
f2 =~ x2 + x3 + x4
x1 ~ 0*1
x2 ~ 0*1
x3 ~ NA*1
x4 ~ NA*1
x5 ~ NA*1
f1 ~ NA*1
f2 ~ NA*1
'
outmm2 <- cfa(mm2, data=dat)
summary(outmm2)

mt1 <- '
f1 =~ a*x1 + a*x5
f2 =~ b*x2 + b*x3 + b*x4
'
outmt1 <- cfa(mt1, data=dat, std.lv=TRUE)
summary(outmt1)
anova(out2, outmt1)

mt2 <- '
f1 =~ 1*x1 + 1*x5
f2 =~ 1*x2 + 1*x3 + 1*x4
'
outmt2 <- cfa(mt2, data=dat, std.lv=FALSE)
summary(outmt2)

datadhd <- read.table("lecture8ex2.csv", sep=",", header=TRUE)

madhd1 <- '
f1m =~ x1m + x2m + x3m
f2m =~ x4m + x5m + x6m
f1d =~ x1d + x2d + x3d
f2d =~ x4d + x5d + x6d
x1m ~~ x1d
x2m ~~ x2d
x3m ~~ x3d
x4m ~~ x4d
x5m ~~ x5d
x6m ~~ x6d
'
outadhd1 <- cfa(madhd1, data=datadhd, meanstructure=TRUE)
summary(outadhd1, fit=TRUE, std=TRUE)

madhd2 <- '
f1m =~ x1m + c2*x2m + c3*x3m
f2m =~ x4m + c5*x5m + c6*x6m
f1d =~ x1d + c2*x2d + c3*x3d
f2d =~ x4d + c5*x5d + c6*x6d
x1m ~~ x1d
x2m ~~ x2d
x3m ~~ x3d
x4m ~~ x4d
x5m ~~ x5d
x6m ~~ x6d
'
outadhd2 <- cfa(madhd2, data=datadhd, meanstructure=TRUE)
summary(outadhd2, fit=TRUE, std=TRUE)
anova(outadhd1, outadhd2)

madhd3 <- '
f1m =~ x1m + c2*x2m + c3*x3m
f2m =~ x4m + c5*x5m + c6*x6m
f1d =~ x1d + c2*x2d + c3*x3d
f2d =~ x4d + c5*x5d + c6*x6d
x1m ~~ x1d
x2m ~~ x2d
x3m ~~ x3d
x4m ~~ x4d
x5m ~~ x5d
x6m ~~ x6d
f1m ~ NA*1
f2m ~ NA*1
f1d ~ NA*1
f2d ~ NA*1
x1m ~ 0*1
x2m ~ i2*1
x3m ~ i3*1
x4m ~ 0*1
x5m ~ i5*1
x6m ~ i6*1
x1d ~ 0*1
x2d ~ i2*1
x3d ~ i3*1
x4d ~ 0*1
x5d ~ i5*1
x6d ~ i6*1
'
outadhd3 <- cfa(madhd3, data=datadhd, meanstructure=TRUE)
summary(outadhd3, std=TRUE)
anova(outadhd3, outadhd2)

madhd4 <- '
f1m =~ x1m + c2*x2m + c3*x3m
f2m =~ x4m + c5*x5m + c6*x6m
f1d =~ x1d + c2*x2d + c3*x3d
f2d =~ x4d + c5*x5d + c6*x6d
x1m ~~ x1d
x2m ~~ x2d
x3m ~~ x3d
x4m ~~ x4d
x5m ~~ x5d
x6m ~~ x6d
f1m ~ NA*1
f2m ~ NA*1
f1d ~ NA*1
f2d ~ NA*1
x1m ~ 0*1
x2m ~ i2*1
x3m ~ i3*1
x4m ~ 0*1
x5m ~ i5*1
x6m ~ i6*1
x1d ~ 0*1
x2d ~ i2*1
x3d ~ i3*1
x4d ~ 0*1
x5d ~ i5*1
x6d ~ i6*1
x1m ~~ e1*x1m
x2m ~~ e2*x2m
x3m ~~ e3*x3m
x4m ~~ e4*x4m
x5m ~~ e5*x5m
x6m ~~ e6*x6m
x1d ~~ e1*x1d
x2d ~~ e2*x2d
x3d ~~ e3*x3d
x4d ~~ e4*x4d
x5d ~~ e5*x5d
x6d ~~ e6*x6d
'
outadhd4 <- cfa(madhd4, data=datadhd, meanstructure=TRUE)
summary(outadhd4, std=TRUE)
anova(outadhd4, outadhd3)

mt1a <- '
f1 =~ a1*x1 + a2*x5
f2 =~ b1*x2 + b2*x3 + b3*x4
a1 == a2
b1 == b2
b1 == b3
'
outmt1a <- cfa(mt1a, data=dat, std.lv=TRUE)
summary(outmt1a)

m1eff <- '
f1 =~ a1*x1 + NA*x1 + a2*x5
f2 =~ b1*x2 + NA*x2 + b2*x3 + b3*x4
f1 ~~ v1*f1
f2 ~~ v2*f2
x1 ~ i1*1
x5 ~ i2*1
x2 ~ j1*1
x3 ~ j2*1
x4 ~ j3*1
f1 ~ m1*1 + NA*1
f2 ~ m2*1 + NA*1
# (a1 + a2)/2 == 1
# a1 + a2 == 2
a1 == 2 - a2
# (b1 + b2 + b3)/3 == 1
# b1 + b2 + b3 == 3
b1 == 3 - b2 - b3
# (i1 + i2)/2 == 0
i1 == 0 - i2
# (j1 + j2 + j3)/3 == 0
j1 == 0 - j2 - j3
'
out1eff <- cfa(m1eff, data=dat)
summary(out1eff)

m1inequal <- '
f1 =~ a1*x1 + a2*x5
f2 =~ b1*x2 + b2*x3 + b3*x4
x1 ~~ e1*x1
a1 > 0.01
b1 > 0.01
e1 > 0.01
b3 > b2
'
outm1inequal<- cfa(m1inequal, data=dat, std.lv=TRUE)
summary(outm1inequal)

m2eff <- '
f1 =~ a1*x1 + NA*x1 + a2*x5
f2 =~ b1*x2 + NA*x2 + b2*x3 + b3*x4
f1 ~~ v1*f1
f2 ~~ v2*f2
x1 ~ i1*1
x5 ~ i2*1
x2 ~ j1*1
x3 ~ j2*1
x4 ~ j3*1
f1 ~ m1*1 + NA*1
f2 ~ m2*1 + NA*1
# (a1 + a2)/2 == 1
# a1 + a2 == 2
a1 == 2 - a2
# (b1 + b2 + b3)/3 == 1
# b1 + b2 + b3 == 3
b1 == 3 - b2 - b3
# (i1 + i2)/2 == 0
i1 == 0 - i2
# (j1 + j2 + j3)/3 == 0
j1 == 0 - j2 - j3
mdiff := m1 - m2
vratio := v1/v2 - 1
'
out2eff <- cfa(m2eff, data=dat)
summary(out2eff)

out2eff2 <- cfa(m2eff, data=dat, se="boot", bootstrap=500)
parameterEstimates(out2eff2, boot.ci.type = "bca")

datneed <- read.table("lecture8manifestneed.csv", sep=",", header=TRUE)
mneed <- '
naff =~ x2 + x6 + x10 + x13
nach =~ x1 + x5 + x9
ndom =~ x4 + x8 + x12 + x15
nauto =~ x3 + x7 + x11 + x14
'
outneed <- cfa(mneed, data=datneed, std.lv=TRUE)
fscores <- lavPredict(outneed)
head(fscores)

fscores2 <- lavPredict(outneed, method="Bartlett")
head(fscores2)

lavInspect(outneed, "cov.lv")
cov(fscores)
cor(fscores)
cov(fscores2)
cor(fscores2)

fscores <- lavPredict(outneed, fsm=TRUE)
wfs <- attr(fscores,"fsm")[[1]]
fscoresa <- as.matrix(scale(datneed[,colnames(wfs)], center=TRUE, scale=FALSE)) %*% t(wfs)


newdat <- datneed
newdat <- newdat[1,]
newdat[,paste0("x", 1:15)] <- 1
newdat
lavPredict(outneed, newdata=newdat)

newdat[,paste0("x", 1:15)] <- 5
newdat
lavPredict(outneed, newdata=newdat)

parameterEstimates(outneed)
datneed[,"naff"] <- with(datneed, x2 + x6 + x10 + x13)
datneed[,"nach"] <- with(datneed, x1 + x5 + x9)
datneed[,"ndom"] <- with(datneed, x4 + x8 + x12 + x15)
datneed[,"nauto"] <- with(datneed, x3 + x7 + x11 + x14)
head(datneed)

library(semTools)
compRelSEM(outneed, obs.var=TRUE)
compRelSEM(outneed)
compRelSEM(outneed, tau.eq=TRUE)

fscores <- lavPredict(outneed, fsm=TRUE)
wfs <- attr(fscores,"fsm")[[1]]
lambda <- inspect(outneed, "est")[["lambda"]]
psi <- inspect(outneed, "est")[["psi"]]
n1 <- wfs %*% lambda %*% psi %*% t(lambda) %*% t(wfs)
d1 <- cov(fscores)
diag(n1) / diag(d1)