rm(list = ls())

#### PACKAGES ####
library(haven)
library(skimr)
library(psych)
library(lavaan)
library(semPlot)
library(semTools)

#### FUNCTIONS ####
cfa_model <- function(model.cfa, data.cfa, estimator.cfa = "MLR", robust = TRUE, decimals = 2, res_cutoff = 0.1) {
  
  message("NOTE: This function only works with Pearson correlations.")
  message(paste("Estimator used:", estimator.cfa))
  message(paste("Fit indices will be rounded to", decimals, "decimals."))
  message(paste("Residual cutoff set to", res_cutoff))
  
  if (robust) {
    message("Robust fit indices are being used.")
  } else {
    message("Non-robust fit indices are being used.")
  }
  
  message("Arguments description:")
  message("$model.cfa: CFA model specification in lavaan syntax.")
  message("$data.cfa: data frame containing observed variables.")
  message("$estimator.cfa: estimator used to fit the CFA model.")
  message("$robust: logical, whether robust fit indices are requested.")
  message("$decimals: number of decimals used to round fit indices.")
  message("$res_cutoff: cutoff for reporting standardized residual correlations.")
  
  all_model <- cfa(
    model = model.cfa,
    data = data.cfa,
    estimator = estimator.cfa
  )
  
  cfa.all <- summary(
    all_model,
    fit.measures = TRUE,
    standardized = TRUE
  )
  
  if (robust) {
    fit.idx <- fitMeasures(
      all_model,
      c(
        "chisq",
        "df",
        "pvalue",
        "cfi.robust",
        "tli.robust",
        "rmsea.robust",
        "srmr"
      )
    )
  } else {
    fit.idx <- fitMeasures(
      all_model,
      c(
        "chisq",
        "df",
        "pvalue",
        "cfi",
        "tli",
        "rmsea",
        "srmr"
      )
    )
  }
  
  usrmr <- lavResiduals(
    all_model,
    type = "cor.bentler"
  )$summary["usrmr", ]
  
  r2 <- inspect(all_model, "r2")
  R2_bar <- mean(unlist(r2))
  usrmr_corr <- usrmr / R2_bar
  
  fit.idx <- c(fit.idx, usrmr = usrmr, usrmr_corr = usrmr_corr)
  fit.idx <- round(fit.idx, decimals)
  
  residuals <- resid(all_model, type = "cor.bentler")$cov
  res <- round(residuals, 2)
  
  idx <- which(
    res > res_cutoff & lower.tri(res),
    arr.ind = TRUE
  )
  
  resid_table <- data.frame(
    item1 = rownames(res)[idx[, 1]],
    item2 = colnames(res)[idx[, 2]],
    residuo = res[idx]
  )
  
  loadings <- standardizedSolution(all_model)
  
  message("Output structure:")
  message("$fit.indices: global model fit indices including USRMR.")
  message("$loadings: standardized factor loadings.")
  message(sprintf(
    "$residuals: standardized residual correlations greater than %s.",
    res_cutoff
  ))
  message("$cfa.all: full lavaan summary with fit measures and standardized estimates.")
  message("$model: fitted lavaan CFA object.")
  
  fit.cfa <- list(
    fit.indices = fit.idx,
    loadings = loadings,
    residuals = resid_table,
    cfa.all = cfa.all,
    model = all_model
  )
  
  return(fit.cfa)
}
invariance <- function(data.cfa, model.cfa, gr, estimator.cfa = "MLR"){
  
  config <- cfa(model.cfa, data.cfa, group = gr, meanstructure = TRUE, estimator = estimator.cfa)
  
  weak <- cfa(model.cfa, data.cfa, group = gr, group.equal="loadings", meanstructure = TRUE, estimator = estimator.cfa)
  
  strong <- cfa(model.cfa, data.cfa, group = gr, group.equal= c("loadings", "intercepts"), meanstructure = TRUE, estimator = estimator.cfa)
  
  strict <- cfa(model.cfa, data.cfa, group = gr, 
                group.equal= c("loadings", "intercepts", "residuals"), meanstructure = TRUE, estimator = estimator.cfa)
  
  return(summary(compareFit(config, weak, strong, strict)))
}

#### LOAD DATA ####
data <- read_sav("data.sav")
data$sexo <- factor(data$sexo - 1, levels = c(0, 1), labels = c("man", "woman"))

#### GD ####
model.gd <- "GD =~ gd1 + gd2 + gd3 + gd4"

gd.cfa <- cfa_model(model.cfa = model.gd, data.cfa = data)

gd.cfa$fit.indices
gd.cfa$loadings
gd.cfa$residuals

#### ASD ####
model.asd <- "
ASD1 =~ asd2 + asd4 + asd5 + asd7 + asd9 + asd10 + asd12 + asd14 + asd16 + asd18 + asd19 + asd20 + asd22 + asd23 + asd25 + asd26
ASD2 =~ asd1 + asd3 + asd6 + asd8 + asd11 + asd13 + asd15 + asd17 + asd21 + asd24"

asd.cfa <- cfa_model(model.cfa = model.asd, data.cfa = data)
asd.cfa$fit.indices
asd.cfa$loadings
asd.cfa$residuals

#### PAQ ####
model.paq <- "
PAQ1 =~ paq2 + paq3 + paq7 + paq8 + paq9 + paq12 + paq15 + paq21 + paq22
PAQ2 =~ paq17 + paq6 + paq10 + paq16 + paq19 + paq20 + paq24"

paq.cfa <- cfa_model(model.cfa = model.paq, data.cfa = data)
paq.cfa$fit.indices
paq.cfa$loadings
paq.cfa$residuals

#### ALL SCALES: GD, ASD, PAQ ####
model.all <-"
GD =~ gd1 + gd2 + gd3 + gd4
ASD1 =~ asd2 + asd4 + asd5 + asd7 + asd9 + asd10 + asd12 + asd14 + asd16 + asd18 + asd19 + asd20 + asd22 + asd23 + asd25 + asd26
ASD2 =~ asd1 + asd3 + asd6 + asd8 + asd11 + asd13 + asd15 + asd17 + asd21 + asd24
PAQ1 =~ paq2 + paq3 + paq7 + paq8 + paq9 + paq12 + paq15 + paq21 + paq22
PAQ2 =~ paq17 + paq6 + paq10 + paq16 + paq19 + paq20 + paq24"

all.cfa <- cfa_model(model.cfa = model.all, data.cfa = data)

all.cfa$fit.indices
all.cfa$loadings
all.cfa$residuals

#### MEASUREMENT INVARIANCE BY SEX ####
invariance(data.cfa = data, model.cfa = model.gd, gr = "sexo")
invariance(data.cfa = data, model.cfa = model.asd, gr = "sexo")
invariance(data.cfa = data, model.cfa = model.paq, gr = "sexo")

#### PARTIAL MEASUREMENT INVARIANCE BY SEX ####
invariance(data.cfa = data, model.cfa = "ASD1 =~ asd2 + asd4 + asd5 + asd7 + asd9 + asd10 + asd12 + asd14 + asd16 + asd18 + asd19 + asd20 + asd22 + asd23 + asd25 + asd26", gr = "sexo")
invariance(data.cfa = data, model.cfa = "ASD2 =~ asd1 + asd3 + asd6 + asd8 + asd11 + asd13 + asd15 + asd17 + asd21 + asd24", gr = "sexo")

invariance(data.cfa = data, model.cfa = "PAQ1 =~ paq2 + paq3 + paq7 + paq8 + paq9 + paq12 + paq15 + paq21 + paq22", gr = "sexo")
invariance(data.cfa = data, model.cfa = "PAQ2 =~ paq17 + paq6 + paq10 + paq16 + paq19 + paq20 + paq24", gr = "sexo")

#### MEASUREMENT INVARIANCE BY AGE ####
invariance(data.cfa = data, model.cfa = model.gd, gr = "edad")
invariance(data.cfa = data, model.cfa = model.asd, gr = "edad")
invariance(data.cfa = data, model.cfa = model.paq, gr = "edad")

#### PARTIAL MEASUREMENT INVARIANCE BY AGE ####
invariance(data.cfa = data, model.cfa = "PAQ1 =~ paq2 + paq3 + paq7 + paq8 + paq9 + paq12 + paq15 + paq21 + paq22", gr = "edad")
invariance(data.cfa = data, model.cfa = "PAQ2 =~ paq17 + paq6 + paq10 + paq16 + paq19 + paq20 + paq24", gr = "edad")

#### CORRELATIONS AMONG ALL FACTORS ####
model.all <- "GD =~ gd1 + gd2 + gd3 + gd4

ASD1 =~ asd2 + asd4 + asd5 + asd7 + asd9 + asd10 + asd12 + asd14 + asd16 + asd18 + asd19 + asd20 + asd22 + asd23 + asd25 + asd26
ASD2 =~ asd1 + asd3 + asd6 + asd8 + asd11 + asd13 + asd15 + asd17 + asd21 + asd24

PAQ1 =~ paq2 + paq3 + paq7 + paq8 + paq9 + paq12 + paq15 + paq21 + paq22
PAQ2 =~ paq17 + paq6 + paq10 + paq16 + paq19 + paq20 + paq24"

all.cfa <- cfa_model(model.cfa = model.all, data.cfa = data)
tail(all.cfa$loadings, 10)[c("lhs", "op", "rhs", "est.std", "pvalue")]