# Item 2
library(simsem)
set.seed(123321)
loading <- matrix(0, 7, 3)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loading[7, 3] <- NA
LY <- bind(loading, 0.7)

VTE <- bind(c(rep(NA, 6), 0), 0.51)
RTE <- binds(diag(7))

RPS <- binds(diag(3))
path <- matrix(0, 3, 3)
path[2, 1] <- NA
path[3, 2] <- NA
pathVal <- matrix(0, 3, 3)
pathVal[2, 1] <- 0.6
pathVal[3, 2] <- 0.7
BE <- bind(path, pathVal)

mod <- model.sem(LY=LY, VTE=VTE, RTE=RTE, RPS=RPS, BE=BE)
dat <- generate(mod, 200)
dat <- data.frame(id = 1:nrow(dat), dat)
dat <- dat[,-ncol(dat)]
last <- dat[,ncol(dat)]
last <- 100 * (last - min(last)) / (max(last) - min(last))
dat[,ncol(dat)] <- last

write.csv(dat, file="hw3i2.csv", row.names=FALSE)

# Solution

colnames(dat) <- c("id", "potency1", "potency2", "potency3", "collect1", "collect2", "collect3", "perform")

script2a <- "
collect =~ collect1 + collect2 + collect3
potency =~ potency1 + potency2 + potency3
per =~ perform
perform ~~ 0*perform
collect ~~ per
potency ~~ per
collect ~~ potency
"
measurement <- cfa(script2a, data=dat)

script2b <- "
collect =~ collect1 + collect2 + collect3
potency =~ potency1 + potency2 + potency3
per =~ perform
perform ~~ 0*perform
per ~ b*collect
per ~ potency
collect ~ a*potency
indirect := a * b
"
partial <- sem(script2b, data=dat)

partialboot <- sem(script2b, data=dat, se="boot")
parameterEstimates(partialboot, boot.ci.type="bca.simple")

script2c <- "
collect =~ collect1 + collect2 + collect3
potency =~ potency1 + potency2 + potency3
per =~ perform
perform ~~ 0*perform
per ~ b*collect
collect ~ a*potency
indirect := a * b
"
full <- sem(script2c, data=dat)

anova(full, partial)

# Item 3

library(simsem)
set.seed(123321)
loading.in <- matrix(0, 9, 3)
loading.in[1:3, 1] <- paste0("load", 1:3)
loading.in[4:6, 2] <- paste0("load", 4:6)
loading.in[7:9, 3] <- paste0("load", 7:9)
LY.in <- bind(loading.in, 0.7)

latent.cor1 <- matrix(NA, 3, 3)
diag(latent.cor1) <- 1
latent.val1 <- matrix(NA, 3, 3)
diag(latent.val1) <- 1
latent.val1[lower.tri(latent.val1)] <- latent.val1[upper.tri(latent.val1)] <- c(0.63, 0.73, 0.64)
RPS1 <- binds(latent.cor1, latent.val1)

latent.cor2 <- matrix(NA, 3, 3)
diag(latent.cor2) <- 1
latent.val2 <- matrix(NA, 3, 3)
diag(latent.val2) <- 1
latent.val2[lower.tri(latent.val2)] <- latent.val2[upper.tri(latent.val2)] <- c(0.58, 0.65, 0.40)
RPS2 <- binds(latent.cor2, latent.val2)

RTE <- binds(diag(9))

VTE <- bind(rep(NA, 9), 0.51)

TY.in <- bind(paste0("int", 1:9), 0)

VPS1 <- bind(rep(1, 3))
VPS2 <- bind(rep(NA, 3), c(1.1, 1.5, 1))

AL1 <- bind(rep(0, 3))
AL2 <- bind(rep(NA, 3), c(0.5, 0, 0.8))

strong <- model(LY = LY.in, RPS = list(RPS1, RPS2), VPS=list(VPS1, VPS2), RTE = RTE, VTE=VTE, TY=TY.in, AL=list(AL1, AL2), ngroups=2, modelType = "CFA")

dat <- generate(strong,list(300, 150))

# create score with seven-point scale

toLikert <- function(vec, start, end) {
	vec <- (end - start) * (vec - min(vec)) / (max(vec) - min(vec)) + start
	round(vec)
}
dat[,1:9] <- lapply(dat[,1:9], toLikert, start = 1, end = 7)

dat <- data.frame(id = 1:nrow(dat), dat)
write.csv(dat, file="hw3i3.csv", row.names=FALSE)

# Solution

colnames(dat) <- c("id", "vs1", "vs2", "vs3", "ss1", "ss2", "ss3", "rm1", "rm2", "rm3", "country")

script3a <- "
vs =~ vs1 + vs2 + vs3
ss =~ ss1 + ss2 + ss3
rm =~ rm1 + rm2 + rm3
"
configural <- cfa(script3a, data=dat, group="country")
weak <- cfa(script3a, data=dat, group="country", group.equal="loadings")
strong <- cfa(script3a, data=dat, group="country", group.equal=c("loadings", "intercepts"))

script3b <- "
vs =~ vs1 + vs2 + vs3
ss =~ ss1 + ss2 + ss3
rm =~ rm1 + rm2 + rm3	  
vs ~~ c(vara1, vara2)*vs
ss ~~ c(varb1, varb2)*ss
rm ~~ c(varc1, varc2)*rm
vs ~~ c(covab1, covab2)*ss
vs ~~ c(covac1, covac2)*rm
ss ~~ c(covbc1, covbc2)*rm
corab1 := covab1 / sqrt(vara1 * varb1)
corab2 := covab2 / sqrt(vara2 * varb2)
corac1 := covac1 / sqrt(vara1 * varc1)
corac2 := covac2 / sqrt(vara2 * varc2)
corbc1 := covbc1 / sqrt(varb1 * varc1)
corbc2 := covbc2 / sqrt(varb2 * varc2)
diffab := corab1 - corab2
diffac := corac1 - corac2
diffbc := corbc1 - corbc2
"
strong2 <- cfa(script3b, data=dat, group="country", group.equal=c("loadings", "intercepts"), se="boot")
parameterEstimates(strong2, boot.ci.type="bca.simple")

# Item 4

library(simsem)
set.seed(123321)
loading <- matrix(0, 21, 7)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loading[7:9, 3] <- NA
loading[10:12, 4] <- NA
loading[13:15, 5] <- NA
loading[16:18, 6] <- NA
loading[19:21, 7] <- NA
LY <- bind(loading, "runif(1, 0.5, 0.9)")

RTE <- binds(diag(21))

ivcor <- matrix(NA, 3, 3)
diag(ivcor) <- 1
latentcor <- diag(7)
latentcor[1:3, 1:3] <- ivcor
latentcor[5:7, 5:7] <- ivcor
RPS <- binds(latentcor, "runif(1, 0.3, 0.7)")

path <- matrix(0, 7, 7)
path[4, 1:3] <- NA
path[5:7, 4] <- NA
pathVal <- matrix(0, 7, 7)
pathVal[4, 1:3] <- c(0.25, 0.45, 0.15)
pathVal[5:7, 4] <- c(0.50, 0.56, -0.53)
BE <- bind(path, pathVal)

mod <- model.sem(LY=LY, RTE=RTE, RPS=RPS, BE=BE)
dat <- generate(mod, 300)
dat <- data.frame(id = 1:nrow(dat), dat)
dat <- dat[,-ncol(dat)]
write.csv(dat, file="hw3i4.csv", row.names=FALSE)

#solution

colnames(dat) <- c("id", "dis1", "dis2", "dis3", "proc1", "proc2", "proc3", "inter1", "inter2", "inter3", "overall1", "overall2", "overall3", "js1", "js2", "js3", "oc1", "oc2", "oc3", "turn1", "turn2", "turn3")

script4a <- "
dis =~ dis1 + dis2 + dis3
proc =~ proc1 + proc2 + proc3
inter =~ inter1 + inter2 + inter3
overall =~ overall1 + overall2 + overall3
"
four <- cfa(script4a, data=dat)

script4b <- "
specific =~ dis1 + dis2 + dis3 + proc1 + proc2 + proc3 + inter1 + inter2 + inter3
overall =~ overall1 + overall2 + overall3
"
two <- cfa(script4b, data=dat)

script4c <- "
overall =~ dis1 + dis2 + dis3 + proc1 + proc2 + proc3 + inter1 + inter2 + inter3+ overall1 + overall2 + overall3
"
one <- cfa(script4c, data=dat)

rbind(fitMeasures(four), fitMeasures(two), fitMeasures(one))

script4d <- "
dis =~ dis1 + dis2 + dis3
proc =~ proc1 + proc2 + proc3
inter =~ inter1 + inter2 + inter3
overall =~ overall1 + overall2 + overall3
js =~ js1 + js2 + js3
oc =~ oc1 + oc2 + oc3
turn =~ turn1 + turn2 + turn3
"
seven <- cfa(script4d, data=dat)

script4e <- "
dis =~ dis1 + dis2 + dis3
proc =~ proc1 + proc2 + proc3
inter =~ inter1 + inter2 + inter3
overall =~ overall1 + overall2 + overall3
js =~ js1 + js2 + js3
oc =~ oc1 + oc2 + oc3
turn =~ turn1 + turn2 + turn3
overall ~ a1*dis + a2*proc + a3*inter
js ~ dis + proc + inter + b1*overall
oc ~ dis + proc + inter + b2*overall
turn ~ dis + proc + inter + b3*overall
js ~~ NA*oc
js ~~ NA*turn
oc ~~ NA*turn
dis ~~ NA*proc
dis ~~ NA*inter
proc ~~ NA*inter
ab11 := a1 * b1
ab12 := a1 * b2
ab13 := a1 * b3
ab21 := a2 * b1
ab22 := a2 * b2
ab23 := a2 * b3
ab31 := a3 * b1
ab32 := a3 * b2
ab33 := a3 * b3
diff12 := a1 - a2
diff13 := a1 - a3
diff23 := a2 - a3
"
partial <- sem(script4e, data=dat)
partialboot <- sem(script4e, data=dat, se="boot")
parameterEstimates(partialboot, boot.ci.type="bca.simple")

script4f <- "
dis =~ dis1 + dis2 + dis3
proc =~ proc1 + proc2 + proc3
inter =~ inter1 + inter2 + inter3
overall =~ overall1 + overall2 + overall3
js =~ js1 + js2 + js3
oc =~ oc1 + oc2 + oc3
turn =~ turn1 + turn2 + turn3
overall ~ a1*dis + a2*proc + a3*inter
js ~ b1*overall
oc ~ b2*overall
turn ~ b3*overall
js ~~ NA*oc
js ~~ NA*turn
oc ~~ NA*turn
dis ~~ NA*proc
dis ~~ NA*inter
proc ~~ NA*inter
ab11 := a1 * b1
ab12 := a1 * b2
ab13 := a1 * b3
ab21 := a2 * b1
ab22 := a2 * b2
ab23 := a2 * b3
ab31 := a3 * b1
ab32 := a3 * b2
ab33 := a3 * b3
"
full <- sem(script4f, data=dat)

anova(full, partial)