lca_confirmatory.R
In tidySEM: Tidy Structural Equation Modeling

## ----setup, include = FALSE---------------------------------------------------
library(yaml)
library(scales)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  results = "hide",
  tidy.opts = list(width.cutoff = 60), tidy = TRUE
)
options(scipen = 1, digits = 2)
eval_results <- FALSE
if(suppressWarnings(tryCatch({isTRUE(as.logical(readLines("pkgdown.txt")))}, error = function(e){FALSE}))){
  eval_results <- TRUE
  knitr::opts_chunk$set(
  results = "markup"
)
}
run_everything = suppressWarnings(tryCatch({isTRUE(as.logical(readLines("run_everything.txt")))}, error = function(e){FALSE}))

## ---- echo = TRUE, eval=TRUE--------------------------------------------------
# Load required packages
library(tidySEM) 
library(ggplot2)
# Load data
df <- zegwaard_carecompass[, c("burdened", "trapped", "negaffect", "loneliness")]

## ----echo = TRUE, eval=FALSE--------------------------------------------------
#  desc <- tidySEM::descriptives(df)
#  desc <- desc[, c("name", "n", "missing", "unique",
#                   "mean", "median", "sd", "min", "max",
#                   "skew_2se", "kurt_2se")]
#  desc

## ----tabdesc, echo = FALSE, eval=TRUE-----------------------------------------
desc <- tidySEM::descriptives(df)
desc <- desc[, c("name", "n", "missing", "unique",
                 "mean", "median", "sd", "min", "max",
                 "skew_2se", "kurt_2se")]
knitr::kable(desc, caption = "Descriptive statistics")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  df_plot <- df
#  names(df_plot) <- paste0("Value.", names(df_plot))
#  df_plot <- reshape(df_plot, varying = names(df_plot), direction = "long",
#                     timevar = "Variable")
#  ggplot(df_plot, aes(x = Value)) +
#    geom_density() +
#    facet_wrap(~Variable)+
#    theme_bw()

## ---- echo = FALSE, eval = run_everything-------------------------------------
#  # mcartest <- mice::mcar(df)
#  df_plot <- df
#  names(df_plot) <- paste0("Value.", names(df_plot))
#  df_plot <- reshape(df_plot, varying = names(df_plot), direction = "long",
#                     timevar = "Variable")
#  p <- ggplot(df_plot, aes(x = Value)) +
#    geom_density() +
#    facet_wrap(~Variable)+
#    theme_bw()
#  ggsave("plot_lpa_desc.png", p, device = "png", width = 100, height = 100, units = "mm")

## ----figdesc, echo = FALSE, eval = eval_results, out.width="80%"--------------
#  knitr::include_graphics("plot_lpa_desc.png")

## ----fitlca, eval = run_everything, echo = FALSE------------------------------
#  set.seed(123) # setting seed
#  res <- mx_profiles(data = df,
#                     classes = 1:5)
#  saveRDS(res, "res_lpa.RData")
#  fit <- table_fit(res)
#  write.csv(fit, "lpatabfit.csv", row.names = FALSE)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  set.seed(123)
#  res <- mx_profiles(data = df,
#                     classes = 1:5)

## ----eval = eval_results, echo = FALSE----------------------------------------
#  fit <- read.csv("lpatabfit.csv", stringsAsFactors = FALSE)
#  class(fit) <- c("tidy_fit", "data.frame")

## ----fit_table, include = TRUE, eval=F----------------------------------------
#  fit <- table_fit(res) # model fit table
#  fit[ , c("Name", "LL", "Parameters", "n",
#           "BIC", "Entropy",
#           "prob_min", "prob_max",
#           "n_min", "n_max",
#           "np_ratio", "np_local")]

## ----tabfit, echo = FALSE, eval = eval_results--------------------------------
#  tbl <- fit[ , c("Name", "LL", "Parameters", "n",
#           "BIC", "Entropy",
#           "prob_min", "prob_max",
#           "n_min", "n_max")]
#  names(tbl) <- c("Name", "LL", "p", "n",
#           "BIC", "Entropy",
#           "p_min", "p_max",
#           "n_min", "n_max")
#  knitr::kable(tbl, caption = "Model fit table")

## ----lmr_table, echo = TRUE, eval=FALSE---------------------------------------
#  lr_lmr(res)

## ----tablmr, echo = FALSE, eval = run_everything------------------------------
#  res <- readRDS("res_lpa.RData")
#  tbl <- lr_lmr(res)
#  write.csv(tbl, "lpatablmr.csv", row.names = FALSE)

## ----echo = FALSE, eval = TRUE------------------------------------------------
tbl <- read.csv("lpatablmr.csv", stringsAsFactors = FALSE)
knitr::kable(tbl, caption = "LMR test table", digits = 2)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  library(future)
#  library(progressr)
#  plan(multisession) # Parallel processing for Windows
#  handlers("progress") # Progress bar
#  set.seed(1)
#  res_blrt <- BLRT(res, replications = 20)

## ----eval = run_everything, echo = FALSE--------------------------------------
#  library(future)
#  library(progressr)
#  plan(multisession) # Parallel processing for Windows
#  handlers("progress") # Progress bar
#  set.seed(1)
#  res_blrt <- BLRT(res, replications = 20)
#  write.csv(res_blrt, "appendixbresblrt.csv", row.names = FALSE)

## ----eval = eval_results, echo = FALSE----------------------------------------
#  res_blrt <- read.csv("appendixbresblrt.csv", stringsAsFactors = FALSE)
#  knitr::kable(res_blrt, caption = "BLRT test table", digits = 2)

## ----eval = run_everything, echo = FALSE--------------------------------------
#  res_alt <- mx_profiles(df, classes = 4, variances = "varying")
#  compare <- list(res[[4]], res_alt)
#  fit_compare <- table_fit(compare)
#  write.csv(fit_compare, "lpa_fit_compare.csv", row.names = FALSE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  res_alt <- mx_profiles(df, classes = 4, variances = "varying")
#  compare <- list(res[[4]], res_alt)
#  table_fit(compare)

## ----tabfitcomp, echo = FALSE, eval = eval_results----------------------------
#  fit_compare <- read.csv("lpa_fit_compare.csv", stringsAsFactors = FALSE)
#  class(fit_compare) <- c("tidy_fit", "data.frame")
#  knitr::kable(fit_compare[ , c("Name", "LL", "Parameters",
#           "BIC", "Entropy",
#           "prob_min", "prob_max",
#           "n_min", "n_max")], caption = "Comparing competing theoretical models")

## ---- echo = TRUE, eval =FALSE------------------------------------------------
#  res_final <- mx_switch_labels(res[[4]])

## ---- echo = FALSE, eval = run_everything-------------------------------------
#  res_final <- mx_switch_labels(res[[4]])
#  cp <- class_prob(res_final)
#  out <- list(counts = cp$sum.posterior$proportion)

## ----echo = TRUE, eval = FALSE------------------------------------------------
#  table_results(res_final, columns = c("label", "est", "se", "confint", "class"))

## ----echo = FALSE, eval = run_everything--------------------------------------
#  tab <- table_results(res_final, columns = c("label", "est", "se", "confint", "class"))
#  write.csv(tab, "lpa_tab_res.csv", row.names = FALSE)

## ---- eval = eval_results, echo=FALSE-----------------------------------------
#  tab <- read.csv("lpa_tab_res.csv", stringsAsFactors = FALSE)
#  knitr::kable(tab, caption = "Four-class model results")

## ----echo = TRUE, eval = FALSE------------------------------------------------
#  plot_bivariate(res_final)

## ----echo = FALSE, eval = run_everything--------------------------------------
#  p <- plot_bivariate(res_final, return_list = TRUE)
#  p[[1]] <- p[[1]] + scale_y_continuous(breaks= c(0, .1, .2, .3), labels = c(".0", ".1", ".2", ".3"))
#  p <- tidySEM:::merge_corplots(p)
#  ggsave("lpa_bivariate.png", p, device = "png", width = 210, height = 100, units = "mm", scale = 1.5)

## ----echo = FALSE, eval = eval_results, fig.cap="Bivariate profile plot", out.width="80%"----
#  knitr::include_graphics("lpa_bivariate.png")

## ----echo = TRUE, eval = FALSE------------------------------------------------
#  plot_profiles(res_final)

## ----echo = FALSE, eval = run_everything--------------------------------------
#  p1 <- plot_profiles(res_final, ci = NULL, sd= FALSE, add_line = TRUE, rawdata = FALSE)
#  p1 <- p1 + theme(legend.position = c(.85, .2))
#  p2 <- plot_profiles(res_final)
#  p2 <- p2 + theme(legend.position = "none")
#  p <- ggpubr::ggarrange(p1, p2)
#  ggsave("lpa_profiles.png", p, device = "png", width = 210, height = 100, units = "mm", scale = 1)

## ----echo = FALSE, eval = eval_results, fig.cap="Bivariate profile plot", out.width="80%"----
#  knitr::include_graphics("lpa_profiles.png")

## ----echo = TRUE, eval=FALSE--------------------------------------------------
#  aux_sex <- BCH(res_final, data = zegwaard_carecompass$sexpatient)

## ----echo = FALSE, eval=run_everything----------------------------------------
#  aux_sex <- BCH(res_final, data = zegwaard_carecompass$sexpatient)
#  saveRDS(aux_sex, "lpa_aux_sex.RData")

## ----echo = FALSE, eval=FALSE-------------------------------------------------
#  aux_sex <- readRDS("lpa_aux_sex.RData")

## ----echo = TRUE, eval=FALSE--------------------------------------------------
#  df_aux <- zegwaard_carecompass[, c("freqvisit", "distance")]
#  df_aux$freqvisit <- as.numeric(df_aux$freqvisit)
#  aux_model <- BCH(res_final, model = "freqvisit ~ distance",
#                   data = df_aux)

## ----echo = FALSE, eval=run_everything----------------------------------------
#  df_aux <- zegwaard_carecompass[, c("freqvisit", "distance")]
#  df_aux$freqvisit <- as.numeric(df_aux$freqvisit)
#  aux_model <- BCH(res_final, model = "freqvisit ~ distance",
#                 data = df_aux)
#  saveRDS(aux_model, "lpa_aux_model.RData")

## ----echo = FALSE, eval=FALSE-------------------------------------------------
#  aux_model <- readRDS("lpa_aux_model.RData")

## ----echo = TRUE, eval = FALSE------------------------------------------------
#  df_new <- data.frame(
#    burdened = 2,
#    trapped = 0.5,
#    negaffect = 1.5,
#    loneliness = 4
#  )
#  predict(res_final, newdata = df_new)

## ----echo = FALSE, eval = TRUE------------------------------------------------
structure(c(0.000808074819286418, 1.35412723878036e-15, 0.999191879285034, 
4.58956785816116e-08, 3), dim = c(1L, 5L), dimnames = list(NULL, 
    c("class1", "class2", "class3", "class4", "predicted")))