library(simsem)

genmodel <- "
f1 =~ 0.7*x1 + 0.7*x2 + 0.7*x3
f2 =~ 0.8*x4 + 0.9*x5
f2 ~~ 0.3*f1 
f1 ~~ 1*f1
f2 ~~ 1*f2
x1 ~~ 0.51*x1
x2 ~~ 0.51*x2
x3 ~~ 0.51*x3
x4 ~~ 0.36*x4
x5 ~~ 0.19*x5
"

analysismodel <- "
f1 =~ x1 + x2 + x3
f2 =~ x4 + x5
"

cfa(analysismodel, data=generate(genmodel, 200))
simout200 <- sim(nRep=1000, model=analysismodel, generate=genmodel, n=200, lavaanfun="cfa", std.lv=TRUE)
summary(simout200)
summaryParam(simout200, detail=TRUE, improper=FALSE)

simout400 <- sim(nRep=100, model=analysismodel, generate=genmodel, n=400, lavaanfun="cfa", std.lv=TRUE)
summary(simout400)

simout600 <- sim(nRep=100, model=analysismodel, generate=genmodel, n=600, lavaanfun="cfa", std.lv=TRUE)
summary(simout600)

simout800 <- sim(nRep=100, model=analysismodel, generate=genmodel, n=800, lavaanfun="cfa", std.lv=TRUE)
summary(simout800)

simout1000 <- sim(nRep=100, model=analysismodel, generate=genmodel, n=1000, lavaanfun="cfa", std.lv=TRUE)
summary(simout1000)

genalternativemodel <- "
f1 =~ 0.7*x1 + 0.7*x2
f2 =~ 0.7*x3 + 0.8*x4 + 0.9*x5
f2 ~~ 0.3*f1 
f1 ~~ 1*f1
f2 ~~ 1*f2
x1 ~~ 0.51*x1
x2 ~~ 0.51*x2
x3 ~~ 0.51*x3
x4 ~~ 0.36*x4
x5 ~~ 0.19*x5
"

simoutalt200 <- sim(nRep=1000, model=analysismodel, generate=genalternativemodel, 
                    n=200, lavaanfun="cfa", std.lv=TRUE)
summary(simoutalt200)
summaryParam(simoutalt200, detail=TRUE)
summaryFit(simoutalt200, improper=FALSE)
summaryConverge(simoutalt200)
ruleThumb <- c(RMSEA=0.05, CFI=0.95, TLI=0.95, SRMR=0.06)
plotPowerFit(simoutalt200, cutoff=ruleThumb, alpha=0.05, 
             usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
getPowerFit(simoutalt200, cutoff=ruleThumb, alpha=0.05, 
            usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
plotCutoff(simout200, alpha=0.05, 
             usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
getCutoff(simout200, alpha=0.05, 
          usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
dfiCutoff <- getCutoff(simout200, alpha=0.05, 
                       usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
getPowerFit(simoutalt200, cutoff=dfiCutoff, alpha=0.05, 
            usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
plotPowerFit(simoutalt200, cutoff=dfiCutoff, alpha=0.05, 
            usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))

##############################
# Example 2

library(simstandard)
set.seed(123321)
genmodel <- "
f1 =~ 0.6*x1 + 0.6*x2 + 0.6*x3 + 0.6*x4
f2 =~ 0.6*x5 + 0.6*x6 + 0.6*x7 + 0.45*x8
f3 =~ 0.45*x8 + 0.6*x9 + 0.6*x10 + 0.6*x11
f2 ~ 0.3*f1 
f3 ~ 0.1*f1 + -0.2*f2
"
popdat <- sim_standardized(genmodel, n = 1000000, observed=TRUE,
                           latent=FALSE, errors=FALSE)
head(popdat)

medmodel <- "
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7 + x8
f3 =~ x8 + x9 + x10 + x11
f2 ~ a*f1
f3 ~ c*f1 + b*f2
ind := a*b
"
paramout <- sem(medmodel, popdat)
parameterEstimates(paramout)

simout200 <- sim(nRep=1000, model=medmodel, n=200, rawData=popdat, lavaanfun="sem")

pv <- simout200@paramValue
pv[,names(coef(paramout, type="all"))] <- coef(paramout, type="all")
simout200@paramValue <- pv
summary(simout200)

summaryParam(simout200)

simout700 <- sim(nRep=100, model=medmodel, n=700, rawData=popdat, lavaanfun="sem")
pv <- simout700@paramValue
pv[,names(coef(paramout, type="all"))] <- coef(paramout, type="all")
simout700@paramValue <- pv
summary(simout700)

simout2000 <- sim(nRep=100, model=medmodel, n=2000, rawData=popdat, lavaanfun="sem")
pv <- simout2000@paramValue
pv[,names(coef(paramout, type="all"))] <- coef(paramout, type="all")
simout2000@paramValue <- pv
summary(simout2000)

summaryParam(simout700, detail=TRUE)[c("a", "b", "c", "ind"),]
summaryParam(simout2000, detail=TRUE)[c("a", "b", "c", "ind"),]

simout4000 <- sim(nRep=100, model=medmodel, n=4000, rawData=popdat, lavaanfun="sem")
pv <- simout4000@paramValue
pv[,names(coef(paramout, type="all"))] <- coef(paramout, type="all")
simout4000@paramValue <- pv
summaryParam(simout4000, detail=TRUE)[c("a", "b", "c", "ind"),]

##############################
# Example 3

library(simsem)

genmodel <- "
f1 =~ 0.6*x1 + 0.6*x2 + 0.6*x3 + 0.6*x4
f2 =~ 0.6*x5 + 0.6*x6 + 0.6*x7 + 0.6*x8
f2 ~~ 0.7*f1 
f1 ~~ 1*f1
f2 ~~ 1.2*f2
x1 ~~ 0.64*x1
x2 ~~ 0.64*x2
x3 ~~ 0.64*x3
x4 ~~ 0.64*x4
x5 ~~ 0.64*x5
x6 ~~ 0.64*x6
x7 ~~ 0.64*x7
x8 ~~ 0.64*x8
f1 ~ 0*1
f2 ~ 0.5*1
x1 ~ 3*1
x2 ~ 3*1
x3 ~ 3*1
x4 ~ 3*1
x5 ~ 3*1
x6 ~ 3*1
x7 ~ 3*1
x8 ~ 3*1
"

analysismodel <- "
f1 =~ x1 + L2*x2 + L3*x3 + L4*x4
f2 =~ x5 + L2*x6 + L3*x7 + L4*x8
f1 ~~ NA*f2
f1 ~~ NA*f1 + V1*f1
f2 ~~ NA*f2 + V2*f2
x1 ~ 0*1
x2 ~ I2*1
x3 ~ I3*1
x4 ~ I4*1
x5 ~ 0*1
x6 ~ I2*1
x7 ~ I3*1
x8 ~ I4*1
f1 ~ M1*1
f2 ~ M2*1
eff := (M2 - M1)/sqrt(V1)
"

simout200 <- sim(nRep=1000, model=analysismodel, generate=genmodel, n=200, lavaanfun="cfa")
summary(simout200)
summaryParam(simout200, detail=TRUE)["eff",]
mean(2*1.96*simout200@se[,"eff"])

simout120 <- sim(nRep=100, model=analysismodel, generate=genmodel, n=120, lavaanfun="cfa")
summaryParam(simout120, detail=TRUE)["eff",]
mean(2*1.96*simout120@se[,"eff"])

library(semTools)

plotRMSEAdist(0.05, n = 200, df = 4, ptile = .95, rmseaScale = TRUE)
findRMSEApower(rmsea0 = 0.05, rmseaA = 0.08, df = 4, n = 200, alpha = 0.05, group = 1)
plotRMSEApower(rmsea0 = 0.05, rmseaA = 0.08, df = 4,
               nlow = 100, nhigh = 2000, steps = 10)
abline(h=0.8)
findRMSEAsamplesize(rmsea0 = .05, rmseaA = .08, df = 4, power = 0.80)

library(MBESS)
ci.rmsea(rmsea=.05, df=4, N=200, conf.level=.95)
ci.rmsea(rmsea=.05, df=4, N=1803, conf.level=.95)

ss.aipe.rmsea(RMSEA=.05, df=4, width=0.04, conf.level=.95)
ss.aipe.rmsea(RMSEA=.05, df=40, width=0.04, conf.level=.90)