BECFA.R

# Bayesian Exploratory/Confirmatory CFA:


# Functions:
# becfa: Runs the Bayesian Exploratory/Confirmatory factor analysis model. And compares it vs regular CFA plus frequentist CFAs
# becfa_mi: Runs the measurement invariance steps using ONLY BECFA

becfa_mi <- function(
  model, # Lavaan model
  data,
  group,
  ..., # Lavaan CFA options
  blavaanOptions = list(),
  tol = 0.05,
  tuning = "~ dunif(0.001,100)",
  verbose=TRUE
){
  library("lavaan")
  library("blavaan")
  library("dplyr")
  library("modeest")
  
  blavaanOptions <- c(blavaanOptions,list( std.lv = TRUE, inits = "simple", convergence = "none"))
  
  if (missing(group)){
    stop("'group' needed for measurement invariance tests")
  }
  
  # Fit CFA model:
  # if (verbose) message("Fitting initial CFA model")
  CFAfit <- cfa(model,data=data,group=group, ...)

  # Obtain parameter table:
  pars <- parTable(CFAfit)
  
  # All indicators:
  allInd <- unique(pars$rhs[pars$op == "=~"])
  
  # All latents:
  allLats <- unique(pars$lhs[pars$op == "=~"])
  
  # Remove measurement model from model:
  modelsplit <- strsplit(model,split="\n")[[1]]
  modelsplit <- modelsplit[!grepl("=~",modelsplit)]
  
  cPrior <- "dnorm(0,0.0001)"
  ePrior <- "ddexp(0, TUNING)"
  
  # Create measurement model for each latent:
  mod <- sapply(seq_along(allLats),function(lat){
    # Indicators:
    Ind <- unique(pars$rhs[pars$lhs == allLats[lat] & pars$op == "=~"])
    ordInd <- c(allInd[allInd%in%Ind],allInd[!allInd%in%Ind])
    paste0(allLats[lat]," =~ ",paste0('prior("',ifelse(ordInd %in% Ind,cPrior,ePrior),'") *',ordInd,collapse=" + "))
  })
  
  # New model:
  newModel <- paste(c(mod,modelsplit),collapse="\n")
  
  # Replace first instance of ddexp with trick to compute tuning:
  # newModel <- sub(ePrior,paste0("ddexp(0, TUNING)\\n TUNING",tuning,"\\n"),newModel,fixed = TRUE)
  
  # Extra line:
  extra <- list(syntax = paste0("TUNING",tuning), monitor = "TUNING")

  #### MODEL COMPUTED, TIME FOR MI ####
  if (verbose) message("Estimating configural invariance BECFA")
  BECFA_configural <- do.call("bcfa",c(list(model=newModel,data=data,group=group,
              jagextra = extra),list(...),blavaanOptions))
 
  if (verbose) message("Estimating weak invariance BECFA")
  BECFA_weak <- do.call("bcfa",c(list(model=newModel,data=data,group=group,
           jagextra = extra, group.equal = "loadings"),list(...),blavaanOptions))
  
  if (verbose) message("Estimating strong invariance BECFA")
  BECFA_strong <- do.call("bcfa",c(list(model=newModel,data=data,group=group,
                                      jagextra = extra, 
                                      group.equal = c("loadings","intercepts")
                                      ),list(...),blavaanOptions))
  
  if (verbose) message("Estimating strict invariance BECFA")
  BECFA_strict <- do.call("bcfa",c(list(model=newModel,data=data,group=group,
                                      jagextra = extra, 
                                      group.equal = c("loadings","intercepts","residuals")
  ),list(...),blavaanOptions))
  
  if (verbose) message("Estimating full invariance BECFA")
  BECFA_full <- do.call("bcfa",c(list(model=newModel,data=data,group=group,
                                      jagextra = extra, 
                                      group.equal = c("loadings","intercepts","residuals","means")
  ),list(...),blavaanOptions))
  
  Res <- list(
    configural = BECFA_configural,
    weak = BECFA_weak,
    strong = BECFA_strong,
    strict = BECFA_strict,
    full = BECFA_full
  )
  
  class(Res) <- "BECFA_MI"
  
  return(
    Res
  ) 
}


#### MI methods ####

### Methods:
print.BECFA_MI <- function(object){
  objName <- deparse(substitute(object))
  cat("Bayesian Exploratory/Confirmatory Factor Analysis Measurement Invariance tests completed.")
  cat("\n\nRun",paste0("summary(",objName,") to obtain a data frame of fit measures."))
  cat("\n\nRun",paste0("plot(",objName,", model, group) to plot the estimated posterior modes as a path diagram. Where 'model' is 'configural', 'weak', 'strong','strict' or 'full' and 'group' is a number indicating the group."))
}
  
summary.BECFA_MI <- function(object){
  library("lavaan")
  objName <- deparse(substitute(object))
 
  message("Computing fit measures")
  
  fits <- list()
  N <- length(object)
  pb <- txtProgressBar(0,N,style = 3)
  for (i in 1:N){
    fits[[i]] <- fitMeasures(object[[i]], fit.measures = c("ppp","dic","waic","looic"))
#     
#     if (i == 1){
#       fits[[i]]$BF <- NA
#     } else {
#       fits[[i]]$BF <- BF(object[[i-1]], object[[i]])
#     }
    
    setTxtProgressBar(pb, i)
  }
  close(pb)
  
  FitTable <- as.data.frame(do.call(rbind,fits))
  FitTable$object <- paste0(objName,"$",names(object))
  rownames(FitTable) <- c(
    "Configural", "Weak","Strong","Strict","Full"
  )
  
  return(FitTable)
}



# Plot method:

plot.BECFA_MI <- function(x,model = "configural", group = 1, ...){
  message(paste0("Plotting '",model,"' invariance model for group ",group))
  
  # Obtain parameter estimates:
  parEsts <- parTable(x[[model]])
  Mode <- x[[model]]@external$runjags$summaries[,"Mode"]
  Model <- left_join(parEsts,data.frame(jlabel=names(Mode),mode=Mode,stringsAsFactors=FALSE),by="jlabel")
  semPlotModel <- semPlotModel_blavaanModel(Model)
  semPaths(semPlotModel, "est", "hide",intercepts=FALSE,include=group,ask=FALSE,title = FALSE,...)
}

# This function does:
# 1. Fit CFA model in lavaan
# 2. Fit CFA model in blavaan with peaked priors at 0 for cross-loadings
# 3. Re-fit CFA model in lavaan

becfa <- function(
  model, # Lavaan model
  data,
  ..., # Lavaan CFA options
  blavaanOptions = list( std.lv = TRUE, #ov.cp = "fa", lv.cp = "fa", #test = "none",
                        inits = "Mplus",
                        convergence = "none"), # Additional blavaan options #jagcontrol=list(method="rjparallel"),
  tol = 0.05,
  tuning = "~ dunif(0.001,100)",
  verbose=TRUE
){
  library("lavaan")
  library("blavaan")
  library("dplyr")
  library("modeest")

  # Fit CFA model:
  if (verbose) message("Fitting CFA model")
  CFAfit <- cfa(model,data, ...)

  if (verbose) message("Fitting Bayesian CFA model")
  BCFAfit <- do.call("bcfa",c(list(model=model,data=data),list(...),blavaanOptions))
  # capture.output(BCFAfit <- do.call("bcfa",c(list(model=model,data=data,inits="simple"),list(...),blavaanOptions)))

  # Obtain parameter table:
  pars <- parTable(CFAfit)
  
  # All indicators:
  allInd <- unique(pars$rhs[pars$op == "=~"])
  
  # All latents:
  allLats <- unique(pars$lhs[pars$op == "=~"])
  
  # Remove measurement model from model:
  modelsplit <- strsplit(model,split="\n")[[1]]
  modelsplit <- modelsplit[!grepl("=~",modelsplit)]
  
#   # Start loop:
#   Results <- lapply(tuning,function(t){
#     
 
    cPrior <- "dnorm(0,0.0001)"
    ePrior <- "ddexp(0, TUNING)"
    
    # Create measurement model for each latent:
    mod <- sapply(seq_along(allLats),function(lat){
      # Indicators:
      Ind <- unique(pars$rhs[pars$lhs == allLats[lat] & pars$op == "=~"])
      ordInd <- c(allInd[allInd%in%Ind],allInd[!allInd%in%Ind])
      paste0(allLats[lat]," =~ ",paste0('prior("',ifelse(ordInd %in% Ind,cPrior,ePrior),'") *',ordInd,collapse=" + "))
    })
  
    # New model:
    newModel <- paste(c(mod,modelsplit),collapse="\n")
    
    # Replace first instance of ddexp with trick to compute tuning:
    newModel <- sub(ePrior,paste0("ddexp(0, TUNING)\\n TUNING",tuning,"\\n"),newModel,fixed = TRUE)
    
    # Compute B-CFA:
    if (verbose) message("Fitting Bayesian E/CFA model")
    BECFAfit <- do.call("bcfa",c(list(model=newModel,data=data),list(...),blavaanOptions))
   
  # Parameter estimates:
  parEsts <- parTable(BECFAfit)
  Mode <- BECFAfit@external$runjags$summaries[,"Mode"]
  
  Crossloadings <- left_join(parEsts,data.frame(jlabel=names(Mode),mode=Mode,stringsAsFactors=FALSE),by="jlabel") %>% 
    filter(op == "=~",grepl("ddexp",prior),mode < (0-tol) | mode > (0+tol))
#   
#   if (nrow(Crossloadings) == 0){
#     message("Original CFA model is best model. Not adding cross-loadings")
#     return(CFAfit)
#   } else {
#     message(paste("Adding",nrow(Crossloadings),"Crossloadings to CFA model."))
#     
    mod <- sapply(seq_along(allLats),function(lat){
      # Indicators:
      Ind <- unique(pars$rhs[pars$lhs == allLats[lat] & pars$op == "=~"])
      add <- Crossloadings$rhs[Crossloadings$lhs == allLats[lat]]
      paste0(allLats[lat]," =~ ",paste0(c(Ind,add),collapse=" + "))
    })
    
    # New model:
    newModel <- paste(c(mod,modelsplit),collapse="\n")
    
    # Compute B-CFA:
    if (verbose) message("Re-fitting CFA model")
    CFArefit <- do.call("cfa",c(list(model=newModel,data=data),list(...)))
    
    Res <- list(
      CFA = CFAfit,
      BCFA = BCFAfit,
      BECFA = BECFAfit,
      ECFA = CFArefit
    )
    
    class(Res) <- "BECFA"
    
   return(
    Res
   ) 
  # }

  
}



### Methods:
print.BECFA <- summary.BECFA <- function(object){
  library("lavaan")
  
  fitBCFA <- fitMeasures(object$BCFA)
  fitBECFA <- fitMeasures(object$BECFA)
  
  df <- rbind(fitBECFA,fitBCFA)
  rownames(df) <- c("object$BECFA","object$BCFA")
  cat("Bayesian CFA vs Bayesian E/CFA:\n")
  print(df[,c("npar","ppp","dic","waic","looic")])
  
  
  cat("\n\n")
  
  
  cat("\n============ BECFA results ============\n\n")
  cat("Original CFA vs re-fitted E/CFA:\n")
  print(as.data.frame(lavTestLRT(object$CFA,object$ECFA)))

  
}

# Changed fun from semPlot:
library("semPlot")
semPlotModel_blavaanModel <- function(object, est = "mode", ...)
{
  
  # Check if parTable, otherwise run lavaanify:
  if (!is.data.frame(object) &  !is.list(object))
  {
    object <- lavaanify(object, ...)
  }
  
  varNames <- lavaanNames(object, type="ov")
  factNames <- lavaanNames(object, type="lv")
  #   rm(Lambda)
  
  factNames <- factNames[!factNames%in%varNames]
  
  # Extract number of variables and factors
  n <- length(varNames)
  k <- length(factNames)
  
  # Extract parameter names:
  if (is.null(object$label)) object$label <- rep("",nrow(object))
  
  semModel <- new("semPlotModel")
  
  # Set estimates to 1 or ustart:
  object$est <- ifelse(is.na(object$ustart),1,object$ustart)
  
  if (is.null(object$group)) object$group <- ""
  
  # Create edges dataframe
  semModel@Pars <- data.frame(
    label = object$label,
    lhs = ifelse(object$op=="~"|object$op=="~1",object$rhs,object$lhs),
    edge = "--",
    rhs = ifelse(object$op=="~"|object$op=="~1",object$lhs,object$rhs),
    est = ifelse(object$free==0,object$est,object[[est]]),
    std = NA,
    group = object$group,
    fixed = object$free==0,
    par = object$free,
    stringsAsFactors=FALSE)
  
  semModel@Pars$edge[object$op=="~~"] <- "<->"  
  semModel@Pars$edge[object$op=="~*~"] <- "<->"  
  semModel@Pars$edge[object$op=="~"] <- "~>"
  semModel@Pars$edge[object$op=="=~"] <- "->"
  semModel@Pars$edge[object$op=="~1"] <- "int"
  semModel@Pars$edge[grepl("\\|",object$op)] <- "|"
  
  # Move thresholds to Thresholds slot:
  semModel@Thresholds <- semModel@Pars[grepl("\\|",semModel@Pars$edge),-(3:4)]
  # Remove thresholds from Pars:
  #   semModel@Pars <- semModel@Pars[!grepl("\\|",semModel@Pars$edge),]
  
  # Remove weird edges:
  semModel@Pars <- semModel@Pars[!object$op%in%c(':=','<','>','==','|','<', '>'),]
  
  
  semModel@Vars <- data.frame(
    name = c(varNames,factNames),
    manifest = c(varNames,factNames)%in%varNames,
    exogenous = NA,
    stringsAsFactors=FALSE)
  
  semModel@ObsCovs <- list()  
  semModel@ImpCovs <- list()
  semModel@Computed <- FALSE
  semModel@Original <- list(object)
  
  return(semModel)
}

plot.BECFA <- function(x,...){
  
  parEsts <- parTable(x$BECFA)
  Mode <- x$BECFA@external$runjags$summaries[,"Mode"]
  Model <- left_join(parEsts,data.frame(jlabel=names(Mode),mode=Mode,stringsAsFactors=FALSE),by="jlabel")
  semPlotModel <- semPlotModel_blavaanModel(Model)
  semPaths(semPlotModel, "est", "hide",intercepts=FALSE,...)
}