library(semTools)

######### Distribution

skew(1:5)
kurtosis(1:5)
mardiaSkew(HolzingerSwineford1939[,paste("x", 1:9, sep="")])
mardiaKurtosis(HolzingerSwineford1939[,paste("x", 1:9, sep="")])

############## EFA

unrotated2 <- efaUnrotate(HolzingerSwineford1939, nf=2, varList=paste0("x", 1:9))
unrotated3 <- efaUnrotate(HolzingerSwineford1939, nf=3, varList=paste0("x", 1:9))
anova(unrotated2, unrotated3)

orthRotate(unrotated3, method="varimax")

dat <- HolzingerSwineford1939[,paste0("x", 1:9)]
miss <- matrix(rbinom(prod(dim(dat)), 1, 0.2), nrow(dat))
dat[miss == 1] <- NA
dat <- data.frame(dat, z=rnorm(nrow(HolzingerSwineford1939), 0, 1))
miss2 <- efaUnrotate(dat, nf=2, varList=paste0("x", 1:9), missing = "fiml", aux = "z")
summary(oblqRotate(miss2, method = "quartimin"))
miss3 <- efaUnrotate(dat, nf=3, varList=paste0("x", 1:9), missing = "fiml", aux = "z")
summary(oblqRotate(miss3, method = "quartimin"))
anova(miss2, miss3)

############# Reliability

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit <- cfa(HS.model, data=HolzingerSwineford1939)
reliability(fit)

fit3 <- cfa(HS.model, data=HolzingerSwineford1939, estimator="MLR", group="school", group.equal="loadings")
reliability(fit3)

library(psych)
dat <- iqitems
for(i in 1:ncol(iqitems)) {
	dat[,i] <- ordered(iqitems[,i])
}
iq.model <- '
reason =~ reason.4 + reason.16 + reason.17 + reason.19
letter =~ letter.7 + letter.33 + letter.34 + letter.58
matrix =~ matrix.45 + matrix.46 + matrix.47 + matrix.55
rotate =~ rotate.3 + rotate.4 + rotate.6 + rotate.8
'
fit4 <- cfa(iq.model, data=dat)
reliability(fit4) # Should provide a warning for coefficient alpha

HS.model3 <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 
			  higher =~ visual + textual + speed'

fit6 <- cfa(HS.model3, data=HolzingerSwineford1939)
reliability(fit6) # Should provide a warning for the endogenous variable
reliabilityL2(fit6, "higher")

################################# Maximal Reliability


mod1 <- ' visual  =~ x1 + x2 + x3'
fit <- cfa(mod1, data=HolzingerSwineford1939)
maximalRelia(fit)

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit2 <- cfa(HS.model, data=HolzingerSwineford1939)
maximalRelia(fit2)

fit3 <- cfa(HS.model, data=HolzingerSwineford1939, group="school")
maximalRelia(fit3)
	
################ Probing interaction ####################

#orthogonalize for residual centering

dat2wayMC <- indProd(PoliticalDemocracy, 1:3, 9:11, match = TRUE)
# dat2wayDMC <- indProd(PoliticalDemocracy, 1:3, 9:11, doubleMC=TRUE, match = TRUE)
# dat2wayRC <-orthogonalize(PoliticalDemocracy,(1:3), (9:11), match = TRUE)

library(lavaan) 
model1 <- "
f1 =~ x1 + x2 + x3
f2 =~ y1 + y2 + y3
f12 =~ y1.x1 + y2.x2 + y3.x3
f3 =~ y4 + y5 + y6 + y7 + y8
f3 ~ f1 + f2 + f12
f12 ~~ f1
f12 ~~ f2
x1 ~ 0*1
y1 ~ 0*1
y1.x1 ~ 0*1
y4 ~ 0*1
f1 ~ NA*1
f2 ~ NA*1
f12 ~ NA*1
f3 ~ NA*1
y1.x1 ~~ y1
y1.x1 ~~ x1
y2.x2 ~~ y2
y2.x2 ~~ x2
y3.x3 ~~ y3
y3.x3 ~~ x3
"

fitMC2way <- sem(model1, data=dat2wayMC, meanstructure=TRUE, std.lv=FALSE)
summary(fitMC2way)
result2wayMC <- probe2WayMC(fitMC2way, c("f1", "f2", "f12"), "f3", "f2", c(-1, 0, 1))
result2wayMC

plotProbe(result2wayMC, xlim=c(-2, 2))

######### measurementInvariance

library(lavaan)
HW.model <- ' visual =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed =~ x7 + x8 + x9 '

measurementInvariance(HW.model, data=HolzingerSwineford1939, group="school", strict = TRUE)

model <- ' f1 =~ u1 + u2 + u3 + u4
           f2 =~ u5 + u6 + u7 + u8'

measurementInvarianceCat(model, data = datCat, group = "g", parameterization="theta", estimator="wlsmv")

############### Longitudinal Measurement Invariance ###################################

model <- ' f1t1 =~ y1t1 + y2t1 + y3t1
              f1t2 =~ y1t2 + y2t2 + y3t2
			  f1t3 =~ y1t3 + y2t3 + y3t3'

var1 <- c("y1t1", "y2t1", "y3t1")
var2 <- c("y1t2", "y2t2", "y3t2")
var3 <- c("y1t3", "y2t3", "y3t3")
constrainedVar <- list(var1, var2, var3)
longInvariance(model, auto=1, constrainAuto=TRUE, varList=constrainedVar, data=exLong)

################################ Partial Invariance Tests


conf <- "
f1 =~ NA*x1 + x2 + x3
f2 =~ NA*x4 + x5 + x6
f1 ~~ c(1, 1)*f1
f2 ~~ c(1, 1)*f2
"

weak <- "
f1 =~ NA*x1 + x2 + x3
f2 =~ NA*x4 + x5 + x6
f1 ~~ c(1, NA)*f1
f2 ~~ c(1, NA)*f2
"

configural <- cfa(conf, data = HolzingerSwineford1939, std.lv = TRUE, group="school")
weak <- cfa(weak, data = HolzingerSwineford1939, group="school", group.equal="loadings")
models <- list(fit.configural = configural, fit.loadings = weak)
partialInvariance(models, "metric")

partialInvariance(models, "metric", free = "x5") # "x5" is free across groups in advance
partialInvariance(models, "metric", fix = "x4") # "x4" is fixed across groups in advance

# Use the result from the measurementInvariance function
HW.model <- ' visual =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed =~ x7 + x8 + x9 '

models2 <- measurementInvariance(HW.model, data=HolzingerSwineford1939, group="school", strict = TRUE)

partialInvariance(models2, "metric")
partialInvariance(models2, "scalar")
partialInvariance(models2, "strict")
partialInvariance(models2, "means")

################################ Partial Invariance Cat Tests

f <- rnorm(1000, 0, 1)
u1 <- 0.9*f + rnorm(1000, 1, sqrt(0.19))
u2 <- 0.8*f + rnorm(1000, 1, sqrt(0.36))
u3 <- 0.6*f + rnorm(1000, 1, sqrt(0.64))
u4 <- 0.7*f + rnorm(1000, 1, sqrt(0.51))
u1 <- as.numeric(cut(u1, breaks = c(-Inf, 0, Inf)))
u2 <- as.numeric(cut(u2, breaks = c(-Inf, 0.5, Inf)))
u3 <- as.numeric(cut(u3, breaks = c(-Inf, 0, Inf)))
u4 <- as.numeric(cut(u4, breaks = c(-Inf, -0.5, Inf)))
g <- rep(c(1, 2), 500)
dat2 <- data.frame(u1, u2, u3, u4, g)

configural2 <- "
f1 =~ NA*u1 + u2 + u3 + u4
u1 | c(t11, t11)*t1 
u2 | c(t21, t21)*t1 
u3 | c(t31, t31)*t1 
u4 | c(t41, t41)*t1 
f1 ~~ c(1, 1)*f1
f1 ~ c(0, NA)*1
u1 ~~ c(1, 1)*u1
u2 ~~ c(1, NA)*u2
u3 ~~ c(1, NA)*u3
u4 ~~ c(1, NA)*u4
"

outConfigural2 <- cfa(configural2, data = dat2, group = "g", parameterization="theta", estimator="wlsmv", ordered = c("u1", "u2", "u3", "u4"))

weak2 <- "
f1 =~ NA*u1 + c(f11, f11)*u1 + c(f21, f21)*u2 + c(f31, f31)*u3 + c(f41, f41)*u4
u1 | c(t11, t11)*t1 
u2 | c(t21, t21)*t1 
u3 | c(t31, t31)*t1 
u4 | c(t41, t41)*t1 
f1 ~~ c(1, NA)*f1
f1 ~ c(0, NA)*1
u1 ~~ c(1, 1)*u1
u2 ~~ c(1, NA)*u2
u3 ~~ c(1, NA)*u3
u4 ~~ c(1, NA)*u4
"

outWeak2 <- cfa(weak2, data = dat2, group = "g", parameterization="theta", estimator="wlsmv", ordered = c("u1", "u2", "u3", "u4"))
modelsCat <- list(fit.configural = outConfigural2, fit.loadings = outWeak2)

partialInvarianceCat(modelsCat, type = "metric") 
partialInvarianceCat(modelsCat, type = "metric", free = "u2") 
partialInvarianceCat(modelsCat, type = "metric", fix = "u3") 

######### moreFitIndices

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit <- cfa(HS.model, data=HolzingerSwineford1939, fixed.x = FALSE)
moreFitIndices(fit)

########### Saris, Satorra, and van der Veld (2009)

library(lavaan)

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit <- cfa(HS.model, data=HolzingerSwineford1939, group="sex", meanstructure=TRUE)
miPowerFit(fit)

model <- ' 
  # latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a*y2 + b*y3 + c*y4
     dem65 =~ y5 + a*y6 + b*y7 + c*y8

  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'
fit2 <- sem(model, data=PoliticalDemocracy, meanstructure=TRUE)
miPowerFit(fit2, stdLoad=0.3, cor=0.2, stdBeta=0.2, intcept=0.5)

################ Power Analysis ################################################

plotRMSEApower(rmsea0=.05, rmseaA=.075, df=23, nlow=100, nhigh=500, steps=10)

plotRMSEAdist(rmsea=c(.05, .08), n=200, df=20, ptile=0.95, rmseaScale = TRUE)
plotRMSEAdist(rmsea=c(.05, .01), n=200, df=20, ptile=0.05, rmseaScale = FALSE)

findRMSEApower(rmsea0=.05, rmseaA=.08, df=20, n=200)

findRMSEAsamplesize(rmsea0=.05, rmseaA=.08, df=20, power=0.80)


#################### outClipboard #####################

library(lavaan)
HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit <- cfa(HS.model, data=HolzingerSwineford1939, group="school", meanstructure=TRUE)

clipboard(fit)
clipboard(fit, "mifit")
clipboard(fit, "coef")
clipboard(fit, "se")
clipboard(fit, "samp")
clipboard(fit, "fit")

########### Auxiliary

library(lavaan)

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '
			  
dat <- data.frame(HolzingerSwineford1939, z=rnorm(nrow(HolzingerSwineford1939), 0, 1))
			  
fit <- cfa(HS.model, data=dat, meanstructure=TRUE) #, group="sex", meanstructure=TRUE)
fitaux <- auxiliary(fit, data=dat, aux="z", fun="cfa", missing="ml")

################### runMI function #########################


HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

HSMiss <- HolzingerSwineford1939[,paste("x", 1:9, sep="")]
randomMiss <- rbinom(prod(dim(HSMiss)), 1, 0.1)
randomMiss <- matrix(as.logical(randomMiss), nrow=nrow(HSMiss))
HSMiss[randomMiss] <- NA

out <- cfa.mi(HS.model, data=HSMiss, m = 3, chi="all")

summary(out)
inspect(out, "fit")
inspect(out, "impute")
standardizedSolution(out)

outscaled <- cfa.mi(HS.model, data=HSMiss, m = 3, chi="all", estimator="mlm")
summary(outscaled)
inspect(outscaled, "fit")
inspect(outscaled, "impute")

HS.model2 <- ' visual  =~ x1 + a*x2 + a*x3
              textual =~ x4 + b*x5 + b*x6
              speed   =~ x7 + c*x8 + c*x9 '

out2 <- cfa.mi(HS.model2, data=HSMiss, m = 3, chi="all")

anova(out, out2)

############# Multivariate Wald Test

HS.model <- ' visual  =~ x1 + con1*x2 + con1*x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + con2*x8 + con2*x9 '

fit <- cfa(HS.model, data=HolzingerSwineford1939)
wald(fit, "con2 - con1")

model.syntax <- '
  # intercept and slope with fixed coefficients
    i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
    s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4

  # regressions
    i ~ x1 + x2
    s ~ x1 + x2

  # time-varying covariates
    t1 ~ c1
    t2 ~ c2
    t3 ~ c3
    t4 ~ c4
'

fit2 <- growth(model.syntax, data=Demo.growth)
wald.syntax <- '
i~x1 - i~x2
1/2*s~x1 - 1/2*s~x2
'
wald(fit2, wald.syntax)

model3 <- ' f1  =~ x1 + p2*x2 + p3*x3 + p4*x4 + p5*x5 + p6*x6
			  p4 == 2*p2'

fit3 <- cfa(model3, data=HolzingerSwineford1939)
wald(fit3, "p3; p6 - 0.5*p5")

################################# Single Parameter Test

library(lavaan)
HW.model <- ' visual =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed =~ x7 + x8 + x9 '

models <- measurementInvariance(HW.model, data=HolzingerSwineford1939, group="school")
singleParamTest(models[[1]], models[[2]])
singleParamTest(models[[3]], models[[4]])