ICSUS LCA

Install packages and load dataset.

## R markdown link: https://rpubs.com/jasohosk/1424635

## Run libraries
packages <- c(
  "haven",
  "tidyverse",
  "psych",
  "ggplot2",
  "tidySEM",
  "OpenMx",
  "gmodels",
  "mice",
  "nnet",
  "future",
  "car",
  "mlogit",
  "dfidx",
  "tidyLPA",
  "broom"
)
installed <- packages %in% rownames(installed.packages())
if (any(!installed)) {
  install.packages(packages[!installed])
}
lapply(packages, library, character.only = TRUE)

## Set working library
setwd(
  "/Users/jasonhoskin/Library/CloudStorage/OneDrive-IndianaUniversity/study-gambling"
)

## Import dataset
url <- "ICSUS_2025_State_18 to 25_CLN.sav"
df1 <- read_sav(url)

## Remove unnecessary variables
rm(installed, packages, url)

Run descriptive statistics and simplify dataset.

## Sample size and number of variables
dim(df1)

## [1] 5402  257

## Prepare demographics
age <- df1 |>
  mutate(age = case_when(
    age < 5 ~ 0,
    ## < 21
    age >= 5 ~ 1 ## 21+
  )) |>
  dplyr::select(autoid, age)
ethnicity <- df1 |>
  mutate(ethnicity = case_when(
    ethnic == 2 ~ 0,
    ## Not Hispanic/Latino
    ethnic == 1 ~ 1 ## Hispanic/Latino
  )) |>
  dplyr::select(autoid, ethnicity)
race <- df1 |>
  mutate(race = case_when(
    race == 1 ~ 0,
    ## White
    race == 2 ~ 1,
    ## Black/African American
    race == 3 ~ 2,
    ## Asian
    race == 4 ~ 4,
    ## Other
    race == 5 ~ 4,
    ## Other
    race == 6 ~ 3,
    ## Multiracial
    race == 7 ~ 4,
    ## Other
  )) |>
  dplyr::select(autoid, race)
gender <- df1 |>
  mutate(
    gender = case_when(
      gender3 == 1 ~ 0,
      # Female
      gender2 == 1 ~ 1,
      # Male
      gender1 == 1 ~ 2,
      # Other
      gender4 == 1 ~ 2,
      # Other
      gender5 == 1 ~ 2,
      # Other
      gender6 == 1 ~ 2,
      # Other
      gender7 == 1 ~ 2,
      # Other
      gender8 == 1 ~ 2 # Other
    )
  ) |>
  dplyr::select(autoid, gender)
greekmem <- df1 |>
  mutate(greekmem = case_when(
    greekmem == 2 ~ 0,
    # No Greek Membership
    greekmem == 1 ~ 1 # Greek Membership
  )) |>
  dplyr::select(autoid, greekmem)
studentstatus <- df1 |>
  mutate(studentstatus = case_when(
    stustat == 1 ~ 0 ## Full-Time
    ,
    stustat == 2 ~ 1 ## Part-Time
  )) |>
  dplyr::select(autoid, studentstatus)
residence <- df1 |>
  mutate(residence = case_when(
    residenc < 4 ~ 0,
    ## Off Campus
    residenc >= 4 ~ 1
  )) |> ## On Campus
  dplyr::select(autoid, residence)
## Combine dems
dfs_to_join <- list(age, ethnicity, race, gender, studentstatus, residence, greekmem)

# Use reduce to join them all by "uid"
combined_demographics <- dfs_to_join %>%
  reduce(full_join, by = "autoid")

## Simplify dataset
df2 <- df1 |>
  dplyr::select(autoid |
                  smokmo:othmo |
                  gmpools:gmoth |
                  gcqsleep:gcqdepress |
                  mhdays)

## Combine prepared dems to new dataset

dems_to_join <- list(combined_demographics, df2)

df2 <- dems_to_join %>%
  reduce(full_join, by = "autoid")

## New sample size and number of variables
dim(df2)

## [1] 5402   46

## Calculate frequencies
freq_list <- lapply(df2[, 2:45], table, useNA = "ifany")
freq_list

## $age
## 
##    0    1 
## 3087 2315 
## 
## $ethnicity
## 
##    0    1 <NA> 
## 4419  670  313 
## 
## $race
## 
##    0    1    2    3    4 <NA> 
## 4148  280  487  245  194   48 
## 
## $gender
## 
##    0    1    2 <NA> 
## 3318 1802  255   27 
## 
## $studentstatus
## 
##    0    1 <NA> 
## 5168  216   18 
## 
## $residence
## 
##    0    1 <NA> 
## 2602 2797    3 
## 
## $greekmem
## 
##    0    1 <NA> 
## 4667  727    8 
## 
## $smokmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 4106  869  238   86   37   33   18   12    3 
## 
## $cigarmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 4538  683  125   21    6    9    1    3   16 
## 
## $snufmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 4993  248   42   22   14   15   13   44   11 
## 
## $hookmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 5027  306   31   11    9    1    2    4   11 
## 
## $ecigmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 3445 1013  242  115   55   55   88  373   16 
## 
## $alcmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 1424 1258  978  708  412  282  115  161   64 
## 
## $marmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 3237 1049  321  170  126  124  161  202   12 
## 
## $cocmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 5275   97    9    6    2    2    1    1    9 
## 
## $hallumo
## 
##    1    2    3    4    5    6    7 <NA> 
## 5074  269   34    9    3    1    2   10 
## 
## $hermo
## 
##    1    2    3    4    5 <NA> 
## 5378   11    1    1    1   10 
## 
## $methmo
## 
##    1    2    3    4    5    6    8 <NA> 
## 5367   20    1    1    1    2    1    9 
## 
## $inhmo
## 
##    1    2    3    4    5    6    7 <NA> 
## 5271   94   13    7    4    2    2    9 
## 
## $rxstmo
## 
##    1    2    3    4    5    6    7    8 <NA> 
## 5152  182   33   11    4    2    5    2   11 
## 
## $rxpkmo
## 
##    1    2    3    4    5    6    7 <NA> 
## 5259  119   10    3    1    1    2    7 
## 
## $rxsedmo
## 
##    1    2    3    4    5    6    7 <NA> 
## 5286   92    7    6    2    1    2    6 
## 
## $othmo
## 
##    1    2    3    4    5    6 <NA> 
## 5272   92   13    6    2    3   14 
## 
## $gmpools
## 
##    1    2    3    4 <NA> 
## 3887  278   77   27 1133 
## 
## $gmfant
## 
##    1    2    3    4 <NA> 
## 4013  162   53   35 1139 
## 
## $gmvideo
## 
##    1    2    3    4 <NA> 
## 3782  342  112   36 1130 
## 
## $gmonsport
## 
##    1    2    3    4 <NA> 
## 3921  167  103   78 1133 
## 
## $gmothsp
## 
##    1    2    3    4 <NA> 
## 4136   64   33   32 1137 
## 
## $gmonline
## 
##    1    2    3    4 <NA> 
## 4076  127   45   21 1133 
## 
## $gmesports
## 
##    1    2    3    4 <NA> 
## 4180   42   24   22 1134 
## 
## $gmhorse
## 
##    1    2    3    4 <NA> 
## 4175   73    8    9 1137 
## 
## $gmcards
## 
##    1    2    3    4 <NA> 
## 3885  262   93   27 1135 
## 
## $gmlott
## 
##    1    2    3    4 <NA> 
## 3274  876  100   18 1134 
## 
## $gmcasino
## 
##    1    2    3    4 <NA> 
## 3930  297   30    7 1138 
## 
## $gmcharit
## 
##    1    2    3    4 <NA> 
## 3831  393   35    8 1135 
## 
## $gmoth
## 
##    1    2    3    4 <NA> 
## 4132   62   16   13 1179 
## 
## $gcqsleep
## 
##    1    2    3 <NA> 
## 1748   32    6 3616 
## 
## $gcqhyg
## 
##    1    2    3 <NA> 
## 1770   15    1 3616 
## 
## $gcqfrnd
## 
##    1    2    3 <NA> 
## 1771   15    1 3615 
## 
## $gcqfmly
## 
##    1    2    3 <NA> 
## 1760   18    7 3617 
## 
## $gcqacad
## 
##    1    2    3 <NA> 
## 1764   15    4 3619 
## 
## $gcqmoney
## 
##    1    2    3 <NA> 
## 1705   64   15 3618 
## 
## $gcqbad
## 
##    1    2    3 <NA> 
## 1650  123   10 3619 
## 
## $gcqdepress
## 
##    1    2    3 <NA> 
## 1743   28    3 3628

Assess NAs and remove all participants with no responses to the gambling behaviors items.

colSums(is.na(df2))

##        autoid           age     ethnicity          race        gender 
##             0             0           313            48            27 
## studentstatus     residence      greekmem        smokmo       cigarmo 
##            18             3             8             3            16 
##        snufmo        hookmo        ecigmo         alcmo         marmo 
##            11            11            16            64            12 
##         cocmo       hallumo         hermo        methmo         inhmo 
##             9            10            10             9             9 
##        rxstmo        rxpkmo       rxsedmo         othmo       gmpools 
##            11             7             6            14          1133 
##        gmfant       gmvideo     gmonsport       gmothsp      gmonline 
##          1139          1130          1133          1137          1133 
##     gmesports       gmhorse       gmcards        gmlott      gmcasino 
##          1134          1137          1135          1134          1138 
##      gmcharit         gmoth      gcqsleep        gcqhyg       gcqfrnd 
##          1135          1179          3616          3616          3615 
##       gcqfmly       gcqacad      gcqmoney        gcqbad    gcqdepress 
##          3617          3619          3618          3619          3628 
##        mhdays 
##          1086

df3 <- df2 |>
  filter(!if_all(
    c(
      gmpools,
      gmfant,
      gmvideo,
      gmonsport,
      gmothsp,
      gmonline,
      gmesports,
      gmhorse,
      gmcards,
      gmlott,
      gmcasino,
      gmcharit,
      gmoth
    ),
    is.na
  ))

## Number of participants removed from dataset:
dim(df2)[1] - dim(df3)[1]

## [1] 1124

colSums(is.na(df3))

##        autoid           age     ethnicity          race        gender 
##             0             0           235            38            24 
## studentstatus     residence      greekmem        smokmo       cigarmo 
##            11             2             5             2            14 
##        snufmo        hookmo        ecigmo         alcmo         marmo 
##             9            11            13            44             8 
##         cocmo       hallumo         hermo        methmo         inhmo 
##             6             8             9             7             6 
##        rxstmo        rxpkmo       rxsedmo         othmo       gmpools 
##            10             6             6             8             9 
##        gmfant       gmvideo     gmonsport       gmothsp      gmonline 
##            15             6             9            13             9 
##     gmesports       gmhorse       gmcards        gmlott      gmcasino 
##            10            13            11            10            14 
##      gmcharit         gmoth      gcqsleep        gcqhyg       gcqfrnd 
##            11            55          2492          2492          2491 
##       gcqfmly       gcqacad      gcqmoney        gcqbad    gcqdepress 
##          2493          2495          2494          2495          2504 
##        mhdays 
##            34

Discriminate skip logic NAs from missing data NAs.

## Create a gambling participation indicator
gm_vars <- names(df3)[25:37] # gmpools:gmoth

df3$any_gamble <- apply(df3[gm_vars], 1, function(x) {
  if (all(is.na(x)))
    return(NA)
  if (all(x == 1, na.rm = TRUE))
    return(0)  # never gambled
  return(1)  # gambled at least once
})

## Create a function to separate skip logic NAs from missing data NAs
classify_gcq_na <- function(data, gate_var, consequence_vars) {
  gate <- data[[gate_var]]
  
  results <- lapply(consequence_vars, function(var) {
    item <- data[[var]]
    
    skip_na <- sum(is.na(item) & gate == 0, na.rm = TRUE)
    true_na <- sum(is.na(item) & gate == 1, na.rm = TRUE)
    observed <- sum(!is.na(item))
    
    data.frame(
      variable = var,
      skip_logic_na = skip_na,
      true_missing_na = true_na,
      observed = observed
    )
  })
  
  do.call(rbind, results)
}

## Apply function to dataset
gcq_vars <- names(df3)[38:45]  # gcqsleep:gcqdepress

na_summary <- classify_gcq_na(data = df3,
                              gate_var = "any_gamble",
                              consequence_vars = gcq_vars)
na_summary

##     variable skip_logic_na true_missing_na observed
## 1   gcqsleep          2462              30     1786
## 2     gcqhyg          2462              30     1786
## 3    gcqfrnd          2462              29     1787
## 4    gcqfmly          2462              31     1785
## 5    gcqacad          2462              33     1783
## 6   gcqmoney          2462              32     1784
## 7     gcqbad          2462              33     1783
## 8 gcqdepress          2462              42     1774

Make the variables dichotomous and convert to class “factor”.

df4 <- df3 |>
  mutate(across(
    .cols  = starts_with("gm"),
    .fns   = ~ case_when(.x == 1 ~ 1, .x > 1 ~ 2, .default = NA),
    .names = "{.col}_dich"
  ))

df5 <- df4 |>
  mutate(across(
    .cols  = starts_with("gcq"),
    .fns   = ~ case_when(.x == 1 ~ 1, .x > 1 ~ 2, .default = NA),
    .names = "{.col}_dich"
  ))

df6 <- df5 |>
  mutate(across(
    .cols  = ends_with("mo"),
    .fns   = ~ case_when(.x == 1 ~ 1, .x > 1 ~ 2, .default = NA),
    .names = "{.col}_dich"
  ))

## Filter out unnecessary variables
df6 <- df6 |>
  dplyr::select(autoid:greekmem | gmpools_dich:othmo_dich | mhdays)

# Convert substances and gambling to factor with 1 as reference
cols_to_factor <- colnames(df6)[9:45]
df6[cols_to_factor] <- lapply(df6[cols_to_factor], function(x) {
  x <- as.factor(x)
  x <- relevel(x, ref = "1")
  return(x)
})

## Convert demographics to class "factor"
df6$age <- factor(df6$age,
                  levels = c(0, 1),
                  labels = c("<21 (ref)", "21+"))

df6$ethnicity <- factor(
  df6$ethnicity,
  levels = c(0, 1),
  labels = c("Not Hispanic/Latino (ref)", "Hispanic/Latino")
)

df6$race <- factor(
  df6$race,
  levels = c(0, 1, 2, 3, 4),
  labels = c(
    "White (ref)",
    "African American/Black",
    "Asian",
    "Multiracial",
    "Other"
  )
)

df6$gender <- factor(
  df6$gender,
  levels = c(0, 1, 2),
  labels = c("Female (ref)", "Male", "Other")
)

df6$studentstatus <- factor(
  df6$studentstatus,
  levels = c(0, 1),
  labels = c("Full Time (ref)", "Part Time")
)

df6$residence <- factor(
  df6$residence,
  levels = c(0, 1),
  labels = c("Off Campus (ref)", "On Campus")
)

df6$greekmem <- factor(df6$greekmem,
                       levels = c(0, 1),
                       labels = c("No Greek (ref)", "Greek"))

# Reference group verification
demo_vars <- c("age",
               "ethnicity",
               "race",
               "gender",
               "studentstatus",
               "residence",
               "greekmem")

for (v in demo_vars) {
  cat(sprintf(
    "%-16s | Reference: %-30s | Levels: %s\n",
    v,
    levels(df6[[v]])[1],
    paste(levels(df6[[v]]), collapse = " / ")
  ))
}

## age              | Reference: <21 (ref)                      | Levels: <21 (ref) / 21+
## ethnicity        | Reference: Not Hispanic/Latino (ref)      | Levels: Not Hispanic/Latino (ref) / Hispanic/Latino
## race             | Reference: White (ref)                    | Levels: White (ref) / African American/Black / Asian / Multiracial / Other
## gender           | Reference: Female (ref)                   | Levels: Female (ref) / Male / Other
## studentstatus    | Reference: Full Time (ref)                | Levels: Full Time (ref) / Part Time
## residence        | Reference: Off Campus (ref)               | Levels: Off Campus (ref) / On Campus
## greekmem         | Reference: No Greek (ref)                 | Levels: No Greek (ref) / Greek

lapply(df6[demo_vars], table, useNA = "ifany")

## $age
## 
## <21 (ref)       21+ 
##      2407      1871 
## 
## $ethnicity
## 
## Not Hispanic/Latino (ref)           Hispanic/Latino                      <NA> 
##                      3525                       518                       235 
## 
## $race
## 
##            White (ref) African American/Black                  Asian 
##                   3309                    222                    379 
##            Multiracial                  Other                   <NA> 
##                    190                    140                     38 
## 
## $gender
## 
## Female (ref)         Male        Other         <NA> 
##         2683         1346          225           24 
## 
## $studentstatus
## 
## Full Time (ref)       Part Time            <NA> 
##            4099             168              11 
## 
## $residence
## 
## Off Campus (ref)        On Campus             <NA> 
##             2045             2231                2 
## 
## $greekmem
## 
## No Greek (ref)          Greek           <NA> 
##           3741            532              5

## Create binary gambling consequences variable (any vs. no consequences)
gcq_dich_cols <- c(
  "gcqsleep_dich", "gcqhyg_dich", "gcqfrnd_dich", "gcqfmly_dich",
  "gcqacad_dich",  "gcqmoney_dich", "gcqbad_dich", "gcqdepress_dich"
)

df6 <- df6 |>
  mutate(gcq_any_dich = case_when(
    if_all(all_of(gcq_dich_cols), ~ is.na(.x))              ~ NA_integer_,
    if_any(all_of(gcq_dich_cols), ~ as.character(.x) == "2") ~ 2L,
    .default = 1L
  ))

df6$gcq_any_dich <- ordered(df6$gcq_any_dich, levels = c(1, 2))

table(df6$gcq_any_dich, useNA = "ifany")

## 
##    1    2 <NA> 
## 1602  185 2491

Conduct the first LCA using the tidySEM package.

set.seed(123)

df6_lca1 <- df6 |>
  dplyr::select(gmpools_dich:gmoth_dich) |>
  as.data.frame() |>
  mutate(across(gmpools_dich:gmoth_dich, as.ordered))

## Run first LCA
## lca1 <- tidySEM::mx_lca(data = df6_lca1,
## classes = 1:4,
## missing = "fiml")

## Save LCA for future knitting
## lca1 <- saveRDS(lca1, file = "lca1_results.rds")

## Load saved LCA
lca1 <- readRDS("lca1_results.rds")

Assess latent classes for first LCA.

fit1 <- table_fit(lca1)
fit1

##   Name Classes        LL    n Parameters      AIC      BIC    saBIC   Entropy
## 1    1       1 -13815.68 4278         13 27657.36 27740.06 27698.75 1.0000000
## 2    2       2 -12071.01 4278         27 24196.01 24367.77 24281.97 0.8538182
## 3    3       3 -11805.73 4278         41 23693.46 23954.27 23823.99 0.7098784
## 4    4       4 -11697.07 4278         55 23504.14 23854.00 23679.24 0.7785529
##    prob_min  prob_max      n_min     n_max  np_ratio   np_local
## 1 1.0000000 1.0000000 1.00000000 1.0000000 329.07692 329.076923
## 2 0.8703959 0.9827206 0.16432913 0.8356709 158.44444  54.076923
## 3 0.7706319 0.9206315 0.10402057 0.6280972 104.34146  34.230769
## 4 0.5839871 0.9831005 0.01496026 0.7515194  77.78182   4.923077

# Lo-Mendell-Rubin Adjusted Likelihood Ratio Test (LMR-ALRT)
## 1-class vs 2-class
calc_lrt(
  n = as.numeric(fit1[1, 4]),
  null_ll = as.numeric(fit1[1, 3]),
  null_param = as.numeric(fit1[1, 5]),
  null_classes = as.numeric(fit1[1, 2]),
  alt_ll = as.numeric(fit1[2, 3]),
  alt_param = as.numeric(fit1[2, 5]),
  alt_classes = as.numeric(fit1[2, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 3489.348, LMR LR (df = 14) = 3355.573, p < 0.001

## 1-class vs 3-class
calc_lrt(
  n = as.numeric(fit1[1, 4]),
  null_ll = as.numeric(fit1[1, 3]),
  null_param = as.numeric(fit1[1, 5]),
  null_classes = as.numeric(fit1[1, 2]),
  alt_ll = as.numeric(fit1[3, 3]),
  alt_param = as.numeric(fit1[3, 5]),
  alt_classes = as.numeric(fit1[3, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 4019.897, LMR LR (df = 28) = 3941.333, p < 0.001

## 1-class vs 4-class
calc_lrt(
  n = as.numeric(fit1[1, 4]),
  null_ll = as.numeric(fit1[1, 3]),
  null_param = as.numeric(fit1[1, 5]),
  null_classes = as.numeric(fit1[1, 2]),
  alt_ll = as.numeric(fit1[4, 3]),
  alt_param = as.numeric(fit1[4, 5]),
  alt_classes = as.numeric(fit1[4, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 4237.225, LMR LR (df = 42) = 4181.656, p < 0.001

## 2-class vs 3-class
calc_lrt(
  n = as.numeric(fit1[2, 4]),
  null_ll = as.numeric(fit1[2, 3]),
  null_param = as.numeric(fit1[2, 5]),
  null_classes = as.numeric(fit1[2, 2]),
  alt_ll = as.numeric(fit1[3, 3]),
  alt_param = as.numeric(fit1[3, 5]),
  alt_classes = as.numeric(fit1[3, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 530.549, LMR LR (df = 14) = 510.208, p < 0.001

## 2-class vs 4-class
calc_lrt(
  n = as.numeric(fit1[2, 4]),
  null_ll = as.numeric(fit1[2, 3]),
  null_param = as.numeric(fit1[2, 5]),
  null_classes = as.numeric(fit1[2, 2]),
  alt_ll = as.numeric(fit1[4, 3]),
  alt_param = as.numeric(fit1[4, 5]),
  alt_classes = as.numeric(fit1[4, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 747.877, LMR LR (df = 28) = 733.261, p < 0.001

## 3-class vs 4-class
calc_lrt(
  n = as.numeric(fit1[3, 4]),
  null_ll = as.numeric(fit1[3, 3]),
  null_param = as.numeric(fit1[3, 5]),
  null_classes = as.numeric(fit1[3, 2]),
  alt_ll = as.numeric(fit1[4, 3]),
  alt_param = as.numeric(fit1[4, 5]),
  alt_classes = as.numeric(fit1[4, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 217.328, LMR LR (df = 14) = 208.996, p < 0.001

# Visualize item probability profiles
plot_prob(lca1[[3]]) +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    size = 7
  ))

plot_prob(lca1[[4]]) +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    size = 7
  ))

# Minimum sample size for 3 classes
fit1[3, 4] * fit1[3, 12]

## [1] 445

# Minimum sample size for 4 classes
fit1[4, 4] * fit1[4, 12]

## [1] 64

# Get class proportions from the 3-class LCA2 model
probs <- class_prob(lca1[[3]], type = "sum.posterior")
probs

## $sum.posterior
##    class     count proportion
## 1 class1 1181.3268   0.276140
## 2 class2  467.9448   0.109384
## 3 class3 2628.7283   0.614476

Conduct the second latent class analysis.

set.seed(123)

## Refilter dataset to only include appropriate variables
gm_behavior_vars <- c(
  "gmpools_dich",
  "gmfant_dich",
  "gmvideo_dich",
  "gmonsport_dich",
  "gmothsp_dich",
  "gmonline_dich",
  "gmesports_dich",
  "gmhorse_dich",
  "gmcards_dich",
  "gmlott_dich",
  "gmcasino_dich",
  "gmcharit_dich",
  "gmoth_dich"
)

gcq_vars <- c("gcq_any_dich")

lca2_vars <- c(gm_behavior_vars, gcq_vars)

diagnostic_lca2 <- df6 |>
  mutate(
    any_behavior = if_any(all_of(gm_behavior_vars), ~ as.character(.x) == "2"),
    any_conseq   = if_any(all_of(gcq_vars), ~ as.character(.x) == "2"),
    current_rule = if_any(all_of(lca2_vars), ~ as.character(.x) == "2")
  )

table(diagnostic_lca2$any_behavior,
      diagnostic_lca2$any_conseq,
      useNA = "ifany")

##        
##         FALSE TRUE <NA>
##   FALSE     0    0 2381
##   TRUE   1602  185   29
##   <NA>      0    0   81

table(diagnostic_lca2$any_behavior,
      diagnostic_lca2$current_rule,
      useNA = "ifany")

##        
##         TRUE <NA>
##   FALSE    0 2381
##   TRUE  1816    0
##   <NA>     0   81

df6_lca2 <- diagnostic_lca2 |>
  filter(any_behavior) |>
  dplyr::select(autoid, all_of(lca2_vars)) |>
  mutate(across(all_of(lca2_vars), ~ ordered(.x, levels = c(1, 2))))

## lca2 <- tidySEM::mx_lca(data = df6_lca2[2:15],
##                          classes = 1:4,
##                          missing = "fiml")

## Save LCA for future knitting
## lca2 <- saveRDS(lca2, file = "lca2_results.rds")

## Run second LCA
lca2 <- readRDS("lca2_results.rds")

Assess latent classes for second LCA.

fit2 <- table_fit(lca2)
fit2

##   Name Classes        LL    n Parameters      AIC      BIC    saBIC   Entropy
## 1    1       1 -11354.87 1816         21 22751.74 22867.34 22800.62 1.0000000
## 2    2       2 -10549.84 1816         43 21185.68 21422.37 21285.76 0.8525082
## 3    3       3 -10346.39 1816         65 20822.79 21180.57 20974.07 0.8818748
## 4    4       4 -10115.17 1816         87 20404.34 20883.23 20606.83 0.7868395
##    prob_min  prob_max      n_min     n_max np_ratio  np_local
## 1 1.0000000 1.0000000 1.00000000 1.0000000 86.47619 86.476190
## 2 0.8752273 0.9790906 0.18116740 0.8188326 42.23256 15.666667
## 3 0.8765232 0.9788917 0.06442731 0.7637665 27.93846  5.571429
## 4 0.8066107 0.9354766 0.01762115 0.5952643 20.87356  1.523810

# Lo-Mendell-Rubin Adjusted Likelihood Ratio Test (LMR-ALRT)
## 1-class vs 2-class
calc_lrt(
  n = as.numeric(fit2[1, 4]),
  null_ll = as.numeric(fit2[1, 3]),
  null_param = as.numeric(fit2[1, 5]),
  null_classes = as.numeric(fit2[1, 2]),
  alt_ll = as.numeric(fit2[2, 3]),
  alt_param = as.numeric(fit2[2, 5]),
  alt_classes = as.numeric(fit2[2, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 1610.067, LMR LR (df = 22) = 1541.592, p < 0.001

## 2-class vs 3-class
calc_lrt(
  n = as.numeric(fit2[2, 4]),
  null_ll = as.numeric(fit2[2, 3]),
  null_param = as.numeric(fit2[2, 5]),
  null_classes = as.numeric(fit2[2, 2]),
  alt_ll = as.numeric(fit2[3, 3]),
  alt_param = as.numeric(fit2[3, 5]),
  alt_classes = as.numeric(fit2[3, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 406.890, LMR LR (df = 22) = 389.586, p < 0.001

## 3-class vs 4-class
calc_lrt(
  n = as.numeric(fit2[3, 4]),
  null_ll = as.numeric(fit2[3, 3]),
  null_param = as.numeric(fit2[3, 5]),
  null_classes = as.numeric(fit2[3, 2]),
  alt_ll = as.numeric(fit2[4, 3]),
  alt_param = as.numeric(fit2[4, 5]),
  alt_classes = as.numeric(fit2[4, 2])
)

## Lo-Mendell-Rubin ad-hoc adjusted likelihood ratio rest:
## 
## LR = 462.442, LMR LR (df = 22) = 442.775, p < 0.001

# Visualize item probability profiles
plot_prob(lca2[[3]]) +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    size = 7
  ))

plot_prob(lca2[[4]]) +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    size = 7
  ))

# Minimum sample size for 3 classes
fit2[3, 4] * fit2[3, 12]

## [1] 117

# Minimum sample size for 4 classes
fit2[4, 4] * fit2[4, 12]

## [1] 32

# Get class proportions from the 3-class LCA2 model
probs <- class_prob(lca2[[4]], type = "sum.posterior")
probs

## $sum.posterior
##    class      count proportion
## 1 class1   33.97541 0.01870893
## 2 class2  330.13765 0.18179386
## 3 class3  406.84396 0.22403302
## 4 class4 1045.04298 0.57546420

Conduct multi-start LCA to locate global maxima.

# 20 runs per model with different random seeds
n_runs <- 20

# Pre-generate seeds for full reproducibility
set.seed(3817)
seeds_m1 <- sample(1:99999, n_runs)
seeds_m2 <- sample(1:99999, n_runs)

cat("Model 1 seeds:", seeds_m1, "\n")

## Model 1 seeds: 29402 51485 42324 21712 33101 15891 17118 83342 37078 12567 75927 64547 34252 28848 5843 26407 37628 4639 52898 3359

cat("Model 2 seeds:", seeds_m2, "\n")

## Model 2 seeds: 82904 48667 5422 31219 98765 19818 45640 13290 60001 12316 9610 63358 96461 75508 94495 65501 52021 41393 78085 62407

# Summarise a single multi-start run into a one-row-per-k data frame
summarise_run <- function(lca_run, seed, run_idx) {
  fit <- table_fit(lca_run)
  n_k <- nrow(fit)

  props_str <- vapply(seq_len(n_k), function(k) {
    cp <- tryCatch(
      class_prob(lca_run[[k]], type = "sum.posterior"),
      error = function(e) NULL
    )
    if (is.null(cp)) return(NA_character_)
    pct <- round(cp$sum.posterior$proportion * 100, 1)
    paste0(pct, "%", collapse = " / ")
  }, character(1))

  data.frame(
    run         = run_idx,
    seed        = seed,
    k           = as.integer(fit[[2]]),
    LL          = as.numeric(fit[[3]]),
    AIC         = as.numeric(fit$AIC),
    BIC         = as.numeric(fit$BIC),
    class_props = props_str,
    stringsAsFactors = FALSE
  )
}

# ---- Model 1 ----
if (file.exists("multistart_lca1_raw.rds") && file.exists("multistart_lca1_summary.rds")) {
  cat("\n=== Model 1: Loading saved results ===\n")
  raw_runs_m1 <- readRDS("multistart_lca1_raw.rds")
  summary_m1  <- readRDS("multistart_lca1_summary.rds")
} else {
  cat("\n=== Model 1: 20 runs ===\n")
  results_m1  <- vector("list", n_runs)
  raw_runs_m1 <- vector("list", n_runs)

  for (i in seq_len(n_runs)) {
    cat(sprintf("  Run %02d/%d  seed=%d  ...", i, n_runs, seeds_m1[i]))
    set.seed(seeds_m1[i])
    lca_run <- tryCatch(
      tidySEM::mx_lca(data = df6_lca1, classes = 1:4, missing = "fiml"),
      error = function(e) { cat(" ERROR:", conditionMessage(e)); NULL }
    )
    if (!is.null(lca_run)) {
      raw_runs_m1[[i]] <- lca_run
      results_m1[[i]]  <- summarise_run(lca_run, seeds_m1[i], i)
      cat(" done\n")
    }
  }

  summary_m1 <- do.call(rbind, Filter(Negate(is.null), results_m1))
  saveRDS(raw_runs_m1, "multistart_lca1_raw.rds")
  saveRDS(summary_m1,  "multistart_lca1_summary.rds")
}

## 
## === Model 1: 20 runs ===
##   Run 01/20  seed=29402  ...

##  done
##   Run 02/20  seed=51485  ...

##  done
##   Run 03/20  seed=42324  ...

##  done
##   Run 04/20  seed=21712  ...

##  done
##   Run 05/20  seed=33101  ...

##  done
##   Run 06/20  seed=15891  ...

##  done
##   Run 07/20  seed=17118  ...

##  done
##   Run 08/20  seed=83342  ...

##  done
##   Run 09/20  seed=37078  ...

##  done
##   Run 10/20  seed=12567  ...

##  done
##   Run 11/20  seed=75927  ...

##  done
##   Run 12/20  seed=64547  ...

##  done
##   Run 13/20  seed=34252  ...

##  done
##   Run 14/20  seed=28848  ...

##  done
##   Run 15/20  seed=5843  ...

##  done
##   Run 16/20  seed=26407  ...

##  done
##   Run 17/20  seed=37628  ...

##  done
##   Run 18/20  seed=4639  ...

##  done
##   Run 19/20  seed=52898  ...

##  done
##   Run 20/20  seed=3359  ...

##  done

# ---- Model 2 ----
if (file.exists("multistart_lca2_raw.rds") && file.exists("multistart_lca2_summary.rds")) {
  cat("\n=== Model 2: Loading saved results ===\n")
  raw_runs_m2 <- readRDS("multistart_lca2_raw.rds")
  summary_m2  <- readRDS("multistart_lca2_summary.rds")
} else {
  cat("\n=== Model 2: 20 runs ===\n")
  results_m2  <- vector("list", n_runs)
  raw_runs_m2 <- vector("list", n_runs)

  for (i in seq_len(n_runs)) {
    cat(sprintf("  Run %02d/%d  seed=%d  ...", i, n_runs, seeds_m2[i]))
    set.seed(seeds_m2[i])
    lca_run <- tryCatch(
      tidySEM::mx_lca(data = df6_lca2[2:15], classes = 1:4, missing = "fiml"),
      error = function(e) { cat(" ERROR:", conditionMessage(e)); NULL }
    )
    if (!is.null(lca_run)) {
      raw_runs_m2[[i]] <- lca_run
      results_m2[[i]]  <- summarise_run(lca_run, seeds_m2[i], i)
      cat(" done\n")
    }
  }

  summary_m2 <- do.call(rbind, Filter(Negate(is.null), results_m2))
  saveRDS(raw_runs_m2, "multistart_lca2_raw.rds")
  saveRDS(summary_m2,  "multistart_lca2_summary.rds")
}

## 
## === Model 2: 20 runs ===
##   Run 01/20  seed=82904  ...

##  done
##   Run 02/20  seed=48667  ...

##  done
##   Run 03/20  seed=5422  ...

##  done
##   Run 04/20  seed=31219  ...

##  done
##   Run 05/20  seed=98765  ...

##  done
##   Run 06/20  seed=19818  ...

##  done
##   Run 07/20  seed=45640  ...

##  done
##   Run 08/20  seed=13290  ...

##  done
##   Run 09/20  seed=60001  ...

##  done
##   Run 10/20  seed=12316  ...

##  done
##   Run 11/20  seed=9610  ...

##  done
##   Run 12/20  seed=63358  ...

##  done
##   Run 13/20  seed=96461  ...

##  done
##   Run 14/20  seed=75508  ...

##  done
##   Run 15/20  seed=94495  ...

##  done
##   Run 16/20  seed=65501  ...

##  done
##   Run 17/20  seed=52021  ...

##  done
##   Run 18/20  seed=41393  ...

##  done
##   Run 19/20  seed=78085  ...

##  done
##   Run 20/20  seed=62407  ...

##  done

# ---- Results: sorted by k, then best LL ----
cat("\n===== MODEL 1 RESULTS =====\n")

## 
## ===== MODEL 1 RESULTS =====

summary_m1 |> dplyr::arrange(k, dplyr::desc(LL))

##    run  seed k        LL      AIC      BIC                  class_props
## 1   18  4639 1 -13815.68 27657.36 27740.06                         100%
## 2    8 83342 1 -13815.68 27657.36 27740.06                         100%
## 3   11 75927 1 -13815.68 27657.36 27740.06                         100%
## 4    7 17118 1 -13815.68 27657.36 27740.06                         100%
## 5   12 64547 1 -13815.68 27657.36 27740.06                         100%
## 6    1 29402 1 -13815.68 27657.36 27740.06                         100%
## 7   16 26407 1 -13815.68 27657.36 27740.06                         100%
## 8    4 21712 1 -13815.68 27657.36 27740.06                         100%
## 9   13 34252 1 -13815.68 27657.36 27740.06                         100%
## 10   6 15891 1 -13815.68 27657.36 27740.06                         100%
## 11  19 52898 1 -13815.68 27657.36 27740.06                         100%
## 12   9 37078 1 -13815.68 27657.36 27740.06                         100%
## 13  17 37628 1 -13815.68 27657.36 27740.06                         100%
## 14  20  3359 1 -13815.68 27657.36 27740.06                         100%
## 15  15  5843 1 -13815.68 27657.36 27740.06                         100%
## 16   5 33101 1 -13815.68 27657.36 27740.06                         100%
## 17   3 42324 1 -13815.68 27657.36 27740.06                         100%
## 18  14 28848 1 -13815.68 27657.36 27740.06                         100%
## 19   2 51485 1 -13815.68 27657.36 27740.06                         100%
## 20  10 12567 1 -13815.68 27657.36 27740.06                         100%
## 21  15  5843 2 -12050.73 24155.46 24327.21                84.6% / 15.4%
## 22   3 42324 2 -12050.73 24155.46 24327.21                84.6% / 15.4%
## 23   2 51485 2 -12050.73 24155.46 24327.21                15.4% / 84.6%
## 24  14 28848 2 -12050.73 24155.46 24327.21                15.4% / 84.6%
## 25  12 64547 2 -12050.73 24155.46 24327.21                84.6% / 15.4%
## 26  13 34252 2 -12050.73 24155.46 24327.21                15.4% / 84.6%
## 27  10 12567 2 -12050.73 24155.46 24327.21                84.6% / 15.4%
## 28  18  4639 2 -12055.86 24165.72 24337.48                16.2% / 83.8%
## 29   4 21712 2 -12055.86 24165.72 24337.48                16.2% / 83.8%
## 30  19 52898 2 -12055.86 24165.72 24337.48                16.2% / 83.8%
## 31   1 29402 2 -12055.86 24165.73 24337.48                84.1% / 15.9%
## 32   6 15891 2 -12055.86 24165.73 24337.48                15.9% / 84.1%
## 33   9 37078 2 -12071.01 24196.01 24367.77                17.2% / 82.8%
## 34   5 33101 2 -12073.31 24200.61 24372.37                82.2% / 17.8%
## 35   7 17118 2 -12080.02 24214.04 24385.79                18.4% / 81.6%
## 36  20  3359 2 -12080.08 24214.17 24385.92                82.2% / 17.8%
## 37  11 75927 2 -12084.59 24223.19 24394.94                82.3% / 17.7%
## 38  16 26407 2 -12092.58 24239.15 24410.91                19.2% / 80.8%
## 39  17 37628 2 -12168.72 24391.44 24563.19                23.1% / 76.9%
## 40   8 83342 2 -12423.93 24901.85 25073.61                36.9% / 63.1%
## 41  15  5843 3 -11758.02 23598.04 23858.85         7.6% / 24.9% / 67.5%
## 42  20  3359 3 -11767.91 23617.83 23878.64           8% / 30.4% / 61.6%
## 43   3 42324 3 -11768.26 23618.52 23879.33         62.7% / 30.3% / 7.1%
## 44   6 15891 3 -11769.38 23620.76 23881.57           61% / 31.2% / 7.8%
## 45  12 64547 3 -11769.47 23620.94 23881.75         31.2% / 7.4% / 61.5%
## 46   2 51485 3 -11769.95 23621.90 23882.71         6.7% / 63.5% / 29.8%
## 47  18  4639 3 -11770.54 23623.07 23883.88         62.9% / 6.8% / 30.4%
## 48   5 33101 3 -11770.79 23623.58 23884.40         62.6% / 30.6% / 6.8%
## 49   8 83342 3 -11776.81 23635.61 23896.42         56.5% / 8.4% / 35.1%
## 50   1 29402 3 -11777.28 23636.55 23897.37         8.4% / 35.6% / 55.9%
## 51   7 17118 3 -11778.99 23639.98 23900.79         8.5% / 54.2% / 37.3%
## 52  16 26407 3 -11783.60 23649.19 23910.00         8.3% / 38.8% / 52.9%
## 53  13 34252 3 -11788.62 23659.23 23920.04           21.2% / 69% / 9.8%
## 54  10 12567 3 -11796.01 23674.01 23934.82           22% / 69.4% / 8.7%
## 55  17 37628 3 -11809.49 23700.99 23961.80        26.6% / 10.7% / 62.7%
## 56  19 52898 3 -11809.91 23701.81 23962.62         53.3% / 37.1% / 9.6%
## 57  11 75927 3 -11821.99 23725.97 23986.78           4.8% / 29.2% / 66%
## 58  14 28848 3 -11911.55 23905.09 24165.90              13% / 75% / 12%
## 59   9 37078 3 -11923.58 23929.16 24189.97         24.6% / 5.2% / 70.2%
## 60   4 21712 3 -12065.18 24212.36 24473.18           17.3% / 0% / 82.7%
## 61   2 51485 4 -11652.11 23414.22 23764.09  63.7% / 1.4% / 7.2% / 27.8%
## 62  20  3359 4 -11688.65 23487.31 23837.18      63% / 4.6% / 28.4% / 4%
## 63  11 75927 4 -11694.46 23498.93 23848.80    8% / 2.4% / 21.4% / 68.2%
## 64  13 34252 4 -11699.36 23508.71 23858.58  67.1% / 2.2% / 7.5% / 23.2%
## 65  12 64547 4 -11700.19 23510.38 23860.25  2.6% / 21.6% / 68.1% / 7.7%
## 66  17 37628 4 -11703.63 23517.27 23867.13  7.8% / 8.5% / 19.6% / 64.1%
## 67   7 17118 4 -11704.17 23518.33 23868.20  22.3% / 67.4% / 4.8% / 5.4%
## 68  16 26407 4 -11704.36 23518.72 23868.59    2.2% / 7.3% / 64% / 26.5%
## 69   8 83342 4 -11709.42 23528.83 23878.70  63.1% / 7.8% / 21.3% / 7.8%
## 70  19 52898 4 -11714.90 23539.80 23889.67 22.5% / 10.2% / 59.5% / 7.8%
## 71   4 21712 4 -11715.96 23541.93 23891.80  31.8% / 8.1% / 2.4% / 57.7%
## 72   3 42324 4 -11717.50 23545.00 23894.87  7.2% / 1.5% / 63.8% / 27.5%
## 73  15  5843 4 -11717.56 23545.13 23895.00          8% / 4% / 66% / 22%
## 74  18  4639 4 -11722.55 23555.10 23904.97 52.5% / 15.3% / 1.3% / 30.9%
## 75   5 33101 4 -11726.16 23562.32 23912.19  6.9% / 67.6% / 17.1% / 8.4%
## 76   6 15891 4 -11728.26 23566.51 23916.38  8.2% / 54.6% / 7.3% / 29.9%
## 77   9 37078 4 -11730.84 23571.67 23921.54    1.5% / 70% / 9.9% / 18.6%
## 78  10 12567 4 -11734.32 23578.64 23928.51  27.3% / 7.6% / 60.2% / 4.9%
## 79   1 29402 4 -11743.79 23597.58 23947.45 14.5% / 63.2% / 7.4% / 14.9%
## 80  14 28848 4 -11751.46 23612.92 23962.78    20.2% / 61.4% / 9% / 9.4%

cat("\n===== MODEL 2 RESULTS =====\n")

## 
## ===== MODEL 2 RESULTS =====

summary_m2 |> dplyr::arrange(k, dplyr::desc(LL))

##    run  seed k         LL      AIC      BIC                   class_props
## 1   10 12316 1 -10393.342 20814.68 20891.75                          100%
## 2   13 96461 1 -10393.342 20814.68 20891.75                          100%
## 3    3  5422 1 -10393.342 20814.68 20891.75                          100%
## 4   16 65501 1 -10393.342 20814.68 20891.75                          100%
## 5    2 48667 1 -10393.342 20814.68 20891.75                          100%
## 6    4 31219 1 -10393.342 20814.68 20891.75                          100%
## 7   20 62407 1 -10393.342 20814.68 20891.75                          100%
## 8    5 98765 1 -10393.342 20814.68 20891.75                          100%
## 9    1 82904 1 -10393.342 20814.68 20891.75                          100%
## 10  12 63358 1 -10393.342 20814.68 20891.75                          100%
## 11  19 78085 1 -10393.342 20814.68 20891.75                          100%
## 12  17 52021 1 -10393.342 20814.68 20891.75                          100%
## 13  18 41393 1 -10393.342 20814.68 20891.75                          100%
## 14   6 19818 1 -10393.342 20814.68 20891.75                          100%
## 15   8 13290 1 -10393.342 20814.68 20891.75                          100%
## 16  15 94495 1 -10393.342 20814.68 20891.75                          100%
## 17   9 60001 1 -10393.342 20814.68 20891.75                          100%
## 18  14 75508 1 -10393.342 20814.68 20891.75                          100%
## 19  11  9610 1 -10393.342 20814.68 20891.75                          100%
## 20   7 45640 1 -10393.342 20814.68 20891.75                          100%
## 21   5 98765 2  -9721.796 19501.59 19661.22                 82.5% / 17.5%
## 22   3  5422 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 23  20 62407 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 24  12 63358 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 25   2 48667 2  -9721.796 19501.59 19661.22                 82.5% / 17.5%
## 26  11  9610 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 27  10 12316 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 28  14 75508 2  -9721.796 19501.59 19661.22                 82.5% / 17.5%
## 29  19 78085 2  -9721.796 19501.59 19661.22                 82.5% / 17.5%
## 30   4 31219 2  -9721.796 19501.59 19661.22                 82.5% / 17.5%
## 31  13 96461 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 32   6 19818 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 33   1 82904 2  -9721.796 19501.59 19661.22                 17.5% / 82.5%
## 34  17 52021 2  -9735.595 19529.19 19688.82                 79.9% / 20.1%
## 35   8 13290 2  -9793.943 19645.89 19805.51                 27.4% / 72.6%
## 36  18 41393 2  -9795.191 19648.38 19808.01                 71.9% / 28.1%
## 37   7 45640 2  -9825.554 19709.11 19868.74                 73.6% / 26.4%
## 38  15 94495 2  -9828.369 19714.74 19874.37                 26.5% / 73.5%
## 39   9 60001 2  -9873.683 19805.37 19964.99                 65.5% / 34.5%
## 40  16 65501 2  -9875.289 19808.58 19968.21                     33% / 67%
## 41   9 60001 3  -9577.875 19243.75 19485.94         13.2% / 17.6% / 69.2%
## 42  20 62407 3  -9581.671 19251.34 19493.54         12.5% / 69.8% / 17.6%
## 43  19 78085 3  -9581.986 19251.97 19494.17         17.7% / 64.3% / 18.1%
## 44   7 45640 3  -9583.223 19254.45 19496.64           12.3% / 70% / 17.7%
## 45  11  9610 3  -9584.820 19257.64 19499.83         16.9% / 18.4% / 64.6%
## 46  13 96461 3  -9590.146 19268.29 19510.49           19% / 19.6% / 61.4%
## 47  17 52021 3  -9590.146 19268.29 19510.49           61.4% / 19.6% / 19%
## 48  16 65501 3  -9593.420 19274.84 19517.03          17.8% / 3.7% / 78.5%
## 49  14 75508 3  -9599.803 19287.61 19529.80         33.5% / 48.7% / 17.9%
## 50   1 82904 3  -9602.648 19293.30 19535.49            4.5% / 19% / 76.6%
## 51   3  5422 3  -9606.091 19300.18 19542.37         17.3% / 28.2% / 54.6%
## 52   5 98765 3  -9606.091 19300.18 19542.37         17.3% / 54.6% / 28.2%
## 53   4 31219 3  -9606.091 19300.18 19542.37         17.3% / 28.2% / 54.6%
## 54   8 13290 3  -9606.091 19300.18 19542.37         17.3% / 28.2% / 54.6%
## 55   6 19818 3  -9606.147 19300.29 19542.49           28% / 17.2% / 54.8%
## 56  18 41393 3  -9606.455 19300.91 19543.10          21.6% / 1.8% / 76.6%
## 57  12 63358 3  -9608.334 19304.67 19546.86         12.5% / 22.3% / 65.2%
## 58  15 94495 3  -9622.315 19332.63 19574.82         20.2% / 51.6% / 28.1%
## 59   2 48667 3  -9628.516 19345.03 19587.23         15.8% / 20.7% / 63.5%
## 60  10 12316 3  -9657.001 19402.00 19644.19         27.9% / 59.9% / 12.2%
## 61   9 60001 4  -9445.345 19008.69 19333.45  4.1% / 63.7% / 18.7% / 13.5%
## 62  15 94495 4  -9445.345 19008.69 19333.45  13.5% / 18.7% / 63.7% / 4.1%
## 63   4 31219 4  -9445.345 19008.69 19333.45  4.1% / 18.7% / 13.5% / 63.7%
## 64   6 19818 4  -9449.821 19017.64 19342.40      18.6% / 65.4% / 4% / 12%
## 65   2 48667 4  -9467.452 19052.90 19377.66    19.3% / 60.9% / 13.8% / 6%
## 66  20 62407 4  -9476.018 19070.04 19394.79    48.9% / 3.8% / 28.2% / 19%
## 67   5 98765 4  -9488.861 19095.72 19420.48    23.9% / 52% / 3.6% / 20.5%
## 68  13 96461 4  -9491.341 19100.68 19425.44 51.6% / 21.2% / 13.7% / 13.5%
## 69  16 65501 4  -9491.341 19100.68 19425.44 21.2% / 51.6% / 13.7% / 13.5%
## 70   8 13290 4  -9492.442 19102.88 19427.64    13.6% / 10.1% / 68.3% / 8%
## 71   1 82904 4  -9493.985 19105.97 19430.73 12.9% / 21.3% / 52.5% / 13.4%
## 72   7 45640 4  -9495.632 19109.26 19434.02   13% / 20.6% / 52.8% / 13.5%
## 73  18 41393 4  -9503.698 19125.40 19450.16    11.8% / 12.2% / 8.1% / 68%
## 74  10 12316 4  -9505.682 19129.36 19454.12 13.8% / 15.3% / 55.8% / 15.1%
## 75  14 75508 4  -9518.022 19154.04 19478.80 41.5% / 24.7% / 19.1% / 14.7%
## 76   3  5422 4  -9532.395 19182.79 19507.55   10.2% / 10% / 19.3% / 60.5%
## 77  11  9610 4  -9540.726 19199.45 19524.21 16.4% / 31.5% / 31.2% / 20.8%
## 78  12 63358 4  -9554.309 19226.62 19551.38 40.8% / 25.8% / 16.8% / 16.5%
## 79  17 52021 4  -9554.977 19227.95 19552.71  16.4% / 10.7% / 67.6% / 5.2%
## 80  19 78085 4  -9571.413 19260.83 19585.58  2.5% / 67.2% / 17.8% / 12.6%

write_csv(summary_m1, file = "summary_m1.csv")
write_csv(summary_m2, file = "summary_m2.csv")

Export LCA results.

## Export LCAs
save_path <- "/Users/jasonhoskin/Library/CloudStorage/OneDrive-IndianaUniversity/study-gambling"

lca_theme <- theme(axis.text.x = element_text(
  angle = 45,
  hjust = 1,
  size = 7
))

ggsave(
  "lca1_3class.png",
  plot = plot_prob(lca1[[3]]) + lca_theme,
  width = 10,
  height = 6,
  dpi = 300,
  path = save_path
)

ggsave(
  "lca2_3class.png",
  plot = plot_prob(lca2[[3]]) + lca_theme,
  width = 10,
  height = 6,
  dpi = 300,
  path = save_path
)

Apply latent class assignments to dataframes.

## Assign latent classes to dataframe
# Extract most likely class assignment for each participant
lca1_classes <- class_prob(lca1[[3]], type = "individual")$individual[, "predicted"]
lca2_classes <- class_prob(lca2[[3]], type = "individual")$individual[, "predicted"]

# Add LCA 1 classes to gamblers-only subset
df6$lca1_class <- lca1_classes

# Add LCA 2 classes to gamblers-only subset
df6_lca2$lca2_class <- lca2_classes

# Then merge back to full dataset by row identifier
df6 <- df6 |>
  left_join(df6_lca2 |> dplyr::select(autoid, lca2_class), by = "autoid")

## Verify that it worked
table(df6$lca1_class, useNA = "always")

## 
##    1    2    3 <NA> 
## 1146  445 2687    0

# Should show counts for classes 1, 2, 3 with no NAs

table(df6$lca2_class, useNA = "always")

## 
##    1    2    3 <NA> 
##  312  117 1387 2462

# Should show counts for classes 1, 2, 3 with 2462 NAs

Prepare dataframes for regression modeling.

## Filter out any unnecessary variables
df6_model1 <- df6 |>
  dplyr::select(autoid:greekmem | smokmo_dich:lca1_class | mhdays)

## Filter out non-gamblers to only include participants in 2nd LCA
df6_model2 <- df6 |>
  filter(!is.na(lca2_class)) |>
  dplyr::select(autoid:greekmem |
                  smokmo_dich:othmo_dich | lca2_class | mhdays)

## Make DVs class "factor"
df6_model1 <- df6_model1 |>
  mutate(lca1_class = as.factor(lca1_class))

df6_model2 <- df6_model2 |>
  mutate(lca2_class = as.factor(lca2_class))

### 4. Relevel DVs so no gambling is the reference group
df6_model1$lca1_class <- relevel(df6_model1$lca1_class, ref = 3)

df6_model2$lca2_class <- relevel(df6_model2$lca2_class, ref = 3)

Treat missing data using MICE and evaluate the imputations.

plan(multisession, workers = parallel::detectCores() - 1)

# Strip all haven labels across the whole data frame
df6_model1 <- haven::zap_labels(df6_model1)
df6_model1 <- haven::zap_label(df6_model1)   # variable labels
df6_model1 <- haven::zap_labels(df6_model1)  # value labels

# Strip all haven labels across the whole data frame
df6_model2 <- haven::zap_labels(df6_model2)
df6_model2 <- haven::zap_label(df6_model2)   # variable labels
df6_model2 <- haven::zap_labels(df6_model2)  # value labels

# Define variables
outcome1    <- "lca1_class"
predictors <- c(
  "age",
  "ethnicity",
  "race",
  "gender",
  "studentstatus",
  "residence",
  "greekmem",
  "smokmo_dich",
  "cigarmo_dich",
  "snufmo_dich",
  "hookmo_dich",
  "ecigmo_dich",
  "alcmo_dich",
  "marmo_dich",
  "cocmo_dich",
  "hallumo_dich",
  "hermo_dich",
  "methmo_dich",
  "inhmo_dich",
  "rxstmo_dich",
  "rxpkmo_dich",
  "rxsedmo_dich",
  "othmo_dich",
  "mhdays"
)
all_vars1 <- c(outcome1, predictors)

# Establish imputation method
meth1 <- make.method(df6_model1[, all_vars1])
meth1[] <- "cart"

# Trimmed predictor matrix
pred1 <- quickpred(df6_model1[, all_vars1], mincor = 0.1, minpuc = 0.1)

# Impute in parallel
imp1 <- futuremice(
  df6_model1[, all_vars1],
  method = meth1,
  predictorMatrix = pred1,
  m = 20,
  maxit = 10,
  parallelseed = 42
)

## Evaluate imputations
plot(imp1)

stripplot(imp1, pch = 20, cex = 1.2)

# Define variables for model 2
outcome2    <- "lca2_class"
all_vars2 <- c(outcome2, predictors)

# Create imputation method for method 2
meth2 <- make.method(df6_model2[, all_vars2])
meth2[] <- "cart"

# Trimmed predictor matrix
pred2 <- quickpred(df6_model2[, all_vars2], mincor = 0.1, minpuc = 0.1)

# Impute in parallel
imp2 <- futuremice(
  df6_model2[, all_vars2],
  method = meth2,
  predictorMatrix = pred2,
  m = 20,
  maxit = 25,
  parallelseed = 42
)

## Evaluate imputations
plot(imp2)

stripplot(imp2, pch = 20, cex = 1.2)

Run multinomial logistic regression models.

# Fit model 1a: null model
fit1a <- with(imp1, multinom(lca1_class ~ 1, trace = FALSE))

# Pool model 1a
pooled1a <- pool(fit1a)
summary(pooled1a)

##          term   estimate  std.error statistic       df      p.value
## 1 (Intercept) -0.8521477 0.03528119 -24.15303 4274.001 6.31298e-121
## 2 (Intercept) -1.7981063 0.05117962 -35.13325 4274.001 8.75161e-238

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled1a)$term,
  fmi  = pooled1a$pooled$fmi,
  lambda = pooled1a$pooled$lambda
)

##          term          fmi lambda
## 1 (Intercept) 0.0004676173      0
## 2 (Intercept) 0.0004676173      0

# Fit model 1b: covariate model
fit1b <- with(
  imp1,
  multinom(
    lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
    trace = FALSE
  )
)

# Pool model 1b
pooled1b <- pool(fit1b)
summary(pooled1b)

##                          term    estimate  std.error   statistic       df
## 1                 (Intercept) -0.90589370 0.06895316 -13.1378134 4222.197
## 2                      age21+  0.13306035 0.07609206   1.7486760 4250.994
## 3                  genderMale  0.18582202 0.08006103   2.3210046 4225.336
## 4                 genderOther -0.24475314 0.16314979  -1.5001744 4210.950
## 5  raceAfrican American/Black -0.69905149 0.18305327  -3.8188419 4225.060
## 6                   raceAsian -0.99866581 0.15136043  -6.5979317 4132.906
## 7             raceMultiracial -0.11295103 0.16959683  -0.6659973 4178.481
## 8                   raceOther -0.60972873 0.23403602  -2.6052773 3729.321
## 9    ethnicityHispanic/Latino -0.09066452 0.11770694  -0.7702564 1673.410
## 10     studentstatusPart Time  0.15355311 0.17843027   0.8605777 4249.292
## 11         residenceOn Campus  0.13080162 0.07842130   1.6679348 4247.587
## 12              greekmemGreek  0.22880358 0.11132463   2.0552827 4240.974
## 13                (Intercept) -2.65619141 0.12052671 -22.0381977 4241.826
## 14                     age21+  0.40536820 0.11374212   3.5639234 4247.007
## 15                 genderMale  1.53422031 0.11085642  13.8397064 4243.662
## 16                genderOther -0.86215574 0.42535208  -2.0269226 4251.041
## 17 raceAfrican American/Black -0.38196541 0.24817418  -1.5391022 4192.615
## 18                  raceAsian -1.43758928 0.24934764  -5.7654017 4227.241
## 19            raceMultiracial -0.70892306 0.32876012  -2.1563536 4242.650
## 20                  raceOther -0.51982702 0.35942112  -1.4462896 4091.741
## 21   ethnicityHispanic/Latino -0.35837174 0.19004298  -1.8857405 2760.591
## 22     studentstatusPart Time -0.03624918 0.28284965  -0.1281571 3957.880
## 23         residenceOn Campus  0.22305855 0.12004676   1.8580972 4247.832
## 24              greekmemGreek  0.69567410 0.14634098   4.7537887 4246.278
##          p.value
## 1   1.135958e-38
## 2   8.041922e-02
## 3   2.033395e-02
## 4   1.336442e-01
## 5   1.360027e-04
## 6   4.696323e-11
## 7   5.054496e-01
## 8   9.216466e-03
## 9   4.412565e-01
## 10  3.895192e-01
## 11  9.540238e-02
## 12  3.991281e-02
## 13 5.370301e-102
## 14  3.693699e-04
## 15  1.223502e-42
## 16  4.273262e-02
## 17  1.238548e-01
## 18  8.727182e-09
## 19  3.111209e-02
## 20  1.481726e-01
## 21  5.943484e-02
## 22  8.980312e-01
## 23  6.322435e-02
## 24  2.062493e-06

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled1b)$term,
  fmi  = pooled1b$pooled$fmi,
  lambda = pooled1b$pooled$lambda
)

##                          term          fmi       lambda
## 1                 (Intercept) 0.0042763522 0.0038048025
## 2                      age21+ 0.0006954888 0.0002254484
## 3                  genderMale 0.0039934525 0.0035221192
## 4                 genderOther 0.0051989819 0.0047266111
## 5  raceAfrican American/Black 0.0040188480 0.0035474959
## 6                   raceAsian 0.0097932627 0.0093141967
## 7             raceMultiracial 0.0073731378 0.0068981377
## 8                   raceOther 0.0233669097 0.0228432909
## 9    ethnicityHispanic/Latino 0.0816817745 0.0805848888
## 10     studentstatusPart Time 0.0010357771 0.0005657086
## 11         residenceOn Campus 0.0013394518 0.0008693376
## 12              greekmemGreek 0.0023099407 0.0018395512
## 13                (Intercept) 0.0021981255 0.0017277778
## 14                     age21+ 0.0014362132 0.0009660804
## 15                 genderMale 0.0019459690 0.0014757058
## 16                genderOther 0.0006856007 0.0002155608
## 17 raceAfrican American/Black 0.0064932222 0.0060194034
## 18                  raceAsian 0.0038150296 0.0033438243
## 19            raceMultiracial 0.0020869910 0.0016166821
## 20                  raceOther 0.0116444383 0.0111614585
## 21   ethnicityHispanic/Latino 0.0475351126 0.0468453183
## 22     studentstatusPart Time 0.0165896313 0.0160928189
## 23         residenceOn Campus 0.0012977249 0.0008276182
## 24              greekmemGreek 0.0015540090 0.0010838510

# Fit model 1c: covariates + alcohol + marijuana + electronic vapor products + smoking + cigars
fit1c <- with(
  imp1,
  nnet::multinom(
    lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
    trace = FALSE
  )
)

# Pool model 1c
pooled1c <- pool(fit1c)
summary(pooled1c)

##                          term    estimate  std.error    statistic       df
## 1                 (Intercept) -1.09154408 0.08989229 -12.14280030 4218.020
## 2                      age21+  0.02011902 0.08136534   0.24726775 4238.161
## 3                  genderMale  0.16396216 0.08358490   1.96162414 4197.964
## 4                 genderOther -0.25366598 0.16448201  -1.54221103 4192.297
## 5  raceAfrican American/Black -0.62903913 0.18428640  -3.41337784 4217.341
## 6                   raceAsian -0.87261799 0.15350840  -5.68449656 4127.941
## 7             raceMultiracial -0.10969586 0.17055078  -0.64318591 4175.164
## 8                   raceOther -0.59016747 0.23558213  -2.50514535 3700.793
## 9    ethnicityHispanic/Latino -0.08377707 0.11825246  -0.70845941 1727.078
## 10     studentstatusPart Time  0.15065425 0.17938153   0.83985376 4237.961
## 11         residenceOn Campus  0.12152907 0.07931347   1.53226270 4238.151
## 12              greekmemGreek  0.13550127 0.11353846   1.19343939 4230.245
## 13                alcmo_dich2  0.13055184 0.09674665   1.34941967 4051.547
## 14                marmo_dich2  0.02249661 0.09435374   0.23842842 4217.679
## 15               ecigmo_dich2  0.26832379 0.10238481   2.62073827 4225.562
## 16               smokmo_dich2  0.06647428 0.11237047   0.59156356 4240.042
## 17              cigarmo_dich2  0.35381440 0.12484946   2.83392810 4214.268
## 18                (Intercept) -3.06352545 0.16267745 -18.83190017 4112.636
## 19                     age21+  0.11132027 0.12303231   0.90480521 4232.897
## 20                 genderMale  1.39008094 0.11932325  11.64970713 4226.966
## 21                genderOther -0.91393693 0.42802923  -2.13522081 4238.844
## 22 raceAfrican American/Black -0.20185349 0.25570957  -0.78938572 4181.855
## 23                  raceAsian -1.07772902 0.25609393  -4.20833492 4210.264
## 24            raceMultiracial -0.65374507 0.33271511  -1.96487943 4233.022
## 25                  raceOther -0.50693987 0.37071110  -1.36747962 4098.519
## 26   ethnicityHispanic/Latino -0.33344840 0.19557964  -1.70492387 2642.131
## 27     studentstatusPart Time -0.01395012 0.28929991  -0.04822027 3868.349
## 28         residenceOn Campus  0.17664522 0.12434734   1.42057902 4235.915
## 29              greekmemGreek  0.37281996 0.15549838   2.39758104 4236.964
## 30                alcmo_dich2  0.20147127 0.16420064   1.22698224 3793.614
## 31                marmo_dich2  0.21243362 0.14218671   1.49404687 4231.948
## 32               ecigmo_dich2  0.43949356 0.15407969   2.85237833 4231.190
## 33               smokmo_dich2  0.20409079 0.15788542   1.29265129 4238.424
## 34              cigarmo_dich2  0.91068264 0.15338533   5.93722130 4234.560
##         p.value
## 1  2.245069e-33
## 2  8.047130e-01
## 3  4.987225e-02
## 4  1.230978e-01
## 5  6.476743e-04
## 6  1.402563e-08
## 7  5.201388e-01
## 8  1.228275e-02
## 9  4.787556e-01
## 10 4.010378e-01
## 11 1.255323e-01
## 12 2.327643e-01
## 13 1.772776e-01
## 14 8.115604e-01
## 15 8.805379e-03
## 16 5.541744e-01
## 17 4.619838e-03
## 18 5.964358e-76
## 19 3.656201e-01
## 20 6.795938e-31
## 21 3.280015e-02
## 22 4.299313e-01
## 23 2.626218e-05
## 24 4.949350e-02
## 25 1.715500e-01
## 26 8.832618e-02
## 27 9.615432e-01
## 28 1.555128e-01
## 29 1.654676e-02
## 30 2.199054e-01
## 31 1.352379e-01
## 32 4.360388e-03
## 33 1.962022e-01
## 34 3.131136e-09

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled1c)$term,
  fmi  = pooled1c$pooled$fmi,
  lambda = pooled1c$pooled$lambda
)

##                          term          fmi       lambda
## 1                 (Intercept) 0.0037492849 0.0032770185
## 2                      age21+ 0.0012435431 0.0007723383
## 3                  genderMale 0.0054376743 0.0049639564
## 4                 genderOther 0.0058497010 0.0053755394
## 5  raceAfrican American/Black 0.0038145353 0.0033422238
## 6                   raceAsian 0.0095717584 0.0090920091
## 7             raceMultiracial 0.0069839693 0.0065084056
## 8                   raceOther 0.0239299904 0.0234026406
## 9    ethnicityHispanic/Latino 0.0794323276 0.0783669036
## 10     studentstatusPart Time 0.0012784717 0.0008072612
## 11         residenceOn Campus 0.0012453741 0.0007741690
## 12              greekmemGreek 0.0024089342 0.0019373987
## 13                alcmo_dich2 0.0129059248 0.0124187772
## 14                marmo_dich2 0.0037821189 0.0033098299
## 15               ecigmo_dich2 0.0029665253 0.0024947312
## 16               smokmo_dich2 0.0008933556 0.0004221946
## 17              cigarmo_dich2 0.0041008830 0.0036283631
## 18                (Intercept) 0.0103027270 0.0098215482
## 19                     age21+ 0.0020578822 0.0015864763
## 20                 genderMale 0.0028063544 0.0023346412
## 21                genderOther 0.0011212891 0.0006501025
## 22 raceAfrican American/Black 0.0065583159 0.0060833094
## 23                  raceAsian 0.0044548514 0.0039820502
## 24            raceMultiracial 0.0020405103 0.0015691100
## 25                  raceOther 0.0109443159 0.0104617931
## 26   ethnicityHispanic/Latino 0.0504932405 0.0497747692
## 27     studentstatusPart Time 0.0191497039 0.0186427192
## 28         residenceOn Campus 0.0016148245 0.0011435452
## 29              greekmemGreek 0.0014468288 0.0009755868
## 30                alcmo_dich2 0.0213507027 0.0208348931
## 31                marmo_dich2 0.0021870322 0.0017155816
## 32               ecigmo_dich2 0.0022872570 0.0018157693
## 33               smokmo_dich2 0.0011970616 0.0007258641
## 34              cigarmo_dich2 0.0018203883 0.0013490553

# Fit model 1d: model c + rx stimulants + rx painkillers + rx sedatives + mhdays
fit1d <- with(
  imp1,
  nnet::multinom(
    lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
    trace = FALSE
  )
)

# Pool model 1d
pooled1d <- pool(fit1d)
summary(pooled1d)

##                          term     estimate   std.error    statistic       df
## 1                 (Intercept) -1.148945052 0.097721667 -11.75732149 4209.515
## 2                      age21+  0.020530723 0.081515972   0.25186135 4230.182
## 3                  genderMale  0.177819171 0.084351817   2.10806570 4194.805
## 4                 genderOther -0.293691002 0.166425336  -1.76470127 4183.700
## 5  raceAfrican American/Black -0.612081521 0.184500297  -3.31750967 4210.346
## 6                   raceAsian -0.875214766 0.153753540  -5.69232272 4117.890
## 7             raceMultiracial -0.108520354 0.170826869  -0.63526514 4167.518
## 8                   raceOther -0.575159715 0.235922732  -2.43791563 3688.826
## 9    ethnicityHispanic/Latino -0.089359049 0.118421566  -0.75458426 1763.913
## 10     studentstatusPart Time  0.166567013 0.179623233   0.92731330 4229.316
## 11         residenceOn Campus  0.117083211 0.079419660   1.47423460 4230.306
## 12              greekmemGreek  0.138256260 0.114123612   1.21146061 4222.191
## 13                alcmo_dich2  0.130103892 0.096870262   1.34307361 4040.249
## 14                marmo_dich2 -0.006249400 0.095462174  -0.06546467 4202.611
## 15               ecigmo_dich2  0.256384099 0.102743644   2.49537674 4214.765
## 16               smokmo_dich2  0.048971772 0.113110964   0.43295337 4231.389
## 17              cigarmo_dich2  0.341582120 0.125922778   2.71263171 4204.409
## 18               rxstmo_dich2  0.385136493 0.217146744   1.77362316 4165.283
## 19               rxpkmo_dich2  0.443037537 0.269833991   1.64188928 4232.371
## 20              rxsedmo_dich2 -0.476287419 0.321305439  -1.48235094 4230.021
## 21                     mhdays  0.006811424 0.004274753   1.59340742 4126.659
## 22                (Intercept) -2.866103961 0.171602238 -16.70201971 4127.506
## 23                     age21+  0.082238018 0.123791441   0.66432717 4224.301
## 24                 genderMale  1.324559262 0.120845754  10.96074312 4219.110
## 25                genderOther -0.813988469 0.430625302  -1.89024766 4230.822
## 26 raceAfrican American/Black -0.209299486 0.256163247  -0.81705509 4169.629
## 27                  raceAsian -1.103763760 0.257448472  -4.28731913 4201.251
## 28            raceMultiracial -0.622804678 0.332755645  -1.87165774 4224.998
## 29                  raceOther -0.445765742 0.369938674  -1.20497200 4087.471
## 30   ethnicityHispanic/Latino -0.355214940 0.196531201  -1.80742263 2640.938
## 31     studentstatusPart Time  0.021964334 0.290288402   0.07566384 3855.467
## 32         residenceOn Campus  0.193992021 0.124965210   1.55236822 4228.368
## 33              greekmemGreek  0.317825043 0.157235560   2.02133056 4227.704
## 34                alcmo_dich2  0.209037118 0.164483874   1.27086694 3845.761
## 35                marmo_dich2  0.219968449 0.144541033   1.52184086 4222.280
## 36               ecigmo_dich2  0.440227608 0.154942912   2.84122456 4223.836
## 37               smokmo_dich2  0.183682084 0.159708129   1.15011105 4229.597
## 38              cigarmo_dich2  0.870838020 0.154654327   5.63086749 4224.778
## 39               rxstmo_dich2  0.620928646 0.265148111   2.34181810 4219.977
## 40               rxpkmo_dich2  0.331613291 0.363789158   0.91155353 4211.242
## 41              rxsedmo_dich2 -0.540195757 0.429776539  -1.25692240 4230.902
## 42                     mhdays -0.024503872 0.007591798  -3.22767690 4183.457
##         p.value
## 1  1.998041e-31
## 2  8.011605e-01
## 3  3.508447e-02
## 4  7.768694e-02
## 5  9.159464e-04
## 6  1.340445e-08
## 7  5.252904e-01
## 8  1.481913e-02
## 9  4.505992e-01
## 10 3.538168e-01
## 11 1.404928e-01
## 12 2.257867e-01
## 13 1.793236e-01
## 14 9.478071e-01
## 15 1.262030e-02
## 16 6.650707e-01
## 17 6.702304e-03
## 18 7.619855e-02
## 19 1.006872e-01
## 20 1.383214e-01
## 21 1.111454e-01
## 22 1.191630e-60
## 23 5.065172e-01
## 24 1.385268e-27
## 25 5.879313e-02
## 26 4.139436e-01
## 27 1.848998e-05
## 28 6.132301e-02
## 29 2.282839e-01
## 30 7.081015e-02
## 31 9.396905e-01
## 32 1.206490e-01
## 33 4.330845e-02
## 34 2.038529e-01
## 35 1.281238e-01
## 36 4.515593e-03
## 37 2.501632e-01
## 38 1.909430e-08
## 39 1.923633e-02
## 40 3.620560e-01
## 41 2.088511e-01
## 42 1.257627e-03

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled1d)$term,
  fmi  = pooled1d$pooled$fmi,
  lambda = pooled1d$pooled$lambda
)

##                          term          fmi       lambda
## 1                 (Intercept) 0.0038050563 0.0033318625
## 2                      age21+ 0.0012422950 0.0007702010
## 3                  genderMale 0.0050766429 0.0046023962
## 4                 genderOther 0.0059030730 0.0054279628
## 5  raceAfrican American/Black 0.0037247745 0.0032516359
## 6                   raceAsian 0.0096906308 0.0092097686
## 7             raceMultiracial 0.0069751777 0.0064987375
## 8                   raceOther 0.0240862790 0.0235573035
## 9    ethnicityHispanic/Latino 0.0779234002 0.0768785026
## 10     studentstatusPart Time 0.0013909272 0.0009188070
## 11         residenceOn Campus 0.0012203145 0.0007482240
## 12              greekmemGreek 0.0024203853 0.0019479559
## 13                alcmo_dich2 0.0130607641 0.0125723312
## 14                marmo_dich2 0.0044333124 0.0039596402
## 15               ecigmo_dich2 0.0032768713 0.0028040161
## 16               smokmo_dich2 0.0010217474 0.0005496839
## 17              cigarmo_dich2 0.0042758206 0.0038022761
## 18               rxstmo_dich2 0.0071137155 0.0066370862
## 19               rxpkmo_dich2 0.0008288107 0.0003567655
## 20              rxsedmo_dich2 0.0012704648 0.0007983662
## 21                     mhdays 0.0092551049 0.0087750531
## 22                (Intercept) 0.0092122352 0.0087322612
## 23                     age21+ 0.0021435298 0.0016712053
## 24                 genderMale 0.0027949284 0.0023223317
## 25                genderOther 0.0011275084 0.0006554316
## 26 raceAfrican American/Black 0.0068424704 0.0063662076
## 27                  raceAsian 0.0045499649 0.0040761949
## 28            raceMultiracial 0.0020477496 0.0015754577
## 29                  raceOther 0.0111005435 0.0106167933
## 30   ethnicityHispanic/Latino 0.0504394294 0.0497205931
## 31     studentstatusPart Time 0.0193361487 0.0188275671
## 32         residenceOn Campus 0.0015461473 0.0010739946
## 33              greekmemGreek 0.0016507126 0.0011785352
## 34                alcmo_dich2 0.0196291573 0.0191194450
## 35                marmo_dich2 0.0024089518 0.0019365271
## 36               ecigmo_dich2 0.0022061273 0.0017337805
## 37               smokmo_dich2 0.0013435393 0.0008714279
## 38              cigarmo_dich2 0.0020782390 0.0016059370
## 39               rxstmo_dich2 0.0026925757 0.0022200276
## 40               rxpkmo_dich2 0.0036369353 0.0031638556
## 41              rxsedmo_dich2 0.0011127069 0.0006406320
## 42                     mhdays 0.0059202287 0.0054450990

# Fit model 2a: null model
fit2a <- with(imp2, multinom(lca2_class ~ 1, trace = FALSE))

# Pool model 2a
pooled2a <- pool(fit2a)
summary(pooled2a)

##          term  estimate  std.error statistic       df       p.value
## 1 (Intercept) -1.491889 0.06265856 -23.80981 1812.003 2.997536e-109
## 2 (Intercept) -2.472710 0.09626986 -25.68520 1812.003 2.462857e-124

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled2a)$term,
  fmi  = pooled2a$pooled$fmi,
  lambda = pooled2a$pooled$lambda
)

##          term         fmi lambda
## 1 (Intercept) 0.001101926      0
## 2 (Intercept) 0.001101926      0

# Fit model 2b: covariate model
fit2b <- with(
  imp2,
  multinom(
    lca2_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
    trace = FALSE
  )
)

# Pool model 2b
pooled2b <- pool(fit2b)
summary(pooled2b)

##                          term    estimate std.error    statistic       df
## 1                 (Intercept) -2.76273929 0.1657924 -16.66385144 1787.848
## 2                      age21+  0.38169973 0.1412967   2.70140590 1788.907
## 3                  genderMale  1.66724094 0.1429523  11.66291550 1789.085
## 4                 genderOther -0.31554357 0.5303322  -0.59499234 1789.374
## 5  raceAfrican American/Black  0.43000215 0.3103759   1.38542376 1756.446
## 6                   raceAsian -0.59085706 0.3296026  -1.79263446 1749.851
## 7             raceMultiracial -0.61816958 0.4217739  -1.46564225 1787.481
## 8                   raceOther  0.23928498 0.4463803   0.53605634 1743.758
## 9    ethnicityHispanic/Latino -0.13584229 0.2304524  -0.58945928 1518.932
## 10     studentstatusPart Time -0.38798574 0.3700409  -1.04849413 1745.102
## 11         residenceOn Campus  0.28196878 0.1507771   1.87010384 1789.163
## 12              greekmemGreek  0.63126771 0.1771202   3.56406457 1789.719
## 13                (Intercept) -3.18055071 0.2201537 -14.44695666 1780.429
## 14                     age21+  0.16448092 0.2069757   0.79468713 1788.948
## 15                 genderMale  0.69293428 0.2080955   3.32988633 1788.833
## 16                genderOther  1.28412015 0.3445001   3.72748892 1786.920
## 17 raceAfrican American/Black  0.59784323 0.4509119   1.32585381 1762.921
## 18                  raceAsian  1.19429165 0.2872053   4.15832091 1782.215
## 19            raceMultiracial  1.08562995 0.3503753   3.09847721 1780.983
## 20                  raceOther -0.04641665 0.7651108  -0.06066658 1737.787
## 21   ethnicityHispanic/Latino  0.14339073 0.3162269   0.45344258 1279.555
## 22     studentstatusPart Time  0.18123821 0.4532748   0.39984179 1789.458
## 23         residenceOn Campus  0.02850135 0.2146861   0.13275823 1787.322
## 24              greekmemGreek -0.05754110 0.3182681  -0.18079443 1737.947
##         p.value
## 1  4.567605e-58
## 2  6.969879e-03
## 3  2.397174e-30
## 4  5.519238e-01
## 5  1.660988e-01
## 6  7.320407e-02
## 7  1.429217e-01
## 8  5.919880e-01
## 9  5.556409e-01
## 10 2.945562e-01
## 11 6.163258e-02
## 12 3.747316e-04
## 13 8.088838e-45
## 14 4.269009e-01
## 15 8.864743e-04
## 16 1.994174e-04
## 17 1.850599e-01
## 18 3.358235e-05
## 19 1.975577e-03
## 20 9.516317e-01
## 21 6.503070e-01
## 22 6.893208e-01
## 23 8.943995e-01
## 24 8.565500e-01

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled2b)$term,
  fmi  = pooled2b$pooled$fmi,
  lambda = pooled2b$pooled$lambda
)

##                          term         fmi       lambda
## 1                 (Intercept) 0.002207549 0.0010919793
## 2                      age21+ 0.001696057 0.0005805758
## 3                  genderMale 0.001606051 0.0004905803
## 4                 genderOther 0.001456299 0.0003408410
## 5  raceAfrican American/Black 0.010981383 0.0098558645
## 6                   raceAsian 0.012264851 0.0111365593
## 7             raceMultiracial 0.002375567 0.0012599560
## 8                   raceOther 0.013375341 0.0122443823
## 9    ethnicityHispanic/Latino 0.039616682 0.0383529634
## 10     studentstatusPart Time 0.013135734 0.0120053712
## 11         residenceOn Campus 0.001565897 0.0004504292
## 12              greekmemGreek 0.001272234 0.0001567853
## 13                (Intercept) 0.005033754 0.0039167113
## 14                     age21+ 0.001675696 0.0005602168
## 15                 genderMale 0.001733140 0.0006176544
## 16                genderOther 0.002624359 0.0015086766
## 17 raceAfrican American/Black 0.009616325 0.0084933903
## 18                  raceAsian 0.004436810 0.0033202171
## 19            raceMultiracial 0.004852932 0.0037360334
## 20                  raceOther 0.014406291 0.0132726340
## 21   ethnicityHispanic/Latino 0.061138935 0.0596726003
## 22     studentstatusPart Time 0.001411913 0.0002964570
## 23         residenceOn Campus 0.002447044 0.0013314139
## 24              greekmemGreek 0.014379385 0.0132458010

# Fit model 2c: covariates + alcohol + marijuana + vaping + smoking + cigars
fit2c <- with(
  imp2,
  multinom(
    lca2_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
    trace = FALSE
  )
)

# Pool model 2c
pooled2c <- pool(fit2c)
summary(pooled2c)

##                          term    estimate std.error    statistic       df
## 1                 (Intercept) -2.99170865 0.2180998 -13.71715559 1755.336
## 2                      age21+  0.18823308 0.1531333   1.22921088 1775.814
## 3                  genderMale  1.53668781 0.1521037  10.10289823 1777.695
## 4                 genderOther -0.30542061 0.5342014  -0.57173307 1779.012
## 5  raceAfrican American/Black  0.59742198 0.3185121   1.87566504 1745.619
## 6                   raceAsian -0.30086057 0.3385103  -0.88877816 1727.256
## 7             raceMultiracial -0.59926141 0.4257510  -1.40753956 1777.120
## 8                   raceOther  0.22459372 0.4516864   0.49723372 1728.935
## 9    ethnicityHispanic/Latino -0.09950893 0.2344836  -0.42437482 1503.581
## 10     studentstatusPart Time -0.38058537 0.3779161  -1.00706315 1732.655
## 11         residenceOn Campus  0.20978793 0.1549535   1.35387648 1778.580
## 12              greekmemGreek  0.42243084 0.1869261   2.25988128 1779.673
## 13                alcmo_dich2  0.07435978 0.2148771   0.34605731 1713.300
## 14                marmo_dich2 -0.10355684 0.1756919  -0.58942307 1776.826
## 15               ecigmo_dich2  0.52117445 0.1904485   2.73656307 1774.811
## 16               smokmo_dich2  0.01596974 0.1917674   0.08327665 1776.895
## 17              cigarmo_dich2  0.67815226 0.1810900   3.74483533 1770.996
## 18                (Intercept) -3.00242964 0.2615740 -11.47831967 1769.610
## 19                     age21+  0.29481923 0.2240508   1.31585856 1779.378
## 20                 genderMale  0.71654658 0.2159980   3.31737614 1778.200
## 21                genderOther  1.30066332 0.3480394   3.73711521 1777.387
## 22 raceAfrican American/Black  0.53760501 0.4546578   1.18243884 1750.757
## 23                  raceAsian  1.09789730 0.2973273   3.69255447 1770.163
## 24            raceMultiracial  1.07114234 0.3514848   3.04747799 1770.088
## 25                  raceOther -0.06684505 0.7678959  -0.08704963 1727.260
## 26   ethnicityHispanic/Latino  0.11708531 0.3172681   0.36904222 1287.160
## 27     studentstatusPart Time  0.20100135 0.4537840   0.44294500 1779.521
## 28         residenceOn Campus  0.01463019 0.2150364   0.06803585 1777.718
## 29              greekmemGreek  0.03863924 0.3240909   0.11922348 1727.588
## 30                alcmo_dich2 -0.30224805 0.2631828  -1.14843409 1776.151
## 31                marmo_dich2 -0.10641517 0.2661881  -0.39977426 1777.283
## 32               ecigmo_dich2  0.25618261 0.2863814   0.89455043 1759.767
## 33               smokmo_dich2 -0.21601402 0.3152955  -0.68511616 1776.821
## 34              cigarmo_dich2 -0.14919143 0.3344989  -0.44601470 1765.218
##         p.value
## 1  9.352440e-41
## 2  2.191556e-01
## 3  2.261683e-23
## 4  5.675751e-01
## 5  6.086807e-02
## 6  3.742461e-01
## 7  1.594423e-01
## 8  6.190875e-01
## 9  6.713532e-01
## 10 3.140450e-01
## 11 1.759478e-01
## 12 2.394906e-02
## 13 7.293421e-01
## 14 5.556524e-01
## 15 6.270357e-03
## 16 9.336410e-01
## 17 1.862835e-04
## 18 1.818721e-29
## 19 1.883910e-01
## 20 9.269938e-04
## 21 1.920270e-04
## 22 2.371922e-01
## 23 2.287446e-04
## 24 2.341838e-03
## 25 9.306422e-01
## 26 7.121569e-01
## 27 6.578594e-01
## 28 9.457647e-01
## 29 9.051122e-01
## 30 2.509441e-01
## 31 6.893708e-01
## 32 3.711498e-01
## 33 4.933600e-01
## 34 6.556413e-01

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled2c)$term,
  fmi  = pooled2c$pooled$fmi,
  lambda = pooled2c$pooled$lambda
)

##                          term         fmi       lambda
## 1                 (Intercept) 0.009125404 0.0079970605
## 2                      age21+ 0.003107807 0.0019856951
## 3                  genderMale 0.002290826 0.0011689820
## 4                 genderOther 0.001652295 0.0005305640
## 5  raceAfrican American/Black 0.011210781 0.0100785488
## 6                   raceAsian 0.014578222 0.0134378570
## 7             raceMultiracial 0.002550856 0.0014289420
## 8                   raceOther 0.014293053 0.0131534645
## 9    ethnicityHispanic/Latino 0.040366773 0.0390911581
## 10     studentstatusPart Time 0.013646981 0.0125090922
## 11         residenceOn Campus 0.001869033 0.0007472727
## 12              greekmemGreek 0.001304156 0.0001824503
## 13                alcmo_dich2 0.016819344 0.0156723092
## 14                marmo_dich2 0.002680207 0.0015582524
## 15               ecigmo_dich2 0.003508554 0.0023862601
## 16               smokmo_dich2 0.002649936 0.0015279911
## 17              cigarmo_dich2 0.004875920 0.0037527521
## 18                (Intercept) 0.005325229 0.0042016899
## 19                     age21+ 0.001462241 0.0003405263
## 20                 genderMale 0.002053482 0.0009316903
## 21                genderOther 0.002431641 0.0013097608
## 22 raceAfrican American/Black 0.010142943 0.0090128131
## 23                  raceAsian 0.005148232 0.0040248444
## 24            raceMultiracial 0.005172716 0.0040493075
## 25                  raceOther 0.014577602 0.0134372385
## 26   ethnicityHispanic/Latino 0.059971580 0.0585120894
## 27     studentstatusPart Time 0.001386354 0.0002646442
## 28         residenceOn Campus 0.002280422 0.0011585808
## 29              greekmemGreek 0.014522194 0.0133819830
## 30                alcmo_dich2 0.002968361 0.0018463050
## 31                marmo_dich2 0.002478179 0.0013562857
## 32               ecigmo_dich2 0.008066063 0.0069393571
## 33               smokmo_dich2 0.002682457 0.0015605016
## 34              cigarmo_dich2 0.006630487 0.0055056318

# Fit model 2d: model c + rx stimulants + rx painkillers + rx sedatives
fit2d <- with(
  imp2,
  nnet::multinom(
    lca2_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
    trace = FALSE
  )
)

# Pool model 2d
pooled2d <- pool(fit2d)
summary(pooled2d)

##                          term     estimate   std.error    statistic       df
## 1                 (Intercept) -2.775468375 0.229346088 -12.10165997 1751.363
## 2                      age21+  0.176325103 0.153728471   1.14699054 1768.201
## 3                  genderMale  1.475549933 0.154273678   9.56449570 1768.662
## 4                 genderOther -0.168553032 0.536922083  -0.31392457 1770.425
## 5  raceAfrican American/Black  0.572967365 0.320992225   1.78498830 1735.070
## 6                   raceAsian -0.305695077 0.339840188  -0.89952598 1713.484
## 7             raceMultiracial -0.579492852 0.425970417  -1.36040633 1768.322
## 8                   raceOther  0.259817549 0.455708221   0.57014014 1718.368
## 9    ethnicityHispanic/Latino -0.100868097 0.236061246  -0.42729630 1499.560
## 10     studentstatusPart Time -0.387142757 0.381550685  -1.01465617 1722.938
## 11         residenceOn Campus  0.230010815 0.155979064   1.47462620 1770.364
## 12              greekmemGreek  0.399695815 0.188886955   2.11605833 1771.748
## 13                alcmo_dich2  0.063535258 0.215820111   0.29438989 1705.022
## 14                marmo_dich2 -0.061888236 0.177772774  -0.34813113 1768.252
## 15               ecigmo_dich2  0.561420726 0.191722779   2.92829432 1765.147
## 16               smokmo_dich2  0.047592250 0.193599195   0.24582876 1769.074
## 17              cigarmo_dich2  0.634928853 0.182521681   3.47864895 1761.106
## 18               rxstmo_dich2 -0.069210317 0.312937393  -0.22116346 1763.836
## 19               rxpkmo_dich2 -0.065359394 0.424443792  -0.15398834 1622.004
## 20              rxsedmo_dich2 -0.034077464 0.519587456  -0.06558562 1762.215
## 21                     mhdays -0.028381889 0.009525474  -2.97957751 1729.866
## 22                (Intercept) -3.202766564 0.285855031 -11.20416372 1748.015
## 23                     age21+  0.296685125 0.225315977   1.31675138 1771.507
## 24                 genderMale  0.777716699 0.219600895   3.54150058 1768.073
## 25                genderOther  1.201198882 0.354549845   3.38795489 1769.676
## 26 raceAfrican American/Black  0.517115976 0.457232753   1.13096880 1741.119
## 27                  raceAsian  1.041543042 0.300020741   3.47157013 1760.738
## 28            raceMultiracial  1.077203056 0.354416060   3.03937428 1763.606
## 29                  raceOther -0.103118834 0.770845507  -0.13377367 1719.675
## 30   ethnicityHispanic/Latino  0.083935603 0.320601586   0.26180657 1269.052
## 31     studentstatusPart Time  0.253618649 0.456317474   0.55579430 1771.284
## 32         residenceOn Campus -0.004130798 0.215677045  -0.01915270 1769.697
## 33              greekmemGreek  0.037070643 0.326470409   0.11354978 1712.458
## 34                alcmo_dich2 -0.310727076 0.264839737  -1.17326455 1767.926
## 35                marmo_dich2 -0.150749185 0.273424038  -0.55133845 1768.967
## 36               ecigmo_dich2  0.211640466 0.289701009   0.73054791 1746.646
## 37               smokmo_dich2 -0.287840997 0.326013190  -0.88291212 1769.670
## 38              cigarmo_dich2 -0.182583957 0.344187376  -0.53047837 1756.110
## 39               rxstmo_dich2  0.444819249 0.496531033   0.89585387 1769.542
## 40               rxpkmo_dich2  0.756445492 0.500116560   1.51253838 1766.803
## 41              rxsedmo_dich2 -0.632284946 0.773647429  -0.81727790 1770.545
## 42                     mhdays  0.022033504 0.010872888   2.02646291 1711.976
##         p.value
## 1  1.959585e-32
## 2  2.515407e-01
## 3  3.628084e-21
## 4  7.536153e-01
## 5  7.443786e-02
## 6  3.684990e-01
## 7  1.738748e-01
## 8  5.686572e-01
## 9  6.692249e-01
## 10 3.104122e-01
## 11 1.404909e-01
## 12 3.447899e-02
## 13 7.684958e-01
## 14 7.277831e-01
## 15 3.451949e-03
## 16 8.058433e-01
## 17 5.162063e-04
## 18 8.249908e-01
## 19 8.776381e-01
## 20 9.477152e-01
## 21 2.926684e-03
## 22 3.461742e-28
## 23 1.880923e-01
## 24 4.081697e-04
## 25 7.195937e-04
## 26 2.582240e-01
## 27 5.299165e-04
## 28 2.405553e-03
## 29 8.935972e-01
## 30 7.935131e-01
## 31 5.784217e-01
## 32 9.847214e-01
## 33 9.096080e-01
## 34 2.408478e-01
## 35 5.814713e-01
## 36 4.651533e-01
## 37 3.774037e-01
## 38 5.958474e-01
## 39 3.704527e-01
## 40 1.305758e-01
## 41 4.138796e-01
## 42 4.287213e-02

## View fraction of missing information (FMI)
data.frame(
  term = summary(pooled2d)$term,
  fmi  = pooled2d$pooled$fmi,
  lambda = pooled2d$pooled$lambda
)

##                          term         fmi       lambda
## 1                 (Intercept) 0.008203442 0.0070714880
## 2                      age21+ 0.002960819 0.0018337121
## 3                  genderMale 0.002763642 0.0016366062
## 4                 genderOther 0.001954053 0.0008272245
## 5  raceAfrican American/Black 0.011769786 0.0106313177
## 6                   raceAsian 0.015602909 0.0144545785
## 7             raceMultiracial 0.002909400 0.0017823127
## 8                   raceOther 0.014795103 0.0136490954
## 9    ethnicityHispanic/Latino 0.040188946 0.0389096758
## 10     studentstatusPart Time 0.014011583 0.0128677034
## 11         residenceOn Campus 0.001983738 0.0008569041
## 12              greekmemGreek 0.001268798 0.0001420372
## 13                alcmo_dich2 0.016939654 0.0157871946
## 14                marmo_dich2 0.002939188 0.0018120892
## 15               ecigmo_dich2 0.004150565 0.0030228570
## 16               smokmo_dich2 0.002582862 0.0014558847
## 17              cigarmo_dich2 0.005508401 0.0043796473
## 18               rxstmo_dich2 0.004613393 0.0034853714
## 19               rxpkmo_dich2 0.027658538 0.0264603388
## 20              rxsedmo_dich2 0.005154856 0.0040264117
## 21                     mhdays 0.012763577 0.0116228344
## 22                (Intercept) 0.009009079 0.0078758801
## 23                     age21+ 0.001399767 0.0002730010
## 24                 genderMale 0.003014375 0.0018872473
## 25                genderOther 0.002309874 0.0011829709
## 26 raceAfrican American/Black 0.010538991 0.0094030627
## 27                  raceAsian 0.005623038 0.0044941786
## 28            raceMultiracial 0.004692142 0.0035640621
## 29                  raceOther 0.014573893 0.0134284980
## 30   ethnicityHispanic/Latino 0.061213522 0.0597351781
## 31     studentstatusPart Time 0.001518626 0.0003918521
## 32         residenceOn Campus 0.002300522 0.0011736213
## 33              greekmemGreek 0.015768990 0.0146201659
## 34                alcmo_dich2 0.003075724 0.0019485720
## 35                marmo_dich2 0.002630603 0.0015036106
## 36               ecigmo_dich2 0.009325543 0.0081918185
## 37               smokmo_dich2 0.002312949 0.0011860456
## 38              cigarmo_dich2 0.006967987 0.0058376855
## 39               rxstmo_dich2 0.002371903 0.0012449842
## 40               rxpkmo_dich2 0.003527471 0.0024001142
## 41              rxsedmo_dich2 0.001894760 0.0007679414
## 42                     mhdays 0.015846495 0.0146974386

Test for outliers.

# Check on a single imputed dataset
single_imp1 <- complete(imp1, 1)
## Null model
model_check1a <- multinom(lca1_class ~ 1, data = single_imp1, trace = FALSE)
## Covariate model
model_check1b <- multinom(
  lca1_class ~
    age + gender + race + ethnicity + studentstatus + residence + greekmem,
  data = single_imp1,
  trace = FALSE
)
## Covariate model + alcohol + marijuana + vaping + smoking + cigars
model_check1c <- multinom(
  lca1_class ~
    age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
  data = single_imp1,
  trace = FALSE
)
## model c + rx stimulants + rx painkillers + rx sedatives
model_check1d <- multinom(
  lca1_class ~
    age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
  data = single_imp1,
  trace = FALSE
)

# Residuals and influence
plot(fitted(model_check1a))

plot(fitted(model_check1b))

plot(fitted(model_check1c))

plot(fitted(model_check1d))

# Check on a single imputed dataset
single_imp2 <- complete(imp2, 1)
## Null model
model_check2a <- multinom(lca2_class ~ 1, data = single_imp2, trace = FALSE)
## Covariate model
model_check2b <- multinom(
  lca2_class ~
    age + gender + race + ethnicity + studentstatus + residence + greekmem,
  data = single_imp2,
  trace = FALSE
)
## Covariate model + alcohol + marijuana + vaping
model_check2c <- multinom(
  lca2_class ~
    age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
  data = single_imp2,
  trace = FALSE
)
## model c + rx stimulants + rx painkillers + rx sedatives
model_check2d <- multinom(
  lca2_class ~
    age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
  data = single_imp2,
  trace = FALSE
)

# Residuals and influence
plot(fitted(model_check2a))

plot(fitted(model_check2b))

plot(fitted(model_check2c))

plot(fitted(model_check2d))

Test for multicollinearity among predictors.

### Test for multicollinearity among predictors
# Fit a simple linear model with same predictors to extract VIF

vif_check1b <- lm(
  as.numeric(lca1_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
  data = df6_model1
)
vif(vif_check1b)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.132489  1        1.064185
## gender        1.016385  2        1.004071
## race          1.174237  4        1.020280
## ethnicity     1.150410  1        1.072572
## studentstatus 1.029457  1        1.014622
## residence     1.205837  1        1.098106
## greekmem      1.069431  1        1.034133

vif_check1c <- lm(
  as.numeric(lca1_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
  data = df6_model1
)

vif(vif_check1c)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.281054  1        1.131837
## gender        1.117460  2        1.028154
## race          1.237150  4        1.026958
## ethnicity     1.150534  1        1.072629
## studentstatus 1.031830  1        1.015791
## residence     1.227470  1        1.107913
## greekmem      1.107798  1        1.052520
## alcmo_dich    1.434667  1        1.197776
## marmo_dich    1.657801  1        1.287556
## ecigmo_dich   1.897877  1        1.377635
## smokmo_dich   1.760500  1        1.326838
## cigarmo_dich  1.451277  1        1.204690

vif_check1d <- lm(
  as.numeric(lca1_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
  data = df6_model1
)
vif(vif_check1d)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.280975  1        1.131801
## gender        1.168048  2        1.039597
## race          1.242775  4        1.027541
## ethnicity     1.153518  1        1.074020
## studentstatus 1.035161  1        1.017429
## residence     1.228151  1        1.108220
## greekmem      1.116074  1        1.056444
## alcmo_dich    1.434321  1        1.197631
## marmo_dich    1.687643  1        1.299093
## ecigmo_dich   1.899183  1        1.378109
## smokmo_dich   1.788407  1        1.337314
## cigarmo_dich  1.476468  1        1.215100
## rxstmo_dich   1.366379  1        1.168922
## rxpkmo_dich   1.366153  1        1.168826
## rxsedmo_dich  1.374381  1        1.172340
## mhdays        1.078938  1        1.038720

# Fit a simple linear model with same predictors to extract VIF
vif_check2b <- lm(
  as.numeric(lca2_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
  data = df6_model2
)
vif(vif_check2b)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.113499  1        1.055225
## gender        1.019402  2        1.004816
## race          1.127217  4        1.015082
## ethnicity     1.112375  1        1.054692
## studentstatus 1.029068  1        1.014430
## residence     1.206640  1        1.098471
## greekmem      1.102474  1        1.049988

vif_check2c <- lm(
  as.numeric(lca2_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
  data = df6_model2
)
vif(vif_check2c)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.292036  1        1.136678
## gender        1.148279  2        1.035170
## race          1.190469  4        1.022033
## ethnicity     1.117230  1        1.056991
## studentstatus 1.030227  1        1.015001
## residence     1.226182  1        1.107331
## greekmem      1.148999  1        1.071914
## alcmo_dich    1.461152  1        1.208781
## marmo_dich    1.651463  1        1.285093
## ecigmo_dich   1.954846  1        1.398158
## smokmo_dich   1.831376  1        1.353284
## cigarmo_dich  1.558880  1        1.248551

vif_check2d <- lm(
  as.numeric(lca2_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
  data = df6_model2
)
vif(vif_check2d)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.288033  1        1.134916
## gender        1.209071  2        1.048608
## race          1.201590  4        1.023221
## ethnicity     1.124182  1        1.060275
## studentstatus 1.035204  1        1.017450
## residence     1.228619  1        1.108431
## greekmem      1.162621  1        1.078249
## alcmo_dich    1.456158  1        1.206714
## marmo_dich    1.682050  1        1.296939
## ecigmo_dich   1.954040  1        1.397870
## smokmo_dich   1.867710  1        1.366642
## cigarmo_dich  1.590546  1        1.261169
## rxstmo_dich   1.513676  1        1.230315
## rxpkmo_dich   1.528437  1        1.236300
## rxsedmo_dich  1.536385  1        1.239510
## mhdays        1.094832  1        1.046342

### Absence of complete separation
## Check for separation in model 1
table(df6_model1$age, df6_model1$lca1_class)

##            
##                3    1    2
##   <21 (ref) 1570  625  212
##   21+       1117  521  233

table(df6_model1$gender, df6_model1$lca1_class)

##               
##                   3    1    2
##   Female (ref) 1780  744  159
##   Male          722  344  280
##   Other         164   55    6

table(df6_model1$race, df6_model1$lca1_class)

##                         
##                             3    1    2
##   White (ref)            1964  962  383
##   African American/Black  162   39   21
##   Asian                   303   57   19
##   Multiracial             125   54   11
##   Other                   104   26   10

table(df6_model1$ethnicity, df6_model1$lca1_class)

##                            
##                                3    1    2
##   Not Hispanic/Latino (ref) 2182  968  375
##   Hispanic/Latino            344  136   38

table(df6_model1$studentstatus, df6_model1$lca1_class)

##                  
##                      3    1    2
##   Full Time (ref) 2582 1092  425
##   Part Time         98   53   17

table(df6_model1$residence, df6_model1$lca1_class)

##                   
##                       3    1    2
##   Off Campus (ref) 1316  524  205
##   On Campus        1370  621  240

table(df6_model1$greekmem, df6_model1$lca1_class)

##                 
##                     3    1    2
##   No Greek (ref) 2396  992  353
##   Greek           287  153   92

table(df6_model1$alcmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1  894  283   75
##   2 1770  850  362

table(df6_model1$marmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 1816  677  208
##   2  865  467  237

table(df6_model1$ecigmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 1956  709  218
##   2  721  434  227

table(df6_model1$smokmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 2230  861  267
##   2  456  285  177

table(df6_model1$cigarmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 2453  961  271
##   2  227  180  172

table(df6_model1$rxstmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 2619 1083  399
##   2   64   58   45

table(df6_model1$rxpkmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 2639 1109  424
##   2   44   36   20

table(df6_model1$rxsedmo_dich, df6_model1$lca1_class)

##    
##        3    1    2
##   1 2642 1123  432
##   2   39   23   13

## Check for separation in model 2
table(df6_model2$age, df6_model2$lca2_class)

##            
##               3   1   2
##   <21 (ref) 762 141  60
##   21+       625 171  57

table(df6_model2$gender, df6_model2$lca2_class)

##               
##                  3   1   2
##   Female (ref) 867  84  50
##   Male         451 224  54
##   Other         64   4  13

table(df6_model2$race, df6_model2$lca2_class)

##                         
##                             3    1    2
##   White (ref)            1155  267   77
##   African American/Black   54   17    6
##   Asian                    81   12   19
##   Multiracial              57    7   12
##   Other                    32    8    2

table(df6_model2$ethnicity, df6_model2$lca2_class)

##                            
##                                3    1    2
##   Not Hispanic/Latino (ref) 1170  261   95
##   Hispanic/Latino            155   30   15

table(df6_model2$studentstatus, df6_model2$lca2_class)

##                  
##                      3    1    2
##   Full Time (ref) 1324  301  111
##   Part Time         60   10    6

table(df6_model2$residence, df6_model2$lca2_class)

##                   
##                      3   1   2
##   Off Campus (ref) 648 140  53
##   On Campus        738 172  64

table(df6_model2$greekmem, df6_model2$lca2_class)

##                 
##                     3    1    2
##   No Greek (ref) 1204  239  103
##   Greek           181   73   13

table(df6_model2$alcmo_dich, df6_model2$lca2_class)

##    
##        3    1    2
##   1  342   47   38
##   2 1031  257   79

table(df6_model2$marmo_dich, df6_model2$lca2_class)

##    
##       3   1   2
##   1 812 146  75
##   2 573 166  42

table(df6_model2$ecigmo_dich, df6_model2$lca2_class)

##    
##       3   1   2
##   1 870 141  77
##   2 514 171  39

table(df6_model2$smokmo_dich, df6_model2$lca2_class)

##    
##        3    1    2
##   1 1038  181   94
##   2  349  130   23

table(df6_model2$cigarmo_dich, df6_model2$lca2_class)

##    
##        3    1    2
##   1 1159  172  100
##   2  223  138   16

table(df6_model2$rxstmo_dich, df6_model2$lca2_class)

##    
##        3    1    2
##   1 1311  282  109
##   2   71   29    8

table(df6_model2$rxpkmo_dich, df6_model2$lca2_class)

##    
##        3    1    2
##   1 1343  298  108
##   2   43   13    8

table(df6_model2$rxsedmo_dich, df6_model2$lca2_class)

##    
##        3    1    2
##   1 1358  304  114
##   2   29    8    3

Conduct demographic statistics for the two separate LCAs.

demo_vars <- c(
  "age", "gender", "race", "ethnicity", "studentstatus", "residence", "greekmem",
  "alcmo_dich", "marmo_dich", "ecigmo_dich", "smokmo_dich", "cigarmo_dich",
  "rxstmo_dich", "rxpkmo_dich", "rxsedmo_dich"
)

var_labels <- c(
  age = "Age", gender = "Gender", race = "Race", ethnicity = "Ethnicity",
  studentstatus = "Student Status", residence = "Residence",
  greekmem = "Greek Membership", alcmo_dich = "Alcohol Use",
  marmo_dich = "Marijuana Use", ecigmo_dich = "E-Cigarette Use",
  smokmo_dich = "Cigarette Smoking", cigarmo_dich = "Cigar Use",
  rxstmo_dich = "Rx Stimulant Use", rxpkmo_dich = "Rx Painkiller Use",
  rxsedmo_dich = "Rx Sedative Use", mhdays = "Mental Health Days"
)

subst_vars <- c(
  "alcmo_dich", "marmo_dich", "ecigmo_dich", "smokmo_dich", "cigarmo_dich",
  "rxstmo_dich", "rxpkmo_dich", "rxsedmo_dich"
)

recode_level <- function(var, level) {
  if (var %in% subst_vars) {
    case_when(level == "1" ~ "None", level == "2" ~ "Any", .default = level)
  } else {
    level
  }
}

make_freq_table <- function(data, class_col) {
  class_sym <- sym(class_col)

  totals <- data |>
    count(class = paste0("Class ", as.character(!!class_sym))) |>
    mutate(variable = "Total N", level = "", `n (%)` = as.character(n)) |>
    select(variable, level, class, `n (%)`) |>
    pivot_wider(names_from = class, values_from = `n (%)`, values_fill = "0")

  freq_rows <- map_dfr(demo_vars, function(var) {
    data |>
      count(class = !!class_sym, level = as.character(.data[[var]])) |>
      group_by(class) |>
      mutate(prop = n / sum(n)) |>
      ungroup() |>
      mutate(
        variable = var_labels[var],
        level = recode_level(var, level),
        class = paste0("Class ", class),
        `n (%)` = sprintf("%d (%.1f%%)", n, prop * 100)
      ) |>
      select(variable, level, class, `n (%)`) |>
      pivot_wider(names_from = class, values_from = `n (%)`, values_fill = "0 (0.0%)")
  })

  mhdays_row <- data |>
    group_by(class = paste0("Class ", as.character(!!class_sym))) |>
    summarise(
      `n (%)` = sprintf("%.2f (%.2f)", mean(mhdays, na.rm = TRUE), sd(mhdays, na.rm = TRUE)),
      .groups = "drop"
    ) |>
    pivot_wider(names_from = class, values_from = `n (%)`) |>
    mutate(variable = "Mental Health Days", level = "Mean (SD)")

  bind_rows(totals, freq_rows, mhdays_row) |>
    select(variable, level, sort(tidyselect::peek_vars()))
}

cat("=== LCA1: Descriptive Statistics by Class ===\n")

## === LCA1: Descriptive Statistics by Class ===

lca1_desc <- make_freq_table(single_imp1, "lca1_class")
print(lca1_desc, n = Inf)

## # A tibble: 36 × 5
##    variable           level                       `Class 1`  `Class 2` `Class 3`
##    <chr>              <chr>                       <chr>      <chr>     <chr>    
##  1 Total N            ""                          1146       445       2687     
##  2 Age                "21+"                       521 (45.5… 233 (52.… 1117 (41…
##  3 Age                "<21 (ref)"                 625 (54.5… 212 (47.… 1570 (58…
##  4 Gender             "Female (ref)"              746 (65.1… 159 (35.… 1796 (66…
##  5 Gender             "Male"                      345 (30.1… 280 (62.… 727 (27.…
##  6 Gender             "Other"                     55 (4.8%)  6 (1.3%)  164 (6.1…
##  7 Race               "African American/Black"    39 (3.4%)  21 (4.7%) 163 (6.1…
##  8 Race               "Asian"                     57 (5.0%)  19 (4.3%) 306 (11.…
##  9 Race               "Multiracial"               54 (4.7%)  11 (2.5%) 127 (4.7…
## 10 Race               "Other"                     26 (2.3%)  11 (2.5%) 106 (3.9…
## 11 Race               "White (ref)"               970 (84.6… 383 (86.… 1985 (73…
## 12 Ethnicity          "Hispanic/Latino"           141 (12.3… 41 (9.2%) 359 (13.…
## 13 Ethnicity          "Not Hispanic/Latino (ref)" 1005 (87.… 404 (90.… 2328 (86…
## 14 Student Status     "Full Time (ref)"           1093 (95.… 428 (96.… 2589 (96…
## 15 Student Status     "Part Time"                 53 (4.6%)  17 (3.8%) 98 (3.6%)
## 16 Residence          "Off Campus (ref)"          525 (45.8… 205 (46.… 1317 (49…
## 17 Residence          "On Campus"                 621 (54.2… 240 (53.… 1370 (51…
## 18 Greek Membership   "Greek"                     153 (13.4… 92 (20.7… 287 (10.…
## 19 Greek Membership   "No Greek (ref)"            993 (86.6… 353 (79.… 2400 (89…
## 20 Alcohol Use        "None"                      283 (24.7… 75 (16.9… 901 (33.…
## 21 Alcohol Use        "Any"                       863 (75.3… 370 (83.… 1786 (66…
## 22 Marijuana Use      "None"                      678 (59.2… 208 (46.… 1820 (67…
## 23 Marijuana Use      "Any"                       468 (40.8… 237 (53.… 867 (32.…
## 24 E-Cigarette Use    "None"                      711 (62.0… 218 (49.… 1963 (73…
## 25 E-Cigarette Use    "Any"                       435 (38.0… 227 (51.… 724 (26.…
## 26 Cigarette Smoking  "None"                      861 (75.1… 268 (60.… 2231 (83…
## 27 Cigarette Smoking  "Any"                       285 (24.9… 177 (39.… 456 (17.…
## 28 Cigar Use          "None"                      966 (84.3… 273 (61.… 2460 (91…
## 29 Cigar Use          "Any"                       180 (15.7… 172 (38.… 227 (8.4…
## 30 Rx Stimulant Use   "None"                      1087 (94.… 400 (89.… 2623 (97…
## 31 Rx Stimulant Use   "Any"                       59 (5.1%)  45 (10.1… 64 (2.4%)
## 32 Rx Painkiller Use  "None"                      1110 (96.… 424 (95.… 2643 (98…
## 33 Rx Painkiller Use  "Any"                       36 (3.1%)  21 (4.7%) 44 (1.6%)
## 34 Rx Sedative Use    "None"                      1123 (98.… 432 (97.… 2648 (98…
## 35 Rx Sedative Use    "Any"                       23 (2.0%)  13 (2.9%) 39 (1.5%)
## 36 Mental Health Days "Mean (SD)"                 9.15 (8.5… 6.38 (7.… 8.59 (8.…

cat("\n=== LCA2: Descriptive Statistics by Class ===\n")

## 
## === LCA2: Descriptive Statistics by Class ===

lca2_desc <- make_freq_table(single_imp2, "lca2_class")
print(lca2_desc, n = Inf)

## # A tibble: 36 × 5
##    variable           level                       `Class 1`  `Class 2` `Class 3`
##    <chr>              <chr>                       <chr>      <chr>     <chr>    
##  1 Total N            ""                          312        117       1387     
##  2 Age                "21+"                       171 (54.8… 57 (48.7… 625 (45.…
##  3 Age                "<21 (ref)"                 141 (45.2… 60 (51.3… 762 (54.…
##  4 Gender             "Female (ref)"              84 (26.9%) 50 (42.7… 871 (62.…
##  5 Gender             "Male"                      224 (71.8… 54 (46.2… 452 (32.…
##  6 Gender             "Other"                     4 (1.3%)   13 (11.1… 64 (4.6%)
##  7 Race               "African American/Black"    17 (5.4%)  6 (5.1%)  55 (4.0%)
##  8 Race               "Asian"                     12 (3.8%)  19 (16.2… 81 (5.8%)
##  9 Race               "Multiracial"               7 (2.2%)   12 (10.3… 57 (4.1%)
## 10 Race               "Other"                     8 (2.6%)   2 (1.7%)  34 (2.5%)
## 11 Race               "White (ref)"               268 (85.9… 78 (66.7… 1160 (83…
## 12 Ethnicity          "Hispanic/Latino"           33 (10.6%) 15 (12.8… 167 (12.…
## 13 Ethnicity          "Not Hispanic/Latino (ref)" 279 (89.4… 102 (87.… 1220 (88…
## 14 Student Status     "Full Time (ref)"           302 (96.8… 111 (94.… 1327 (95…
## 15 Student Status     "Part Time"                 10 (3.2%)  6 (5.1%)  60 (4.3%)
## 16 Residence          "Off Campus (ref)"          140 (44.9… 53 (45.3… 648 (46.…
## 17 Residence          "On Campus"                 172 (55.1… 64 (54.7… 739 (53.…
## 18 Greek Membership   "Greek"                     73 (23.4%) 13 (11.1… 181 (13.…
## 19 Greek Membership   "No Greek (ref)"            239 (76.6… 104 (88.… 1206 (87…
## 20 Alcohol Use        "None"                      47 (15.1%) 38 (32.5… 342 (24.…
## 21 Alcohol Use        "Any"                       265 (84.9… 79 (67.5… 1045 (75…
## 22 Marijuana Use      "None"                      146 (46.8… 75 (64.1… 814 (58.…
## 23 Marijuana Use      "Any"                       166 (53.2… 42 (35.9… 573 (41.…
## 24 E-Cigarette Use    "None"                      141 (45.2… 78 (66.7… 870 (62.…
## 25 E-Cigarette Use    "Any"                       171 (54.8… 39 (33.3… 517 (37.…
## 26 Cigarette Smoking  "None"                      182 (58.3… 94 (80.3… 1038 (74…
## 27 Cigarette Smoking  "Any"                       130 (41.7… 23 (19.7… 349 (25.…
## 28 Cigar Use          "None"                      174 (55.8… 100 (85.… 1164 (83…
## 29 Cigar Use          "Any"                       138 (44.2… 17 (14.5… 223 (16.…
## 30 Rx Stimulant Use   "None"                      283 (90.7… 109 (93.… 1315 (94…
## 31 Rx Stimulant Use   "Any"                       29 (9.3%)  8 (6.8%)  72 (5.2%)
## 32 Rx Painkiller Use  "None"                      298 (95.5… 109 (93.… 1344 (96…
## 33 Rx Painkiller Use  "Any"                       14 (4.5%)  8 (6.8%)  43 (3.1%)
## 34 Rx Sedative Use    "None"                      304 (97.4… 114 (97.… 1358 (97…
## 35 Rx Sedative Use    "Any"                       8 (2.6%)   3 (2.6%)  29 (2.1%)
## 36 Mental Health Days "Mean (SD)"                 6.06 (7.9… 10.47 (9… 8.91 (8.…

Compare log-likelihoods per models.

## --- Model 1: Gambling Participation LCA  ---
pooled_lrt <- function(larger_fit, smaller_fit, label) {
  # Extract per-imputation log-likelihoods
  ll_large <- sapply(larger_fit$analyses, function(m)
    as.numeric(logLik(m)))
  ll_small <- sapply(smaller_fit$analyses, function(m)
    as.numeric(logLik(m)))
  
  # LRT statistic per imputation
  lrt_per_imp <- -2 * (ll_small - ll_large)
  
  # Degrees of freedom: difference in number of estimated parameters
  # attr(logLik(model), "df") returns total params; we want the difference
  df_large <- attr(logLik(larger_fit$analyses[[1]]), "df")
  df_small <- attr(logLik(smaller_fit$analyses[[1]]), "df")
  df_diff  <- df_large - df_small
  
  # Pool: simple mean of LRT statistics across imputations
  mean_lrt <- mean(lrt_per_imp)
  sd_lrt   <- sd(lrt_per_imp)
  
  # P-value from chi-square distribution
  p_val <- pchisq(mean_lrt, df = df_diff, lower.tail = FALSE)
  
  # Report
  cat(
    sprintf(
      "\n--- %s ---\n  Mean LRT (chi-sq): %.3f | df: %d | p: %.4f\n  SD across imputations: %.3f\n",
      label,
      mean_lrt,
      df_diff,
      p_val,
      sd_lrt
    )
  )
  
  data.frame(
    comparison    = label,
    mean_LRT_chisq = round(mean_lrt, 3),
    df            = df_diff,
    p_value       = round(p_val, 4),
    sd_across_imp = round(sd_lrt, 3)
  )
}

lrt_table <- bind_rows(
  pooled_lrt(fit1b, fit1a, "Model 1: A → B"),
  pooled_lrt(fit1c, fit1b, "Model 1: B → C"),
  pooled_lrt(fit1d, fit1c, "Model 1: C → D"),
  pooled_lrt(fit1d, fit1b, "Model 1: B → D"),
  pooled_lrt(fit1c, fit1a, "Model 1: A → C"),
  pooled_lrt(fit1d, fit1a, "Model 1: A → D"),
  pooled_lrt(fit2b, fit2a, "Model 2: A → B"),
  pooled_lrt(fit2c, fit2b, "Model 2: B → C"),
  pooled_lrt(fit2c, fit2a, "Model 2: A → C"),
  pooled_lrt(fit2d, fit2c, "Model 2: C → D"),
  pooled_lrt(fit2d, fit2b, "Model 2: B → D"),
  pooled_lrt(fit2d, fit2a, "Model 2: A → D")
)

## 
## --- Model 1: A → B ---
##   Mean LRT (chi-sq): 379.035 | df: 22 | p: 0.0000
##   SD across imputations: 2.037
## 
## --- Model 1: B → C ---
##   Mean LRT (chi-sq): 155.214 | df: 10 | p: 0.0000
##   SD across imputations: 0.766
## 
## --- Model 1: C → D ---
##   Mean LRT (chi-sq): 28.257 | df: 8 | p: 0.0004
##   SD across imputations: 0.821
## 
## --- Model 1: B → D ---
##   Mean LRT (chi-sq): 183.470 | df: 18 | p: 0.0000
##   SD across imputations: 1.012
## 
## --- Model 1: A → C ---
##   Mean LRT (chi-sq): 534.248 | df: 32 | p: 0.0000
##   SD across imputations: 2.163
## 
## --- Model 1: A → D ---
##   Mean LRT (chi-sq): 562.505 | df: 40 | p: 0.0000
##   SD across imputations: 2.274
## 
## --- Model 2: A → B ---
##   Mean LRT (chi-sq): 239.897 | df: 22 | p: 0.0000
##   SD across imputations: 1.153
## 
## --- Model 2: B → C ---
##   Mean LRT (chi-sq): 51.004 | df: 10 | p: 0.0000
##   SD across imputations: 0.544
## 
## --- Model 2: A → C ---
##   Mean LRT (chi-sq): 290.900 | df: 32 | p: 0.0000
##   SD across imputations: 0.964
## 
## --- Model 2: C → D ---
##   Mean LRT (chi-sq): 19.858 | df: 8 | p: 0.0109
##   SD across imputations: 0.897
## 
## --- Model 2: B → D ---
##   Mean LRT (chi-sq): 70.862 | df: 18 | p: 0.0000
##   SD across imputations: 0.987
## 
## --- Model 2: A → D ---
##   Mean LRT (chi-sq): 310.758 | df: 40 | p: 0.0000
##   SD across imputations: 1.146

print(lrt_table)

##        comparison mean_LRT_chisq df p_value sd_across_imp
## 1  Model 1: A → B        379.035 22  0.0000         2.037
## 2  Model 1: B → C        155.214 10  0.0000         0.766
## 3  Model 1: C → D         28.257  8  0.0004         0.821
## 4  Model 1: B → D        183.470 18  0.0000         1.012
## 5  Model 1: A → C        534.248 32  0.0000         2.163
## 6  Model 1: A → D        562.505 40  0.0000         2.274
## 7  Model 2: A → B        239.897 22  0.0000         1.153
## 8  Model 2: B → C         51.004 10  0.0000         0.544
## 9  Model 2: A → C        290.900 32  0.0000         0.964
## 10 Model 2: C → D         19.858  8  0.0109         0.897
## 11 Model 2: B → D         70.862 18  0.0000         0.987
## 12 Model 2: A → D        310.758 40  0.0000         1.146

Calculate odds ratios.

compute_or_table <- function(pooled_model, model_label, ref_class = 3) {
  s <- summary(pooled_model)
  
  # number of rows per outcome contrast
  n_params_per_contrast <- length(unique(s$term))
  n_contrasts <- nrow(s) / n_params_per_contrast
  
  # use the actual non-reference class labels, not ref_class + 1, +2, ...
  all_classes <- seq_len(n_contrasts + 1)
  nonref_classes <- all_classes[all_classes != ref_class]
  
  if (length(nonref_classes) != n_contrasts) {
    stop("Mismatch between detected contrasts and class labels.")
  }
  
  s$contrast <- rep(paste0("Class ", nonref_classes, " vs. Class ", ref_class),
                    each = n_params_per_contrast)
  
  s |>
    dplyr::mutate(
      OR = round(exp(estimate), 3),
      CI_lower = round(exp(estimate - 1.96 * std.error), 3),
      CI_upper = round(exp(estimate + 1.96 * std.error), 3),
      p_value = round(p.value, 4),
      sig = dplyr::case_when(
        p.value < .001 ~ "***",
        p.value < .01  ~ "**",
        p.value < .05  ~ "*",
        p.value < .10  ~ ".",
        TRUE           ~ ""
      ),
      model = model_label
    ) |>
    dplyr::select(model, contrast, term, OR, CI_lower, CI_upper, p_value, sig)
}

or_1b <- compute_or_table(pooled1b, "Model 1b")
or_1c <- compute_or_table(pooled1c, "Model 1c")
or_1d <- compute_or_table(pooled1d, "Model 1d")
or_2b <- compute_or_table(pooled2b, "Model 2b")
or_2c <- compute_or_table(pooled2c, "Model 2c")
or_2d <- compute_or_table(pooled2d, "Model 2d")

print(or_1b, row.names = FALSE)

##     model            contrast                       term    OR CI_lower
##  Model 1b Class 1 vs. Class 3                (Intercept) 0.404    0.353
##  Model 1b Class 1 vs. Class 3                     age21+ 1.142    0.984
##  Model 1b Class 1 vs. Class 3                 genderMale 1.204    1.029
##  Model 1b Class 1 vs. Class 3                genderOther 0.783    0.569
##  Model 1b Class 1 vs. Class 3 raceAfrican American/Black 0.497    0.347
##  Model 1b Class 1 vs. Class 3                  raceAsian 0.368    0.274
##  Model 1b Class 1 vs. Class 3            raceMultiracial 0.893    0.641
##  Model 1b Class 1 vs. Class 3                  raceOther 0.543    0.344
##  Model 1b Class 1 vs. Class 3   ethnicityHispanic/Latino 0.913    0.725
##  Model 1b Class 1 vs. Class 3     studentstatusPart Time 1.166    0.822
##  Model 1b Class 1 vs. Class 3         residenceOn Campus 1.140    0.977
##  Model 1b Class 1 vs. Class 3              greekmemGreek 1.257    1.011
##  Model 1b Class 2 vs. Class 3                (Intercept) 0.070    0.055
##  Model 1b Class 2 vs. Class 3                     age21+ 1.500    1.200
##  Model 1b Class 2 vs. Class 3                 genderMale 4.638    3.732
##  Model 1b Class 2 vs. Class 3                genderOther 0.422    0.183
##  Model 1b Class 2 vs. Class 3 raceAfrican American/Black 0.683    0.420
##  Model 1b Class 2 vs. Class 3                  raceAsian 0.237    0.146
##  Model 1b Class 2 vs. Class 3            raceMultiracial 0.492    0.258
##  Model 1b Class 2 vs. Class 3                  raceOther 0.595    0.294
##  Model 1b Class 2 vs. Class 3   ethnicityHispanic/Latino 0.699    0.481
##  Model 1b Class 2 vs. Class 3     studentstatusPart Time 0.964    0.554
##  Model 1b Class 2 vs. Class 3         residenceOn Campus 1.250    0.988
##  Model 1b Class 2 vs. Class 3              greekmemGreek 2.005    1.505
##  CI_upper p_value sig
##     0.463  0.0000 ***
##     1.326  0.0804   .
##     1.409  0.0203   *
##     1.078  0.1336    
##     0.712  0.0001 ***
##     0.496  0.0000 ***
##     1.245  0.5054    
##     0.860  0.0092  **
##     1.150  0.4413    
##     1.654  0.3895    
##     1.329  0.0954   .
##     1.564  0.0399   *
##     0.089  0.0000 ***
##     1.874  0.0004 ***
##     5.763  0.0000 ***
##     0.972  0.0427   *
##     1.110  0.1239    
##     0.387  0.0000 ***
##     0.937  0.0311   *
##     1.203  0.1482    
##     1.014  0.0594   .
##     1.679  0.8980    
##     1.581  0.0632   .
##     2.671  0.0000 ***

print(or_1c, row.names = FALSE)

##     model            contrast                       term    OR CI_lower
##  Model 1c Class 1 vs. Class 3                (Intercept) 0.336    0.281
##  Model 1c Class 1 vs. Class 3                     age21+ 1.020    0.870
##  Model 1c Class 1 vs. Class 3                 genderMale 1.178    1.000
##  Model 1c Class 1 vs. Class 3                genderOther 0.776    0.562
##  Model 1c Class 1 vs. Class 3 raceAfrican American/Black 0.533    0.371
##  Model 1c Class 1 vs. Class 3                  raceAsian 0.418    0.309
##  Model 1c Class 1 vs. Class 3            raceMultiracial 0.896    0.641
##  Model 1c Class 1 vs. Class 3                  raceOther 0.554    0.349
##  Model 1c Class 1 vs. Class 3   ethnicityHispanic/Latino 0.920    0.729
##  Model 1c Class 1 vs. Class 3     studentstatusPart Time 1.163    0.818
##  Model 1c Class 1 vs. Class 3         residenceOn Campus 1.129    0.967
##  Model 1c Class 1 vs. Class 3              greekmemGreek 1.145    0.917
##  Model 1c Class 1 vs. Class 3                alcmo_dich2 1.139    0.943
##  Model 1c Class 1 vs. Class 3                marmo_dich2 1.023    0.850
##  Model 1c Class 1 vs. Class 3               ecigmo_dich2 1.308    1.070
##  Model 1c Class 1 vs. Class 3               smokmo_dich2 1.069    0.857
##  Model 1c Class 1 vs. Class 3              cigarmo_dich2 1.424    1.115
##  Model 1c Class 2 vs. Class 3                (Intercept) 0.047    0.034
##  Model 1c Class 2 vs. Class 3                     age21+ 1.118    0.878
##  Model 1c Class 2 vs. Class 3                 genderMale 4.015    3.178
##  Model 1c Class 2 vs. Class 3                genderOther 0.401    0.173
##  Model 1c Class 2 vs. Class 3 raceAfrican American/Black 0.817    0.495
##  Model 1c Class 2 vs. Class 3                  raceAsian 0.340    0.206
##  Model 1c Class 2 vs. Class 3            raceMultiracial 0.520    0.271
##  Model 1c Class 2 vs. Class 3                  raceOther 0.602    0.291
##  Model 1c Class 2 vs. Class 3   ethnicityHispanic/Latino 0.716    0.488
##  Model 1c Class 2 vs. Class 3     studentstatusPart Time 0.986    0.559
##  Model 1c Class 2 vs. Class 3         residenceOn Campus 1.193    0.935
##  Model 1c Class 2 vs. Class 3              greekmemGreek 1.452    1.070
##  Model 1c Class 2 vs. Class 3                alcmo_dich2 1.223    0.887
##  Model 1c Class 2 vs. Class 3                marmo_dich2 1.237    0.936
##  Model 1c Class 2 vs. Class 3               ecigmo_dich2 1.552    1.147
##  Model 1c Class 2 vs. Class 3               smokmo_dich2 1.226    0.900
##  Model 1c Class 2 vs. Class 3              cigarmo_dich2 2.486    1.841
##  CI_upper p_value sig
##     0.400  0.0000 ***
##     1.197  0.8047    
##     1.388  0.0499   *
##     1.071  0.1231    
##     0.765  0.0006 ***
##     0.565  0.0000 ***
##     1.252  0.5201    
##     0.879  0.0123   *
##     1.160  0.4788    
##     1.652  0.4010    
##     1.319  0.1255    
##     1.431  0.2328    
##     1.377  0.1773    
##     1.231  0.8116    
##     1.598  0.0088  **
##     1.332  0.5542    
##     1.819  0.0046  **
##     0.064  0.0000 ***
##     1.423  0.3656    
##     5.073  0.0000 ***
##     0.928  0.0328   *
##     1.349  0.4299    
##     0.562  0.0000 ***
##     0.998  0.0495   *
##     1.246  0.1716    
##     1.051  0.0883   .
##     1.739  0.9615    
##     1.523  0.1555    
##     1.969  0.0165   *
##     1.688  0.2199    
##     1.634  0.1352    
##     2.099  0.0044  **
##     1.671  0.1962    
##     3.358  0.0000 ***

print(or_1d, row.names = FALSE)

##     model            contrast                       term    OR CI_lower
##  Model 1d Class 1 vs. Class 3                (Intercept) 0.317    0.262
##  Model 1d Class 1 vs. Class 3                     age21+ 1.021    0.870
##  Model 1d Class 1 vs. Class 3                 genderMale 1.195    1.013
##  Model 1d Class 1 vs. Class 3                genderOther 0.746    0.538
##  Model 1d Class 1 vs. Class 3 raceAfrican American/Black 0.542    0.378
##  Model 1d Class 1 vs. Class 3                  raceAsian 0.417    0.308
##  Model 1d Class 1 vs. Class 3            raceMultiracial 0.897    0.642
##  Model 1d Class 1 vs. Class 3                  raceOther 0.563    0.354
##  Model 1d Class 1 vs. Class 3   ethnicityHispanic/Latino 0.915    0.725
##  Model 1d Class 1 vs. Class 3     studentstatusPart Time 1.181    0.831
##  Model 1d Class 1 vs. Class 3         residenceOn Campus 1.124    0.962
##  Model 1d Class 1 vs. Class 3              greekmemGreek 1.148    0.918
##  Model 1d Class 1 vs. Class 3                alcmo_dich2 1.139    0.942
##  Model 1d Class 1 vs. Class 3                marmo_dich2 0.994    0.824
##  Model 1d Class 1 vs. Class 3               ecigmo_dich2 1.292    1.057
##  Model 1d Class 1 vs. Class 3               smokmo_dich2 1.050    0.841
##  Model 1d Class 1 vs. Class 3              cigarmo_dich2 1.407    1.099
##  Model 1d Class 1 vs. Class 3               rxstmo_dich2 1.470    0.960
##  Model 1d Class 1 vs. Class 3               rxpkmo_dich2 1.557    0.918
##  Model 1d Class 1 vs. Class 3              rxsedmo_dich2 0.621    0.331
##  Model 1d Class 1 vs. Class 3                     mhdays 1.007    0.998
##  Model 1d Class 2 vs. Class 3                (Intercept) 0.057    0.041
##  Model 1d Class 2 vs. Class 3                     age21+ 1.086    0.852
##  Model 1d Class 2 vs. Class 3                 genderMale 3.761    2.967
##  Model 1d Class 2 vs. Class 3                genderOther 0.443    0.191
##  Model 1d Class 2 vs. Class 3 raceAfrican American/Black 0.811    0.491
##  Model 1d Class 2 vs. Class 3                  raceAsian 0.332    0.200
##  Model 1d Class 2 vs. Class 3            raceMultiracial 0.536    0.279
##  Model 1d Class 2 vs. Class 3                  raceOther 0.640    0.310
##  Model 1d Class 2 vs. Class 3   ethnicityHispanic/Latino 0.701    0.477
##  Model 1d Class 2 vs. Class 3     studentstatusPart Time 1.022    0.579
##  Model 1d Class 2 vs. Class 3         residenceOn Campus 1.214    0.950
##  Model 1d Class 2 vs. Class 3              greekmemGreek 1.374    1.010
##  Model 1d Class 2 vs. Class 3                alcmo_dich2 1.232    0.893
##  Model 1d Class 2 vs. Class 3                marmo_dich2 1.246    0.939
##  Model 1d Class 2 vs. Class 3               ecigmo_dich2 1.553    1.146
##  Model 1d Class 2 vs. Class 3               smokmo_dich2 1.202    0.879
##  Model 1d Class 2 vs. Class 3              cigarmo_dich2 2.389    1.764
##  Model 1d Class 2 vs. Class 3               rxstmo_dich2 1.861    1.107
##  Model 1d Class 2 vs. Class 3               rxpkmo_dich2 1.393    0.683
##  Model 1d Class 2 vs. Class 3              rxsedmo_dich2 0.583    0.251
##  Model 1d Class 2 vs. Class 3                     mhdays 0.976    0.961
##  CI_upper p_value sig
##     0.384  0.0000 ***
##     1.198  0.8012    
##     1.409  0.0351   *
##     1.033  0.0777   .
##     0.778  0.0009 ***
##     0.563  0.0000 ***
##     1.254  0.5253    
##     0.893  0.0148   *
##     1.153  0.4506    
##     1.680  0.3538    
##     1.314  0.1405    
##     1.436  0.2258    
##     1.377  0.1793    
##     1.198  0.9478    
##     1.581  0.0126   *
##     1.311  0.6651    
##     1.801  0.0067  **
##     2.250  0.0762   .
##     2.643  0.1007    
##     1.166  0.1383    
##     1.015  0.1111    
##     0.080  0.0000 ***
##     1.384  0.5065    
##     4.766  0.0000 ***
##     1.030  0.0588   .
##     1.340  0.4139    
##     0.549  0.0000 ***
##     1.030  0.0613   .
##     1.322  0.2283    
##     1.030  0.0708   .
##     1.806  0.9397    
##     1.551  0.1206    
##     1.870  0.0433   *
##     1.701  0.2039    
##     1.654  0.1281    
##     2.104  0.0045  **
##     1.643  0.2502    
##     3.235  0.0000 ***
##     3.129  0.0192   *
##     2.842  0.3621    
##     1.353  0.2089    
##     0.990  0.0013  **

print(or_2b, row.names = FALSE)

##     model            contrast                       term    OR CI_lower
##  Model 2b Class 1 vs. Class 3                (Intercept) 0.063    0.046
##  Model 2b Class 1 vs. Class 3                     age21+ 1.465    1.110
##  Model 2b Class 1 vs. Class 3                 genderMale 5.298    4.003
##  Model 2b Class 1 vs. Class 3                genderOther 0.729    0.258
##  Model 2b Class 1 vs. Class 3 raceAfrican American/Black 1.537    0.837
##  Model 2b Class 1 vs. Class 3                  raceAsian 0.554    0.290
##  Model 2b Class 1 vs. Class 3            raceMultiracial 0.539    0.236
##  Model 2b Class 1 vs. Class 3                  raceOther 1.270    0.530
##  Model 2b Class 1 vs. Class 3   ethnicityHispanic/Latino 0.873    0.556
##  Model 2b Class 1 vs. Class 3     studentstatusPart Time 0.678    0.328
##  Model 2b Class 1 vs. Class 3         residenceOn Campus 1.326    0.987
##  Model 2b Class 1 vs. Class 3              greekmemGreek 1.880    1.329
##  Model 2b Class 2 vs. Class 3                (Intercept) 0.042    0.027
##  Model 2b Class 2 vs. Class 3                     age21+ 1.179    0.786
##  Model 2b Class 2 vs. Class 3                 genderMale 2.000    1.330
##  Model 2b Class 2 vs. Class 3                genderOther 3.611    1.838
##  Model 2b Class 2 vs. Class 3 raceAfrican American/Black 1.818    0.751
##  Model 2b Class 2 vs. Class 3                  raceAsian 3.301    1.880
##  Model 2b Class 2 vs. Class 3            raceMultiracial 2.961    1.490
##  Model 2b Class 2 vs. Class 3                  raceOther 0.955    0.213
##  Model 2b Class 2 vs. Class 3   ethnicityHispanic/Latino 1.154    0.621
##  Model 2b Class 2 vs. Class 3     studentstatusPart Time 1.199    0.493
##  Model 2b Class 2 vs. Class 3         residenceOn Campus 1.029    0.676
##  Model 2b Class 2 vs. Class 3              greekmemGreek 0.944    0.506
##  CI_upper p_value sig
##     0.087  0.0000 ***
##     1.932  0.0070  **
##     7.011  0.0000 ***
##     2.062  0.5519    
##     2.825  0.1661    
##     1.057  0.0732   .
##     1.232  0.1429    
##     3.047  0.5920    
##     1.371  0.5556    
##     1.401  0.2946    
##     1.782  0.0616   .
##     2.660  0.0004 ***
##     0.064  0.0000 ***
##     1.769  0.4269    
##     3.007  0.0009 ***
##     7.095  0.0002 ***
##     4.400  0.1851    
##     5.796  0.0000 ***
##     5.885  0.0020  **
##     4.277  0.9516    
##     2.145  0.6503    
##     2.914  0.6893    
##     1.567  0.8944    
##     1.762  0.8566

print(or_2c, row.names = FALSE)

##     model            contrast                       term    OR CI_lower
##  Model 2c Class 1 vs. Class 3                (Intercept) 0.050    0.033
##  Model 2c Class 1 vs. Class 3                     age21+ 1.207    0.894
##  Model 2c Class 1 vs. Class 3                 genderMale 4.649    3.451
##  Model 2c Class 1 vs. Class 3                genderOther 0.737    0.259
##  Model 2c Class 1 vs. Class 3 raceAfrican American/Black 1.817    0.973
##  Model 2c Class 1 vs. Class 3                  raceAsian 0.740    0.381
##  Model 2c Class 1 vs. Class 3            raceMultiracial 0.549    0.238
##  Model 2c Class 1 vs. Class 3                  raceOther 1.252    0.516
##  Model 2c Class 1 vs. Class 3   ethnicityHispanic/Latino 0.905    0.572
##  Model 2c Class 1 vs. Class 3     studentstatusPart Time 0.683    0.326
##  Model 2c Class 1 vs. Class 3         residenceOn Campus 1.233    0.910
##  Model 2c Class 1 vs. Class 3              greekmemGreek 1.526    1.058
##  Model 2c Class 1 vs. Class 3                alcmo_dich2 1.077    0.707
##  Model 2c Class 1 vs. Class 3                marmo_dich2 0.902    0.639
##  Model 2c Class 1 vs. Class 3               ecigmo_dich2 1.684    1.159
##  Model 2c Class 1 vs. Class 3               smokmo_dich2 1.016    0.698
##  Model 2c Class 1 vs. Class 3              cigarmo_dich2 1.970    1.382
##  Model 2c Class 2 vs. Class 3                (Intercept) 0.050    0.030
##  Model 2c Class 2 vs. Class 3                     age21+ 1.343    0.866
##  Model 2c Class 2 vs. Class 3                 genderMale 2.047    1.341
##  Model 2c Class 2 vs. Class 3                genderOther 3.672    1.856
##  Model 2c Class 2 vs. Class 3 raceAfrican American/Black 1.712    0.702
##  Model 2c Class 2 vs. Class 3                  raceAsian 2.998    1.674
##  Model 2c Class 2 vs. Class 3            raceMultiracial 2.919    1.466
##  Model 2c Class 2 vs. Class 3                  raceOther 0.935    0.208
##  Model 2c Class 2 vs. Class 3   ethnicityHispanic/Latino 1.124    0.604
##  Model 2c Class 2 vs. Class 3     studentstatusPart Time 1.223    0.502
##  Model 2c Class 2 vs. Class 3         residenceOn Campus 1.015    0.666
##  Model 2c Class 2 vs. Class 3              greekmemGreek 1.039    0.551
##  Model 2c Class 2 vs. Class 3                alcmo_dich2 0.739    0.441
##  Model 2c Class 2 vs. Class 3                marmo_dich2 0.899    0.534
##  Model 2c Class 2 vs. Class 3               ecigmo_dich2 1.292    0.737
##  Model 2c Class 2 vs. Class 3               smokmo_dich2 0.806    0.434
##  Model 2c Class 2 vs. Class 3              cigarmo_dich2 0.861    0.447
##  CI_upper p_value sig
##     0.077  0.0000 ***
##     1.630  0.2192    
##     6.264  0.0000 ***
##     2.099  0.5676    
##     3.393  0.0609   .
##     1.437  0.3742    
##     1.265  0.1594    
##     3.034  0.6191    
##     1.433  0.6714    
##     1.434  0.3140    
##     1.671  0.1759    
##     2.201  0.0239   *
##     1.641  0.7293    
##     1.272  0.5557    
##     2.446  0.0063  **
##     1.480  0.9336    
##     2.810  0.0002 ***
##     0.083  0.0000 ***
##     2.083  0.1884    
##     3.126  0.0009 ***
##     7.263  0.0002 ***
##     4.173  0.2372    
##     5.369  0.0002 ***
##     5.813  0.0023  **
##     4.213  0.9306    
##     2.094  0.7122    
##     2.976  0.6579    
##     1.547  0.9458    
##     1.962  0.9051    
##     1.238  0.2509    
##     1.515  0.6894    
##     2.265  0.3711    
##     1.495  0.4934    
##     1.659  0.6556

print(or_2d, row.names = FALSE)

##     model            contrast                       term    OR CI_lower
##  Model 2d Class 1 vs. Class 3                (Intercept) 0.062    0.040
##  Model 2d Class 1 vs. Class 3                     age21+ 1.193    0.883
##  Model 2d Class 1 vs. Class 3                 genderMale 4.373    3.232
##  Model 2d Class 1 vs. Class 3                genderOther 0.845    0.295
##  Model 2d Class 1 vs. Class 3 raceAfrican American/Black 1.774    0.945
##  Model 2d Class 1 vs. Class 3                  raceAsian 0.737    0.378
##  Model 2d Class 1 vs. Class 3            raceMultiracial 0.560    0.243
##  Model 2d Class 1 vs. Class 3                  raceOther 1.297    0.531
##  Model 2d Class 1 vs. Class 3   ethnicityHispanic/Latino 0.904    0.569
##  Model 2d Class 1 vs. Class 3     studentstatusPart Time 0.679    0.321
##  Model 2d Class 1 vs. Class 3         residenceOn Campus 1.259    0.927
##  Model 2d Class 1 vs. Class 3              greekmemGreek 1.491    1.030
##  Model 2d Class 1 vs. Class 3                alcmo_dich2 1.066    0.698
##  Model 2d Class 1 vs. Class 3                marmo_dich2 0.940    0.663
##  Model 2d Class 1 vs. Class 3               ecigmo_dich2 1.753    1.204
##  Model 2d Class 1 vs. Class 3               smokmo_dich2 1.049    0.718
##  Model 2d Class 1 vs. Class 3              cigarmo_dich2 1.887    1.319
##  Model 2d Class 1 vs. Class 3               rxstmo_dich2 0.933    0.505
##  Model 2d Class 1 vs. Class 3               rxpkmo_dich2 0.937    0.408
##  Model 2d Class 1 vs. Class 3              rxsedmo_dich2 0.966    0.349
##  Model 2d Class 1 vs. Class 3                     mhdays 0.972    0.954
##  Model 2d Class 2 vs. Class 3                (Intercept) 0.041    0.023
##  Model 2d Class 2 vs. Class 3                     age21+ 1.345    0.865
##  Model 2d Class 2 vs. Class 3                 genderMale 2.176    1.415
##  Model 2d Class 2 vs. Class 3                genderOther 3.324    1.659
##  Model 2d Class 2 vs. Class 3 raceAfrican American/Black 1.677    0.685
##  Model 2d Class 2 vs. Class 3                  raceAsian 2.834    1.574
##  Model 2d Class 2 vs. Class 3            raceMultiracial 2.936    1.466
##  Model 2d Class 2 vs. Class 3                  raceOther 0.902    0.199
##  Model 2d Class 2 vs. Class 3   ethnicityHispanic/Latino 1.088    0.580
##  Model 2d Class 2 vs. Class 3     studentstatusPart Time 1.289    0.527
##  Model 2d Class 2 vs. Class 3         residenceOn Campus 0.996    0.653
##  Model 2d Class 2 vs. Class 3              greekmemGreek 1.038    0.547
##  Model 2d Class 2 vs. Class 3                alcmo_dich2 0.733    0.436
##  Model 2d Class 2 vs. Class 3                marmo_dich2 0.860    0.503
##  Model 2d Class 2 vs. Class 3               ecigmo_dich2 1.236    0.700
##  Model 2d Class 2 vs. Class 3               smokmo_dich2 0.750    0.396
##  Model 2d Class 2 vs. Class 3              cigarmo_dich2 0.833    0.424
##  Model 2d Class 2 vs. Class 3               rxstmo_dich2 1.560    0.590
##  Model 2d Class 2 vs. Class 3               rxpkmo_dich2 2.131    0.799
##  Model 2d Class 2 vs. Class 3              rxsedmo_dich2 0.531    0.117
##  Model 2d Class 2 vs. Class 3                     mhdays 1.022    1.001
##  CI_upper p_value sig
##     0.098  0.0000 ***
##     1.612  0.2515    
##     5.918  0.0000 ***
##     2.420  0.7536    
##     3.327  0.0744   .
##     1.434  0.3685    
##     1.291  0.1739    
##     3.168  0.5687    
##     1.436  0.6692    
##     1.434  0.3104    
##     1.709  0.1405    
##     2.160  0.0345   *
##     1.627  0.7685    
##     1.332  0.7278    
##     2.553  0.0035  **
##     1.533  0.8058    
##     2.698  0.0005 ***
##     1.723  0.8250    
##     2.152  0.8776    
##     2.676  0.9477    
##     0.990  0.0029  **
##     0.071  0.0000 ***
##     2.092  0.1881    
##     3.347  0.0004 ***
##     6.660  0.0007 ***
##     4.109  0.2582    
##     5.102  0.0005 ***
##     5.882  0.0024  **
##     4.087  0.8936    
##     2.039  0.7935    
##     3.152  0.5784    
##     1.520  0.9847    
##     1.968  0.9096    
##     1.232  0.2408    
##     1.470  0.5815    
##     2.180  0.4652    
##     1.421  0.3774    
##     1.636  0.5958    
##     4.129  0.3705    
##     5.678  0.1306    
##     2.421  0.4139    
##     1.044  0.0429   *

Evaluate approximate Pseudo-R^2 effect sizes.

#   McFadden's R²   = 1 - (LL_full / LL_null)
#                     Strict; values of .2–.4 indicate excellent fit.
#                     Best for comparing models with the same outcome.
#
#   Cox & Snell R²  = 1 - exp((2/n) * (LL_null - LL_full))
#                     Analogous to OLS R² but cannot reach 1.0.
#
#   Nagelkerke R²   = Cox-Snell R² / (1 - exp((2/n) * LL_null))
#                     Rescales Cox-Snell to [0, 1].

compute_pseudo_r2 <- function(null_model, full_model, n, model_label) {
  ll_null <- as.numeric(logLik(null_model))
  ll_full <- as.numeric(logLik(full_model))
  
  mcfadden    <- round(1 - (ll_full / ll_null), 6)
  cox_snell   <- round(1 - exp((2 / n) * (ll_null - ll_full)), 4)
  nagelkerke_max <- 1 - exp((2 / n) * ll_null)
  nagelkerke  <- round(cox_snell / nagelkerke_max, 4)
  
  data.frame(
    model       = model_label,
    LL_null     = round(ll_null, 2),
    LL_full     = round(ll_full, 2),
    McFadden    = mcfadden,
    Cox_Snell   = cox_snell,
    Nagelkerke  = nagelkerke
  )
}

## Model 1 pseudo-R² (single imputation approximation)
n1 <- nrow(single_imp1)

model1a_si <- multinom(lca1_class ~ 1, data = single_imp1, trace = FALSE)
model1b_si <- multinom(
  lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
  data = single_imp1,
  trace = FALSE
)
model1c_si <- multinom(
  lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
  data = single_imp1,
  trace = FALSE
)
model1d_si <- multinom(
  lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
  data = single_imp1,
  trace = FALSE
)

## Model 2 pseudo-R² (single imputation approximation)
n2 <- nrow(single_imp2)
model2a_si <- multinom(lca2_class ~ 1, data = single_imp2, trace = FALSE)
model2b_si <- multinom(
  lca2_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
  data = single_imp2,
  trace = FALSE
)
model2c_si <- multinom(
  lca2_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
  data = single_imp2,
  trace = FALSE
)
model2d_si <- multinom(
  lca2_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
  data = single_imp2,
  trace = FALSE
)

pseudo_r2_table <- bind_rows(
  compute_pseudo_r2(model1a_si, model1b_si, n1, "Model 1b"),
  compute_pseudo_r2(model1a_si, model1c_si, n1, "Model 1c"),
  compute_pseudo_r2(model1a_si, model1d_si, n1, "Model 1d"),
  compute_pseudo_r2(model2a_si, model2b_si, n2, "Model 2b"),
  compute_pseudo_r2(model2a_si, model2c_si, n2, "Model 2c"),
  compute_pseudo_r2(model2a_si, model2d_si, n2, "Model 2d")
)
print(pseudo_r2_table, row.names = FALSE)

##     model  LL_null  LL_full McFadden Cox_Snell Nagelkerke
##  Model 1b -3766.25 -3575.90 0.050540    0.0851     0.1028
##  Model 1c -3766.25 -3498.08 0.071204    0.1178     0.1423
##  Model 1d -3766.25 -3484.11 0.074911    0.1236     0.1493
##  Model 2b -1244.18 -1124.35 0.096315    0.1236     0.1657
##  Model 2c -1244.18 -1099.10 0.116608    0.1477     0.1980
##  Model 2d -1244.18 -1089.69 0.124172    0.1565     0.2098

# Incremental pseudo-R² (how much variance each block adds)
incr_r2 <- data.frame(
  comparison = c("1b - 1a", "1c - 1b", "1d - 1c", "2b - 2a", "2c - 2b", "2d - 2c"),
  delta_McFadden  = round(
    c(
      pseudo_r2_table$McFadden[1] - 0,
      pseudo_r2_table$McFadden[2] - pseudo_r2_table$McFadden[1],
      pseudo_r2_table$McFadden[3] - pseudo_r2_table$McFadden[2],
      pseudo_r2_table$McFadden[4] - 0,
      pseudo_r2_table$McFadden[5] - pseudo_r2_table$McFadden[4],
      pseudo_r2_table$McFadden[6] - pseudo_r2_table$McFadden[5]
    ),
    4
  ),
  delta_Nagelkerke = round(
    c(
      pseudo_r2_table$Nagelkerke[1] - 0,
      pseudo_r2_table$Nagelkerke[2] - pseudo_r2_table$Nagelkerke[1],
      pseudo_r2_table$Nagelkerke[3] - pseudo_r2_table$Nagelkerke[2],
      pseudo_r2_table$Nagelkerke[4] - 0,
      pseudo_r2_table$Nagelkerke[5] - pseudo_r2_table$Nagelkerke[4],
      pseudo_r2_table$Nagelkerke[6] - pseudo_r2_table$Nagelkerke[5]
    ),
    4
  )
)
print(incr_r2, row.names = FALSE)

##  comparison delta_McFadden delta_Nagelkerke
##     1b - 1a         0.0505           0.1028
##     1c - 1b         0.0207           0.0395
##     1d - 1c         0.0037           0.0070
##     2b - 2a         0.0963           0.1657
##     2c - 2b         0.0203           0.0323
##     2d - 2c         0.0076           0.0118

Generate a combined results table.

build_results_table <- function(pooled_model,
                                r2_row,
                                model_label,
                                ref_class = 3) {
  s <- summary(pooled_model)
  
  # number of rows per outcome contrast
  n_params_per_contrast <- length(unique(s$term))
  n_contrasts <- nrow(s) / n_params_per_contrast
  
  # correct contrast labels
  all_classes <- seq_len(n_contrasts + 1)
  nonref_classes <- all_classes[all_classes != ref_class]
  
  if (length(nonref_classes) != n_contrasts) {
    stop("Mismatch between detected contrasts and class labels.")
  }
  
  s$contrast <- rep(paste0("Class ", nonref_classes, " vs. Class ", ref_class),
                    each = n_params_per_contrast)
  
  results <- s |>
    dplyr::mutate(
      model = model_label,
      OR = exp(estimate),
      CI_lower = exp(estimate - 1.96 * std.error),
      CI_upper = exp(estimate + 1.96 * std.error),
      OR_CI = sprintf("%.3f [%.3f, %.3f]", OR, CI_lower, CI_upper),
      p.value = round(p.value, 4),
      sig = dplyr::case_when(
        p.value < .001 ~ "***",
        p.value < .01  ~ "**",
        p.value < .05  ~ "*",
        p.value < .10  ~ ".",
        TRUE           ~ ""
      ),
      fmi = round(pooled_model$pooled$fmi, 3),
      McFadden = round(r2_row$McFadden, 3),
      Nagelkerke = round(r2_row$Nagelkerke, 3)
    ) |>
    dplyr::select(model,
                  contrast,
                  term,
                  OR_CI,
                  p.value,
                  sig,
                  fmi,
                  McFadden,
                  Nagelkerke)
  
  return(results)
}
full_table <- dplyr::bind_rows(
  build_results_table(pooled1b, pseudo_r2_table[1, ], "1b"),
  build_results_table(pooled1c, pseudo_r2_table[2, ], "1c"),
  build_results_table(pooled1d, pseudo_r2_table[3, ], "1d"),
  build_results_table(pooled2b, pseudo_r2_table[4, ], "2b"),
  build_results_table(pooled2c, pseudo_r2_table[5, ], "2c"),
  build_results_table(pooled2d, pseudo_r2_table[6, ], "2d")
)

write_csv(full_table, file = "full_table.csv")
print(full_table, width = getOption("width"))

##     model            contrast                       term                OR_CI
## 1      1b Class 1 vs. Class 3                (Intercept) 0.404 [0.353, 0.463]
## 2      1b Class 1 vs. Class 3                     age21+ 1.142 [0.984, 1.326]
## 3      1b Class 1 vs. Class 3                 genderMale 1.204 [1.029, 1.409]
## 4      1b Class 1 vs. Class 3                genderOther 0.783 [0.569, 1.078]
## 5      1b Class 1 vs. Class 3 raceAfrican American/Black 0.497 [0.347, 0.712]
## 6      1b Class 1 vs. Class 3                  raceAsian 0.368 [0.274, 0.496]
## 7      1b Class 1 vs. Class 3            raceMultiracial 0.893 [0.641, 1.245]
## 8      1b Class 1 vs. Class 3                  raceOther 0.543 [0.344, 0.860]
## 9      1b Class 1 vs. Class 3   ethnicityHispanic/Latino 0.913 [0.725, 1.150]
## 10     1b Class 1 vs. Class 3     studentstatusPart Time 1.166 [0.822, 1.654]
## 11     1b Class 1 vs. Class 3         residenceOn Campus 1.140 [0.977, 1.329]
## 12     1b Class 1 vs. Class 3              greekmemGreek 1.257 [1.011, 1.564]
## 13     1b Class 2 vs. Class 3                (Intercept) 0.070 [0.055, 0.089]
## 14     1b Class 2 vs. Class 3                     age21+ 1.500 [1.200, 1.874]
## 15     1b Class 2 vs. Class 3                 genderMale 4.638 [3.732, 5.763]
## 16     1b Class 2 vs. Class 3                genderOther 0.422 [0.183, 0.972]
## 17     1b Class 2 vs. Class 3 raceAfrican American/Black 0.683 [0.420, 1.110]
## 18     1b Class 2 vs. Class 3                  raceAsian 0.237 [0.146, 0.387]
## 19     1b Class 2 vs. Class 3            raceMultiracial 0.492 [0.258, 0.937]
## 20     1b Class 2 vs. Class 3                  raceOther 0.595 [0.294, 1.203]
## 21     1b Class 2 vs. Class 3   ethnicityHispanic/Latino 0.699 [0.481, 1.014]
## 22     1b Class 2 vs. Class 3     studentstatusPart Time 0.964 [0.554, 1.679]
## 23     1b Class 2 vs. Class 3         residenceOn Campus 1.250 [0.988, 1.581]
## 24     1b Class 2 vs. Class 3              greekmemGreek 2.005 [1.505, 2.671]
## 25     1c Class 1 vs. Class 3                (Intercept) 0.336 [0.281, 0.400]
## 26     1c Class 1 vs. Class 3                     age21+ 1.020 [0.870, 1.197]
## 27     1c Class 1 vs. Class 3                 genderMale 1.178 [1.000, 1.388]
## 28     1c Class 1 vs. Class 3                genderOther 0.776 [0.562, 1.071]
## 29     1c Class 1 vs. Class 3 raceAfrican American/Black 0.533 [0.371, 0.765]
## 30     1c Class 1 vs. Class 3                  raceAsian 0.418 [0.309, 0.565]
## 31     1c Class 1 vs. Class 3            raceMultiracial 0.896 [0.641, 1.252]
## 32     1c Class 1 vs. Class 3                  raceOther 0.554 [0.349, 0.879]
## 33     1c Class 1 vs. Class 3   ethnicityHispanic/Latino 0.920 [0.729, 1.160]
## 34     1c Class 1 vs. Class 3     studentstatusPart Time 1.163 [0.818, 1.652]
## 35     1c Class 1 vs. Class 3         residenceOn Campus 1.129 [0.967, 1.319]
## 36     1c Class 1 vs. Class 3              greekmemGreek 1.145 [0.917, 1.431]
## 37     1c Class 1 vs. Class 3                alcmo_dich2 1.139 [0.943, 1.377]
## 38     1c Class 1 vs. Class 3                marmo_dich2 1.023 [0.850, 1.231]
## 39     1c Class 1 vs. Class 3               ecigmo_dich2 1.308 [1.070, 1.598]
## 40     1c Class 1 vs. Class 3               smokmo_dich2 1.069 [0.857, 1.332]
## 41     1c Class 1 vs. Class 3              cigarmo_dich2 1.424 [1.115, 1.819]
## 42     1c Class 2 vs. Class 3                (Intercept) 0.047 [0.034, 0.064]
## 43     1c Class 2 vs. Class 3                     age21+ 1.118 [0.878, 1.423]
## 44     1c Class 2 vs. Class 3                 genderMale 4.015 [3.178, 5.073]
## 45     1c Class 2 vs. Class 3                genderOther 0.401 [0.173, 0.928]
## 46     1c Class 2 vs. Class 3 raceAfrican American/Black 0.817 [0.495, 1.349]
## 47     1c Class 2 vs. Class 3                  raceAsian 0.340 [0.206, 0.562]
## 48     1c Class 2 vs. Class 3            raceMultiracial 0.520 [0.271, 0.998]
## 49     1c Class 2 vs. Class 3                  raceOther 0.602 [0.291, 1.246]
## 50     1c Class 2 vs. Class 3   ethnicityHispanic/Latino 0.716 [0.488, 1.051]
## 51     1c Class 2 vs. Class 3     studentstatusPart Time 0.986 [0.559, 1.739]
## 52     1c Class 2 vs. Class 3         residenceOn Campus 1.193 [0.935, 1.523]
## 53     1c Class 2 vs. Class 3              greekmemGreek 1.452 [1.070, 1.969]
## 54     1c Class 2 vs. Class 3                alcmo_dich2 1.223 [0.887, 1.688]
## 55     1c Class 2 vs. Class 3                marmo_dich2 1.237 [0.936, 1.634]
## 56     1c Class 2 vs. Class 3               ecigmo_dich2 1.552 [1.147, 2.099]
## 57     1c Class 2 vs. Class 3               smokmo_dich2 1.226 [0.900, 1.671]
## 58     1c Class 2 vs. Class 3              cigarmo_dich2 2.486 [1.841, 3.358]
## 59     1d Class 1 vs. Class 3                (Intercept) 0.317 [0.262, 0.384]
## 60     1d Class 1 vs. Class 3                     age21+ 1.021 [0.870, 1.198]
## 61     1d Class 1 vs. Class 3                 genderMale 1.195 [1.013, 1.409]
## 62     1d Class 1 vs. Class 3                genderOther 0.746 [0.538, 1.033]
## 63     1d Class 1 vs. Class 3 raceAfrican American/Black 0.542 [0.378, 0.778]
## 64     1d Class 1 vs. Class 3                  raceAsian 0.417 [0.308, 0.563]
## 65     1d Class 1 vs. Class 3            raceMultiracial 0.897 [0.642, 1.254]
## 66     1d Class 1 vs. Class 3                  raceOther 0.563 [0.354, 0.893]
## 67     1d Class 1 vs. Class 3   ethnicityHispanic/Latino 0.915 [0.725, 1.153]
## 68     1d Class 1 vs. Class 3     studentstatusPart Time 1.181 [0.831, 1.680]
## 69     1d Class 1 vs. Class 3         residenceOn Campus 1.124 [0.962, 1.314]
## 70     1d Class 1 vs. Class 3              greekmemGreek 1.148 [0.918, 1.436]
## 71     1d Class 1 vs. Class 3                alcmo_dich2 1.139 [0.942, 1.377]
## 72     1d Class 1 vs. Class 3                marmo_dich2 0.994 [0.824, 1.198]
## 73     1d Class 1 vs. Class 3               ecigmo_dich2 1.292 [1.057, 1.581]
## 74     1d Class 1 vs. Class 3               smokmo_dich2 1.050 [0.841, 1.311]
## 75     1d Class 1 vs. Class 3              cigarmo_dich2 1.407 [1.099, 1.801]
## 76     1d Class 1 vs. Class 3               rxstmo_dich2 1.470 [0.960, 2.250]
## 77     1d Class 1 vs. Class 3               rxpkmo_dich2 1.557 [0.918, 2.643]
## 78     1d Class 1 vs. Class 3              rxsedmo_dich2 0.621 [0.331, 1.166]
## 79     1d Class 1 vs. Class 3                     mhdays 1.007 [0.998, 1.015]
## 80     1d Class 2 vs. Class 3                (Intercept) 0.057 [0.041, 0.080]
## 81     1d Class 2 vs. Class 3                     age21+ 1.086 [0.852, 1.384]
## 82     1d Class 2 vs. Class 3                 genderMale 3.761 [2.967, 4.766]
## 83     1d Class 2 vs. Class 3                genderOther 0.443 [0.191, 1.030]
## 84     1d Class 2 vs. Class 3 raceAfrican American/Black 0.811 [0.491, 1.340]
## 85     1d Class 2 vs. Class 3                  raceAsian 0.332 [0.200, 0.549]
## 86     1d Class 2 vs. Class 3            raceMultiracial 0.536 [0.279, 1.030]
## 87     1d Class 2 vs. Class 3                  raceOther 0.640 [0.310, 1.322]
## 88     1d Class 2 vs. Class 3   ethnicityHispanic/Latino 0.701 [0.477, 1.030]
## 89     1d Class 2 vs. Class 3     studentstatusPart Time 1.022 [0.579, 1.806]
## 90     1d Class 2 vs. Class 3         residenceOn Campus 1.214 [0.950, 1.551]
## 91     1d Class 2 vs. Class 3              greekmemGreek 1.374 [1.010, 1.870]
## 92     1d Class 2 vs. Class 3                alcmo_dich2 1.232 [0.893, 1.701]
## 93     1d Class 2 vs. Class 3                marmo_dich2 1.246 [0.939, 1.654]
## 94     1d Class 2 vs. Class 3               ecigmo_dich2 1.553 [1.146, 2.104]
## 95     1d Class 2 vs. Class 3               smokmo_dich2 1.202 [0.879, 1.643]
## 96     1d Class 2 vs. Class 3              cigarmo_dich2 2.389 [1.764, 3.235]
## 97     1d Class 2 vs. Class 3               rxstmo_dich2 1.861 [1.107, 3.129]
## 98     1d Class 2 vs. Class 3               rxpkmo_dich2 1.393 [0.683, 2.842]
## 99     1d Class 2 vs. Class 3              rxsedmo_dich2 0.583 [0.251, 1.353]
## 100    1d Class 2 vs. Class 3                     mhdays 0.976 [0.961, 0.990]
## 101    2b Class 1 vs. Class 3                (Intercept) 0.063 [0.046, 0.087]
## 102    2b Class 1 vs. Class 3                     age21+ 1.465 [1.110, 1.932]
## 103    2b Class 1 vs. Class 3                 genderMale 5.298 [4.003, 7.011]
## 104    2b Class 1 vs. Class 3                genderOther 0.729 [0.258, 2.062]
## 105    2b Class 1 vs. Class 3 raceAfrican American/Black 1.537 [0.837, 2.825]
## 106    2b Class 1 vs. Class 3                  raceAsian 0.554 [0.290, 1.057]
## 107    2b Class 1 vs. Class 3            raceMultiracial 0.539 [0.236, 1.232]
## 108    2b Class 1 vs. Class 3                  raceOther 1.270 [0.530, 3.047]
## 109    2b Class 1 vs. Class 3   ethnicityHispanic/Latino 0.873 [0.556, 1.371]
## 110    2b Class 1 vs. Class 3     studentstatusPart Time 0.678 [0.328, 1.401]
## 111    2b Class 1 vs. Class 3         residenceOn Campus 1.326 [0.987, 1.782]
## 112    2b Class 1 vs. Class 3              greekmemGreek 1.880 [1.329, 2.660]
## 113    2b Class 2 vs. Class 3                (Intercept) 0.042 [0.027, 0.064]
## 114    2b Class 2 vs. Class 3                     age21+ 1.179 [0.786, 1.769]
## 115    2b Class 2 vs. Class 3                 genderMale 2.000 [1.330, 3.007]
## 116    2b Class 2 vs. Class 3                genderOther 3.611 [1.838, 7.095]
## 117    2b Class 2 vs. Class 3 raceAfrican American/Black 1.818 [0.751, 4.400]
## 118    2b Class 2 vs. Class 3                  raceAsian 3.301 [1.880, 5.796]
## 119    2b Class 2 vs. Class 3            raceMultiracial 2.961 [1.490, 5.885]
## 120    2b Class 2 vs. Class 3                  raceOther 0.955 [0.213, 4.277]
## 121    2b Class 2 vs. Class 3   ethnicityHispanic/Latino 1.154 [0.621, 2.145]
## 122    2b Class 2 vs. Class 3     studentstatusPart Time 1.199 [0.493, 2.914]
## 123    2b Class 2 vs. Class 3         residenceOn Campus 1.029 [0.676, 1.567]
## 124    2b Class 2 vs. Class 3              greekmemGreek 0.944 [0.506, 1.762]
## 125    2c Class 1 vs. Class 3                (Intercept) 0.050 [0.033, 0.077]
## 126    2c Class 1 vs. Class 3                     age21+ 1.207 [0.894, 1.630]
## 127    2c Class 1 vs. Class 3                 genderMale 4.649 [3.451, 6.264]
## 128    2c Class 1 vs. Class 3                genderOther 0.737 [0.259, 2.099]
## 129    2c Class 1 vs. Class 3 raceAfrican American/Black 1.817 [0.973, 3.393]
## 130    2c Class 1 vs. Class 3                  raceAsian 0.740 [0.381, 1.437]
## 131    2c Class 1 vs. Class 3            raceMultiracial 0.549 [0.238, 1.265]
## 132    2c Class 1 vs. Class 3                  raceOther 1.252 [0.516, 3.034]
## 133    2c Class 1 vs. Class 3   ethnicityHispanic/Latino 0.905 [0.572, 1.433]
## 134    2c Class 1 vs. Class 3     studentstatusPart Time 0.683 [0.326, 1.434]
## 135    2c Class 1 vs. Class 3         residenceOn Campus 1.233 [0.910, 1.671]
## 136    2c Class 1 vs. Class 3              greekmemGreek 1.526 [1.058, 2.201]
## 137    2c Class 1 vs. Class 3                alcmo_dich2 1.077 [0.707, 1.641]
## 138    2c Class 1 vs. Class 3                marmo_dich2 0.902 [0.639, 1.272]
## 139    2c Class 1 vs. Class 3               ecigmo_dich2 1.684 [1.159, 2.446]
## 140    2c Class 1 vs. Class 3               smokmo_dich2 1.016 [0.698, 1.480]
## 141    2c Class 1 vs. Class 3              cigarmo_dich2 1.970 [1.382, 2.810]
## 142    2c Class 2 vs. Class 3                (Intercept) 0.050 [0.030, 0.083]
## 143    2c Class 2 vs. Class 3                     age21+ 1.343 [0.866, 2.083]
## 144    2c Class 2 vs. Class 3                 genderMale 2.047 [1.341, 3.126]
## 145    2c Class 2 vs. Class 3                genderOther 3.672 [1.856, 7.263]
## 146    2c Class 2 vs. Class 3 raceAfrican American/Black 1.712 [0.702, 4.173]
## 147    2c Class 2 vs. Class 3                  raceAsian 2.998 [1.674, 5.369]
## 148    2c Class 2 vs. Class 3            raceMultiracial 2.919 [1.466, 5.813]
## 149    2c Class 2 vs. Class 3                  raceOther 0.935 [0.208, 4.213]
## 150    2c Class 2 vs. Class 3   ethnicityHispanic/Latino 1.124 [0.604, 2.094]
## 151    2c Class 2 vs. Class 3     studentstatusPart Time 1.223 [0.502, 2.976]
## 152    2c Class 2 vs. Class 3         residenceOn Campus 1.015 [0.666, 1.547]
## 153    2c Class 2 vs. Class 3              greekmemGreek 1.039 [0.551, 1.962]
## 154    2c Class 2 vs. Class 3                alcmo_dich2 0.739 [0.441, 1.238]
## 155    2c Class 2 vs. Class 3                marmo_dich2 0.899 [0.534, 1.515]
## 156    2c Class 2 vs. Class 3               ecigmo_dich2 1.292 [0.737, 2.265]
## 157    2c Class 2 vs. Class 3               smokmo_dich2 0.806 [0.434, 1.495]
## 158    2c Class 2 vs. Class 3              cigarmo_dich2 0.861 [0.447, 1.659]
## 159    2d Class 1 vs. Class 3                (Intercept) 0.062 [0.040, 0.098]
## 160    2d Class 1 vs. Class 3                     age21+ 1.193 [0.883, 1.612]
## 161    2d Class 1 vs. Class 3                 genderMale 4.373 [3.232, 5.918]
## 162    2d Class 1 vs. Class 3                genderOther 0.845 [0.295, 2.420]
## 163    2d Class 1 vs. Class 3 raceAfrican American/Black 1.774 [0.945, 3.327]
## 164    2d Class 1 vs. Class 3                  raceAsian 0.737 [0.378, 1.434]
## 165    2d Class 1 vs. Class 3            raceMultiracial 0.560 [0.243, 1.291]
## 166    2d Class 1 vs. Class 3                  raceOther 1.297 [0.531, 3.168]
## 167    2d Class 1 vs. Class 3   ethnicityHispanic/Latino 0.904 [0.569, 1.436]
## 168    2d Class 1 vs. Class 3     studentstatusPart Time 0.679 [0.321, 1.434]
## 169    2d Class 1 vs. Class 3         residenceOn Campus 1.259 [0.927, 1.709]
## 170    2d Class 1 vs. Class 3              greekmemGreek 1.491 [1.030, 2.160]
## 171    2d Class 1 vs. Class 3                alcmo_dich2 1.066 [0.698, 1.627]
## 172    2d Class 1 vs. Class 3                marmo_dich2 0.940 [0.663, 1.332]
## 173    2d Class 1 vs. Class 3               ecigmo_dich2 1.753 [1.204, 2.553]
## 174    2d Class 1 vs. Class 3               smokmo_dich2 1.049 [0.718, 1.533]
## 175    2d Class 1 vs. Class 3              cigarmo_dich2 1.887 [1.319, 2.698]
## 176    2d Class 1 vs. Class 3               rxstmo_dich2 0.933 [0.505, 1.723]
## 177    2d Class 1 vs. Class 3               rxpkmo_dich2 0.937 [0.408, 2.152]
## 178    2d Class 1 vs. Class 3              rxsedmo_dich2 0.966 [0.349, 2.676]
## 179    2d Class 1 vs. Class 3                     mhdays 0.972 [0.954, 0.990]
## 180    2d Class 2 vs. Class 3                (Intercept) 0.041 [0.023, 0.071]
## 181    2d Class 2 vs. Class 3                     age21+ 1.345 [0.865, 2.092]
## 182    2d Class 2 vs. Class 3                 genderMale 2.176 [1.415, 3.347]
## 183    2d Class 2 vs. Class 3                genderOther 3.324 [1.659, 6.660]
## 184    2d Class 2 vs. Class 3 raceAfrican American/Black 1.677 [0.685, 4.109]
## 185    2d Class 2 vs. Class 3                  raceAsian 2.834 [1.574, 5.102]
## 186    2d Class 2 vs. Class 3            raceMultiracial 2.936 [1.466, 5.882]
## 187    2d Class 2 vs. Class 3                  raceOther 0.902 [0.199, 4.087]
## 188    2d Class 2 vs. Class 3   ethnicityHispanic/Latino 1.088 [0.580, 2.039]
## 189    2d Class 2 vs. Class 3     studentstatusPart Time 1.289 [0.527, 3.152]
## 190    2d Class 2 vs. Class 3         residenceOn Campus 0.996 [0.653, 1.520]
## 191    2d Class 2 vs. Class 3              greekmemGreek 1.038 [0.547, 1.968]
## 192    2d Class 2 vs. Class 3                alcmo_dich2 0.733 [0.436, 1.232]
## 193    2d Class 2 vs. Class 3                marmo_dich2 0.860 [0.503, 1.470]
## 194    2d Class 2 vs. Class 3               ecigmo_dich2 1.236 [0.700, 2.180]
## 195    2d Class 2 vs. Class 3               smokmo_dich2 0.750 [0.396, 1.421]
## 196    2d Class 2 vs. Class 3              cigarmo_dich2 0.833 [0.424, 1.636]
## 197    2d Class 2 vs. Class 3               rxstmo_dich2 1.560 [0.590, 4.129]
## 198    2d Class 2 vs. Class 3               rxpkmo_dich2 2.131 [0.799, 5.678]
## 199    2d Class 2 vs. Class 3              rxsedmo_dich2 0.531 [0.117, 2.421]
## 200    2d Class 2 vs. Class 3                     mhdays 1.022 [1.001, 1.044]
##     p.value sig   fmi McFadden Nagelkerke
## 1    0.0000 *** 0.004    0.051      0.103
## 2    0.0804   . 0.001    0.051      0.103
## 3    0.0203   * 0.004    0.051      0.103
## 4    0.1336     0.005    0.051      0.103
## 5    0.0001 *** 0.004    0.051      0.103
## 6    0.0000 *** 0.010    0.051      0.103
## 7    0.5054     0.007    0.051      0.103
## 8    0.0092  ** 0.023    0.051      0.103
## 9    0.4413     0.082    0.051      0.103
## 10   0.3895     0.001    0.051      0.103
## 11   0.0954   . 0.001    0.051      0.103
## 12   0.0399   * 0.002    0.051      0.103
## 13   0.0000 *** 0.002    0.051      0.103
## 14   0.0004 *** 0.001    0.051      0.103
## 15   0.0000 *** 0.002    0.051      0.103
## 16   0.0427   * 0.001    0.051      0.103
## 17   0.1239     0.006    0.051      0.103
## 18   0.0000 *** 0.004    0.051      0.103
## 19   0.0311   * 0.002    0.051      0.103
## 20   0.1482     0.012    0.051      0.103
## 21   0.0594   . 0.048    0.051      0.103
## 22   0.8980     0.017    0.051      0.103
## 23   0.0632   . 0.001    0.051      0.103
## 24   0.0000 *** 0.002    0.051      0.103
## 25   0.0000 *** 0.004    0.071      0.142
## 26   0.8047     0.001    0.071      0.142
## 27   0.0499   * 0.005    0.071      0.142
## 28   0.1231     0.006    0.071      0.142
## 29   0.0006 *** 0.004    0.071      0.142
## 30   0.0000 *** 0.010    0.071      0.142
## 31   0.5201     0.007    0.071      0.142
## 32   0.0123   * 0.024    0.071      0.142
## 33   0.4788     0.079    0.071      0.142
## 34   0.4010     0.001    0.071      0.142
## 35   0.1255     0.001    0.071      0.142
## 36   0.2328     0.002    0.071      0.142
## 37   0.1773     0.013    0.071      0.142
## 38   0.8116     0.004    0.071      0.142
## 39   0.0088  ** 0.003    0.071      0.142
## 40   0.5542     0.001    0.071      0.142
## 41   0.0046  ** 0.004    0.071      0.142
## 42   0.0000 *** 0.010    0.071      0.142
## 43   0.3656     0.002    0.071      0.142
## 44   0.0000 *** 0.003    0.071      0.142
## 45   0.0328   * 0.001    0.071      0.142
## 46   0.4299     0.007    0.071      0.142
## 47   0.0000 *** 0.004    0.071      0.142
## 48   0.0495   * 0.002    0.071      0.142
## 49   0.1716     0.011    0.071      0.142
## 50   0.0883   . 0.050    0.071      0.142
## 51   0.9615     0.019    0.071      0.142
## 52   0.1555     0.002    0.071      0.142
## 53   0.0165   * 0.001    0.071      0.142
## 54   0.2199     0.021    0.071      0.142
## 55   0.1352     0.002    0.071      0.142
## 56   0.0044  ** 0.002    0.071      0.142
## 57   0.1962     0.001    0.071      0.142
## 58   0.0000 *** 0.002    0.071      0.142
## 59   0.0000 *** 0.004    0.075      0.149
## 60   0.8012     0.001    0.075      0.149
## 61   0.0351   * 0.005    0.075      0.149
## 62   0.0777   . 0.006    0.075      0.149
## 63   0.0009 *** 0.004    0.075      0.149
## 64   0.0000 *** 0.010    0.075      0.149
## 65   0.5253     0.007    0.075      0.149
## 66   0.0148   * 0.024    0.075      0.149
## 67   0.4506     0.078    0.075      0.149
## 68   0.3538     0.001    0.075      0.149
## 69   0.1405     0.001    0.075      0.149
## 70   0.2258     0.002    0.075      0.149
## 71   0.1793     0.013    0.075      0.149
## 72   0.9478     0.004    0.075      0.149
## 73   0.0126   * 0.003    0.075      0.149
## 74   0.6651     0.001    0.075      0.149
## 75   0.0067  ** 0.004    0.075      0.149
## 76   0.0762   . 0.007    0.075      0.149
## 77   0.1007     0.001    0.075      0.149
## 78   0.1383     0.001    0.075      0.149
## 79   0.1111     0.009    0.075      0.149
## 80   0.0000 *** 0.009    0.075      0.149
## 81   0.5065     0.002    0.075      0.149
## 82   0.0000 *** 0.003    0.075      0.149
## 83   0.0588   . 0.001    0.075      0.149
## 84   0.4139     0.007    0.075      0.149
## 85   0.0000 *** 0.005    0.075      0.149
## 86   0.0613   . 0.002    0.075      0.149
## 87   0.2283     0.011    0.075      0.149
## 88   0.0708   . 0.050    0.075      0.149
## 89   0.9397     0.019    0.075      0.149
## 90   0.1206     0.002    0.075      0.149
## 91   0.0433   * 0.002    0.075      0.149
## 92   0.2039     0.020    0.075      0.149
## 93   0.1281     0.002    0.075      0.149
## 94   0.0045  ** 0.002    0.075      0.149
## 95   0.2502     0.001    0.075      0.149
## 96   0.0000 *** 0.002    0.075      0.149
## 97   0.0192   * 0.003    0.075      0.149
## 98   0.3621     0.004    0.075      0.149
## 99   0.2089     0.001    0.075      0.149
## 100  0.0013  ** 0.006    0.075      0.149
## 101  0.0000 *** 0.002    0.096      0.166
## 102  0.0070  ** 0.002    0.096      0.166
## 103  0.0000 *** 0.002    0.096      0.166
## 104  0.5519     0.001    0.096      0.166
## 105  0.1661     0.011    0.096      0.166
## 106  0.0732   . 0.012    0.096      0.166
## 107  0.1429     0.002    0.096      0.166
## 108  0.5920     0.013    0.096      0.166
## 109  0.5556     0.040    0.096      0.166
## 110  0.2946     0.013    0.096      0.166
## 111  0.0616   . 0.002    0.096      0.166
## 112  0.0004 *** 0.001    0.096      0.166
## 113  0.0000 *** 0.005    0.096      0.166
## 114  0.4269     0.002    0.096      0.166
## 115  0.0009 *** 0.002    0.096      0.166
## 116  0.0002 *** 0.003    0.096      0.166
## 117  0.1851     0.010    0.096      0.166
## 118  0.0000 *** 0.004    0.096      0.166
## 119  0.0020  ** 0.005    0.096      0.166
## 120  0.9516     0.014    0.096      0.166
## 121  0.6503     0.061    0.096      0.166
## 122  0.6893     0.001    0.096      0.166
## 123  0.8944     0.002    0.096      0.166
## 124  0.8566     0.014    0.096      0.166
## 125  0.0000 *** 0.009    0.117      0.198
## 126  0.2192     0.003    0.117      0.198
## 127  0.0000 *** 0.002    0.117      0.198
## 128  0.5676     0.002    0.117      0.198
## 129  0.0609   . 0.011    0.117      0.198
## 130  0.3742     0.015    0.117      0.198
## 131  0.1594     0.003    0.117      0.198
## 132  0.6191     0.014    0.117      0.198
## 133  0.6714     0.040    0.117      0.198
## 134  0.3140     0.014    0.117      0.198
## 135  0.1759     0.002    0.117      0.198
## 136  0.0239   * 0.001    0.117      0.198
## 137  0.7293     0.017    0.117      0.198
## 138  0.5557     0.003    0.117      0.198
## 139  0.0063  ** 0.004    0.117      0.198
## 140  0.9336     0.003    0.117      0.198
## 141  0.0002 *** 0.005    0.117      0.198
## 142  0.0000 *** 0.005    0.117      0.198
## 143  0.1884     0.001    0.117      0.198
## 144  0.0009 *** 0.002    0.117      0.198
## 145  0.0002 *** 0.002    0.117      0.198
## 146  0.2372     0.010    0.117      0.198
## 147  0.0002 *** 0.005    0.117      0.198
## 148  0.0023  ** 0.005    0.117      0.198
## 149  0.9306     0.015    0.117      0.198
## 150  0.7122     0.060    0.117      0.198
## 151  0.6579     0.001    0.117      0.198
## 152  0.9458     0.002    0.117      0.198
## 153  0.9051     0.015    0.117      0.198
## 154  0.2509     0.003    0.117      0.198
## 155  0.6894     0.002    0.117      0.198
## 156  0.3711     0.008    0.117      0.198
## 157  0.4934     0.003    0.117      0.198
## 158  0.6556     0.007    0.117      0.198
## 159  0.0000 *** 0.008    0.124      0.210
## 160  0.2515     0.003    0.124      0.210
## 161  0.0000 *** 0.003    0.124      0.210
## 162  0.7536     0.002    0.124      0.210
## 163  0.0744   . 0.012    0.124      0.210
## 164  0.3685     0.016    0.124      0.210
## 165  0.1739     0.003    0.124      0.210
## 166  0.5687     0.015    0.124      0.210
## 167  0.6692     0.040    0.124      0.210
## 168  0.3104     0.014    0.124      0.210
## 169  0.1405     0.002    0.124      0.210
## 170  0.0345   * 0.001    0.124      0.210
## 171  0.7685     0.017    0.124      0.210
## 172  0.7278     0.003    0.124      0.210
## 173  0.0035  ** 0.004    0.124      0.210
## 174  0.8058     0.003    0.124      0.210
## 175  0.0005 *** 0.006    0.124      0.210
## 176  0.8250     0.005    0.124      0.210
## 177  0.8776     0.028    0.124      0.210
## 178  0.9477     0.005    0.124      0.210
## 179  0.0029  ** 0.013    0.124      0.210
## 180  0.0000 *** 0.009    0.124      0.210
## 181  0.1881     0.001    0.124      0.210
## 182  0.0004 *** 0.003    0.124      0.210
## 183  0.0007 *** 0.002    0.124      0.210
## 184  0.2582     0.011    0.124      0.210
## 185  0.0005 *** 0.006    0.124      0.210
## 186  0.0024  ** 0.005    0.124      0.210
## 187  0.8936     0.015    0.124      0.210
## 188  0.7935     0.061    0.124      0.210
## 189  0.5784     0.002    0.124      0.210
## 190  0.9847     0.002    0.124      0.210
## 191  0.9096     0.016    0.124      0.210
## 192  0.2408     0.003    0.124      0.210
## 193  0.5815     0.003    0.124      0.210
## 194  0.4652     0.009    0.124      0.210
## 195  0.3774     0.002    0.124      0.210
## 196  0.5958     0.007    0.124      0.210
## 197  0.3705     0.002    0.124      0.210
## 198  0.1306     0.004    0.124      0.210
## 199  0.4139     0.002    0.124      0.210
## 200  0.0429   * 0.016    0.124      0.210

Generate forest plots.

all_or <- bind_rows(or_1b, or_1c, or_1d, or_2b, or_2c, or_2d) |>
  filter(term != "(Intercept)", p_value < 0.05) |>
  mutate(
    term_label = term |>
      gsub("^gender", "Gender: ", x = _) |>
      gsub("^race", "Race: ", x = _) |>
      gsub("^age", "Age: ", x = _) |>
      gsub("^greekmem", "Greek Mem: ", x = _) |>
      gsub("ecigmo_dich2", "E-cig Use", x = _) |>
      gsub("cigarmo_dich2", "Cigar Use", x = _) |>
      gsub("rxstmo_dich2", "Rx Stimulant Use", x = _) |>
      gsub("mhdays", "Mental Health Days", x = _),
    series = if_else(grepl("1", model), "Model Series 1", "Model Series 2")
  )

make_forest_plot <- function(data, title) {
  data |>
    mutate(
      panel = paste0(model, " | ", contrast),
      term_label = forcats::fct_reorder(term_label, OR)
    ) |>
    ggplot(aes(x = OR, y = term_label)) +
    geom_vline(xintercept = 1,
               linetype = "dashed",
               color = "gray50") +
    geom_pointrange(aes(xmin = CI_lower, xmax = CI_upper)) +
    scale_x_log10() +
    facet_wrap(~ panel, scales = "free_y", ncol = 2) +
    labs(title = title, x = "Odds Ratio (log scale)", y = NULL) +
    theme_bw(base_size = 11) +
    theme(
      strip.text = element_text(size = 8.5),
      axis.text.y = element_text(size = 9),
      panel.spacing = unit(0.8, "cm")
    )
}

ggsave(
  "forest_plot_series1.png",
  plot = make_forest_plot(
    filter(all_or, series == "Model Series 1"),
    "Forest Plot: Model Series 1 — Significant Odds Ratios"
  ),
  width = 10,
  height = 6,
  dpi = 300,
  path = "/Users/jasonhoskin/Library/CloudStorage/OneDrive-IndianaUniversity/study-gambling"
)

ggsave(
  "forest_plot_series2.png",
  plot = make_forest_plot(
    filter(all_or, series == "Model Series 2"),
    "Forest Plot: Model Series 2 — Significant Odds Ratios"
  ),
  width = 10,
  height = 6,
  dpi = 300,
  path = "/Users/jasonhoskin/Library/CloudStorage/OneDrive-IndianaUniversity/study-gambling"
)