# Exercise 13.1

set.seed(123321)
S <- matrix(runif(64, -0.2, 0.2), 8, 8)
S[1:3, 1:3] <- runif(9, 0.5, 0.7)
S[4:8, 4:8] <- runif(25, 0.5, 0.7)
S[lower.tri(S)] <- t(S)[lower.tri(S)]
diag(S) <- 1
S

dat <- MASS::mvrnorm(1000000, rep(0, 8), S, empirical = TRUE)
dat <- data.frame(dat)
colnames(dat) <- paste0("x", 1:8)

out1 <- lm(x1 ~ x2 + x3 + x4 + x5 + x6 + x7 + x8, data = dat)
out2 <- lm(x2 ~ x1 + x3 + x4 + x5 + x6 + x7 + x8, data = dat)
out3 <- lm(x3 ~ x1 + x2 + x4 + x5 + x6 + x7 + x8, data = dat)
out4 <- lm(x4 ~ x1 + x2 + x3 + x5 + x6 + x7 + x8, data = dat)
out5 <- lm(x5 ~ x1 + x2 + x3 + x4 + x6 + x7 + x8, data = dat)
out6 <- lm(x6 ~ x1 + x2 + x3 + x4 + x5 + x7 + x8, data = dat)
out7 <- lm(x7 ~ x1 + x2 + x3 + x4 + x5 + x6 + x8, data = dat)
out8 <- lm(x8 ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = dat)
C <- S
C[1,1] <- summary(out1)$r.squared
C[2,2] <- summary(out2)$r.squared
C[3,3] <- summary(out3)$r.squared
C[4,4] <- summary(out4)$r.squared
C[5,5] <- summary(out5)$r.squared
C[6,6] <- summary(out6)$r.squared
C[7,7] <- summary(out7)$r.squared
C[8,8] <- summary(out8)$r.squared

eigen(C)

vec <- eigen(C)$vectors[,1:2,drop = FALSE]
val <- eigen(C)$values[1:2]
newc <- vec %*% diag(val) %*% t(vec)

rotmat1 <- matrix(c(cos(pi/6), -sin(pi/6), sin(pi/6), cos(pi/6)), 2, 2, byrow = TRUE)
rotmat2 <- matrix(c(cos(pi/3), -sin(pi/3), sin(pi/3), cos(pi/3)), 2, 2, byrow = TRUE)
round(vec %*% solve(rotmat1), 3)
round(vec %*% solve(rotmat2), 3)

# Exercise 13.2

library(psych)
library(GPArotation)
library(lavaan)

out <- factanal(covmat = cov(PoliticalDemocracy), factors = 1, n.obs = 75, rotation = "quartimin")
out2 <- factanal(covmat = cov(PoliticalDemocracy), factors = 2, n.obs = 75, rotation = "quartimin")
out3 <- factanal(covmat = cov(PoliticalDemocracy), factors = 3, n.obs = 75, rotation = "quartimin")
out4 <- factanal(covmat = cov(PoliticalDemocracy), factors = 4, n.obs = 75, rotation = "quartimin")

fa.parallel(PoliticalDemocracy)

out11 <- fa(PoliticalDemocracy, nfactors = 1, fm = "ml")
out21 <- fa(PoliticalDemocracy, nfactors = 2, fm = "ml")
out31 <- fa(PoliticalDemocracy, nfactors = 3, fm = "ml")
out41 <- fa(PoliticalDemocracy, nfactors = 4, fm = "ml")

J <- ncol(PoliticalDemocracy)
result <- list(out11, out21, out31, out41)
sumout <- sapply(result, function(u) c(chi = u$STATISTIC, df = u$dof, u$RMSEA[1:3], TLI = u$TLI))
sumout <- rbind(sumout, p = pchisq(sumout[1,], sumout[2,], lower.tail = FALSE), aic = sumout[1,] + J*(J+1) - 2*sumout[2,])
round(sumout, 3)

chidiff <- sumout[1, 1:3] - sumout[1, 2:4]
dfdiff <- sumout[2, 1:3] - sumout[2, 2:4]
pdiff <- pchisq(chidiff, dfdiff, lower.tail = FALSE)

out32 <- fa(PoliticalDemocracy, nfactors = 3, fm = "ml", rotate = "varimax")
out33 <- fa(PoliticalDemocracy, nfactors = 3, fm = "ml", rotate = "quartimin")
out34 <- fa(PoliticalDemocracy, nfactors = 3, fm = "ml", rotate = "bifactor")

library(MBESS)
data(HS.data)

dat <- HS.data[HS.data$school == "Pasteur",7:25]

fa.parallel(dat)

hs1 <- fa(dat, nfactors = 1, fm = "ml")
hs2 <- fa(dat, nfactors = 2, fm = "ml")
hs3 <- fa(dat, nfactors = 3, fm = "ml")
hs4 <- fa(dat, nfactors = 4, fm = "ml")
hs5 <- fa(dat, nfactors = 5, fm = "ml")
hs6 <- fa(dat, nfactors = 6, fm = "ml")
hs7 <- fa(dat, nfactors = 7, fm = "ml")

J <- ncol(dat)
result <- list(hs1, hs2, hs3, hs4, hs5, hs6, hs7)
sumout <- sapply(result, function(u) c(chi = u$STATISTIC, df = u$dof, u$RMSEA[1:3], TLI = u$TLI))
sumout <- rbind(sumout, p = pchisq(sumout[1,], sumout[2,], lower.tail = FALSE), aic = sumout[1,] + J*(J+1) - 2*sumout[2,])
round(sumout, 3)

# Pick four solutions

hs41 <- fa(dat, nfactors = 4, fm = "ml", rotate = "quartimin")

# Very simple structure with four factors (Loadings > .4 indicate salient items)

dat2 <- HS.data[HS.data$school == "Grant-White",7:25]

syntax <- '
spatial =~ visual + cubes + paper + flags
verbal =~ general + paragrap + sentence + wordc + wordm 
speed =~ addition + code + counting + straight  
memory =~ wordr + numberr + figurer + object + numberf + figurew
'
fit <- cfa(syntax, data = dat2)
summary(fit, fit = TRUE)