- Install packages and load dataset.
## R markdown link: http://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/spring-2026/research-in-health-behavior"
)
## 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
## 1086
df3 <- df2 |>
filter(!if_all(
c(
gmpools,
gmfant,
gmvideo,
gmonsport,
gmothsp,
gmonline,
gmesports,
gmhorse,
gmcards,
gmlott,
gmcasino,
gmcharit,
gmoth
),
is.na
))
## Number of participants removed from dataset:
dim(df2)[1] - dim(df3)[1]
## [1] 1124
colSums(is.na(df3))
## autoid age ethnicity race gender
## 0 0 235 38 24
## studentstatus residence greekmem smokmo cigarmo
## 11 2 5 2 14
## snufmo hookmo ecigmo alcmo marmo
## 9 11 13 44 8
## cocmo hallumo hermo methmo inhmo
## 6 8 9 7 6
## rxstmo rxpkmo rxsedmo othmo gmpools
## 10 6 6 8 9
## gmfant gmvideo gmonsport gmothsp gmonline
## 15 6 9 13 9
## gmesports gmhorse gmcards gmlott gmcasino
## 10 13 11 10 14
## gmcharit gmoth gcqsleep gcqhyg gcqfrnd
## 11 55 2492 2492 2491
## gcqfmly gcqacad gcqmoney gcqbad gcqdepress
## 2493 2495 2494 2495 2504
## mhdays
## 34
- Discriminate skip logic NAs from missing data NAs.
## Create a gambling participation indicator
gm_vars <- names(df3)[25:37] # gmpools:gmoth
df3$any_gamble <- apply(df3[gm_vars], 1, function(x) {
if (all(is.na(x)))
return(NA)
if (all(x == 1, na.rm = TRUE))
return(0) # never gambled
return(1) # gambled at least once
})
## Create a function to separate skip logic NAs from missing data NAs
classify_gcq_na <- function(data, gate_var, consequence_vars) {
gate <- data[[gate_var]]
results <- lapply(consequence_vars, function(var) {
item <- data[[var]]
skip_na <- sum(is.na(item) & gate == 0, na.rm = TRUE)
true_na <- sum(is.na(item) & gate == 1, na.rm = TRUE)
observed <- sum(!is.na(item))
data.frame(
variable = var,
skip_logic_na = skip_na,
true_missing_na = true_na,
observed = observed
)
})
do.call(rbind, results)
}
## Apply function to dataset
gcq_vars <- names(df3)[38:45] # gcqsleep:gcqdepress
na_summary <- classify_gcq_na(data = df3,
gate_var = "any_gamble",
consequence_vars = gcq_vars)
na_summary
## variable skip_logic_na true_missing_na observed
## 1 gcqsleep 2462 30 1786
## 2 gcqhyg 2462 30 1786
## 3 gcqfrnd 2462 29 1787
## 4 gcqfmly 2462 31 1785
## 5 gcqacad 2462 33 1783
## 6 gcqmoney 2462 32 1784
## 7 gcqbad 2462 33 1783
## 8 gcqdepress 2462 42 1774
- Make the variables dichotomous and convert to class “factor”.
df4 <- df3 |>
mutate(across(
.cols = starts_with("gm"),
.fns = ~ case_when(.x == 1 ~ 1, .x > 1 ~ 2, .default = NA),
.names = "{.col}_dich"
))
df5 <- df4 |>
mutate(across(
.cols = starts_with("gcq"),
.fns = ~ case_when(.x == 1 ~ 1, .x > 1 ~ 2, .default = NA),
.names = "{.col}_dich"
))
df6 <- df5 |>
mutate(across(
.cols = ends_with("mo"),
.fns = ~ case_when(.x == 1 ~ 1, .x > 1 ~ 2, .default = NA),
.names = "{.col}_dich"
))
## Filter out unnecessary variables
df6 <- df6 |>
dplyr::select(autoid:greekmem | gmpools_dich:othmo_dich | mhdays)
# Convert substances and gambling to factor with 1 as reference.
cols_to_factor <- colnames(df6)[9:45]
df6[cols_to_factor] <- lapply(df6[cols_to_factor], function(x) {
x <- as.factor(x)
x <- relevel(x, ref = "1")
return(x)
})
## Convert demographics to class "factor"
df6$age <- factor(df6$age,
levels = c(0, 1),
labels = c("<21 (ref)", "21+"))
df6$ethnicity <- factor(
df6$ethnicity,
levels = c(0, 1),
labels = c("Not Hispanic/Latino (ref)", "Hispanic/Latino")
)
df6$race <- factor(
df6$race,
levels = c(0, 1, 2, 3, 4),
labels = c(
"White (ref)",
"African American/Black",
"Asian",
"Multiracial",
"Other"
)
)
df6$gender <- factor(
df6$gender,
levels = c(0, 1, 2),
labels = c("Female (ref)", "Male", "Other")
)
df6$studentstatus <- factor(
df6$studentstatus,
levels = c(0, 1),
labels = c("Full Time (ref)", "Part Time")
)
df6$residence <- factor(
df6$residence,
levels = c(0, 1),
labels = c("Off Campus (ref)", "On Campus")
)
df6$greekmem <- factor(df6$greekmem,
levels = c(0, 1),
labels = c("No Greek (ref)", "Greek"))
# Reference group verification
demo_vars <- c("age",
"ethnicity",
"race",
"gender",
"studentstatus",
"residence",
"greekmem")
for (v in demo_vars) {
cat(sprintf(
"%-16s | Reference: %-30s | Levels: %s\n",
v,
levels(df6[[v]])[1],
paste(levels(df6[[v]]), collapse = " / ")
))
}
## age | Reference: <21 (ref) | Levels: <21 (ref) / 21+
## ethnicity | Reference: Not Hispanic/Latino (ref) | Levels: Not Hispanic/Latino (ref) / Hispanic/Latino
## race | Reference: White (ref) | Levels: White (ref) / African American/Black / Asian / Multiracial / Other
## gender | Reference: Female (ref) | Levels: Female (ref) / Male / Other
## studentstatus | Reference: Full Time (ref) | Levels: Full Time (ref) / Part Time
## residence | Reference: Off Campus (ref) | Levels: Off Campus (ref) / On Campus
## greekmem | Reference: No Greek (ref) | Levels: No Greek (ref) / Greek
lapply(df6[demo_vars], table, useNA = "ifany")
## $age
##
## <21 (ref) 21+
## 2407 1871
##
## $ethnicity
##
## Not Hispanic/Latino (ref) Hispanic/Latino <NA>
## 3525 518 235
##
## $race
##
## White (ref) African American/Black Asian
## 3309 222 379
## Multiracial Other <NA>
## 190 140 38
##
## $gender
##
## Female (ref) Male Other <NA>
## 2683 1346 225 24
##
## $studentstatus
##
## Full Time (ref) Part Time <NA>
## 4099 168 11
##
## $residence
##
## Off Campus (ref) On Campus <NA>
## 2045 2231 2
##
## $greekmem
##
## No Greek (ref) Greek <NA>
## 3741 532 5
- Conduct the first LCA using the
tidySEM package.
## 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")
## 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.
set.seed(123)
## 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)
## Run second LCA (code use for knitting)
lca2 <- readRDS("lca2_results.rds")
## 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")
- 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 -10766.31 1816 43 21618.61 21855.30 21718.69 0.8706046
## 3 3 3 -10407.67 1816 65 20945.34 21303.12 21096.62 0.7554996
## 4 4 4 -10137.62 1816 87 20449.24 20928.12 20651.72 0.7906488
## 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.8765436 0.9804656 0.15583700 0.8441630 42.23256 13.476190
## 3 0.8232718 0.9289432 0.18447137 0.6200441 27.93846 15.952381
## 4 0.8263384 0.9779214 0.01817181 0.4614537 20.87356 1.571429
# 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 = 1177.133, LMR LR (df = 22) = 1127.071, 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 = 717.276, LMR LR (df = 22) = 686.771, 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 = 540.100, LMR LR (df = 22) = 517.129, 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] 335
# Minimum sample size for 4 classes
fit2[4, 4] * fit2[4, 12]
## [1] 33
- 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>
## 355 335 1126 2462
# Should show counts for classes 1, 2, 3 with 2462 NAs
- Prepare dataframes for regression modeling.
## Filter out any unnecessary variables
df6_model1 <- df6 |>
dplyr::select(autoid:greekmem | smokmo_dich:lca1_class | mhdays)
## Filter out non-gamblers to only include participants in 2nd LCA
df6_model2 <- df6 |>
filter(!is.na(lca2_class)) |>
dplyr::select(autoid:greekmem |
smokmo_dich:othmo_dich | lca2_class | mhdays)
## Make DVs class "factor"
df6_model1 <- df6_model1 |>
mutate(lca1_class = as.factor(lca1_class))
df6_model2 <- df6_model2 |>
mutate(lca2_class = as.factor(lca2_class))
### 4. Relevel DVs so no gambling is the reference group
df6_model1$lca1_class <- relevel(df6_model1$lca1_class, ref = 3)
df6_model2$lca2_class <- relevel(df6_model2$lca2_class, ref = 3)
- Treat missing data using
MICE and evaluate the
imputations.
plan(multisession, workers = parallel::detectCores() - 1)
# Strip all haven labels across the whole data frame
df6_model1 <- haven::zap_labels(df6_model1)
df6_model1 <- haven::zap_label(df6_model1) # variable labels
df6_model1 <- haven::zap_labels(df6_model1) # value labels
# Strip all haven labels across the whole data frame
df6_model2 <- haven::zap_labels(df6_model2)
df6_model2 <- haven::zap_label(df6_model2) # variable labels
df6_model2 <- haven::zap_labels(df6_model2) # value labels
# Define variables
outcome1 <- "lca1_class"
predictors <- c(
"age",
"ethnicity",
"race",
"gender",
"studentstatus",
"residence",
"greekmem",
"smokmo_dich",
"cigarmo_dich",
"snufmo_dich",
"hookmo_dich",
"ecigmo_dich",
"alcmo_dich",
"marmo_dich",
"cocmo_dich",
"hallumo_dich",
"hermo_dich",
"methmo_dich",
"inhmo_dich",
"rxstmo_dich",
"rxpkmo_dich",
"rxsedmo_dich",
"othmo_dich",
"mhdays"
)
all_vars1 <- c(outcome1, predictors)
# Establish imputation method
meth1 <- make.method(df6_model1[, all_vars1])
meth1[] <- "cart"
# Trimmed predictor matrix
pred1 <- quickpred(df6_model1[, all_vars1], mincor = 0.1, minpuc = 0.1)
# Impute in parallel
imp1 <- futuremice(
df6_model1[, all_vars1],
method = meth1,
predictorMatrix = pred1,
m = 20,
maxit = 10,
parallelseed = 42
)
## Evaluate imputations
plot(imp1)
stripplot(imp1, pch = 20, cex = 1.2)
# Define variables for model 2
outcome2 <- "lca2_class"
all_vars2 <- c(outcome2, predictors)
# Create imputation method for method 2
meth2 <- make.method(df6_model2[, all_vars2])
meth2[] <- "cart"
# Trimmed predictor matrix
pred2 <- quickpred(df6_model2[, all_vars2], mincor = 0.1, minpuc = 0.1)
# Impute in parallel
imp2 <- futuremice(
df6_model2[, all_vars2],
method = meth2,
predictorMatrix = pred2,
m = 20,
maxit = 25,
parallelseed = 42
)
## Evaluate imputations
plot(imp2)
stripplot(imp2, pch = 20, cex = 1.2)
- Run multinomial logistic regression models.
# Fit model 1a: null model
fit1a <- with(imp1, multinom(lca1_class ~ 1, trace = FALSE))
# Pool model 1a
pooled1a <- pool(fit1a)
summary(pooled1a)
## term estimate std.error statistic df p.value
## 1 (Intercept) -0.8521477 0.03528119 -24.15303 4274.001 6.31298e-121
## 2 (Intercept) -1.7981063 0.05117962 -35.13325 4274.001 8.75161e-238
## View fraction of missing information (FMI)
data.frame(
term = summary(pooled1a)$term,
fmi = pooled1a$pooled$fmi,
lambda = pooled1a$pooled$lambda
)
## term fmi lambda
## 1 (Intercept) 0.0004676173 0
## 2 (Intercept) 0.0004676173 0
# Fit model 1b: covariate model
fit1b <- with(
imp1,
multinom(
lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
trace = FALSE
)
)
# Pool model 1b
pooled1b <- pool(fit1b)
summary(pooled1b)
## term estimate std.error statistic df
## 1 (Intercept) -0.90589370 0.06895316 -13.1378134 4222.197
## 2 age21+ 0.13306035 0.07609206 1.7486760 4250.994
## 3 genderMale 0.18582202 0.08006103 2.3210046 4225.336
## 4 genderOther -0.24475314 0.16314979 -1.5001744 4210.950
## 5 raceAfrican American/Black -0.69905149 0.18305327 -3.8188419 4225.060
## 6 raceAsian -0.99866581 0.15136043 -6.5979317 4132.906
## 7 raceMultiracial -0.11295103 0.16959683 -0.6659973 4178.481
## 8 raceOther -0.60972873 0.23403602 -2.6052773 3729.321
## 9 ethnicityHispanic/Latino -0.09066452 0.11770694 -0.7702564 1673.410
## 10 studentstatusPart Time 0.15355311 0.17843027 0.8605777 4249.292
## 11 residenceOn Campus 0.13080162 0.07842130 1.6679348 4247.587
## 12 greekmemGreek 0.22880358 0.11132463 2.0552827 4240.974
## 13 (Intercept) -2.65619141 0.12052671 -22.0381977 4241.826
## 14 age21+ 0.40536820 0.11374212 3.5639234 4247.007
## 15 genderMale 1.53422031 0.11085642 13.8397064 4243.662
## 16 genderOther -0.86215574 0.42535208 -2.0269226 4251.041
## 17 raceAfrican American/Black -0.38196541 0.24817418 -1.5391022 4192.615
## 18 raceAsian -1.43758928 0.24934764 -5.7654017 4227.241
## 19 raceMultiracial -0.70892306 0.32876012 -2.1563536 4242.650
## 20 raceOther -0.51982702 0.35942112 -1.4462896 4091.741
## 21 ethnicityHispanic/Latino -0.35837174 0.19004298 -1.8857405 2760.591
## 22 studentstatusPart Time -0.03624918 0.28284965 -0.1281571 3957.880
## 23 residenceOn Campus 0.22305855 0.12004676 1.8580972 4247.832
## 24 greekmemGreek 0.69567410 0.14634098 4.7537887 4246.278
## p.value
## 1 1.135958e-38
## 2 8.041922e-02
## 3 2.033395e-02
## 4 1.336442e-01
## 5 1.360027e-04
## 6 4.696323e-11
## 7 5.054496e-01
## 8 9.216466e-03
## 9 4.412565e-01
## 10 3.895192e-01
## 11 9.540238e-02
## 12 3.991281e-02
## 13 5.370301e-102
## 14 3.693699e-04
## 15 1.223502e-42
## 16 4.273262e-02
## 17 1.238548e-01
## 18 8.727182e-09
## 19 3.111209e-02
## 20 1.481726e-01
## 21 5.943484e-02
## 22 8.980312e-01
## 23 6.322435e-02
## 24 2.062493e-06
## View fraction of missing information (FMI)
data.frame(
term = summary(pooled1b)$term,
fmi = pooled1b$pooled$fmi,
lambda = pooled1b$pooled$lambda
)
## term fmi lambda
## 1 (Intercept) 0.0042763522 0.0038048025
## 2 age21+ 0.0006954888 0.0002254484
## 3 genderMale 0.0039934525 0.0035221192
## 4 genderOther 0.0051989819 0.0047266111
## 5 raceAfrican American/Black 0.0040188480 0.0035474959
## 6 raceAsian 0.0097932627 0.0093141967
## 7 raceMultiracial 0.0073731378 0.0068981377
## 8 raceOther 0.0233669097 0.0228432909
## 9 ethnicityHispanic/Latino 0.0816817745 0.0805848888
## 10 studentstatusPart Time 0.0010357771 0.0005657086
## 11 residenceOn Campus 0.0013394518 0.0008693376
## 12 greekmemGreek 0.0023099407 0.0018395512
## 13 (Intercept) 0.0021981255 0.0017277778
## 14 age21+ 0.0014362132 0.0009660804
## 15 genderMale 0.0019459690 0.0014757058
## 16 genderOther 0.0006856007 0.0002155608
## 17 raceAfrican American/Black 0.0064932222 0.0060194034
## 18 raceAsian 0.0038150296 0.0033438243
## 19 raceMultiracial 0.0020869910 0.0016166821
## 20 raceOther 0.0116444383 0.0111614585
## 21 ethnicityHispanic/Latino 0.0475351126 0.0468453183
## 22 studentstatusPart Time 0.0165896313 0.0160928189
## 23 residenceOn Campus 0.0012977249 0.0008276182
## 24 greekmemGreek 0.0015540090 0.0010838510
# Fit model 1c: covariates + alcohol + marijuana + electronic vapor products + smoking + cigars
fit1c <- with(
imp1,
nnet::multinom(
lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
trace = FALSE
)
)
# Pool model 1c
pooled1c <- pool(fit1c)
summary(pooled1c)
## term estimate std.error statistic df
## 1 (Intercept) -1.09154408 0.08989229 -12.14280030 4218.020
## 2 age21+ 0.02011902 0.08136534 0.24726775 4238.161
## 3 genderMale 0.16396216 0.08358490 1.96162414 4197.964
## 4 genderOther -0.25366598 0.16448201 -1.54221103 4192.297
## 5 raceAfrican American/Black -0.62903913 0.18428640 -3.41337784 4217.341
## 6 raceAsian -0.87261799 0.15350840 -5.68449656 4127.941
## 7 raceMultiracial -0.10969586 0.17055078 -0.64318591 4175.164
## 8 raceOther -0.59016747 0.23558213 -2.50514535 3700.793
## 9 ethnicityHispanic/Latino -0.08377707 0.11825246 -0.70845941 1727.078
## 10 studentstatusPart Time 0.15065425 0.17938153 0.83985376 4237.961
## 11 residenceOn Campus 0.12152907 0.07931347 1.53226270 4238.151
## 12 greekmemGreek 0.13550127 0.11353846 1.19343939 4230.245
## 13 alcmo_dich2 0.13055184 0.09674665 1.34941967 4051.547
## 14 marmo_dich2 0.02249661 0.09435374 0.23842842 4217.679
## 15 ecigmo_dich2 0.26832379 0.10238481 2.62073827 4225.562
## 16 smokmo_dich2 0.06647428 0.11237047 0.59156356 4240.042
## 17 cigarmo_dich2 0.35381440 0.12484946 2.83392810 4214.268
## 18 (Intercept) -3.06352545 0.16267745 -18.83190017 4112.636
## 19 age21+ 0.11132027 0.12303231 0.90480521 4232.897
## 20 genderMale 1.39008094 0.11932325 11.64970713 4226.966
## 21 genderOther -0.91393693 0.42802923 -2.13522081 4238.844
## 22 raceAfrican American/Black -0.20185349 0.25570957 -0.78938572 4181.855
## 23 raceAsian -1.07772902 0.25609393 -4.20833492 4210.264
## 24 raceMultiracial -0.65374507 0.33271511 -1.96487943 4233.022
## 25 raceOther -0.50693987 0.37071110 -1.36747962 4098.519
## 26 ethnicityHispanic/Latino -0.33344840 0.19557964 -1.70492387 2642.131
## 27 studentstatusPart Time -0.01395012 0.28929991 -0.04822027 3868.349
## 28 residenceOn Campus 0.17664522 0.12434734 1.42057902 4235.915
## 29 greekmemGreek 0.37281996 0.15549838 2.39758104 4236.964
## 30 alcmo_dich2 0.20147127 0.16420064 1.22698224 3793.614
## 31 marmo_dich2 0.21243362 0.14218671 1.49404687 4231.948
## 32 ecigmo_dich2 0.43949356 0.15407969 2.85237833 4231.190
## 33 smokmo_dich2 0.20409079 0.15788542 1.29265129 4238.424
## 34 cigarmo_dich2 0.91068264 0.15338533 5.93722130 4234.560
## p.value
## 1 2.245069e-33
## 2 8.047130e-01
## 3 4.987225e-02
## 4 1.230978e-01
## 5 6.476743e-04
## 6 1.402563e-08
## 7 5.201388e-01
## 8 1.228275e-02
## 9 4.787556e-01
## 10 4.010378e-01
## 11 1.255323e-01
## 12 2.327643e-01
## 13 1.772776e-01
## 14 8.115604e-01
## 15 8.805379e-03
## 16 5.541744e-01
## 17 4.619838e-03
## 18 5.964358e-76
## 19 3.656201e-01
## 20 6.795938e-31
## 21 3.280015e-02
## 22 4.299313e-01
## 23 2.626218e-05
## 24 4.949350e-02
## 25 1.715500e-01
## 26 8.832618e-02
## 27 9.615432e-01
## 28 1.555128e-01
## 29 1.654676e-02
## 30 2.199054e-01
## 31 1.352379e-01
## 32 4.360388e-03
## 33 1.962022e-01
## 34 3.131136e-09
## View fraction of missing information (FMI)
data.frame(
term = summary(pooled1c)$term,
fmi = pooled1c$pooled$fmi,
lambda = pooled1c$pooled$lambda
)
## term fmi lambda
## 1 (Intercept) 0.0037492849 0.0032770185
## 2 age21+ 0.0012435431 0.0007723383
## 3 genderMale 0.0054376743 0.0049639564
## 4 genderOther 0.0058497010 0.0053755394
## 5 raceAfrican American/Black 0.0038145353 0.0033422238
## 6 raceAsian 0.0095717584 0.0090920091
## 7 raceMultiracial 0.0069839693 0.0065084056
## 8 raceOther 0.0239299904 0.0234026406
## 9 ethnicityHispanic/Latino 0.0794323276 0.0783669036
## 10 studentstatusPart Time 0.0012784717 0.0008072612
## 11 residenceOn Campus 0.0012453741 0.0007741690
## 12 greekmemGreek 0.0024089342 0.0019373987
## 13 alcmo_dich2 0.0129059248 0.0124187772
## 14 marmo_dich2 0.0037821189 0.0033098299
## 15 ecigmo_dich2 0.0029665253 0.0024947312
## 16 smokmo_dich2 0.0008933556 0.0004221946
## 17 cigarmo_dich2 0.0041008830 0.0036283631
## 18 (Intercept) 0.0103027270 0.0098215482
## 19 age21+ 0.0020578822 0.0015864763
## 20 genderMale 0.0028063544 0.0023346412
## 21 genderOther 0.0011212891 0.0006501025
## 22 raceAfrican American/Black 0.0065583159 0.0060833094
## 23 raceAsian 0.0044548514 0.0039820502
## 24 raceMultiracial 0.0020405103 0.0015691100
## 25 raceOther 0.0109443159 0.0104617931
## 26 ethnicityHispanic/Latino 0.0504932405 0.0497747692
## 27 studentstatusPart Time 0.0191497039 0.0186427192
## 28 residenceOn Campus 0.0016148245 0.0011435452
## 29 greekmemGreek 0.0014468288 0.0009755868
## 30 alcmo_dich2 0.0213507027 0.0208348931
## 31 marmo_dich2 0.0021870322 0.0017155816
## 32 ecigmo_dich2 0.0022872570 0.0018157693
## 33 smokmo_dich2 0.0011970616 0.0007258641
## 34 cigarmo_dich2 0.0018203883 0.0013490553
# Fit model 1d: model c + rx stimulants + rx painkillers + rx sedatives + mhdays
fit1d <- with(
imp1,
nnet::multinom(
lca1_class ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
trace = FALSE
)
)
# Pool model 1d
pooled1d <- pool(fit1d)
summary(pooled1d)
## term estimate std.error statistic df
## 1 (Intercept) -1.148945052 0.097721667 -11.75732149 4209.515
## 2 age21+ 0.020530723 0.081515972 0.25186135 4230.182
## 3 genderMale 0.177819171 0.084351817 2.10806570 4194.805
## 4 genderOther -0.293691002 0.166425336 -1.76470127 4183.700
## 5 raceAfrican American/Black -0.612081521 0.184500297 -3.31750967 4210.346
## 6 raceAsian -0.875214766 0.153753540 -5.69232272 4117.890
## 7 raceMultiracial -0.108520354 0.170826869 -0.63526514 4167.518
## 8 raceOther -0.575159715 0.235922732 -2.43791563 3688.826
## 9 ethnicityHispanic/Latino -0.089359049 0.118421566 -0.75458426 1763.913
## 10 studentstatusPart Time 0.166567013 0.179623233 0.92731330 4229.316
## 11 residenceOn Campus 0.117083211 0.079419660 1.47423460 4230.306
## 12 greekmemGreek 0.138256260 0.114123612 1.21146061 4222.191
## 13 alcmo_dich2 0.130103892 0.096870262 1.34307361 4040.249
## 14 marmo_dich2 -0.006249400 0.095462174 -0.06546467 4202.611
## 15 ecigmo_dich2 0.256384099 0.102743644 2.49537674 4214.765
## 16 smokmo_dich2 0.048971772 0.113110964 0.43295337 4231.389
## 17 cigarmo_dich2 0.341582120 0.125922778 2.71263171 4204.409
## 18 rxstmo_dich2 0.385136493 0.217146744 1.77362316 4165.283
## 19 rxpkmo_dich2 0.443037537 0.269833991 1.64188928 4232.371
## 20 rxsedmo_dich2 -0.476287419 0.321305439 -1.48235094 4230.021
## 21 mhdays 0.006811424 0.004274753 1.59340742 4126.659
## 22 (Intercept) -2.866103961 0.171602238 -16.70201971 4127.506
## 23 age21+ 0.082238018 0.123791441 0.66432717 4224.301
## 24 genderMale 1.324559262 0.120845754 10.96074312 4219.110
## 25 genderOther -0.813988469 0.430625302 -1.89024766 4230.822
## 26 raceAfrican American/Black -0.209299486 0.256163247 -0.81705509 4169.629
## 27 raceAsian -1.103763760 0.257448472 -4.28731913 4201.251
## 28 raceMultiracial -0.622804678 0.332755645 -1.87165774 4224.998
## 29 raceOther -0.445765742 0.369938674 -1.20497200 4087.471
## 30 ethnicityHispanic/Latino -0.355214940 0.196531201 -1.80742263 2640.938
## 31 studentstatusPart Time 0.021964334 0.290288402 0.07566384 3855.467
## 32 residenceOn Campus 0.193992021 0.124965210 1.55236822 4228.368
## 33 greekmemGreek 0.317825043 0.157235560 2.02133056 4227.704
## 34 alcmo_dich2 0.209037118 0.164483874 1.27086694 3845.761
## 35 marmo_dich2 0.219968449 0.144541033 1.52184086 4222.280
## 36 ecigmo_dich2 0.440227608 0.154942912 2.84122456 4223.836
## 37 smokmo_dich2 0.183682084 0.159708129 1.15011105 4229.597
## 38 cigarmo_dich2 0.870838020 0.154654327 5.63086749 4224.778
## 39 rxstmo_dich2 0.620928646 0.265148111 2.34181810 4219.977
## 40 rxpkmo_dich2 0.331613291 0.363789158 0.91155353 4211.242
## 41 rxsedmo_dich2 -0.540195757 0.429776539 -1.25692240 4230.902
## 42 mhdays -0.024503872 0.007591798 -3.22767690 4183.457
## p.value
## 1 1.998041e-31
## 2 8.011605e-01
## 3 3.508447e-02
## 4 7.768694e-02
## 5 9.159464e-04
## 6 1.340445e-08
## 7 5.252904e-01
## 8 1.481913e-02
## 9 4.505992e-01
## 10 3.538168e-01
## 11 1.404928e-01
## 12 2.257867e-01
## 13 1.793236e-01
## 14 9.478071e-01
## 15 1.262030e-02
## 16 6.650707e-01
## 17 6.702304e-03
## 18 7.619855e-02
## 19 1.006872e-01
## 20 1.383214e-01
## 21 1.111454e-01
## 22 1.191630e-60
## 23 5.065172e-01
## 24 1.385268e-27
## 25 5.879313e-02
## 26 4.139436e-01
## 27 1.848998e-05
## 28 6.132301e-02
## 29 2.282839e-01
## 30 7.081015e-02
## 31 9.396905e-01
## 32 1.206490e-01
## 33 4.330845e-02
## 34 2.038529e-01
## 35 1.281238e-01
## 36 4.515593e-03
## 37 2.501632e-01
## 38 1.909430e-08
## 39 1.923633e-02
## 40 3.620560e-01
## 41 2.088511e-01
## 42 1.257627e-03
## View fraction of missing information (FMI)
data.frame(
term = summary(pooled1d)$term,
fmi = pooled1d$pooled$fmi,
lambda = pooled1d$pooled$lambda
)
## term fmi lambda
## 1 (Intercept) 0.0038050563 0.0033318625
## 2 age21+ 0.0012422950 0.0007702010
## 3 genderMale 0.0050766429 0.0046023962
## 4 genderOther 0.0059030730 0.0054279628
## 5 raceAfrican American/Black 0.0037247745 0.0032516359
## 6 raceAsian 0.0096906308 0.0092097686
## 7 raceMultiracial 0.0069751777 0.0064987375
## 8 raceOther 0.0240862790 0.0235573035
## 9 ethnicityHispanic/Latino 0.0779234002 0.0768785026
## 10 studentstatusPart Time 0.0013909272 0.0009188070
## 11 residenceOn Campus 0.0012203145 0.0007482240
## 12 greekmemGreek 0.0024203853 0.0019479559
## 13 alcmo_dich2 0.0130607641 0.0125723312
## 14 marmo_dich2 0.0044333124 0.0039596402
## 15 ecigmo_dich2 0.0032768713 0.0028040161
## 16 smokmo_dich2 0.0010217474 0.0005496839
## 17 cigarmo_dich2 0.0042758206 0.0038022761
## 18 rxstmo_dich2 0.0071137155 0.0066370862
## 19 rxpkmo_dich2 0.0008288107 0.0003567655
## 20 rxsedmo_dich2 0.0012704648 0.0007983662
## 21 mhdays 0.0092551049 0.0087750531
## 22 (Intercept) 0.0092122352 0.0087322612
## 23 age21+ 0.0021435298 0.0016712053
## 24 genderMale 0.0027949284 0.0023223317
## 25 genderOther 0.0011275084 0.0006554316
## 26 raceAfrican American/Black 0.0068424704 0.0063662076
## 27 raceAsian 0.0045499649 0.0040761949
## 28 raceMultiracial 0.0020477496 0.0015754577
## 29 raceOther 0.0111005435 0.0106167933
## 30 ethnicityHispanic/Latino 0.0504394294 0.0497205931
## 31 studentstatusPart Time 0.0193361487 0.0188275671
## 32 residenceOn Campus 0.0015461473 0.0010739946
## 33 greekmemGreek 0.0016507126 0.0011785352
## 34 alcmo_dich2 0.0196291573 0.0191194450
## 35 marmo_dich2 0.0024089518 0.0019365271
## 36 ecigmo_dich2 0.0022061273 0.0017337805
## 37 smokmo_dich2 0.0013435393 0.0008714279
## 38 cigarmo_dich2 0.0020782390 0.0016059370
## 39 rxstmo_dich2 0.0026925757 0.0022200276
## 40 rxpkmo_dich2 0.0036369353 0.0031638556
## 41 rxsedmo_dich2 0.0011127069 0.0006406320
## 42 mhdays 0.0059202287 0.0054450990
# Fit model 2a: null model
fit2a <- with(imp2, multinom(lca2_class ~ 1, trace = FALSE))
# Pool model 2a
pooled2a <- pool(fit2a)
summary(pooled2a)
## term estimate std.error statistic df p.value
## 1 (Intercept) -1.154309 0.06086872 -18.96391 1812.003 2.663345e-73
## 2 (Intercept) -1.212296 0.06223483 -19.47939 1812.003 6.943915e-77
## 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) -1.69799209 0.1338409 -12.68664947 1781.002
## 2 age21+ -0.46173310 0.1365747 -3.38081054 1788.805
## 3 genderMale 0.94266193 0.1341224 7.02836930 1774.229
## 4 genderOther 1.58449153 0.2567450 6.17146016 1770.721
## 5 raceAfrican American/Black 0.80943093 0.2933297 2.75945800 1776.636
## 6 raceAsian 1.31784346 0.2279526 5.78121777 1776.389
## 7 raceMultiracial 0.36363376 0.2897738 1.25488812 1771.096
## 8 raceOther 0.85728641 0.3779526 2.26823813 1732.698
## 9 ethnicityHispanic/Latino 0.30412511 0.2022799 1.50348686 1411.081
## 10 studentstatusPart Time 0.21175644 0.3080844 0.68733253 1784.660
## 11 residenceOn Campus 0.14478930 0.1387769 1.04332402 1785.544
## 12 greekmemGreek -0.24134059 0.2098181 -1.15023719 1778.803
## 13 (Intercept) -2.46703055 0.1606329 -15.35818939 1787.159
## 14 age21+ 0.21040409 0.1411508 1.49063354 1788.957
## 15 genderMale 1.85634668 0.1428555 12.99457319 1789.042
## 16 genderOther 0.69507787 0.4070274 1.70769298 1789.892
## 17 raceAfrican American/Black 0.69796300 0.3145385 2.21900668 1764.266
## 18 raceAsian 0.25977498 0.2964208 0.87637242 1766.773
## 19 raceMultiracial -0.57252245 0.4058470 -1.41068537 1786.142
## 20 raceOther 0.09054606 0.5096306 0.17766999 1746.384
## 21 ethnicityHispanic/Latino -0.01640983 0.2301250 -0.07130831 1528.832
## 22 studentstatusPart Time -0.25664427 0.3645885 -0.70392852 1754.221
## 23 residenceOn Campus 0.16850838 0.1494115 1.12781400 1788.237
## 24 greekmemGreek 0.50451088 0.1788882 2.82025877 1789.545
## p.value
## 1 2.263035e-35
## 2 7.382719e-04
## 3 2.968837e-12
## 4 8.367722e-10
## 5 5.849146e-03
## 6 8.741975e-09
## 7 2.096849e-01
## 8 2.343730e-02
## 9 1.329373e-01
## 10 4.919625e-01
## 11 2.969395e-01
## 12 2.502008e-01
## 13 4.270937e-50
## 14 1.362340e-01
## 15 5.842851e-37
## 16 8.786674e-02
## 17 2.661307e-02
## 18 3.809467e-01
## 19 1.585115e-01
## 20 8.590028e-01
## 21 9.431617e-01
## 22 4.815707e-01
## 23 2.595498e-01
## 24 4.851391e-03
## 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.004846517 3.729623e-03
## 2 age21+ 0.001747511 6.320238e-04
## 3 genderMale 0.006865403 5.746523e-03
## 4 genderOther 0.007782829 6.662769e-03
## 5 raceAfrican American/Black 0.006191489 5.073365e-03
## 6 raceAsian 0.006262428 5.144228e-03
## 7 raceMultiracial 0.007687946 6.568016e-03
## 8 raceOther 0.015247144 1.411113e-02
## 9 ethnicityHispanic/Latino 0.049449872 4.810356e-02
## 10 studentstatusPart Time 0.003547305 2.431244e-03
## 11 residenceOn Campus 0.003200031 2.084134e-03
## 12 greekmemGreek 0.005546254 4.428767e-03
## 13 (Intercept) 0.002519369 1.403719e-03
## 14 age21+ 0.001671038 5.555603e-04
## 15 genderMale 0.001628036 5.125627e-04
## 16 genderOther 0.001177035 6.158767e-05
## 17 raceAfrican American/Black 0.009316755 8.194337e-03
## 18 raceAsian 0.008741385 7.619907e-03
## 19 raceMultiracial 0.002955595 1.839797e-03
## 20 raceOther 0.012904652 1.177485e-02
## 21 ethnicityHispanic/Latino 0.038681319 3.742455e-02
## 22 studentstatusPart Time 0.011425022 1.029858e-02
## 23 residenceOn Campus 0.002024407 9.088757e-04
## 24 greekmemGreek 0.001365416 2.499634e-04
# 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) -1.54871685 0.1650662 -9.38239969 1757.975
## 2 age21+ -0.34016543 0.1465450 -2.32123567 1775.877
## 3 genderMale 0.98850353 0.1393690 7.09270675 1759.316
## 4 genderOther 1.59758385 0.2599879 6.14483889 1757.067
## 5 raceAfrican American/Black 0.72168113 0.2952697 2.44414220 1766.819
## 6 raceAsian 1.19552970 0.2326831 5.13801622 1767.458
## 7 raceMultiracial 0.35158468 0.2913588 1.20670682 1758.537
## 8 raceOther 0.89757155 0.3822535 2.34810554 1719.494
## 9 ethnicityHispanic/Latino 0.26266809 0.2040925 1.28700517 1419.661
## 10 studentstatusPart Time 0.23471064 0.3076426 0.76293274 1775.316
## 11 residenceOn Campus 0.14844440 0.1394382 1.06458893 1775.932
## 12 greekmemGreek -0.10854398 0.2134621 -0.50849298 1770.579
## 13 alcmo_dich2 -0.13700386 0.1709849 -0.80126297 1661.140
## 14 marmo_dich2 -0.06751389 0.1699816 -0.39718353 1760.296
## 15 ecigmo_dich2 0.05000586 0.1854190 0.26969114 1739.011
## 16 smokmo_dich2 -0.27077048 0.2077158 -1.30356247 1770.531
## 17 cigarmo_dich2 -0.36770739 0.2229386 -1.64936657 1768.844
## 18 (Intercept) -2.55258907 0.2070104 -12.33072621 1759.523
## 19 age21+ 0.08099796 0.1527501 0.53026435 1776.061
## 20 genderMale 1.73079805 0.1512836 11.44075110 1777.825
## 21 genderOther 0.66217056 0.4110492 1.61092776 1779.595
## 22 raceAfrican American/Black 0.82380269 0.3208733 2.56737673 1754.715
## 23 raceAsian 0.49738686 0.3046162 1.63283124 1749.531
## 24 raceMultiracial -0.55038864 0.4091703 -1.34513346 1776.285
## 25 raceOther 0.01482968 0.5140167 0.02885058 1732.895
## 26 ethnicityHispanic/Latino 0.01955432 0.2334742 0.08375365 1493.185
## 27 studentstatusPart Time -0.23941634 0.3701371 -0.64683148 1742.753
## 28 residenceOn Campus 0.09580748 0.1531652 0.62551722 1777.484
## 29 greekmemGreek 0.31336783 0.1882648 1.66450537 1779.619
## 30 alcmo_dich2 -0.18452058 0.2083786 -0.88550642 1737.454
## 31 marmo_dich2 0.07649457 0.1763341 0.43380484 1777.476
## 32 ecigmo_dich2 0.42129325 0.1911498 2.20399564 1775.470
## 33 smokmo_dich2 0.04784519 0.1922617 0.24885457 1778.142
## 34 cigarmo_dich2 0.59000166 0.1829290 3.22530404 1770.321
## p.value
## 1 1.919670e-20
## 2 2.038683e-02
## 3 1.898047e-12
## 4 9.880002e-10
## 5 1.461689e-02
## 6 3.083655e-07
## 7 2.277074e-01
## 8 1.898194e-02
## 9 1.983023e-01
## 10 4.456049e-01
## 11 2.872067e-01
## 12 6.111710e-01
## 13 4.230940e-01
## 14 6.912803e-01
## 15 7.874299e-01
## 16 1.925523e-01
## 17 9.925012e-02
## 18 1.441702e-33
## 19 5.959949e-01
## 20 2.695510e-29
## 21 1.073728e-01
## 22 1.032932e-02
## 23 1.026844e-01
## 24 1.787538e-01
## 25 9.769871e-01
## 26 9.332635e-01
## 27 5.178262e-01
## 28 5.317119e-01
## 29 9.618752e-02
## 30 3.760061e-01
## 31 6.644828e-01
## 32 2.765257e-02
## 33 8.035020e-01
## 34 1.281341e-03
## 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.008504563 0.0073772071
## 2 age21+ 0.003082039 0.0019599380
## 3 genderMale 0.008177789 0.0070509206
## 4 genderOther 0.008721360 0.0075936683
## 5 raceAfrican American/Black 0.006173094 0.0050487399
## 6 raceAsian 0.005984834 0.0048606743
## 7 raceMultiracial 0.008368629 0.0072414783
## 8 raceOther 0.015851088 0.0147070569
## 9 ethnicityHispanic/Latino 0.048092791 0.0467527008
## 10 studentstatusPart Time 0.003309446 0.0021872467
## 11 residenceOn Campus 0.003059284 0.0019371914
## 12 greekmemGreek 0.005013234 0.0038899573
## 13 alcmo_dich2 0.023884624 0.0227100956
## 14 marmo_dich2 0.007933941 0.0068074224
## 15 ecigmo_dich2 0.012492429 0.0113573703
## 16 smokmo_dich2 0.005028995 0.0039057051
## 17 cigarmo_dich2 0.005565024 0.0044412694
## 18 (Intercept) 0.008126634 0.0069998405
## 19 age21+ 0.003005675 0.0018836042
## 20 genderMale 0.002230661 0.0011088317
## 21 genderOther 0.001346148 0.0002244403
## 22 raceAfrican American/Black 0.009267627 0.0081390469
## 23 raceAsian 0.010404322 0.0092736984
## 24 raceMultiracial 0.002911854 0.0017898189
## 25 raceOther 0.013604589 0.0124668087
## 26 ethnicityHispanic/Latino 0.041344146 0.0400609635
## 27 studentstatusPart Time 0.011778185 0.0106447418
## 28 residenceOn Campus 0.002387701 0.0012658331
## 29 greekmemGreek 0.001333471 0.0002117641
## 30 alcmo_dich2 0.012781739 0.0116459958
## 31 marmo_dich2 0.002391242 0.0012693728
## 32 ecigmo_dich2 0.003247818 0.0021256462
## 33 smokmo_dich2 0.002081521 0.0009597229
## 34 cigarmo_dich2 0.005097288 0.0039739423
# 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) -1.75593714 0.182250041 -9.63476948 1742.669
## 2 age21+ -0.33568545 0.147192773 -2.28058377 1767.897
## 3 genderMale 1.05689657 0.142343673 7.42496347 1748.328
## 4 genderOther 1.50156775 0.262331001 5.72394321 1747.470
## 5 raceAfrican American/Black 0.71558540 0.297941492 2.40176485 1758.934
## 6 raceAsian 1.18717531 0.235279252 5.04581383 1756.116
## 7 raceMultiracial 0.32500154 0.292090501 1.11267410 1752.126
## 8 raceOther 0.84074172 0.384656764 2.18569331 1705.424
## 9 ethnicityHispanic/Latino 0.26999056 0.205183150 1.31585150 1386.785
## 10 studentstatusPart Time 0.19037277 0.310451102 0.61321337 1766.130
## 11 residenceOn Campus 0.13246308 0.139840236 0.94724580 1767.804
## 12 greekmemGreek -0.06154009 0.214443480 -0.28697578 1761.419
## 13 alcmo_dich2 -0.12587195 0.171792414 -0.73269796 1650.931
## 14 marmo_dich2 -0.05293592 0.172205316 -0.30740005 1752.408
## 15 ecigmo_dich2 0.01528989 0.186867515 0.08182209 1727.873
## 16 smokmo_dich2 -0.29867336 0.209769695 -1.42381557 1763.864
## 17 cigarmo_dich2 -0.31281089 0.224297826 -1.39462292 1762.223
## 18 rxstmo_dich2 -0.62268922 0.441424590 -1.41063555 1770.233
## 19 rxpkmo_dich2 0.11825321 0.432842902 0.27320122 1767.549
## 20 rxsedmo_dich2 0.61759743 0.554868200 1.11305248 1771.413
## 21 mhdays 0.02143039 0.007409292 2.89236694 1735.252
## 22 (Intercept) -2.45078385 0.219514016 -11.16458936 1754.955
## 23 age21+ 0.07305186 0.153096674 0.47716165 1768.130
## 24 genderMale 1.69351260 0.153445135 11.03660018 1769.290
## 25 genderOther 0.72420275 0.413839645 1.74995982 1770.568
## 26 raceAfrican American/Black 0.82379571 0.320990779 2.56641551 1746.488
## 27 raceAsian 0.48852621 0.304715735 1.60321951 1739.983
## 28 raceMultiracial -0.53880889 0.408976700 -1.31745621 1768.024
## 29 raceOther 0.03914116 0.514640841 0.07605529 1723.788
## 30 ethnicityHispanic/Latino 0.01081878 0.234851442 0.04606650 1432.760
## 31 studentstatusPart Time -0.22299840 0.371621413 -0.60006875 1733.655
## 32 residenceOn Campus 0.09982964 0.153697773 0.64951912 1768.992
## 33 greekmemGreek 0.29411482 0.189814594 1.54948474 1771.555
## 34 alcmo_dich2 -0.18389105 0.208641058 -0.88137520 1734.047
## 35 marmo_dich2 0.08552738 0.178386432 0.47945002 1769.056
## 36 ecigmo_dich2 0.42643349 0.191780272 2.22355243 1767.077
## 37 smokmo_dich2 0.05307844 0.194200590 0.27331760 1769.979
## 38 cigarmo_dich2 0.56924698 0.184603804 3.08361455 1760.079
## 39 rxstmo_dich2 0.06777027 0.302891861 0.22374411 1765.355
## 40 rxpkmo_dich2 0.24573057 0.400045787 0.61425611 1650.338
## 41 rxsedmo_dich2 -0.27960707 0.523598204 -0.53401076 1760.696
## 42 mhdays -0.01209008 0.009056721 -1.33492875 1683.953
## p.value
## 1 1.930990e-21
## 2 2.269163e-02
## 3 1.753513e-13
## 4 1.222426e-08
## 5 1.641935e-02
## 6 4.983354e-07
## 7 2.660011e-01
## 8 2.897365e-02
## 9 1.884414e-01
## 10 5.398142e-01
## 11 3.436430e-01
## 12 7.741646e-01
## 13 4.638467e-01
## 14 7.585754e-01
## 15 9.347977e-01
## 16 1.546767e-01
## 17 1.633054e-01
## 18 1.585278e-01
## 19 7.847305e-01
## 20 2.658369e-01
## 21 3.871161e-03
## 22 5.206606e-28
## 23 6.333060e-01
## 24 1.953368e-27
## 25 8.029847e-02
## 26 1.035827e-02
## 27 1.090677e-01
## 28 1.878565e-01
## 29 9.393839e-01
## 30 9.632637e-01
## 31 5.485388e-01
## 32 5.160872e-01
## 33 1.214439e-01
## 34 3.782369e-01
## 35 6.316778e-01
## 36 2.630459e-02
## 37 7.846410e-01
## 38 2.076794e-03
## 39 8.229823e-01
## 40 5.391308e-01
## 41 5.934015e-01
## 42 1.820801e-01
## 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.010208203 0.0090729045
## 2 age21+ 0.003087708 0.0019605513
## 3 genderMale 0.008935554 0.0078024743
## 4 genderOther 0.009135683 0.0080022764
## 5 raceAfrican American/Black 0.006167135 0.0050377371
## 6 raceAsian 0.006966427 0.0058361272
## 7 raceMultiracial 0.008012720 0.0068810422
## 8 raceOther 0.016877644 0.0157253838
## 9 ethnicityHispanic/Latino 0.050594943 0.0492267123
## 10 studentstatusPart Time 0.003786426 0.0026589325
## 11 residenceOn Campus 0.003126095 0.0019989222
## 12 greekmemGreek 0.005409915 0.0042812505
## 13 alcmo_dich2 0.024266980 0.0230856556
## 14 marmo_dich2 0.007941586 0.0068100090
## 15 ecigmo_dich2 0.013130711 0.0119890786
## 16 smokmo_dich2 0.004603713 0.0034756987
## 17 cigarmo_dich2 0.005152147 0.0040237050
## 18 rxstmo_dich2 0.002047283 0.0009204376
## 19 rxpkmo_dich2 0.003230404 0.0021031864
## 20 rxsedmo_dich2 0.001450273 0.0003235041
## 21 mhdays 0.011733969 0.0105955788
## 22 (Intercept) 0.007280221 0.0061495316
## 23 age21+ 0.002990417 0.0018632989
## 24 genderMale 0.002486284 0.0013593343
## 25 genderOther 0.001883441 0.0007566244
## 26 raceAfrican American/Black 0.009361426 0.0082276401
## 27 raceAsian 0.010777073 0.0096406770
## 28 raceMultiracial 0.003034670 0.0019075342
## 29 raceOther 0.013862555 0.0127190669
## 30 ethnicityHispanic/Latino 0.046414190 0.0450840004
## 31 studentstatusPart Time 0.012045366 0.0109062872
## 32 residenceOn Campus 0.002619514 0.0014925249
## 33 greekmemGreek 0.001374273 0.0002475087
## 34 alcmo_dich2 0.011969330 0.0108304208
## 35 marmo_dich2 0.002590935 0.0014639544
## 36 ecigmo_dich2 0.003419913 0.0022926088
## 37 smokmo_dich2 0.002168532 0.0010416614
## 38 cigarmo_dich2 0.005825202 0.0046961503
## 39 rxstmo_dich2 0.004074805 0.0029471437
## 40 rxpkmo_dich2 0.024339244 0.0231575834
## 41 rxsedmo_dich2 0.005636080 0.0045072085
## 42 mhdays 0.019999499 0.0188362608
- 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.034133
vif_check1c <- lm(
as.numeric(lca1_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
data = df6_model1
)
vif(vif_check1c)
## GVIF Df GVIF^(1/(2*Df))
## age 1.281054 1 1.131837
## gender 1.117460 2 1.028154
## race 1.237150 4 1.026958
## ethnicity 1.150534 1 1.072629
## studentstatus 1.031830 1 1.015791
## residence 1.227470 1 1.107913
## greekmem 1.107798 1 1.052520
## alcmo_dich 1.434667 1 1.197776
## marmo_dich 1.657801 1 1.287556
## ecigmo_dich 1.897877 1 1.377635
## smokmo_dich 1.760500 1 1.326838
## cigarmo_dich 1.451277 1 1.204690
vif_check1d <- lm(
as.numeric(lca1_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
data = df6_model1
)
vif(vif_check1d)
## GVIF Df GVIF^(1/(2*Df))
## age 1.280975 1 1.131801
## gender 1.168048 2 1.039597
## race 1.242775 4 1.027541
## ethnicity 1.153518 1 1.074020
## studentstatus 1.035161 1 1.017429
## residence 1.228151 1 1.108220
## greekmem 1.116074 1 1.056444
## alcmo_dich 1.434321 1 1.197631
## marmo_dich 1.687643 1 1.299093
## ecigmo_dich 1.899183 1 1.378109
## smokmo_dich 1.788407 1 1.337314
## cigarmo_dich 1.476468 1 1.215100
## rxstmo_dich 1.366379 1 1.168922
## rxpkmo_dich 1.366153 1 1.168826
## rxsedmo_dich 1.374381 1 1.172340
## mhdays 1.078938 1 1.038720
# Fit a simple linear model with same predictors to extract VIF
vif_check2b <- lm(
as.numeric(lca2_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem,
data = df6_model2
)
vif(vif_check2b)
## GVIF Df GVIF^(1/(2*Df))
## age 1.113499 1 1.055225
## gender 1.019402 2 1.004816
## race 1.127217 4 1.015082
## ethnicity 1.112375 1 1.054692
## studentstatus 1.029068 1 1.014430
## residence 1.206640 1 1.098471
## greekmem 1.102474 1 1.049988
vif_check2c <- lm(
as.numeric(lca2_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich,
data = df6_model2
)
vif(vif_check2c)
## GVIF Df GVIF^(1/(2*Df))
## age 1.292036 1 1.136678
## gender 1.148279 2 1.035170
## race 1.190469 4 1.022033
## ethnicity 1.117230 1 1.056991
## studentstatus 1.030227 1 1.015001
## residence 1.226182 1 1.107331
## greekmem 1.148999 1 1.071914
## alcmo_dich 1.461152 1 1.208781
## marmo_dich 1.651463 1 1.285093
## ecigmo_dich 1.954846 1 1.398158
## smokmo_dich 1.831376 1 1.353284
## cigarmo_dich 1.558880 1 1.248551
vif_check2d <- lm(
as.numeric(lca2_class) ~ age + gender + race + ethnicity + studentstatus + residence + greekmem + alcmo_dich + marmo_dich + ecigmo_dich + smokmo_dich + cigarmo_dich + rxstmo_dich + rxpkmo_dich + rxsedmo_dich + mhdays,
data = df6_model2
)
vif(vif_check2d)
## GVIF Df GVIF^(1/(2*Df))
## age 1.288033 1 1.134916
## gender 1.209071 2 1.048608
## race 1.201590 4 1.023221
## ethnicity 1.124182 1 1.060275
## studentstatus 1.035204 1 1.017450
## residence 1.228619 1 1.108431
## greekmem 1.162621 1 1.078249
## alcmo_dich 1.456158 1 1.206714
## marmo_dich 1.682050 1 1.296939
## ecigmo_dich 1.954040 1 1.397870
## smokmo_dich 1.867710 1 1.366642
## cigarmo_dich 1.590546 1 1.261169
## rxstmo_dich 1.513676 1 1.230315
## rxpkmo_dich 1.528437 1 1.236300
## rxsedmo_dich 1.536385 1 1.239510
## mhdays 1.094832 1 1.046342
### Absence of complete separation
## Check for separation in model 1
table(df6_model1$age, df6_model1$lca1_class)
##
## 3 1 2
## <21 (ref) 1570 625 212
## 21+ 1117 521 233
table(df6_model1$gender, df6_model1$lca1_class)
##
## 3 1 2
## Female (ref) 1780 744 159
## Male 722 344 280
## Other 164 55 6
table(df6_model1$race, df6_model1$lca1_class)
##
## 3 1 2
## White (ref) 1964 962 383
## African American/Black 162 39 21
## Asian 303 57 19
## Multiracial 125 54 11
## Other 104 26 10
table(df6_model1$ethnicity, df6_model1$lca1_class)
##
## 3 1 2
## Not Hispanic/Latino (ref) 2182 968 375
## Hispanic/Latino 344 136 38
table(df6_model1$studentstatus, df6_model1$lca1_class)
##
## 3 1 2
## Full Time (ref) 2582 1092 425
## Part Time 98 53 17
table(df6_model1$residence, df6_model1$lca1_class)
##
## 3 1 2
## Off Campus (ref) 1316 524 205
## On Campus 1370 621 240
table(df6_model1$greekmem, df6_model1$lca1_class)
##
## 3 1 2
## No Greek (ref) 2396 992 353
## Greek 287 153 92
table(df6_model1$alcmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 894 283 75
## 2 1770 850 362
table(df6_model1$marmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 1816 677 208
## 2 865 467 237
table(df6_model1$ecigmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 1956 709 218
## 2 721 434 227
table(df6_model1$smokmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 2230 861 267
## 2 456 285 177
table(df6_model1$cigarmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 2453 961 271
## 2 227 180 172
table(df6_model1$rxstmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 2619 1083 399
## 2 64 58 45
table(df6_model1$rxpkmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 2639 1109 424
## 2 44 36 20
table(df6_model1$rxsedmo_dich, df6_model1$lca1_class)
##
## 3 1 2
## 1 2642 1123 432
## 2 39 23 13
## Check for separation in model 2
table(df6_model2$age, df6_model2$lca2_class)
##
## 3 1 2
## <21 (ref) 584 223 156
## 21+ 542 132 179
table(df6_model2$gender, df6_model2$lca2_class)
##
## 3 1 2
## Female (ref) 762 151 88
## Male 324 166 239
## Other 38 35 8
table(df6_model2$race, df6_model2$lca2_class)
##
## 3 1 2
## White (ref) 966 251 282
## African American/Black 38 20 19
## Asian 47 46 19
## Multiracial 48 20 8
## Other 21 15 6
table(df6_model2$ethnicity, df6_model2$lca2_class)
##
## 3 1 2
## Not Hispanic/Latino (ref) 969 280 277
## Hispanic/Latino 118 51 31
table(df6_model2$studentstatus, df6_model2$lca2_class)
##
## 3 1 2
## Full Time (ref) 1076 337 323
## Part Time 48 17 11
table(df6_model2$residence, df6_model2$lca2_class)
##
## 3 1 2
## Off Campus (ref) 519 164 158
## On Campus 606 191 177
table(df6_model2$greekmem, df6_model2$lca2_class)
##
## 3 1 2
## No Greek (ref) 969 318 259
## Greek 156 35 76
table(df6_model2$alcmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 248 120 59
## 2 867 231 269
table(df6_model2$marmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 642 236 155
## 2 482 119 180
table(df6_model2$ecigmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 681 250 157
## 2 443 103 178
table(df6_model2$smokmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 822 295 196
## 2 304 60 138
table(df6_model2$cigarmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 927 314 190
## 2 194 40 143
table(df6_model2$rxstmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 1060 344 298
## 2 62 10 36
table(df6_model2$rxpkmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 1092 342 315
## 2 34 12 18
table(df6_model2$rxsedmo_dich, df6_model2$lca2_class)
##
## 3 1 2
## 1 1103 347 326
## 2 23 8 9
- 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): 329.129 | df: 22 | p: 0.0000
## SD across imputations: 2.144
##
## --- Model 2: B → C ---
## Mean LRT (chi-sq): 60.197 | df: 10 | p: 0.0000
## SD across imputations: 0.626
##
## --- Model 2: A → C ---
## Mean LRT (chi-sq): 389.326 | df: 32 | p: 0.0000
## SD across imputations: 2.156
##
## --- Model 2: C → D ---
## Mean LRT (chi-sq): 16.607 | df: 8 | p: 0.0345
## SD across imputations: 0.735
##
## --- Model 2: B → D ---
## Mean LRT (chi-sq): 76.804 | df: 18 | p: 0.0000
## SD across imputations: 0.865
##
## --- Model 2: A → D ---
## Mean LRT (chi-sq): 405.933 | df: 40 | p: 0.0000
## SD across imputations: 2.092
print(lrt_table)
## comparison mean_LRT_chisq df p_value sd_across_imp
## 1 Model 1: A → B 379.035 22 0.0000 2.037
## 2 Model 1: B → C 155.214 10 0.0000 0.766
## 3 Model 1: C → D 28.257 8 0.0004 0.821
## 4 Model 1: B → D 183.470 18 0.0000 1.012
## 5 Model 1: A → C 534.248 32 0.0000 2.163
## 6 Model 1: A → D 562.505 40 0.0000 2.274
## 7 Model 2: A → B 329.129 22 0.0000 2.144
## 8 Model 2: B → C 60.197 10 0.0000 0.626
## 9 Model 2: A → C 389.326 32 0.0000 2.156
## 10 Model 2: C → D 16.607 8 0.0345 0.735
## 11 Model 2: B → D 76.804 18 0.0000 0.865
## 12 Model 2: A → D 405.933 40 0.0000 2.092
- 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.183 0.141
## Model 2b Class 1 vs. Class 3 age21+ 0.630 0.482
## Model 2b Class 1 vs. Class 3 genderMale 2.567 1.973
## Model 2b Class 1 vs. Class 3 genderOther 4.877 2.948
## Model 2b Class 1 vs. Class 3 raceAfrican American/Black 2.247 1.264
## Model 2b Class 1 vs. Class 3 raceAsian 3.735 2.389
## Model 2b Class 1 vs. Class 3 raceMultiracial 1.439 0.815
## Model 2b Class 1 vs. Class 3 raceOther 2.357 1.124
## Model 2b Class 1 vs. Class 3 ethnicityHispanic/Latino 1.355 0.912
## Model 2b Class 1 vs. Class 3 studentstatusPart Time 1.236 0.676
## Model 2b Class 1 vs. Class 3 residenceOn Campus 1.156 0.881
## Model 2b Class 1 vs. Class 3 greekmemGreek 0.786 0.521
## Model 2b Class 2 vs. Class 3 (Intercept) 0.085 0.062
## Model 2b Class 2 vs. Class 3 age21+ 1.234 0.936
## Model 2b Class 2 vs. Class 3 genderMale 6.400 4.837
## Model 2b Class 2 vs. Class 3 genderOther 2.004 0.902
## Model 2b Class 2 vs. Class 3 raceAfrican American/Black 2.010 1.085
## Model 2b Class 2 vs. Class 3 raceAsian 1.297 0.725
## Model 2b Class 2 vs. Class 3 raceMultiracial 0.564 0.255
## Model 2b Class 2 vs. Class 3 raceOther 1.095 0.403
## Model 2b Class 2 vs. Class 3 ethnicityHispanic/Latino 0.984 0.627
## Model 2b Class 2 vs. Class 3 studentstatusPart Time 0.774 0.379
## Model 2b Class 2 vs. Class 3 residenceOn Campus 1.184 0.883
## Model 2b Class 2 vs. Class 3 greekmemGreek 1.656 1.166
## CI_upper p_value sig
## 0.238 0.0000 ***
## 0.824 0.0007 ***
## 3.339 0.0000 ***
## 8.066 0.0000 ***
## 3.992 0.0058 **
## 5.839 0.0000 ***
## 2.539 0.2097
## 4.943 0.0234 *
## 2.015 0.1329
## 2.261 0.4920
## 1.517 0.2969
## 1.185 0.2502
## 0.116 0.0000 ***
## 1.628 0.1362
## 8.468 0.0000 ***
## 4.450 0.0879 .
## 3.723 0.0266 *
## 2.318 0.3809
## 1.250 0.1585
## 2.973 0.8590
## 1.544 0.9432
## 1.581 0.4816
## 1.586 0.2595
## 2.352 0.0049 **
print(or_2c, row.names = FALSE)
## model contrast term OR CI_lower
## Model 2c Class 1 vs. Class 3 (Intercept) 0.213 0.154
## Model 2c Class 1 vs. Class 3 age21+ 0.712 0.534
## Model 2c Class 1 vs. Class 3 genderMale 2.687 2.045
## Model 2c Class 1 vs. Class 3 genderOther 4.941 2.968
## Model 2c Class 1 vs. Class 3 raceAfrican American/Black 2.058 1.154
## Model 2c Class 1 vs. Class 3 raceAsian 3.305 2.095
## Model 2c Class 1 vs. Class 3 raceMultiracial 1.421 0.803
## Model 2c Class 1 vs. Class 3 raceOther 2.454 1.160
## Model 2c Class 1 vs. Class 3 ethnicityHispanic/Latino 1.300 0.872
## Model 2c Class 1 vs. Class 3 studentstatusPart Time 1.265 0.692
## Model 2c Class 1 vs. Class 3 residenceOn Campus 1.160 0.883
## Model 2c Class 1 vs. Class 3 greekmemGreek 0.897 0.590
## Model 2c Class 1 vs. Class 3 alcmo_dich2 0.872 0.624
## Model 2c Class 1 vs. Class 3 marmo_dich2 0.935 0.670
## Model 2c Class 1 vs. Class 3 ecigmo_dich2 1.051 0.731
## Model 2c Class 1 vs. Class 3 smokmo_dich2 0.763 0.508
## Model 2c Class 1 vs. Class 3 cigarmo_dich2 0.692 0.447
## Model 2c Class 2 vs. Class 3 (Intercept) 0.078 0.052
## Model 2c Class 2 vs. Class 3 age21+ 1.084 0.804
## Model 2c Class 2 vs. Class 3 genderMale 5.645 4.197
## Model 2c Class 2 vs. Class 3 genderOther 1.939 0.866
## Model 2c Class 2 vs. Class 3 raceAfrican American/Black 2.279 1.215
## Model 2c Class 2 vs. Class 3 raceAsian 1.644 0.905
## Model 2c Class 2 vs. Class 3 raceMultiracial 0.577 0.259
## Model 2c Class 2 vs. Class 3 raceOther 1.015 0.371
## Model 2c Class 2 vs. Class 3 ethnicityHispanic/Latino 1.020 0.645
## Model 2c Class 2 vs. Class 3 studentstatusPart Time 0.787 0.381
## Model 2c Class 2 vs. Class 3 residenceOn Campus 1.101 0.815
## Model 2c Class 2 vs. Class 3 greekmemGreek 1.368 0.946
## Model 2c Class 2 vs. Class 3 alcmo_dich2 0.832 0.553
## Model 2c Class 2 vs. Class 3 marmo_dich2 1.079 0.764
## Model 2c Class 2 vs. Class 3 ecigmo_dich2 1.524 1.048
## Model 2c Class 2 vs. Class 3 smokmo_dich2 1.049 0.720
## Model 2c Class 2 vs. Class 3 cigarmo_dich2 1.804 1.260
## CI_upper p_value sig
## 0.294 0.0000 ***
## 0.948 0.0204 *
## 3.531 0.0000 ***
## 8.225 0.0000 ***
## 3.671 0.0146 *
## 5.215 0.0000 ***
## 2.516 0.2277
## 5.190 0.0190 *
## 1.940 0.1983
## 2.311 0.4456
## 1.525 0.2872
## 1.363 0.6112
## 1.219 0.4231
## 1.304 0.6913
## 1.512 0.7874
## 1.146 0.1926
## 1.072 0.0993 .
## 0.117 0.0000 ***
## 1.463 0.5960
## 7.594 0.0000 ***
## 4.340 0.1074
## 4.275 0.0103 *
## 2.987 0.1027
## 1.286 0.1788
## 2.780 0.9770
## 1.611 0.9333
## 1.626 0.5178
## 1.486 0.5317
## 1.979 0.0962 .
## 1.251 0.3760
## 1.525 0.6645
## 2.217 0.0277 *
## 1.529 0.8035
## 2.582 0.0013 **
print(or_2d, row.names = FALSE)
## model contrast term OR CI_lower
## Model 2d Class 1 vs. Class 3 (Intercept) 0.173 0.121
## Model 2d Class 1 vs. Class 3 age21+ 0.715 0.536
## Model 2d Class 1 vs. Class 3 genderMale 2.877 2.177
## Model 2d Class 1 vs. Class 3 genderOther 4.489 2.684
## Model 2d Class 1 vs. Class 3 raceAfrican American/Black 2.045 1.141
## Model 2d Class 1 vs. Class 3 raceAsian 3.278 2.067
## Model 2d Class 1 vs. Class 3 raceMultiracial 1.384 0.781
## Model 2d Class 1 vs. Class 3 raceOther 2.318 1.091
## Model 2d Class 1 vs. Class 3 ethnicityHispanic/Latino 1.310 0.876
## Model 2d Class 1 vs. Class 3 studentstatusPart Time 1.210 0.658
## Model 2d Class 1 vs. Class 3 residenceOn Campus 1.142 0.868
## Model 2d Class 1 vs. Class 3 greekmemGreek 0.940 0.618
## Model 2d Class 1 vs. Class 3 alcmo_dich2 0.882 0.630
## Model 2d Class 1 vs. Class 3 marmo_dich2 0.948 0.677
## Model 2d Class 1 vs. Class 3 ecigmo_dich2 1.015 0.704
## Model 2d Class 1 vs. Class 3 smokmo_dich2 0.742 0.492
## Model 2d Class 1 vs. Class 3 cigarmo_dich2 0.731 0.471
## Model 2d Class 1 vs. Class 3 rxstmo_dich2 0.536 0.226
## Model 2d Class 1 vs. Class 3 rxpkmo_dich2 1.126 0.482
## Model 2d Class 1 vs. Class 3 rxsedmo_dich2 1.854 0.625
## Model 2d Class 1 vs. Class 3 mhdays 1.022 1.007
## Model 2d Class 2 vs. Class 3 (Intercept) 0.086 0.056
## Model 2d Class 2 vs. Class 3 age21+ 1.076 0.797
## Model 2d Class 2 vs. Class 3 genderMale 5.439 4.026
## Model 2d Class 2 vs. Class 3 genderOther 2.063 0.917
## Model 2d Class 2 vs. Class 3 raceAfrican American/Black 2.279 1.215
## Model 2d Class 2 vs. Class 3 raceAsian 1.630 0.897
## Model 2d Class 2 vs. Class 3 raceMultiracial 0.583 0.262
## Model 2d Class 2 vs. Class 3 raceOther 1.040 0.379
## Model 2d Class 2 vs. Class 3 ethnicityHispanic/Latino 1.011 0.638
## Model 2d Class 2 vs. Class 3 studentstatusPart Time 0.800 0.386
## Model 2d Class 2 vs. Class 3 residenceOn Campus 1.105 0.818
## Model 2d Class 2 vs. Class 3 greekmemGreek 1.342 0.925
## Model 2d Class 2 vs. Class 3 alcmo_dich2 0.832 0.553
## Model 2d Class 2 vs. Class 3 marmo_dich2 1.089 0.768
## Model 2d Class 2 vs. Class 3 ecigmo_dich2 1.532 1.052
## Model 2d Class 2 vs. Class 3 smokmo_dich2 1.055 0.721
## Model 2d Class 2 vs. Class 3 cigarmo_dich2 1.767 1.231
## Model 2d Class 2 vs. Class 3 rxstmo_dich2 1.070 0.591
## Model 2d Class 2 vs. Class 3 rxpkmo_dich2 1.279 0.584
## Model 2d Class 2 vs. Class 3 rxsedmo_dich2 0.756 0.271
## Model 2d Class 2 vs. Class 3 mhdays 0.988 0.971
## CI_upper p_value sig
## 0.247 0.0000 ***
## 0.954 0.0227 *
## 3.803 0.0000 ***
## 7.506 0.0000 ***
## 3.668 0.0164 *
## 5.198 0.0000 ***
## 2.453 0.2660
## 4.927 0.0290 *
## 1.958 0.1884
## 2.223 0.5398
## 1.502 0.3436
## 1.432 0.7742
## 1.235 0.4638
## 1.329 0.7586
## 1.465 0.9348
## 1.119 0.1547
## 1.135 0.1633
## 1.274 0.1585
## 2.629 0.7847
## 5.502 0.2658
## 1.037 0.0039 **
## 0.133 0.0000 ***
## 1.452 0.6333
## 7.347 0.0000 ***
## 4.643 0.0803 .
## 4.276 0.0104 *
## 2.962 0.1091
## 1.301 0.1879
## 2.851 0.9394
## 1.602 0.9633
## 1.658 0.5485
## 1.493 0.5161
## 1.947 0.1214
## 1.252 0.3782
## 1.545 0.6317
## 2.231 0.0263 *
## 1.543 0.7846
## 2.537 0.0021 **
## 1.938 0.8230
## 2.801 0.5391
## 2.110 0.5934
## 1.006 0.1821
- 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 -1683.88 -1520.04 0.097298 0.1651 0.1957
## Model 2c -1683.88 -1490.52 0.114834 0.1918 0.2274
## Model 2d -1683.88 -1481.95 0.119921 0.1994 0.2364
# 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.0973 0.1957
## 2c - 2b 0.0175 0.0317
## 2d - 2c 0.0051 0.0090
- 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.183 [0.141, 0.238]
## 102 2b Class 1 vs. Class 3 age21+ 0.630 [0.482, 0.824]
## 103 2b Class 1 vs. Class 3 genderMale 2.567 [1.973, 3.339]
## 104 2b Class 1 vs. Class 3 genderOther 4.877 [2.948, 8.066]
## 105 2b Class 1 vs. Class 3 raceAfrican American/Black 2.247 [1.264, 3.992]
## 106 2b Class 1 vs. Class 3 raceAsian 3.735 [2.389, 5.839]
## 107 2b Class 1 vs. Class 3 raceMultiracial 1.439 [0.815, 2.539]
## 108 2b Class 1 vs. Class 3 raceOther 2.357 [1.124, 4.943]
## 109 2b Class 1 vs. Class 3 ethnicityHispanic/Latino 1.355 [0.912, 2.015]
## 110 2b Class 1 vs. Class 3 studentstatusPart Time 1.236 [0.676, 2.261]
## 111 2b Class 1 vs. Class 3 residenceOn Campus 1.156 [0.881, 1.517]
## 112 2b Class 1 vs. Class 3 greekmemGreek 0.786 [0.521, 1.185]
## 113 2b Class 2 vs. Class 3 (Intercept) 0.085 [0.062, 0.116]
## 114 2b Class 2 vs. Class 3 age21+ 1.234 [0.936, 1.628]
## 115 2b Class 2 vs. Class 3 genderMale 6.400 [4.837, 8.468]
## 116 2b Class 2 vs. Class 3 genderOther 2.004 [0.902, 4.450]
## 117 2b Class 2 vs. Class 3 raceAfrican American/Black 2.010 [1.085, 3.723]
## 118 2b Class 2 vs. Class 3 raceAsian 1.297 [0.725, 2.318]
## 119 2b Class 2 vs. Class 3 raceMultiracial 0.564 [0.255, 1.250]
## 120 2b Class 2 vs. Class 3 raceOther 1.095 [0.403, 2.973]
## 121 2b Class 2 vs. Class 3 ethnicityHispanic/Latino 0.984 [0.627, 1.544]
## 122 2b Class 2 vs. Class 3 studentstatusPart Time 0.774 [0.379, 1.581]
## 123 2b Class 2 vs. Class 3 residenceOn Campus 1.184 [0.883, 1.586]
## 124 2b Class 2 vs. Class 3 greekmemGreek 1.656 [1.166, 2.352]
## 125 2c Class 1 vs. Class 3 (Intercept) 0.213 [0.154, 0.294]
## 126 2c Class 1 vs. Class 3 age21+ 0.712 [0.534, 0.948]
## 127 2c Class 1 vs. Class 3 genderMale 2.687 [2.045, 3.531]
## 128 2c Class 1 vs. Class 3 genderOther 4.941 [2.968, 8.225]
## 129 2c Class 1 vs. Class 3 raceAfrican American/Black 2.058 [1.154, 3.671]
## 130 2c Class 1 vs. Class 3 raceAsian 3.305 [2.095, 5.215]
## 131 2c Class 1 vs. Class 3 raceMultiracial 1.421 [0.803, 2.516]
## 132 2c Class 1 vs. Class 3 raceOther 2.454 [1.160, 5.190]
## 133 2c Class 1 vs. Class 3 ethnicityHispanic/Latino 1.300 [0.872, 1.940]
## 134 2c Class 1 vs. Class 3 studentstatusPart Time 1.265 [0.692, 2.311]
## 135 2c Class 1 vs. Class 3 residenceOn Campus 1.160 [0.883, 1.525]
## 136 2c Class 1 vs. Class 3 greekmemGreek 0.897 [0.590, 1.363]
## 137 2c Class 1 vs. Class 3 alcmo_dich2 0.872 [0.624, 1.219]
## 138 2c Class 1 vs. Class 3 marmo_dich2 0.935 [0.670, 1.304]
## 139 2c Class 1 vs. Class 3 ecigmo_dich2 1.051 [0.731, 1.512]
## 140 2c Class 1 vs. Class 3 smokmo_dich2 0.763 [0.508, 1.146]
## 141 2c Class 1 vs. Class 3 cigarmo_dich2 0.692 [0.447, 1.072]
## 142 2c Class 2 vs. Class 3 (Intercept) 0.078 [0.052, 0.117]
## 143 2c Class 2 vs. Class 3 age21+ 1.084 [0.804, 1.463]
## 144 2c Class 2 vs. Class 3 genderMale 5.645 [4.197, 7.594]
## 145 2c Class 2 vs. Class 3 genderOther 1.939 [0.866, 4.340]
## 146 2c Class 2 vs. Class 3 raceAfrican American/Black 2.279 [1.215, 4.275]
## 147 2c Class 2 vs. Class 3 raceAsian 1.644 [0.905, 2.987]
## 148 2c Class 2 vs. Class 3 raceMultiracial 0.577 [0.259, 1.286]
## 149 2c Class 2 vs. Class 3 raceOther 1.015 [0.371, 2.780]
## 150 2c Class 2 vs. Class 3 ethnicityHispanic/Latino 1.020 [0.645, 1.611]
## 151 2c Class 2 vs. Class 3 studentstatusPart Time 0.787 [0.381, 1.626]
## 152 2c Class 2 vs. Class 3 residenceOn Campus 1.101 [0.815, 1.486]
## 153 2c Class 2 vs. Class 3 greekmemGreek 1.368 [0.946, 1.979]
## 154 2c Class 2 vs. Class 3 alcmo_dich2 0.832 [0.553, 1.251]
## 155 2c Class 2 vs. Class 3 marmo_dich2 1.079 [0.764, 1.525]
## 156 2c Class 2 vs. Class 3 ecigmo_dich2 1.524 [1.048, 2.217]
## 157 2c Class 2 vs. Class 3 smokmo_dich2 1.049 [0.720, 1.529]
## 158 2c Class 2 vs. Class 3 cigarmo_dich2 1.804 [1.260, 2.582]
## 159 2d Class 1 vs. Class 3 (Intercept) 0.173 [0.121, 0.247]
## 160 2d Class 1 vs. Class 3 age21+ 0.715 [0.536, 0.954]
## 161 2d Class 1 vs. Class 3 genderMale 2.877 [2.177, 3.803]
## 162 2d Class 1 vs. Class 3 genderOther 4.489 [2.684, 7.506]
## 163 2d Class 1 vs. Class 3 raceAfrican American/Black 2.045 [1.141, 3.668]
## 164 2d Class 1 vs. Class 3 raceAsian 3.278 [2.067, 5.198]
## 165 2d Class 1 vs. Class 3 raceMultiracial 1.384 [0.781, 2.453]
## 166 2d Class 1 vs. Class 3 raceOther 2.318 [1.091, 4.927]
## 167 2d Class 1 vs. Class 3 ethnicityHispanic/Latino 1.310 [0.876, 1.958]
## 168 2d Class 1 vs. Class 3 studentstatusPart Time 1.210 [0.658, 2.223]
## 169 2d Class 1 vs. Class 3 residenceOn Campus 1.142 [0.868, 1.502]
## 170 2d Class 1 vs. Class 3 greekmemGreek 0.940 [0.618, 1.432]
## 171 2d Class 1 vs. Class 3 alcmo_dich2 0.882 [0.630, 1.235]
## 172 2d Class 1 vs. Class 3 marmo_dich2 0.948 [0.677, 1.329]
## 173 2d Class 1 vs. Class 3 ecigmo_dich2 1.015 [0.704, 1.465]
## 174 2d Class 1 vs. Class 3 smokmo_dich2 0.742 [0.492, 1.119]
## 175 2d Class 1 vs. Class 3 cigarmo_dich2 0.731 [0.471, 1.135]
## 176 2d Class 1 vs. Class 3 rxstmo_dich2 0.536 [0.226, 1.274]
## 177 2d Class 1 vs. Class 3 rxpkmo_dich2 1.126 [0.482, 2.629]
## 178 2d Class 1 vs. Class 3 rxsedmo_dich2 1.854 [0.625, 5.502]
## 179 2d Class 1 vs. Class 3 mhdays 1.022 [1.007, 1.037]
## 180 2d Class 2 vs. Class 3 (Intercept) 0.086 [0.056, 0.133]
## 181 2d Class 2 vs. Class 3 age21+ 1.076 [0.797, 1.452]
## 182 2d Class 2 vs. Class 3 genderMale 5.439 [4.026, 7.347]
## 183 2d Class 2 vs. Class 3 genderOther 2.063 [0.917, 4.643]
## 184 2d Class 2 vs. Class 3 raceAfrican American/Black 2.279 [1.215, 4.276]
## 185 2d Class 2 vs. Class 3 raceAsian 1.630 [0.897, 2.962]
## 186 2d Class 2 vs. Class 3 raceMultiracial 0.583 [0.262, 1.301]
## 187 2d Class 2 vs. Class 3 raceOther 1.040 [0.379, 2.851]
## 188 2d Class 2 vs. Class 3 ethnicityHispanic/Latino 1.011 [0.638, 1.602]
## 189 2d Class 2 vs. Class 3 studentstatusPart Time 0.800 [0.386, 1.658]
## 190 2d Class 2 vs. Class 3 residenceOn Campus 1.105 [0.818, 1.493]
## 191 2d Class 2 vs. Class 3 greekmemGreek 1.342 [0.925, 1.947]
## 192 2d Class 2 vs. Class 3 alcmo_dich2 0.832 [0.553, 1.252]
## 193 2d Class 2 vs. Class 3 marmo_dich2 1.089 [0.768, 1.545]
## 194 2d Class 2 vs. Class 3 ecigmo_dich2 1.532 [1.052, 2.231]
## 195 2d Class 2 vs. Class 3 smokmo_dich2 1.055 [0.721, 1.543]
## 196 2d Class 2 vs. Class 3 cigarmo_dich2 1.767 [1.231, 2.537]
## 197 2d Class 2 vs. Class 3 rxstmo_dich2 1.070 [0.591, 1.938]
## 198 2d Class 2 vs. Class 3 rxpkmo_dich2 1.279 [0.584, 2.801]
## 199 2d Class 2 vs. Class 3 rxsedmo_dich2 0.756 [0.271, 2.110]
## 200 2d Class 2 vs. Class 3 mhdays 0.988 [0.971, 1.006]
## 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.005 0.097 0.196
## 102 0.0007 *** 0.002 0.097 0.196
## 103 0.0000 *** 0.007 0.097 0.196
## 104 0.0000 *** 0.008 0.097 0.196
## 105 0.0058 ** 0.006 0.097 0.196
## 106 0.0000 *** 0.006 0.097 0.196
## 107 0.2097 0.008 0.097 0.196
## 108 0.0234 * 0.015 0.097 0.196
## 109 0.1329 0.049 0.097 0.196
## 110 0.4920 0.004 0.097 0.196
## 111 0.2969 0.003 0.097 0.196
## 112 0.2502 0.006 0.097 0.196
## 113 0.0000 *** 0.003 0.097 0.196
## 114 0.1362 0.002 0.097 0.196
## 115 0.0000 *** 0.002 0.097 0.196
## 116 0.0879 . 0.001 0.097 0.196
## 117 0.0266 * 0.009 0.097 0.196
## 118 0.3809 0.009 0.097 0.196
## 119 0.1585 0.003 0.097 0.196
## 120 0.8590 0.013 0.097 0.196
## 121 0.9432 0.039 0.097 0.196
## 122 0.4816 0.011 0.097 0.196
## 123 0.2595 0.002 0.097 0.196
## 124 0.0049 ** 0.001 0.097 0.196
## 125 0.0000 *** 0.009 0.115 0.227
## 126 0.0204 * 0.003 0.115 0.227
## 127 0.0000 *** 0.008 0.115 0.227
## 128 0.0000 *** 0.009 0.115 0.227
## 129 0.0146 * 0.006 0.115 0.227
## 130 0.0000 *** 0.006 0.115 0.227
## 131 0.2277 0.008 0.115 0.227
## 132 0.0190 * 0.016 0.115 0.227
## 133 0.1983 0.048 0.115 0.227
## 134 0.4456 0.003 0.115 0.227
## 135 0.2872 0.003 0.115 0.227
## 136 0.6112 0.005 0.115 0.227
## 137 0.4231 0.024 0.115 0.227
## 138 0.6913 0.008 0.115 0.227
## 139 0.7874 0.012 0.115 0.227
## 140 0.1926 0.005 0.115 0.227
## 141 0.0993 . 0.006 0.115 0.227
## 142 0.0000 *** 0.008 0.115 0.227
## 143 0.5960 0.003 0.115 0.227
## 144 0.0000 *** 0.002 0.115 0.227
## 145 0.1074 0.001 0.115 0.227
## 146 0.0103 * 0.009 0.115 0.227
## 147 0.1027 0.010 0.115 0.227
## 148 0.1788 0.003 0.115 0.227
## 149 0.9770 0.014 0.115 0.227
## 150 0.9333 0.041 0.115 0.227
## 151 0.5178 0.012 0.115 0.227
## 152 0.5317 0.002 0.115 0.227
## 153 0.0962 . 0.001 0.115 0.227
## 154 0.3760 0.013 0.115 0.227
## 155 0.6645 0.002 0.115 0.227
## 156 0.0277 * 0.003 0.115 0.227
## 157 0.8035 0.002 0.115 0.227
## 158 0.0013 ** 0.005 0.115 0.227
## 159 0.0000 *** 0.010 0.120 0.236
## 160 0.0227 * 0.003 0.120 0.236
## 161 0.0000 *** 0.009 0.120 0.236
## 162 0.0000 *** 0.009 0.120 0.236
## 163 0.0164 * 0.006 0.120 0.236
## 164 0.0000 *** 0.007 0.120 0.236
## 165 0.2660 0.008 0.120 0.236
## 166 0.0290 * 0.017 0.120 0.236
## 167 0.1884 0.051 0.120 0.236
## 168 0.5398 0.004 0.120 0.236
## 169 0.3436 0.003 0.120 0.236
## 170 0.7742 0.005 0.120 0.236
## 171 0.4638 0.024 0.120 0.236
## 172 0.7586 0.008 0.120 0.236
## 173 0.9348 0.013 0.120 0.236
## 174 0.1547 0.005 0.120 0.236
## 175 0.1633 0.005 0.120 0.236
## 176 0.1585 0.002 0.120 0.236
## 177 0.7847 0.003 0.120 0.236
## 178 0.2658 0.001 0.120 0.236
## 179 0.0039 ** 0.012 0.120 0.236
## 180 0.0000 *** 0.007 0.120 0.236
## 181 0.6333 0.003 0.120 0.236
## 182 0.0000 *** 0.002 0.120 0.236
## 183 0.0803 . 0.002 0.120 0.236
## 184 0.0104 * 0.009 0.120 0.236
## 185 0.1091 0.011 0.120 0.236
## 186 0.1879 0.003 0.120 0.236
## 187 0.9394 0.014 0.120 0.236
## 188 0.9633 0.046 0.120 0.236
## 189 0.5485 0.012 0.120 0.236
## 190 0.5161 0.003 0.120 0.236
## 191 0.1214 0.001 0.120 0.236
## 192 0.3782 0.012 0.120 0.236
## 193 0.6317 0.003 0.120 0.236
## 194 0.0263 * 0.003 0.120 0.236
## 195 0.7846 0.002 0.120 0.236
## 196 0.0021 ** 0.006 0.120 0.236
## 197 0.8230 0.004 0.120 0.236
## 198 0.5391 0.024 0.120 0.236
## 199 0.5934 0.006 0.120 0.236
## 200 0.1821 0.020 0.120 0.236
library(dplyr)
library(ggplot2)
all_or <- bind_rows(or_1b, or_1c, or_1d, or_2b, or_2c, or_2d) |>
filter(term != "(Intercept)", p_value < 0.05) |>
mutate(
term_label = term |>
gsub("^gender", "Gender: ", x = _) |>
gsub("^race", "Race: ", x = _) |>
gsub("^age", "Age: ", x = _) |>
gsub("^greekmem", "Greek Mem: ", x = _) |>
gsub("ecigmo_dich2", "E-cig Use", x = _) |>
gsub("cigarmo_dich2", "Cigar Use", x = _) |>
gsub("rxstmo_dich2", "Rx Stimulant Use", x = _) |>
gsub("mhdays", "Mental Health Days", x = _),
series = if_else(grepl("1", model), "Model Series 1", "Model Series 2")
)
make_forest_plot <- function(data, title) {
data |>
mutate(
panel = paste0(model, " | ", contrast),
term_label = forcats::fct_reorder(term_label, OR)
) |>
ggplot(aes(x = OR, y = term_label)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "gray50") +
geom_pointrange(aes(xmin = CI_lower, xmax = CI_upper)) +
scale_x_log10() +
facet_wrap(~ panel, scales = "free_y", ncol = 2) +
labs(
title = title,
x = "Odds Ratio (log scale)",
y = NULL
) +
theme_bw(base_size = 11) +
theme(
strip.text = element_text(size = 8.5),
axis.text.y = element_text(size = 9),
panel.spacing = unit(0.8, "cm")
)
}
make_forest_plot(
filter(all_or, series == "Model Series 1"),
"Forest Plot: Model Series 1 — Significant Odds Ratios"
)
make_forest_plot(
filter(all_or, series == "Model Series 2"),
"Forest Plot: Model Series 2 — Significant Odds Ratios"
)
