ICSUS LCA Gambling: Tables and Figures

library(haven)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(future)
library(mice)

Attaching package: 'mice'

The following object is masked from 'package:stats':

    filter

The following objects are masked from 'package:base':

    cbind, rbind

library(nnet)
library(ggtext)

data_raw <- read_spss(
  "ICSUS_2025_State_18 to 25_CLN.sav"
)

df <- data_raw |>
  select(autoid, starts_with("gm"))
df <- df[, 1:14]

# Factor and label variables
df <- df |>
  mutate(across(
    .cols = starts_with("gm"),
    .fns = ~ case_when(.x == 1 ~ "None", .x > 1 ~ "Any", .default = NA)
  ))

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

df <- df |>
  mutate(across(
    .cols = starts_with("gm"),
    .fns = as.factor
  ))

## gmpools
df |>
  group_by(gmpools) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmpools     n   prop
  <fct>   <int>  <dbl>
1 Any       382  8.93 
2 None     3887 90.9  
3 <NA>        9  0.210

## gmfant
df |>
  group_by(gmfant) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmfant     n   prop
  <fct>  <int>  <dbl>
1 Any      250  5.84 
2 None    4013 93.8  
3 <NA>      15  0.351

## gmvideo
df |>
  group_by(gmvideo) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmvideo     n   prop
  <fct>   <int>  <dbl>
1 Any       490 11.5  
2 None     3782 88.4  
3 <NA>        6  0.140

## gmonsport
df |>
  group_by(gmonsport) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmonsport     n   prop
  <fct>     <int>  <dbl>
1 Any         348  8.13 
2 None       3921 91.7  
3 <NA>          9  0.210

## gmothsp
df |>
  group_by(gmothsp) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmothsp     n   prop
  <fct>   <int>  <dbl>
1 Any       129  3.02 
2 None     4136 96.7  
3 <NA>       13  0.304

## gmonline
df |>
  group_by(gmonline) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmonline     n   prop
  <fct>    <int>  <dbl>
1 Any        193  4.51 
2 None      4076 95.3  
3 <NA>         9  0.210

## gmesports
df |>
  group_by(gmesports) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmesports     n   prop
  <fct>     <int>  <dbl>
1 Any          88  2.06 
2 None       4180 97.7  
3 <NA>         10  0.234

## gmhorse
df |>
  group_by(gmhorse) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmhorse     n   prop
  <fct>   <int>  <dbl>
1 Any        90  2.10 
2 None     4175 97.6  
3 <NA>       13  0.304

## gmcards
df |>
  group_by(gmcards) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmcards     n   prop
  <fct>   <int>  <dbl>
1 Any       382  8.93 
2 None     3885 90.8  
3 <NA>       11  0.257

## gmlott
df |>
  group_by(gmlott) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmlott     n   prop
  <fct>  <int>  <dbl>
1 Any      994 23.2  
2 None    3274 76.5  
3 <NA>      10  0.234

## gmcasino
df |>
  group_by(gmcasino) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmcasino     n   prop
  <fct>    <int>  <dbl>
1 Any        334  7.81 
2 None      3930 91.9  
3 <NA>        14  0.327

## gmcharit
df |>
  group_by(gmcharit) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmcharit     n   prop
  <fct>    <int>  <dbl>
1 Any        436 10.2  
2 None      3831 89.6  
3 <NA>        11  0.257

## gmoth
df |>
  group_by(gmoth) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  gmoth     n  prop
  <fct> <int> <dbl>
1 Any      91  2.13
2 None   4132 96.6 
3 <NA>     55  1.29

remove(df, data_raw)

data_raw <- read.csv(
  "/Users/jasonhoskin/Documents/study-gambling/results_final/lca1_k3_regression_pre_imputation.csv"
)
df <- data_raw |>
  select(
    autoid,
    mhdays,
    lca_class,
    age_cat,
    ethnicity_cat,
    race_cat,
    gender_cat,
    studentstatus_cat,
    residence_cat,
    greekmem_cat,
    alcmo_dich,
    marmo_dich,
    ecigmo_dich,
    smokmo_dich,
    cigarmo_dich,
    rxstmo_dich,
    rxpkmo_dich,
    rxsedmo_dich,
    gcq_items_missing,
    gcq_count
  )

# Rename Classes
df <- df |>
  mutate(
    lca_class = case_when(
      lca_class == "Class 1" ~ 1,
      lca_class == "Class 2" ~ 2,
      lca_class == "Class 3: low/no gambling" ~ 3
    )
  )

# Ensure all missing data is marked as NA
df <- df |>
  mutate(
    age_cat = case_when(
      age_cat == "<21" ~ "<21",
      age_cat == "21+" ~ "21+",
      .default = NA
    ),
    ethnicity_cat = case_when(
      ethnicity_cat == "Not Hispanic/Latino" ~ "Not Hispanic/Latino",
      ethnicity_cat == "Hispanic/Latino" ~ "Hispanic/Latino",
      .default = NA
    ),
    race_cat = case_when(
      race_cat == "White" ~ "White",
      race_cat == "Asian" ~ "Asian",
      race_cat == "Black/African American" ~ "Black/African American",
      race_cat == "Multiracial" ~ "Multiracial",
      race_cat == "Other" ~ "Other",
      .default = NA
    ),
    gender_cat = case_when(
      gender_cat == "Male" ~ "Male",
      gender_cat == "Female" ~ "Female",
      gender_cat == "Other" ~ "Other",
      .default = NA
    ),
    studentstatus_cat = case_when(
      studentstatus_cat == "Full time" ~ "Full time",
      studentstatus_cat == "Part time" ~ "Part time",
      .default = NA
    ),
    residence_cat = case_when(
      residence_cat == "Off campus" ~ "Off campus",
      residence_cat == "On campus" ~ "On campus",
      .default = NA
    ),
    greekmem_cat = case_when(
      greekmem_cat == "No Greek membership" ~ "No Greek membership",
      greekmem_cat == "Greek membership" ~ "Greek membership",
      .default = NA
    ),
    alcmo_dich = case_when(
      alcmo_dich == 0 ~ "None",
      alcmo_dich == 1 ~ "Any",
      .default = NA
    ),
    marmo_dich = case_when(
      marmo_dich == 0 ~ "None",
      marmo_dich == 1 ~ "Any",
      .default = NA
    ),
    ecigmo_dich = case_when(
      ecigmo_dich == 0 ~ "None",
      ecigmo_dich == 1 ~ "Any",
      .default = NA
    ),
    smokmo_dich = case_when(
      smokmo_dich == 0 ~ "None",
      smokmo_dich == 1 ~ "Any",
      .default = NA
    ),
    cigarmo_dich = case_when(
      cigarmo_dich == 0 ~ "None",
      cigarmo_dich == 1 ~ "Any",
      .default = NA
    ),
    rxstmo_dich = case_when(
      rxstmo_dich == 0 ~ "None",
      rxstmo_dich == 1 ~ "Any",
      .default = NA
    ),
    rxpkmo_dich = case_when(
      rxpkmo_dich == 0 ~ "None",
      rxpkmo_dich == 1 ~ "Any",
      .default = NA
    ),
    rxsedmo_dich = case_when(
      rxsedmo_dich == 0 ~ "None",
      rxsedmo_dich == 1 ~ "Any",
      .default = NA
    )
  )

## Convert dataframe classes and set reference groups
df$lca_class <- as.factor(df$lca_class)
df$lca_class <- relevel(df$lca_class, ref = "3")

df$age_cat <- as.factor(df$age_cat)
df$age_cat <- relevel(df$age_cat, ref = "21+")

df$ethnicity_cat <- as.factor(df$ethnicity_cat)
df$ethnicity_cat <- relevel(df$ethnicity_cat, ref = "Not Hispanic/Latino")

df$race_cat <- as.factor(df$race_cat)
df$race_cat <- relevel(df$race_cat, ref = "White")

df$gender_cat <- as.factor(df$gender_cat)
df$gender_cat <- relevel(df$gender_cat, ref = "Female")

df$studentstatus_cat <- as.factor(df$studentstatus_cat)
df$studentstatus_cat <- relevel(df$studentstatus_cat, ref = "Full time")

df$residence_cat <- as.factor(df$residence_cat)
df$residence_cat <- relevel(df$residence_cat, ref = "On campus")

df$greekmem_cat <- as.factor(df$greekmem_cat)
df$greekmem_cat <- relevel(df$greekmem_cat, ref = "No Greek membership")

df$alcmo_dich <- as.factor(df$alcmo_dich)
df$alcmo_dich <- relevel(df$alcmo_dich, ref = "None")

df$marmo_dich <- as.factor(df$marmo_dich)
df$marmo_dich <- relevel(df$marmo_dich, ref = "None")

df$ecigmo_dich <- as.factor(df$ecigmo_dich)
df$ecigmo_dich <- relevel(df$ecigmo_dich, ref = "None")

df$smokmo_dich <- as.factor(df$smokmo_dich)
df$smokmo_dich <- relevel(df$smokmo_dich, ref = "None")

df$cigarmo_dich <- as.factor(df$cigarmo_dich)
df$cigarmo_dich <- relevel(df$cigarmo_dich, ref = "None")

df$rxstmo_dich <- as.factor(df$rxstmo_dich)
df$rxstmo_dich <- relevel(df$rxstmo_dich, ref = "None")

df$rxpkmo_dich <- as.factor(df$rxpkmo_dich)
df$rxpkmo_dich <- relevel(df$rxpkmo_dich, ref = "None")

df$rxsedmo_dich <- as.factor(df$rxsedmo_dich)
df$rxsedmo_dich <- relevel(df$rxsedmo_dich, ref = "None")

df$gcq_items_missing <- as.integer(df$gcq_items_missing)
df$gcq_count <- as.integer(df$gcq_count)
df$mhdays <- as.integer(df$mhdays)

Warning: NAs introduced by coercion

## Confirm classes are correct
sapply(df, class)

           autoid            mhdays         lca_class           age_cat 
        "integer"         "integer"          "factor"          "factor" 
    ethnicity_cat          race_cat        gender_cat studentstatus_cat 
         "factor"          "factor"          "factor"          "factor" 
    residence_cat      greekmem_cat        alcmo_dich        marmo_dich 
         "factor"          "factor"          "factor"          "factor" 
      ecigmo_dich       smokmo_dich      cigarmo_dich       rxstmo_dich 
         "factor"          "factor"          "factor"          "factor" 
      rxpkmo_dich      rxsedmo_dich gcq_items_missing         gcq_count 
         "factor"          "factor"         "integer"         "integer" 

# Calculate general frequencies and proportions
### Bad mental health days
df |>
  filter(!is.na(mhdays)) |>
  summarize(
    median = median(mhdays),
    mean = mean(mhdays),
    sd = sd(mhdays),
    mad = mad(mhdays),
    rangemin = min(mhdays),
    rangemax = max(mhdays),
    IQR = IQR(mhdays)
  )

  median     mean       sd   mad rangemin rangemax IQR
1      5 8.254532 8.235722 7.413        0       30  10

df |>
  filter(!is.na(mhdays)) |>
  ggplot(aes(mhdays)) +
  geom_histogram(bins = 80)

### Gambling consequences (omitted no participants with no reported gambling and NAs)
df |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count)) |>
  summarize(
    n = n(),
    median = median(gcq_count),
    mean = mean(gcq_count),
    sd = sd(gcq_count),
    mad = mad(gcq_count),
    rangemin = min(gcq_count),
    rangemax = max(gcq_count),
    IQR = IQR(gcq_count)
  )

     n median      mean        sd mad rangemin rangemax IQR
1 1765      0 0.2016997 0.8046756   0        0        8   0

df |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count) & gcq_count > 0) |>
  summarize(
    n = n(),
    median = median(gcq_count),
    mean = mean(gcq_count),
    sd = sd(gcq_count),
    mad = mad(gcq_count),
    rangemin = min(gcq_count),
    rangemax = max(gcq_count),
    IQR = IQR(gcq_count)
  )

    n median     mean       sd mad rangemin rangemax IQR
1 184      1 1.934783 1.694119   0        1        8   1

### Gambling consequences (omitted no participants with no reported gambling, NAs, and remaining 0 reported consequences)
df |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count)) |>
  ggplot(aes(gcq_count)) +
  geom_histogram(bins = 80)

df |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count) & gcq_count > 0) |>
  ggplot(aes(gcq_count)) +
  geom_histogram(bins = 50)

### Proportion of all sample who reported any gambling
sum(df$gcq_items_missing != 8) / nrow(df)

[1] 0.4177186

### Proportion of all sample who reported any gambling consequences
sum(df$gcq_count != 0, na.rm = TRUE) / nrow(df)

[1] 0.04301075

### Proportion all sample who reported any gambling behaviors who reported any gambling consequences
sum(df$gcq_count != 0, na.rm = TRUE) / sum(df$gcq_items_missing != 8)

[1] 0.1029659

## LCA class membership
df |>
  group_by(lca_class) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  lca_class     n  prop
  <fct>     <int> <dbl>
1 3          3239 75.7 
2 1           319  7.46
3 2           720 16.8 

## Age
df |>
  group_by(age_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 2 × 3
  age_cat     n  prop
  <fct>   <int> <dbl>
1 21+      1871  43.7
2 <21      2407  56.3

## Ethnicity
df |>
  group_by(ethnicity_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  ethnicity_cat           n  prop
  <fct>               <int> <dbl>
1 Not Hispanic/Latino  3525 82.4 
2 Hispanic/Latino       518 12.1 
3 <NA>                  235  5.49

## Race
df |>
  group_by(race_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 6 × 3
  race_cat                   n   prop
  <fct>                  <int>  <dbl>
1 White                   3309 77.3  
2 Asian                    379  8.86 
3 Black/African American   222  5.19 
4 Multiracial              190  4.44 
5 Other                    140  3.27 
6 <NA>                      38  0.888

## Gender
df |>
  group_by(gender_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 4 × 3
  gender_cat     n   prop
  <fct>      <int>  <dbl>
1 Female      2683 62.7  
2 Male        1346 31.5  
3 Other        225  5.26 
4 <NA>          24  0.561

## Student Status
df |>
  group_by(studentstatus_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  studentstatus_cat     n   prop
  <fct>             <int>  <dbl>
1 Full time          4099 95.8  
2 Part time           168  3.93 
3 <NA>                 11  0.257

## Residence Status
df |>
  group_by(residence_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  residence_cat     n    prop
  <fct>         <int>   <dbl>
1 On campus      2231 52.2   
2 Off campus     2045 47.8   
3 <NA>              2  0.0468

## Greek life membership
df |>
  group_by(greekmem_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  greekmem_cat            n   prop
  <fct>               <int>  <dbl>
1 No Greek membership  3741 87.4  
2 Greek membership      532 12.4  
3 <NA>                    5  0.117

## Any alcohol use in the month
df |>
  group_by(alcmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  alcmo_dich     n  prop
  <fct>      <int> <dbl>
1 None        1252 29.3 
2 Any         2982 69.7 
3 <NA>          44  1.03

## Any marijuana use in the month
df |>
  group_by(marmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  marmo_dich     n   prop
  <fct>      <int>  <dbl>
1 None        2701 63.1  
2 Any         1569 36.7  
3 <NA>           8  0.187

## Any e-cigarette use in the month
df |>
  group_by(ecigmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  ecigmo_dich     n   prop
  <fct>       <int>  <dbl>
1 None         2883 67.4  
2 Any          1382 32.3  
3 <NA>           13  0.304

## Any smoking use in the month
df |>
  group_by(smokmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  smokmo_dich     n    prop
  <fct>       <int>   <dbl>
1 None         3358 78.5   
2 Any           918 21.5   
3 <NA>            2  0.0468

## Any cigar use in the month
df |>
  group_by(cigarmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  cigarmo_dich     n   prop
  <fct>        <int>  <dbl>
1 None          3685 86.1  
2 Any            579 13.5  
3 <NA>            14  0.327

## Any prescription stimulant use in the month
df |>
  group_by(rxstmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  rxstmo_dich     n   prop
  <fct>       <int>  <dbl>
1 None         4101 95.9  
2 Any           167  3.90 
3 <NA>           10  0.234

## Any prescription painkiller use in the month
df |>
  group_by(rxpkmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  rxpkmo_dich     n   prop
  <fct>       <int>  <dbl>
1 None         4172 97.5  
2 Any           100  2.34 
3 <NA>            6  0.140

## Any prescription sedative use in the month
df |>
  group_by(rxsedmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

# A tibble: 3 × 3
  rxsedmo_dich     n   prop
  <fct>        <int>  <dbl>
1 None          4197 98.1  
2 Any             75  1.75 
3 <NA>             6  0.140

# Calculate frequencies and proportions by latent class
### Bad mental health days
df |>
  group_by(lca_class) |>
  filter(!is.na(mhdays)) |>
  summarize(
    median = median(mhdays),
    mean = mean(mhdays),
    sd = sd(mhdays),
    mad = mad(mhdays),
    rangemin = min(mhdays),
    rangemax = max(mhdays),
    IQR = IQR(mhdays)
  )

# A tibble: 3 × 8
  lca_class median  mean    sd   mad rangemin rangemax   IQR
  <fct>      <dbl> <dbl> <dbl> <dbl>    <int>    <int> <dbl>
1 3              5  8.44  8.27  7.41        0       30    12
2 1              3  5.63  7.38  4.45        0       30     7
3 2              5  8.60  8.27  5.93        0       30    10

df |>
  filter(!is.na(mhdays)) |>
  ggplot(aes(mhdays, fill = lca_class)) +
  geom_histogram(bins = 80)

### Gambling consequences (omitted participants with no reported gambling and NAs)
df |>
  group_by(lca_class) |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count)) |>
  summarize(
    n = n(),
    median = median(gcq_count),
    mean = mean(gcq_count),
    sd = sd(gcq_count),
    mad = mad(gcq_count),
    rangemin = min(gcq_count),
    rangemax = max(gcq_count),
    IQR = IQR(gcq_count)
  )

# A tibble: 3 × 9
  lca_class     n median  mean    sd   mad rangemin rangemax   IQR
  <fct>     <int>  <int> <dbl> <dbl> <dbl>    <int>    <int> <dbl>
1 3           751      0 0.104 0.516     0        0        6     0
2 1           303      0 0.422 1.18      0        0        8     0
3 2           711      0 0.211 0.836     0        0        8     0

df |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count)) |>
  ggplot(aes(gcq_count, fill = lca_class)) +
  geom_histogram(bins = 80)

### Gambling consequences (omitted participants with no reported gambling, NAs, and remaining 0 reported consequences)
df |>
  group_by(lca_class) |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count) & gcq_count > 0) |>
  summarize(
    n = n(),
    median = median(gcq_count),
    mean = mean(gcq_count),
    sd = sd(gcq_count),
    mad = mad(gcq_count),
    rangemin = min(gcq_count),
    rangemax = max(gcq_count),
    IQR = IQR(gcq_count)
  )

# A tibble: 3 × 9
  lca_class     n median  mean    sd   mad rangemin rangemax   IQR
  <fct>     <int>  <dbl> <dbl> <dbl> <dbl>    <int>    <int> <dbl>
1 3            45      1  1.73  1.29     0        1        6     1
2 1            62      1  2.06  1.87     0        1        8     1
3 2            77      1  1.95  1.76     0        1        8     1

df |>
  filter(gcq_items_missing != 8 & !is.na(gcq_count) & gcq_count > 0) |>
  ggplot(aes(gcq_count, fill = lca_class)) +
  geom_histogram(bins = 50)

### Proportion of each latent class who reported any gambling
total_latent_class_membership <- df |>
  group_by(lca_class) |>
  summarize(n())
total_latent_class_membership <- pull(total_latent_class_membership)
df |>
  group_by(lca_class) |>
  filter(gcq_items_missing != 8) |>
  summarize(n_any_gambling = n()) |>
  mutate(
    prop_latent_class = (n_any_gambling / total_latent_class_membership) * 100
  )

# A tibble: 3 × 3
  lca_class n_any_gambling prop_latent_class
  <fct>              <int>             <dbl>
1 3                    763              23.6
2 1                    309              96.9
3 2                    715              99.3

### Proportion of each latent class who reported any gambling consequences
df |>
  group_by(lca_class) |>
  filter(gcq_count != 0, na.rm = TRUE) |>
  summarize(n_any_gcq = n()) |>
  mutate(prop_latent_class = (n_any_gcq / total_latent_class_membership) * 100)

# A tibble: 3 × 3
  lca_class n_any_gcq prop_latent_class
  <fct>         <int>             <dbl>
1 3                45              1.39
2 1                62             19.4 
3 2                77             10.7 

### Proportion each latent class that reported any gambling behaviors who also reported any gambling consequences
total_latent_class_membership_gamblers <- df |>
  group_by(lca_class) |>
  filter(gcq_items_missing != 8, na.rm = TRUE) |>
  summarize(n_gamblers = n())
total_latent_class_membership_gamblers <- pull(
  total_latent_class_membership_gamblers
)
df |>
  group_by(lca_class) |>
  filter(gcq_items_missing != 8 & gcq_count != 0) |>
  summarize(n = n()) |>
  mutate(prop_latent_class = (n / total_latent_class_membership_gamblers) * 100)

# A tibble: 3 × 3
  lca_class     n prop_latent_class
  <fct>     <int>             <dbl>
1 3            45              5.90
2 1            62             20.1 
3 2            77             10.8 

## Age by latent class
df |>
  group_by(lca_class, age_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and age_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, age_cat))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 6 × 4
# Groups:   lca_class [3]
  lca_class age_cat     n  prop
  <fct>     <fct>   <int> <dbl>
1 3         21+      1351  41.7
2 3         <21      1888  58.3
3 1         21+       173  54.2
4 1         <21       146  45.8
5 2         21+       347  48.2
6 2         <21       373  51.8

## Ethnicity by latent class
df |>
  group_by(lca_class, ethnicity_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and ethnicity_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, ethnicity_cat))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 9 × 4
# Groups:   lca_class [3]
  lca_class ethnicity_cat           n  prop
  <fct>     <fct>               <int> <dbl>
1 3         Not Hispanic/Latino  2652 81.9 
2 3         Hispanic/Latino       404 12.5 
3 3         <NA>                  183  5.65
4 1         Not Hispanic/Latino   266 83.4 
5 1         Hispanic/Latino        28  8.78
6 1         <NA>                   25  7.84
7 2         Not Hispanic/Latino   607 84.3 
8 2         Hispanic/Latino        86 11.9 
9 2         <NA>                   27  3.75

## Race by latent class
df |>
  group_by(lca_class, race_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and race_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, race_cat))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 18 × 4
# Groups:   lca_class [3]
   lca_class race_cat                   n   prop
   <fct>     <fct>                  <int>  <dbl>
 1 3         White                   2426 74.9  
 2 3         Asian                    336 10.4  
 3 3         Black/African American   179  5.53 
 4 3         Multiracial              149  4.60 
 5 3         Other                    117  3.61 
 6 3         <NA>                      32  0.988
 7 1         White                    274 85.9  
 8 1         Asian                     13  4.08 
 9 1         Black/African American    18  5.64 
10 1         Multiracial                7  2.19 
11 1         Other                      6  1.88 
12 1         <NA>                       1  0.313
13 2         White                    609 84.6  
14 2         Asian                     30  4.17 
15 2         Black/African American    25  3.47 
16 2         Multiracial               34  4.72 
17 2         Other                     17  2.36 
18 2         <NA>                       5  0.694

## Gender
df |>
  group_by(lca_class, gender_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and gender_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, gender_cat))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 11 × 4
# Groups:   lca_class [3]
   lca_class gender_cat     n   prop
   <fct>     <fct>      <int>  <dbl>
 1 3         Female      2165 66.8  
 2 3         Male         862 26.6  
 3 3         Other        189  5.84 
 4 3         <NA>          23  0.710
 5 1         Female        82 25.7  
 6 1         Male         232 72.7  
 7 1         Other          5  1.57 
 8 2         Female       436 60.6  
 9 2         Male         252 35    
10 2         Other         31  4.31 
11 2         <NA>           1  0.139

## Student Status
df |>
  group_by(lca_class, studentstatus_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and studentstatus_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, studentstatus_cat))` for per-operation
  grouping (`?dplyr::dplyr_by`) instead.

# A tibble: 9 × 4
# Groups:   lca_class [3]
  lca_class studentstatus_cat     n   prop
  <fct>     <fct>             <int>  <dbl>
1 3         Full time          3109 96.0  
2 3         Part time           122  3.77 
3 3         <NA>                  8  0.247
4 1         Full time           306 95.9  
5 1         Part time            12  3.76 
6 1         <NA>                  1  0.313
7 2         Full time           684 95    
8 2         Part time            34  4.72 
9 2         <NA>                  2  0.278

## Residence Status
df |>
  group_by(lca_class, residence_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and residence_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, residence_cat))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 8 × 4
# Groups:   lca_class [3]
  lca_class residence_cat     n    prop
  <fct>     <fct>         <int>   <dbl>
1 3         On campus      1661 51.3   
2 3         Off campus     1577 48.7   
3 3         <NA>              1  0.0309
4 1         On campus       173 54.2   
5 1         Off campus      146 45.8   
6 2         On campus       397 55.1   
7 2         Off campus      322 44.7   
8 2         <NA>              1  0.139 

## Greek life membership
df |>
  group_by(lca_class, greekmem_cat) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and greekmem_cat.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, greekmem_cat))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 7 × 4
# Groups:   lca_class [3]
  lca_class greekmem_cat            n   prop
  <fct>     <fct>               <int>  <dbl>
1 3         No Greek membership  2889 89.2  
2 3         Greek membership      345 10.7  
3 3         <NA>                    5  0.154
4 1         No Greek membership   243 76.2  
5 1         Greek membership       76 23.8  
6 2         No Greek membership   609 84.6  
7 2         Greek membership      111 15.4  

## Any alcohol use in the month
df |>
  group_by(lca_class, alcmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and alcmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, alcmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 9 × 4
# Groups:   lca_class [3]
  lca_class alcmo_dich     n   prop
  <fct>     <fct>      <int>  <dbl>
1 3         None        1039 32.1  
2 3         Any         2172 67.1  
3 3         <NA>          28  0.864
4 1         None          46 14.4  
5 1         Any          265 83.1  
6 1         <NA>           8  2.51 
7 2         None         167 23.2  
8 2         Any          545 75.7  
9 2         <NA>           8  1.11 

## Any marijuana use in the month
df |>
  group_by(lca_class, marmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and marmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, marmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 8 × 4
# Groups:   lca_class [3]
  lca_class marmo_dich     n   prop
  <fct>     <fct>      <int>  <dbl>
1 3         None        2155 66.5  
2 3         Any         1078 33.3  
3 3         <NA>           6  0.185
4 1         None         145 45.5  
5 1         Any          174 54.5  
6 2         None         401 55.7  
7 2         Any          317 44.0  
8 2         <NA>           2  0.278

## Any e-cigarette use in the month
df |>
  group_by(lca_class, ecigmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and ecigmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, ecigmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 8 × 4
# Groups:   lca_class [3]
  lca_class ecigmo_dich     n   prop
  <fct>     <fct>       <int>  <dbl>
1 3         None         2314 71.4  
2 3         Any           914 28.2  
3 3         <NA>           11  0.340
4 1         None          142 44.5  
5 1         Any           177 55.5  
6 2         None          427 59.3  
7 2         Any           291 40.4  
8 2         <NA>            2  0.278

## Any smoking use in the month
df |>
  group_by(lca_class, smokmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and smokmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, smokmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 8 × 4
# Groups:   lca_class [3]
  lca_class smokmo_dich     n    prop
  <fct>     <fct>       <int>   <dbl>
1 3         None         2660 82.1   
2 3         Any           578 17.8   
3 3         <NA>            1  0.0309
4 1         None          183 57.4   
5 1         Any           135 42.3   
6 1         <NA>            1  0.313 
7 2         None          515 71.5   
8 2         Any           205 28.5   

## Any cigar use in the month
df |>
  group_by(lca_class, cigarmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and cigarmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, cigarmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 9 × 4
# Groups:   lca_class [3]
  lca_class cigarmo_dich     n   prop
  <fct>     <fct>        <int>  <dbl>
1 3         None          2929 90.4  
2 3         Any            299  9.23 
3 3         <NA>            11  0.340
4 1         None           177 55.5  
5 1         Any            140 43.9  
6 1         <NA>             2  0.627
7 2         None           579 80.4  
8 2         Any            140 19.4  
9 2         <NA>             1  0.139

## Any prescription stimulant use in the month
df |>
  group_by(lca_class, rxstmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and rxstmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, rxstmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 9 × 4
# Groups:   lca_class [3]
  lca_class rxstmo_dich     n   prop
  <fct>     <fct>       <int>  <dbl>
1 3         None         3153 97.3  
2 3         Any            81  2.50 
3 3         <NA>            5  0.154
4 1         None          283 88.7  
5 1         Any            35 11.0  
6 1         <NA>            1  0.313
7 2         None          665 92.4  
8 2         Any            51  7.08 
9 2         <NA>            4  0.556

## Any prescription painkiller use in the month
df |>
  group_by(lca_class, rxpkmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and rxpkmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, rxpkmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 9 × 4
# Groups:   lca_class [3]
  lca_class rxpkmo_dich     n   prop
  <fct>     <fct>       <int>  <dbl>
1 3         None         3183 98.3  
2 3         Any            52  1.61 
3 3         <NA>            4  0.123
4 1         None          302 94.7  
5 1         Any            16  5.02 
6 1         <NA>            1  0.313
7 2         None          687 95.4  
8 2         Any            32  4.44 
9 2         <NA>            1  0.139

## Any prescription sedative use in the month
df |>
  group_by(lca_class, rxsedmo_dich) |>
  summarize(n = n()) |>
  mutate(prop = (n / sum(n)) * 100)

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by lca_class and rxsedmo_dich.
ℹ Output is grouped by lca_class.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(lca_class, rxsedmo_dich))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 7 × 4
# Groups:   lca_class [3]
  lca_class rxsedmo_dich     n   prop
  <fct>     <fct>        <int>  <dbl>
1 3         None          3189 98.5  
2 3         Any             44  1.36 
3 3         <NA>             6  0.185
4 1         None           310 97.2  
5 1         Any              9  2.82 
6 2         None           698 96.9  
7 2         Any             22  3.06 

## Compare class membership proportions with theoretical class memberships
class_membership <- df |>
  group_by(lca_class) |>
  summarize(n = n()) |>
  mutate(sample_proportion = round(n / sum(n), digits = 3))
est_proportion <- data.frame(
  lca_class = as.factor(1:3),
  est_proportion = c(0.076, 0.249, 0.675)
)
prob <- dplyr::left_join(class_membership, est_proportion)

Joining with `by = join_by(lca_class)`

prob <- prob |>
  mutate(
    diff = est_proportion - sample_proportion,
    lca_class_integer = c(3, 1, 2)
  )
prob <- arrange(prob, lca_class_integer)
prob$diff <- round(prob$diff, digits = 3)
prob |>
  select(lca_class, n, sample_proportion, est_proportion, diff)

# A tibble: 3 × 5
  lca_class     n sample_proportion est_proportion   diff
  <fct>     <int>             <dbl>          <dbl>  <dbl>
1 1           319             0.075          0.076  0.001
2 2           720             0.168          0.249  0.081
3 3          3239             0.757          0.675 -0.082

# Create LCA profile plot
lca <- data.frame(
  class = as.factor(c(rep(1, 13), rep(2, 13), rep(3, 13))),
  item = as.factor(c(rep(seq(1, 13, 1), 3))),
  cond_prob = c(
    0.78,
    0.66,
    0.33,
    0.75,
    0.36,
    0.31,
    0.11,
    0.10,
    0.39,
    0.49,
    0.33,
    0.23,
    0.11,
    0.11,
    0.03,
    0.22,
    0.08,
    0.01,
    0.09,
    0.04,
    0.05,
    0.20,
    0.58,
    0.19,
    0.28,
    0.05,
    0.01,
    0.00,
    0.05,
    0.01,
    0.00,
    0.00,
    0.00,
    0.00,
    0.01,
    0.08,
    0.01,
    0.02,
    0.00
  )
)
plot_base <- ggplot(data = lca, mapping = aes(x = item, y = cond_prob)) +
  geom_col(aes(fill = class), position = "dodge")
plot_labs <- plot_base +
  labs(x = "Gambling Type", y = "Conditional Probability")
plot_legend <- plot_labs +
  scale_fill_hue(
    labels = c(
      "**Class 1**<br>Broad sports/event and multi-activity<br>Estimated prevalence: 7.6%",
      "**Class 2**<br>Lottery-centered moderate<br>Estimated prevalence: 24.9%",
      "**Class 3**<br>Low/no gambling<br>Estimated prevalence: 67.5%"
    )
  ) +
  theme(
    legend.position = "left",
    legend.title = element_blank(),
    legend.key.spacing.y = unit(1, 'cm'),
    legend.text = element_markdown()
  )
plot_axis <- plot_legend +
  scale_x_discrete(
    labels = c(
      "Sports pools",
      "Fantasy sports",
      "Video game wagering",
      "Non-sports events",
      "Other sports",
      "Online gambling",
      "E-sports",
      "Horse racing",
      "Card games",
      "Lottery",
      "Casino",
      "Charitable gambling",
      "Other gambling"
    )
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_axis

plan(multisession, workers = parallel::detectCores() - 1)
set.seed(123)

# Define variables
outcome <- "lca_class"
predictors <- c(
  "mhdays",
  "age_cat",
  "ethnicity_cat",
  "race_cat",
  "gender_cat",
  "studentstatus_cat",
  "residence_cat",
  "greekmem_cat",
  "alcmo_dich",
  "marmo_dich",
  "ecigmo_dich",
  "smokmo_dich",
  "cigarmo_dich",
  "rxstmo_dich",
  "rxpkmo_dich",
  "rxsedmo_dich",
  "gcq_count"
)
all_vars <- c(outcome, predictors)

# Establish imputation method
meth <- make.method(df[, all_vars])
meth[] <- "cart"

# Trimmed predictor matrix
pred <- quickpred(df[, all_vars], mincor = 0.1, minpuc = 0.1)

# Impute in parallel
imp <- futuremice(
  df[, all_vars],
  method = meth,
  predictorMatrix = pred,
  m = 20,
  maxit = 10,
  parallelseed = 42
)

## Evaluate imputations
plot(imp)

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

# Fit model 1: null model
fit1 <- with(imp, multinom(lca_class ~ 1, trace = FALSE))

# Pool model 1
pooled1 <- pool(fit1)
summary(pooled1)

         term  estimate  std.error statistic       df       p.value
1 (Intercept) -2.317828 0.05868162 -39.49837 4274.001 3.770743e-291
2 (Intercept) -1.503768 0.04120226 -36.49724 4274.001 4.092233e-254

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

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

# Fit model 2: covariate model
fit2 <- with(
  imp,
  multinom(
    lca_class ~ age_cat +
      gender_cat +
      race_cat +
      ethnicity_cat +
      studentstatus_cat +
      residence_cat +
      greekmem_cat,
    trace = FALSE
  )
)

# Pool model 2
pooled2 <- pool(fit2)
summary(pooled2)

                             term    estimate  std.error   statistic       df
1                     (Intercept) -2.91104072 0.13893687 -20.9522555 4234.427
2                      age_cat<21 -0.46869064 0.13174721  -3.5574996 4243.048
3                  gender_catMale  2.03581180 0.13691191  14.8695013 4247.904
4                 gender_catOther -0.28216028 0.46884752  -0.6018167 4250.701
5                   race_catAsian -1.41102464 0.29775169  -4.7389307 3868.101
6  race_catBlack/African American -0.06480914 0.27084137  -0.2392882 4242.315
7             race_catMultiracial -0.81099068 0.40459032  -2.0044738 4240.897
8                   race_catOther -0.70307137 0.45087338  -1.5593543 4158.815
9    ethnicity_catHispanic/Latino -0.21759305 0.22184513  -0.9808331 1205.353
10     studentstatus_catPart time -0.09365494 0.32753757  -0.2859365 4140.213
11        residence_catOff campus -0.29991340 0.14029115  -2.1377927 4247.204
12   greekmem_catGreek membership  0.93195586 0.16206820   5.7503930 4243.151
13                    (Intercept) -1.34876069 0.07937811 -16.9915953 4245.660
14                     age_cat<21 -0.23248489 0.08848287  -2.6274564 4249.216
15                 gender_catMale  0.43604600 0.09021694   4.8333052 4237.950
16                gender_catOther -0.20421789 0.20213106  -1.0103241 4236.162
17                  race_catAsian -1.12018032 0.19802498  -5.6567628 4118.087
18 race_catBlack/African American -0.54421678 0.21966392  -2.4774974 3944.813
19            race_catMultiracial -0.08406774 0.19861420  -0.4232716 4203.275
20                  race_catOther -0.48739166 0.27692830  -1.7599923 3758.570
21   ethnicity_catHispanic/Latino -0.07129572 0.13374731  -0.5330629 3845.336
22     studentstatus_catPart time  0.12965935 0.20338426   0.6375093 4227.734
23        residence_catOff campus -0.17152432 0.09240007  -1.8563223 4241.293
24   greekmem_catGreek membership  0.41495073 0.12353055   3.3590940 4246.364
        p.value
1  7.917009e-93
2  3.784944e-04
3  8.605090e-49
4  5.473282e-01
5  2.225724e-06
6  8.108937e-01
7  4.508280e-02
8  1.189886e-01
9  3.268719e-01
10 7.749411e-01
11 3.259050e-02
12 9.529412e-09
13 1.073113e-62
14 8.633502e-03
15 1.390373e-06
16 3.123977e-01
17 1.647234e-08
18 1.327231e-02
19 6.721188e-01
20 7.849045e-02
21 5.940209e-01
22 5.238277e-01
23 6.347693e-02
24 7.888772e-04

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

                             term          fmi       lambda
1                     (Intercept) 0.0030849454 0.0026141947
2                      age_cat<21 0.0020321016 0.0015618110
3                  gender_catMale 0.0012852819 0.0008151773
4                 gender_catOther 0.0007574639 0.0002874202
5                   race_catAsian 0.0194115025 0.0189046206
6  race_catBlack/African American 0.0021325210 0.0016621966
7             race_catMultiracial 0.0023199377 0.0018495443
8                   race_catOther 0.0084810530 0.0080043400
9    ethnicity_catHispanic/Latino 0.1052952228 0.1038119016
10     studentstatus_catPart time 0.0094370752 0.0089586826
11        residence_catOff campus 0.0014036733 0.0009335469
12   greekmem_catGreek membership 0.0020177260 0.0015474401
13                    (Intercept) 0.0016505895 0.0011804085
14                     age_cat<21 0.0010499579 0.0005798877
15                 gender_catMale 0.0026843928 0.0022138444
16                gender_catOther 0.0028918577 0.0024212087
17                  race_catAsian 0.0104874883 0.0100070359
18 race_catBlack/African American 0.0170181731 0.0165199327
19            race_catMultiracial 0.0057664552 0.0052934921
20                  race_catOther 0.0225633284 0.0220433558
21   ethnicity_catHispanic/Latino 0.0200883962 0.0195788663
22     studentstatus_catPart time 0.0037679389 0.0032967663
23        residence_catOff campus 0.0022684378 0.0017980641
24   greekmem_catGreek membership 0.0015403687 0.0010702136

# Fit model 3: covariates + alcohol + marijuana + electronic vapor products + smoking + cigars
fit3 <- with(
  imp,
  nnet::multinom(
    lca_class ~ age_cat +
      gender_cat +
      race_cat +
      ethnicity_cat +
      studentstatus_cat +
      residence_cat +
      greekmem_cat +
      alcmo_dich +
      marmo_dich +
      ecigmo_dich +
      smokmo_dich +
      cigarmo_dich,
    trace = FALSE
  )
)

# Pool model 3
pooled3 <- pool(fit3)
summary(pooled3)

                             term    estimate  std.error   statistic       df
1                     (Intercept) -3.80029240 0.22513803 -16.8798333 4031.020
2                      age_cat<21 -0.14523696 0.14221040  -1.0212822 4231.295
3                  gender_catMale  1.91045093 0.14663797  13.0283511 4231.571
4                 gender_catOther -0.30414846 0.47169446  -0.6447997 4238.829
5                   race_catAsian -1.01973842 0.30718765  -3.3195944 3807.735
6  race_catBlack/African American  0.13690439 0.28143100   0.4864581 4232.956
7             race_catMultiracial -0.75020302 0.40947002  -1.8321318 4228.630
8                   race_catOther -0.69139753 0.46507723  -1.4866295 4158.820
9    ethnicity_catHispanic/Latino -0.18436139 0.22777479  -0.8094021 1343.062
10     studentstatus_catPart time -0.05556481 0.33467344  -0.1660270 4127.719
11        residence_catOff campus -0.23897560 0.14569710  -1.6402221 4235.237
12   greekmem_catGreek membership  0.60576808 0.17245113   3.5126942 4234.147
13                  alcmo_dichAny  0.28619014 0.19936064   1.4355398 3795.203
14                  marmo_dichAny  0.08350941 0.16357919   0.5105136 4236.224
15                 ecigmo_dichAny  0.67484841 0.17666501   3.8199325 4236.109
16                 smokmo_dichAny  0.15426231 0.17786916   0.8672797 4228.953
17                cigarmo_dichAny  0.84434934 0.16937564   4.9850695 4228.588
18                    (Intercept) -1.65775402 0.12472720 -13.2910384 4192.823
19                     age_cat<21 -0.11356871 0.09468804  -1.1993987 4234.724
20                 gender_catMale  0.40302816 0.09483783   4.2496560 4222.296
21                gender_catOther -0.22109558 0.20350915  -1.0864159 4225.866
22                  race_catAsian -0.97775041 0.20059261  -4.8743092 4109.368
23 race_catBlack/African American -0.46062392 0.22122610  -2.0821410 3939.413
24            race_catMultiracial -0.07457687 0.19976184  -0.3733289 4199.165
25                  race_catOther -0.48138704 0.27874817  -1.7269604 3764.309
26   ethnicity_catHispanic/Latino -0.05809635 0.13451641  -0.4318904 3857.993
27     studentstatus_catPart time  0.13271596 0.20449262   0.6490012 4214.150
28        residence_catOff campus -0.15767992 0.09361463  -1.6843513 4231.826
29   greekmem_catGreek membership  0.30455171 0.12645188   2.4084396 4237.555
30                  alcmo_dichAny  0.06052670 0.11661286   0.5190397 4096.294
31                  marmo_dichAny  0.09499402 0.10965591   0.8662919 4201.689
32                 ecigmo_dichAny  0.22886962 0.11901632   1.9230104 4172.118
33                 smokmo_dichAny  0.15351897 0.12686985   1.2100508 4233.421
34                cigarmo_dichAny  0.38300883 0.13394062   2.8595419 4233.989
        p.value
1  8.038255e-62
2  3.071791e-01
3  4.518147e-38
4  5.190919e-01
5  9.099655e-04
6  6.266676e-01
7  6.700209e-02
8  1.371885e-01
9  4.184273e-01
10 8.681439e-01
11 1.010333e-01
12 4.482215e-04
13 1.512157e-01
14 6.097183e-01
15 1.353998e-04
16 3.858380e-01
17 6.441930e-07
18 1.627564e-39
19 2.304401e-01
20 2.187191e-05
21 2.773570e-01
22 1.133109e-06
23 3.739398e-02
24 7.089225e-01
25 8.425681e-02
26 6.658452e-01
27 5.163729e-01
28 9.218753e-02
29 1.606331e-02
30 6.037611e-01
31 3.863796e-01
32 5.454673e-02
33 2.263270e-01
34 4.263284e-03

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

                             term         fmi       lambda
1                     (Intercept) 0.013698117 0.0132088826
2                      age_cat<21 0.002273580 0.0018020975
3                  gender_catMale 0.002237200 0.0017657306
4                 gender_catOther 0.001124035 0.0006528479
5                   race_catAsian 0.020944658 0.0204305476
6  race_catBlack/African American 0.002049703 0.0015783000
7             race_catMultiracial 0.002609144 0.0021375234
8                   race_catOther 0.007950791 0.0074738236
9    ethnicity_catHispanic/Latino 0.097222303 0.0958789464
10     studentstatus_catPart time 0.009582637 0.0091028670
11        residence_catOff campus 0.001719198 0.0012478925
12   greekmem_catGreek membership 0.001880834 0.0014094832
13                  alcmo_dichAny 0.021305238 0.0207896198
14                  marmo_dichAny 0.001566134 0.0010948662
15                 ecigmo_dichAny 0.001584346 0.0011130744
16                 smokmo_dichAny 0.002569821 0.0020982176
17                cigarmo_dichAny 0.002614154 0.0021425307
18                    (Intercept) 0.005812364 0.0053382438
19                     age_cat<21 0.001796171 0.0013248452
20                 gender_catMale 0.003319993 0.0028480015
21                gender_catOther 0.002932248 0.0024604712
22                  race_catAsian 0.010453834 0.0099723458
23 race_catBlack/African American 0.016909523 0.0164105446
24            race_catMultiracial 0.005347462 0.0048738361
25                  race_catOther 0.022180836 0.0216614524
26   ethnicity_catHispanic/Latino 0.019462952 0.0189547692
27     studentstatus_catPart time 0.004111611 0.0036390832
28        residence_catOff campus 0.002203296 0.0017318397
29   greekmem_catGreek membership 0.001348035 0.0008768126
30                  alcmo_dichAny 0.011042888 0.0105601515
31                  marmo_dichAny 0.005154177 0.0046807444
32                 ecigmo_dichAny 0.007171463 0.0066956424
33                 smokmo_dichAny 0.001984606 0.0015132243
34                cigarmo_dichAny 0.001903646 0.0014322890

# Fit model 4: model c + rx stimulants + rx painkillers + rx sedatives + mhdays
fit4 <- with(
  imp,
  nnet::multinom(
    lca_class ~ age_cat +
      gender_cat +
      race_cat +
      ethnicity_cat +
      studentstatus_cat +
      residence_cat +
      greekmem_cat +
      alcmo_dich +
      marmo_dich +
      ecigmo_dich +
      smokmo_dich +
      cigarmo_dich +
      rxstmo_dich +
      rxpkmo_dich +
      rxsedmo_dich +
      mhdays +
      gcq_count,
    trace = FALSE
  )
)

# Pool model 4
pooled4 <- pool(fit4)
summary(pooled4)

                             term     estimate   std.error   statistic       df
1                     (Intercept) -3.677409125 0.240783717 -15.2726653 4071.490
2                      age_cat<21 -0.074643704 0.146580630  -0.5092331 4217.780
3                  gender_catMale  1.775977727 0.151309903  11.7373529 4195.055
4                 gender_catOther -0.255263870 0.484916822  -0.5264075 4215.600
5                   race_catAsian -1.339031800 0.336285141  -3.9818346 3425.427
6  race_catBlack/African American  0.071608959 0.291442290   0.2457054 4121.409
7             race_catMultiracial -0.810782723 0.418080587  -1.9392977 4202.607
8                   race_catOther -0.650917530 0.482528845  -1.3489712 4090.947
9    ethnicity_catHispanic/Latino -0.251951372 0.231871179  -1.0866006 1746.853
10     studentstatus_catPart time -0.097811699 0.348858164  -0.2803767 4143.346
11        residence_catOff campus -0.285110055 0.149585564  -1.9059998 4206.896
12   greekmem_catGreek membership  0.532052398 0.179377819   2.9660992 4210.219
13                  alcmo_dichAny  0.368302500 0.206647966   1.7822701 3697.646
14                  marmo_dichAny  0.108016407 0.170153910   0.6348159 4219.044
15                 ecigmo_dichAny  0.678181609 0.182026962   3.7257206 4214.823
16                 smokmo_dichAny  0.147987191 0.184994742   0.7999535 4198.341
17                cigarmo_dichAny  0.801323351 0.174945153   4.5804261 4214.585
18                 rxstmo_dichAny  0.556824030 0.297018296   1.8747129 4197.820
19                 rxpkmo_dichAny  0.293001559 0.405186670   0.7231273 3756.986
20                rxsedmo_dichAny -0.516316927 0.508386390  -1.0155994 4202.587
21                         mhdays -0.038360301 0.009886238  -3.8801717 3488.765
22                      gcq_count  1.357002658 0.145177193   9.3472165 3527.110
23                    (Intercept) -1.722070084 0.135031003 -12.7531459 4126.559
24                     age_cat<21 -0.091722485 0.096467675  -0.9508106 4218.638
25                 gender_catMale  0.369809749 0.097643988   3.7873274 4181.573
26                gender_catOther -0.302853546 0.210108991  -1.4414116 4196.421
27                  race_catAsian -1.148459410 0.211631750  -5.4266877 4077.236
28 race_catBlack/African American -0.506706116 0.226976590  -2.2324158 3885.021
29            race_catMultiracial -0.139463160 0.204651675  -0.6814660 4210.933
30                  race_catOther -0.436642861 0.280997650  -1.5539022 3826.868
31   ethnicity_catHispanic/Latino -0.089660635 0.137245881  -0.6532847 3921.655
32     studentstatus_catPart time  0.136123134 0.208660112   0.6523678 4197.988
33        residence_catOff campus -0.161412644 0.095279486  -1.6940965 4203.529
34   greekmem_catGreek membership  0.301628349 0.129149238   2.3355024 4217.221
35                  alcmo_dichAny  0.078278404 0.119002147   0.6577898 4048.110
36                  marmo_dichAny  0.058686614 0.113348990   0.5177515 4166.790
37                 ecigmo_dichAny  0.210375370 0.121786111   1.7274168 4092.734
38                 smokmo_dichAny  0.127003605 0.130277515   0.9748697 4217.831
39                cigarmo_dichAny  0.341077384 0.137182861   2.4862973 4223.140
40                 rxstmo_dichAny  0.500187742 0.226245048   2.2108229 4162.297
41                 rxpkmo_dichAny  0.492892419 0.285964969   1.7236112 4203.786
42                rxsedmo_dichAny -0.090450531 0.335190839  -0.2698479 4221.247
43                         mhdays -0.001377652 0.005443280  -0.2530921 3395.315
44                      gcq_count  1.140693745 0.135834482   8.3976743 3899.151
        p.value
1  2.988577e-51
2  6.106155e-01
3  2.521105e-31
4  5.986328e-01
5  6.980532e-05
6  8.059225e-01
7  5.253194e-02
8  1.774209e-01
9  2.773633e-01
10 7.792025e-01
11 5.671832e-02
12 3.032980e-03
13 7.478722e-02
14 5.255829e-01
15 1.973085e-04
16 4.237830e-01
17 4.774274e-06
18 6.090161e-02
19 4.696467e-01
20 3.098786e-01
21 1.063109e-04
22 1.551049e-20
23 1.455030e-36
24 3.417550e-01
25 1.544151e-04
26 1.495430e-01
27 6.074473e-08
28 2.564427e-02
29 4.956141e-01
30 1.202904e-01
31 5.136111e-01
32 5.141997e-01
33 9.032103e-02
34 1.956384e-02
35 5.107106e-01
36 6.046591e-01
37 8.416829e-02
38 3.296808e-01
39 1.294637e-02
40 2.710228e-02
41 8.485159e-02
42 7.872905e-01
43 8.002123e-01
44 6.296536e-17

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

                             term         fmi      lambda
1                     (Intercept) 0.011711324 0.011225976
2                      age_cat<21 0.002717347 0.002244565
3                  gender_catMale 0.004899713 0.004425410
4                 gender_catOther 0.002969251 0.002496344
5                   race_catAsian 0.030846964 0.030281271
6  race_catBlack/African American 0.009422791 0.008942209
7             race_catMultiracial 0.004260205 0.003786450
8                   race_catOther 0.010862048 0.010378592
9    ethnicity_catHispanic/Latino 0.078577294 0.077522946
10     studentstatus_catPart time 0.008275143 0.007796552
11        residence_catOff campus 0.003865981 0.003392522
12   greekmem_catGreek membership 0.003541123 0.003067883
13                  alcmo_dichAny 0.023804648 0.023276782
14                  marmo_dichAny 0.002564660 0.002091947
15                 ecigmo_dichAny 0.003055920 0.002582967
16                 smokmo_dichAny 0.004629005 0.004154945
17                cigarmo_dichAny 0.003082216 0.002609248
18                 rxstmo_dichAny 0.004672623 0.004198525
19                 rxpkmo_dichAny 0.022161462 0.021641057
20                rxsedmo_dichAny 0.004262042 0.003788286
21                         mhdays 0.029257980 0.028701643
22                      gcq_count 0.028284902 0.027734061
23                    (Intercept) 0.009163101 0.008682993
24                     age_cat<21 0.002614309 0.002141573
25                 gender_catMale 0.005914881 0.005439534
26                gender_catOther 0.004788513 0.004314312
27                  race_catAsian 0.011465339 0.010980554
28 race_catBlack/African American 0.018374008 0.017868799
29            race_catMultiracial 0.003468780 0.002995585
30                  race_catOther 0.020143375 0.019631416
31   ethnicity_catHispanic/Latino 0.017205099 0.016704012
32     studentstatus_catPart time 0.004658588 0.004184501
33        residence_catOff campus 0.004177660 0.003703970
34   greekmem_catGreek membership 0.002783196 0.002310382
35                  alcmo_dichAny 0.012675830 0.012188155
36                  marmo_dichAny 0.006898702 0.006422142
37                 ecigmo_dichAny 0.010781612 0.010298328
38                 smokmo_dichAny 0.002711327 0.002238547
39                cigarmo_dichAny 0.002028798 0.001556289
40                 rxstmo_dichAny 0.007177641 0.006700700
41                 rxpkmo_dichAny 0.004154378 0.003680705
42                rxsedmo_dichAny 0.002285228 0.001812629
43                         mhdays 0.031595924 0.031025657
44                      gcq_count 0.017928717 0.017425111

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_1 <- compute_or_table(pooled1, "Model 1")
or_2 <- compute_or_table(pooled2, "Model 2")
or_3 <- compute_or_table(pooled3, "Model 3")
or_4 <- compute_or_table(pooled4, "Model 4")

all_or <- bind_rows(
  "Model 1" = or_1,
  "Model 2" = or_2,
  "Model 3" = or_3,
  "Model 4" = or_4,
)
all_or

      model            contrast                           term    OR CI_lower
1   Model 1 Class 1 vs. Class 3                    (Intercept) 0.098    0.088
2   Model 1 Class 2 vs. Class 3                    (Intercept) 0.222    0.205
3   Model 2 Class 1 vs. Class 3                    (Intercept) 0.054    0.041
4   Model 2 Class 1 vs. Class 3                     age_cat<21 0.626    0.483
5   Model 2 Class 1 vs. Class 3                 gender_catMale 7.658    5.856
6   Model 2 Class 1 vs. Class 3                gender_catOther 0.754    0.301
7   Model 2 Class 1 vs. Class 3                  race_catAsian 0.244    0.136
8   Model 2 Class 1 vs. Class 3 race_catBlack/African American 0.937    0.551
9   Model 2 Class 1 vs. Class 3            race_catMultiracial 0.444    0.201
10  Model 2 Class 1 vs. Class 3                  race_catOther 0.495    0.205
11  Model 2 Class 1 vs. Class 3   ethnicity_catHispanic/Latino 0.804    0.521
12  Model 2 Class 1 vs. Class 3     studentstatus_catPart time 0.911    0.479
13  Model 2 Class 1 vs. Class 3        residence_catOff campus 0.741    0.563
14  Model 2 Class 1 vs. Class 3   greekmem_catGreek membership 2.539    1.848
15  Model 2 Class 2 vs. Class 3                    (Intercept) 0.260    0.222
16  Model 2 Class 2 vs. Class 3                     age_cat<21 0.793    0.666
17  Model 2 Class 2 vs. Class 3                 gender_catMale 1.547    1.296
18  Model 2 Class 2 vs. Class 3                gender_catOther 0.815    0.549
19  Model 2 Class 2 vs. Class 3                  race_catAsian 0.326    0.221
20  Model 2 Class 2 vs. Class 3 race_catBlack/African American 0.580    0.377
21  Model 2 Class 2 vs. Class 3            race_catMultiracial 0.919    0.623
22  Model 2 Class 2 vs. Class 3                  race_catOther 0.614    0.357
23  Model 2 Class 2 vs. Class 3   ethnicity_catHispanic/Latino 0.931    0.716
24  Model 2 Class 2 vs. Class 3     studentstatus_catPart time 1.138    0.764
25  Model 2 Class 2 vs. Class 3        residence_catOff campus 0.842    0.703
26  Model 2 Class 2 vs. Class 3   greekmem_catGreek membership 1.514    1.189
27  Model 3 Class 1 vs. Class 3                    (Intercept) 0.022    0.014
28  Model 3 Class 1 vs. Class 3                     age_cat<21 0.865    0.654
29  Model 3 Class 1 vs. Class 3                 gender_catMale 6.756    5.068
30  Model 3 Class 1 vs. Class 3                gender_catOther 0.738    0.293
31  Model 3 Class 1 vs. Class 3                  race_catAsian 0.361    0.198
32  Model 3 Class 1 vs. Class 3 race_catBlack/African American 1.147    0.661
33  Model 3 Class 1 vs. Class 3            race_catMultiracial 0.472    0.212
34  Model 3 Class 1 vs. Class 3                  race_catOther 0.501    0.201
35  Model 3 Class 1 vs. Class 3   ethnicity_catHispanic/Latino 0.832    0.532
36  Model 3 Class 1 vs. Class 3     studentstatus_catPart time 0.946    0.491
37  Model 3 Class 1 vs. Class 3        residence_catOff campus 0.787    0.592
38  Model 3 Class 1 vs. Class 3   greekmem_catGreek membership 1.833    1.307
39  Model 3 Class 1 vs. Class 3                  alcmo_dichAny 1.331    0.901
40  Model 3 Class 1 vs. Class 3                  marmo_dichAny 1.087    0.789
41  Model 3 Class 1 vs. Class 3                 ecigmo_dichAny 1.964    1.389
42  Model 3 Class 1 vs. Class 3                 smokmo_dichAny 1.167    0.823
43  Model 3 Class 1 vs. Class 3                cigarmo_dichAny 2.326    1.669
44  Model 3 Class 2 vs. Class 3                    (Intercept) 0.191    0.149
45  Model 3 Class 2 vs. Class 3                     age_cat<21 0.893    0.741
46  Model 3 Class 2 vs. Class 3                 gender_catMale 1.496    1.243
47  Model 3 Class 2 vs. Class 3                gender_catOther 0.802    0.538
48  Model 3 Class 2 vs. Class 3                  race_catAsian 0.376    0.254
49  Model 3 Class 2 vs. Class 3 race_catBlack/African American 0.631    0.409
50  Model 3 Class 2 vs. Class 3            race_catMultiracial 0.928    0.627
51  Model 3 Class 2 vs. Class 3                  race_catOther 0.618    0.358
52  Model 3 Class 2 vs. Class 3   ethnicity_catHispanic/Latino 0.944    0.725
53  Model 3 Class 2 vs. Class 3     studentstatus_catPart time 1.142    0.765
54  Model 3 Class 2 vs. Class 3        residence_catOff campus 0.854    0.711
55  Model 3 Class 2 vs. Class 3   greekmem_catGreek membership 1.356    1.058
56  Model 3 Class 2 vs. Class 3                  alcmo_dichAny 1.062    0.845
57  Model 3 Class 2 vs. Class 3                  marmo_dichAny 1.100    0.887
58  Model 3 Class 2 vs. Class 3                 ecigmo_dichAny 1.257    0.996
59  Model 3 Class 2 vs. Class 3                 smokmo_dichAny 1.166    0.909
60  Model 3 Class 2 vs. Class 3                cigarmo_dichAny 1.467    1.128
61  Model 4 Class 1 vs. Class 3                    (Intercept) 0.025    0.016
62  Model 4 Class 1 vs. Class 3                     age_cat<21 0.928    0.696
63  Model 4 Class 1 vs. Class 3                 gender_catMale 5.906    4.390
64  Model 4 Class 1 vs. Class 3                gender_catOther 0.775    0.299
65  Model 4 Class 1 vs. Class 3                  race_catAsian 0.262    0.136
66  Model 4 Class 1 vs. Class 3 race_catBlack/African American 1.074    0.607
67  Model 4 Class 1 vs. Class 3            race_catMultiracial 0.445    0.196
68  Model 4 Class 1 vs. Class 3                  race_catOther 0.522    0.203
69  Model 4 Class 1 vs. Class 3   ethnicity_catHispanic/Latino 0.777    0.493
70  Model 4 Class 1 vs. Class 3     studentstatus_catPart time 0.907    0.458
71  Model 4 Class 1 vs. Class 3        residence_catOff campus 0.752    0.561
72  Model 4 Class 1 vs. Class 3   greekmem_catGreek membership 1.702    1.198
73  Model 4 Class 1 vs. Class 3                  alcmo_dichAny 1.445    0.964
74  Model 4 Class 1 vs. Class 3                  marmo_dichAny 1.114    0.798
75  Model 4 Class 1 vs. Class 3                 ecigmo_dichAny 1.970    1.379
76  Model 4 Class 1 vs. Class 3                 smokmo_dichAny 1.159    0.807
77  Model 4 Class 1 vs. Class 3                cigarmo_dichAny 2.228    1.582
78  Model 4 Class 1 vs. Class 3                 rxstmo_dichAny 1.745    0.975
79  Model 4 Class 1 vs. Class 3                 rxpkmo_dichAny 1.340    0.606
80  Model 4 Class 1 vs. Class 3                rxsedmo_dichAny 0.597    0.220
81  Model 4 Class 1 vs. Class 3                         mhdays 0.962    0.944
82  Model 4 Class 1 vs. Class 3                      gcq_count 3.885    2.923
83  Model 4 Class 2 vs. Class 3                    (Intercept) 0.179    0.137
84  Model 4 Class 2 vs. Class 3                     age_cat<21 0.912    0.755
85  Model 4 Class 2 vs. Class 3                 gender_catMale 1.447    1.195
86  Model 4 Class 2 vs. Class 3                gender_catOther 0.739    0.489
87  Model 4 Class 2 vs. Class 3                  race_catAsian 0.317    0.209
88  Model 4 Class 2 vs. Class 3 race_catBlack/African American 0.602    0.386
89  Model 4 Class 2 vs. Class 3            race_catMultiracial 0.870    0.582
90  Model 4 Class 2 vs. Class 3                  race_catOther 0.646    0.373
91  Model 4 Class 2 vs. Class 3   ethnicity_catHispanic/Latino 0.914    0.699
92  Model 4 Class 2 vs. Class 3     studentstatus_catPart time 1.146    0.761
93  Model 4 Class 2 vs. Class 3        residence_catOff campus 0.851    0.706
94  Model 4 Class 2 vs. Class 3   greekmem_catGreek membership 1.352    1.050
95  Model 4 Class 2 vs. Class 3                  alcmo_dichAny 1.081    0.856
96  Model 4 Class 2 vs. Class 3                  marmo_dichAny 1.060    0.849
97  Model 4 Class 2 vs. Class 3                 ecigmo_dichAny 1.234    0.972
98  Model 4 Class 2 vs. Class 3                 smokmo_dichAny 1.135    0.880
99  Model 4 Class 2 vs. Class 3                cigarmo_dichAny 1.406    1.075
100 Model 4 Class 2 vs. Class 3                 rxstmo_dichAny 1.649    1.058
101 Model 4 Class 2 vs. Class 3                 rxpkmo_dichAny 1.637    0.935
102 Model 4 Class 2 vs. Class 3                rxsedmo_dichAny 0.914    0.474
103 Model 4 Class 2 vs. Class 3                         mhdays 0.999    0.988
104 Model 4 Class 2 vs. Class 3                      gcq_count 3.129    2.398
    CI_upper p_value sig
1      0.110  0.0000 ***
2      0.241  0.0000 ***
3      0.071  0.0000 ***
4      0.810  0.0004 ***
5     10.016  0.0000 ***
6      1.890  0.5473    
7      0.437  0.0000 ***
8      1.594  0.8109    
9      0.982  0.0451   *
10     1.198  0.1190    
11     1.243  0.3269    
12     1.730  0.7749    
13     0.975  0.0326   *
14     3.489  0.0000 ***
15     0.303  0.0000 ***
16     0.943  0.0086  **
17     1.846  0.0000 ***
18     1.212  0.3124    
19     0.481  0.0000 ***
20     0.893  0.0133   *
21     1.357  0.6721    
22     1.057  0.0785   .
23     1.210  0.5940    
24     1.696  0.5238    
25     1.010  0.0635   .
26     1.929  0.0008 ***
27     0.035  0.0000 ***
28     1.143  0.3072    
29     9.006  0.0000 ***
30     1.860  0.5191    
31     0.659  0.0009 ***
32     1.991  0.6267    
33     1.054  0.0670   .
34     1.246  0.1372    
35     1.300  0.4184    
36     1.823  0.8681    
37     1.048  0.1010    
38     2.570  0.0004 ***
39     1.968  0.1512    
40     1.498  0.6097    
41     2.776  0.0001 ***
42     1.653  0.3858    
43     3.242  0.0000 ***
44     0.243  0.0000 ***
45     1.075  0.2304    
46     1.802  0.0000 ***
47     1.195  0.2774    
48     0.557  0.0000 ***
49     0.973  0.0374   *
50     1.373  0.7089    
51     1.067  0.0843   .
52     1.228  0.6658    
53     1.705  0.5164    
54     1.026  0.0922   .
55     1.737  0.0161   *
56     1.335  0.6038    
57     1.363  0.3864    
58     1.587  0.0545   .
59     1.495  0.2263    
60     1.907  0.0043  **
61     0.041  0.0000 ***
62     1.237  0.6106    
63     7.945  0.0000 ***
64     2.004  0.5986    
65     0.507  0.0001 ***
66     1.902  0.8059    
67     1.009  0.0525   .
68     1.343  0.1774    
69     1.224  0.2774    
70     1.797  0.7792    
71     1.008  0.0567   .
72     2.420  0.0030  **
73     2.167  0.0748   .
74     1.555  0.5256    
75     2.815  0.0002 ***
76     1.666  0.4238    
77     3.140  0.0000 ***
78     3.124  0.0609   .
79     2.966  0.4696    
80     1.616  0.3099    
81     0.981  0.0001 ***
82     5.163  0.0000 ***
83     0.233  0.0000 ***
84     1.102  0.3418    
85     1.753  0.0002 ***
86     1.115  0.1495    
87     0.480  0.0000 ***
88     0.940  0.0256   *
89     1.299  0.4956    
90     1.121  0.1203    
91     1.196  0.5136    
92     1.725  0.5142    
93     1.026  0.0903   .
94     1.742  0.0196   *
95     1.366  0.5107    
96     1.324  0.6047    
97     1.567  0.0842   .
98     1.466  0.3297    
99     1.840  0.0129   *
100    2.569  0.0271   *
101    2.867  0.0849   .
102    1.762  0.7873    
103    1.009  0.8002    
104    4.083  0.0000 ***

# Make forest plot for final model
or_4$labels <- dplyr::recode(
  or_4$term,
  "(Intercept)" = "(Intercept)",
  "age_cat<21" = "Age: <21",
  "gender_catMale" = "Gender: Male",
  "gender_catOther" = "Gender: Other",
  "race_catAsian" = "Race: Asian",
  "race_catBlack/African American" = "Race: Black/African American",
  "race_catMultiracial" = "Race: Multiracial",
  "race_catOther" = "Race: Other",
  "ethnicity_catHispanic/Latino" = "Ethnicity: Hispanic/Latino",
  "studentstatus_catPart time" = "Student Status: Part Time",
  "residence_catOff campus" = "Residence: Off Campus",
  "greekmem_catGreek membership" = "Greek Life Membership",
  "alcmo_dichAny" = "Alcohol Use: Any",
  "marmo_dichAny" = "Marijuana Use: Any",
  "ecigmo_dichAny" = "E-Cigarette Use: Any",
  "smokmo_dichAny" = "Smoking Use: Any",
  "cigarmo_dichAny" = "Cigar Use: Any",
  "rxstmo_dichAny" = "Prescription Stimulant Use: Any",
  "rxpkmo_dichAny" = "Prescription Painkillers: Any",
  "rxsedmo_dichAny" = "Prescription Sedative Use: Any",
  "mhdays" = "Number of Poor Mental Health Days",
  "gcq_count" = "Number of Gambling Consequences"
)

or_4 |>
  mutate(labels = fct_reorder(labels, desc(OR))) |>
  ggplot(aes(x = OR, y = labels, color = contrast)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray50") +
  geom_pointrange(
    aes(xmin = CI_lower, xmax = CI_upper),
    position = position_dodge(width = 0.5)
  ) +
  ## scale_x_log10() +
  labs(
    title = "Forest Plot of Odds Ratios",
    x = "Odds Ratio",
    y = NULL,
    color = "Outcome Level"
  ) +
  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")
  )