### DATA PREPRATION ###

#uncomment two lines below if packages are not installed
#install.packages("lavaan", dependencies=TRUE)
#install.packages("lme4", dependencies=TRUE)

#load libraries
library(lavaan)
library(lme4)

#load datasets
widedat <- read.csv("https://stats.idre.ucla.edu/wp-content/uploads/2021/08/gpa_wide.csv")
longdat <- read.csv("https://stats.idre.ucla.edu/wp-content/uploads/2021/08/gpa_long.csv")
hsbdemo <- read.csv("https://stats.idre.ucla.edu/wp-content/uploads/2021/07/hsbdemo.csv", stringsAsFactors=TRUE)

### PART 1: GROWTH MODEL ###

#The latent growth model (LGM)
m1 <-   'i =~ 1*gpa0 + 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 
         s =~ 0*gpa0 + 1*gpa1 + 2*gpa2 + 3*gpa3 + 4*gpa4'
fit_m1 <- growth(m1, data=widedat)
summary(fit_m1)

#Ordinal versus measured time in an LGM
m2 <-   'i =~ 1*gpa0 + 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 
         s =~ 0*gpa0 + 3*gpa1 + 6*gpa2 + 9*gpa3 + 12*gpa4'
fit_m2 <- growth(m2, data=widedat)
summary(fit_m2)

#Equivalence of the LGM to the hierarchical linear model (HLM)
m3 <- lmer(gpa ~ time + (time|student),dat=longdat)
summary(m3)
#variance covariance matrix
as.data.frame(VarCorr(m3))

#LGM with homogenous residual variance
m4 <-   'i =~ 1*gpa0 + 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4  
          s =~ 0*gpa0 + 1*gpa1 + 2*gpa2 + 3*gpa3 + 4*gpa4  
          gpa0 ~~ a*gpa0
          gpa1 ~~ a*gpa1
          gpa2 ~~ a*gpa2
          gpa3 ~~ a*gpa3
          gpa4 ~~ a*gpa4'
fit_m4 <- growth(m4, data=widedat)
summary(fit_m4)

#Adding a predictor to the latent growth model (LGM)
m5 <-   'i =~ 1*gpa0 + 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4  
         s =~ 0*gpa0 + 1*gpa1 + 2*gpa2 + 3*gpa3 + 4*gpa4
         i ~ sex
         s ~ sex
         gpa0 ~~ a*gpa0
         gpa1 ~~ a*gpa1
         gpa2 ~~ a*gpa2
         gpa3 ~~ a*gpa3
         gpa4 ~~ a*gpa4'
fit_m5 <- growth(m5, data=widedat)
summary(fit_m5)

#Exercise 1 (challenge), try to match m5 with cfa() instead of growth() in lavaan
#what is wrong with this output?
m6 <-    'i =~ 1*gpa0 + 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4 
          s =~ 0*gpa0 + 1*gpa1 + 2*gpa2 + 3*gpa3 + 4*gpa4
          i ~ sex
          s ~ sex
          #constrain the residual variances
          gpa0 ~~ a*gpa0
          gpa1 ~~ a*gpa1
          gpa2 ~~ a*gpa2
          gpa3 ~~ a*gpa3
          gpa4 ~~ a*gpa4'
fit_m6 <- cfa(m6, data=widedat,meanstructure=TRUE)
summary(fit_m6)

#correct, matches m5
m7 <-    'i =~ 1*gpa0 + 1*gpa1 + 1*gpa2 + 1*gpa3 + 1*gpa4  
          s =~ 0*gpa0 + 1*gpa1 + 2*gpa2 + 3*gpa3 + 4*gpa4
          #observed intercepts constrained to 0
          gpa0 ~ 0 
          gpa1 ~ 0 
          gpa2 ~ 0 
          gpa3 ~ 0 
          gpa4 ~ 0 
          #estimate intercept and slope means, add gender as a predictor
          i ~ 1 + sex
          s ~ 1 + sex
          #constrain the residual variances
          gpa0 ~~ a*gpa0
          gpa1 ~~ a*gpa1
          gpa2 ~~ a*gpa2
          gpa3 ~~ a*gpa3
          gpa4 ~~ a*gpa4'
fit_m7 <- cfa(m7, data=widedat, meanstructure=TRUE)
summary(fit_m7)

m8 <- lmer(gpa ~ time + sex + time:sex + (time|student),dat=longdat)
summary(m8)
#variance covariance matrix
as.data.frame(VarCorr(m8))

### PART2: MEASUREMENT INVARIANCE ###

#separate cfa models
onefac <- 'f1  =~ read + write + math + science'

#equivalent meanstructure
onefac_b <- 'f1  =~ read + write + math + science
        read ~ 1
        write ~ 1
        math ~ 1
        science ~ 1'

#split dataset by gender
femaledat<- subset(hsbdemo,female=="female")
maledat<- subset(hsbdemo,female=="male")
dim(femaledat)
dim(maledat)

fitfemale <- cfa(onefac, data = femaledat, meanstructure = TRUE)
summary(fitfemale, standardized=TRUE)

fitmale <- cfa(onefac, data = maledat, meanstructure = TRUE)
summary(fitmale, standardized=TRUE)

# configural invariance
fit.Configural <- cfa(onefac, data = hsbdemo, group = "female", meanstructure = TRUE)
summary(fit.Configural, standardized=TRUE)
fitMeasures(fit.Configural, "cfi")

#  metric (weak) invariance 
fit.Metric <- cfa(onefac, data = hsbdemo, group = "female", group.equal = c("loadings"), meanstructure = TRUE)
summary(fit.Metric, standardized=TRUE)
fitMeasures(fit.Metric, "cfi")

#Exercise 2A
#one factor, variance std method (WRONG) 
onefac_varstd <- 'f1  =~ NA*read + write + math + science
           f1 ~~ c(1,1)*f1'

# metric invariance (WRONG)
fit.MetricWRONG <- cfa(onefac_varstd , data = hsbdemo, group = "female", group.equal=c("loadings"),meanstructure = TRUE)
summary(fit.MetricWRONG)

#one factor, variance std method (CORRECT) 
onefac_varstd <- 'f1  =~ NA*read + write + math + science
           f1 ~~ c(1,NA)*f1'
# metric invariance (CORRECT)
fit.MetricCORRECT <- cfa(onefac_varstd , data = hsbdemo, group = "female", group.equal=c("loadings"),meanstructure = TRUE)
summary(fit.MetricCORRECT)

#Exercise 2B
#metric (weak) invariance, manual variance std method 
onefac_varstd2 <- 'f1  =~ c(NA,NA)*read + c(a,a)*read + c(b,b)*write + c(c,c)*math + c(d,d)*science
                   f1 ~~ c(1,NA)*f1'

#metric (weak) invariance, manual variance std method
fit.MetricVS <- cfa(onefac_varstd2 , data = hsbdemo, group = "female", meanstructure = TRUE)
summary(fit.MetricVS)

# scalar (strong) invariance
fit.Scalar <- cfa(onefac, data = hsbdemo, group = "female", group.equal = c("loadings","intercepts"), meanstructure = TRUE)
summary(fit.Scalar, standardized=TRUE)
fitMeasures(fit.Scalar, "cfi")

# scalar (strong) invariance manual method
onefac_strong <- 'f1  =~ c(1,1)*read + c(a,a)*write + c(b,b)*math + c(c,c)*science
                  read ~ c(d,d)*1
                  write ~ c(e,e)*1
                  math ~ c(f,f)*1
                  science ~ c(g,g)*1
                  f1 ~ c(0,NA)*1'
fit.Residual <- cfa(onefac_strong, data = hsbdemo, group = "female", meanstructure = TRUE)
summary(fit.Residual, standardized=TRUE, fit.measures=TRUE)

# residual invariance
fit.Residual <- cfa(onefac, data = hsbdemo, group = "female", group.equal = c("loadings","intercepts","residuals"), meanstructure = TRUE)
summary(fit.Residual, standardized=TRUE)
fitMeasures(fit.Residual, "cfi")

#residual invariance manual method
onefac_resid <- 'f1  =~ c(1,1)*read + c(a,a)*write + c(b,b)*math + c(c,c)*science
                  #strong invariance
                  read ~ c(d,d)*1
                  write ~ c(e,e)*1
                  math ~ c(f,f)*1
                  science ~ c(g,g)*1
                  f1 ~ c(0,NA)*1
                  #residual invariance
                  read ~~ c(h,h)*read
                  write ~~ c(i,i)*write
                  math ~~ c(j,j)*math
                  science ~~ c(k,k)*science'
fit.Residual2 <- cfa(onefac_resid, data = hsbdemo, group = "female", meanstructure = TRUE)
summary(fit.Residual2, standardized=TRUE)

# model comparison tests
lavTestLRT(fit.Configural, fit.Metric, fit.Scalar, fit.Residual)

#manual chi-square difference test
pchisq(q=4.179,df=3,lower.tail=FALSE)
#alternatively
1-pchisq(q=4.179,df=3,lower.tail=TRUE)

# partial measurement invariance
onefac_d <- 'f1  =~ read + write + c(a,b)*math + science
        read ~ 1
        write ~ 1
        math ~ 1
        science ~ 1'
fit.Partial <- cfa(onefac_d, data = hsbdemo, group = "female", group.equal = c("loadings"), meanstructure = TRUE)
summary(fit.Partial, standardized=TRUE)
lavTestLRT(fit.Configural, fit.Partial)

#Exercise 3A
onefac_c <- 'f1  =~ c(a,b)*read + write + math + science
        read ~ 1
        write ~ 1
        math ~ 1
        science ~ 1'
fit.Partial2 <- cfa(onefac_c, data = hsbdemo, group = "female", group.equal = c("loadings"), meanstructure = TRUE)
summary(fit.Partial2, standardized=TRUE)

#Exercise 3B (challenge)
onefac_e <- 'f1  =~ c(NA,NA)*read + c(1,1)*write + math + science
        read ~ 1
        write ~ 1
        math ~ 1
        science ~ 1'
fit.Partial3 <- cfa(onefac_e, data = hsbdemo, group = "female", group.equal = c("loadings"), meanstructure = TRUE, auto.fix.first = FALSE)
summary(fit.Partial3, standardized=TRUE)
lavTestLRT(fit.Configural, fit.Partial3)