Medication Prescriptions

Load Packages

library(tidyverse)
library(codebookr)
library(summarytools)
library(broom)
library(performance)
library(gt)
library(gtsummary)
library(janitor)
library(forcats)
library(here)
library(yaml)
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)
save(med_rx, file = "med_rx.rda")

Delete NA values for meds data

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

Make medication presribed variable (yes/no)

med_rx_drop %>%  
mutate(mesalamine_2 = case_when(MESALAMINE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(prednisone_2 = case_when(PREDNISONE>= 1 ~ '1',TRUE ~ "0")) %>%  
mutate(infliximab_2 = case_when(INFLIXIMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(vedolizumab_2 = case_when(VEDOLIZUMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(adalimumab_2 = case_when(ADALIMUMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(tofacitinib_2 = case_when(TOFACITINIB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(ustekinumab_2 = case_when(USTEKINUMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(methylpred_2 = case_when(METHYLPREDNISOLONE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(dexamethasone_2 = case_when(DEXAMETHASONE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(budesonide_2 = case_when(BUDESONIDE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(methotrexate_2 = case_when(METHOTREXATE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(azathioprine_2 = case_when(AZATHIOPRINE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(mercaptopurine_2 = case_when(MERCAPTOPURINE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(sulfasalazine_2 = case_when(SULFASALAZINE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(balsalazide_2 = case_when(BALSALAZIDE>= 1 ~ '1',TRUE ~ "0")) %>%   
mutate(cyclosporine_2 = case_when(CYCLOSPORINE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(certolizumab_2 = case_when(CERTOLIZUMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(prednisolone_2 = case_when(PREDNISOLONE>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(upadacitinib_2 = case_when(UPADACITINIB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(golimumab_2 = case_when(GOLIMUMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(natalizumab_2 = case_when(NATALIZUMAB>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(ozanimod_2 = case_when(OZANIMOD>= 1 ~ '1',TRUE ~ "0")) -> med_rx_2
View(med_rx_2)

Make numeric

med_rx_2$mesalamine_2 = as.numeric(med_rx_2$mesalamine_2)
class(med_rx_2$mesalamine_2)

[1] "numeric"

med_rx_2$prednisone_2 = as.numeric(med_rx_2$prednisone_2)
class(med_rx_2$prednisone_2)

[1] "numeric"

med_rx_2$infliximab_2 = as.numeric(med_rx_2$infliximab_2) 

med_rx_2$vedolizumab_2 = as.numeric(med_rx_2$vedolizumab_2)

med_rx_2$adalimumab_2 = as.numeric(med_rx_2$adalimumab_2)

med_rx_2$tofacitinib_2 = as.numeric(med_rx_2$tofacitinib_2)

med_rx_2$ustekinumab_2 = as.numeric(med_rx_2$ustekinumab_2)

med_rx_2$methylpred_2 = as.numeric(med_rx_2$methylpred_2)

med_rx_2$dexamethasone_2 = as.numeric(med_rx_2$dexamethasone_2)

med_rx_2$budesonide_2 = as.numeric(med_rx_2$budesonide_2)

med_rx_2$methotrexate_2 = as.numeric(med_rx_2$methotrexate_2)

med_rx_2$azathioprine_2 = as.numeric(med_rx_2$azathioprine_2)

med_rx_2$mercaptopurine_2 = as.numeric(med_rx_2$mercaptopurine_2)

med_rx_2$sulfasalazine_2 = as.numeric(med_rx_2$sulfasalazine_2)

med_rx_2$balsalazide_2 = as.numeric(med_rx_2$balsalazide_2)

med_rx_2$cyclosporine_2 = as.numeric(med_rx_2$cyclosporine_2)

med_rx_2$certolizumab_2 = as.numeric(med_rx_2$certolizumab_2)

med_rx_2$prednisolone_2 = as.numeric(med_rx_2$prednisolone_2)

med_rx_2$upadacitinib_2 = as.numeric(med_rx_2$upadacitinib_2)

med_rx_2$golimumab_2 = as.numeric(med_rx_2$golimumab_2)

med_rx_2$natalizumab_2 = as.numeric(med_rx_2$natalizumab_2)

med_rx_2$ozanimod_2 = as.numeric(med_rx_2$ozanimod_2)

Meds by class - count

med_rx_2 %>% 
  mutate(ASA = mesalamine_2 + sulfasalazine_2 + balsalazide_2) %>% 
  mutate(biologic = infliximab_2 + vedolizumab_2 + adalimumab_2 + ustekinumab_2 + certolizumab_2 +
           golimumab_2 + natalizumab_2) %>% 
  mutate(steroids = prednisone_2 + methylpred_2 + dexamethasone_2 + budesonide_2) %>% 
  mutate(immuno = methotrexate_2 + azathioprine_2 + mercaptopurine_2) %>% 
mutate(small = tofacitinib_2 + upadacitinib_2 + ozanimod_2) -> rx_class

View(rx_class)
  

Meds by class - yes/no

rx_class %>% 
mutate(ASA_2 = case_when(ASA>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(biologic_2 = case_when(biologic>= 1 ~ '1',TRUE ~ "0")) %>%  
mutate(steroids_2 = case_when(steroids>= 1 ~ '1',TRUE ~ "0")) %>% 
mutate(immuno_2 = case_when(immuno>= 1 ~ '1',TRUE ~ "0")) %>%
mutate(small_2 = case_when(small>= 1 ~ '1',TRUE ~ "0")) -> rx_class_2

View(rx_class_2)

Create UC only cohort

UC_meds <- rx_class_2[ which(rx_class_2$ibd_3=='UC'),]
View(UC_meds)

save(rx_class_2, file = "rx_class_2.rda")
save(UC_meds, file = "UC_meds.rda")

Baseline Characteristics

rx_class_2 %>% 
  dplyr::select(ibd_3, age_yrs, gender, race_5, ethnic_3, lang_3, max_ch, ASA_2, immuno_2, biologic_2, small_2, steroids_2, 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_2 ~ "5-ASA", immuno_2 ~ "Immunomodulator", biologic_2 ~ "Biologic", small_2 ~ "Small Molecule", steroids_2 ~ "Corticosteroids"),
        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)
5-ASA
0	4,889 (62%)
1	2,947 (38%)
Immunomodulator
0	5,626 (72%)
1	2,210 (28%)
Biologic
0	3,796 (48%)
1	4,040 (52%)
Small Molecule
0	7,661 (98%)
1	175 (2.2%)
Corticosteroids
0	3,307 (42%)
1	4,529 (58%)
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)

Baseline characteristics 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", ASA_2 ~ "5-ASA", immuno_2 ~ "Immunomodulator", biologic_2 ~ "Biologic", small_2 ~ "Small Molecule", steroids_2 ~ "Corticosteroids"),
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
5-ASA					<0.001
0	1,932 (60%)	1,469 (61%)	958 (65%)	530 (69%)
1	1,267 (40%)	923 (39%)	522 (35%)	235 (31%)
Immunomodulator					0.002
0	2,340 (73%)	1,735 (73%)	1,041 (70%)	510 (67%)
1	859 (27%)	657 (27%)	439 (30%)	255 (33%)
Biologic					0.5
0	1,569 (49%)	1,167 (49%)	706 (48%)	354 (46%)
1	1,630 (51%)	1,225 (51%)	774 (52%)	411 (54%)
Small Molecule					0.4
0	3,119 (97%)	2,339 (98%)	1,455 (98%)	748 (98%)
1	80 (2.5%)	53 (2.2%)	25 (1.7%)	17 (2.2%)
Corticosteroids					<0.001
0	1,413 (44%)	1,004 (42%)	622 (42%)	268 (35%)
1	1,786 (56%)	1,388 (58%)	858 (58%)	497 (65%)
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

baseline$ASA_2 = as.numeric(baseline$ASA_2)
class(baseline$ASA_2)

[1] "numeric"

baseline$immuno_2 = as.numeric(baseline$immuno_2)
class(baseline$immuno_2)

[1] "numeric"

baseline$biologic_2 = as.numeric(baseline$biologic_2)
class(baseline$biologic_2)

[1] "numeric"

baseline$small_2 = as.numeric(baseline$small_2)
class(baseline$small_2)

[1] "numeric"

baseline$steroids_2 = as.numeric(baseline$steroids_2)
class(baseline$steroids_2)

[1] "numeric"

5-ASA limited to UC patients

UC_meds$ASA_2 = as.numeric(UC_meds$ASA_2)
class(UC_meds$ASA_2)

[1] "numeric"

tbl_uv_asa<-
  tbl_uvregression(
    UC_meds[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_2")],
    method = glm,
    y = ASA_2,
    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, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	1.00	0.99, 1.00	0.004
Gender	3,784
Male		—	—
Female		1.11	0.98, 1.27	0.10
Race	3,784
White		—	—
Black		1.05	0.80, 1.39	0.7
Asian		2.36	1.64, 3.48	<0.001
Native		0.74	0.29, 1.90	0.5
Other		1.25	0.90, 1.75	0.2
Ethnicity	3,784
NonHispanic		—	—
Hispanic		1.24	0.84, 1.84	0.3
Primary Language	3,784
English		—	—
Other		1.62	0.89, 3.07	0.12
Charlson Comorbidity Index	3,784	0.94	0.92, 0.95	<0.001
Total SVI	3,784	0.62	0.48, 0.81	<0.001
SVI Quartiles	3,784
First		—	—
Second		0.90	0.77, 1.05	0.2
Third		0.84	0.70, 1.01	0.067
Fourth		0.71	0.55, 0.92	0.008
Soceioeconomic Status	3,784	0.53	0.41, 0.69	<0.001
Household Composition	3,784	0.45	0.35, 0.57	<0.001
Minority Status and Language	3,784	1.87	1.49, 2.34	<0.001
Housing and Transportation	3,784	0.80	0.64, 1.01	0.061
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Immunomodulators

tbl_uv_immuno<-
  tbl_uvregression(
    baseline[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "immuno_2")],
    method = glm,
    y = immuno_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_immuno, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		0.48	0.43, 0.53	<0.001
Unspecified		0.32	0.12, 0.71	0.011
Age	7,836	0.99	0.98, 0.99	<0.001
Gender	7,836
Male		—	—
Female		0.93	0.84, 1.03	0.2
Race	7,836
White		—	—
Black		0.94	0.77, 1.14	0.5
Asian		0.73	0.53, 0.99	0.049
Native		1.52	0.69, 3.18	0.3
Other		0.69	0.51, 0.93	0.015
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.01	0.72, 1.40	>0.9
Primary Language	7,836
English		—	—
Other		0.66	0.38, 1.11	0.13
Charlson Comorbidity Index	7,836	0.96	0.95, 0.98	<0.001
Total SVI	7,836	1.48	1.23, 1.80	<0.001
SVI Quartiles	7,836
First		—	—
Second		1.03	0.92, 1.16	0.6
Third		1.15	1.00, 1.32	0.046
Fourth		1.36	1.15, 1.61	<0.001
Soceioeconomic Status	7,836	1.62	1.34, 1.96	<0.001
Household Composition	7,836	1.46	1.21, 1.75	<0.001
Minority Status and Language	7,836	0.98	0.82, 1.16	0.8
Housing and Transportation	7,836	1.14	0.96, 1.36	0.12
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Biologics

tbl_uv_biologic<-
  tbl_uvregression(
    baseline[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "biologic_2")],
    method = glm,
    y = biologic_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_biologic, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		0.24	0.22, 0.27	<0.001
Unspecified		0.18	0.08, 0.35	<0.001
Age	7,836	0.97	0.97, 0.97	<0.001
Gender	7,836
Male		—	—
Female		0.78	0.72, 0.86	<0.001
Race	7,836
White		—	—
Black		1.41	1.19, 1.69	<0.001
Asian		0.71	0.55, 0.92	0.010
Native		0.78	0.37, 1.62	0.5
Other		1.15	0.90, 1.47	0.3
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.10	0.82, 1.49	0.5
Primary Language	7,836
English		—	—
Other		0.63	0.40, 0.98	0.041
Charlson Comorbidity Index	7,836	0.93	0.92, 0.94	<0.001
Total SVI	7,836	1.17	0.99, 1.40	0.070
SVI Quartiles	7,836
First		—	—
Second		1.01	0.91, 1.12	0.8
Third		1.06	0.93, 1.19	0.4
Fourth		1.12	0.95, 1.31	0.2
Soceioeconomic Status	7,836	1.20	1.01, 1.43	0.040
Household Composition	7,836	1.11	0.94, 1.31	0.2
Minority Status and Language	7,836	1.23	1.06, 1.44	0.008
Housing and Transportation	7,836	0.95	0.81, 1.11	0.5
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Infliximab

tbl_uv_infliximab<-
  tbl_uvregression(
    med_rx_2[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "infliximab_2")],
    method = glm,
    y = infliximab_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_infliximab, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		0.44	0.40, 0.49	<0.001
Unspecified		0.15	0.04, 0.42	0.002
Age	7,836	0.96	0.96, 0.97	<0.001
Gender	7,836
Male		—	—
Female		0.72	0.66, 0.80	<0.001
Race	7,836
White		—	—
Black		1.41	1.17, 1.70	<0.001
Asian		0.92	0.68, 1.24	0.6
Native		0.73	0.27, 1.68	0.5
Other		1.18	0.90, 1.54	0.2
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.04	0.74, 1.45	0.8
Primary Language	7,836
English		—	—
Other		0.87	0.51, 1.42	0.6
Charlson Comorbidity Index	7,836	0.93	0.92, 0.94	<0.001
Total SVI	7,836	0.92	0.76, 1.12	0.4
SVI Quartiles	7,836
First		—	—
Second		0.88	0.78, 0.99	0.036
Third		0.94	0.82, 1.08	0.4
Fourth		0.95	0.79, 1.13	0.6
Soceioeconomic Status	7,836	0.96	0.79, 1.17	0.7
Household Composition	7,836	0.92	0.76, 1.11	0.4
Minority Status and Language	7,836	1.22	1.03, 1.45	0.023
Housing and Transportation	7,836	0.78	0.66, 0.93	0.006
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Adalimumab

tbl_uv_ada<-
  tbl_uvregression(
    med_rx_2[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "adalimumab_2")],
    method = glm,
    y = adalimumab_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_ada, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		0.35	0.31, 0.40	<0.001
Unspecified		0.67	0.27, 1.44	0.3
Age	7,836	0.99	0.99, 0.99	<0.001
Gender	7,836
Male		—	—
Female		0.96	0.85, 1.08	0.5
Race	7,836
White		—	—
Black		1.26	1.01, 1.56	0.034
Asian		0.67	0.44, 0.98	0.049
Native		0.77	0.23, 1.99	0.6
Other		0.90	0.63, 1.24	0.5
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.03	0.68, 1.50	0.9
Primary Language	7,836
English		—	—
Other		0.67	0.32, 1.23	0.2
Charlson Comorbidity Index	7,836	0.94	0.93, 0.96	<0.001
Total SVI	7,836	1.48	1.18, 1.85	<0.001
SVI Quartiles	7,836
First		—	—
Second		1.25	1.08, 1.44	0.002
Third		1.25	1.06, 1.47	0.008
Fourth		1.27	1.03, 1.55	0.024
Soceioeconomic Status	7,836	1.41	1.12, 1.76	0.003
Household Composition	7,836	1.31	1.05, 1.63	0.017
Minority Status and Language	7,836	1.26	1.03, 1.54	0.028
Housing and Transportation	7,836	1.18	0.96, 1.45	0.11
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Vedolizumab

tbl_uv_vedo<-
  tbl_uvregression(
    med_rx_2[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "vedolizumab_2")],
    method = glm,
    y = vedolizumab_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_vedo, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		1.00	0.86, 1.16	>0.9
Unspecified		0.00	0.00, 0.05	>0.9
Age	7,836	1.00	0.99, 1.00	0.2
Gender	7,836
Male		—	—
Female		1.03	0.88, 1.20	0.7
Race	7,836
White		—	—
Black		1.00	0.73, 1.33	>0.9
Asian		0.94	0.58, 1.44	0.8
Native		1.11	0.27, 3.17	0.9
Other		1.59	1.09, 2.24	0.011
Ethnicity	7,836
NonHispanic		—	—
Hispanic		1.44	0.90, 2.20	0.11
Primary Language	7,836
English		—	—
Other		1.17	0.54, 2.23	0.7
Charlson Comorbidity Index	7,836	1.03	1.02, 1.05	<0.001
Total SVI	7,836	0.69	0.51, 0.93	0.016
SVI Quartiles	7,836
First		—	—
Second		0.99	0.83, 1.18	>0.9
Third		0.87	0.70, 1.07	0.2
Fourth		0.72	0.53, 0.96	0.029
Soceioeconomic Status	7,836	0.76	0.56, 1.02	0.073
Household Composition	7,836	0.61	0.45, 0.82	<0.001
Minority Status and Language	7,836	1.23	0.95, 1.60	0.12
Housing and Transportation	7,836	0.73	0.56, 0.95	0.020
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Ustekinumab

tbl_uv_uste<-
  tbl_uvregression(
    med_rx_2[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "ustekinumab_2")],
    method = glm,
    y = ustekinumab_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_uste, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		0.19	0.16, 0.22	<0.001
Unspecified		0.10	0.01, 0.48	0.026
Age	7,836	0.99	0.99, 0.99	<0.001
Gender	7,836
Male		—	—
Female		1.17	1.02, 1.34	0.023
Race	7,836
White		—	—
Black		1.29	1.01, 1.63	0.040
Asian		0.52	0.30, 0.83	0.011
Native		0.53	0.09, 1.76	0.4
Other		0.65	0.41, 0.99	0.059
Ethnicity	7,836
NonHispanic		—	—
Hispanic		0.57	0.31, 0.96	0.051
Primary Language	7,836
English		—	—
Other		0.67	0.28, 1.36	0.3
Charlson Comorbidity Index	7,836	0.98	0.96, 0.99	0.002
Total SVI	7,836	1.41	1.09, 1.83	0.009
SVI Quartiles	7,836
First		—	—
Second		0.98	0.83, 1.15	0.8
Third		1.13	0.93, 1.35	0.2
Fourth		1.27	1.01, 1.60	0.037
Soceioeconomic Status	7,836	1.51	1.17, 1.96	0.002
Household Composition	7,836	1.42	1.10, 1.82	0.007
Minority Status and Language	7,836	1.03	0.81, 1.30	0.8
Housing and Transportation	7,836	1.06	0.84, 1.35	0.6
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Small molecules

UC_meds$small_2 = as.numeric(UC_meds$small_2)
class(UC_meds$small_2)

[1] "numeric"

tbl_uv_small<-
  tbl_uvregression(
    UC_meds[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_2")],
    method = glm,
    y = small_2,
    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, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
Age	3,784	0.99	0.98, 1.00	0.062
Gender	3,784
Male		—	—
Female		0.59	0.42, 0.82	0.002
Race	3,784
White		—	—
Black		0.68	0.26, 1.42	0.4
Asian		0.64	0.19, 1.54	0.4
Native		0.00		>0.9
Other		1.30	0.58, 2.54	0.5
Ethnicity	3,784
NonHispanic		—	—
Hispanic		0.44	0.07, 1.41	0.3
Primary Language	3,784
English		—	—
Other		0.51	0.03, 2.36	0.5
Charlson Comorbidity Index	3,784	0.96	0.92, 1.00	0.036
Total SVI	3,784	0.64	0.31, 1.28	0.2
SVI Quartiles	3,784
First		—	—
Second		0.83	0.56, 1.21	0.3
Third		0.77	0.46, 1.23	0.3
Fourth		0.72	0.33, 1.38	0.4
Soceioeconomic Status	3,784	0.57	0.28, 1.14	0.12
Household Composition	3,784	1.13	0.59, 2.12	0.7
Minority Status and Language	3,784	0.79	0.45, 1.40	0.4
Housing and Transportation	3,784	0.76	0.41, 1.36	0.4
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Steroids

tbl_uv_steroids<-
  tbl_uvregression(
    baseline[c("ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "max_ch", "RPL_THEMES", "RPL_4", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4", "steroids_2")],
    method = glm,
    y = steroids_2,
    method.args = list(family = binomial),
    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"),
    exponentiate = TRUE
  )
print(tbl_uv_steroids, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	7,836
CD		—	—
UC		0.99	0.91, 1.09	0.9
Unspecified		0.73	0.39, 1.36	0.3
Age	7,836	1.01	1.00, 1.01	<0.001
Gender	7,836
Male		—	—
Female		1.23	1.12, 1.34	<0.001
Race	7,836
White		—	—
Black		1.29	1.08, 1.55	0.005
Asian		0.61	0.47, 0.79	<0.001
Native		1.03	0.49, 2.21	>0.9
Other		0.82	0.64, 1.05	0.12
Ethnicity	7,836
NonHispanic		—	—
Hispanic		0.96	0.71, 1.30	0.8
Primary Language	7,836
English		—	—
Other		0.60	0.38, 0.92	0.021
Charlson Comorbidity Index	7,836	1.10	1.09, 1.11	<0.001
Total SVI	7,836	1.46	1.23, 1.74	<0.001
SVI Quartiles	7,836
First		—	—
Second		1.09	0.98, 1.22	0.10
Third		1.09	0.96, 1.24	0.2
Fourth		1.47	1.25, 1.73	<0.001
Soceioeconomic Status	7,836	1.56	1.31, 1.86	<0.001
Household Composition	7,836	1.55	1.31, 1.84	<0.001
Minority Status and Language	7,836	0.80	0.68, 0.93	0.004
Housing and Transportation	7,836	1.20	1.03, 1.41	0.021
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Multivariable Analysis with RPL_THEMES

Mesalamine Limited to UC patients

mesalamine_rx <- glm(ASA_2 ~  age_yrs + gender + race_5 + ethnic_3 + lang_3 
                     + max_ch + RPL_THEMES,
              family = "binomial",
              data = UC_meds)
summary(mesalamine_rx )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8853  -1.2976   0.8957   0.9973   1.7327  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.428249   0.112030   3.823 0.000132 ***
age_yrs           0.005396   0.002012   2.682 0.007329 ** 
genderFemale      0.090846   0.067877   1.338 0.180769    
race_5Black       0.184954   0.148659   1.244 0.213443    
race_5Asian       0.732002   0.195427   3.746 0.000180 ***
race_5Native     -0.346942   0.478125  -0.726 0.468065    
race_5Other       0.136301   0.184547   0.739 0.460167    
ethnic_3Hispanic  0.147884   0.216588   0.683 0.494739    
lang_3Other       0.215863   0.328184   0.658 0.510699    
max_ch           -0.073720   0.007417  -9.939  < 2e-16 ***
RPL_THEMES       -0.510102   0.142516  -3.579 0.000345 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5135.6  on 3783  degrees of freedom
Residual deviance: 4983.5  on 3773  degrees of freedom
AIC: 5005.5

Number of Fisher Scoring iterations: 4

broom::glance(mesalamine_rx )

broom::tidy(mesalamine_rx , exponentiate = TRUE)

tbl_regression(mesalamine_rx, 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	1.01	1.00, 1.01	0.007
Gender
Male	—	—
Female	1.10	0.96, 1.25	0.2
Race
White	—	—
Black	1.20	0.90, 1.61	0.2
Asian	2.08	1.43, 3.09	<0.001
Native	0.71	0.27, 1.83	0.5
Other	1.15	0.80, 1.65	0.5
Ethnicity
NonHispanic	—	—
Hispanic	1.16	0.76, 1.78	0.5
Primary Language
English	—	—
Other	1.24	0.66, 2.42	0.5
Charlson Comorbidity Index	0.93	0.92, 0.94	<0.001
Total SVI	0.60	0.45, 0.79	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(mesalamine_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log |   PCP
------------------------------------------------------------------------------
5005.480 | 5074.104 |     0.040 | 0.483 | 1.149 |    0.658 |      -Inf | 0.534

performance::check_model(mesalamine_rx, panel = TRUE)

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

Immunomodulators

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

Call:
glm(formula = immuno_2 ~ 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.2656  -0.8500  -0.6844   1.3160   2.1605  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -0.165390   0.081804  -2.022  0.04320 *  
ibd_3UC          -0.671254   0.052915 -12.685  < 2e-16 ***
ibd_3Unspecified -1.137523   0.446188  -2.549  0.01079 *  
age_yrs          -0.011312   0.001579  -7.165 7.77e-13 ***
genderFemale     -0.037702   0.051638  -0.730  0.46531    
race_5Black      -0.248140   0.104272  -2.380  0.01733 *  
race_5Asian      -0.257248   0.162348  -1.585  0.11307    
race_5Native      0.486987   0.390056   1.249  0.21185    
race_5Other      -0.409378   0.162112  -2.525  0.01156 *  
ethnic_3Hispanic  0.147020   0.182114   0.807  0.41950    
lang_3Other      -0.239689   0.284557  -0.842  0.39961    
max_ch           -0.011465   0.006284  -1.824  0.06809 .  
RPL_THEMES        0.326741   0.102755   3.180  0.00147 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9322.7  on 7835  degrees of freedom
Residual deviance: 8995.9  on 7823  degrees of freedom
AIC: 9021.9

Number of Fisher Scoring iterations: 4

broom::glance(immuno_rx )

broom::tidy(immuno_rx , exponentiate = TRUE)

tbl_regression(immuno_rx, 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.51	0.46, 0.57	<0.001
Unspecified	0.32	0.12, 0.72	0.011
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	0.96	0.87, 1.07	0.5
Race
White	—	—
Black	0.78	0.63, 0.96	0.017
Asian	0.77	0.56, 1.06	0.11
Native	1.63	0.73, 3.45	0.2
Other	0.66	0.48, 0.91	0.012
Ethnicity
NonHispanic	—	—
Hispanic	1.16	0.81, 1.65	0.4
Primary Language
English	—	—
Other	0.79	0.44, 1.34	0.4
Charlson Comorbidity Index	0.99	0.98, 1.00	0.068
Total SVI	1.39	1.13, 1.70	0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(immuno_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
9021.907 | 9112.471 |     0.041 | 0.441 | 1.072 |    0.574 |      -Inf |       1.276e-04 | 0.612

performance::check_model(immuno_rx, panel = TRUE)

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

Biologics

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

Call:
glm(formula = biologic_2 ~ 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  
-2.0403  -1.0047   0.5697   0.9587   2.1631  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.234195   0.087984  25.393  < 2e-16 ***
ibd_3UC          -1.377142   0.050702 -27.161  < 2e-16 ***
ibd_3Unspecified -1.768512   0.365780  -4.835 1.33e-06 ***
age_yrs          -0.027214   0.001543 -17.637  < 2e-16 ***
genderFemale     -0.211840   0.050407  -4.203 2.64e-05 ***
race_5Black       0.236898   0.101743   2.328   0.0199 *  
race_5Asian      -0.378743   0.148013  -2.559   0.0105 *  
race_5Native     -0.118172   0.398078  -0.297   0.7666    
race_5Other       0.135652   0.146829   0.924   0.3555    
ethnic_3Hispanic  0.114051   0.177789   0.641   0.5212    
lang_3Other      -0.237877   0.253556  -0.938   0.3482    
max_ch           -0.024948   0.005834  -4.277 1.90e-05 ***
RPL_THEMES       -0.209637   0.101487  -2.066   0.0389 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 10855.4  on 7835  degrees of freedom
Residual deviance:  9331.5  on 7823  degrees of freedom
AIC: 9357.5

Number of Fisher Scoring iterations: 4

broom::glance(biologic_rx )

broom::tidy(biologic_rx , exponentiate = TRUE)

tbl_regression(biologic_rx, 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.25	0.23, 0.28	<0.001
Unspecified	0.17	0.08, 0.34	<0.001
Age	0.97	0.97, 0.98	<0.001
Gender
Male	—	—
Female	0.81	0.73, 0.89	<0.001
Race
White	—	—
Black	1.27	1.04, 1.55	0.020
Asian	0.68	0.51, 0.91	0.011
Native	0.89	0.40, 1.94	0.8
Other	1.15	0.86, 1.53	0.4
Ethnicity
NonHispanic	—	—
Hispanic	1.12	0.79, 1.59	0.5
Primary Language
English	—	—
Other	0.79	0.48, 1.29	0.3
Charlson Comorbidity Index	0.98	0.96, 0.99	<0.001
Total SVI	0.81	0.66, 0.99	0.039
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(biologic_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
9357.459 | 9448.023 |     0.182 | 0.452 | 1.092 |    0.595 |      -Inf |       1.276e-04 | 0.591

performance::check_model(biologic_rx, panel = TRUE)

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

Infliximab

infliximab_rx <- glm(infliximab_2 ~  ibd_3 + age_yrs + gender + race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEMES,
              family = "binomial",
              data = med_rx_2)
summary(infliximab_rx )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5020  -0.8092  -0.5677   1.0415   2.6500  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       1.0754985  0.0869036  12.376  < 2e-16 ***
ibd_3UC          -0.7226124  0.0558051 -12.949  < 2e-16 ***
ibd_3Unspecified -1.8481913  0.6078547  -3.041 0.002362 ** 
age_yrs          -0.0358048  0.0017789 -20.128  < 2e-16 ***
genderFemale     -0.2507976  0.0542993  -4.619 3.86e-06 ***
race_5Black       0.3112473  0.1033024   3.013 0.002587 ** 
race_5Asian      -0.2187263  0.1634072  -1.339 0.180722    
race_5Native     -0.1884973  0.4760557  -0.396 0.692137    
race_5Other       0.0791040  0.1537462   0.515 0.606895    
ethnic_3Hispanic -0.0828051  0.1890478  -0.438 0.661378    
lang_3Other       0.1452288  0.2826756   0.514 0.607416    
max_ch           -0.0005387  0.0070684  -0.076 0.939247    
RPL_THEMES       -0.3611875  0.1096050  -3.295 0.000983 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9141.4  on 7835  degrees of freedom
Residual deviance: 8227.7  on 7823  degrees of freedom
AIC: 8253.7

Number of Fisher Scoring iterations: 4

broom::glance(infliximab_rx )

broom::tidy(infliximab_rx , exponentiate = TRUE)

tbl_regression(infliximab_rx, 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.49	0.44, 0.54	<0.001
Unspecified	0.16	0.04, 0.45	0.002
Age	0.96	0.96, 0.97	<0.001
Gender
Male	—	—
Female	0.78	0.70, 0.87	<0.001
Race
White	—	—
Black	1.37	1.11, 1.67	0.003
Asian	0.80	0.58, 1.10	0.2
Native	0.83	0.30, 1.99	0.7
Other	1.08	0.80, 1.46	0.6
Ethnicity
NonHispanic	—	—
Hispanic	0.92	0.63, 1.33	0.7
Primary Language
English	—	—
Other	1.16	0.65, 1.98	0.6
Charlson Comorbidity Index	1.00	0.99, 1.01	>0.9
Total SVI	0.70	0.56, 0.86	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(infliximab_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8253.689 | 8344.254 |     0.118 | 0.416 | 1.026 |    0.525 |      -Inf |       2.530e-04 | 0.652

performance::check_model(infliximab_rx, panel = TRUE)

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

Adalimumab

ada_rx <- glm(adalimumab_2 ~  ibd_3 + age_yrs + gender + race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEMES,
              family = "binomial",
              data = med_rx_2)
summary(ada_rx )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8745  -0.7441  -0.4942  -0.3953   2.4445  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.0701423  0.0963809 -11.103  < 2e-16 ***
ibd_3UC          -0.9982753  0.0665945 -14.990  < 2e-16 ***
ibd_3Unspecified -0.3967320  0.4196983  -0.945   0.3445    
age_yrs          -0.0002828  0.0018672  -0.151   0.8796    
genderFemale     -0.0300197  0.0616693  -0.487   0.6264    
race_5Black       0.0957953  0.1162978   0.824   0.4101    
race_5Asian      -0.2725446  0.2069800  -1.317   0.1879    
race_5Native     -0.1381513  0.5477863  -0.252   0.8009    
race_5Other      -0.0836809  0.1852752  -0.452   0.6515    
ethnic_3Hispanic  0.1604023  0.2151993   0.745   0.4561    
lang_3Other      -0.3079545  0.3507759  -0.878   0.3800    
max_ch           -0.0524628  0.0083435  -6.288 3.22e-10 ***
RPL_THEMES        0.2061805  0.1214483   1.698   0.0896 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7194.8  on 7835  degrees of freedom
Residual deviance: 6855.5  on 7823  degrees of freedom
AIC: 6881.5

Number of Fisher Scoring iterations: 5

broom::glance(ada_rx )

broom::tidy(ada_rx , exponentiate = TRUE)

tbl_regression(ada_rx, 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.37	0.32, 0.42	<0.001
Unspecified	0.67	0.27, 1.44	0.3
Age	1.00	1.00, 1.00	0.9
Gender
Male	—	—
Female	0.97	0.86, 1.10	0.6
Race
White	—	—
Black	1.10	0.87, 1.38	0.4
Asian	0.76	0.50, 1.12	0.2
Native	0.87	0.25, 2.29	0.8
Other	0.92	0.63, 1.31	0.7
Ethnicity
NonHispanic	—	—
Hispanic	1.17	0.76, 1.77	0.5
Primary Language
English	—	—
Other	0.73	0.35, 1.40	0.4
Charlson Comorbidity Index	0.95	0.93, 0.96	<0.001
Total SVI	1.23	0.97, 1.56	0.090
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ada_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
6881.500 | 6972.065 |     0.042 | 0.370 | 0.936 |    0.437 |  -259.442 |       1.906e-04 | 0.727

performance::check_model(ada_rx, panel = TRUE)

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

Vedolizumab

vedo_rx <- glm(vedolizumab_2 ~  ibd_3 + age_yrs + gender + race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEMES,
              family = "binomial",
              data = med_rx_2)
summary(vedo_rx )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8089  -0.4668  -0.4297  -0.3932   2.4687  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -1.898003   0.124262 -15.274  < 2e-16 ***
ibd_3UC           -0.021442   0.078563  -0.273  0.78490    
ibd_3Unspecified -13.275341 228.869785  -0.058  0.95375    
age_yrs           -0.009499   0.002374  -4.001 6.31e-05 ***
genderFemale       0.059036   0.077869   0.758  0.44836    
race_5Black        0.074423   0.157383   0.473  0.63630    
race_5Asian       -0.059343   0.237333  -0.250  0.80255    
race_5Native       0.114435   0.613432   0.187  0.85201    
race_5Other        0.425644   0.196141   2.170  0.03000 *  
ethnic_3Hispanic   0.221249   0.244324   0.906  0.36517    
lang_3Other        0.194876   0.369305   0.528  0.59772    
max_ch             0.050666   0.008144   6.221 4.94e-10 ***
RPL_THEMES        -0.433789   0.160511  -2.703  0.00688 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4931.9  on 7835  degrees of freedom
Residual deviance: 4872.7  on 7823  degrees of freedom
AIC: 4898.7

Number of Fisher Scoring iterations: 14

broom::glance(vedo_rx )

broom::tidy(vedo_rx , exponentiate = TRUE)

tbl_regression(vedo_rx, 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.98	0.84, 1.14	0.8
Unspecified	0.00	0.00, 0.05	>0.9
Age	0.99	0.99, 1.00	<0.001
Gender
Male	—	—
Female	1.06	0.91, 1.24	0.4
Race
White	—	—
Black	1.08	0.78, 1.45	0.6
Asian	0.94	0.58, 1.47	0.8
Native	1.12	0.27, 3.21	0.9
Other	1.53	1.03, 2.22	0.030
Ethnicity
NonHispanic	—	—
Hispanic	1.25	0.76, 1.98	0.4
Primary Language
English	—	—
Other	1.22	0.55, 2.38	0.6
Charlson Comorbidity Index	1.05	1.04, 1.07	<0.001
Total SVI	0.65	0.47, 0.89	0.007
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(vedo_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
4898.719 | 4989.283 |     0.007 | 0.293 | 0.789 |    0.311 |   -74.722 |       8.071e-04 | 0.829

performance::check_model(vedo_rx, panel = TRUE)

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

Ustekinumab

uste_rx <- glm(ustekinumab_2 ~  ibd_3 + age_yrs + gender + race_5 + ethnic_3 + lang_3 + max_ch + RPL_THEMES,
              family = "binomial",
              data = med_rx_2)
summary(uste_rx )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7837  -0.6511  -0.3221  -0.2815   2.7254  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.303528   0.110871 -11.757  < 2e-16 ***
ibd_3UC          -1.641980   0.089797 -18.286  < 2e-16 ***
ibd_3Unspecified -2.251208   1.013849  -2.220  0.02639 *  
age_yrs          -0.004644   0.002194  -2.116  0.03430 *  
genderFemale      0.193877   0.072194   2.685  0.00724 ** 
race_5Black       0.094076   0.132245   0.711  0.47685    
race_5Asian      -0.447005   0.271876  -1.644  0.10015    
race_5Native     -0.488673   0.750917  -0.651  0.51520    
race_5Other      -0.231723   0.238942  -0.970  0.33215    
ethnic_3Hispanic -0.306426   0.306306  -1.000  0.31712    
lang_3Other      -0.041524   0.416924  -0.100  0.92067    
max_ch           -0.007596   0.008814  -0.862  0.38883    
RPL_THEMES        0.056831   0.140893   0.403  0.68668    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5812.6  on 7835  degrees of freedom
Residual deviance: 5323.8  on 7823  degrees of freedom
AIC: 5349.8

Number of Fisher Scoring iterations: 6

broom::glance(uste_rx )

broom::tidy(uste_rx , exponentiate = TRUE)

tbl_regression(uste_rx, 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.19	0.16, 0.23	<0.001
Unspecified	0.11	0.01, 0.49	0.026
Age	1.00	0.99, 1.00	0.034
Gender
Male	—	—
Female	1.21	1.05, 1.40	0.007
Race
White	—	—
Black	1.10	0.84, 1.42	0.5
Asian	0.64	0.36, 1.06	0.10
Native	0.61	0.10, 2.15	0.5
Other	0.79	0.48, 1.24	0.3
Ethnicity
NonHispanic	—	—
Hispanic	0.74	0.39, 1.29	0.3
Primary Language
English	—	—
Other	0.96	0.39, 2.04	>0.9
Charlson Comorbidity Index	0.99	0.98, 1.01	0.4
Total SVI	1.06	0.80, 1.39	0.7
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(uste_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
5349.790 | 5440.355 |     0.059 | 0.318 | 0.825 |    0.340 |  -127.674 |       3.025e-04 | 0.798

performance::check_model(uste_rx, panel = TRUE)

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

Small molecules limited to UC patients only

small_rx <- glm(small_2 ~  age_yrs + gender + race_5 + ethnic_3 + lang_3 
                     + max_ch + RPL_THEMES,
              family = "binomial",
              data = UC_meds)
summary(small_rx )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.4594  -0.3168  -0.2678  -0.2338   2.9530  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -2.467573   0.265755  -9.285  < 2e-16 ***
age_yrs           -0.004787   0.005022  -0.953  0.34053    
genderFemale      -0.516495   0.171773  -3.007  0.00264 ** 
race_5Black       -0.258714   0.433124  -0.597  0.55029    
race_5Asian       -0.539520   0.519949  -1.038  0.29944    
race_5Native     -12.286357 338.759401  -0.036  0.97107    
race_5Other        0.381766   0.392848   0.972  0.33116    
ethnic_3Hispanic  -0.996527   0.747999  -1.332  0.18278    
lang_3Other       -0.454824   1.037457  -0.438  0.66109    
max_ch            -0.032679   0.021027  -1.554  0.12014    
RPL_THEMES        -0.391026   0.368234  -1.062  0.28828    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1243.2  on 3783  degrees of freedom
Residual deviance: 1219.7  on 3773  degrees of freedom
AIC: 1241.7

Number of Fisher Scoring iterations: 14

broom::glance(small_rx )

broom::tidy(small_rx , exponentiate = TRUE)

tbl_regression(small_rx, 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	1.00	0.99, 1.01	0.3
Gender
Male	—	—
Female	0.60	0.42, 0.83	0.003
Race
White	—	—
Black	0.77	0.30, 1.66	0.6
Asian	0.58	0.18, 1.42	0.3
Native	0.00		>0.9
Other	1.46	0.63, 2.98	0.3
Ethnicity
NonHispanic	—	—
Hispanic	0.37	0.06, 1.27	0.2
Primary Language
English	—	—
Other	0.63	0.03, 3.15	0.7
Charlson Comorbidity Index	0.97	0.93, 1.01	0.12
Total SVI	0.68	0.32, 1.38	0.3
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(small_rx)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1241.731 | 1310.355 |     0.006 | 0.193 | 0.569 |    0.161 |    -5.770 |           0.004 | 0.926

performance::check_model(small_rx, panel = TRUE)

# Margins 
cplot(small_rx, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Small Molecule Rx Given SVI")

Steroids

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

Call:
glm(formula = steroids_2 ~ 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  
-2.2118  -1.2069   0.7211   1.0913   1.4961  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.109851   0.076117   1.443 0.148966    
ibd_3UC          -0.004122   0.047891  -0.086 0.931406    
ibd_3Unspecified -0.312056   0.327583  -0.953 0.340792    
age_yrs          -0.007646   0.001427  -5.360 8.34e-08 ***
genderFemale      0.217705   0.047318   4.601 4.21e-06 ***
race_5Black       0.163414   0.096696   1.690 0.091032 .  
race_5Asian      -0.321734   0.136980  -2.349 0.018835 *  
race_5Native     -0.012846   0.385404  -0.033 0.973411    
race_5Other      -0.112555   0.135982  -0.828 0.407830    
ethnic_3Hispanic  0.101630   0.165234   0.615 0.538511    
lang_3Other      -0.320614   0.233509  -1.373 0.169745    
max_ch            0.109799   0.006368  17.241  < 2e-16 ***
RPL_THEMES        0.326855   0.095571   3.420 0.000626 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 10672  on 7835  degrees of freedom
Residual deviance: 10244  on 7823  degrees of freedom
AIC: 10270

Number of Fisher Scoring iterations: 4

broom::glance(steroids_rx )

broom::tidy(steroids_rx , exponentiate = TRUE)

tbl_regression(steroids_rx, 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.91, 1.09	>0.9
Unspecified	0.73	0.38, 1.40	0.3
Age	0.99	0.99, 1.00	<0.001
Gender
Male	—	—
Female	1.24	1.13, 1.36	<0.001
Race
White	—	—
Black	1.18	0.98, 1.42	0.091
Asian	0.72	0.55, 0.95	0.019
Native	0.99	0.47, 2.15	>0.9
Other	0.89	0.68, 1.17	0.4
Ethnicity
NonHispanic	—	—
Hispanic	1.11	0.80, 1.53	0.5
Primary Language
English	—	—
Other	0.73	0.46, 1.15	0.2
Charlson Comorbidity Index	1.12	1.10, 1.13	<0.001
Total SVI	1.39	1.15, 1.67	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(steroids_rx)

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log |   PCP
--------------------------------------------------------------------------------
10270.410 | 10360.974 |     0.051 | 0.481 | 1.144 |    0.654 |      -Inf | 0.537

performance::check_model(steroids_rx, panel = TRUE)

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

Multivariable Analysis with all RPL_Themes

Mesalamine (linited to UC patients only)

mesalamine_rx2 <- glm(ASA_2 ~  age_yrs + gender + race_5 + ethnic_3 + lang_3 
                     + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = UC_meds)
summary(mesalamine_rx2 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8639  -1.2590   0.8501   1.0000   1.8527  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.285543   0.133923   2.132 0.032996 *  
age_yrs           0.005968   0.002025   2.948 0.003203 ** 
genderFemale      0.085650   0.068193   1.256 0.209118    
race_5Black       0.087034   0.152011   0.573 0.566949    
race_5Asian       0.516324   0.199233   2.592 0.009554 ** 
race_5Native     -0.325855   0.482127  -0.676 0.499123    
race_5Other       0.063642   0.185888   0.342 0.732075    
ethnic_3Hispanic  0.142460   0.217695   0.654 0.512853    
lang_3Other       0.099579   0.329558   0.302 0.762530    
max_ch           -0.073121   0.007452  -9.813  < 2e-16 ***
RPL_THEME1       -0.425063   0.193172  -2.200 0.027776 *  
RPL_THEME2       -0.430123   0.169483  -2.538 0.011154 *  
RPL_THEME3        0.471579   0.125033   3.772 0.000162 ***
RPL_THEME4        0.051994   0.142250   0.366 0.714730    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5135.6  on 3783  degrees of freedom
Residual deviance: 4949.7  on 3770  degrees of freedom
AIC: 4977.7

Number of Fisher Scoring iterations: 4

broom::glance(mesalamine_rx2 )

broom::tidy(mesalamine_rx2 , exponentiate = TRUE)

tbl_regression(mesalamine_rx2, 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	1.01	1.00, 1.01	0.003
Gender
Male	—	—
Female	1.09	0.95, 1.25	0.2
Race
White	—	—
Black	1.09	0.81, 1.47	0.6
Asian	1.68	1.14, 2.51	0.010
Native	0.72	0.28, 1.89	0.5
Other	1.07	0.74, 1.54	0.7
Ethnicity
NonHispanic	—	—
Hispanic	1.15	0.76, 1.78	0.5
Primary Language
English	—	—
Other	1.10	0.59, 2.16	0.8
Charlson Comorbidity Index	0.93	0.92, 0.94	<0.001
Soceioeconomic Status	0.65	0.45, 0.95	0.028
Household Composition	0.65	0.47, 0.91	0.011
Minority Status and Language	1.60	1.25, 2.05	<0.001
Housing and Transportation	1.05	0.80, 1.39	0.7
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(mesalamine_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log |   PCP
------------------------------------------------------------------------------
4977.653 | 5064.993 |     0.049 | 0.481 | 1.146 |    0.654 |      -Inf | 0.538

performance::check_model(mesalamine_rx2, panel = TRUE)

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

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

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

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

Immunomodulators

immuno_rx2 <- glm(immuno_2 ~  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_rx2 )

Call:
glm(formula = immuno_2 ~ 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.2288  -0.8516  -0.6837   1.3151   2.1681  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -0.165606   0.097976  -1.690  0.09097 .  
ibd_3UC          -0.666095   0.052970 -12.575  < 2e-16 ***
ibd_3Unspecified -1.147217   0.446291  -2.571  0.01015 *  
age_yrs          -0.011379   0.001581  -7.196 6.22e-13 ***
genderFemale     -0.037675   0.051669  -0.729  0.46591    
race_5Black      -0.259139   0.105923  -2.446  0.01443 *  
race_5Asian      -0.210733   0.164730  -1.279  0.20080    
race_5Native      0.486267   0.390944   1.244  0.21356    
race_5Other      -0.405689   0.162522  -2.496  0.01255 *  
ethnic_3Hispanic  0.157174   0.182154   0.863  0.38821    
lang_3Other      -0.227904   0.285465  -0.798  0.42466    
max_ch           -0.011545   0.006289  -1.836  0.06640 .  
RPL_THEME1        0.415114   0.145530   2.852  0.00434 ** 
RPL_THEME2        0.131502   0.128313   1.025  0.30543    
RPL_THEME3       -0.013149   0.094282  -0.139  0.88908    
RPL_THEME4       -0.156370   0.107524  -1.454  0.14587    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9322.7  on 7835  degrees of freedom
Residual deviance: 8987.2  on 7820  degrees of freedom
AIC: 9019.2

Number of Fisher Scoring iterations: 4

broom::glance(immuno_rx2 )

broom::tidy(immuno_rx2 , exponentiate = TRUE)

tbl_regression(immuno_rx2, 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.51	0.46, 0.57	<0.001
Unspecified	0.32	0.12, 0.71	0.010
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	0.96	0.87, 1.07	0.5
Race
White	—	—
Black	0.77	0.63, 0.95	0.014
Asian	0.81	0.58, 1.11	0.2
Native	1.63	0.73, 3.45	0.2
Other	0.67	0.48, 0.91	0.013
Ethnicity
NonHispanic	—	—
Hispanic	1.17	0.81, 1.66	0.4
Primary Language
English	—	—
Other	0.80	0.44, 1.36	0.4
Charlson Comorbidity Index	0.99	0.98, 1.00	0.066
Soceioeconomic Status	1.51	1.14, 2.01	0.004
Household Composition	1.14	0.89, 1.47	0.3
Minority Status and Language	0.99	0.82, 1.19	0.9
Housing and Transportation	0.86	0.69, 1.06	0.15
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(immuno_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
9019.220 | 9130.683 |     0.042 | 0.440 | 1.072 |    0.573 |      -Inf |       1.276e-04 | 0.612

performance::check_model(immuno_rx2, panel = TRUE)

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

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

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

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

Biologics

biologic_rx2 <- glm(biologic_2 ~  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_rx2 )

Call:
glm(formula = biologic_2 ~ 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  
-2.0386  -1.0013   0.5658   0.9579   2.1592  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.215922   0.102812  21.553  < 2e-16 ***
ibd_3UC          -1.376958   0.050766 -27.124  < 2e-16 ***
ibd_3Unspecified -1.763538   0.365518  -4.825 1.40e-06 ***
age_yrs          -0.027289   0.001547 -17.639  < 2e-16 ***
genderFemale     -0.210630   0.050438  -4.176 2.97e-05 ***
race_5Black       0.204174   0.103259   1.977  0.04801 *  
race_5Asian      -0.400663   0.150623  -2.660  0.00781 ** 
race_5Native     -0.111078   0.399392  -0.278  0.78092    
race_5Other       0.122775   0.147286   0.834  0.40452    
ethnic_3Hispanic  0.117524   0.177906   0.661  0.50887    
lang_3Other      -0.246638   0.254094  -0.971  0.33172    
max_ch           -0.024830   0.005839  -4.252 2.12e-05 ***
RPL_THEME1       -0.118698   0.142579  -0.833  0.40512    
RPL_THEME2        0.097153   0.125504   0.774  0.43887    
RPL_THEME3        0.129201   0.091917   1.406  0.15984    
RPL_THEME4       -0.260460   0.105051  -2.479  0.01316 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 10855.4  on 7835  degrees of freedom
Residual deviance:  9323.7  on 7820  degrees of freedom
AIC: 9355.7

Number of Fisher Scoring iterations: 4

broom::glance(biologic_rx2 )

broom::tidy(biologic_rx2 , exponentiate = TRUE)

tbl_regression(biologic_rx2, 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.25	0.23, 0.28	<0.001
Unspecified	0.17	0.08, 0.34	<0.001
Age	0.97	0.97, 0.98	<0.001
Gender
Male	—	—
Female	0.81	0.73, 0.89	<0.001
Race
White	—	—
Black	1.23	1.00, 1.50	0.048
Asian	0.67	0.50, 0.90	0.008
Native	0.89	0.40, 1.95	0.8
Other	1.13	0.85, 1.51	0.4
Ethnicity
NonHispanic	—	—
Hispanic	1.12	0.79, 1.60	0.5
Primary Language
English	—	—
Other	0.78	0.47, 1.28	0.3
Charlson Comorbidity Index	0.98	0.96, 0.99	<0.001
Soceioeconomic Status	0.89	0.67, 1.17	0.4
Household Composition	1.10	0.86, 1.41	0.4
Minority Status and Language	1.14	0.95, 1.36	0.2
Housing and Transportation	0.77	0.63, 0.95	0.013
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(biologic_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
9355.689 | 9467.153 |     0.183 | 0.452 | 1.092 |    0.595 |      -Inf |       1.276e-04 | 0.592

performance::check_model(biologic_rx2, panel = TRUE)

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

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

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

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

Infliximab

infliximab_rx2 <- glm(infliximab_2 ~  ibd_3 + age_yrs + gender + race_5 + 
                  ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                  RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = med_rx_2)
summary(infliximab_rx2 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5087  -0.8075  -0.5685   1.0313   2.6720  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       1.0923163  0.1040109  10.502  < 2e-16 ***
ibd_3UC          -0.7201915  0.0558795 -12.888  < 2e-16 ***
ibd_3Unspecified -1.8408739  0.6076414  -3.030  0.00245 ** 
age_yrs          -0.0358916  0.0017815 -20.147  < 2e-16 ***
genderFemale     -0.2503739  0.0543408  -4.607 4.08e-06 ***
race_5Black       0.2781475  0.1048624   2.652  0.00799 ** 
race_5Asian      -0.2334231  0.1663057  -1.404  0.16044    
race_5Native     -0.1779649  0.4757576  -0.374  0.70835    
race_5Other       0.0638614  0.1544066   0.414  0.67917    
ethnic_3Hispanic -0.0776895  0.1891648  -0.411  0.68129    
lang_3Other       0.1414929  0.2835130   0.499  0.61773    
max_ch           -0.0004843  0.0070779  -0.068  0.94545    
RPL_THEME1       -0.1417644  0.1538962  -0.921  0.35696    
RPL_THEME2        0.0590500  0.1354933   0.436  0.66297    
RPL_THEME3        0.0831059  0.0996033   0.834  0.40407    
RPL_THEME4       -0.3695278  0.1135716  -3.254  0.00114 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9141.4  on 7835  degrees of freedom
Residual deviance: 8218.8  on 7820  degrees of freedom
AIC: 8250.8

Number of Fisher Scoring iterations: 4

broom::glance(infliximab_rx2 )

broom::tidy(infliximab_rx2 , exponentiate = TRUE)

tbl_regression(infliximab_rx2, 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.49	0.44, 0.54	<0.001
Unspecified	0.16	0.04, 0.45	0.002
Age	0.96	0.96, 0.97	<0.001
Gender
Male	—	—
Female	0.78	0.70, 0.87	<0.001
Race
White	—	—
Black	1.32	1.07, 1.62	0.008
Asian	0.79	0.57, 1.09	0.2
Native	0.84	0.30, 2.01	0.7
Other	1.07	0.78, 1.44	0.7
Ethnicity
NonHispanic	—	—
Hispanic	0.93	0.63, 1.33	0.7
Primary Language
English	—	—
Other	1.15	0.65, 1.98	0.6
Charlson Comorbidity Index	1.00	0.99, 1.01	>0.9
Soceioeconomic Status	0.87	0.64, 1.17	0.4
Household Composition	1.06	0.81, 1.38	0.7
Minority Status and Language	1.09	0.89, 1.32	0.4
Housing and Transportation	0.69	0.55, 0.86	0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(infliximab_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8250.794 | 8362.258 |     0.119 | 0.416 | 1.025 |    0.524 |      -Inf |       2.110e-04 | 0.653

performance::check_model(infliximab_rx2, panel = TRUE)

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

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

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

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

Adalimumab

ada_rx2 <- glm(adalimumab_2 ~  ibd_3 + age_yrs + gender + race_5 + 
                  ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                  RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = med_rx_2)
summary(ada_rx2 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8941  -0.7388  -0.4964  -0.3946   2.4648  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.1651375  0.1161478 -10.032  < 2e-16 ***
ibd_3UC          -1.0000948  0.0666511 -15.005  < 2e-16 ***
ibd_3Unspecified -0.3861828  0.4198106  -0.920   0.3576    
age_yrs          -0.0002222  0.0018700  -0.119   0.9054    
genderFemale     -0.0281759  0.0616888  -0.457   0.6479    
race_5Black       0.0701357  0.1180682   0.594   0.5525    
race_5Asian      -0.3201678  0.2095517  -1.528   0.1265    
race_5Native     -0.1384152  0.5484260  -0.252   0.8007    
race_5Other      -0.0994968  0.1855083  -0.536   0.5917    
ethnic_3Hispanic  0.1581851  0.2152103   0.735   0.4623    
lang_3Other      -0.3258789  0.3517139  -0.927   0.3542    
max_ch           -0.0523842  0.0083487  -6.275 3.51e-10 ***
RPL_THEME1        0.0468808  0.1733148   0.270   0.7868    
RPL_THEME2        0.1154303  0.1528099   0.755   0.4500    
RPL_THEME3        0.2259481  0.1126143   2.006   0.0448 *  
RPL_THEME4       -0.0016623  0.1282836  -0.013   0.9897    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7194.8  on 7835  degrees of freedom
Residual deviance: 6853.0  on 7820  degrees of freedom
AIC: 6885

Number of Fisher Scoring iterations: 5

broom::glance(ada_rx2 )

broom::tidy(ada_rx2 , exponentiate = TRUE)

tbl_regression(ada_rx2, 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.37	0.32, 0.42	<0.001
Unspecified	0.68	0.27, 1.46	0.4
Age	1.00	1.00, 1.00	>0.9
Gender
Male	—	—
Female	0.97	0.86, 1.10	0.6
Race
White	—	—
Black	1.07	0.85, 1.35	0.6
Asian	0.73	0.47, 1.08	0.13
Native	0.87	0.25, 2.30	0.8
Other	0.91	0.62, 1.29	0.6
Ethnicity
NonHispanic	—	—
Hispanic	1.17	0.76, 1.76	0.5
Primary Language
English	—	—
Other	0.72	0.34, 1.38	0.4
Charlson Comorbidity Index	0.95	0.93, 0.96	<0.001
Soceioeconomic Status	1.05	0.75, 1.47	0.8
Household Composition	1.12	0.83, 1.51	0.5
Minority Status and Language	1.25	1.01, 1.56	0.045
Housing and Transportation	1.00	0.78, 1.28	>0.9
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(ada_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
6884.961 | 6996.424 |     0.042 | 0.369 | 0.936 |    0.437 |  -259.496 |       1.361e-04 | 0.727

performance::check_model(ada_rx2, panel = TRUE)

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

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

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

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

Vedolizumab

vedo_rx2 <- glm(vedolizumab_2 ~  ibd_3 + age_yrs + gender + race_5 + 
                  ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                  RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = med_rx_2)
summary(vedo_rx2 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8163  -0.4700  -0.4291  -0.3872   2.5053  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.893e+00  1.493e-01 -12.682  < 2e-16 ***
ibd_3UC          -2.286e-02  7.864e-02  -0.291 0.771296    
ibd_3Unspecified -1.326e+01  2.283e+02  -0.058 0.953688    
age_yrs          -9.245e-03  2.380e-03  -3.885 0.000102 ***
genderFemale      5.942e-02  7.793e-02   0.763 0.445757    
race_5Black      -3.445e-04  1.594e-01  -0.002 0.998275    
race_5Asian      -1.448e-01  2.406e-01  -0.602 0.547208    
race_5Native      1.128e-01  6.142e-01   0.184 0.854308    
race_5Other       3.825e-01  1.965e-01   1.947 0.051567 .  
ethnic_3Hispanic  2.124e-01  2.438e-01   0.871 0.383514    
lang_3Other       1.017e-01  3.697e-01   0.275 0.783272    
max_ch            5.122e-02  8.152e-03   6.284  3.3e-10 ***
RPL_THEME1        1.343e-01  2.200e-01   0.610 0.541610    
RPL_THEME2       -4.819e-01  1.962e-01  -2.456 0.014055 *  
RPL_THEME3        1.624e-01  1.420e-01   1.143 0.252904    
RPL_THEME4       -2.608e-01  1.629e-01  -1.601 0.109399    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4931.9  on 7835  degrees of freedom
Residual deviance: 4864.4  on 7820  degrees of freedom
AIC: 4896.4

Number of Fisher Scoring iterations: 14

broom::glance(vedo_rx2 )

broom::tidy(vedo_rx2 , exponentiate = TRUE)

tbl_regression(vedo_rx2, 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.98	0.84, 1.14	0.8
Unspecified	0.00	0.00, 0.05	>0.9
Age	0.99	0.99, 1.00	<0.001
Gender
Male	—	—
Female	1.06	0.91, 1.24	0.4
Race
White	—	—
Black	1.00	0.72, 1.35	>0.9
Asian	0.87	0.53, 1.35	0.5
Native	1.12	0.27, 3.21	0.9
Other	1.47	0.98, 2.13	0.052
Ethnicity
NonHispanic	—	—
Hispanic	1.24	0.75, 1.96	0.4
Primary Language
English	—	—
Other	1.11	0.50, 2.17	0.8
Charlson Comorbidity Index	1.05	1.04, 1.07	<0.001
Soceioeconomic Status	1.14	0.74, 1.76	0.5
Household Composition	0.62	0.42, 0.91	0.014
Minority Status and Language	1.18	0.89, 1.55	0.3
Housing and Transportation	0.77	0.56, 1.06	0.11
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(vedo_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
4896.418 | 5007.882 |     0.009 | 0.292 | 0.789 |    0.310 |   -74.766 |       8.071e-04 | 0.829

performance::check_model(vedo_rx2, panel = TRUE)

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

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

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

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

Ustekinumab

uste_rx2 <- glm(ustekinumab_2 ~  ibd_3 + age_yrs + gender + race_5 + 
                  ethnic_3 + lang_3 + max_ch + RPL_THEME1 + 
                  RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = med_rx_2)
summary(uste_rx2 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8146  -0.6480  -0.3222  -0.2805   2.7189  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.303695   0.133113  -9.794  < 2e-16 ***
ibd_3UC          -1.638761   0.089837 -18.241  < 2e-16 ***
ibd_3Unspecified -2.254963   1.013906  -2.224  0.02615 *  
age_yrs          -0.004644   0.002196  -2.115  0.03444 *  
genderFemale      0.193955   0.072218   2.686  0.00724 ** 
race_5Black       0.075631   0.134427   0.563  0.57369    
race_5Asian      -0.430816   0.274326  -1.570  0.11631    
race_5Native     -0.493730   0.751551  -0.657  0.51121    
race_5Other      -0.235413   0.239346  -0.984  0.32533    
ethnic_3Hispanic -0.297271   0.306236  -0.971  0.33169    
lang_3Other      -0.047504   0.418311  -0.114  0.90958    
max_ch           -0.007527   0.008816  -0.854  0.39322    
RPL_THEME1        0.186582   0.202523   0.921  0.35690    
RPL_THEME2        0.070568   0.179227   0.394  0.69378    
RPL_THEME3        0.030294   0.130832   0.232  0.81689    
RPL_THEME4       -0.200830   0.149735  -1.341  0.17984    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5812.6  on 7835  degrees of freedom
Residual deviance: 5321.4  on 7820  degrees of freedom
AIC: 5353.4

Number of Fisher Scoring iterations: 6

broom::glance(uste_rx2 )

broom::tidy(uste_rx2 , exponentiate = TRUE)

tbl_regression(uste_rx2, 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.19	0.16, 0.23	<0.001
Unspecified	0.10	0.01, 0.48	0.026
Age	1.00	0.99, 1.00	0.034
Gender
Male	—	—
Female	1.21	1.05, 1.40	0.007
Race
White	—	—
Black	1.08	0.82, 1.40	0.6
Asian	0.65	0.37, 1.08	0.12
Native	0.61	0.10, 2.14	0.5
Other	0.79	0.48, 1.24	0.3
Ethnicity
NonHispanic	—	—
Hispanic	0.74	0.39, 1.31	0.3
Primary Language
English	—	—
Other	0.95	0.39, 2.03	>0.9
Charlson Comorbidity Index	0.99	0.98, 1.01	0.4
Soceioeconomic Status	1.21	0.81, 1.79	0.4
Household Composition	1.07	0.76, 1.53	0.7
Minority Status and Language	1.03	0.80, 1.33	0.8
Housing and Transportation	0.82	0.61, 1.10	0.2
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(uste_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
5353.404 | 5464.867 |     0.059 | 0.317 | 0.825 |    0.340 |  -127.706 |       3.822e-04 | 0.798

performance::check_model(uste_rx2, panel = TRUE)

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

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

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

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

Small molecules (limited to UC patients)

small_rx2 <- glm(small_2 ~  age_yrs + gender + race_5 + ethnic_3 + lang_3 
                     + max_ch + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = UC_meds)
summary(small_rx2 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5209  -0.3113  -0.2681  -0.2286   3.0166  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -2.471234   0.323463  -7.640 2.17e-14 ***
age_yrs           -0.005255   0.005021  -1.047  0.29526    
genderFemale      -0.521588   0.171900  -3.034  0.00241 ** 
race_5Black       -0.195435   0.438533  -0.446  0.65585    
race_5Asian       -0.469641   0.529691  -0.887  0.37528    
race_5Native     -12.268729 337.483036  -0.036  0.97100    
race_5Other        0.408990   0.395401   1.034  0.30096    
ethnic_3Hispanic  -0.994946   0.748889  -1.329  0.18399    
lang_3Other       -0.336405   1.041295  -0.323  0.74665    
max_ch            -0.033320   0.021051  -1.583  0.11346    
RPL_THEME1        -0.933555   0.507655  -1.839  0.06592 .  
RPL_THEME2         0.684031   0.433985   1.576  0.11499    
RPL_THEME3        -0.104698   0.311944  -0.336  0.73715    
RPL_THEME4        -0.056519   0.360053  -0.157  0.87526    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1243.2  on 3783  degrees of freedom
Residual deviance: 1215.8  on 3770  degrees of freedom
AIC: 1243.8

Number of Fisher Scoring iterations: 14

broom::glance(small_rx2 )

broom::tidy(small_rx2 , exponentiate = TRUE)

tbl_regression(small_rx2, 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	1.0	0.98, 1.00	0.3
Gender
Male	—	—
Female	0.59	0.42, 0.83	0.002
Race
White	—	—
Black	0.82	0.31, 1.79	0.7
Asian	0.63	0.19, 1.57	0.4
Native	0.00		>0.9
Other	1.51	0.64, 3.08	0.3
Ethnicity
NonHispanic	—	—
Hispanic	0.37	0.06, 1.28	0.2
Primary Language
English	—	—
Other	0.71	0.04, 3.60	0.7
Charlson Comorbidity Index	0.97	0.93, 1.01	0.11
Soceioeconomic Status	0.39	0.14, 1.05	0.066
Household Composition	1.98	0.85, 4.65	0.11
Minority Status and Language	0.90	0.49, 1.66	0.7
Housing and Transportation	0.95	0.46, 1.91	0.9
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(small_rx2)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1243.834 | 1331.174 |     0.008 | 0.192 | 0.568 |    0.161 |    -5.773 |           0.005 | 0.926

performance::check_model(small_rx2, panel = TRUE)

# Margins 
cplot(small_rx2, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Small Molecule Rx Given Theme 1")

cplot(small_rx2, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Small Molecule Rx Given THEME2")

cplot(small_rx2, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Small Molecule Rx Given THEME3")

cplot(small_rx2, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Small Molecule Rx Given THEME4")

Steroids

steroids_rx2 <- glm(steroids_2 ~  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(steroids_rx2 )

Call:
glm(formula = steroids_2 ~ 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  
-2.1862  -1.2026   0.7198   1.0884   1.5242  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.174584   0.091044   1.918  0.05516 .  
ibd_3UC           0.001206   0.047968   0.025  0.97994    
ibd_3Unspecified -0.331217   0.327634  -1.011  0.31205    
age_yrs          -0.007808   0.001430  -5.462 4.70e-08 ***
genderFemale      0.216942   0.047362   4.581 4.64e-06 ***
race_5Black       0.196869   0.098189   2.005  0.04496 *  
race_5Asian      -0.233178   0.139356  -1.673  0.09428 .  
race_5Native     -0.019671   0.384849  -0.051  0.95924    
race_5Other      -0.087439   0.136379  -0.641  0.52143    
ethnic_3Hispanic  0.114002   0.165343   0.689  0.49052    
lang_3Other      -0.279139   0.234392  -1.191  0.23369    
max_ch            0.109789   0.006375  17.223  < 2e-16 ***
RPL_THEME1        0.358779   0.134015   2.677  0.00742 ** 
RPL_THEME2        0.117920   0.117862   1.000  0.31707    
RPL_THEME3       -0.197243   0.086499  -2.280  0.02259 *  
RPL_THEME4       -0.038961   0.098837  -0.394  0.69344    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 10672  on 7835  degrees of freedom
Residual deviance: 10231  on 7820  degrees of freedom
AIC: 10263

Number of Fisher Scoring iterations: 4

broom::glance(steroids_rx2 )

broom::tidy(steroids_rx2 , exponentiate = TRUE)

tbl_regression(steroids_rx2, 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.91, 1.10	>0.9
Unspecified	0.72	0.38, 1.37	0.3
Age	0.99	0.99, 1.00	<0.001
Gender
Male	—	—
Female	1.24	1.13, 1.36	<0.001
Race
White	—	—
Black	1.22	1.01, 1.48	0.045
Asian	0.79	0.60, 1.04	0.094
Native	0.98	0.46, 2.13	>0.9
Other	0.92	0.70, 1.20	0.5
Ethnicity
NonHispanic	—	—
Hispanic	1.12	0.81, 1.55	0.5
Primary Language
English	—	—
Other	0.76	0.48, 1.20	0.2
Charlson Comorbidity Index	1.12	1.10, 1.13	<0.001
Soceioeconomic Status	1.43	1.10, 1.86	0.007
Household Composition	1.13	0.89, 1.42	0.3
Minority Status and Language	0.82	0.69, 0.97	0.023
Housing and Transportation	0.96	0.79, 1.17	0.7
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(steroids_rx2)

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log |   PCP
--------------------------------------------------------------------------------
10262.916 | 10374.380 |     0.053 | 0.481 | 1.144 |    0.653 |      -Inf | 0.538

performance::check_model(steroids_rx2, panel = TRUE)

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

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

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

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