ICSUS LCA
Jason Hoskin
2026-05-19
- 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)
## Set seed
set.seed(123)- 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
## 1086df3 <- 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] 1124colSums(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) / Greeklapply(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- Conduct the first LCA using the
tidySEMpackage.
## Set seed
set.seed(123)
df6_lca1 <- df6 |>
dplyr::select(gmpools_dich:gmoth_dich) |>
as.data.frame()
## Run first LCA (code used for knitting)
lca1 <- readRDS("lca1_results.rds")
## Run first LCA (code used for running first time)
## lca1 <- tidySEM::mx_lca(data = df6_lca1,
## classes = 1:4,
## missing = "fiml",
## mutate(across(gmpools_dich:gcqdepress_dich, as.ordered))
## Save LCA for future knitting
## lca1 <- saveRDS(lca1, file = "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- Conduct the second latent class analysis.
random_seed <- sample(1:10000, 1)
cat("Seed used:", random_seed, "\n")## Seed used: 2463set.seed(random_seed)
## Refilter dataset to only include appropriate variables
df6_lca2 <- df6 |>
filter(if_any(gmpools_dich:gcqdepress_dich, ~ . == 2)) |>
dplyr::select(autoid | gmpools_dich:gcqdepress_dich) |>
mutate(across(gmpools_dich:gcqdepress_dich, as.ordered))
## lca2 <- tidySEM::mx_lca(data = df6_lca2[2:22],
## classes = 1:4,
## missing = "fiml")
## Save LCA for future knitting
## lca2 <- saveRDS(lca2, file = "lca2_results.rds")
## Run second LCA (code use for knitting)
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[[3]], type = "sum.posterior")
probs## $sum.posterior
## class count proportion
## 1 class1 317.5183 0.17484485
## 2 class2 120.9421 0.06659809
## 3 class3 1377.5396 0.75855705- 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
MICEand 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 + gender + 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.034133vif_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.204690vif_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.049988vif_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.248551vif_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 233table(df6_model1$gender, df6_model1$lca1_class)##
## 3 1 2
## Female (ref) 1780 744 159
## Male 722 344 280
## Other 164 55 6table(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 10table(df6_model1$ethnicity, df6_model1$lca1_class)##
## 3 1 2
## Not Hispanic/Latino (ref) 2182 968 375
## Hispanic/Latino 344 136 38table(df6_model1$studentstatus, df6_model1$lca1_class)##
## 3 1 2
## Full Time (ref) 2582 1092 425
## Part Time 98 53 17table(df6_model1$residence, df6_model1$lca1_class)##
## 3 1 2
## Off Campus (ref) 1316 524 205
## On Campus 1370 621 240table(df6_model1$greekmem, df6_model1$lca1_class)##
## 3 1 2
## No Greek (ref) 2396 992 353
## Greek 287 153 92table(df6_model1$alcmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 894 283 75
## 2 1770 850 362table(df6_model1$marmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 1816 677 208
## 2 865 467 237table(df6_model1$ecigmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 1956 709 218
## 2 721 434 227table(df6_model1$smokmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 2230 861 267
## 2 456 285 177table(df6_model1$cigarmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 2453 961 271
## 2 227 180 172table(df6_model1$rxstmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 2619 1083 399
## 2 64 58 45table(df6_model1$rxpkmo_dich, df6_model1$lca1_class)##
## 3 1 2
## 1 2639 1109 424
## 2 44 36 20table(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 57table(df6_model2$gender, df6_model2$lca2_class)##
## 3 1 2
## Female (ref) 867 84 50
## Male 451 224 54
## Other 64 4 13table(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 2table(df6_model2$ethnicity, df6_model2$lca2_class)##
## 3 1 2
## Not Hispanic/Latino (ref) 1170 261 95
## Hispanic/Latino 155 30 15table(df6_model2$studentstatus, df6_model2$lca2_class)##
## 3 1 2
## Full Time (ref) 1324 301 111
## Part Time 60 10 6table(df6_model2$residence, df6_model2$lca2_class)##
## 3 1 2
## Off Campus (ref) 648 140 53
## On Campus 738 172 64table(df6_model2$greekmem, df6_model2$lca2_class)##
## 3 1 2
## No Greek (ref) 1204 239 103
## Greek 181 73 13table(df6_model2$alcmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 342 47 38
## 2 1031 257 79table(df6_model2$marmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 812 146 75
## 2 573 166 42table(df6_model2$ecigmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 870 141 77
## 2 514 171 39table(df6_model2$smokmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 1038 181 94
## 2 349 130 23table(df6_model2$cigarmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 1159 172 100
## 2 223 138 16table(df6_model2$rxstmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 1311 282 109
## 2 71 29 8table(df6_model2$rxpkmo_dich, df6_model2$lca2_class)##
## 3 1 2
## 1 1343 298 108
## 2 43 13 8table(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.146print(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.8566print(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.6556print(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.210all_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"
)
## 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
)