Medication Access

Load Packages

library(tidyverse)
library(codebookr)
library(summarytools)
library(broom)
library(performance)
library(gt)
library(gtsummary)
library(janitor)
library(forcats)
library(here)
library(margins)

Import Data

load(file = "~/Desktop/R-Code/SDOH_ALL/med_access_all.rda")
View(med_access_all)

Data Cleaning

Select variables

med_access_all %>%  
  dplyr::select(sbj_id:lang_3, RPL_THEMES:RPL_4, max_ch, steroid_dep, MESALAMINE:OZANIMOD) -> med_rx 

View(med_rx)

Drop NA values for meds

med_rx %>% 
  na.omit(MESALAMINE:OZANIMOD) -> med_rx_drop

Make med access variable yes/no

med_rx_drop %>%  
mutate(mesalamine_3 = case_when(MESALAMINE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(prednisone_3 = case_when(PREDNISONE>= 3 ~ '1',TRUE ~ "0")) %>%  
mutate(infliximab_3 = case_when(INFLIXIMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(vedolizumab_3 = case_when(VEDOLIZUMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(adalimumab_3 = case_when(ADALIMUMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(tofacitinib_3 = case_when(TOFACITINIB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(ustekinumab_3 = case_when(USTEKINUMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(methylpred_3 = case_when(METHYLPREDNISOLONE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(dexamethasone_3 = case_when(DEXAMETHASONE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(budesonide_3 = case_when(BUDESONIDE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(methotrexate_3 = case_when(METHOTREXATE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(azathioprine_3 = case_when(AZATHIOPRINE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(mercaptopurine_3 = case_when(MERCAPTOPURINE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(sulfasalazine_3 = case_when(SULFASALAZINE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(balsalazide_3 = case_when(BALSALAZIDE>= 3 ~ '1',TRUE ~ "0")) %>%   
mutate(cyclosporine_3 = case_when(CYCLOSPORINE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(certolizumab_3 = case_when(CERTOLIZUMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(prednisolone_3 = case_when(PREDNISOLONE>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(upadacitinib_3 = case_when(UPADACITINIB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(golimumab_3 = case_when(GOLIMUMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(natalizumab_3 = case_when(NATALIZUMAB>= 3 ~ '1',TRUE ~ "0")) %>% 
mutate(ozanimod_3 = case_when(OZANIMOD>= 3 ~ '1',TRUE ~ "0")) -> med_rx_3
View(med_rx_3)

Make numeric

med_rx_3$mesalamine_3 = as.numeric(med_rx_3$mesalamine_3)
med_rx_3$prednisone_3 = as.numeric(med_rx_3$prednisone_3)
med_rx_3$infliximab_3 = as.numeric(med_rx_3$infliximab_3)
med_rx_3$vedolizumab_3 = as.numeric(med_rx_3$vedolizumab_3)
med_rx_3$adalimumab_3 = as.numeric(med_rx_3$adalimumab_3)
med_rx_3$tofacitinib_3 = as.numeric(med_rx_3$tofacitinib_3)
med_rx_3$ustekinumab_3 = as.numeric(med_rx_3$ustekinumab_3)
med_rx_3$methylpred_3 = as.numeric(med_rx_3$methylpred_3)
med_rx_3$dexamethasone_3 = as.numeric(med_rx_3$dexamethasone_3)
med_rx_3$budesonide_3 = as.numeric(med_rx_3$budesonide_3)
med_rx_3$methotrexate_3 = as.numeric(med_rx_3$methotrexate_3)
med_rx_3$azathioprine_3 = as.numeric(med_rx_3$azathioprine_3)
med_rx_3$mercaptopurine_3 = as.numeric(med_rx_3$mercaptopurine_3)
med_rx_3$sulfasalazine_3 = as.numeric(med_rx_3$sulfasalazine_3)
med_rx_3$balsalazide_3 = as.numeric(med_rx_3$balsalazide_3)
med_rx_3$cyclosporine_3 = as.numeric(med_rx_3$cyclosporine_3)
med_rx_3$certolizumab_3 = as.numeric(med_rx_3$certolizumab_3)
med_rx_3$prednisolone_3 = as.numeric(med_rx_3$prednisolone_3)
med_rx_3$upadacitinib_3 = as.numeric(med_rx_3$upadacitinib_3)
med_rx_3$golimumab_3 = as.numeric(med_rx_3$golimumab_3)
med_rx_3$natalizumab_3 = as.numeric(med_rx_3$natalizumab_3)
med_rx_3$ozanimod_3 = as.numeric(med_rx_3$ozanimod_3) 

Meds by class

med_rx_3 %>% 
  mutate(ASA = mesalamine_3 + sulfasalazine_3 + balsalazide_3) %>% 
  mutate(biologic = infliximab_3 + vedolizumab_3 + adalimumab_3 + ustekinumab_3 + certolizumab_3 +
           golimumab_3 + natalizumab_3) %>% 
  mutate(steroids = prednisone_3 + prednisolone_3 + methylpred_3 + dexamethasone_3 + budesonide_3) %>% 
  mutate(immuno = methotrexate_3 + azathioprine_3 + mercaptopurine_3) %>% 
mutate(small = tofacitinib_3 + upadacitinib_3 + ozanimod_3) %>% 

mutate(steroid_count = PREDNISONE + PREDNISOLONE + METHYLPREDNISOLONE + DEXAMETHASONE + BUDESONIDE) -> rx_class_final


View(rx_class_final)

Med access yes/no by class

rx_class_final %>%  
mutate(ASA_3 = case_when(ASA>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(immuno_3 = case_when(immuno>= 1 ~ '1',TRUE ~ "0")) %>%  
mutate(biologic_3 = case_when(biologic>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(small_3 = case_when(small>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(steroids_3 = case_when(steroids>= 1 ~ '1',TRUE ~ "0")) -> access_class_dich

Make numeric

access_class_dich$ASA_3 = as.numeric(access_class_dich$ASA_3)
access_class_dich$immuno_3 = as.numeric(access_class_dich$immuno_3)
access_class_dich$biologic_3 = as.numeric(access_class_dich$biologic_3)
access_class_dich$small_3 = as.numeric(access_class_dich$small_3)
access_class_dich$steroids_3 = as.numeric(access_class_dich$steroids_3)

Make UC only cohort

UC_meds_3 <- access_class_dich[ which(access_class_dich$ibd_3=='UC'),]
View(UC_meds_3)

Baseline Characteristics

By medication class

access_class_dich %>% 
  dplyr::select(ibd_3, age_yrs, gender, race_5, ethnic_3, lang_3, max_ch, ASA_3, immuno_3, biologic_3, small_3, steroids_3, RPL_THEMES, RPL_4, RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4) -> baseline
baseline %>% tbl_summary(label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index", ASA_3 ~ "ASA Accessed",immuno_3 ~ "Immunomodulator Accessed", biologic_3 ~ "Biologic Accessed", small_3 ~ "Small Molecule Accessed", steroids_3 ~ "Steroids Accessed" ),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)")

Characteristic	N = 7,8361
IBD Diagnosis
CD	4,012 (51%)
UC	3,784 (48%)
Unspecified	40 (0.5%)
Age	45 (19)
Gender
Male	3,647 (47%)
Female	4,189 (53%)
Race
White	6,753 (86%)
Black	555 (7.1%)
Asian	237 (3.0%)
Native	29 (0.4%)
Other	262 (3.3%)
Ethnicity
NonHispanic	7,660 (98%)
UNKNOWN	0 (0%)
Hispanic	176 (2.2%)
CHOOSE NOT TO DISCLOSE	0 (0%)
Primary Language
English	7,754 (99%)
Other	82 (1.0%)
Charlson Comorbidity Index	3.2 (5.0)
ASA Accessed	504 (6.4%)
Immunomodulator Accessed	204 (2.6%)
Biologic Accessed	2,490 (32%)
Small Molecule Accessed	30 (0.4%)
Steroids Accessed	1,303 (17%)
Total SVI	0.36 (0.26)
SVI Quartiles
First	3,199 (41%)
Second	2,392 (31%)
Third	1,480 (19%)
Fourth	765 (9.8%)
Soceioeconomic Status	0.34 (0.26)
Household Composition	0.38 (0.26)
Minority Status and Language	0.49 (0.29)
Housing and Transportation	0.43 (0.28)
1 n (%); Mean (SD)

By individual medication

rx_class_final %>% 
  dplyr::select(ibd_3, age_yrs, gender, race_5, ethnic_3, lang_3, max_ch, mesalamine_3, azathioprine_3, mercaptopurine_3, methotrexate_3, infliximab_3, adalimumab_3, certolizumab_3, golimumab_3, ustekinumab_3, vedolizumab_3, prednisone_3, budesonide_3, steroid_count, RPL_THEMES, RPL_4, RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4) -> baseline2
baseline2 %>% tbl_summary(label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index", mesalamine_3 ~ "Mesalamine", azathioprine_3 ~ "AZA", mercaptopurine_3 ~ "6MP", methotrexate_3 ~ "Methotrexate", infliximab_3 ~ "Infliximab", adalimumab_3 ~ "Adalimumab", certolizumab_3 ~ "Certolizumab", golimumab_3 ~  "Golimumab", ustekinumab_3 ~ "Ustekinumab", vedolizumab_3 ~ "Vedolizumab", prednisone_3 ~ "Prednisone", budesonide_3 ~ "Budesonide", steroid_count ~ "Total Steroid Prescriptions"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)")

Characteristic	N = 7,8361
IBD Diagnosis
CD	4,012 (51%)
UC	3,784 (48%)
Unspecified	40 (0.5%)
Age	45 (19)
Gender
Male	3,647 (47%)
Female	4,189 (53%)
Race
White	6,753 (86%)
Black	555 (7.1%)
Asian	237 (3.0%)
Native	29 (0.4%)
Other	262 (3.3%)
Ethnicity
NonHispanic	7,660 (98%)
UNKNOWN	0 (0%)
Hispanic	176 (2.2%)
CHOOSE NOT TO DISCLOSE	0 (0%)
Primary Language
English	7,754 (99%)
Other	82 (1.0%)
Charlson Comorbidity Index	3.2 (5.0)
Mesalamine	487 (6.2%)
AZA	86 (1.1%)
6MP	22 (0.3%)
Methotrexate	96 (1.2%)
Infliximab	1,774 (23%)
Adalimumab	191 (2.4%)
Certolizumab	16 (0.2%)
Golimumab	8 (0.1%)
Ustekinumab	211 (2.7%)
Vedolizumab	667 (8.5%)
Prednisone	783 (10.0%)
Budesonide	167 (2.1%)
Total Steroid Prescriptions	2.08 (4.89)
Total SVI	0.36 (0.26)
SVI Quartiles
First	3,199 (41%)
Second	2,392 (31%)
Third	1,480 (19%)
Fourth	765 (9.8%)
Soceioeconomic Status	0.34 (0.26)
Household Composition	0.38 (0.26)
Minority Status and Language	0.49 (0.29)
Housing and Transportation	0.43 (0.28)
1 n (%); Mean (SD)

By SVI Quartile

baseline %>%  
tbl_summary(by = RPL_4,
         label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index", mesalamine_3 ~ "Mesalamine", azathioprine_3 ~ "AZA", mercaptopurine_3 ~ "6MP", methotrexate_3 ~ "Methotrexate", infliximab_3 ~ "Infliximab", adalimumab_3 ~ "Adalimumab", certolizumab_3 ~ "Certolizumab", golimumab_3 ~  "Golimumab", ustekinumab_3 ~ "Ustekinumab", vedolizumab_3 ~ "Vedolizumab", prednisone_3 ~ "Prednisone", budesonide_3 ~ "Budesonide", steroid_count ~ "Total Steroid Prescriptions"),
statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()

There was an error in 'add_p()/add_difference()' for variable 'ibd_3', p-value omitted:
Error in stats::fisher.test(structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, : FEXACT[f3xact()] error: hash key 2e+10 > INT_MAX, kyy=3200, it[i (= nco = 4)]= 0.
Rather set 'simulate.p.value=TRUE'

There was an error in 'add_p()/add_difference()' for variable 'race_5', p-value omitted:
Error in stats::fisher.test(structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, : FEXACT error 5.
The hash table key cannot be computed because the largest key
is larger than the largest representable int.
The algorithm cannot proceed.
Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.

Characteristic	First, N = 3,1991	Second, N = 2,3921	Third, N = 1,4801	Fourth, N = 7651	p-value2
IBD Diagnosis
CD	1,537 (48%)	1,175 (49%)	821 (55%)	479 (63%)
UC	1,647 (51%)	1,206 (50%)	648 (44%)	283 (37%)
Unspecified	15 (0.5%)	11 (0.5%)	11 (0.7%)	3 (0.4%)
Age	46 (20)	45 (19)	45 (19)	43 (18)	0.007
Gender					0.2
Male	1,516 (47%)	1,130 (47%)	663 (45%)	338 (44%)
Female	1,683 (53%)	1,262 (53%)	817 (55%)	427 (56%)
Race
White	2,887 (90%)	2,092 (87%)	1,252 (85%)	522 (68%)
Black	83 (2.6%)	131 (5.5%)	144 (9.7%)	197 (26%)
Asian	103 (3.2%)	93 (3.9%)	35 (2.4%)	6 (0.8%)
Native	9 (0.3%)	14 (0.6%)	2 (0.1%)	4 (0.5%)
Other	117 (3.7%)	62 (2.6%)	47 (3.2%)	36 (4.7%)
Ethnicity					0.022
NonHispanic	3,129 (98%)	2,349 (98%)	1,446 (98%)	736 (96%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Hispanic	70 (2.2%)	43 (1.8%)	34 (2.3%)	29 (3.8%)
CHOOSE NOT TO DISCLOSE	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Primary Language					<0.001
English	3,181 (99%)	2,366 (99%)	1,458 (99%)	749 (98%)
Other	18 (0.6%)	26 (1.1%)	22 (1.5%)	16 (2.1%)
Charlson Comorbidity Index	3.1 (4.9)	3.4 (5.2)	3.3 (5.0)	3.0 (4.8)	0.14
Mesalamine	213 (6.7%)	155 (6.5%)	83 (5.6%)	36 (4.7%)	0.2
AZA	32 (1.0%)	29 (1.2%)	16 (1.1%)	9 (1.2%)	0.9
6MP	11 (0.3%)	3 (0.1%)	5 (0.3%)	3 (0.4%)	0.3
Methotrexate	42 (1.3%)	33 (1.4%)	12 (0.8%)	9 (1.2%)	0.4
Infliximab	775 (24%)	511 (21%)	324 (22%)	164 (21%)	0.048
Adalimumab	69 (2.2%)	63 (2.6%)	37 (2.5%)	22 (2.9%)	0.6
Certolizumab	7 (0.2%)	5 (0.2%)	2 (0.1%)	2 (0.3%)	0.9
Golimumab	2 (<0.1%)	2 (<0.1%)	3 (0.2%)	1 (0.1%)	0.4
Ustekinumab	89 (2.8%)	57 (2.4%)	37 (2.5%)	28 (3.7%)	0.3
Vedolizumab	295 (9.2%)	209 (8.7%)	116 (7.8%)	47 (6.1%)	0.035
Prednisone	308 (9.6%)	231 (9.7%)	142 (9.6%)	102 (13%)	0.015
Budesonide	57 (1.8%)	62 (2.6%)	27 (1.8%)	21 (2.7%)	0.10
Total Steroid Prescriptions	1.91 (4.33)	2.21 (5.30)	2.00 (4.77)	2.58 (5.86)	<0.001
Total SVI	0.12 (0.07)	0.37 (0.07)	0.62 (0.07)	0.86 (0.07)	<0.001
Soceioeconomic Status	0.13 (0.11)	0.33 (0.16)	0.57 (0.16)	0.79 (0.11)	<0.001
Household Composition	0.21 (0.14)	0.36 (0.21)	0.57 (0.25)	0.76 (0.20)	<0.001
Minority Status and Language	0.42 (0.27)	0.49 (0.29)	0.54 (0.29)	0.67 (0.25)	<0.001
Housing and Transportation	0.19 (0.16)	0.50 (0.20)	0.63 (0.22)	0.81 (0.17)	<0.001
1 n (%); Mean (SD)
2 Kruskal-Wallis rank sum test; Pearson's Chi-squared test; Fisher's exact test

Bivariate analysis by class

ASA

UC_meds_3$ASA_3 = as.numeric(UC_meds_3$ASA_3)
tbl_uv_ASA_3<-
  tbl_uvregression(
    UC_meds_3[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "ASA_3")],
    method = glm,
    y = ASA_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_ASA_3, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	0.99	0.99, 1.00	<0.001
Gender	3,784
Male		—	—
Female		1.39	1.14, 1.70	0.001
Race	3,784
White		—	—
Black		1.18	0.77, 1.75	0.4
Asian		1.79	1.16, 2.68	0.006
Native		0.97	0.15, 3.44	>0.9
Other		1.93	1.26, 2.87	0.002
Ethnicity	3,784
NonHispanic		—	—
Hispanic		1.31	0.75, 2.15	0.3
Primary Language	3,784
English		—	—
Other		1.88	0.88, 3.64	0.079
Charlson Comorbidity Index	3,784	0.96	0.94, 0.98	<0.001
Total SVI	3,784	0.93	0.62, 1.38	0.7
SVI Quartiles	3,784
First		—	—
Second		1.01	0.81, 1.27	>0.9
Third		1.05	0.79, 1.37	0.7
Fourth		0.83	0.54, 1.23	0.4
Soceioeconomic Status	3,784	0.90	0.60, 1.33	0.6
Household Composition	3,784	0.68	0.46, 1.00	0.051
Minority Status and Language	3,784	1.36	0.97, 1.90	0.075
Housing and Transportation	3,784	0.97	0.69, 1.38	0.9
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Immunomodulators

tbl_uv_immuno<-
  tbl_uvregression(
    access_class_dich[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "immuno_3")],
    method = glm,
    y = immuno_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_immuno, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.98	0.98, 0.99	<0.001
Gender	7,836
Male		—	—
Female		0.65	0.49, 0.86	0.003
Race	7,836
White		—	—
Black		0.67	0.33, 1.20	0.2
Asian		1.43	0.67, 2.68	0.3
Native		1.30	0.07, 6.12	0.8
Other		0.42	0.10, 1.12	0.14
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.33	0.52, 2.79	0.5
Primary Language	7,836
English		—	—
Other		1.43	0.35, 3.86	0.5
Charlson Comorbidity Index	7,836	1.04	1.02, 1.07	<0.001
Total SVI	7,836	0.93	0.53, 1.59	0.8
SVI Quartiles	7,836
First		—	—
Second		1.02	0.74, 1.42	0.9
Third		0.84	0.55, 1.24	0.4
Fourth		1.03	0.62, 1.65	0.9
Soceioeconomic Status	7,836	1.34	0.78, 2.27	0.3
Household Composition	7,836	0.90	0.53, 1.52	0.7
Minority Status and Language	7,836	0.84	0.52, 1.36	0.5
Housing and Transportation	7,836	0.71	0.43, 1.16	0.2
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Biologics

tbl_uv_biologics<-
  tbl_uvregression(
    access_class_dich[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "biologic_3")],
    method = glm,
    y = biologic_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_biologics, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.97	0.96, 0.97	<0.001
Gender	7,836
Male		—	—
Female		0.74	0.67, 0.81	<0.001
Race	7,836
White		—	—
Black		1.37	1.14, 1.63	<0.001
Asian		0.95	0.71, 1.25	0.7
Native		0.70	0.28, 1.57	0.4
Other		1.32	1.02, 1.70	0.032
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.14	0.83, 1.56	0.4
Primary Language	7,836
English		—	—
Other		0.89	0.54, 1.41	0.6
Charlson Comorbidity Index	7,836	0.95	0.94, 0.96	<0.001
Total SVI	7,836	0.82	0.68, 0.98	0.034
SVI Quartiles	7,836
First		—	—
Second		0.85	0.76, 0.95	0.006
Third		0.88	0.77, 1.00	0.053
Fourth		0.87	0.74, 1.03	0.12
Soceioeconomic Status	7,836	0.85	0.70, 1.02	0.085
Household Composition	7,836	0.77	0.65, 0.93	0.006
Minority Status and Language	7,836	1.49	1.26, 1.75	<0.001
Housing and Transportation	7,836	0.69	0.58, 0.81	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Small Molecules (UC only)

UC_meds_3$small = as.numeric(UC_meds_3$small_3)
tbl_uv_small_3<-
  tbl_uvregression(
    UC_meds_3[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "small_3")],
    method = glm,
    y = small_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_small_3, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	0.99	0.97, 1.01	0.5
Gender	3,784
Male		—	—
Female		0.73	0.34, 1.55	0.4
Race	3,784
White		—	—
Black		1.28	0.20, 4.37	0.7
Asian		0.91	0.05, 4.35	>0.9
Native		0.00		>0.9
Other		1.81	0.29, 6.20	0.4
Ethnicity	3,784
NonHispanic		—	—
Hispanic		0.00		>0.9
Primary Language	3,784
English		—	—
Other		0.00		>0.9
Charlson Comorbidity Index	3,784	1.00	0.92, 1.06	>0.9
Total SVI	3,784	0.42	0.07, 2.05	0.3
SVI Quartiles	3,784
First		—	—
Second		0.82	0.34, 1.84	0.6
Third		0.51	0.12, 1.54	0.3
Fourth		0.39	0.02, 1.91	0.4
Soceioeconomic Status	3,784	0.64	0.12, 2.92	0.6
Household Composition	3,784	0.46	0.09, 2.05	0.3
Minority Status and Language	3,784	1.35	0.37, 4.95	0.6
Housing and Transportation	3,784	0.39	0.09, 1.51	0.2
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Steroids

tbl_uv_steroids<-
  tbl_uvregression(
    access_class_dich[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "steroids_3")],
    method = glm,
    y = steroids_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_steroids, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	1.00	1.00, 1.01	0.028
Gender	7,836
Male		—	—
Female		1.04	0.93, 1.18	0.5
Race	7,836
White		—	—
Black		1.68	1.37, 2.05	<0.001
Asian		0.58	0.37, 0.87	0.012
Native		1.34	0.49, 3.09	0.5
Other		0.76	0.52, 1.09	0.2
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.21	0.81, 1.74	0.3
Primary Language	7,836
English		—	—
Other		0.54	0.24, 1.05	0.10
Charlson Comorbidity Index	7,836	1.10	1.09, 1.11	<0.001
Total SVI	7,836	1.29	1.03, 1.63	0.028
SVI Quartiles	7,836
First		—	—
Second		1.09	0.95, 1.26	0.2
Third		0.94	0.79, 1.11	0.5
Fourth		1.41	1.15, 1.71	<0.001
Soceioeconomic Status	7,836	1.44	1.14, 1.81	0.002
Household Composition	7,836	1.45	1.16, 1.81	0.001
Minority Status and Language	7,836	0.83	0.67, 1.02	0.076
Housing and Transportation	7,836	1.04	0.84, 1.28	0.7
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Bivariate Analysis by Medication

Mesalamine limited to UC patients only

tbl_uv_mes<-
  tbl_uvregression(
    UC_meds_3[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "mesalamine_3")],
    method = glm,
    y = mesalamine_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_mes, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	0.99	0.99, 1.00	<0.001
Gender	3,784
Male		—	—
Female		1.43	1.17, 1.76	<0.001
Race	3,784
White		—	—
Black		1.23	0.81, 1.82	0.3
Asian		1.87	1.21, 2.80	0.003
Native		1.01	0.16, 3.58	>0.9
Other		2.01	1.32, 2.99	<0.001
Ethnicity	3,784
NonHispanic		—	—
Hispanic		1.26	0.71, 2.09	0.4
Primary Language	3,784
English		—	—
Other		1.70	0.77, 3.37	0.2
Charlson Comorbidity Index	3,784	0.96	0.94, 0.98	<0.001
Total SVI	3,784	0.93	0.62, 1.40	0.7
SVI Quartiles	3,784
First		—	—
Second		0.99	0.78, 1.24	>0.9
Third		1.03	0.77, 1.35	0.8
Fourth		0.85	0.55, 1.26	0.4
Soceioeconomic Status	3,784	0.91	0.60, 1.36	0.6
Household Composition	3,784	0.69	0.46, 1.02	0.063
Minority Status and Language	3,784	1.38	0.98, 1.94	0.063
Housing and Transportation	3,784	0.97	0.68, 1.37	0.8
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Immunomodulators

## MTX
tbl_uv_mtx<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "methotrexate_3")],
    method = glm,
    y = methotrexate_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_mtx, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.98	0.97, 0.99	<0.001
Gender	7,836
Male		—	—
Female		0.67	0.45, 1.01	0.057
Race	7,836
White		—	—
Black		0.42	0.10, 1.11	0.14
Asian		1.32	0.40, 3.18	0.6
Native		2.74	0.15, 13.0	0.3
Other		0.29	0.02, 1.33	0.2
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.41	0.34, 3.80	0.6
Primary Language	7,836
English		—	—
Other		3.13	0.76, 8.57	0.056
Charlson Comorbidity Index	7,836	1.09	1.06, 1.12	<0.001
Total SVI	7,836	0.68	0.30, 1.51	0.4
SVI Quartiles	7,836
First		—	—
Second		1.05	0.66, 1.66	0.8
Third		0.61	0.31, 1.13	0.14
Fourth		0.89	0.41, 1.76	0.8
Soceioeconomic Status	7,836	1.10	0.50, 2.37	0.8
Household Composition	7,836	0.58	0.26, 1.26	0.2
Minority Status and Language	7,836	0.74	0.37, 1.49	0.4
Housing and Transportation	7,836	0.64	0.31, 1.32	0.2
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## 6MP
tbl_uv_mp<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "mercaptopurine_3")],
    method = glm,
    y = mercaptopurine_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_mp, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	1.00	0.97, 1.02	0.7
Gender	7,836
Male		—	—
Female		0.41	0.15, 0.96	0.049
Race	7,836
White		—	—
Black		0.58	0.03, 2.77	0.6
Asian		0.00		>0.9
Native		0.00		>0.9
Other		0.00		>0.9
Ethnicity	7,836
NonHispanic		—	—
Hispanic		0.00		>0.9
Primary Language	7,836
English		—	—
Other		0.00		>0.9
Charlson Comorbidity Index	7,836	0.92	0.78, 1.02	0.2
Total SVI	7,836	0.90	0.16, 4.45	>0.9
SVI Quartiles	7,836
First		—	—
Second		0.36	0.08, 1.17	0.12
Third		0.98	0.31, 2.71	>0.9
Fourth		1.14	0.26, 3.66	0.8
Soceioeconomic Status	7,836	0.98	0.18, 4.79	>0.9
Household Composition	7,836	1.76	0.36, 7.94	0.5
Minority Status and Language	7,836	0.82	0.19, 3.50	0.8
Housing and Transportation	7,836	0.56	0.12, 2.45	0.4
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## AZA 
tbl_uv_aza<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "azathioprine_3")],
    method = glm,
    y = azathioprine_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_aza, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.99	0.98, 1.00	0.030
Gender	7,836
Male		—	—
Female		0.72	0.47, 1.10	0.13
Race	7,836
White		—	—
Black		1.00	0.39, 2.12	>0.9
Asian		1.97	0.69, 4.46	0.15
Native		0.00	0.00, 31,531	>0.9
Other		0.70	0.12, 2.25	0.6
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.58	0.39, 4.28	0.4
Primary Language	7,836
English		—	—
Other		0.00	0.00, 7.01	>0.9
Charlson Comorbidity Index	7,836	0.99	0.95, 1.03	0.7
Total SVI	7,836	1.30	0.56, 2.90	0.5
SVI Quartiles	7,836
First		—	—
Second		1.21	0.73, 2.01	0.5
Third		1.08	0.58, 1.95	0.8
Fourth		1.18	0.53, 2.38	0.7
Soceioeconomic Status	7,836	1.76	0.78, 3.88	0.2
Household Composition	7,836	1.21	0.54, 2.64	0.6
Minority Status and Language	7,836	0.97	0.46, 2.03	>0.9
Housing and Transportation	7,836	0.85	0.39, 1.78	0.7
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Biologics

## IFX
tbl_uv_ifx<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "infliximab_3")],
    method = glm,
    y = infliximab_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_ifx, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.96	0.96, 0.96	<0.001
Gender	7,836
Male		—	—
Female		0.71	0.64, 0.79	<0.001
Race	7,836
White		—	—
Black		1.48	1.22, 1.78	<0.001
Asian		1.05	0.76, 1.42	0.8
Native		0.57	0.17, 1.47	0.3
Other		1.27	0.95, 1.67	0.10
Ethnicity	7,836
NonHispanic		—	—
Hispanic		0.91	0.62, 1.29	0.6
Primary Language	7,836
English		—	—
Other		1.03	0.60, 1.69	>0.9
Charlson Comorbidity Index	7,836	0.92	0.91, 0.93	<0.001
Total SVI	7,836	0.81	0.66, 1.00	0.046
SVI Quartiles	7,836
First		—	—
Second		0.85	0.75, 0.96	0.012
Third		0.88	0.76, 1.02	0.080
Fourth		0.85	0.70, 1.03	0.10
Soceioeconomic Status	7,836	0.84	0.68, 1.04	0.11
Household Composition	7,836	0.79	0.65, 0.97	0.026
Minority Status and Language	7,836	1.39	1.16, 1.67	<0.001
Housing and Transportation	7,836	0.68	0.57, 0.82	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## 6MP
tbl_uv_ada<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "adalimumab_3")],
    method = glm,
    y = adalimumab_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_ada, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.96	0.95, 0.97	<0.001
Gender	7,836
Male		—	—
Female		0.94	0.70, 1.25	0.6
Race	7,836
White		—	—
Black		1.57	0.95, 2.46	0.062
Asian		1.09	0.43, 2.28	0.8
Native		1.50	0.08, 7.09	0.7
Other		1.15	0.49, 2.31	0.7
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.17	0.41, 2.61	0.7
Primary Language	7,836
English		—	—
Other		0.00	0.00, 3.99	>0.9
Charlson Comorbidity Index	7,836	0.91	0.86, 0.95	<0.001
Total SVI	7,836	1.38	0.79, 2.38	0.2
SVI Quartiles	7,836
First		—	—
Second		1.23	0.87, 1.73	0.2
Third		1.16	0.77, 1.73	0.5
Fourth		1.34	0.81, 2.15	0.2
Soceioeconomic Status	7,836	1.49	0.86, 2.56	0.2
Household Composition	7,836	1.28	0.74, 2.17	0.4
Minority Status and Language	7,836	1.68	1.02, 2.78	0.041
Housing and Transportation	7,836	0.88	0.53, 1.45	0.6
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## Vedo 
tbl_uv_vedo<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "vedolizumab_3")],
    method = glm,
    y = vedolizumab_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_vedo, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	1.00	0.99, 1.00	0.4
Gender	7,836
Male		—	—
Female		1.00	0.85, 1.17	>0.9
Race	7,836
White		—	—
Black		0.91	0.65, 1.25	0.6
Asian		0.89	0.53, 1.41	0.7
Native		0.81	0.13, 2.69	0.8
Other		1.68	1.14, 2.38	0.006
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.47	0.90, 2.28	0.10
Primary Language	7,836
English		—	—
Other		1.33	0.62, 2.53	0.4
Charlson Comorbidity Index	7,836	1.03	1.02, 1.05	<0.001
Total SVI	7,836	0.61	0.44, 0.84	0.003
SVI Quartiles	7,836
First		—	—
Second		0.94	0.78, 1.13	0.5
Third		0.84	0.67, 1.04	0.12
Fourth		0.64	0.46, 0.88	0.007
Soceioeconomic Status	7,836	0.68	0.49, 0.94	0.019
Household Composition	7,836	0.55	0.40, 0.75	<0.001
Minority Status and Language	7,836	1.30	0.99, 1.71	0.063
Housing and Transportation	7,836	0.65	0.49, 0.87	0.003
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## uste
tbl_uv_uste<-
  tbl_uvregression(
    baseline[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "ustekinumab_3")],
    method = glm,
    y = ustekinumab_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_uste, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	7,836	0.98	0.97, 0.99	<0.001
Gender	7,836
Male		—	—
Female		0.91	0.69, 1.20	0.5
Race	7,836
White		—	—
Black		1.50	0.93, 2.30	0.079
Asian		0.62	0.19, 1.49	0.4
Native		1.30	0.07, 6.12	0.8
Other		0.42	0.10, 1.12	0.14
Ethnicity	7,836
NonHispanic		—	—
Hispanic		0.84	0.26, 2.00	0.7
Primary Language	7,836
English		—	—
Other		0.90	0.15, 2.89	0.9
Charlson Comorbidity Index	7,836	0.98	0.95, 1.00	0.12
Total SVI	7,836	1.24	0.73, 2.09	0.4
SVI Quartiles	7,836
First		—	—
Second		0.85	0.61, 1.19	0.4
Third		0.90	0.60, 1.31	0.6
Fourth		1.33	0.85, 2.02	0.2
Soceioeconomic Status	7,836	1.14	0.67, 1.92	0.6
Household Composition	7,836	1.24	0.74, 2.06	0.4
Minority Status and Language	7,836	1.36	0.85, 2.20	0.2
Housing and Transportation	7,836	0.88	0.54, 1.43	0.6
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Small molecules limited to UC only patients

## Tofa 
tbl_uv_tofa<-
  tbl_uvregression(
    UC_meds_3[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "tofacitinib_3")],
    method = glm,
    y = tofacitinib_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )
print(tbl_uv_tofa, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	0.99	0.97, 1.01	0.5
Gender	3,784
Male		—	—
Female		0.73	0.34, 1.55	0.4
Race	3,784
White		—	—
Black		1.28	0.20, 4.37	0.7
Asian		0.91	0.05, 4.35	>0.9
Native		0.00		>0.9
Other		1.81	0.29, 6.20	0.4
Ethnicity	3,784
NonHispanic		—	—
Hispanic		0.00		>0.9
Primary Language	3,784
English		—	—
Other		0.00		>0.9
Charlson Comorbidity Index	3,784	1.00	0.92, 1.06	>0.9
Total SVI	3,784	0.42	0.07, 2.05	0.3
SVI Quartiles	3,784
First		—	—
Second		0.82	0.34, 1.84	0.6
Third		0.51	0.12, 1.54	0.3
Fourth		0.39	0.02, 1.91	0.4
Soceioeconomic Status	3,784	0.64	0.12, 2.92	0.6
Household Composition	3,784	0.46	0.09, 2.05	0.3
Minority Status and Language	3,784	1.35	0.37, 4.95	0.6
Housing and Transportation	3,784	0.39	0.09, 1.51	0.2
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## Upa 
tbl_uv_upa<-
  tbl_uvregression(
    UC_meds_3[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "upadacitinib_3")],
    method = glm,
    y = upadacitinib_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )

Warning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not converge! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

print(tbl_uv_upa, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	1.00	0.00, 33,234,572,863,151,126,253,589,885,910,374,792,622,009,002,930,739,352,140,257,320,052,974,679,648,356,012,118,124,517,226,347,799,556,268,642,517,016,984,103,244,776,589,737,823,977,045,904,041,344,183,447,140,882,591,323,890,986,127,607,033,250,480,587,244,620,166,054,550,460,586,521,182,482,216,256,275,726,589,547,969,092,765,351,936	>0.9
Gender	3,784
Male		—	—
Female		1.00	0.00, Inf	>0.9
Race	3,784
White		—	—
Black		1.00	0.00, Inf	>0.9
Asian		1.00	0.00, Inf	>0.9
Native		1.00	0.00, Inf	>0.9
Other		1.00	0.00, Inf	>0.9
Ethnicity	3,784
NonHispanic		—	—
Hispanic		1.00	0.00, Inf	>0.9
Primary Language	3,784
English		—	—
Other		1.00	0.00, Inf	>0.9
Charlson Comorbidity Index	3,784	1.00	0.00, Inf	>0.9
Total SVI	3,784	1.00	0.00, Inf	>0.9
SVI Quartiles	3,784
First		—	—
Second		1.00	0.00, Inf	>0.9
Third		1.00	0.00, Inf	>0.9
Fourth		1.00	0.00, Inf	>0.9
Soceioeconomic Status	3,784	1.00	0.00, Inf	>0.9
Household Composition	3,784	1.00	0.00, Inf	>0.9
Minority Status and Language	3,784	1.00	0.00, Inf	>0.9
Housing and Transportation	3,784	1.00	0.00, Inf	>0.9
1 OR = Odds Ratio, CI = Confidence Interval

NULL

## Ozanimod
tbl_uv_ozan<-
  tbl_uvregression(
    UC_meds_3[c("age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "ozanimod_3")],
    method = glm,
    y = ozanimod_3,
    method.args = list(family = binomial),
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartiles", max_ch ~ "Charlson Comorbidity Index"),
    exponentiate = TRUE
  )

Warning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not convergeWarning: glm.fit: algorithm did not converge! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.
! `broom::tidy()` failed to tidy the model.
✖ need at least two non-NA values to interpolateapprox(sp$y, sp$x, xout = cutoff)
✔ `tidy_parameters()` used instead.
ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

print(tbl_uv_ozan, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	1.00	0.00, 33,234,572,863,151,126,253,589,885,910,374,792,622,009,002,930,739,352,140,257,320,052,974,679,648,356,012,118,124,517,226,347,799,556,268,642,517,016,984,103,244,776,589,737,823,977,045,904,041,344,183,447,140,882,591,323,890,986,127,607,033,250,480,587,244,620,166,054,550,460,586,521,182,482,216,256,275,726,589,547,969,092,765,351,936	>0.9
Gender	3,784
Male		—	—
Female		1.00	0.00, Inf	>0.9
Race	3,784
White		—	—
Black		1.00	0.00, Inf	>0.9
Asian		1.00	0.00, Inf	>0.9
Native		1.00	0.00, Inf	>0.9
Other		1.00	0.00, Inf	>0.9
Ethnicity	3,784
NonHispanic		—	—
Hispanic		1.00	0.00, Inf	>0.9
Primary Language	3,784
English		—	—
Other		1.00	0.00, Inf	>0.9
Charlson Comorbidity Index	3,784	1.00	0.00, Inf	>0.9
Total SVI	3,784	1.00	0.00, Inf	>0.9
SVI Quartiles	3,784
First		—	—
Second		1.00	0.00, Inf	>0.9
Third		1.00	0.00, Inf	>0.9
Fourth		1.00	0.00, Inf	>0.9
Soceioeconomic Status	3,784	1.00	0.00, Inf	>0.9
Household Composition	3,784	1.00	0.00, Inf	>0.9
Minority Status and Language	3,784	1.00	0.00, Inf	>0.9
Housing and Transportation	3,784	1.00	0.00, Inf	>0.9
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Multivariable Models by Class (RPL_THEMES)

ASA

## Total SVI 
ASA_access <- glm(ASA_3 ~  age_yrs + gender + race_5 + ethnic_3 + lang_3 
                     + max_ch + RPL_THEMES,
              family = "binomial",
              data = UC_meds_3)
summary(ASA_access )

Call:
glm(formula = ASA_3 ~ age_yrs + gender + race_5 + ethnic_3 + 
    lang_3 + max_ch + RPL_THEMES, family = "binomial", data = UC_meds_3)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8679  -0.5441  -0.4834  -0.4223   2.4125  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.812443   0.166056 -10.915  < 2e-16 ***
age_yrs          -0.006741   0.003023  -2.230 0.025753 *  
genderFemale      0.356928   0.103059   3.463 0.000534 ***
race_5Black       0.128121   0.214084   0.598 0.549534    
race_5Asian       0.469260   0.218976   2.143 0.032115 *  
race_5Native     -0.078508   0.755811  -0.104 0.917271    
race_5Other       0.643663   0.229272   2.807 0.004994 ** 
ethnic_3Hispanic -0.083076   0.295729  -0.281 0.778773    
lang_3Other       0.337591   0.378089   0.893 0.371917    
max_ch           -0.023294   0.012248  -1.902 0.057187 .  
RPL_THEMES       -0.158865   0.213469  -0.744 0.456752    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2796.4  on 3783  degrees of freedom
Residual deviance: 2751.5  on 3773  degrees of freedom
AIC: 2773.5

Number of Fisher Scoring iterations: 5

broom::glance(ASA_access )

broom::tidy(ASA_access , exponentiate = TRUE)

tbl_regression(ASA_access, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
Age	0.99	0.99, 1.00	0.026
Gender
Male	—	—
Female	1.43	1.17, 1.75	<0.001
Race
White	—	—
Black	1.14	0.73, 1.70	0.5
Asian	1.60	1.02, 2.42	0.032
Native	0.92	0.15, 3.30	>0.9
Other	1.90	1.20, 2.95	0.005
Ethnicity
NonHispanic	—	—
Hispanic	0.92	0.50, 1.60	0.8
Primary Language
English	—	—
Other	1.40	0.64, 2.84	0.4
Charlson Comorbidity Index	0.98	0.95, 1.00	0.057
Total SVI	0.85	0.56, 1.29	0.5
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ASA_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
2773.513 | 2842.137 |     0.012 | 0.324 | 0.854 |    0.364 |   -59.248 |       5.288e-04 | 0.789

performance::check_model(ASA_access, panel = TRUE)

# Margins 
cplot(ASA_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Mesalamine Rx Given SVI")

## All Themes 
ASA_access2 <- glm(ASA_3 ~  age_yrs + gender + race_5 + ethnic_3 + lang_3 
                     + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = UC_meds_3)
summary(ASA_access2 )

Call:
glm(formula = ASA_3 ~ age_yrs + gender + race_5 + ethnic_3 + 
    lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = UC_meds_3)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8540  -0.5411  -0.4829  -0.4196   2.4255  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.800114   0.198825  -9.054  < 2e-16 ***
age_yrs          -0.006562   0.003030  -2.165 0.030355 *  
genderFemale      0.356414   0.103097   3.457 0.000546 ***
race_5Black       0.101615   0.217633   0.467 0.640564    
race_5Asian       0.411965   0.225647   1.826 0.067895 .  
race_5Native     -0.080771   0.756433  -0.107 0.914964    
race_5Other       0.623641   0.230490   2.706 0.006816 ** 
ethnic_3Hispanic -0.083917   0.295595  -0.284 0.776493    
lang_3Other       0.297934   0.379240   0.786 0.432096    
max_ch           -0.022898   0.012254  -1.869 0.061674 .  
RPL_THEME1        0.035467   0.283107   0.125 0.900303    
RPL_THEME2       -0.305557   0.250045  -1.222 0.221704    
RPL_THEME3        0.063012   0.185924   0.339 0.734677    
RPL_THEME4       -0.012995   0.209754  -0.062 0.950598    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2796.4  on 3783  degrees of freedom
Residual deviance: 2749.7  on 3770  degrees of freedom
AIC: 2777.7

Number of Fisher Scoring iterations: 5

broom::glance(ASA_access2 )

broom::tidy(ASA_access2 , exponentiate = TRUE)

tbl_regression(ASA_access2, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
Age	0.99	0.99, 1.00	0.030
Gender
Male	—	—
Female	1.43	1.17, 1.75	<0.001
Race
White	—	—
Black	1.11	0.71, 1.67	0.6
Asian	1.51	0.95, 2.32	0.068
Native	0.92	0.14, 3.30	>0.9
Other	1.87	1.17, 2.90	0.007
Ethnicity
NonHispanic	—	—
Hispanic	0.92	0.50, 1.60	0.8
Primary Language
English	—	—
Other	1.35	0.61, 2.73	0.4
Charlson Comorbidity Index	0.98	0.95, 1.00	0.062
Soceioeconomic Status	1.04	0.59, 1.80	>0.9
Household Composition	0.74	0.45, 1.20	0.2
Minority Status and Language	1.07	0.74, 1.53	0.7
Housing and Transportation	0.99	0.65, 1.49	>0.9
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ASA_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
2777.703 | 2865.042 |     0.013 | 0.324 | 0.854 |    0.363 |   -59.264 |       4.113e-04 | 0.790

performance::check_model(ASA_access2, panel = TRUE)

# Margins 
cplot(ASA_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of ASA Rx Given Theme 1")

cplot(ASA_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of ASA Rx Given THEME2")

cplot(ASA_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of ASA Rx Given THEME3")

cplot(ASA_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of ASA Rx Given THEME4")

Immunomodulator

# Total SVI 
immuno_access <- glm(immuno_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline)
summary(immuno_access )

Call:
glm(formula = immuno_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8345  -0.2581  -0.2061  -0.1596   3.2574  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -2.344969   0.212946 -11.012  < 2e-16 ***
ibd_3UC          -0.583910   0.152993  -3.817 0.000135 ***
ibd_3Unspecified -0.238358   1.022199  -0.233 0.815620    
age_yrs          -0.027577   0.004616  -5.974 2.32e-09 ***
genderFemale     -0.355211   0.144626  -2.456 0.014047 *  
race_5Black      -0.446640   0.338691  -1.319 0.187261    
race_5Asian       0.419083   0.364312   1.150 0.250004    
race_5Native      0.404722   1.028292   0.394 0.693886    
race_5Other      -1.041447   0.606626  -1.717 0.086018 .  
ethnic_3Hispanic  0.557254   0.442458   1.259 0.207868    
lang_3Other       0.582826   0.623475   0.935 0.349890    
max_ch            0.096174   0.014359   6.698 2.11e-11 ***
RPL_THEMES       -0.122577   0.291634  -0.420 0.674258    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1891.2  on 7835  degrees of freedom
Residual deviance: 1805.6  on 7823  degrees of freedom
AIC: 1831.6

Number of Fisher Scoring iterations: 7

broom::glance(immuno_access )

broom::tidy(immuno_access , exponentiate = TRUE)

tbl_regression(immuno_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.56	0.41, 0.75	<0.001
Unspecified	0.79	0.04, 3.73	0.8
Age	0.97	0.96, 0.98	<0.001
Gender
Male	—	—
Female	0.70	0.53, 0.93	0.014
Race
White	—	—
Black	0.64	0.31, 1.18	0.2
Asian	1.52	0.69, 2.94	0.3
Native	1.50	0.08, 7.25	0.7
Other	0.35	0.08, 0.98	0.086
Ethnicity
NonHispanic	—	—
Hispanic	1.75	0.66, 3.84	0.2
Primary Language
English	—	—
Other	1.79	0.42, 5.24	0.3
Charlson Comorbidity Index	1.10	1.07, 1.13	<0.001
Total SVI	0.88	0.50, 1.56	0.7
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(immuno_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1831.643 | 1922.207 |     0.015 | 0.158 | 0.480 |    0.115 |    -5.369 |           0.006 | 0.950

performance::check_model(immuno_access, panel = TRUE)

# Margins 
cplot(immuno_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given SVI")

## All themes
immuno_access2 <- glm(immuno_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline)
summary(immuno_access2 )

Call:
glm(formula = immuno_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9226  -0.2583  -0.2027  -0.1569   3.2314  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -2.149293   0.256840  -8.368  < 2e-16 ***
ibd_3UC          -0.575934   0.153338  -3.756 0.000173 ***
ibd_3Unspecified -0.282443   1.025334  -0.275 0.782960    
age_yrs          -0.027375   0.004619  -5.926 3.10e-09 ***
genderFemale     -0.354884   0.144766  -2.451 0.014229 *  
race_5Black      -0.492708   0.342152  -1.440 0.149860    
race_5Asian       0.526077   0.371515   1.416 0.156766    
race_5Native      0.388969   1.031450   0.377 0.706093    
race_5Other      -1.053039   0.606773  -1.735 0.082657 .  
ethnic_3Hispanic  0.567800   0.442271   1.284 0.199202    
lang_3Other       0.570205   0.622552   0.916 0.359711    
max_ch            0.096060   0.014351   6.693 2.18e-11 ***
RPL_THEME1        1.027842   0.400195   2.568 0.010218 *  
RPL_THEME2       -0.508736   0.360880  -1.410 0.158626    
RPL_THEME3       -0.275267   0.262561  -1.048 0.294458    
RPL_THEME4       -0.686773   0.300819  -2.283 0.022430 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1891.2  on 7835  degrees of freedom
Residual deviance: 1796.4  on 7820  degrees of freedom
AIC: 1828.4

Number of Fisher Scoring iterations: 7

broom::glance(immuno_access2 )

broom::tidy(immuno_access2 , exponentiate = TRUE)

tbl_regression(immuno_access2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.56	0.41, 0.76	<0.001
Unspecified	0.75	0.04, 3.60	0.8
Age	0.97	0.96, 0.98	<0.001
Gender
Male	—	—
Female	0.70	0.53, 0.93	0.014
Race
White	—	—
Black	0.61	0.29, 1.14	0.15
Asian	1.69	0.76, 3.33	0.2
Native	1.48	0.08, 7.19	0.7
Other	0.35	0.08, 0.97	0.083
Ethnicity
NonHispanic	—	—
Hispanic	1.76	0.67, 3.88	0.2
Primary Language
English	—	—
Other	1.77	0.41, 5.17	0.4
Charlson Comorbidity Index	1.10	1.07, 1.13	<0.001
Soceioeconomic Status	2.80	1.27, 6.09	0.010
Household Composition	0.60	0.30, 1.22	0.2
Minority Status and Language	0.76	0.45, 1.27	0.3
Housing and Transportation	0.50	0.28, 0.90	0.022
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(immuno_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1828.448 | 1939.911 |     0.016 | 0.158 | 0.479 |    0.115 |    -5.374 |           0.006 | 0.950

performance::check_model(immuno_access2, panel = TRUE)

# Margins 
cplot(immuno_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given Theme 1")

cplot(immuno_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given THEME2")

cplot(immuno_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given THEME3")

cplot(immuno_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given THEME4")

Biologics

# Total SVI 
biologic_access <- glm(biologic_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline)
summary(biologic_access )

Call:
glm(formula = biologic_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6682  -0.8710  -0.6122   1.0796   2.7391  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       1.284955   0.084680  15.174  < 2e-16 ***
ibd_3UC          -0.836952   0.053124 -15.755  < 2e-16 ***
ibd_3Unspecified -2.565503   0.731322  -3.508 0.000451 ***
age_yrs          -0.033454   0.001672 -20.007  < 2e-16 ***
genderFemale     -0.234009   0.051862  -4.512 6.42e-06 ***
race_5Black       0.309697   0.100180   3.091 0.001992 ** 
race_5Asian      -0.144650   0.154975  -0.933 0.350629    
race_5Native     -0.211205   0.450785  -0.469 0.639407    
race_5Other       0.232594   0.146557   1.587 0.112500    
ethnic_3Hispanic  0.024677   0.179770   0.137 0.890817    
lang_3Other       0.143213   0.267571   0.535 0.592488    
max_ch            0.015833   0.006367   2.487 0.012895 *  
RPL_THEMES       -0.521087   0.105155  -4.955 7.22e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9797.7  on 7835  degrees of freedom
Residual deviance: 8852.7  on 7823  degrees of freedom
AIC: 8878.7

Number of Fisher Scoring iterations: 5

broom::glance(biologic_access )

broom::tidy(biologic_access , exponentiate = TRUE)

tbl_regression(biologic_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.43	0.39, 0.48	<0.001
Unspecified	0.08	0.01, 0.25	<0.001
Age	0.97	0.96, 0.97	<0.001
Gender
Male	—	—
Female	0.79	0.71, 0.88	<0.001
Race
White	—	—
Black	1.36	1.12, 1.66	0.002
Asian	0.87	0.64, 1.17	0.4
Native	0.81	0.31, 1.87	0.6
Other	1.26	0.94, 1.68	0.11
Ethnicity
NonHispanic	—	—
Hispanic	1.02	0.72, 1.45	0.9
Primary Language
English	—	—
Other	1.15	0.67, 1.93	0.6
Charlson Comorbidity Index	1.02	1.00, 1.03	0.013
Total SVI	0.59	0.48, 0.73	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(biologic_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8878.714 | 8969.278 |     0.119 | 0.437 | 1.064 |    0.565 |      -Inf |       1.320e-04 | 0.618

performance::check_model(biologic_access, panel = TRUE)

# Margins 
cplot(biologic_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given SVI")

## All themes
biologic_access2 <- glm(biologic_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline)
summary(biologic_access2 )

Call:
glm(formula = biologic_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6595  -0.8662  -0.6093   1.0725   2.7284  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       1.234932   0.100777  12.254  < 2e-16 ***
ibd_3UC          -0.838950   0.053275 -15.747  < 2e-16 ***
ibd_3Unspecified -2.542845   0.730976  -3.479 0.000504 ***
age_yrs          -0.033506   0.001678 -19.963  < 2e-16 ***
genderFemale     -0.231555   0.051973  -4.455 8.38e-06 ***
race_5Black       0.222303   0.101595   2.188 0.028661 *  
race_5Asian      -0.238173   0.157965  -1.508 0.131616    
race_5Native     -0.198160   0.451025  -0.439 0.660404    
race_5Other       0.185677   0.147454   1.259 0.207949    
ethnic_3Hispanic  0.022071   0.180266   0.122 0.902556    
lang_3Other       0.084938   0.268362   0.317 0.751620    
max_ch            0.016404   0.006387   2.568 0.010222 *  
RPL_THEME1       -0.160785   0.147274  -1.092 0.274945    
RPL_THEME2       -0.068216   0.129846  -0.525 0.599331    
RPL_THEME3        0.327757   0.095222   3.442 0.000577 ***
RPL_THEME4       -0.492909   0.108826  -4.529 5.92e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9797.7  on 7835  degrees of freedom
Residual deviance: 8823.1  on 7820  degrees of freedom
AIC: 8855.1

Number of Fisher Scoring iterations: 5

broom::glance(biologic_access2 )

broom::tidy(biologic_access2 , exponentiate = TRUE)

tbl_regression(biologic_access2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.43	0.39, 0.48	<0.001
Unspecified	0.08	0.01, 0.26	<0.001
Age	0.97	0.96, 0.97	<0.001
Gender
Male	—	—
Female	0.79	0.72, 0.88	<0.001
Race
White	—	—
Black	1.25	1.02, 1.52	0.029
Asian	0.79	0.58, 1.07	0.13
Native	0.82	0.32, 1.90	0.7
Other	1.20	0.90, 1.60	0.2
Ethnicity
NonHispanic	—	—
Hispanic	1.02	0.71, 1.45	>0.9
Primary Language
English	—	—
Other	1.09	0.63, 1.82	0.8
Charlson Comorbidity Index	1.02	1.00, 1.03	0.010
Soceioeconomic Status	0.85	0.64, 1.14	0.3
Household Composition	0.93	0.72, 1.20	0.6
Minority Status and Language	1.39	1.15, 1.67	<0.001
Housing and Transportation	0.61	0.49, 0.76	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(biologic_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8855.065 | 8966.529 |     0.122 | 0.436 | 1.062 |    0.563 |      -Inf |       1.576e-04 | 0.619

performance::check_model(biologic_access2, panel = TRUE)

# Margins 
cplot(biologic_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Biologic Access Given Theme 1")

cplot(biologic_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Biologic Access Given THEME2")

cplot(biologic_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Biologic Access Given THEME3")

cplot(biologic_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Biologic Access Given THEME4")

Infliximab

infliximab_access <- glm(infliximab_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline2)
summary(infliximab_access )

Call:
glm(formula = infliximab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5056  -0.7419  -0.5205  -0.2813   2.7958  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.935834   0.090771  10.310  < 2e-16 ***
ibd_3UC          -0.784993   0.059865 -13.113  < 2e-16 ***
ibd_3Unspecified -2.771162   1.018461  -2.721 0.006510 ** 
age_yrs          -0.036700   0.001912 -19.196  < 2e-16 ***
genderFemale     -0.258076   0.057623  -4.479 7.51e-06 ***
race_5Black       0.396978   0.107807   3.682 0.000231 ***
race_5Asian      -0.102515   0.169525  -0.605 0.545364    
race_5Native     -0.395531   0.554894  -0.713 0.475967    
race_5Other       0.195414   0.159993   1.221 0.221939    
ethnic_3Hispanic -0.278084   0.205920  -1.350 0.176874    
lang_3Other       0.329616   0.290337   1.135 0.256255    
max_ch           -0.006600   0.007836  -0.842 0.399604    
RPL_THEMES       -0.525140   0.116888  -4.493 7.03e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8382.6  on 7835  degrees of freedom
Residual deviance: 7481.6  on 7823  degrees of freedom
AIC: 7507.6

Number of Fisher Scoring iterations: 6

broom::glance(infliximab_access )

broom::tidy(infliximab_access , exponentiate = TRUE)

tbl_regression(infliximab_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.46	0.41, 0.51	<0.001
Unspecified	0.06	0.00, 0.29	0.007
Age	0.96	0.96, 0.97	<0.001
Gender
Male	—	—
Female	0.77	0.69, 0.86	<0.001
Race
White	—	—
Black	1.49	1.20, 1.83	<0.001
Asian	0.90	0.64, 1.25	0.5
Native	0.67	0.19, 1.80	0.5
Other	1.22	0.88, 1.66	0.2
Ethnicity
NonHispanic	—	—
Hispanic	0.76	0.50, 1.12	0.2
Primary Language
English	—	—
Other	1.39	0.77, 2.42	0.3
Charlson Comorbidity Index	0.99	0.98, 1.01	0.4
Total SVI	0.59	0.47, 0.74	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(infliximab_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7507.605 | 7598.170 |     0.119 | 0.392 | 0.978 |    0.477 |      -Inf |       2.618e-04 | 0.691

performance::check_model(infliximab_access, panel = TRUE)

# Margins 
cplot(infliximab_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given SVI")

## All themes
infliximab_acess2 <- glm(infliximab_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline2)
summary(infliximab_acess2 )

Call:
glm(formula = infliximab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4394  -0.7397  -0.5198  -0.2766   2.7925  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.920437   0.109122   8.435  < 2e-16 ***
ibd_3UC          -0.784223   0.059985 -13.074  < 2e-16 ***
ibd_3Unspecified -2.750677   1.018104  -2.702   0.0069 ** 
age_yrs          -0.036760   0.001916 -19.191  < 2e-16 ***
genderFemale     -0.256613   0.057703  -4.447 8.70e-06 ***
race_5Black       0.333399   0.109365   3.049   0.0023 ** 
race_5Asian      -0.162541   0.172818  -0.941   0.3469    
race_5Native     -0.377090   0.554109  -0.681   0.4962    
race_5Other       0.159484   0.160876   0.991   0.3215    
ethnic_3Hispanic -0.277365   0.206120  -1.346   0.1784    
lang_3Other       0.299447   0.291442   1.027   0.3042    
max_ch           -0.006377   0.007852  -0.812   0.4167    
RPL_THEME1       -0.203731   0.163372  -1.247   0.2124    
RPL_THEME2       -0.004339   0.144107  -0.030   0.9760    
RPL_THEME3        0.216992   0.105877   2.049   0.0404 *  
RPL_THEME4       -0.479469   0.120844  -3.968 7.26e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8382.6  on 7835  degrees of freedom
Residual deviance: 7464.9  on 7820  degrees of freedom
AIC: 7496.9

Number of Fisher Scoring iterations: 6

broom::glance(infliximab_acess2 )

broom::tidy(infliximab_acess2 , exponentiate = TRUE)

tbl_regression(infliximab_acess2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.46	0.41, 0.51	<0.001
Unspecified	0.06	0.00, 0.30	0.007
Age	0.96	0.96, 0.97	<0.001
Gender
Male	—	—
Female	0.77	0.69, 0.87	<0.001
Race
White	—	—
Black	1.40	1.12, 1.73	0.002
Asian	0.85	0.60, 1.18	0.3
Native	0.69	0.20, 1.83	0.5
Other	1.17	0.85, 1.60	0.3
Ethnicity
NonHispanic	—	—
Hispanic	0.76	0.50, 1.12	0.2
Primary Language
English	—	—
Other	1.35	0.75, 2.35	0.3
Charlson Comorbidity Index	0.99	0.98, 1.01	0.4
Soceioeconomic Status	0.82	0.59, 1.12	0.2
Household Composition	1.00	0.75, 1.32	>0.9
Minority Status and Language	1.24	1.01, 1.53	0.040
Housing and Transportation	0.62	0.49, 0.78	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(infliximab_acess2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7496.929 | 7608.392 |     0.121 | 0.391 | 0.977 |    0.476 |      -Inf |       1.500e-04 | 0.692

performance::check_model(infliximab_acess2, panel = TRUE)

# Margins 
cplot(infliximab_acess2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Infliximab Access Given Theme 1")

cplot(infliximab_acess2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Infliximab Access Given THEME2")

cplot(infliximab_acess2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Infliximab Access Given THEME3")

cplot(infliximab_acess2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Infliximab Access Given THEME4")

Adalimumab

ada_access <- glm(adalimumab_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline2)
summary(ada_access )

Call:
glm(formula = adalimumab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5183  -0.2587  -0.1784  -0.1193   3.5359  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -1.929610   0.217872  -8.857  < 2e-16 ***
ibd_3UC           -1.069394   0.178884  -5.978 2.26e-09 ***
ibd_3Unspecified  -0.244499   1.022196  -0.239    0.811    
age_yrs           -0.038065   0.005325  -7.148 8.82e-13 ***
genderFemale       0.067628   0.149447   0.453    0.651    
race_5Black        0.264959   0.256382   1.033    0.301    
race_5Asian        0.192940   0.428658   0.450    0.653    
race_5Native       0.589577   1.035143   0.570    0.569    
race_5Other        0.094096   0.409763   0.230    0.818    
ethnic_3Hispanic   0.123315   0.482531   0.256    0.798    
lang_3Other      -13.672672 415.596155  -0.033    0.974    
max_ch            -0.011168   0.023844  -0.468    0.640    
RPL_THEMES         0.091152   0.292490   0.312    0.755    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1796.1  on 7835  degrees of freedom
Residual deviance: 1652.0  on 7823  degrees of freedom
AIC: 1678

Number of Fisher Scoring iterations: 16

broom::glance(ada_access )

broom::tidy(ada_access , exponentiate = TRUE)

tbl_regression(ada_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.34	0.24, 0.48	<0.001
Unspecified	0.78	0.04, 3.71	0.8
Age	0.96	0.95, 0.97	<0.001
Gender
Male	—	—
Female	1.07	0.80, 1.44	0.7
Race
White	—	—
Black	1.30	0.77, 2.11	0.3
Asian	1.21	0.47, 2.58	0.7
Native	1.80	0.10, 8.91	0.6
Other	1.10	0.45, 2.29	0.8
Ethnicity
NonHispanic	—	—
Hispanic	1.13	0.38, 2.65	0.8
Primary Language
English	—	—
Other	0.00	0.00, 58.9	>0.9
Charlson Comorbidity Index	0.99	0.94, 1.03	0.6
Total SVI	1.10	0.61, 1.93	0.8
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ada_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1677.965 | 1768.529 |     0.021 | 0.152 | 0.460 |    0.105 |    -4.715 |           0.007 | 0.953

performance::check_model(ada_access, panel = TRUE)

# Margins 
cplot(ada_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Immunomodulator Rx Given SVI")

## All themes
ada_access2 <- glm(adalimumab_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline2)
summary(ada_access2 )

Call:
glm(formula = adalimumab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5968  -0.2574  -0.1769  -0.1180   3.5978  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -2.08787    0.26869  -7.771 7.81e-15 ***
ibd_3UC           -1.06208    0.17903  -5.932 2.99e-09 ***
ibd_3Unspecified  -0.21461    1.02218  -0.210    0.834    
age_yrs           -0.03800    0.00533  -7.130 1.01e-12 ***
genderFemale       0.06853    0.14951   0.458    0.647    
race_5Black        0.16609    0.26023   0.638    0.523    
race_5Asian        0.10786    0.43688   0.247    0.805    
race_5Native       0.59851    1.03495   0.578    0.563    
race_5Other        0.03852    0.41004   0.094    0.925    
ethnic_3Hispanic   0.12329    0.48145   0.256    0.798    
lang_3Other      -13.71344  414.93287  -0.033    0.974    
max_ch            -0.01141    0.02387  -0.478    0.633    
RPL_THEME1         0.27396    0.41451   0.661    0.509    
RPL_THEME2         0.20887    0.36915   0.566    0.572    
RPL_THEME3         0.44372    0.27543   1.611    0.107    
RPL_THEME4        -0.48034    0.31119  -1.544    0.123    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1796.1  on 7835  degrees of freedom
Residual deviance: 1647.0  on 7820  degrees of freedom
AIC: 1679

Number of Fisher Scoring iterations: 16

broom::glance(ada_access2 )

broom::tidy(ada_access2 , exponentiate = TRUE)

tbl_regression(ada_access2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.35	0.24, 0.49	<0.001
Unspecified	0.81	0.05, 3.83	0.8
Age	0.96	0.95, 0.97	<0.001
Gender
Male	—	—
Female	1.07	0.80, 1.44	0.6
Race
White	—	—
Black	1.18	0.69, 1.93	0.5
Asian	1.11	0.42, 2.42	0.8
Native	1.82	0.10, 8.99	0.6
Other	1.04	0.42, 2.17	>0.9
Ethnicity
NonHispanic	—	—
Hispanic	1.13	0.39, 2.64	0.8
Primary Language
English	—	—
Other	0.00	0.00, 53.7	>0.9
Charlson Comorbidity Index	0.99	0.94, 1.03	0.6
Soceioeconomic Status	1.32	0.58, 2.95	0.5
Household Composition	1.23	0.60, 2.55	0.6
Minority Status and Language	1.56	0.91, 2.68	0.11
Housing and Transportation	0.62	0.33, 1.13	0.12
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ada_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1678.966 | 1790.429 |     0.023 | 0.152 | 0.459 |    0.105 |    -4.718 |           0.007 | 0.954

performance::check_model(ada_access2, panel = TRUE)

# Margins 
cplot(ada_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Adalimumab Access Given Theme 1")

cplot(ada_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Adalimumab Access Given THEME2")

cplot(ada_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Adalimumab Access Given THEME3")

cplot(ada_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Adalimumab Access Given THEME4")

Vedolizumab

vedolizumab_access <- glm(vedolizumab_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline2)
summary(vedolizumab_access )

Call:
glm(formula = vedolizumab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8052  -0.4416  -0.4036  -0.3679   2.5038  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -1.991454   0.130418 -15.270  < 2e-16 ***
ibd_3UC           -0.003793   0.082659  -0.046 0.963397    
ibd_3Unspecified -13.143216 228.805833  -0.057 0.954193    
age_yrs           -0.009210   0.002495  -3.691 0.000224 ***
genderFemale       0.030180   0.081863   0.369 0.712378    
race_5Black        0.017017   0.170682   0.100 0.920585    
race_5Asian       -0.129602   0.254667  -0.509 0.610816    
race_5Native      -0.198861   0.736158  -0.270 0.787057    
race_5Other        0.472693   0.201626   2.344 0.019058 *  
ethnic_3Hispanic   0.211379   0.255009   0.829 0.407156    
lang_3Other        0.337822   0.371132   0.910 0.362693    
max_ch             0.051149   0.008522   6.002 1.95e-09 ***
RPL_THEMES        -0.542136   0.170285  -3.184 0.001454 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4562.1  on 7835  degrees of freedom
Residual deviance: 4502.4  on 7823  degrees of freedom
AIC: 4528.4

Number of Fisher Scoring iterations: 14

broom::glance(vedolizumab_access )

broom::tidy(vedolizumab_access , exponentiate = TRUE)

tbl_regression(vedolizumab_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	1.00	0.85, 1.17	>0.9
Unspecified	0.00	0.00, 0.05	>0.9
Age	0.99	0.99, 1.00	<0.001
Gender
Male	—	—
Female	1.03	0.88, 1.21	0.7
Race
White	—	—
Black	1.02	0.72, 1.41	>0.9
Asian	0.88	0.52, 1.41	0.6
Native	0.82	0.13, 2.76	0.8
Other	1.60	1.06, 2.35	0.019
Ethnicity
NonHispanic	—	—
Hispanic	1.24	0.73, 1.99	0.4
Primary Language
English	—	—
Other	1.40	0.63, 2.76	0.4
Charlson Comorbidity Index	1.05	1.03, 1.07	<0.001
Total SVI	0.58	0.42, 0.81	0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(vedolizumab_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
4528.378 | 4618.942 |     0.008 | 0.278 | 0.759 |    0.287 |   -59.229 |       8.071e-04 | 0.845

performance::check_model(vedolizumab_access, panel = TRUE)

# Margins 
cplot(vedolizumab_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Vedolizumab Access Given SVI")

## All themes
vedolizumab_access2 <- glm(vedolizumab_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline2)
summary(vedolizumab_access2 )

Call:
glm(formula = vedolizumab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8060  -0.4454  -0.4035  -0.3598   2.5837  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -2.006522   0.157146 -12.769  < 2e-16 ***
ibd_3UC           -0.004880   0.082745  -0.059 0.952974    
ibd_3Unspecified -13.124175 228.176438  -0.058 0.954133    
age_yrs           -0.008918   0.002502  -3.565 0.000364 ***
genderFemale       0.030951   0.081939   0.378 0.705631    
race_5Black       -0.078350   0.172688  -0.454 0.650040    
race_5Asian       -0.235960   0.257898  -0.915 0.360225    
race_5Native      -0.199177   0.736960  -0.270 0.786954    
race_5Other        0.418920   0.201921   2.075 0.038017 *  
ethnic_3Hispanic   0.203365   0.254227   0.800 0.423749    
lang_3Other        0.226669   0.371530   0.610 0.541799    
max_ch             0.051796   0.008532   6.071 1.27e-09 ***
RPL_THEME1         0.104277   0.231879   0.450 0.652924    
RPL_THEME2        -0.516225   0.207644  -2.486 0.012915 *  
RPL_THEME3         0.234396   0.149475   1.568 0.116851    
RPL_THEME4        -0.338792   0.172100  -1.969 0.049002 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4562.1  on 7835  degrees of freedom
Residual deviance: 4491.9  on 7820  degrees of freedom
AIC: 4523.9

Number of Fisher Scoring iterations: 14

broom::glance(vedolizumab_access2 )

broom::tidy(vedolizumab_access2 , exponentiate = TRUE)

tbl_regression(vedolizumab_access2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	1.00	0.85, 1.17	>0.9
Unspecified	0.00	0.00, 0.05	>0.9
Age	0.99	0.99, 1.00	<0.001
Gender
Male	—	—
Female	1.03	0.88, 1.21	0.7
Race
White	—	—
Black	0.92	0.65, 1.28	0.7
Asian	0.79	0.46, 1.27	0.4
Native	0.82	0.13, 2.76	0.8
Other	1.52	1.01, 2.23	0.038
Ethnicity
NonHispanic	—	—
Hispanic	1.23	0.73, 1.97	0.4
Primary Language
English	—	—
Other	1.25	0.57, 2.47	0.5
Charlson Comorbidity Index	1.05	1.04, 1.07	<0.001
Soceioeconomic Status	1.11	0.70, 1.74	0.7
Household Composition	0.60	0.40, 0.90	0.013
Minority Status and Language	1.26	0.94, 1.70	0.12
Housing and Transportation	0.71	0.51, 1.00	0.049
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(vedolizumab_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
4523.881 | 4635.345 |     0.009 | 0.278 | 0.758 |    0.287 |   -59.273 |       8.071e-04 | 0.846

performance::check_model(vedolizumab_access2, panel = TRUE)

# Margins 
cplot(vedolizumab_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Vedolizumab Access Given Theme 1")

cplot(vedolizumab_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Vedolizumab Access Given THEME2")

cplot(vedolizumab_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Vedolizumab Access Given THEME3")

cplot(vedolizumab_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Vedolizumab Access Given THEME4")

Ustekinumab

ustekinumab_access <- glm(ustekinumab_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline2)
summary(ustekinumab_access )

Call:
glm(formula = ustekinumab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.4721  -0.2990  -0.1602  -0.1133   3.3177  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -2.260003   0.208089 -10.861  < 2e-16 ***
ibd_3UC           -1.815656   0.208943  -8.690  < 2e-16 ***
ibd_3Unspecified -13.446163 377.303942  -0.036    0.972    
age_yrs           -0.019362   0.004572  -4.235 2.29e-05 ***
genderFemale      -0.032177   0.142052  -0.227    0.821    
race_5Black        0.292511   0.244053   1.199    0.231    
race_5Asian       -0.344452   0.527465  -0.653    0.514    
race_5Native       0.549679   1.035561   0.531    0.596    
race_5Other       -0.828149   0.600155  -1.380    0.168    
ethnic_3Hispanic   0.146688   0.529568   0.277    0.782    
lang_3Other        0.403149   0.749602   0.538    0.591    
max_ch             0.022765   0.017566   1.296    0.195    
RPL_THEMES        -0.144867   0.280467  -0.517    0.605    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1941.6  on 7835  degrees of freedom
Residual deviance: 1789.9  on 7823  degrees of freedom
AIC: 1815.9

Number of Fisher Scoring iterations: 15

broom::glance(ustekinumab_access )

broom::tidy(ustekinumab_access , exponentiate = TRUE)

tbl_regression(ustekinumab_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.16	0.11, 0.24	<0.001
Unspecified	0.00	0.00, 30.7	>0.9
Age	0.98	0.97, 0.99	<0.001
Gender
Male	—	—
Female	0.97	0.73, 1.28	0.8
Race
White	—	—
Black	1.34	0.81, 2.12	0.2
Asian	0.71	0.21, 1.76	0.5
Native	1.73	0.10, 8.58	0.6
Other	0.44	0.11, 1.20	0.2
Ethnicity
NonHispanic	—	—
Hispanic	1.16	0.34, 2.90	0.8
Primary Language
English	—	—
Other	1.50	0.24, 5.20	0.6
Charlson Comorbidity Index	1.02	0.99, 1.06	0.2
Total SVI	0.87	0.50, 1.49	0.6
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ustekinumab_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1815.861 | 1906.426 |     0.019 | 0.160 | 0.478 |    0.114 |    -5.758 |           0.007 | 0.949

performance::check_model(ustekinumab_access, panel = TRUE)

# Margins 
cplot(ustekinumab_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Ustekinumab Access Given SVI")

## All themes
ustekinumab_access2 <- glm(ustekinumab_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline2)
summary(ustekinumab_access2 )

Call:
glm(formula = ustekinumab_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.4893  -0.2965  -0.1638  -0.1133   3.3611  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -2.354405   0.255513  -9.214  < 2e-16 ***
ibd_3UC           -1.815700   0.209015  -8.687  < 2e-16 ***
ibd_3Unspecified -13.427856 377.218570  -0.036    0.972    
age_yrs           -0.019396   0.004573  -4.242 2.22e-05 ***
genderFemale      -0.028041   0.142103  -0.197    0.844    
race_5Black        0.258051   0.248040   1.040    0.298    
race_5Asian       -0.420739   0.532830  -0.790    0.430    
race_5Native       0.543270   1.036255   0.524    0.600    
race_5Other       -0.850198   0.599960  -1.417    0.156    
ethnic_3Hispanic   0.140827   0.529723   0.266    0.790    
lang_3Other        0.392588   0.752152   0.522    0.602    
max_ch             0.022845   0.017592   1.299    0.194    
RPL_THEME1        -0.297369   0.404645  -0.735    0.462    
RPL_THEME2         0.286072   0.358570   0.798    0.425    
RPL_THEME3         0.312155   0.261017   1.196    0.232    
RPL_THEME4        -0.282064   0.297702  -0.947    0.343    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1941.6  on 7835  degrees of freedom
Residual deviance: 1786.8  on 7820  degrees of freedom
AIC: 1818.8

Number of Fisher Scoring iterations: 15

broom::glance(ustekinumab_access2 )

broom::tidy(ustekinumab_access2 , exponentiate = TRUE)

tbl_regression(ustekinumab_access2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.16	0.11, 0.24	<0.001
Unspecified	0.00	0.00, 31.0	>0.9
Age	0.98	0.97, 0.99	<0.001
Gender
Male	—	—
Female	0.97	0.74, 1.29	0.8
Race
White	—	—
Black	1.29	0.78, 2.06	0.3
Asian	0.66	0.19, 1.65	0.4
Native	1.72	0.09, 8.55	0.6
Other	0.43	0.10, 1.18	0.2
Ethnicity
NonHispanic	—	—
Hispanic	1.15	0.34, 2.88	0.8
Primary Language
English	—	—
Other	1.48	0.23, 5.19	0.6
Charlson Comorbidity Index	1.02	0.99, 1.06	0.2
Soceioeconomic Status	0.74	0.33, 1.63	0.5
Household Composition	1.33	0.66, 2.70	0.4
Minority Status and Language	1.37	0.82, 2.28	0.2
Housing and Transportation	0.75	0.42, 1.35	0.3
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ustekinumab_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1818.783 | 1930.247 |     0.019 | 0.160 | 0.478 |    0.114 |    -5.760 |           0.007 | 0.949

performance::check_model(ustekinumab_access2, panel = TRUE)

# Margins 
cplot(ustekinumab_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Ustekinumab Access Given Theme 1")

cplot(ustekinumab_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Ustekinumab Access Given THEME2")

cplot(ustekinumab_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Ustekinumab Access Given THEME3")

cplot(ustekinumab_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Ustekinumab Access Given THEME4")

Small Molecules

small_access <- glm(small_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = UC_meds_3)

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

Steroids

steroid_access <- glm(steroids_3 ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              family = "binomial",
              data = baseline)
summary(steroid_access )

Call:
glm(formula = steroids_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "binomial", 
    data = baseline)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6652  -0.5836  -0.5113  -0.4218   2.5477  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.473599   0.103000 -14.307  < 2e-16 ***
ibd_3UC           0.007799   0.064015   0.122   0.9030    
ibd_3Unspecified -1.329186   0.737475  -1.802   0.0715 .  
age_yrs          -0.016670   0.001979  -8.424  < 2e-16 ***
genderFemale      0.065334   0.063407   1.030   0.3028    
race_5Black       0.516818   0.112741   4.584 4.56e-06 ***
race_5Asian      -0.346263   0.223215  -1.551   0.1208    
race_5Native      0.297234   0.471472   0.630   0.5284    
race_5Other      -0.278414   0.202178  -1.377   0.1685    
ethnic_3Hispanic  0.426193   0.209469   2.035   0.0419 *  
lang_3Other      -0.280633   0.386769  -0.726   0.4681    
max_ch            0.129497   0.006495  19.938  < 2e-16 ***
RPL_THEMES        0.064212   0.127388   0.504   0.6142    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7051.5  on 7835  degrees of freedom
Residual deviance: 6593.4  on 7823  degrees of freedom
AIC: 6619.4

Number of Fisher Scoring iterations: 5

broom::glance(steroid_access )

broom::tidy(steroid_access , exponentiate = TRUE)

tbl_regression(steroid_access, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	1.01	0.89, 1.14	>0.9
Unspecified	0.26	0.04, 0.89	0.071
Age	0.98	0.98, 0.99	<0.001
Gender
Male	—	—
Female	1.07	0.94, 1.21	0.3
Race
White	—	—
Black	1.68	1.34, 2.09	<0.001
Asian	0.71	0.45, 1.07	0.12
Native	1.35	0.49, 3.19	0.5
Other	0.76	0.50, 1.11	0.2
Ethnicity
NonHispanic	—	—
Hispanic	1.53	1.00, 2.28	0.042
Primary Language
English	—	—
Other	0.76	0.33, 1.52	0.5
Charlson Comorbidity Index	1.14	1.12, 1.15	<0.001
Total SVI	1.07	0.83, 1.37	0.6
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(steroid_access)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
6619.371 | 6709.935 |     0.067 | 0.360 | 0.918 |    0.421 |      -Inf |       1.277e-04 | 0.741

performance::check_model(steroid_access, panel = TRUE)

# Margins 
cplot(steroid_access, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Steroid Access Given SVI")

## All themes
steroid_access2 <- glm(steroids_3 ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = baseline)
summary(steroid_access2 )

Call:
glm(formula = steroids_3 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = baseline)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7286  -0.5858  -0.5058  -0.4198   2.5506  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.385706   0.122075 -11.351  < 2e-16 ***
ibd_3UC           0.012909   0.064099   0.201   0.8404    
ibd_3Unspecified -1.349920   0.737672  -1.830   0.0673 .  
age_yrs          -0.016780   0.001981  -8.470  < 2e-16 ***
genderFemale      0.066016   0.063451   1.040   0.2981    
race_5Black       0.534237   0.115158   4.639  3.5e-06 ***
race_5Asian      -0.249034   0.225960  -1.102   0.2704    
race_5Native      0.301842   0.470878   0.641   0.5215    
race_5Other      -0.255745   0.202846  -1.261   0.2074    
ethnic_3Hispanic  0.440742   0.209710   2.102   0.0356 *  
lang_3Other      -0.240422   0.387393  -0.621   0.5349    
max_ch            0.129349   0.006500  19.900  < 2e-16 ***
RPL_THEME1        0.272533   0.179848   1.515   0.1297    
RPL_THEME2        0.094717   0.158879   0.596   0.5511    
RPL_THEME3       -0.216656   0.116173  -1.865   0.0622 .  
RPL_THEME4       -0.216059   0.132036  -1.636   0.1018    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7051.5  on 7835  degrees of freedom
Residual deviance: 6583.9  on 7820  degrees of freedom
AIC: 6615.9

Number of Fisher Scoring iterations: 5

broom::glance(steroid_access2 )

broom::tidy(steroid_access2 , exponentiate = TRUE)

tbl_regression(steroid_access2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	1.01	0.89, 1.15	0.8
Unspecified	0.26	0.04, 0.87	0.067
Age	0.98	0.98, 0.99	<0.001
Gender
Male	—	—
Female	1.07	0.94, 1.21	0.3
Race
White	—	—
Black	1.71	1.36, 2.13	<0.001
Asian	0.78	0.49, 1.19	0.3
Native	1.35	0.49, 3.20	0.5
Other	0.77	0.51, 1.14	0.2
Ethnicity
NonHispanic	—	—
Hispanic	1.55	1.02, 2.32	0.036
Primary Language
English	—	—
Other	0.79	0.34, 1.59	0.5
Charlson Comorbidity Index	1.14	1.12, 1.15	<0.001
Soceioeconomic Status	1.31	0.92, 1.87	0.13
Household Composition	1.10	0.81, 1.50	0.6
Minority Status and Language	0.81	0.64, 1.01	0.062
Housing and Transportation	0.81	0.62, 1.04	0.10
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(steroid_access2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
6615.868 | 6727.332 |     0.068 | 0.360 | 0.918 |    0.420 |      -Inf |       1.276e-04 | 0.742

performance::check_model(steroid_access2, panel = TRUE)

# Margins 
cplot(steroid_access2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Steroid Access Given Theme 1")

cplot(steroid_access2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of SSteroid Access Given THEME2")

cplot(steroid_access2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Steroide Access Given THEME3")

cplot(steroid_access2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Steroid Access Given THEME4")

Steroids - Poisson

steroid_count.p <- glm(steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
                          family = "poisson", 
               data = access_class_dich) 
summary(steroid_count.p)

Call:
glm(formula = steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, family = "poisson", 
    data = access_class_dich)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.0896  -1.7688  -0.9456   0.1930  25.7563  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.9536445  0.0254459  37.477  < 2e-16 ***
ibd_3UC          -0.1853562  0.0160509 -11.548  < 2e-16 ***
ibd_3Unspecified -0.8984805  0.1693805  -5.305 1.13e-07 ***
age_yrs          -0.0141004  0.0004826 -29.219  < 2e-16 ***
genderFemale      0.0493611  0.0157686   3.130 0.001746 ** 
race_5Black       0.3485969  0.0270175  12.903  < 2e-16 ***
race_5Asian      -0.4316746  0.0638844  -6.757 1.41e-11 ***
race_5Native     -0.1207442  0.1391074  -0.868 0.385398    
race_5Other       0.1285500  0.0440705   2.917 0.003535 ** 
ethnic_3Hispanic  0.4941812  0.0470679  10.499  < 2e-16 ***
lang_3Other      -0.4091470  0.1052017  -3.889 0.000101 ***
max_ch            0.0937966  0.0014200  66.052  < 2e-16 ***
RPL_THEMES        0.1100711  0.0314838   3.496 0.000472 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 38056  on 7835  degrees of freedom
Residual deviance: 33571  on 7823  degrees of freedom
AIC: 45936

Number of Fisher Scoring iterations: 6

broom::glance(steroid_count.p)

broom::tidy(steroid_count.p, exponentiate = TRUE)

tbl_regression(steroid_count.p, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.83	0.81, 0.86	<0.001
Unspecified	0.41	0.29, 0.56	<0.001
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.05	1.02, 1.08	0.002
Race
White	—	—
Black	1.42	1.34, 1.49	<0.001
Asian	0.65	0.57, 0.73	<0.001
Native	0.89	0.67, 1.15	0.4
Other	1.14	1.04, 1.24	0.004
Ethnicity
NonHispanic	—	—
Hispanic	1.64	1.49, 1.80	<0.001
Primary Language
English	—	—
Other	0.66	0.54, 0.81	<0.001
Charlson Comorbidity Index	1.10	1.10, 1.10	<0.001
Total SVI	1.12	1.05, 1.19	<0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
steroid_count.p.res <- resid(steroid_count.p)
plot(fitted(steroid_count.p), steroid_count.p.res, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Negative Binomial')
abline(0,0)

# NB regression more appropriate because residuals of the model are smaller 

# Likelihood ratio test 
pchisq(2 * (logLik(steroid_count.p) - logLik(steroid_count.p)), df = 1, lower.tail = FALSE)

'log Lik.' 1 (df=13)

# p-value of loglik is < 0.05 so NB regression is the more appropriate model 

# Model performance 
model_performance(steroid_count.p)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
45936.287 | 46026.851 |           0.439 | 4.748 | 2.072 |    -2.929 |           0.010

performance::check_model(steroid_count.p, panel = TRUE)

# Margins 
cplot(steroid_count.nb, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Steroid Rx Given SVI")

## All themes 
steroid_count.p2 <- glm(steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEME1
                           + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
                           family = "poisson", 
               data = access_class_dich) 
summary(steroid_count.p2)

Call:
glm(formula = steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "poisson", data = access_class_dich)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.2760  -1.7647  -0.9420   0.1924  25.5761  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.9587554  0.0303375  31.603  < 2e-16 ***
ibd_3UC          -0.1814101  0.0160637 -11.293  < 2e-16 ***
ibd_3Unspecified -0.9042024  0.1693905  -5.338 9.40e-08 ***
age_yrs          -0.0141059  0.0004831 -29.197  < 2e-16 ***
genderFemale      0.0511812  0.0157720   3.245 0.001174 ** 
race_5Black       0.3305843  0.0275427  12.003  < 2e-16 ***
race_5Asian      -0.3987125  0.0644696  -6.184 6.23e-10 ***
race_5Native     -0.1119765  0.1391298  -0.805 0.420915    
race_5Other       0.1298473  0.0441305   2.942 0.003257 ** 
ethnic_3Hispanic  0.5035655  0.0470240  10.709  < 2e-16 ***
lang_3Other      -0.4029164  0.1053751  -3.824 0.000131 ***
max_ch            0.0936875  0.0014200  65.976  < 2e-16 ***
RPL_THEME1        0.2455678  0.0444972   5.519 3.41e-08 ***
RPL_THEME2        0.1151459  0.0396422   2.905 0.003677 ** 
RPL_THEME3        0.0100464  0.0288857   0.348 0.727991    
RPL_THEME4       -0.2372832  0.0326819  -7.260 3.86e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 38056  on 7835  degrees of freedom
Residual deviance: 33488  on 7820  degrees of freedom
AIC: 45859

Number of Fisher Scoring iterations: 6

broom::glance(steroid_count.p2)

broom::tidy(steroid_count.p2, exponentiate = TRUE)

tbl_regression(steroid_count.p2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.83	0.81, 0.86	<0.001
Unspecified	0.40	0.28, 0.55	<0.001
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.05	1.02, 1.09	0.001
Race
White	—	—
Black	1.39	1.32, 1.47	<0.001
Asian	0.67	0.59, 0.76	<0.001
Native	0.89	0.67, 1.16	0.4
Other	1.14	1.04, 1.24	0.003
Ethnicity
NonHispanic	—	—
Hispanic	1.65	1.51, 1.81	<0.001
Primary Language
English	—	—
Other	0.67	0.54, 0.82	<0.001
Charlson Comorbidity Index	1.10	1.10, 1.10	<0.001
Soceioeconomic Status	1.28	1.17, 1.39	<0.001
Household Composition	1.12	1.04, 1.21	0.004
Minority Status and Language	1.01	0.95, 1.07	0.7
Housing and Transportation	0.79	0.74, 0.84	<0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
steroid_count.p2.res2 <- resid(steroid_count.p2)
plot(fitted(steroid_count.p2), steroid_count.p2.res2, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Negative Binomial')
abline(0,0)

# NB regression more appropriate because residuals of the model are smaller 

# Likelihood ratio test 
pchisq(2 * (logLik(steroid_count.p2) - logLik(steroid_count.p2)), df = 1, lower.tail = FALSE)

'log Lik.' 1 (df=16)

# p-value of loglik is < 0.05 so NB regression is the more appropriate model 

# Model performance 
model_performance(steroid_count.p2)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
45859.294 | 45970.758 |           0.445 | 4.745 | 2.069 |    -2.924 |           0.010

performance::check_model(steroid_count.p2, panel = TRUE)

# Margins 
cplot(steroid_count.p2, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Steroid Rx Given SVI")

Steroids - negative binomial

library(MASS)
steroid_count.nb <- glm.nb(steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
               data = access_class_dich) 
summary(steroid_count.nb)

Call:
glm.nb(formula = steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, data = access_class_dich, 
    init.theta = 0.5501620579, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8370  -1.1945  -0.4441   0.1016   7.4741  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.892189   0.057075  15.632  < 2e-16 ***
ibd_3UC          -0.195863   0.035686  -5.488 4.05e-08 ***
ibd_3Unspecified -0.883169   0.279659  -3.158 0.001588 ** 
age_yrs          -0.013840   0.001072 -12.916  < 2e-16 ***
genderFemale      0.059275   0.035267   1.681 0.092815 .  
race_5Black       0.388900   0.068057   5.714 1.10e-08 ***
race_5Asian      -0.435657   0.113149  -3.850 0.000118 ***
race_5Native     -0.108531   0.291855  -0.372 0.709992    
race_5Other       0.198292   0.101145   1.960 0.049939 *  
ethnic_3Hispanic  0.476181   0.119680   3.979 6.93e-05 ***
lang_3Other      -0.410571   0.191915  -2.139 0.032408 *  
max_ch            0.101363   0.003880  26.127  < 2e-16 ***
RPL_THEMES        0.139797   0.070812   1.974 0.048360 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.5502) family taken to be 1)

    Null deviance: 8544.3  on 7835  degrees of freedom
Residual deviance: 7655.0  on 7823  degrees of freedom
AIC: 28499

Number of Fisher Scoring iterations: 1

              Theta:  0.5502 
          Std. Err.:  0.0132 

 2 x log-likelihood:  -28471.0690 

broom::glance(steroid_count.nb)

broom::tidy(steroid_count.nb, exponentiate = TRUE)

tbl_regression(steroid_count.nb, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.82	0.77, 0.88	<0.001
Unspecified	0.41	0.24, 0.72	0.002
Age	0.99	0.98, 0.99	<0.001
Gender
Male	—	—
Female	1.06	0.99, 1.14	0.093
Race
White	—	—
Black	1.48	1.29, 1.69	<0.001
Asian	0.65	0.52, 0.81	<0.001
Native	0.90	0.52, 1.64	0.7
Other	1.22	1.01, 1.49	0.050
Ethnicity
NonHispanic	—	—
Hispanic	1.61	1.28, 2.04	<0.001
Primary Language
English	—	—
Other	0.66	0.46, 0.97	0.032
Charlson Comorbidity Index	1.11	1.10, 1.12	<0.001
Total SVI	1.15	1.00, 1.32	0.048
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
steroid_count.nb.res <- resid(steroid_count.nb)
plot(fitted(steroid_count.nb), steroid_count.nb.res, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Negative Binomial')
abline(0,0)

# NB regression more appropriate because residuals of the model are smaller 

# Likelihood ratio test 
pchisq(2 * (logLik(steroid_count.nb) - logLik(steroid_count.nb)), df = 1, lower.tail = FALSE)

'log Lik.' 1 (df=14)

# p-value of loglik is < 0.05 so NB regression is the more appropriate model 

# Model performance 
model_performance(steroid_count.nb)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
28499.069 | 28596.600 |           0.162 | 4.760 | 0.989 |    -2.304 |           0.010

performance::check_model(steroid_count.nb, panel = TRUE)

# Margins 
cplot(steroid_count.nb, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Steroid Rx Given SVI")

## All themes 
steroid_count.nb2 <- glm.nb(steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEME1
                           + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
               data = access_class_dich) 
summary(steroid_count.nb2)

Call:
glm.nb(formula = steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, data = access_class_dich, init.theta = 0.5521851452, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8434  -1.1929  -0.4419   0.1005   7.4426  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.914142   0.067923  13.459  < 2e-16 ***
ibd_3UC          -0.188313   0.035671  -5.279 1.30e-07 ***
ibd_3Unspecified -0.881386   0.278868  -3.161 0.001575 ** 
age_yrs          -0.014098   0.001072 -13.152  < 2e-16 ***
genderFemale      0.056295   0.035232   1.598 0.110074    
race_5Black       0.372441   0.069117   5.389 7.10e-08 ***
race_5Asian      -0.389933   0.114625  -3.402 0.000669 ***
race_5Native     -0.111786   0.291823  -0.383 0.701673    
race_5Other       0.176373   0.101537   1.737 0.082381 .  
ethnic_3Hispanic  0.471661   0.119752   3.939 8.19e-05 ***
lang_3Other      -0.389512   0.192407  -2.024 0.042928 *  
max_ch            0.101670   0.003876  26.232  < 2e-16 ***
RPL_THEME1        0.232069   0.099580   2.330 0.019781 *  
RPL_THEME2        0.165145   0.087896   1.879 0.060263 .  
RPL_THEME3       -0.023504   0.064346  -0.365 0.714907    
RPL_THEME4       -0.220805   0.073376  -3.009 0.002619 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.5522) family taken to be 1)

    Null deviance: 8563.5  on 7835  degrees of freedom
Residual deviance: 7654.7  on 7820  degrees of freedom
AIC: 28488

Number of Fisher Scoring iterations: 1

              Theta:  0.5522 
          Std. Err.:  0.0132 

 2 x log-likelihood:  -28454.1120 

broom::glance(steroid_count.nb2)

broom::tidy(steroid_count.nb2, exponentiate = TRUE)

tbl_regression(steroid_count.nb2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.83	0.77, 0.89	<0.001
Unspecified	0.41	0.24, 0.73	0.002
Age	0.99	0.98, 0.99	<0.001
Gender
Male	—	—
Female	1.06	0.99, 1.13	0.11
Race
White	—	—
Black	1.45	1.27, 1.67	<0.001
Asian	0.68	0.54, 0.85	<0.001
Native	0.89	0.52, 1.64	0.7
Other	1.19	0.98, 1.46	0.082
Ethnicity
NonHispanic	—	—
Hispanic	1.60	1.28, 2.03	<0.001
Primary Language
English	—	—
Other	0.68	0.47, 0.99	0.043
Charlson Comorbidity Index	1.11	1.10, 1.12	<0.001
Soceioeconomic Status	1.26	1.03, 1.54	0.020
Household Composition	1.18	0.99, 1.41	0.060
Minority Status and Language	0.98	0.86, 1.11	0.7
Housing and Transportation	0.80	0.69, 0.93	0.003
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
steroid_count.nb.res2 <- resid(steroid_count.nb2)
plot(fitted(steroid_count.nb2), steroid_count.nb.res2, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Negative Binomial')
abline(0,0)

# NB regression more appropriate because residuals of the model are smaller 

# Likelihood ratio test 
pchisq(2 * (logLik(steroid_count.nb2) - logLik(steroid_count.nb2)), df = 1, lower.tail = FALSE)

'log Lik.' 1 (df=17)

# p-value of loglik is < 0.05 so NB regression is the more appropriate model 

# Model performance 
model_performance(steroid_count.nb2)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
28488.112 | 28606.542 |           0.165 | 4.757 | 0.989 |    -2.294 |           0.010

performance::check_model(steroid_count.nb2, panel = TRUE)

# Margins 
cplot(steroid_count.nb, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Steroid Rx Given SVI")

Steroids - linear model

steroid_count <- lm(steroid_count ~  ibd_3 + age_yrs + gender + race_5 + 
                           ethnic_3 + lang_3  + max_ch + RPL_THEMES,
              data = access_class_dich)
summary(steroid_count)

Call:
lm(formula = steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEMES, data = access_class_dich)

Residuals:
    Min      1Q  Median      3Q     Max 
 -6.806  -1.814  -0.922   0.317 109.536 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       2.606611   0.175835  14.824  < 2e-16 ***
ibd_3UC          -0.359485   0.109398  -3.286  0.00102 ** 
ibd_3Unspecified -1.260233   0.754885  -1.669  0.09507 .  
age_yrs          -0.030943   0.003247  -9.529  < 2e-16 ***
genderFemale      0.105613   0.108114   0.977  0.32866    
race_5Black       0.898529   0.217129   4.138 3.54e-05 ***
race_5Asian      -0.558102   0.319979  -1.744  0.08117 .  
race_5Native     -0.276690   0.884767  -0.313  0.75450    
race_5Other       0.280816   0.315741   0.889  0.37382    
ethnic_3Hispanic  1.172306   0.382145   3.068  0.00216 ** 
lang_3Other      -0.565005   0.541629  -1.043  0.29691    
max_ch            0.254954   0.012364  20.621  < 2e-16 ***
RPL_THEMES        0.275780   0.218063   1.265  0.20602    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.749 on 7823 degrees of freedom
Multiple R-squared:  0.05882,   Adjusted R-squared:  0.05737 
F-statistic: 40.74 on 12 and 7823 DF,  p-value: < 2.2e-16

broom::glance(steroid_count)

broom::tidy(steroid_count)

tbl_regression(steroid_count, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index"))

Characteristic	Beta	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	-0.36	-0.57, -0.15	0.001
Unspecified	-1.3	-2.7, 0.22	0.10
Age	-0.03	-0.04, -0.02	<0.001
Gender
Male	—	—
Female	0.11	-0.11, 0.32	0.3
Race
White	—	—
Black	0.90	0.47, 1.3	<0.001
Asian	-0.56	-1.2, 0.07	0.081
Native	-0.28	-2.0, 1.5	0.8
Other	0.28	-0.34, 0.90	0.4
Ethnicity
NonHispanic	—	—
Hispanic	1.2	0.42, 1.9	0.002
Primary Language
English	—	—
Other	-0.57	-1.6, 0.50	0.3
Charlson Comorbidity Index	0.25	0.23, 0.28	<0.001
Total SVI	0.28	-0.15, 0.70	0.2
1 CI = Confidence Interval

# Model performance 
model_performance(steroid_count)

# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------
46667.357 | 46764.888 | 0.059 |     0.057 | 4.745 | 4.749

performance::check_model(steroid_count, panel = TRUE)

# Margins 
cplot(steroid_count, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Steroids Given SVI")

## All themes
steroid_count2 <- lm(steroid_count ~  ibd_3 + age_yrs + gender + 
                    race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                      RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              data = access_class_dich)
summary(steroid_count2 )

Call:
lm(formula = steroid_count ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, data = access_class_dich)

Residuals:
    Min      1Q  Median      3Q     Max 
 -6.963  -1.808  -0.921   0.329 109.405 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       2.61934    0.20940  12.509  < 2e-16 ***
ibd_3UC          -0.34813    0.10943  -3.181  0.00147 ** 
ibd_3Unspecified -1.27753    0.75473  -1.693  0.09055 .  
age_yrs          -0.03107    0.00325  -9.560  < 2e-16 ***
genderFemale      0.10645    0.10808   0.985  0.32471    
race_5Black       0.85986    0.22030   3.903 9.57e-05 ***
race_5Asian      -0.48405    0.32519  -1.489  0.13665    
race_5Native     -0.26893    0.88452  -0.304  0.76110    
race_5Other       0.28135    0.31640   0.889  0.37390    
ethnic_3Hispanic  1.19139    0.38210   3.118  0.00183 ** 
lang_3Other      -0.55244    0.54298  -1.017  0.30898    
max_ch            0.25486    0.01236  20.616  < 2e-16 ***
RPL_THEME1        0.58927    0.30596   1.926  0.05414 .  
RPL_THEME2        0.25186    0.26960   0.934  0.35024    
RPL_THEME3        0.01470    0.19745   0.074  0.94064    
RPL_THEME4       -0.50171    0.22520  -2.228  0.02592 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.747 on 7820 degrees of freedom
Multiple R-squared:  0.05981,   Adjusted R-squared:  0.05801 
F-statistic: 33.17 on 15 and 7820 DF,  p-value: < 2.2e-16

broom::glance(steroid_count2 )

broom::tidy(steroid_count2)

tbl_regression(steroid_count2, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", max_ch ~ "Charlson Comorbidity Index", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation"))

Characteristic	Beta	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	-0.35	-0.56, -0.13	0.001
Unspecified	-1.3	-2.8, 0.20	0.091
Age	-0.03	-0.04, -0.02	<0.001
Gender
Male	—	—
Female	0.11	-0.11, 0.32	0.3
Race
White	—	—
Black	0.86	0.43, 1.3	<0.001
Asian	-0.48	-1.1, 0.15	0.14
Native	-0.27	-2.0, 1.5	0.8
Other	0.28	-0.34, 0.90	0.4
Ethnicity
NonHispanic	—	—
Hispanic	1.2	0.44, 1.9	0.002
Primary Language
English	—	—
Other	-0.55	-1.6, 0.51	0.3
Charlson Comorbidity Index	0.25	0.23, 0.28	<0.001
Soceioeconomic Status	0.59	-0.01, 1.2	0.054
Household Composition	0.25	-0.28, 0.78	0.4
Minority Status and Language	0.01	-0.37, 0.40	>0.9
Housing and Transportation	-0.50	-0.94, -0.06	0.026
1 CI = Confidence Interval

# Model performance 
model_performance(steroid_count2)

# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------
46665.054 | 46783.484 | 0.060 |     0.058 | 4.742 | 4.747

performance::check_model(steroid_count2, panel = TRUE)

# Margins 
cplot(steroid_count2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Steroid Access Given Theme 1")

cplot(steroid_count2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of SSteroid Access Given THEME2")

cplot(steroid_count2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Steroide Access Given THEME3")

cplot(steroid_count2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Steroid Access Given THEME4")