- Install packages and load dataset.
## 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)
## 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 45
## 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
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
- 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)
# 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)
## 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)
## 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)
# 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"
)
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.90676011 0.06891914 -13.1568683 4227.598
## 2 age21+ 0.13332643 0.07610440 1.7518885 4249.725
## 3 genderMale 0.18877788 0.08007697 2.3574554 4223.170
## 4 genderOther -0.24874855 0.16333732 -1.5229131 4198.657
## 5 raceAfrican American/Black -0.70047632 0.18323989 -3.8227283 4176.365
## 6 raceAsian -0.99553834 0.15127139 -6.5811409 4081.378
## 7 raceMultiracial -0.10858234 0.16961133 -0.6401833 4216.146
## 8 raceOther -0.60957057 0.23239650 -2.6229766 4155.811
## 9 ethnicityHispanic/Latino -0.08861526 0.11457586 -0.7734200 3565.931
## 10 studentstatusPart Time 0.15063741 0.17844730 0.8441563 4233.486
## 11 residenceOn Campus 0.13054777 0.07846147 1.6638456 4240.290
## 12 greekmemGreek 0.22861437 0.11132248 2.0536226 4243.389
## 13 (Intercept) -2.65844908 0.12054919 -22.0528152 4241.146
## 14 age21+ 0.40639876 0.11370971 3.5740022 4249.840
## 15 genderMale 1.53721204 0.11087306 13.8646130 4243.201
## 16 genderOther -0.86410415 0.42535722 -2.0314788 4250.241
## 17 raceAfrican American/Black -0.39001078 0.24831054 -1.5706573 4193.262
## 18 raceAsian -1.43340074 0.24933573 -5.7488781 4172.402
## 19 raceMultiracial -0.70259695 0.32847623 -2.1389583 4232.548
## 20 raceOther -0.53396175 0.35990297 -1.4836270 4195.849
## 21 ethnicityHispanic/Latino -0.34353075 0.19146648 -1.7942083 2034.473
## 22 studentstatusPart Time -0.05058380 0.28244369 -0.1790934 4169.033
## 23 residenceOn Campus 0.22379630 0.12002649 1.8645575 4246.924
## 24 greekmemGreek 0.69609817 0.14631077 4.7576687 4249.674
## p.value
## 1 8.898142e-39
## 2 7.986511e-02
## 3 1.844604e-02
## 4 1.278557e-01
## 5 1.339057e-04
## 6 5.260185e-11
## 7 5.220882e-01
## 8 8.748361e-03
## 9 4.393251e-01
## 10 3.986298e-01
## 11 9.621719e-02
## 12 4.007338e-02
## 13 4.027352e-102
## 14 3.554702e-04
## 15 8.782323e-43
## 16 4.226855e-02
## 17 1.163377e-01
## 18 9.625515e-09
## 19 3.249611e-02
## 20 1.379830e-01
## 21 7.292835e-02
## 22 8.578731e-01
## 23 6.231239e-02
## 24 2.023408e-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.0037809548 0.0033097732
## 2 age21+ 0.0009531585 0.0004830991
## 3 genderMale 0.0041900452 0.0037185633
## 4 genderOther 0.0060890688 0.0056157392
## 5 raceAfrican American/Black 0.0074981903 0.0070230095
## 6 raceAsian 0.0120769156 0.0115929217
## 7 raceMultiracial 0.0047882815 0.0043162980
## 8 raceOther 0.0086407528 0.0081637721
## 9 ethnicityHispanic/Latino 0.0276594412 0.0271142439
## 10 studentstatusPart Time 0.0031864818 0.0027156745
## 11 residenceOn Campus 0.0023975448 0.0019271207
## 12 greekmemGreek 0.0019845348 0.0015142596
## 13 (Intercept) 0.0022876462 0.0018172653
## 14 age21+ 0.0009308272 0.0004607699
## 15 genderMale 0.0020107628 0.0015404792
## 16 genderOther 0.0008516007 0.0003815505
## 17 raceAfrican American/Black 0.0064508229 0.0059770570
## 18 raceAsian 0.0077282949 0.0072527731
## 19 raceMultiracial 0.0032857681 0.0028149034
## 20 raceOther 0.0062793275 0.0058057719
## 21 ethnicityHispanic/Latino 0.0683499910 0.0674345774
## 22 studentstatusPart Time 0.0079199042 0.0074440902
## 23 residenceOn Campus 0.0014498680 0.0009797325
## 24 greekmemGreek 0.0009630307 0.0004929703
# 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.09391018 0.08998732 -12.1562704 4201.510
## 2 age21+ 0.02003767 0.08137022 0.2462531 4238.042
## 3 genderMale 0.16758606 0.08356717 2.0054055 4207.383
## 4 genderOther -0.25667621 0.16467136 -1.5587180 4178.642
## 5 raceAfrican American/Black -0.63010436 0.18451679 -3.4148890 4163.961
## 6 raceAsian -0.86960628 0.15351589 -5.6646010 4041.144
## 7 raceMultiracial -0.10521236 0.17059296 -0.6167450 4206.496
## 8 raceOther -0.58998470 0.23392665 -2.5220927 4136.594
## 9 ethnicityHispanic/Latino -0.08086116 0.11504255 -0.7028804 3683.831
## 10 studentstatusPart Time 0.14695666 0.17937740 0.8192597 4223.068
## 11 residenceOn Campus 0.12134636 0.07935377 1.5291821 4231.958
## 12 greekmemGreek 0.13531343 0.11352811 1.1918936 4233.434
## 13 alcmo_dich2 0.13295256 0.09649157 1.3778671 4177.230
## 14 marmo_dich2 0.02193542 0.09431869 0.2325671 4227.380
## 15 ecigmo_dich2 0.26862104 0.10250190 2.6206445 4211.965
## 16 smokmo_dich2 0.06653760 0.11243493 0.5917877 4235.996
## 17 cigarmo_dich2 0.35080667 0.12486257 2.8095422 4221.731
## 18 (Intercept) -3.07272667 0.16247921 -18.9115064 4220.502
## 19 age21+ 0.11018085 0.12295552 0.8961033 4239.989
## 20 genderMale 1.39385017 0.11935335 11.6783497 4227.015
## 21 genderOther -0.91366287 0.42800749 -2.1346890 4239.090
## 22 raceAfrican American/Black -0.20757871 0.25573237 -0.8117029 4192.403
## 23 raceAsian -1.07171844 0.25603159 -4.1858836 4158.296
## 24 raceMultiracial -0.64779115 0.33237044 -1.9490035 4227.408
## 25 raceOther -0.52073313 0.37141396 -1.4020290 4190.209
## 26 ethnicityHispanic/Latino -0.31465395 0.19601469 -1.6052570 2275.599
## 27 studentstatusPart Time -0.03086130 0.28873954 -0.1068828 4129.931
## 28 residenceOn Campus 0.17825784 0.12433378 1.4337041 4236.558
## 29 greekmemGreek 0.37317876 0.15543836 2.4008150 4240.694
## 30 alcmo_dich2 0.21236628 0.16324436 1.3009103 4187.090
## 31 marmo_dich2 0.21045493 0.14211359 1.4808923 4236.237
## 32 ecigmo_dich2 0.43809140 0.15401281 2.8445127 4239.724
## 33 smokmo_dich2 0.20439206 0.15790530 1.2943965 4235.720
## 34 cigarmo_dich2 0.90764210 0.15345361 5.9147654 4231.433
## p.value
## 1 1.924143e-33
## 2 8.054982e-01
## 3 4.498364e-02
## 4 1.191389e-01
## 5 6.441792e-04
## 6 1.576208e-08
## 7 5.374363e-01
## 8 1.170313e-02
## 9 4.821747e-01
## 10 4.126845e-01
## 11 1.262941e-01
## 12 2.333698e-01
## 13 1.683181e-01
## 14 8.161089e-01
## 15 8.807900e-03
## 16 5.540244e-01
## 17 4.984032e-03
## 18 1.246520e-76
## 19 3.702485e-01
## 20 4.905211e-31
## 21 3.284362e-02
## 22 4.170082e-01
## 23 2.899327e-05
## 24 5.136106e-02
## 25 1.609807e-01
## 26 1.085760e-01
## 27 9.148871e-01
## 28 1.517305e-01
## 29 1.640142e-02
## 30 1.933607e-01
## 31 1.387096e-01
## 32 4.469257e-03
## 33 1.955991e-01
## 34 3.585862e-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.0051680500 0.0046946034
## 2 age21+ 0.0012643350 0.0007931268
## 3 genderMale 0.0046978703 0.0042248610
## 4 genderOther 0.0067651621 0.0062898894
## 5 raceAfrican American/Black 0.0076565626 0.0071800427
## 6 raceAsian 0.0133115637 0.0128233632
## 7 raceMultiracial 0.0047710319 0.0042979576
## 8 raceOther 0.0091394334 0.0086604785
## 9 ethnicityHispanic/Latino 0.0243871350 0.0238576058
## 10 studentstatusPart Time 0.0032386493 0.0027667055
## 11 residenceOn Campus 0.0021856481 0.0017141979
## 12 greekmemGreek 0.0019828094 0.0015114278
## 13 alcmo_dich2 0.0068546313 0.0063792408
## 14 marmo_dich2 0.0027580648 0.0022863750
## 15 ecigmo_dich2 0.0043068662 0.0038341858
## 16 smokmo_dich2 0.0016021540 0.0011308777
## 17 cigarmo_dich2 0.0033787000 0.0029066732
## 18 (Intercept) 0.0035042386 0.0030321337
## 19 age21+ 0.0009037642 0.0004326023
## 20 genderMale 0.0028006313 0.0023289209
## 21 genderOther 0.0010759136 0.0006047331
## 22 raceAfrican American/Black 0.0058421747 0.0053680215
## 23 raceAsian 0.0079803217 0.0075033085
## 24 raceMultiracial 0.0027548137 0.0022831254
## 25 raceOther 0.0059963010 0.0055219730
## 26 ethnicityHispanic/Latino 0.0607014166 0.0598762397
## 27 studentstatusPart Time 0.0094736049 0.0089940392
## 28 residenceOn Campus 0.0015128570 0.0010416010
## 29 greekmemGreek 0.0007605937 0.0002894426
## 30 alcmo_dich2 0.0062104623 0.0057358835
## 31 marmo_dich2 0.0015641560 0.0010928886
## 32 ecigmo_dich2 0.0009556550 0.0004844881
## 33 smokmo_dich2 0.0016451309 0.0011738442
## 34 cigarmo_dich2 0.0022553695 0.0017838938
# Fit model 1d: model c + rx stimulants + rx painkillers + rx sedatives
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,
trace = FALSE
)
)
# Pool model 1d
pooled1d <- pool(fit1d)
summary(pooled1d)
## term estimate std.error statistic df
## 1 (Intercept) -1.090908708 0.09002810 -12.11742516 4196.241
## 2 age21+ 0.015494403 0.08149472 0.19012770 4232.063
## 3 genderMale 0.164376447 0.08364893 1.96507532 4203.827
## 4 genderOther -0.266380280 0.16541511 -1.61037457 4172.674
## 5 raceAfrican American/Black -0.621954853 0.18464321 -3.36841447 4156.666
## 6 raceAsian -0.876627763 0.15367691 -5.70435580 4036.985
## 7 raceMultiracial -0.098626851 0.17074527 -0.57762566 4202.190
## 8 raceOther -0.573230846 0.23414445 -2.44819314 4128.876
## 9 ethnicityHispanic/Latino -0.090302150 0.11525829 -0.78347637 3702.948
## 10 studentstatusPart Time 0.165394495 0.17962248 0.92078951 4215.443
## 11 residenceOn Campus 0.121332255 0.07941679 1.52779099 4225.757
## 12 greekmemGreek 0.127828646 0.11383344 1.12294456 4226.632
## 13 alcmo_dich2 0.132115844 0.09658518 1.36786875 4173.174
## 14 marmo_dich2 0.003592934 0.09515458 0.03775892 4222.446
## 15 ecigmo_dich2 0.262486438 0.10272808 2.55515765 4205.747
## 16 smokmo_dich2 0.049901681 0.11312589 0.44111635 4229.745
## 17 cigarmo_dich2 0.332929337 0.12581169 2.64625116 4212.676
## 18 rxstmo_dich2 0.389472785 0.21716797 1.79341722 4224.035
## 19 rxpkmo_dich2 0.435981142 0.26961321 1.61706150 4232.912
## 20 rxsedmo_dich2 -0.453263070 0.32139251 -1.41031001 4233.689
## 21 (Intercept) -3.064394194 0.16251017 -18.85663055 4214.505
## 22 age21+ 0.098509807 0.12330302 0.79892453 4233.954
## 23 genderMale 1.386777539 0.11956738 11.59829282 4220.736
## 24 genderOther -0.928943838 0.42938063 -2.16345074 4233.208
## 25 raceAfrican American/Black -0.192481202 0.25575653 -0.75259545 4185.923
## 26 raceAsian -1.085223209 0.25710114 -4.22099722 4148.150
## 27 raceMultiracial -0.632921787 0.33229672 -1.90468865 4221.487
## 28 raceOther -0.473650780 0.37055825 -1.27820872 4178.660
## 29 ethnicityHispanic/Latino -0.322998148 0.19650531 -1.64371207 2207.853
## 30 studentstatusPart Time 0.004254594 0.28871527 0.01473630 4109.103
## 31 residenceOn Campus 0.178678793 0.12455909 1.43449015 4230.478
## 32 greekmemGreek 0.349343150 0.15665751 2.22998020 4234.699
## 33 alcmo_dich2 0.220230031 0.16335915 1.34813404 4186.809
## 34 marmo_dich2 0.176256002 0.14366624 1.22684357 4231.507
## 35 ecigmo_dich2 0.426986634 0.15459517 2.76196619 4233.839
## 36 smokmo_dich2 0.179646721 0.15942470 1.12684371 4228.856
## 37 cigarmo_dich2 0.885643729 0.15442946 5.73494011 4225.242
## 38 rxstmo_dich2 0.594528880 0.26333122 2.25772269 4211.259
## 39 rxpkmo_dich2 0.336992361 0.36360395 0.92681161 4133.457
## 40 rxsedmo_dich2 -0.598452439 0.42436781 -1.41022108 4225.325
## p.value
## 1 3.048289e-33
## 2 8.492182e-01
## 3 4.947127e-02
## 4 1.073917e-01
## 5 7.628971e-04
## 6 1.251496e-08
## 7 5.635478e-01
## 8 1.439880e-02
## 9 4.333975e-01
## 10 3.572130e-01
## 11 1.266393e-01
## 12 2.615248e-01
## 13 1.714268e-01
## 14 9.698817e-01
## 15 1.064887e-02
## 16 6.591513e-01
## 17 8.169200e-03
## 18 7.297775e-02
## 19 1.059395e-01
## 20 1.585216e-01
## 21 3.281104e-76
## 22 4.243790e-01
## 23 1.219729e-30
## 24 3.056231e-02
## 25 4.517354e-01
## 26 2.484300e-05
## 27 5.688845e-02
## 28 2.012468e-01
## 29 1.003781e-01
## 30 9.882433e-01
## 31 1.515064e-01
## 32 2.580098e-02
## 33 1.776882e-01
## 34 2.199496e-01
## 35 5.770230e-03
## 36 2.598725e-01
## 37 1.043588e-08
## 38 2.401395e-02
## 39 3.540785e-01
## 40 1.585480e-01
## 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.0051185078 0.0046444433
## 2 age21+ 0.0012624553 0.0007905807
## 3 genderMale 0.0044986640 0.0040251600
## 4 genderOther 0.0067728025 0.0062968538
## 5 raceAfrican American/Black 0.0077425356 0.0072652210
## 6 raceAsian 0.0132596845 0.0127709555
## 7 raceMultiracial 0.0046372693 0.0041636468
## 8 raceOther 0.0092396893 0.0087598878
## 9 ethnicityHispanic/Latino 0.0237450033 0.0232178603
## 10 studentstatusPart Time 0.0034132684 0.0029405540
## 11 residenceOn Campus 0.0022156456 0.0017435180
## 12 greekmemGreek 0.0020973156 0.0016252298
## 13 alcmo_dich2 0.0067408480 0.0062649411
## 14 marmo_dich2 0.0026350314 0.0021627323
## 15 ecigmo_dich2 0.0043321602 0.0038587932
## 16 smokmo_dich2 0.0016436533 0.0011717003
## 17 cigarmo_dich2 0.0036908028 0.0032179097
## 18 rxstmo_dich2 0.0024389806 0.0019667663
## 19 rxpkmo_dich2 0.0011103901 0.0006385382
## 20 rxsedmo_dich2 0.0009638270 0.0004919925
## 21 (Intercept) 0.0035089477 0.0030361734
## 22 age21+ 0.0009120390 0.0004402096
## 23 genderMale 0.0028369516 0.0023645568
## 24 genderOther 0.0010553926 0.0005835477
## 25 raceAfrican American/Black 0.0058845283 0.0054096615
## 26 raceAsian 0.0082229296 0.0077448669
## 27 raceMultiracial 0.0027494286 0.0022770763
## 28 raceOther 0.0063827480 0.0059072944
## 29 ethnicityHispanic/Latino 0.0626920993 0.0618434163
## 30 studentstatusPart Time 0.0102022086 0.0097205673
## 31 residenceOn Campus 0.0015277999 0.0010558740
## 32 greekmemGreek 0.0007606555 0.0002888376
## 33 alcmo_dich2 0.0058216248 0.0053468284
## 34 marmo_dich2 0.0013582003 0.0008863088
## 35 ecigmo_dich2 0.0009346238 0.0004627922
## 36 smokmo_dich2 0.0017790454 0.0013070572
## 37 cigarmo_dich2 0.0022837041 0.0018115511
## 38 rxstmo_dich2 0.0038278288 0.0033548417
## 39 rxpkmo_dich2 0.0090059047 0.0085265217
## 40 rxsedmo_dich2 0.0022728234 0.0018006746
# 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.69553038 0.1337449 -12.67735039 1782.421
## 2 age21+ -0.46065146 0.1365097 -3.37449669 1789.292
## 3 genderMale 0.93751927 0.1339963 6.99660593 1777.955
## 4 genderOther 1.57424610 0.2564040 6.13970977 1785.308
## 5 raceAfrican American/Black 0.80407365 0.2935423 2.73920860 1775.599
## 6 raceAsian 1.31705202 0.2281164 5.77359630 1772.676
## 7 raceMultiracial 0.36168976 0.2894691 1.24949355 1782.796
## 8 raceOther 0.85104517 0.3770052 2.25738294 1754.645
## 9 ethnicityHispanic/Latino 0.31515295 0.2034335 1.54916958 1226.246
## 10 studentstatusPart Time 0.20969294 0.3079477 0.68093691 1789.626
## 11 residenceOn Campus 0.14385465 0.1387134 1.03706368 1785.184
## 12 greekmemGreek -0.24553465 0.2097869 -1.17040023 1785.367
## 13 (Intercept) -2.46410271 0.1606267 -15.34055497 1785.661
## 14 age21+ 0.21131531 0.1410668 1.49798063 1789.337
## 15 genderMale 1.85355336 0.1428080 12.97933557 1789.361
## 16 genderOther 0.69285263 0.4070286 1.70222094 1789.548
## 17 raceAfrican American/Black 0.68314536 0.3143029 2.17352567 1784.085
## 18 raceAsian 0.24900578 0.2967702 0.83905248 1767.315
## 19 raceMultiracial -0.56845650 0.4056802 -1.40124299 1787.676
## 20 raceOther 0.08277533 0.5082090 0.16287656 1779.402
## 21 ethnicityHispanic/Latino -0.02080925 0.2338592 -0.08898196 1159.716
## 22 studentstatusPart Time -0.26025792 0.3642008 -0.71460009 1777.236
## 23 residenceOn Campus 0.16750155 0.1493371 1.12163355 1788.154
## 24 greekmemGreek 0.50500639 0.1788348 2.82387058 1789.559
## p.value
## 1 2.516654e-35
## 2 7.553022e-04
## 3 3.700370e-12
## 4 1.016407e-09
## 5 6.220307e-03
## 6 9.142378e-09
## 7 2.116486e-01
## 8 2.410665e-02
## 9 1.215990e-01
## 10 4.959995e-01
## 11 2.998466e-01
## 12 2.419961e-01
## 13 5.460516e-50
## 14 1.343147e-01
## 15 7.004011e-37
## 16 8.888763e-02
## 17 2.987191e-02
## 18 4.015533e-01
## 19 1.613150e-01
## 20 8.706341e-01
## 21 9.291116e-01
## 22 4.749500e-01
## 23 2.621689e-01
## 24 4.797302e-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.004365621 0.0032490762
## 2 age21+ 0.001498888 0.0003834264
## 3 genderMale 0.005803563 0.0046858317
## 4 genderOther 0.003294070 0.0021781301
## 5 raceAfrican American/Black 0.006486926 0.0053684820
## 6 raceAsian 0.007280021 0.0061606282
## 7 raceMultiracial 0.004233731 0.0031172735
## 8 raceOther 0.011341280 0.0102150176
## 9 ethnicityHispanic/Latino 0.065936098 0.0644138862
## 10 studentstatusPart Time 0.001322362 0.0002069111
## 11 residenceOn Campus 0.003343179 0.0022272170
## 12 greekmemGreek 0.003270892 0.0021549631
## 13 (Intercept) 0.003152515 0.0020366379
## 14 age21+ 0.001475355 0.0003598956
## 15 genderMale 0.001462737 0.0003472782
## 16 genderOther 0.001363929 0.0002484759
## 17 raceAfrican American/Black 0.003765182 0.0026490056
## 18 raceAsian 0.008613834 0.0074925561
## 19 raceMultiracial 0.002286947 0.0011713588
## 20 raceOther 0.005360528 0.0042432081
## 21 ethnicityHispanic/Latino 0.072056412 0.0704574965
## 22 studentstatusPart Time 0.006016721 0.0048987780
## 23 residenceOn Campus 0.002064044 0.0009485047
## 24 greekmemGreek 0.001358029 0.0002425761
# 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.552980180 0.1648584 -9.42008303 1766.883
## 2 age21+ -0.342592416 0.1464198 -2.33979575 1776.886
## 3 genderMale 0.984652973 0.1392662 7.07029276 1761.494
## 4 genderOther 1.588015330 0.2595714 6.11783706 1773.903
## 5 raceAfrican American/Black 0.717451329 0.2955320 2.42766063 1764.223
## 6 raceAsian 1.196452968 0.2329341 5.13644354 1762.144
## 7 raceMultiracial 0.349819477 0.2908594 1.20270979 1773.911
## 8 raceOther 0.891630241 0.3811312 2.33943122 1746.207
## 9 ethnicityHispanic/Latino 0.272927560 0.2055564 1.32775011 1206.413
## 10 studentstatusPart Time 0.232145037 0.3075425 0.75483893 1779.712
## 11 residenceOn Campus 0.147941592 0.1393729 1.06148001 1775.433
## 12 greekmemGreek -0.113418788 0.2134077 -0.53146520 1776.536
## 13 alcmo_dich2 -0.126063327 0.1698746 -0.74209621 1746.022
## 14 marmo_dich2 -0.069911425 0.1695462 -0.41234432 1776.765
## 15 ecigmo_dich2 0.052653478 0.1850924 0.28447134 1751.810
## 16 smokmo_dich2 -0.271071767 0.2072913 -1.30768522 1776.358
## 17 cigarmo_dich2 -0.370148703 0.2229399 -1.66030740 1768.251
## 18 (Intercept) -2.553407637 0.2068613 -12.34357625 1765.346
## 19 age21+ 0.080813894 0.1525908 0.52961192 1777.330
## 20 genderMale 1.728192269 0.1512220 11.42818284 1778.299
## 21 genderOther 0.660486924 0.4110541 1.60681264 1779.586
## 22 raceAfrican American/Black 0.809971527 0.3206124 2.52632595 1773.981
## 23 raceAsian 0.487939580 0.3049762 1.59992695 1751.964
## 24 raceMultiracial -0.546955936 0.4090845 -1.33702424 1777.281
## 25 raceOther 0.007090684 0.5124699 0.01383629 1767.246
## 26 ethnicityHispanic/Latino 0.017626355 0.2375866 0.07418918 1097.070
## 27 studentstatusPart Time -0.243019702 0.3696966 -0.65734896 1766.747
## 28 residenceOn Campus 0.094671504 0.1530654 0.61850351 1778.122
## 29 greekmemGreek 0.313165879 0.1881889 1.66410429 1779.305
## 30 alcmo_dich2 -0.176395479 0.2078345 -0.84873047 1762.447
## 31 marmo_dich2 0.072930322 0.1762710 0.41373968 1779.039
## 32 ecigmo_dich2 0.418618211 0.1910554 2.19108307 1776.885
## 33 smokmo_dich2 0.051787149 0.1921068 0.26957473 1777.289
## 34 cigarmo_dich2 0.587393115 0.1826658 3.21567165 1775.921
## p.value
## 1 1.357579e-20
## 2 1.940449e-02
## 3 2.220369e-12
## 4 1.164374e-09
## 5 1.529603e-02
## 6 3.110150e-07
## 7 2.292491e-01
## 8 1.942532e-02
## 9 1.845118e-01
## 10 4.504455e-01
## 11 2.886162e-01
## 12 5.951629e-01
## 13 4.581289e-01
## 14 6.801368e-01
## 15 7.760828e-01
## 16 1.911493e-01
## 17 9.702993e-02
## 18 1.232567e-33
## 19 5.964472e-01
## 20 3.082899e-29
## 21 1.082729e-01
## 22 1.161261e-02
## 23 1.097951e-01
## 24 1.813858e-01
## 25 9.889621e-01
## 26 9.408734e-01
## 27 5.110422e-01
## 28 5.363227e-01
## 29 9.626767e-02
## 30 3.961465e-01
## 31 6.791146e-01
## 32 2.857524e-02
## 33 7.875187e-01
## 34 1.324795e-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.006154264 0.0050299300
## 2 age21+ 0.002654251 0.0015323049
## 3 genderMale 0.007629159 0.0065030604
## 4 genderOther 0.003854638 0.0027321593
## 5 raceAfrican American/Black 0.006905726 0.0057805487
## 6 raceAsian 0.007460718 0.0063348439
## 7 raceMultiracial 0.003851593 0.0027291157
## 8 raceOther 0.011092104 0.0099601174
## 9 ethnicityHispanic/Latino 0.067300198 0.0657552417
## 10 studentstatusPart Time 0.001282998 0.0001612930
## 11 residenceOn Campus 0.003262688 0.0021405096
## 12 greekmemGreek 0.002805605 0.0016836084
## 13 alcmo_dich2 0.011129574 0.0099975100
## 14 marmo_dich2 0.002706803 0.0015848405
## 15 ecigmo_dich2 0.009914926 0.0087852147
## 16 smokmo_dich2 0.002881203 0.0017591798
## 17 cigarmo_dich2 0.005746643 0.0046227176
## 18 (Intercept) 0.006594531 0.0054697173
## 19 age21+ 0.002457364 0.0013354764
## 20 genderMale 0.002005945 0.0008841620
## 21 genderOther 0.001351254 0.0002295463
## 22 raceAfrican American/Black 0.003825487 0.0027030243
## 23 raceAsian 0.009881361 0.0087517105
## 24 raceMultiracial 0.002479276 0.0013573824
## 25 raceOther 0.006047691 0.0049234674
## 26 ethnicityHispanic/Latino 0.077623815 0.0759438200
## 27 studentstatusPart Time 0.006193953 0.0050695776
## 28 residenceOn Campus 0.002091042 0.0009692421
## 29 greekmemGreek 0.001500679 0.0003789623
## 30 alcmo_dich2 0.007381402 0.0062556312
## 31 marmo_dich2 0.001638372 0.0005166418
## 32 ecigmo_dich2 0.002654439 0.0015324927
## 33 smokmo_dich2 0.002475750 0.0013538576
## 34 cigarmo_dich2 0.003063901 0.0019418072
# 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,
trace = FALSE
)
)
# Pool model 2d
pooled2d <- pool(fit2d)
summary(pooled2d)
## term estimate std.error statistic df
## 1 (Intercept) -1.550232652 0.1648570 -9.40350115 1761.103
## 2 age21+ -0.337691880 0.1465558 -2.30418684 1770.905
## 3 genderMale 0.983158560 0.1393331 7.05617313 1755.327
## 4 genderOther 1.571454531 0.2595088 6.05549670 1767.706
## 5 raceAfrican American/Black 0.712204452 0.2955735 2.40956800 1758.379
## 6 raceAsian 1.188301636 0.2347900 5.06112450 1755.630
## 7 raceMultiracial 0.329695131 0.2911105 1.13254304 1768.074
## 8 raceOther 0.877350784 0.3806562 2.30483753 1739.390
## 9 ethnicityHispanic/Latino 0.274563733 0.2059396 1.33322443 1203.313
## 10 studentstatusPart Time 0.213572490 0.3077810 0.69391048 1773.697
## 11 residenceOn Campus 0.142557858 0.1394099 1.02258042 1769.530
## 12 greekmemGreek -0.098216421 0.2137263 -0.45954292 1770.599
## 13 alcmo_dich2 -0.129523578 0.1700773 -0.76155700 1740.775
## 14 marmo_dich2 -0.049588196 0.1709694 -0.29004139 1770.879
## 15 ecigmo_dich2 0.049101772 0.1854019 0.26483964 1746.788
## 16 smokmo_dich2 -0.276990288 0.2086662 -1.32743269 1770.681
## 17 cigarmo_dich2 -0.348541297 0.2232032 -1.56154252 1762.118
## 18 rxstmo_dich2 -0.564511233 0.4425343 -1.27563280 1773.179
## 19 rxpkmo_dich2 0.081113232 0.4323575 0.18760685 1757.426
## 20 rxsedmo_dich2 0.700750782 0.5534528 1.26614377 1772.390
## 21 (Intercept) -2.546982931 0.2069674 -12.30620203 1759.834
## 22 age21+ 0.080367911 0.1528104 0.52593205 1771.396
## 23 genderMale 1.722360157 0.1514674 11.37116006 1772.420
## 24 genderOther 0.666576279 0.4121117 1.61746523 1773.593
## 25 raceAfrican American/Black 0.812620687 0.3202248 2.53765712 1768.310
## 26 raceAsian 0.476600677 0.3050378 1.56243142 1744.806
## 27 raceMultiracial -0.543701671 0.4093625 -1.32816690 1770.923
## 28 raceOther 0.016503964 0.5123685 0.03221112 1760.896
## 29 ethnicityHispanic/Latino 0.007607376 0.2385131 0.03189501 1072.391
## 30 studentstatusPart Time -0.228576926 0.3702659 -0.61733179 1760.806
## 31 residenceOn Campus 0.091558465 0.1532576 0.59741533 1772.383
## 32 greekmemGreek 0.306658472 0.1892395 1.62047771 1773.498
## 33 alcmo_dich2 -0.175858593 0.2080161 -0.84540862 1757.014
## 34 marmo_dich2 0.067878931 0.1780357 0.38126589 1772.757
## 35 ecigmo_dich2 0.413714173 0.1912927 2.16272878 1770.487
## 36 smokmo_dich2 0.047935689 0.1937949 0.24735275 1771.211
## 37 cigarmo_dich2 0.586194264 0.1836836 3.19132633 1769.846
## 38 rxstmo_dich2 0.044887378 0.3010419 0.14910674 1765.213
## 39 rxpkmo_dich2 0.234720210 0.4022244 0.58355531 1540.834
## 40 rxsedmo_dich2 -0.295769190 0.5222601 -0.56632548 1759.125
## p.value
## 1 1.583182e-20
## 2 2.132741e-02
## 3 2.454116e-12
## 4 1.705895e-09
## 5 1.607359e-02
## 6 4.604719e-07
## 7 2.575598e-01
## 8 2.129293e-02
## 9 1.827105e-01
## 10 4.878292e-01
## 11 3.066460e-01
## 12 6.459008e-01
## 13 4.464276e-01
## 14 7.718185e-01
## 15 7.911642e-01
## 16 1.845367e-01
## 17 1.185754e-01
## 18 2.022524e-01
## 19 8.512065e-01
## 20 2.056280e-01
## 21 1.906415e-33
## 22 5.990012e-01
## 23 5.716972e-29
## 24 1.059557e-01
## 25 1.124496e-02
## 26 1.183677e-01
## 27 1.842941e-01
## 28 9.743073e-01
## 29 9.745617e-01
## 30 5.370957e-01
## 31 5.503064e-01
## 32 1.053075e-01
## 33 3.979978e-01
## 34 7.030516e-01
## 35 3.069555e-02
## 36 8.046639e-01
## 37 1.441053e-03
## 38 8.814864e-01
## 39 5.596049e-01
## 40 5.712448e-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.006110406 0.0049823339
## 2 age21+ 0.002654715 0.0015289830
## 3 genderMale 0.007698205 0.0065682311
## 4 genderOther 0.003941110 0.0028147964
## 5 raceAfrican American/Black 0.006886230 0.0057572934
## 6 raceAsian 0.007620050 0.0064901823
## 7 raceMultiracial 0.003803488 0.0026772531
## 8 raceOther 0.011295210 0.0101590229
## 9 ethnicityHispanic/Latino 0.067317579 0.0657686745
## 10 studentstatusPart Time 0.001295336 0.0001698430
## 11 residenceOn Campus 0.003234398 0.0021084464
## 12 greekmemGreek 0.002787882 0.0016621059
## 13 alcmo_dich2 0.011014132 0.0098785248
## 14 marmo_dich2 0.002666277 0.0015405408
## 15 ecigmo_dich2 0.009733189 0.0086000237
## 16 smokmo_dich2 0.002752344 0.0016265795
## 17 cigarmo_dich2 0.005806542 0.0046787749
## 18 rxstmo_dich2 0.001571766 0.0004462556
## 19 rxpkmo_dich2 0.007145869 0.0060166158
## 20 rxsedmo_dich2 0.001968635 0.0008430722
## 21 (Intercept) 0.006478607 0.0053501404
## 22 age21+ 0.002435819 0.0013101512
## 23 genderMale 0.001953892 0.0008283309
## 24 genderOther 0.001352250 0.0002267553
## 25 raceAfrican American/Black 0.003713989 0.0025878029
## 26 raceAsian 0.010166918 0.0090329627
## 27 raceMultiracial 0.002646666 0.0015209360
## 28 raceOther 0.006171461 0.0050433260
## 29 ethnicityHispanic/Latino 0.079815985 0.0781014480
## 30 studentstatusPart Time 0.006197700 0.0050695366
## 31 residenceOn Campus 0.001971951 0.0008463876
## 32 greekmemGreek 0.001403018 0.0002775202
## 33 alcmo_dich2 0.007256219 0.0061268268
## 34 marmo_dich2 0.001787364 0.0006618287
## 35 ecigmo_dich2 0.002835980 0.0017101868
## 36 smokmo_dich2 0.002518926 0.0013932354
## 37 cigarmo_dich2 0.003105188 0.0019792911
## 38 rxstmo_dich2 0.004819638 0.0036927292
## 39 rxpkmo_dich2 0.036339355 0.0350893370
## 40 rxsedmo_dich2 0.006679030 0.0055503366
- 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,
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,
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,
data = df6_model1
)
vif(vif_check1d)
## GVIF Df GVIF^(1/(2*Df))
## age 1.281642 1 1.132096
## gender 1.122064 2 1.029211
## race 1.238655 4 1.027114
## ethnicity 1.151341 1 1.073006
## studentstatus 1.034848 1 1.017275
## residence 1.226269 1 1.107370
## greekmem 1.110929 1 1.054006
## alcmo_dich 1.435078 1 1.197947
## marmo_dich 1.683316 1 1.297427
## ecigmo_dich 1.899071 1 1.378068
## smokmo_dich 1.785934 1 1.336388
## cigarmo_dich 1.472210 1 1.213347
## rxstmo_dich 1.365110 1 1.168379
## rxpkmo_dich 1.364695 1 1.168202
## rxsedmo_dich 1.369998 1 1.170469
# 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,
data = df6_model2
)
vif(vif_check2d)
## GVIF Df GVIF^(1/(2*Df))
## age 1.290361 1 1.135941
## gender 1.153041 2 1.036242
## race 1.199863 4 1.023037
## ethnicity 1.123367 1 1.059890
## studentstatus 1.034097 1 1.016905
## residence 1.225814 1 1.107165
## greekmem 1.158049 1 1.076127
## alcmo_dich 1.459840 1 1.208238
## marmo_dich 1.685732 1 1.298357
## ecigmo_dich 1.954899 1 1.398177
## smokmo_dich 1.863934 1 1.365260
## cigarmo_dich 1.576686 1 1.255662
## rxstmo_dich 1.507410 1 1.227766
## rxpkmo_dich 1.526068 1 1.235341
## rxsedmo_dich 1.531515 1 1.237544
### 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.122 | df: 22 | p: 0.0000
## SD across imputations: 2.097
##
## --- Model 1: B → C ---
## Mean LRT (chi-sq): 155.148 | df: 10 | p: 0.0000
## SD across imputations: 0.827
##
## --- Model 1: C → D ---
## Mean LRT (chi-sq): 11.350 | df: 6 | p: 0.0781
## SD across imputations: 0.294
##
## --- Model 1: B → D ---
## Mean LRT (chi-sq): 166.498 | df: 16 | p: 0.0000
## SD across imputations: 0.874
##
## --- Model 1: A → C ---
## Mean LRT (chi-sq): 534.270 | df: 32 | p: 0.0000
## SD across imputations: 2.195
##
## --- Model 1: A → D ---
## Mean LRT (chi-sq): 545.620 | df: 38 | p: 0.0000
## SD across imputations: 2.190
##
## --- Model 2: A → B ---
## Mean LRT (chi-sq): 328.288 | df: 22 | p: 0.0000
## SD across imputations: 1.395
##
## --- Model 2: B → C ---
## Mean LRT (chi-sq): 59.780 | df: 10 | p: 0.0000
## SD across imputations: 0.694
##
## --- Model 2: A → C ---
## Mean LRT (chi-sq): 388.068 | df: 32 | p: 0.0000
## SD across imputations: 1.454
##
## --- Model 2: C → D ---
## Mean LRT (chi-sq): 3.686 | df: 6 | p: 0.7191
## SD across imputations: 0.207
##
## --- Model 2: B → D ---
## Mean LRT (chi-sq): 63.466 | df: 16 | p: 0.0000
## SD across imputations: 0.686
##
## --- Model 2: A → D ---
## Mean LRT (chi-sq): 391.754 | df: 38 | p: 0.0000
## SD across imputations: 1.424
print(lrt_table)
## comparison mean_LRT_chisq df p_value sd_across_imp
## 1 Model 1: A → B 379.122 22 0.0000 2.097
## 2 Model 1: B → C 155.148 10 0.0000 0.827
## 3 Model 1: C → D 11.350 6 0.0781 0.294
## 4 Model 1: B → D 166.498 16 0.0000 0.874
## 5 Model 1: A → C 534.270 32 0.0000 2.195
## 6 Model 1: A → D 545.620 38 0.0000 2.190
## 7 Model 2: A → B 328.288 22 0.0000 1.395
## 8 Model 2: B → C 59.780 10 0.0000 0.694
## 9 Model 2: A → C 388.068 32 0.0000 1.454
## 10 Model 2: C → D 3.686 6 0.7191 0.207
## 11 Model 2: B → D 63.466 16 0.0000 0.686
## 12 Model 2: A → D 391.754 38 0.0000 1.424
- 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.143 0.984
## Model 1b Class 1 vs. Class 3 genderMale 1.208 1.032
## Model 1b Class 1 vs. Class 3 genderOther 0.780 0.566
## Model 1b Class 1 vs. Class 3 raceAfrican American/Black 0.496 0.347
## Model 1b Class 1 vs. Class 3 raceAsian 0.370 0.275
## Model 1b Class 1 vs. Class 3 raceMultiracial 0.897 0.643
## Model 1b Class 1 vs. Class 3 raceOther 0.544 0.345
## Model 1b Class 1 vs. Class 3 ethnicityHispanic/Latino 0.915 0.731
## Model 1b Class 1 vs. Class 3 studentstatusPart Time 1.163 0.819
## Model 1b Class 1 vs. Class 3 residenceOn Campus 1.139 0.977
## Model 1b Class 1 vs. Class 3 greekmemGreek 1.257 1.010
## Model 1b Class 2 vs. Class 3 (Intercept) 0.070 0.055
## Model 1b Class 2 vs. Class 3 age21+ 1.501 1.201
## Model 1b Class 2 vs. Class 3 genderMale 4.652 3.743
## Model 1b Class 2 vs. Class 3 genderOther 0.421 0.183
## Model 1b Class 2 vs. Class 3 raceAfrican American/Black 0.677 0.416
## Model 1b Class 2 vs. Class 3 raceAsian 0.238 0.146
## Model 1b Class 2 vs. Class 3 raceMultiracial 0.495 0.260
## Model 1b Class 2 vs. Class 3 raceOther 0.586 0.290
## Model 1b Class 2 vs. Class 3 ethnicityHispanic/Latino 0.709 0.487
## Model 1b Class 2 vs. Class 3 studentstatusPart Time 0.951 0.547
## Model 1b Class 2 vs. Class 3 residenceOn Campus 1.251 0.989
## Model 1b Class 2 vs. Class 3 greekmemGreek 2.006 1.506
## CI_upper p_value sig
## 0.462 0.0000 ***
## 1.326 0.0799 .
## 1.413 0.0184 *
## 1.074 0.1279
## 0.711 0.0001 ***
## 0.497 0.0000 ***
## 1.251 0.5221
## 0.857 0.0087 **
## 1.146 0.4393
## 1.649 0.3986
## 1.329 0.0962 .
## 1.563 0.0401 *
## 0.089 0.0000 ***
## 1.876 0.0004 ***
## 5.781 0.0000 ***
## 0.970 0.0423 *
## 1.102 0.1163
## 0.389 0.0000 ***
## 0.943 0.0325 *
## 1.187 0.1380
## 1.032 0.0729 .
## 1.654 0.8579
## 1.583 0.0623 .
## 2.672 0.0000 ***
print(or_1c, row.names = FALSE)
## model contrast term OR CI_lower
## Model 1c Class 1 vs. Class 3 (Intercept) 0.335 0.281
## Model 1c Class 1 vs. Class 3 age21+ 1.020 0.870
## Model 1c Class 1 vs. Class 3 genderMale 1.182 1.004
## Model 1c Class 1 vs. Class 3 genderOther 0.774 0.560
## Model 1c Class 1 vs. Class 3 raceAfrican American/Black 0.533 0.371
## Model 1c Class 1 vs. Class 3 raceAsian 0.419 0.310
## Model 1c Class 1 vs. Class 3 raceMultiracial 0.900 0.644
## Model 1c Class 1 vs. Class 3 raceOther 0.554 0.350
## Model 1c Class 1 vs. Class 3 ethnicityHispanic/Latino 0.922 0.736
## Model 1c Class 1 vs. Class 3 studentstatusPart Time 1.158 0.815
## Model 1c Class 1 vs. Class 3 residenceOn Campus 1.129 0.966
## Model 1c Class 1 vs. Class 3 greekmemGreek 1.145 0.916
## Model 1c Class 1 vs. Class 3 alcmo_dich2 1.142 0.945
## Model 1c Class 1 vs. Class 3 marmo_dich2 1.022 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.420 1.112
## Model 1c Class 2 vs. Class 3 (Intercept) 0.046 0.034
## Model 1c Class 2 vs. Class 3 age21+ 1.116 0.877
## Model 1c Class 2 vs. Class 3 genderMale 4.030 3.190
## Model 1c Class 2 vs. Class 3 genderOther 0.401 0.173
## Model 1c Class 2 vs. Class 3 raceAfrican American/Black 0.813 0.492
## Model 1c Class 2 vs. Class 3 raceAsian 0.342 0.207
## Model 1c Class 2 vs. Class 3 raceMultiracial 0.523 0.273
## Model 1c Class 2 vs. Class 3 raceOther 0.594 0.287
## Model 1c Class 2 vs. Class 3 ethnicityHispanic/Latino 0.730 0.497
## Model 1c Class 2 vs. Class 3 studentstatusPart Time 0.970 0.551
## Model 1c Class 2 vs. Class 3 residenceOn Campus 1.195 0.937
## Model 1c Class 2 vs. Class 3 greekmemGreek 1.452 1.071
## Model 1c Class 2 vs. Class 3 alcmo_dich2 1.237 0.898
## Model 1c Class 2 vs. Class 3 marmo_dich2 1.234 0.934
## Model 1c Class 2 vs. Class 3 ecigmo_dich2 1.550 1.146
## Model 1c Class 2 vs. Class 3 smokmo_dich2 1.227 0.900
## Model 1c Class 2 vs. Class 3 cigarmo_dich2 2.478 1.835
## CI_upper p_value sig
## 0.400 0.0000 ***
## 1.197 0.8055
## 1.393 0.0450 *
## 1.068 0.1191
## 0.765 0.0006 ***
## 0.566 0.0000 ***
## 1.258 0.5374
## 0.877 0.0117 *
## 1.156 0.4822
## 1.646 0.4127
## 1.319 0.1263
## 1.430 0.2334
## 1.380 0.1683
## 1.230 0.8161
## 1.599 0.0088 **
## 1.332 0.5540
## 1.814 0.0050 **
## 0.064 0.0000 ***
## 1.421 0.3702
## 5.093 0.0000 ***
## 0.928 0.0328 *
## 1.341 0.4170
## 0.566 0.0000 ***
## 1.004 0.0514 .
## 1.230 0.1610
## 1.072 0.1086
## 1.708 0.9149
## 1.525 0.1517
## 1.970 0.0164 *
## 1.703 0.1934
## 1.631 0.1387
## 2.096 0.0045 **
## 1.672 0.1956
## 3.348 0.0000 ***
print(or_1d, row.names = FALSE)
## model contrast term OR CI_lower
## Model 1d Class 1 vs. Class 3 (Intercept) 0.336 0.282
## Model 1d Class 1 vs. Class 3 age21+ 1.016 0.866
## Model 1d Class 1 vs. Class 3 genderMale 1.179 1.000
## Model 1d Class 1 vs. Class 3 genderOther 0.766 0.554
## Model 1d Class 1 vs. Class 3 raceAfrican American/Black 0.537 0.374
## Model 1d Class 1 vs. Class 3 raceAsian 0.416 0.308
## Model 1d Class 1 vs. Class 3 raceMultiracial 0.906 0.648
## Model 1d Class 1 vs. Class 3 raceOther 0.564 0.356
## Model 1d Class 1 vs. Class 3 ethnicityHispanic/Latino 0.914 0.729
## Model 1d Class 1 vs. Class 3 studentstatusPart Time 1.180 0.830
## Model 1d Class 1 vs. Class 3 residenceOn Campus 1.129 0.966
## Model 1d Class 1 vs. Class 3 greekmemGreek 1.136 0.909
## Model 1d Class 1 vs. Class 3 alcmo_dich2 1.141 0.944
## Model 1d Class 1 vs. Class 3 marmo_dich2 1.004 0.833
## Model 1d Class 1 vs. Class 3 ecigmo_dich2 1.300 1.063
## Model 1d Class 1 vs. Class 3 smokmo_dich2 1.051 0.842
## Model 1d Class 1 vs. Class 3 cigarmo_dich2 1.395 1.090
## Model 1d Class 1 vs. Class 3 rxstmo_dich2 1.476 0.964
## Model 1d Class 1 vs. Class 3 rxpkmo_dich2 1.546 0.912
## Model 1d Class 1 vs. Class 3 rxsedmo_dich2 0.636 0.339
## Model 1d Class 2 vs. Class 3 (Intercept) 0.047 0.034
## Model 1d Class 2 vs. Class 3 age21+ 1.104 0.867
## Model 1d Class 2 vs. Class 3 genderMale 4.002 3.166
## Model 1d Class 2 vs. Class 3 genderOther 0.395 0.170
## Model 1d Class 2 vs. Class 3 raceAfrican American/Black 0.825 0.500
## Model 1d Class 2 vs. Class 3 raceAsian 0.338 0.204
## Model 1d Class 2 vs. Class 3 raceMultiracial 0.531 0.277
## Model 1d Class 2 vs. Class 3 raceOther 0.623 0.301
## Model 1d Class 2 vs. Class 3 ethnicityHispanic/Latino 0.724 0.493
## Model 1d Class 2 vs. Class 3 studentstatusPart Time 1.004 0.570
## Model 1d Class 2 vs. Class 3 residenceOn Campus 1.196 0.937
## Model 1d Class 2 vs. Class 3 greekmemGreek 1.418 1.043
## Model 1d Class 2 vs. Class 3 alcmo_dich2 1.246 0.905
## Model 1d Class 2 vs. Class 3 marmo_dich2 1.193 0.900
## Model 1d Class 2 vs. Class 3 ecigmo_dich2 1.533 1.132
## Model 1d Class 2 vs. Class 3 smokmo_dich2 1.197 0.876
## Model 1d Class 2 vs. Class 3 cigarmo_dich2 2.425 1.791
## Model 1d Class 2 vs. Class 3 rxstmo_dich2 1.812 1.082
## Model 1d Class 2 vs. Class 3 rxpkmo_dich2 1.401 0.687
## Model 1d Class 2 vs. Class 3 rxsedmo_dich2 0.550 0.239
## CI_upper p_value sig
## 0.401 0.0000 ***
## 1.192 0.8492
## 1.389 0.0495 *
## 1.060 0.1074
## 0.771 0.0008 ***
## 0.562 0.0000 ***
## 1.266 0.5635
## 0.892 0.0144 *
## 1.145 0.4334
## 1.678 0.3572
## 1.319 0.1266
## 1.420 0.2615
## 1.379 0.1714
## 1.209 0.9699
## 1.590 0.0106 *
## 1.312 0.6592
## 1.785 0.0082 **
## 2.259 0.0730 .
## 2.623 0.1059
## 1.193 0.1585
## 0.064 0.0000 ***
## 1.405 0.4244
## 5.059 0.0000 ***
## 0.916 0.0306 *
## 1.362 0.4517
## 0.559 0.0000 ***
## 1.019 0.0569 .
## 1.287 0.2012
## 1.064 0.1004
## 1.769 0.9882
## 1.526 0.1515
## 1.928 0.0258 *
## 1.717 0.1777
## 1.581 0.2199
## 2.075 0.0058 **
## 1.636 0.2599
## 3.282 0.0000 ***
## 3.036 0.0240 *
## 2.857 0.3541
## 1.263 0.1585
print(or_2b, row.names = FALSE)
## model contrast term OR CI_lower
## Model 2b Class 1 vs. Class 3 (Intercept) 0.184 0.141
## Model 2b Class 1 vs. Class 3 age21+ 0.631 0.483
## Model 2b Class 1 vs. Class 3 genderMale 2.554 1.964
## Model 2b Class 1 vs. Class 3 genderOther 4.827 2.920
## Model 2b Class 1 vs. Class 3 raceAfrican American/Black 2.235 1.257
## Model 2b Class 1 vs. Class 3 raceAsian 3.732 2.387
## Model 2b Class 1 vs. Class 3 raceMultiracial 1.436 0.814
## Model 2b Class 1 vs. Class 3 raceOther 2.342 1.119
## Model 2b Class 1 vs. Class 3 ethnicityHispanic/Latino 1.370 0.920
## Model 2b Class 1 vs. Class 3 studentstatusPart Time 1.233 0.674
## Model 2b Class 1 vs. Class 3 residenceOn Campus 1.155 0.880
## Model 2b Class 1 vs. Class 3 greekmemGreek 0.782 0.519
## Model 2b Class 2 vs. Class 3 (Intercept) 0.085 0.062
## Model 2b Class 2 vs. Class 3 age21+ 1.235 0.937
## Model 2b Class 2 vs. Class 3 genderMale 6.382 4.824
## Model 2b Class 2 vs. Class 3 genderOther 1.999 0.900
## Model 2b Class 2 vs. Class 3 raceAfrican American/Black 1.980 1.069
## Model 2b Class 2 vs. Class 3 raceAsian 1.283 0.717
## Model 2b Class 2 vs. Class 3 raceMultiracial 0.566 0.256
## Model 2b Class 2 vs. Class 3 raceOther 1.086 0.401
## Model 2b Class 2 vs. Class 3 ethnicityHispanic/Latino 0.979 0.619
## Model 2b Class 2 vs. Class 3 studentstatusPart Time 0.771 0.378
## Model 2b Class 2 vs. Class 3 residenceOn Campus 1.182 0.882
## Model 2b Class 2 vs. Class 3 greekmemGreek 1.657 1.167
## CI_upper p_value sig
## 0.238 0.0000 ***
## 0.824 0.0008 ***
## 3.321 0.0000 ***
## 7.979 0.0000 ***
## 3.973 0.0062 **
## 5.837 0.0000 ***
## 2.532 0.2116
## 4.904 0.0241 *
## 2.042 0.1216
## 2.255 0.4960
## 1.515 0.2998
## 1.180 0.2420
## 0.117 0.0000 ***
## 1.629 0.1343
## 8.444 0.0000 ***
## 4.440 0.0889 .
## 3.666 0.0299 *
## 2.295 0.4016
## 1.254 0.1613
## 2.941 0.8706
## 1.549 0.9291
## 1.574 0.4750
## 1.584 0.2622
## 2.353 0.0048 **
print(or_2c, row.names = FALSE)
## model contrast term OR CI_lower
## Model 2c Class 1 vs. Class 3 (Intercept) 0.212 0.153
## Model 2c Class 1 vs. Class 3 age21+ 0.710 0.533
## Model 2c Class 1 vs. Class 3 genderMale 2.677 2.037
## Model 2c Class 1 vs. Class 3 genderOther 4.894 2.942
## Model 2c Class 1 vs. Class 3 raceAfrican American/Black 2.049 1.148
## Model 2c Class 1 vs. Class 3 raceAsian 3.308 2.096
## Model 2c Class 1 vs. Class 3 raceMultiracial 1.419 0.802
## Model 2c Class 1 vs. Class 3 raceOther 2.439 1.156
## Model 2c Class 1 vs. Class 3 ethnicityHispanic/Latino 1.314 0.878
## Model 2c Class 1 vs. Class 3 studentstatusPart Time 1.261 0.690
## Model 2c Class 1 vs. Class 3 residenceOn Campus 1.159 0.882
## Model 2c Class 1 vs. Class 3 greekmemGreek 0.893 0.588
## Model 2c Class 1 vs. Class 3 alcmo_dich2 0.882 0.632
## Model 2c Class 1 vs. Class 3 marmo_dich2 0.932 0.669
## Model 2c Class 1 vs. Class 3 ecigmo_dich2 1.054 0.733
## Model 2c Class 1 vs. Class 3 smokmo_dich2 0.763 0.508
## Model 2c Class 1 vs. Class 3 cigarmo_dich2 0.691 0.446
## 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.630 4.186
## Model 2c Class 2 vs. Class 3 genderOther 1.936 0.865
## Model 2c Class 2 vs. Class 3 raceAfrican American/Black 2.248 1.199
## Model 2c Class 2 vs. Class 3 raceAsian 1.629 0.896
## Model 2c Class 2 vs. Class 3 raceMultiracial 0.579 0.260
## Model 2c Class 2 vs. Class 3 raceOther 1.007 0.369
## Model 2c Class 2 vs. Class 3 ethnicityHispanic/Latino 1.018 0.639
## Model 2c Class 2 vs. Class 3 studentstatusPart Time 0.784 0.380
## Model 2c Class 2 vs. Class 3 residenceOn Campus 1.099 0.814
## Model 2c Class 2 vs. Class 3 greekmemGreek 1.368 0.946
## Model 2c Class 2 vs. Class 3 alcmo_dich2 0.838 0.558
## Model 2c Class 2 vs. Class 3 marmo_dich2 1.076 0.761
## Model 2c Class 2 vs. Class 3 ecigmo_dich2 1.520 1.045
## Model 2c Class 2 vs. Class 3 smokmo_dich2 1.053 0.723
## Model 2c Class 2 vs. Class 3 cigarmo_dich2 1.799 1.258
## CI_upper p_value sig
## 0.292 0.0000 ***
## 0.946 0.0194 *
## 3.517 0.0000 ***
## 8.140 0.0000 ***
## 3.657 0.0153 *
## 5.223 0.0000 ***
## 2.509 0.2292
## 5.148 0.0194 *
## 1.966 0.1845
## 2.305 0.4504
## 1.524 0.2886
## 1.356 0.5952
## 1.230 0.4581
## 1.300 0.6801
## 1.515 0.7761
## 1.145 0.1911
## 1.069 0.0970 .
## 0.117 0.0000 ***
## 1.462 0.5964
## 7.573 0.0000 ***
## 4.333 0.1083
## 4.214 0.0116 *
## 2.961 0.1098
## 1.290 0.1814
## 2.750 0.9890
## 1.621 0.9409
## 1.619 0.5110
## 1.484 0.5363
## 1.978 0.0963 .
## 1.260 0.3961
## 1.520 0.6791
## 2.210 0.0286 *
## 1.535 0.7875
## 2.574 0.0013 **
print(or_2d, row.names = FALSE)
## model contrast term OR CI_lower
## Model 2d Class 1 vs. Class 3 (Intercept) 0.212 0.154
## Model 2d Class 1 vs. Class 3 age21+ 0.713 0.535
## Model 2d Class 1 vs. Class 3 genderMale 2.673 2.034
## Model 2d Class 1 vs. Class 3 genderOther 4.814 2.895
## Model 2d Class 1 vs. Class 3 raceAfrican American/Black 2.038 1.142
## Model 2d Class 1 vs. Class 3 raceAsian 3.282 2.071
## Model 2d Class 1 vs. Class 3 raceMultiracial 1.391 0.786
## Model 2d Class 1 vs. Class 3 raceOther 2.405 1.140
## Model 2d Class 1 vs. Class 3 ethnicityHispanic/Latino 1.316 0.879
## Model 2d Class 1 vs. Class 3 studentstatusPart Time 1.238 0.677
## Model 2d Class 1 vs. Class 3 residenceOn Campus 1.153 0.877
## Model 2d Class 1 vs. Class 3 greekmemGreek 0.906 0.596
## Model 2d Class 1 vs. Class 3 alcmo_dich2 0.879 0.629
## Model 2d Class 1 vs. Class 3 marmo_dich2 0.952 0.681
## Model 2d Class 1 vs. Class 3 ecigmo_dich2 1.050 0.730
## Model 2d Class 1 vs. Class 3 smokmo_dich2 0.758 0.504
## Model 2d Class 1 vs. Class 3 cigarmo_dich2 0.706 0.456
## Model 2d Class 1 vs. Class 3 rxstmo_dich2 0.569 0.239
## Model 2d Class 1 vs. Class 3 rxpkmo_dich2 1.084 0.465
## Model 2d Class 1 vs. Class 3 rxsedmo_dich2 2.015 0.681
## Model 2d Class 2 vs. Class 3 (Intercept) 0.078 0.052
## Model 2d Class 2 vs. Class 3 age21+ 1.084 0.803
## Model 2d Class 2 vs. Class 3 genderMale 5.598 4.160
## Model 2d Class 2 vs. Class 3 genderOther 1.948 0.868
## Model 2d Class 2 vs. Class 3 raceAfrican American/Black 2.254 1.203
## Model 2d Class 2 vs. Class 3 raceAsian 1.611 0.886
## Model 2d Class 2 vs. Class 3 raceMultiracial 0.581 0.260
## Model 2d Class 2 vs. Class 3 raceOther 1.017 0.372
## Model 2d Class 2 vs. Class 3 ethnicityHispanic/Latino 1.008 0.631
## Model 2d Class 2 vs. Class 3 studentstatusPart Time 0.796 0.385
## Model 2d Class 2 vs. Class 3 residenceOn Campus 1.096 0.812
## Model 2d Class 2 vs. Class 3 greekmemGreek 1.359 0.938
## Model 2d Class 2 vs. Class 3 alcmo_dich2 0.839 0.558
## Model 2d Class 2 vs. Class 3 marmo_dich2 1.070 0.755
## Model 2d Class 2 vs. Class 3 ecigmo_dich2 1.512 1.040
## Model 2d Class 2 vs. Class 3 smokmo_dich2 1.049 0.718
## Model 2d Class 2 vs. Class 3 cigarmo_dich2 1.797 1.254
## Model 2d Class 2 vs. Class 3 rxstmo_dich2 1.046 0.580
## Model 2d Class 2 vs. Class 3 rxpkmo_dich2 1.265 0.575
## Model 2d Class 2 vs. Class 3 rxsedmo_dich2 0.744 0.267
## CI_upper p_value sig
## 0.293 0.0000 ***
## 0.951 0.0213 *
## 3.512 0.0000 ***
## 8.005 0.0000 ***
## 3.638 0.0161 *
## 5.199 0.0000 ***
## 2.460 0.2576
## 5.070 0.0213 *
## 1.970 0.1827
## 2.263 0.4878
## 1.516 0.3066
## 1.378 0.6459
## 1.226 0.4464
## 1.330 0.7718
## 1.511 0.7912
## 1.141 0.1845
## 1.093 0.1186
## 1.354 0.2023
## 2.531 0.8512
## 5.963 0.2056
## 0.117 0.0000 ***
## 1.462 0.5990
## 7.533 0.0000 ***
## 4.368 0.1060
## 4.222 0.0112 *
## 2.928 0.1184
## 1.295 0.1843
## 2.775 0.9743
## 1.608 0.9746
## 1.644 0.5371
## 1.480 0.5503
## 1.969 0.1053
## 1.261 0.3980
## 1.517 0.7031
## 2.200 0.0307 *
## 1.534 0.8047
## 2.576 0.0014 **
## 1.887 0.8815
## 2.782 0.5596
## 2.071 0.5712
- 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,
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,
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 -3578.64 0.049812 0.0840 0.1014
## Model 1c -3766.25 -3500.99 0.070429 0.1166 0.1408
## Model 1d -3766.25 -3495.38 0.071920 0.1189 0.1436
## Model 2b -1683.88 -1520.59 0.096972 0.1646 0.1951
## Model 2c -1683.88 -1490.60 0.114782 0.1917 0.2273
## Model 2d -1683.88 -1488.67 0.115929 0.1934 0.2293
# 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.0498 0.1014
## 1c - 1b 0.0206 0.0394
## 1d - 1c 0.0015 0.0028
## 2b - 2a 0.0970 0.1951
## 2c - 2b 0.0178 0.0322
## 2d - 2c 0.0011 0.0020
- 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.462]
## 2 1b Class 1 vs. Class 3 age21+ 1.143 [0.984, 1.326]
## 3 1b Class 1 vs. Class 3 genderMale 1.208 [1.032, 1.413]
## 4 1b Class 1 vs. Class 3 genderOther 0.780 [0.566, 1.074]
## 5 1b Class 1 vs. Class 3 raceAfrican American/Black 0.496 [0.347, 0.711]
## 6 1b Class 1 vs. Class 3 raceAsian 0.370 [0.275, 0.497]
## 7 1b Class 1 vs. Class 3 raceMultiracial 0.897 [0.643, 1.251]
## 8 1b Class 1 vs. Class 3 raceOther 0.544 [0.345, 0.857]
## 9 1b Class 1 vs. Class 3 ethnicityHispanic/Latino 0.915 [0.731, 1.146]
## 10 1b Class 1 vs. Class 3 studentstatusPart Time 1.163 [0.819, 1.649]
## 11 1b Class 1 vs. Class 3 residenceOn Campus 1.139 [0.977, 1.329]
## 12 1b Class 1 vs. Class 3 greekmemGreek 1.257 [1.010, 1.563]
## 13 1b Class 2 vs. Class 3 (Intercept) 0.070 [0.055, 0.089]
## 14 1b Class 2 vs. Class 3 age21+ 1.501 [1.201, 1.876]
## 15 1b Class 2 vs. Class 3 genderMale 4.652 [3.743, 5.781]
## 16 1b Class 2 vs. Class 3 genderOther 0.421 [0.183, 0.970]
## 17 1b Class 2 vs. Class 3 raceAfrican American/Black 0.677 [0.416, 1.102]
## 18 1b Class 2 vs. Class 3 raceAsian 0.238 [0.146, 0.389]
## 19 1b Class 2 vs. Class 3 raceMultiracial 0.495 [0.260, 0.943]
## 20 1b Class 2 vs. Class 3 raceOther 0.586 [0.290, 1.187]
## 21 1b Class 2 vs. Class 3 ethnicityHispanic/Latino 0.709 [0.487, 1.032]
## 22 1b Class 2 vs. Class 3 studentstatusPart Time 0.951 [0.547, 1.654]
## 23 1b Class 2 vs. Class 3 residenceOn Campus 1.251 [0.989, 1.583]
## 24 1b Class 2 vs. Class 3 greekmemGreek 2.006 [1.506, 2.672]
## 25 1c Class 1 vs. Class 3 (Intercept) 0.335 [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.182 [1.004, 1.393]
## 28 1c Class 1 vs. Class 3 genderOther 0.774 [0.560, 1.068]
## 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.419 [0.310, 0.566]
## 31 1c Class 1 vs. Class 3 raceMultiracial 0.900 [0.644, 1.258]
## 32 1c Class 1 vs. Class 3 raceOther 0.554 [0.350, 0.877]
## 33 1c Class 1 vs. Class 3 ethnicityHispanic/Latino 0.922 [0.736, 1.156]
## 34 1c Class 1 vs. Class 3 studentstatusPart Time 1.158 [0.815, 1.646]
## 35 1c Class 1 vs. Class 3 residenceOn Campus 1.129 [0.966, 1.319]
## 36 1c Class 1 vs. Class 3 greekmemGreek 1.145 [0.916, 1.430]
## 37 1c Class 1 vs. Class 3 alcmo_dich2 1.142 [0.945, 1.380]
## 38 1c Class 1 vs. Class 3 marmo_dich2 1.022 [0.850, 1.230]
## 39 1c Class 1 vs. Class 3 ecigmo_dich2 1.308 [1.070, 1.599]
## 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.420 [1.112, 1.814]
## 42 1c Class 2 vs. Class 3 (Intercept) 0.046 [0.034, 0.064]
## 43 1c Class 2 vs. Class 3 age21+ 1.116 [0.877, 1.421]
## 44 1c Class 2 vs. Class 3 genderMale 4.030 [3.190, 5.093]
## 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.813 [0.492, 1.341]
## 47 1c Class 2 vs. Class 3 raceAsian 0.342 [0.207, 0.566]
## 48 1c Class 2 vs. Class 3 raceMultiracial 0.523 [0.273, 1.004]
## 49 1c Class 2 vs. Class 3 raceOther 0.594 [0.287, 1.230]
## 50 1c Class 2 vs. Class 3 ethnicityHispanic/Latino 0.730 [0.497, 1.072]
## 51 1c Class 2 vs. Class 3 studentstatusPart Time 0.970 [0.551, 1.708]
## 52 1c Class 2 vs. Class 3 residenceOn Campus 1.195 [0.937, 1.525]
## 53 1c Class 2 vs. Class 3 greekmemGreek 1.452 [1.071, 1.970]
## 54 1c Class 2 vs. Class 3 alcmo_dich2 1.237 [0.898, 1.703]
## 55 1c Class 2 vs. Class 3 marmo_dich2 1.234 [0.934, 1.631]
## 56 1c Class 2 vs. Class 3 ecigmo_dich2 1.550 [1.146, 2.096]
## 57 1c Class 2 vs. Class 3 smokmo_dich2 1.227 [0.900, 1.672]
## 58 1c Class 2 vs. Class 3 cigarmo_dich2 2.478 [1.835, 3.348]
## 59 1d Class 1 vs. Class 3 (Intercept) 0.336 [0.282, 0.401]
## 60 1d Class 1 vs. Class 3 age21+ 1.016 [0.866, 1.192]
## 61 1d Class 1 vs. Class 3 genderMale 1.179 [1.000, 1.389]
## 62 1d Class 1 vs. Class 3 genderOther 0.766 [0.554, 1.060]
## 63 1d Class 1 vs. Class 3 raceAfrican American/Black 0.537 [0.374, 0.771]
## 64 1d Class 1 vs. Class 3 raceAsian 0.416 [0.308, 0.562]
## 65 1d Class 1 vs. Class 3 raceMultiracial 0.906 [0.648, 1.266]
## 66 1d Class 1 vs. Class 3 raceOther 0.564 [0.356, 0.892]
## 67 1d Class 1 vs. Class 3 ethnicityHispanic/Latino 0.914 [0.729, 1.145]
## 68 1d Class 1 vs. Class 3 studentstatusPart Time 1.180 [0.830, 1.678]
## 69 1d Class 1 vs. Class 3 residenceOn Campus 1.129 [0.966, 1.319]
## 70 1d Class 1 vs. Class 3 greekmemGreek 1.136 [0.909, 1.420]
## 71 1d Class 1 vs. Class 3 alcmo_dich2 1.141 [0.944, 1.379]
## 72 1d Class 1 vs. Class 3 marmo_dich2 1.004 [0.833, 1.209]
## 73 1d Class 1 vs. Class 3 ecigmo_dich2 1.300 [1.063, 1.590]
## 74 1d Class 1 vs. Class 3 smokmo_dich2 1.051 [0.842, 1.312]
## 75 1d Class 1 vs. Class 3 cigarmo_dich2 1.395 [1.090, 1.785]
## 76 1d Class 1 vs. Class 3 rxstmo_dich2 1.476 [0.964, 2.259]
## 77 1d Class 1 vs. Class 3 rxpkmo_dich2 1.546 [0.912, 2.623]
## 78 1d Class 1 vs. Class 3 rxsedmo_dich2 0.636 [0.339, 1.193]
## 79 1d Class 2 vs. Class 3 (Intercept) 0.047 [0.034, 0.064]
## 80 1d Class 2 vs. Class 3 age21+ 1.104 [0.867, 1.405]
## 81 1d Class 2 vs. Class 3 genderMale 4.002 [3.166, 5.059]
## 82 1d Class 2 vs. Class 3 genderOther 0.395 [0.170, 0.916]
## 83 1d Class 2 vs. Class 3 raceAfrican American/Black 0.825 [0.500, 1.362]
## 84 1d Class 2 vs. Class 3 raceAsian 0.338 [0.204, 0.559]
## 85 1d Class 2 vs. Class 3 raceMultiracial 0.531 [0.277, 1.019]
## 86 1d Class 2 vs. Class 3 raceOther 0.623 [0.301, 1.287]
## 87 1d Class 2 vs. Class 3 ethnicityHispanic/Latino 0.724 [0.493, 1.064]
## 88 1d Class 2 vs. Class 3 studentstatusPart Time 1.004 [0.570, 1.769]
## 89 1d Class 2 vs. Class 3 residenceOn Campus 1.196 [0.937, 1.526]
## 90 1d Class 2 vs. Class 3 greekmemGreek 1.418 [1.043, 1.928]
## 91 1d Class 2 vs. Class 3 alcmo_dich2 1.246 [0.905, 1.717]
## 92 1d Class 2 vs. Class 3 marmo_dich2 1.193 [0.900, 1.581]
## 93 1d Class 2 vs. Class 3 ecigmo_dich2 1.533 [1.132, 2.075]
## 94 1d Class 2 vs. Class 3 smokmo_dich2 1.197 [0.876, 1.636]
## 95 1d Class 2 vs. Class 3 cigarmo_dich2 2.425 [1.791, 3.282]
## 96 1d Class 2 vs. Class 3 rxstmo_dich2 1.812 [1.082, 3.036]
## 97 1d Class 2 vs. Class 3 rxpkmo_dich2 1.401 [0.687, 2.857]
## 98 1d Class 2 vs. Class 3 rxsedmo_dich2 0.550 [0.239, 1.263]
## 99 2b Class 1 vs. Class 3 (Intercept) 0.184 [0.141, 0.238]
## 100 2b Class 1 vs. Class 3 age21+ 0.631 [0.483, 0.824]
## 101 2b Class 1 vs. Class 3 genderMale 2.554 [1.964, 3.321]
## 102 2b Class 1 vs. Class 3 genderOther 4.827 [2.920, 7.979]
## 103 2b Class 1 vs. Class 3 raceAfrican American/Black 2.235 [1.257, 3.973]
## 104 2b Class 1 vs. Class 3 raceAsian 3.732 [2.387, 5.837]
## 105 2b Class 1 vs. Class 3 raceMultiracial 1.436 [0.814, 2.532]
## 106 2b Class 1 vs. Class 3 raceOther 2.342 [1.119, 4.904]
## 107 2b Class 1 vs. Class 3 ethnicityHispanic/Latino 1.370 [0.920, 2.042]
## 108 2b Class 1 vs. Class 3 studentstatusPart Time 1.233 [0.674, 2.255]
## 109 2b Class 1 vs. Class 3 residenceOn Campus 1.155 [0.880, 1.515]
## 110 2b Class 1 vs. Class 3 greekmemGreek 0.782 [0.519, 1.180]
## 111 2b Class 2 vs. Class 3 (Intercept) 0.085 [0.062, 0.117]
## 112 2b Class 2 vs. Class 3 age21+ 1.235 [0.937, 1.629]
## 113 2b Class 2 vs. Class 3 genderMale 6.382 [4.824, 8.444]
## 114 2b Class 2 vs. Class 3 genderOther 1.999 [0.900, 4.440]
## 115 2b Class 2 vs. Class 3 raceAfrican American/Black 1.980 [1.069, 3.666]
## 116 2b Class 2 vs. Class 3 raceAsian 1.283 [0.717, 2.295]
## 117 2b Class 2 vs. Class 3 raceMultiracial 0.566 [0.256, 1.254]
## 118 2b Class 2 vs. Class 3 raceOther 1.086 [0.401, 2.941]
## 119 2b Class 2 vs. Class 3 ethnicityHispanic/Latino 0.979 [0.619, 1.549]
## 120 2b Class 2 vs. Class 3 studentstatusPart Time 0.771 [0.378, 1.574]
## 121 2b Class 2 vs. Class 3 residenceOn Campus 1.182 [0.882, 1.584]
## 122 2b Class 2 vs. Class 3 greekmemGreek 1.657 [1.167, 2.353]
## 123 2c Class 1 vs. Class 3 (Intercept) 0.212 [0.153, 0.292]
## 124 2c Class 1 vs. Class 3 age21+ 0.710 [0.533, 0.946]
## 125 2c Class 1 vs. Class 3 genderMale 2.677 [2.037, 3.517]
## 126 2c Class 1 vs. Class 3 genderOther 4.894 [2.942, 8.140]
## 127 2c Class 1 vs. Class 3 raceAfrican American/Black 2.049 [1.148, 3.657]
## 128 2c Class 1 vs. Class 3 raceAsian 3.308 [2.096, 5.223]
## 129 2c Class 1 vs. Class 3 raceMultiracial 1.419 [0.802, 2.509]
## 130 2c Class 1 vs. Class 3 raceOther 2.439 [1.156, 5.148]
## 131 2c Class 1 vs. Class 3 ethnicityHispanic/Latino 1.314 [0.878, 1.966]
## 132 2c Class 1 vs. Class 3 studentstatusPart Time 1.261 [0.690, 2.305]
## 133 2c Class 1 vs. Class 3 residenceOn Campus 1.159 [0.882, 1.524]
## 134 2c Class 1 vs. Class 3 greekmemGreek 0.893 [0.588, 1.356]
## 135 2c Class 1 vs. Class 3 alcmo_dich2 0.882 [0.632, 1.230]
## 136 2c Class 1 vs. Class 3 marmo_dich2 0.932 [0.669, 1.300]
## 137 2c Class 1 vs. Class 3 ecigmo_dich2 1.054 [0.733, 1.515]
## 138 2c Class 1 vs. Class 3 smokmo_dich2 0.763 [0.508, 1.145]
## 139 2c Class 1 vs. Class 3 cigarmo_dich2 0.691 [0.446, 1.069]
## 140 2c Class 2 vs. Class 3 (Intercept) 0.078 [0.052, 0.117]
## 141 2c Class 2 vs. Class 3 age21+ 1.084 [0.804, 1.462]
## 142 2c Class 2 vs. Class 3 genderMale 5.630 [4.186, 7.573]
## 143 2c Class 2 vs. Class 3 genderOther 1.936 [0.865, 4.333]
## 144 2c Class 2 vs. Class 3 raceAfrican American/Black 2.248 [1.199, 4.214]
## 145 2c Class 2 vs. Class 3 raceAsian 1.629 [0.896, 2.961]
## 146 2c Class 2 vs. Class 3 raceMultiracial 0.579 [0.260, 1.290]
## 147 2c Class 2 vs. Class 3 raceOther 1.007 [0.369, 2.750]
## 148 2c Class 2 vs. Class 3 ethnicityHispanic/Latino 1.018 [0.639, 1.621]
## 149 2c Class 2 vs. Class 3 studentstatusPart Time 0.784 [0.380, 1.619]
## 150 2c Class 2 vs. Class 3 residenceOn Campus 1.099 [0.814, 1.484]
## 151 2c Class 2 vs. Class 3 greekmemGreek 1.368 [0.946, 1.978]
## 152 2c Class 2 vs. Class 3 alcmo_dich2 0.838 [0.558, 1.260]
## 153 2c Class 2 vs. Class 3 marmo_dich2 1.076 [0.761, 1.520]
## 154 2c Class 2 vs. Class 3 ecigmo_dich2 1.520 [1.045, 2.210]
## 155 2c Class 2 vs. Class 3 smokmo_dich2 1.053 [0.723, 1.535]
## 156 2c Class 2 vs. Class 3 cigarmo_dich2 1.799 [1.258, 2.574]
## 157 2d Class 1 vs. Class 3 (Intercept) 0.212 [0.154, 0.293]
## 158 2d Class 1 vs. Class 3 age21+ 0.713 [0.535, 0.951]
## 159 2d Class 1 vs. Class 3 genderMale 2.673 [2.034, 3.512]
## 160 2d Class 1 vs. Class 3 genderOther 4.814 [2.895, 8.005]
## 161 2d Class 1 vs. Class 3 raceAfrican American/Black 2.038 [1.142, 3.638]
## 162 2d Class 1 vs. Class 3 raceAsian 3.282 [2.071, 5.199]
## 163 2d Class 1 vs. Class 3 raceMultiracial 1.391 [0.786, 2.460]
## 164 2d Class 1 vs. Class 3 raceOther 2.405 [1.140, 5.070]
## 165 2d Class 1 vs. Class 3 ethnicityHispanic/Latino 1.316 [0.879, 1.970]
## 166 2d Class 1 vs. Class 3 studentstatusPart Time 1.238 [0.677, 2.263]
## 167 2d Class 1 vs. Class 3 residenceOn Campus 1.153 [0.877, 1.516]
## 168 2d Class 1 vs. Class 3 greekmemGreek 0.906 [0.596, 1.378]
## 169 2d Class 1 vs. Class 3 alcmo_dich2 0.879 [0.629, 1.226]
## 170 2d Class 1 vs. Class 3 marmo_dich2 0.952 [0.681, 1.330]
## 171 2d Class 1 vs. Class 3 ecigmo_dich2 1.050 [0.730, 1.511]
## 172 2d Class 1 vs. Class 3 smokmo_dich2 0.758 [0.504, 1.141]
## 173 2d Class 1 vs. Class 3 cigarmo_dich2 0.706 [0.456, 1.093]
## 174 2d Class 1 vs. Class 3 rxstmo_dich2 0.569 [0.239, 1.354]
## 175 2d Class 1 vs. Class 3 rxpkmo_dich2 1.084 [0.465, 2.531]
## 176 2d Class 1 vs. Class 3 rxsedmo_dich2 2.015 [0.681, 5.963]
## 177 2d Class 2 vs. Class 3 (Intercept) 0.078 [0.052, 0.117]
## 178 2d Class 2 vs. Class 3 age21+ 1.084 [0.803, 1.462]
## 179 2d Class 2 vs. Class 3 genderMale 5.598 [4.160, 7.533]
## 180 2d Class 2 vs. Class 3 genderOther 1.948 [0.868, 4.368]
## 181 2d Class 2 vs. Class 3 raceAfrican American/Black 2.254 [1.203, 4.222]
## 182 2d Class 2 vs. Class 3 raceAsian 1.611 [0.886, 2.928]
## 183 2d Class 2 vs. Class 3 raceMultiracial 0.581 [0.260, 1.295]
## 184 2d Class 2 vs. Class 3 raceOther 1.017 [0.372, 2.775]
## 185 2d Class 2 vs. Class 3 ethnicityHispanic/Latino 1.008 [0.631, 1.608]
## 186 2d Class 2 vs. Class 3 studentstatusPart Time 0.796 [0.385, 1.644]
## 187 2d Class 2 vs. Class 3 residenceOn Campus 1.096 [0.812, 1.480]
## 188 2d Class 2 vs. Class 3 greekmemGreek 1.359 [0.938, 1.969]
## 189 2d Class 2 vs. Class 3 alcmo_dich2 0.839 [0.558, 1.261]
## 190 2d Class 2 vs. Class 3 marmo_dich2 1.070 [0.755, 1.517]
## 191 2d Class 2 vs. Class 3 ecigmo_dich2 1.512 [1.040, 2.200]
## 192 2d Class 2 vs. Class 3 smokmo_dich2 1.049 [0.718, 1.534]
## 193 2d Class 2 vs. Class 3 cigarmo_dich2 1.797 [1.254, 2.576]
## 194 2d Class 2 vs. Class 3 rxstmo_dich2 1.046 [0.580, 1.887]
## 195 2d Class 2 vs. Class 3 rxpkmo_dich2 1.265 [0.575, 2.782]
## 196 2d Class 2 vs. Class 3 rxsedmo_dich2 0.744 [0.267, 2.071]
## p.value sig fmi McFadden Nagelkerke
## 1 0.0000 *** 0.004 0.050 0.101
## 2 0.0799 . 0.001 0.050 0.101
## 3 0.0184 * 0.004 0.050 0.101
## 4 0.1279 0.006 0.050 0.101
## 5 0.0001 *** 0.007 0.050 0.101
## 6 0.0000 *** 0.012 0.050 0.101
## 7 0.5221 0.005 0.050 0.101
## 8 0.0087 ** 0.009 0.050 0.101
## 9 0.4393 0.028 0.050 0.101
## 10 0.3986 0.003 0.050 0.101
## 11 0.0962 . 0.002 0.050 0.101
## 12 0.0401 * 0.002 0.050 0.101
## 13 0.0000 *** 0.002 0.050 0.101
## 14 0.0004 *** 0.001 0.050 0.101
## 15 0.0000 *** 0.002 0.050 0.101
## 16 0.0423 * 0.001 0.050 0.101
## 17 0.1163 0.006 0.050 0.101
## 18 0.0000 *** 0.008 0.050 0.101
## 19 0.0325 * 0.003 0.050 0.101
## 20 0.1380 0.006 0.050 0.101
## 21 0.0729 . 0.068 0.050 0.101
## 22 0.8579 0.008 0.050 0.101
## 23 0.0623 . 0.001 0.050 0.101
## 24 0.0000 *** 0.001 0.050 0.101
## 25 0.0000 *** 0.005 0.070 0.141
## 26 0.8055 0.001 0.070 0.141
## 27 0.0450 * 0.005 0.070 0.141
## 28 0.1191 0.007 0.070 0.141
## 29 0.0006 *** 0.008 0.070 0.141
## 30 0.0000 *** 0.013 0.070 0.141
## 31 0.5374 0.005 0.070 0.141
## 32 0.0117 * 0.009 0.070 0.141
## 33 0.4822 0.024 0.070 0.141
## 34 0.4127 0.003 0.070 0.141
## 35 0.1263 0.002 0.070 0.141
## 36 0.2334 0.002 0.070 0.141
## 37 0.1683 0.007 0.070 0.141
## 38 0.8161 0.003 0.070 0.141
## 39 0.0088 ** 0.004 0.070 0.141
## 40 0.5540 0.002 0.070 0.141
## 41 0.0050 ** 0.003 0.070 0.141
## 42 0.0000 *** 0.004 0.070 0.141
## 43 0.3702 0.001 0.070 0.141
## 44 0.0000 *** 0.003 0.070 0.141
## 45 0.0328 * 0.001 0.070 0.141
## 46 0.4170 0.006 0.070 0.141
## 47 0.0000 *** 0.008 0.070 0.141
## 48 0.0514 . 0.003 0.070 0.141
## 49 0.1610 0.006 0.070 0.141
## 50 0.1086 0.061 0.070 0.141
## 51 0.9149 0.009 0.070 0.141
## 52 0.1517 0.002 0.070 0.141
## 53 0.0164 * 0.001 0.070 0.141
## 54 0.1934 0.006 0.070 0.141
## 55 0.1387 0.002 0.070 0.141
## 56 0.0045 ** 0.001 0.070 0.141
## 57 0.1956 0.002 0.070 0.141
## 58 0.0000 *** 0.002 0.070 0.141
## 59 0.0000 *** 0.005 0.072 0.144
## 60 0.8492 0.001 0.072 0.144
## 61 0.0495 * 0.004 0.072 0.144
## 62 0.1074 0.007 0.072 0.144
## 63 0.0008 *** 0.008 0.072 0.144
## 64 0.0000 *** 0.013 0.072 0.144
## 65 0.5635 0.005 0.072 0.144
## 66 0.0144 * 0.009 0.072 0.144
## 67 0.4334 0.024 0.072 0.144
## 68 0.3572 0.003 0.072 0.144
## 69 0.1266 0.002 0.072 0.144
## 70 0.2615 0.002 0.072 0.144
## 71 0.1714 0.007 0.072 0.144
## 72 0.9699 0.003 0.072 0.144
## 73 0.0106 * 0.004 0.072 0.144
## 74 0.6592 0.002 0.072 0.144
## 75 0.0082 ** 0.004 0.072 0.144
## 76 0.0730 . 0.002 0.072 0.144
## 77 0.1059 0.001 0.072 0.144
## 78 0.1585 0.001 0.072 0.144
## 79 0.0000 *** 0.004 0.072 0.144
## 80 0.4244 0.001 0.072 0.144
## 81 0.0000 *** 0.003 0.072 0.144
## 82 0.0306 * 0.001 0.072 0.144
## 83 0.4517 0.006 0.072 0.144
## 84 0.0000 *** 0.008 0.072 0.144
## 85 0.0569 . 0.003 0.072 0.144
## 86 0.2012 0.006 0.072 0.144
## 87 0.1004 0.063 0.072 0.144
## 88 0.9882 0.010 0.072 0.144
## 89 0.1515 0.002 0.072 0.144
## 90 0.0258 * 0.001 0.072 0.144
## 91 0.1777 0.006 0.072 0.144
## 92 0.2199 0.001 0.072 0.144
## 93 0.0058 ** 0.001 0.072 0.144
## 94 0.2599 0.002 0.072 0.144
## 95 0.0000 *** 0.002 0.072 0.144
## 96 0.0240 * 0.004 0.072 0.144
## 97 0.3541 0.009 0.072 0.144
## 98 0.1585 0.002 0.072 0.144
## 99 0.0000 *** 0.004 0.097 0.195
## 100 0.0008 *** 0.001 0.097 0.195
## 101 0.0000 *** 0.006 0.097 0.195
## 102 0.0000 *** 0.003 0.097 0.195
## 103 0.0062 ** 0.006 0.097 0.195
## 104 0.0000 *** 0.007 0.097 0.195
## 105 0.2116 0.004 0.097 0.195
## 106 0.0241 * 0.011 0.097 0.195
## 107 0.1216 0.066 0.097 0.195
## 108 0.4960 0.001 0.097 0.195
## 109 0.2998 0.003 0.097 0.195
## 110 0.2420 0.003 0.097 0.195
## 111 0.0000 *** 0.003 0.097 0.195
## 112 0.1343 0.001 0.097 0.195
## 113 0.0000 *** 0.001 0.097 0.195
## 114 0.0889 . 0.001 0.097 0.195
## 115 0.0299 * 0.004 0.097 0.195
## 116 0.4016 0.009 0.097 0.195
## 117 0.1613 0.002 0.097 0.195
## 118 0.8706 0.005 0.097 0.195
## 119 0.9291 0.072 0.097 0.195
## 120 0.4750 0.006 0.097 0.195
## 121 0.2622 0.002 0.097 0.195
## 122 0.0048 ** 0.001 0.097 0.195
## 123 0.0000 *** 0.006 0.115 0.227
## 124 0.0194 * 0.003 0.115 0.227
## 125 0.0000 *** 0.008 0.115 0.227
## 126 0.0000 *** 0.004 0.115 0.227
## 127 0.0153 * 0.007 0.115 0.227
## 128 0.0000 *** 0.007 0.115 0.227
## 129 0.2292 0.004 0.115 0.227
## 130 0.0194 * 0.011 0.115 0.227
## 131 0.1845 0.067 0.115 0.227
## 132 0.4504 0.001 0.115 0.227
## 133 0.2886 0.003 0.115 0.227
## 134 0.5952 0.003 0.115 0.227
## 135 0.4581 0.011 0.115 0.227
## 136 0.6801 0.003 0.115 0.227
## 137 0.7761 0.010 0.115 0.227
## 138 0.1911 0.003 0.115 0.227
## 139 0.0970 . 0.006 0.115 0.227
## 140 0.0000 *** 0.007 0.115 0.227
## 141 0.5964 0.002 0.115 0.227
## 142 0.0000 *** 0.002 0.115 0.227
## 143 0.1083 0.001 0.115 0.227
## 144 0.0116 * 0.004 0.115 0.227
## 145 0.1098 0.010 0.115 0.227
## 146 0.1814 0.002 0.115 0.227
## 147 0.9890 0.006 0.115 0.227
## 148 0.9409 0.078 0.115 0.227
## 149 0.5110 0.006 0.115 0.227
## 150 0.5363 0.002 0.115 0.227
## 151 0.0963 . 0.002 0.115 0.227
## 152 0.3961 0.007 0.115 0.227
## 153 0.6791 0.002 0.115 0.227
## 154 0.0286 * 0.003 0.115 0.227
## 155 0.7875 0.002 0.115 0.227
## 156 0.0013 ** 0.003 0.115 0.227
## 157 0.0000 *** 0.006 0.116 0.229
## 158 0.0213 * 0.003 0.116 0.229
## 159 0.0000 *** 0.008 0.116 0.229
## 160 0.0000 *** 0.004 0.116 0.229
## 161 0.0161 * 0.007 0.116 0.229
## 162 0.0000 *** 0.008 0.116 0.229
## 163 0.2576 0.004 0.116 0.229
## 164 0.0213 * 0.011 0.116 0.229
## 165 0.1827 0.067 0.116 0.229
## 166 0.4878 0.001 0.116 0.229
## 167 0.3066 0.003 0.116 0.229
## 168 0.6459 0.003 0.116 0.229
## 169 0.4464 0.011 0.116 0.229
## 170 0.7718 0.003 0.116 0.229
## 171 0.7912 0.010 0.116 0.229
## 172 0.1845 0.003 0.116 0.229
## 173 0.1186 0.006 0.116 0.229
## 174 0.2023 0.002 0.116 0.229
## 175 0.8512 0.007 0.116 0.229
## 176 0.2056 0.002 0.116 0.229
## 177 0.0000 *** 0.006 0.116 0.229
## 178 0.5990 0.002 0.116 0.229
## 179 0.0000 *** 0.002 0.116 0.229
## 180 0.1060 0.001 0.116 0.229
## 181 0.0112 * 0.004 0.116 0.229
## 182 0.1184 0.010 0.116 0.229
## 183 0.1843 0.003 0.116 0.229
## 184 0.9743 0.006 0.116 0.229
## 185 0.9746 0.080 0.116 0.229
## 186 0.5371 0.006 0.116 0.229
## 187 0.5503 0.002 0.116 0.229
## 188 0.1053 0.001 0.116 0.229
## 189 0.3980 0.007 0.116 0.229
## 190 0.7031 0.002 0.116 0.229
## 191 0.0307 * 0.003 0.116 0.229
## 192 0.8047 0.003 0.116 0.229
## 193 0.0014 ** 0.003 0.116 0.229
## 194 0.8815 0.005 0.116 0.229
## 195 0.5596 0.036 0.116 0.229
## 196 0.5712 0.007 0.116 0.229