Flu Models

Load Packages

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

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

Loading required package: mvtnorm
mediation: Causal Mediation Analysis
Version: 4.5.0

Import Data

load("~/Desktop/R-Code/SDOH_ALL/Vax_dataset.rda")

Codebook

vax_dataset %>% 
dplyr::select(ibd_3, age_yrs, gender, race_5, ethnic_3, lang_3, relig_affil, mstat_5, max_ch, STATE, r_pct, IC, pop_dens, RPL_THEMES, RPL_4, RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4, act_tob, r_pct, flu_2015, flu_2016, flu_2017, flu_2018, flu_2019, flu_2020, flu_2021, flu_2022, total_flu, prevnar, pvax, any_pneum, any_pneum, total_cov_vax, total_shingrix) -> vax_clean1
print(dfSummary(vax_clean1), method = 'render')

Data Frame Summary

vax_clean1

Dimensions: 15208 x 34
Duplicates: 0

No	Variable	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	ibd_3 [factor]	1. CD 2. UC 3. Unspecified	7712 ( 50.7% ) 7345 ( 48.3% ) 151 ( 1.0% )		15208 (100.0%)	0 (0.0%)
2	age_yrs [numeric]	Mean (sd) : 49.3 (19.5) min ≤ med ≤ max: 1 ≤ 49.3 ≤ 108.3 IQR (CV) : 31.6 (0.4)	11541 distinct values		15208 (100.0%)	0 (0.0%)
3	gender [factor]	1. Male 2. Female	6962 ( 45.8% ) 8246 ( 54.2% )		15208 (100.0%)	0 (0.0%)
4	race_5 [factor]	1. White 2. Black 3. Asian 4. Native 5. Other	13202 ( 86.8% ) 946 ( 6.2% ) 373 ( 2.5% ) 56 ( 0.4% ) 631 ( 4.1% )		15208 (100.0%)	0 (0.0%)
5	ethnic_3 [factor]	1. NonHispanic 2. UNKNOWN 3. Hispanic 4. CHOOSE NOT TO DISCLOSE	14371 ( 98.0% ) 0 ( 0.0% ) 299 ( 2.0% ) 0 ( 0.0% )		14670 (96.5%)	538 (3.5%)
6	lang_3 [factor]	1. English 2. Other	15046 ( 98.9% ) 162 ( 1.1% )		15208 (100.0%)	0 (0.0%)
7	relig_affil [factor]	1. Yes 2. No 3. UNKNOWN 4. PATIENT REFUSED	8193 ( 57.4% ) 6075 ( 42.6% ) 0 ( 0.0% ) 0 ( 0.0% )		14268 (93.8%)	940 (6.2%)
8	mstat_5 [factor]	1. Married 2. Unknown 3. Unmarried 4. DivorcedSeparated 5. Widow	6178 ( 40.6% ) 3226 ( 21.2% ) 4990 ( 32.8% ) 507 ( 3.3% ) 307 ( 2.0% )		15208 (100.0%)	0 (0.0%)
9	max_ch [numeric]	Mean (sd) : 3.3 (4.9) min ≤ med ≤ max: 0 ≤ 1 ≤ 27 IQR (CV) : 5 (1.5)	28 distinct values		14773 (97.1%)	435 (2.9%)
10	STATE [character]	1. FLORIDA 2. ILLINOIS 3. INDIANA 4. MARYLAND 5. MICHIGAN 6. NEW YORK 7. NORTH CAROLINA 8. OHIO 9. PENNSYLVANIA 10. TEXAS	24 ( 0.2% ) 40 ( 0.3% ) 47 ( 0.3% ) 16 ( 0.1% ) 14069 ( 96.5% ) 38 ( 0.3% ) 10 ( 0.1% ) 314 ( 2.2% ) 14 ( 0.1% ) 7 ( 0.0% )		14579 (95.9%)	629 (4.1%)
11	r_pct [numeric]	Mean (sd) : 44.7 (18) min ≤ med ≤ max: 1.2 ≤ 47.9 ≤ 84.8 IQR (CV) : 24.8 (0.4)	2624 distinct values		12943 (85.1%)	2265 (14.9%)
12	IC [numeric]	Min : 0 Mean : 0.6 Max : 1	0 : 2946 ( 37.6% ) 1 : 4890 ( 62.4% )		7836 (51.5%)	7372 (48.5%)
13	pop_dens [numeric]	Mean (sd) : 2290.3 (4015.6) min ≤ med ≤ max: 2.5 ≤ 1523.4 ≤ 166336.1 IQR (CV) : 2985 (1.8)	2641 distinct values		14579 (95.9%)	629 (4.1%)
14	RPL_THEMES [numeric]	Mean (sd) : 0.4 (0.3) min ≤ med ≤ max: 0 ≤ 0.3 ≤ 1 IQR (CV) : 0.4 (0.7)	2538 distinct values		14920 (98.1%)	288 (1.9%)
15	RPL_4 [factor]	1. First 2. Second 3. Third 4. Fourth	5833 ( 39.1% ) 4492 ( 30.1% ) 3052 ( 20.5% ) 1543 ( 10.3% )		14920 (98.1%)	288 (1.9%)
16	RPL_THEME1 [numeric]	Mean (sd) : 0.4 (0.3) min ≤ med ≤ max: 0 ≤ 0.3 ≤ 1 IQR (CV) : 0.4 (0.7)	2476 distinct values		14870 (97.8%)	338 (2.2%)
17	RPL_THEME2 [numeric]	Mean (sd) : 0.4 (0.3) min ≤ med ≤ max: 0 ≤ 0.3 ≤ 1 IQR (CV) : 0.4 (0.7)	2450 distinct values		14921 (98.1%)	287 (1.9%)
18	RPL_THEME3 [numeric]	Mean (sd) : 0.5 (0.3) min ≤ med ≤ max: 0 ≤ 0.5 ≤ 1 IQR (CV) : 0.5 (0.6)	1831 distinct values		14929 (98.2%)	279 (1.8%)
19	RPL_THEME4 [numeric]	Mean (sd) : 0.4 (0.3) min ≤ med ≤ max: 0 ≤ 0.4 ≤ 1 IQR (CV) : 0.5 (0.7)	2454 distinct values		14898 (98.0%)	310 (2.0%)
20	act_tob [factor]	1. No 2. Yes 3. NOT ASKED	12593 ( 87.0% ) 1878 ( 13.0% ) 0 ( 0.0% )		14471 (95.2%)	737 (4.8%)
21	flu_2015 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 10815 ( 71.1% ) 1 : 4393 ( 28.9% )		15208 (100.0%)	0 (0.0%)
22	flu_2016 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 10616 ( 69.8% ) 1 : 4592 ( 30.2% )		15208 (100.0%)	0 (0.0%)
23	flu_2017 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 10756 ( 70.7% ) 1 : 4452 ( 29.3% )		15208 (100.0%)	0 (0.0%)
24	flu_2018 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 10645 ( 70.0% ) 1 : 4563 ( 30.0% )		15208 (100.0%)	0 (0.0%)
25	flu_2019 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 10679 ( 70.2% ) 1 : 4529 ( 29.8% )		15208 (100.0%)	0 (0.0%)
26	flu_2020 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 10596 ( 69.7% ) 1 : 4612 ( 30.3% )		15208 (100.0%)	0 (0.0%)
27	flu_2021 [numeric]	Min : 0 Mean : 0.2 Max : 1	0 : 12014 ( 79.0% ) 1 : 3194 ( 21.0% )		15208 (100.0%)	0 (0.0%)
28	flu_2022 [numeric]	Min : 0 Mean : 0 Max : 1	0 : 15044 ( 98.9% ) 1 : 164 ( 1.1% )		15208 (100.0%)	0 (0.0%)
29	total_flu [numeric]	Mean (sd) : 2 (2.3) min ≤ med ≤ max: 0 ≤ 1 ≤ 8 IQR (CV) : 4 (1.2)	0 : 6500 ( 42.7% ) 1 : 2012 ( 13.2% ) 2 : 1373 ( 9.0% ) 3 : 1244 ( 8.2% ) 4 : 1160 ( 7.6% ) 5 : 1058 ( 7.0% ) 6 : 955 ( 6.3% ) 7 : 899 ( 5.9% ) 8 : 7 ( 0.0% )		15208 (100.0%)	0 (0.0%)
30	prevnar [numeric]	Mean (sd) : 0.3 (0.5) min ≤ med ≤ max: 0 ≤ 0 ≤ 8 IQR (CV) : 1 (1.7)	0 : 10700 ( 70.4% ) 1 : 4235 ( 27.8% ) 2 : 210 ( 1.4% ) 3 : 32 ( 0.2% ) 4 : 24 ( 0.2% ) 5 : 4 ( 0.0% ) 6 : 1 ( 0.0% ) 7 : 1 ( 0.0% ) 8 : 1 ( 0.0% )		15208 (100.0%)	0 (0.0%)
31	pvax [numeric]	Mean (sd) : 0.3 (0.5) min ≤ med ≤ max: 0 ≤ 0 ≤ 6 IQR (CV) : 1 (1.8)	0 : 11222 ( 73.8% ) 1 : 3692 ( 24.3% ) 2 : 280 ( 1.8% ) 3 : 11 ( 0.1% ) 5 : 1 ( 0.0% ) 6 : 2 ( 0.0% )		15208 (100.0%)	0 (0.0%)
32	any_pneum [numeric]	Mean (sd) : 0.6 (0.9) min ≤ med ≤ max: 0 ≤ 0 ≤ 9 IQR (CV) : 1 (1.5)	0 : 9474 ( 62.3% ) 1 : 2816 ( 18.5% ) 2 : 2517 ( 16.6% ) 3 : 309 ( 2.0% ) 4 : 65 ( 0.4% ) 5 : 16 ( 0.1% ) 6 : 7 ( 0.0% ) 7 : 2 ( 0.0% ) 8 : 1 ( 0.0% ) 9 : 1 ( 0.0% )		15208 (100.0%)	0 (0.0%)
33	total_cov_vax [numeric]	Mean (sd) : 1.4 (1.6) min ≤ med ≤ max: 0 ≤ 0 ≤ 7 IQR (CV) : 3 (1.1)	0 : 7737 ( 50.9% ) 1 : 433 ( 2.8% ) 2 : 2203 ( 14.5% ) 3 : 3139 ( 20.6% ) 4 : 1487 ( 9.8% ) 5 : 161 ( 1.1% ) 6 : 39 ( 0.3% ) 7 : 9 ( 0.1% )		15208 (100.0%)	0 (0.0%)
34	total_shingrix [numeric]	Mean (sd) : 0.2 (0.6) min ≤ med ≤ max: 0 ≤ 0 ≤ 5 IQR (CV) : 0 (3.1)	0 : 13624 ( 89.6% ) 1 : 366 ( 2.4% ) 2 : 1115 ( 7.3% ) 3 : 91 ( 0.6% ) 4 : 10 ( 0.1% ) 5 : 2 ( 0.0% )		15208 (100.0%)	0 (0.0%)

Generated by summarytools 1.0.1 (R version 4.2.1)
2022-11-26

NA

save(vax_clean1, file = "vax_clean1.rda")

Patient Characteristics

Baseline Characteristics

vax_clean1 %>% 
  dplyr::select(ibd_3, age_yrs, gender, race_5, ethnic_3, lang_3, relig_affil, mstat_5, act_tob, max_ch, IC, pop_dens,r_pct, 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", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", pop_dens ~ "Population Density", RPL_4 ~ "SVI Quartiles", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC = "Immunocompromised"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)")

Characteristic	N = 15,2081
IBD Diagnosis
CD	7,712 (51%)
UC	7,345 (48%)
Unspecified	151 (1.0%)
Age	49 (19)
Gender
Male	6,962 (46%)
Female	8,246 (54%)
Race
White	13,202 (87%)
Black	946 (6.2%)
Asian	373 (2.5%)
Native	56 (0.4%)
Other	631 (4.1%)
Ethnicity
NonHispanic	14,371 (98%)
UNKNOWN	0 (0%)
Hispanic	299 (2.0%)
CHOOSE NOT TO DISCLOSE	0 (0%)
(Missing)	538
Primary Language
English	15,046 (99%)
Other	162 (1.1%)
Any Religious Affiliation
Yes	8,193 (57%)
No	6,075 (43%)
UNKNOWN	0 (0%)
PATIENT REFUSED	0 (0%)
(Missing)	940
Marital Status
Married	6,178 (41%)
Unknown	3,226 (21%)
Unmarried	4,990 (33%)
DivorcedSeparated	507 (3.3%)
Widow	307 (2.0%)
Active Tobacco Use
No	12,593 (87%)
Yes	1,878 (13%)
NOT ASKED	0 (0%)
(Missing)	737
Charlson Comorbidity Index	3.3 (4.9)
(Missing)	435
Immunocompromised	4,890 (62%)
(Missing)	7,372
Population Density	2,290 (4,016)
(Missing)	629
Percent Republican	45 (18)
(Missing)	2,265
Total SVI	0.37 (0.26)
(Missing)	288
SVI Quartiles
First	5,833 (39%)
Second	4,492 (30%)
Third	3,052 (20%)
Fourth	1,543 (10%)
(Missing)	288
Soceioeconomic Status	0.35 (0.26)
(Missing)	338
Household Composition	0.40 (0.27)
(Missing)	287
Minority Status and Language	0.48 (0.29)
(Missing)	279
Housing and Transportation	0.44 (0.29)
(Missing)	310
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", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", pop_dens ~ "Population Density", RPL_4 ~ "SVI Quartiles", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()

288 observations missing `RPL_4` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `RPL_4` column before passing to `tbl_summary()`.
There was an error in 'add_p()/add_difference()' for variable 'relig_affil', p-value omitted:
Error in stats::fisher.test(structure(c(2L, NA, 2L, NA, 1L, 1L, 2L, NA, : FEXACT error 6.  LDKEY=534 is too small for this problem,
  (ii := key2[itp=1059] = 45352598, ldstp=16020)
Try increasing the size of the workspace and possibly 'mult'
There was an error in 'add_p()/add_difference()' for variable 'act_tob', p-value omitted:
Error in stats::fisher.test(structure(c(1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, : FEXACT error 6.  LDKEY=533 is too small for this problem,
  (ii := key2[itp=950] = 21479782, ldstp=15990)
Try increasing the size of the workspace and possibly 'mult'

Characteristic	First, N = 5,8331	Second, N = 4,4921	Third, N = 3,0521	Fourth, N = 1,5431	p-value2
IBD Diagnosis					<0.001
CD	2,773 (48%)	2,208 (49%)	1,658 (54%)	914 (59%)
UC	3,009 (52%)	2,244 (50%)	1,356 (44%)	611 (40%)
Unspecified	51 (0.9%)	40 (0.9%)	38 (1.2%)	18 (1.2%)
Age	49 (20)	50 (19)	50 (19)	48 (19)	0.008
Gender					<0.001
Male	2,756 (47%)	2,050 (46%)	1,303 (43%)	691 (45%)
Female	3,077 (53%)	2,442 (54%)	1,749 (57%)	852 (55%)
Race					<0.001
White	5,293 (91%)	3,959 (88%)	2,619 (86%)	1,091 (71%)
Black	139 (2.4%)	204 (4.5%)	239 (7.8%)	345 (22%)
Asian	153 (2.6%)	140 (3.1%)	65 (2.1%)	7 (0.5%)
Native	14 (0.2%)	24 (0.5%)	10 (0.3%)	7 (0.5%)
Other	234 (4.0%)	165 (3.7%)	119 (3.9%)	93 (6.0%)
Ethnicity					0.010
NonHispanic	5,517 (98%)	4,262 (98%)	2,888 (98%)	1,440 (97%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Hispanic	100 (1.8%)	86 (2.0%)	59 (2.0%)	48 (3.2%)
CHOOSE NOT TO DISCLOSE	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	216	144	105	55
Primary Language					<0.001
English	5,799 (99%)	4,441 (99%)	3,009 (99%)	1,511 (98%)
Other	34 (0.6%)	51 (1.1%)	43 (1.4%)	32 (2.1%)
Any Religious Affiliation
Yes	3,232 (59%)	2,397 (57%)	1,651 (58%)	779 (55%)
No	2,284 (41%)	1,810 (43%)	1,204 (42%)	642 (45%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)
PATIENT REFUSED	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	317	285	197	122
Marital Status					<0.001
Married	2,695 (46%)	1,848 (41%)	1,109 (36%)	444 (29%)
Unknown	1,134 (19%)	943 (21%)	702 (23%)	388 (25%)
Unmarried	1,761 (30%)	1,432 (32%)	1,039 (34%)	619 (40%)
DivorcedSeparated	133 (2.3%)	168 (3.7%)	128 (4.2%)	72 (4.7%)
Widow	110 (1.9%)	101 (2.2%)	74 (2.4%)	20 (1.3%)
Active Tobacco Use
No	5,133 (91%)	3,710 (87%)	2,389 (83%)	1,120 (78%)
Yes	478 (8.5%)	566 (13%)	501 (17%)	308 (22%)
NOT ASKED	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	222	216	162	115
Charlson Comorbidity Index	3.2 (4.9)	3.4 (5.1)	3.4 (4.9)	3.3 (4.8)	0.020
(Missing)	135	128	93	62
Immunocompromised	1,955 (61%)	1,485 (62%)	945 (64%)	505 (66%)	0.047
(Missing)	2,634	2,100	1,572	778
Population Density	2,026 (3,289)	2,277 (5,130)	2,431 (3,595)	3,081 (3,469)	<0.001
(Missing)	55	97	116	81
Percent Republican	46 (16)	46 (18)	46 (18)	37 (21)	<0.001
(Missing)	614	829	542	263
Total SVI	0.12 (0.07)	0.37 (0.07)	0.61 (0.07)	0.86 (0.07)	<0.001
Soceioeconomic Status	0.13 (0.11)	0.34 (0.16)	0.56 (0.15)	0.78 (0.11)	<0.001
(Missing)	51	0	0	0
Household Composition	0.21 (0.14)	0.38 (0.22)	0.58 (0.24)	0.76 (0.19)	<0.001
Minority Status and Language	0.41 (0.27)	0.47 (0.29)	0.52 (0.28)	0.66 (0.25)	<0.001
Housing and Transportation	0.19 (0.16)	0.49 (0.20)	0.64 (0.21)	0.81 (0.17)	<0.001
(Missing)	22	0	0	0
1 n (%); Mean (SD)
2 Pearson's Chi-squared test; Kruskal-Wallis rank sum test; Fisher's exact test

Bivariate Analysis

Total Flu

tbl_uv_ex1 <-
  tbl_uvregression(
    vax_clean1[c("total_flu", "ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "mstat_5", "relig_affil", "act_tob", "max_ch", "IC", "pop_dens", "r_pct", "RPL_THEMES", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4")],
    method = glm.nb,
    y = total_flu,
    label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised"),
  exponentiate = TRUE)
print(tbl_uv_ex1, method = render)

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

Characteristic	N	IRR1	95% CI1	p-value
IBD Diagnosis	15,208
CD		—	—
UC		1.07	1.03, 1.12	0.002
Unspecified		0.38	0.29, 0.49	<0.001
Age	15,208	1.00	1.00, 1.00	0.039
Gender	15,208
Male		—	—
Female		1.08	1.03, 1.12	0.001
Race	15,208
White		—	—
Black		0.86	0.78, 0.94	0.001
Asian		1.30	1.13, 1.49	<0.001
Native		1.09	0.76, 1.58	0.7
Other		0.83	0.74, 0.93	0.001
Ethnicity	14,670
NonHispanic		—	—
Hispanic		1.10	0.94, 1.28	0.3
Primary Language	15,208
English		—	—
Other		0.72	0.57, 0.90	0.004
Marital Status	15,208
Married		—	—
Unknown		0.76	0.72, 0.81	<0.001
Unmarried		0.90	0.85, 0.95	<0.001
DivorcedSeparated		0.83	0.73, 0.95	0.005
Widow		0.90	0.77, 1.06	0.2
Any Religious Affiliation	14,268
Yes		—	—
No		0.92	0.88, 0.96	<0.001
Active Tobacco Use	14,471
No		—	—
Yes		0.69	0.65, 0.74	<0.001
Charlson Comorbidity Index	14,773	1.02	1.02, 1.03	<0.001
Immunocompromised	7,836	1.23	1.17, 1.28	<0.001
Population Density	14,579	1.00	1.00, 1.00	0.8
Percent Republican	12,943	0.99	0.99, 0.99	<0.001
Total SVI	14,920	0.49	0.45, 0.54	<0.001
Soceioeconomic Status	14,870	0.48	0.44, 0.52	<0.001
Household Composition	14,921	0.49	0.45, 0.53	<0.001
Minority Status and Language	14,929	1.19	1.10, 1.28	<0.001
Housing and Transportation	14,898	0.66	0.61, 0.71	<0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

NULL

Flu 2019

tbl_uv_ex2 <-
  tbl_uvregression(
    vax_clean1[c("flu_2019", "ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "mstat_5", "relig_affil", "act_tob", "max_ch", "IC", "pop_dens", "r_pct", "RPL_THEMES", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4")],
    method = glm,
    y = flu_2019,
    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", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised"),
    exponentiate = TRUE
  )
print(tbl_uv_ex2, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	15,208
CD		—	—
UC		1.04	0.97, 1.12	0.3
Unspecified		0.30	0.18, 0.49	<0.001
Age	15,208	0.99	0.99, 0.99	<0.001
Gender	15,208
Male		—	—
Female		1.09	1.02, 1.17	0.012
Race	15,208
White		—	—
Black		0.81	0.69, 0.94	0.006
Asian		1.64	1.33, 2.03	<0.001
Native		1.20	0.68, 2.06	0.5
Other		0.78	0.65, 0.94	0.008
Ethnicity	14,670
NonHispanic		—	—
Hispanic		1.40	1.10, 1.77	0.005
Primary Language	15,208
English		—	—
Other		0.62	0.42, 0.90	0.015
Marital Status	15,208
Married		—	—
Unknown		0.88	0.81, 0.97	0.010
Unmarried		0.93	0.86, 1.01	0.084
DivorcedSeparated		0.83	0.67, 1.01	0.064
Widow		0.56	0.42, 0.74	<0.001
Any Religious Affiliation	14,268
Yes		—	—
No		0.93	0.87, 1.00	0.054
Active Tobacco Use	14,471
No		—	—
Yes		0.64	0.57, 0.71	<0.001
Charlson Comorbidity Index	14,773	1.01	1.00, 1.02	0.015
Immunocompromised	7,836	1.74	1.58, 1.91	<0.001
Population Density	14,579	1.00	1.00, 1.00	0.8
Percent Republican	12,943	0.99	0.99, 0.99	<0.001
Total SVI	14,920	0.37	0.32, 0.43	<0.001
Soceioeconomic Status	14,870	0.34	0.30, 0.40	<0.001
Household Composition	14,921	0.38	0.33, 0.43	<0.001
Minority Status and Language	14,929	1.31	1.16, 1.48	<0.001
Housing and Transportation	14,898	0.54	0.48, 0.61	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Flu 2021

tbl_uv_ex3 <-
  tbl_uvregression(
    vax_clean1[c("flu_2021", "ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "mstat_5", "relig_affil", "act_tob", "max_ch", "IC", "pop_dens", "r_pct", "RPL_THEMES", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4")],
    method = glm,
    y = flu_2021,
    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", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised"),
    exponentiate = TRUE
  )
print(tbl_uv_ex3, method = render)

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

Characteristic	N	OR1	95% CI1	p-value
IBD Diagnosis	15,208
CD		—	—
UC		1.18	1.09, 1.28	<0.001
Unspecified		0.23	0.10, 0.43	<0.001
Age	15,208	1.01	1.00, 1.01	<0.001
Gender	15,208
Male		—	—
Female		1.16	1.08, 1.26	<0.001
Race	15,208
White		—	—
Black		0.76	0.64, 0.90	0.002
Asian		1.95	1.56, 2.42	<0.001
Native		1.02	0.51, 1.87	>0.9
Other		0.74	0.59, 0.91	0.005
Ethnicity	14,670
NonHispanic		—	—
Hispanic		1.03	0.78, 1.36	0.8
Primary Language	15,208
English		—	—
Other		0.65	0.41, 0.99	0.054
Marital Status	15,208
Married		—	—
Unknown		0.76	0.69, 0.85	<0.001
Unmarried		0.69	0.63, 0.75	<0.001
DivorcedSeparated		0.79	0.63, 0.99	0.043
Widow		1.07	0.82, 1.39	0.6
Any Religious Affiliation	14,268
Yes		—	—
No		0.90	0.83, 0.98	0.012
Active Tobacco Use	14,471
No		—	—
Yes		0.46	0.40, 0.53	<0.001
Charlson Comorbidity Index	14,773	1.05	1.04, 1.06	<0.001
Immunocompromised	7,836	1.13	1.02, 1.25	0.019
Population Density	14,579	1.00	1.00, 1.00	0.3
Percent Republican	12,943	0.98	0.98, 0.98	<0.001
Total SVI	14,920	0.27	0.22, 0.31	<0.001
Soceioeconomic Status	14,870	0.23	0.19, 0.27	<0.001
Household Composition	14,921	0.24	0.20, 0.28	<0.001
Minority Status and Language	14,929	1.59	1.39, 1.83	<0.001
Housing and Transportation	14,898	0.48	0.42, 0.55	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

NULL

Total Flu Multivariable Model

RPL_THEMES Poisson Model

totalflu_SVI <- glm(total_flu ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch 
                      + pop_dens + r_pct + RPL_THEMES,
               family = poisson,
               data = vax_clean1) 
summary(totalflu_SVI)

Call:
glm(formula = total_flu ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    pop_dens + r_pct + RPL_THEMES, family = poisson, data = vax_clean1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-3.650  -1.861  -0.546   1.044   4.238  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.787e+00  3.622e-02  49.333  < 2e-16 ***
ibd_3UC                  -3.237e-02  1.305e-02  -2.480 0.013137 *  
ibd_3Unspecified         -9.004e-01  1.089e-01  -8.268  < 2e-16 ***
age_yrs                  -6.021e-03  4.370e-04 -13.777  < 2e-16 ***
genderFemale              8.448e-02  1.300e-02   6.500 8.03e-11 ***
race_5Black              -1.494e-01  2.980e-02  -5.013 5.35e-07 ***
race_5Asian               1.294e-01  3.785e-02   3.419 0.000628 ***
race_5Native              8.935e-02  1.048e-01   0.853 0.393748    
race_5Other              -1.646e-01  4.125e-02  -3.989 6.63e-05 ***
ethnic_3Hispanic          1.494e-01  4.600e-02   3.247 0.001166 ** 
lang_3Other              -3.328e-01  7.679e-02  -4.333 1.47e-05 ***
mstat_5Unknown           -2.775e-01  1.897e-02 -14.629  < 2e-16 ***
mstat_5Unmarried         -1.470e-01  1.713e-02  -8.584  < 2e-16 ***
mstat_5DivorcedSeparated -1.346e-01  3.819e-02  -3.523 0.000427 ***
mstat_5Widow             -1.393e-01  4.675e-02  -2.979 0.002894 ** 
relig_affilNo            -6.447e-02  1.340e-02  -4.810 1.51e-06 ***
act_tobYes               -2.459e-01  2.247e-02 -10.945  < 2e-16 ***
max_ch                    2.773e-02  1.333e-03  20.804  < 2e-16 ***
pop_dens                 -1.348e-05  2.370e-06  -5.686 1.30e-08 ***
r_pct                    -9.507e-03  4.188e-04 -22.702  < 2e-16 ***
RPL_THEMES               -5.890e-01  2.716e-02 -21.686  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 31304  on 10973  degrees of freedom
Residual deviance: 29047  on 10953  degrees of freedom
  (4234 observations deleted due to missingness)
AIC: 49531

Number of Fisher Scoring iterations: 6

broom::glance(totalflu_SVI)

broom::tidy(totalflu_SVI, exponentiate = TRUE)

model_performance(totalflu_SVI)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
49530.535 | 49683.904 |           0.197 | 2.258 | 1.628 |    -2.255 |           0.008

tbl_regression(totalflu_SVI, label = list(ibd_3 ~ "IBD Diagnosis", age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Primary Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Diagnosis
CD	—	—
UC	0.97	0.94, 0.99	0.013
Unspecified	0.41	0.33, 0.50	<0.001
Age	0.99	0.99, 1.0	<0.001
Gender
Male	—	—
Female	1.09	1.06, 1.12	<0.001
Race
White	—	—
Black	0.86	0.81, 0.91	<0.001
Asian	1.14	1.06, 1.22	<0.001
Native	1.09	0.88, 1.33	0.4
Other	0.85	0.78, 0.92	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.16	1.06, 1.27	0.001
Primary Language
English	—	—
Other	0.72	0.61, 0.83	<0.001
Marital Status
Married	—	—
Unknown	0.76	0.73, 0.79	<0.001
Unmarried	0.86	0.83, 0.89	<0.001
DivorcedSeparated	0.87	0.81, 0.94	<0.001
Widow	0.87	0.79, 0.95	0.003
Any Religious Affiliation
Yes	—	—
No	0.94	0.91, 0.96	<0.001
Active Tobacco Use
No	—	—
Yes	0.78	0.75, 0.82	<0.001
Charlson Comorbidity Index	1.03	1.03, 1.03	<0.001
pop_dens	1.00	1.00, 1.00	<0.001
r_pct	0.99	0.99, 0.99	<0.001
Total SVI	0.55	0.53, 0.59	<0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# Residual plots 
totalflu_res <- resid(totalflu_SVI)
plot(fitted(totalflu_SVI), totalflu_res, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Poisson')
abline(0,0)

# Poisson model assumes that variance is equal to the mean. If variance is significantly different than the mean, NB regression is better. Here was see that the standardized residuals fall outside the -2 to 2 range, indicating the NB regression is more appropriate. 

library(AER)
dispersiontest(totalflu_SVI, trafo=1)

    Overdispersion test

data:  totalflu_SVI
z = 59.183, p-value < 2.2e-16
alternative hypothesis: true alpha is greater than 0
sample estimates:
   alpha 
1.354895 

RPL_THEMES Negative Binomial Model

library(MASS)
totalflu_SVI.nb <- glm.nb(total_flu ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEMES,
               data = vax_clean1) 
summary(totalflu_SVI.nb)

Call:
glm.nb(formula = total_flu ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_THEMES, data = vax_clean1, init.theta = 2.161938564, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6392  -1.2142  -0.1269   0.5989   2.6103  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.399e+00  7.007e-02  19.966  < 2e-16 ***
ibd_3UC                   1.402e-02  2.455e-02   0.571 0.567929    
ibd_3Unspecified         -7.809e-01  1.980e-01  -3.944 8.02e-05 ***
age_yrs                   2.534e-04  8.060e-04   0.314 0.753184    
genderFemale              1.103e-01  2.314e-02   4.767 1.87e-06 ***
race_5Black              -2.280e-01  5.004e-02  -4.557 5.20e-06 ***
race_5Asian               1.487e-01  6.731e-02   2.209 0.027184 *  
race_5Native              1.217e-01  1.819e-01   0.669 0.503579    
race_5Other              -2.122e-01  6.872e-02  -3.088 0.002014 ** 
ethnic_3Hispanic          2.018e-01  7.985e-02   2.527 0.011496 *  
lang_3Other              -4.398e-01  1.281e-01  -3.433 0.000597 ***
mstat_5Unknown           -2.758e-01  3.334e-02  -8.271  < 2e-16 ***
mstat_5Unmarried         -1.444e-01  3.086e-02  -4.681 2.85e-06 ***
mstat_5DivorcedSeparated -1.106e-01  7.141e-02  -1.549 0.121416    
mstat_5Widow             -2.899e-01  9.143e-02  -3.171 0.001518 ** 
relig_affilNo            -9.577e-02  2.381e-02  -4.022 5.77e-05 ***
act_tobYes               -1.722e-01  3.857e-02  -4.464 8.05e-06 ***
max_ch                    2.148e-02  2.544e-03   8.441  < 2e-16 ***
IC                        3.020e-01  2.635e-02  11.462  < 2e-16 ***
pop_dens                 -1.282e-05  3.503e-06  -3.658 0.000254 ***
r_pct                    -7.897e-03  7.438e-04 -10.617  < 2e-16 ***
RPL_THEMES               -4.568e-01  4.780e-02  -9.556  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

    Null deviance: 8459.2  on 6407  degrees of freedom
Residual deviance: 7817.7  on 6386  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 27267

Number of Fisher Scoring iterations: 1

              Theta:  2.1619 
          Std. Err.:  0.0838 

 2 x log-likelihood:  -27221.3070 

broom::glance(totalflu_SVI.nb)

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

tbl_regression(totalflu_SVI.nb, label = list(age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
ibd_3
CD	—	—
UC	1.01	0.97, 1.06	0.6
Unspecified	0.46	0.31, 0.67	<0.001
Age	1.00	1.00, 1.00	0.8
Gender
Male	—	—
Female	1.12	1.07, 1.17	<0.001
race
White	—	—
Black	0.80	0.72, 0.88	<0.001
Asian	1.16	1.02, 1.32	0.027
Native	1.13	0.79, 1.62	0.5
Other	0.81	0.71, 0.92	0.002
Ethnicity
NonHispanic	—	—
Hispanic	1.22	1.05, 1.43	0.011
Preferred Language
English	—	—
Other	0.64	0.50, 0.83	<0.001
Marital Status
Married	—	—
Unknown	0.76	0.71, 0.81	<0.001
Unmarried	0.87	0.81, 0.92	<0.001
DivorcedSeparated	0.90	0.78, 1.03	0.12
Widow	0.75	0.63, 0.90	0.002
Any Religious Affiliation
Yes	—	—
No	0.91	0.87, 0.95	<0.001
Active Tobacco Use
No	—	—
Yes	0.84	0.78, 0.91	<0.001
Charlson Comorbidity Index	1.02	1.02, 1.03	<0.001
IC	1.35	1.28, 1.43	<0.001
Population density	1.00	1.00, 1.00	<0.001
Percent Rebuplican	0.99	0.99, 0.99	<0.001
Total SVI	0.63	0.58, 0.70	<0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
totalflu_res.nb <- resid(totalflu_SVI.nb)
plot(fitted(totalflu_SVI.nb), totalflu_res.nb, 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(totalflu_SVI.nb) - logLik(totalflu_SVI)), df = 1, lower.tail = FALSE)

'log Lik.' 0 (df=23)

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

# Model performance 
model_performance(totalflu_SVI.nb)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
27267.307 | 27422.909 |           0.130 | 2.247 | 1.106 |    -2.124 |           0.011

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

# Margins 
cplot(totalflu_SVI.nb, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

NA
NA

SVI Quartile Poisson Model

totalflu_4 <- glm(total_flu ~ age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5
                + relig_affil+ act_tob + max_ch + pop_dens + r_pct + RPL_4,
               family = poisson,
               data = vax_clean1) 
summary(totalflu_4)
broom::glance(totalflu_4)
broom::tidy(totalflu_4, exponentiate = TRUE)
model_performance(totalflu_4)
tbl_regression(totalflu_4, label = list(age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_4 ~ "SVI Quartile", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use"), exponentiate = TRUE)

# Residual plots 
totalflu_4_res <- resid(totalflu_4)
plot(fitted(totalflu_4), totalflu_4_res, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Poisson')
abline(0,0)
# Poisson model assumes that variance is equal to the mean. If variance is significantly different than the mean, NB regression is better. Here was see that the standardized residuals fall outside the -2 to 2 range, indicating the NB regression is more appropriate. 

library(AER)
dispersiontest(totalflu_4, trafo=1)

SVI Quartile Negative Binomial Model

totalflu_SVI_4.nb <- glm.nb(total_flu ~ ibd_3 + age_yrs + gender + race_5 + 
                               ethnic_3 + lang_3  + mstat_5 + relig_affil 
                               + act_tob + max_ch + IC + r_pct + pop_dens + RPL_4,
               data = vax_clean1) 
summary(totalflu_SVI_4.nb)

Call:
glm.nb(formula = total_flu ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + r_pct + pop_dens + RPL_4, data = vax_clean1, init.theta = 2.153357942, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5988  -1.2244  -0.1288   0.5968   2.5981  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.354e+00  6.952e-02  19.473  < 2e-16 ***
ibd_3UC                   1.655e-02  2.458e-02   0.673 0.500657    
ibd_3Unspecified         -8.046e-01  1.983e-01  -4.058 4.95e-05 ***
age_yrs                   2.265e-04  8.069e-04   0.281 0.778949    
genderFemale              1.108e-01  2.317e-02   4.782 1.73e-06 ***
race_5Black              -2.394e-01  5.019e-02  -4.771 1.83e-06 ***
race_5Asian               1.535e-01  6.741e-02   2.277 0.022801 *  
race_5Native              1.137e-01  1.823e-01   0.624 0.532809    
race_5Other              -2.166e-01  6.881e-02  -3.147 0.001648 ** 
ethnic_3Hispanic          1.972e-01  7.995e-02   2.467 0.013623 *  
lang_3Other              -4.502e-01  1.282e-01  -3.510 0.000447 ***
mstat_5Unknown           -2.797e-01  3.337e-02  -8.380  < 2e-16 ***
mstat_5Unmarried         -1.466e-01  3.088e-02  -4.748 2.06e-06 ***
mstat_5DivorcedSeparated -1.154e-01  7.144e-02  -1.615 0.106238    
mstat_5Widow             -2.930e-01  9.153e-02  -3.202 0.001366 ** 
relig_affilNo            -9.551e-02  2.384e-02  -4.007 6.15e-05 ***
act_tobYes               -1.764e-01  3.857e-02  -4.575 4.77e-06 ***
max_ch                    2.143e-02  2.548e-03   8.411  < 2e-16 ***
IC                        2.996e-01  2.638e-02  11.358  < 2e-16 ***
r_pct                    -7.863e-03  7.470e-04 -10.526  < 2e-16 ***
pop_dens                 -1.311e-05  3.521e-06  -3.724 0.000196 ***
RPL_4Second              -1.441e-01  2.733e-02  -5.272 1.35e-07 ***
RPL_4Third               -1.959e-01  3.271e-02  -5.989 2.11e-09 ***
RPL_4Fourth              -3.464e-01  4.526e-02  -7.653 1.96e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

    Null deviance: 8445.5  on 6407  degrees of freedom
Residual deviance: 7814.8  on 6384  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 27281

Number of Fisher Scoring iterations: 1

              Theta:  2.1534 
          Std. Err.:  0.0832 

 2 x log-likelihood:  -27230.7460 

broom::glance(totalflu_SVI_4.nb)

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

tbl_regression(totalflu_SVI_4.nb, label = list(RPL_4 ~ "SVI Quartile", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", race_5 ~ "Race", r_pct ~ "Percent Rebuplican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", r_pct ~ "Percent Republican", ibd_3 ~ "IBD Type", IC ~ "Immunocompromised", mstat_5 ~ "Marital Status"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Type
CD	—	—
UC	1.02	0.97, 1.07	0.5
Unspecified	0.45	0.30, 0.66	<0.001
Age	1.00	1.00, 1.00	0.8
Gender
Male	—	—
Female	1.12	1.07, 1.17	<0.001
Race
White	—	—
Black	0.79	0.71, 0.87	<0.001
Asian	1.17	1.02, 1.33	0.023
Native	1.12	0.79, 1.61	0.5
Other	0.81	0.70, 0.92	0.002
Ethnicity
NonHispanic	—	—
Hispanic	1.22	1.04, 1.42	0.014
Preferred Language
English	—	—
Other	0.64	0.50, 0.82	<0.001
Marital Status
Married	—	—
Unknown	0.76	0.71, 0.81	<0.001
Unmarried	0.86	0.81, 0.92	<0.001
DivorcedSeparated	0.89	0.77, 1.03	0.11
Widow	0.75	0.62, 0.89	0.001
Any Religious Affiliation
Yes	—	—
No	0.91	0.87, 0.95	<0.001
Active Tobacco Use
No	—	—
Yes	0.84	0.78, 0.90	<0.001
Charlson Comorbidity Index	1.02	1.02, 1.03	<0.001
Immunocompromised	1.35	1.28, 1.42	<0.001
Percent Republican	0.99	0.99, 0.99	<0.001
Population density	1.00	1.00, 1.00	<0.001
SVI Quartile
First	—	—
Second	0.87	0.82, 0.91	<0.001
Third	0.82	0.77, 0.88	<0.001
Fourth	0.71	0.65, 0.77	<0.001
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
totalflu_SVI_4.nb_res <- resid(totalflu_SVI_4.nb)
plot(fitted(totalflu_SVI_4.nb), totalflu_SVI_4.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(totalflu_SVI_4.nb) - logLik(totalflu_SVI_4.nb)), df = 1, lower.tail = FALSE)

'log Lik.' 1 (df=25)

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

# Model performance 
model_performance(totalflu_SVI_4.nb)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
27280.746 | 27449.879 |           0.128 | 2.249 | 1.106 |    -2.125 |           0.011

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

# Margins 
cplot(totalflu_SVI_4.nb, "RPL_4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI Quartile")

All Themes Poisson Model

totalflu_themes <- glm(total_flu  ~ age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5
                + relig_affil + act_tob + max_ch + pop_dens + r_pct + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
               family = poisson,
               data = vax_clean1) 
summary(totalflu_themes)
broom::glance(totalflu_themes)
broom::tidy(totalflu_themes, exponentiate = TRUE)
model_performance(totalflu_themes)
tbl_regression(totalflu_themes, label = list(age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use"), exponentiate = TRUE)

# Residual plots 
totalflu_res <- resid(totalflu_themes)
plot(fitted(totalflu_themes), totalflu_res, col='steelblue', pch=16,
     xlab='Predicted Vaccines', ylab='Standardized Residuals', main='Poisson')
abline(0,0)
# Poisson model assumes that variance is equal to the mean. If variance is significantly different than the mean, NB regression is better. Here was see that the standardized residuals fall outside the -2 to 2 range, indicating the NB regression is more appropriate. 

library(AER)
dispersiontest(totalflu_SVI, trafo=1)

All Themes Negative Binomial Model

totalflu_themes.nb <- glm.nb(total_flu ~ ibd_3 + age_yrs + gender + race_5 + 
                               ethnic_3 + lang_3  + mstat_5 + relig_affil 
                               + act_tob + max_ch + IC + r_pct + pop_dens +
                               RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4, 
               data = vax_clean1) 
summary(totalflu_themes.nb)

Call:
glm.nb(formula = total_flu ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + r_pct + pop_dens + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, data = vax_clean1, init.theta = 2.181681817, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6324  -1.2325  -0.1297   0.5927   2.6832  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.340e+00  8.292e-02  16.162  < 2e-16 ***
ibd_3UC                   1.203e-02  2.450e-02   0.491 0.623433    
ibd_3Unspecified         -7.880e-01  1.982e-01  -3.977 6.99e-05 ***
age_yrs                   3.360e-04  8.045e-04   0.418 0.676159    
genderFemale              1.118e-01  2.310e-02   4.839 1.30e-06 ***
race_5Black              -2.113e-01  5.028e-02  -4.203 2.64e-05 ***
race_5Asian               1.166e-01  6.768e-02   1.723 0.084897 .  
race_5Native              1.344e-01  1.812e-01   0.742 0.458126    
race_5Other              -2.185e-01  6.864e-02  -3.183 0.001459 ** 
ethnic_3Hispanic          1.926e-01  7.971e-02   2.416 0.015705 *  
lang_3Other              -4.632e-01  1.280e-01  -3.618 0.000297 ***
mstat_5Unknown           -2.756e-01  3.330e-02  -8.278  < 2e-16 ***
mstat_5Unmarried         -1.437e-01  3.081e-02  -4.664 3.09e-06 ***
mstat_5DivorcedSeparated -1.078e-01  7.128e-02  -1.513 0.130259    
mstat_5Widow             -2.897e-01  9.119e-02  -3.177 0.001490 ** 
relig_affilNo            -9.473e-02  2.383e-02  -3.976 7.00e-05 ***
act_tobYes               -1.721e-01  3.854e-02  -4.467 7.95e-06 ***
max_ch                    2.160e-02  2.538e-03   8.512  < 2e-16 ***
IC                        3.010e-01  2.632e-02  11.439  < 2e-16 ***
r_pct                    -6.007e-03  8.895e-04  -6.753 1.44e-11 ***
pop_dens                 -1.232e-05  3.522e-06  -3.498 0.000469 ***
RPL_THEME1               -2.292e-01  6.739e-02  -3.401 0.000671 ***
RPL_THEME2               -2.930e-01  5.944e-02  -4.929 8.27e-07 ***
RPL_THEME3                3.503e-02  4.815e-02   0.727 0.466963    
RPL_THEME4               -5.557e-02  5.035e-02  -1.104 0.269790    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

    Null deviance: 8490.6  on 6407  degrees of freedom
Residual deviance: 7820.8  on 6383  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 27248

Number of Fisher Scoring iterations: 1

              Theta:  2.1817 
          Std. Err.:  0.0849 

 2 x log-likelihood:  -27196.1650 

broom::glance(totalflu_themes.nb)

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

model_performance(totalflu_themes.nb)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
27248.165 | 27424.063 |           0.135 | 2.243 | 1.107 |    -2.122 |           0.011

tbl_regression(totalflu_themes.nb, label = list(age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", race_5 ~ "Race", r_pct ~ "Percent Rebuplican", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", r_pct ~ "Percent Republican", ibd_3 ~ "IBD Type", IC ~ "Immunocompromised", mstat_5 ~ "Marital Status"), exponentiate = TRUE)

Characteristic	IRR1	95% CI1	p-value
IBD Type
CD	—	—
UC	1.01	0.96, 1.06	0.6
Unspecified	0.45	0.31, 0.67	<0.001
Age	1.00	1.00, 1.00	0.7
Gender
Male	—	—
Female	1.12	1.07, 1.17	<0.001
Race
White	—	—
Black	0.81	0.73, 0.89	<0.001
Asian	1.12	0.98, 1.28	0.085
Native	1.14	0.80, 1.64	0.5
Other	0.80	0.70, 0.92	0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.21	1.04, 1.42	0.016
Preferred Language
English	—	—
Other	0.63	0.49, 0.81	<0.001
Marital Status
Married	—	—
Unknown	0.76	0.71, 0.81	<0.001
Unmarried	0.87	0.82, 0.92	<0.001
DivorcedSeparated	0.90	0.78, 1.03	0.13
Widow	0.75	0.63, 0.90	0.001
Any Religious Affiliation
Yes	—	—
No	0.91	0.87, 0.95	<0.001
Active Tobacco Use
No	—	—
Yes	0.84	0.78, 0.91	<0.001
Charlson Comorbidity Index	1.02	1.02, 1.03	<0.001
Immunocompromised	1.35	1.28, 1.42	<0.001
Percent Republican	0.99	0.99, 1.00	<0.001
Population density	1.00	1.00, 1.00	<0.001
Soceioeconomic Status	0.80	0.70, 0.91	<0.001
Household Composition	0.75	0.66, 0.84	<0.001
Minority Status and Language	1.04	0.94, 1.14	0.5
Housing and Transportation	0.95	0.86, 1.04	0.3
1 IRR = Incidence Rate Ratio, CI = Confidence Interval

# NB Residual Plot
totalflu_res.nb <- resid(totalflu_themes.nb)
plot(fitted(totalflu_themes.nb), totalflu_res.nb, 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(totalflu_themes.nb) - logLik(totalflu_themes.nb)), df = 1, lower.tail = FALSE)

'log Lik.' 1 (df=26)

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

# Model performance 
model_performance(totalflu_themes.nb)

# Indices of model performance

AIC       |       BIC | Nagelkerke's R2 |  RMSE | Sigma | Score_log | Score_spherical
-------------------------------------------------------------------------------------
27248.165 | 27424.063 |           0.135 | 2.243 | 1.107 |    -2.122 |           0.011

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

# Margins 
cplot(totalflu_themes.nb, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME1")

cplot(totalflu_themes.nb, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME2")

cplot(totalflu_themes.nb, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME3")

cplot(totalflu_themes.nb, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME4")

Flu 2019 + Total SVI: Multivaraible Model

All variables

Flu2019_SVI <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEMES,
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_SVI )

Call:
glm(formula = flu_2019 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_THEMES, family = "binomial", 
    data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9215  -1.0401  -0.7805   1.1957   2.1263  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               8.637e-01  1.628e-01   5.306 1.12e-07 ***
ibd_3UC                   5.478e-02  5.618e-02   0.975 0.329508    
ibd_3Unspecified         -8.824e-01  4.354e-01  -2.027 0.042713 *  
age_yrs                  -1.024e-02  1.855e-03  -5.521 3.38e-08 ***
genderFemale              1.813e-01  5.281e-02   3.433 0.000596 ***
race_5Black              -5.902e-01  1.150e-01  -5.133 2.85e-07 ***
race_5Asian               3.933e-01  1.587e-01   2.478 0.013197 *  
race_5Native              8.523e-02  4.273e-01   0.199 0.841898    
race_5Other              -4.284e-01  1.552e-01  -2.760 0.005787 ** 
ethnic_3Hispanic          4.062e-01  1.852e-01   2.193 0.028282 *  
lang_3Other              -7.432e-01  2.924e-01  -2.542 0.011031 *  
mstat_5Unknown           -3.358e-01  7.544e-02  -4.452 8.51e-06 ***
mstat_5Unmarried         -3.921e-01  7.128e-02  -5.501 3.77e-08 ***
mstat_5DivorcedSeparated -1.565e-01  1.633e-01  -0.958 0.337854    
mstat_5Widow             -6.860e-01  2.244e-01  -3.057 0.002232 ** 
relig_affilNo            -2.040e-01  5.421e-02  -3.764 0.000167 ***
act_tobYes               -1.906e-01  8.595e-02  -2.218 0.026548 *  
max_ch                    2.549e-02  6.003e-03   4.246 2.18e-05 ***
IC                        6.439e-01  6.050e-02  10.643  < 2e-16 ***
pop_dens                 -2.822e-05  8.967e-06  -3.147 0.001649 ** 
r_pct                    -1.396e-02  1.748e-03  -7.988 1.37e-15 ***
RPL_THEMES               -7.494e-01  1.086e-01  -6.901 5.18e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8347.1  on 6386  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8391.1

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_SVI )

broom::tidy(Flu2019_SVI , exponentiate = TRUE)

tbl_regression(Flu2019_SVI, label = list(ibd_3 ~ "IBD Type", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised", RPL_THEMES ~ "Total SVI"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Type
CD	—	—
UC	1.06	0.95, 1.18	0.3
Unspecified	0.41	0.16, 0.92	0.043
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.20	1.08, 1.33	<0.001
Race
White	—	—
Black	0.55	0.44, 0.69	<0.001
Asian	1.48	1.09, 2.03	0.013
Native	1.09	0.46, 2.52	0.8
Other	0.65	0.48, 0.88	0.006
Ethnicity
NonHispanic	—	—
Hispanic	1.50	1.04, 2.16	0.028
Preferred Language
English	—	—
Other	0.48	0.26, 0.83	0.011
Marital Status
Married	—	—
Unknown	0.71	0.62, 0.83	<0.001
Unmarried	0.68	0.59, 0.78	<0.001
DivorcedSeparated	0.86	0.62, 1.18	0.3
Widow	0.50	0.32, 0.77	0.002
Any Religious Affiliation
Yes	—	—
No	0.82	0.73, 0.91	<0.001
Active Tobacco Use
No	—	—
Yes	0.83	0.70, 0.98	0.027
Charlson Comorbidity Index	1.03	1.01, 1.04	<0.001
Immunocompromised	1.90	1.69, 2.14	<0.001
Population density	1.00	1.00, 1.00	0.002
Percent Rebuplican	0.99	0.98, 0.99	<0.001
Total SVI	0.47	0.38, 0.58	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(Flu2019_SVI)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8391.150 | 8539.986 |     0.059 | 0.479 | 1.143 |    0.651 |      -Inf |       1.561e-04 | 0.540

performance::check_model(Flu2019_SVI, panel = TRUE)

# Margins 
cplot(Flu2019_SVI, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

Age and RPL_THEMES as Quadratic

Flu2019_SVI.2 <- glm(flu_2019 ~ I(age_yrs^2) + ibd_3  + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct +  I(RPL_THEMES^2),
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_SVI.2 )

Call:
glm(formula = flu_2019 ~ I(age_yrs^2) + ibd_3 + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + I(RPL_THEMES^2), family = "binomial", 
    data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9693  -1.0458  -0.7808   1.1944   2.2159  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               5.890e-01  1.397e-01   4.216 2.49e-05 ***
I(age_yrs^2)             -1.201e-04  1.902e-05  -6.314 2.72e-10 ***
ibd_3UC                   5.445e-02  5.620e-02   0.969 0.332605    
ibd_3Unspecified         -9.000e-01  4.344e-01  -2.072 0.038291 *  
genderFemale              1.741e-01  5.278e-02   3.298 0.000974 ***
race_5Black              -5.927e-01  1.157e-01  -5.123 3.01e-07 ***
race_5Asian               3.867e-01  1.586e-01   2.438 0.014776 *  
race_5Native              6.873e-02  4.267e-01   0.161 0.872021    
race_5Other              -4.216e-01  1.552e-01  -2.716 0.006604 ** 
ethnic_3Hispanic          4.143e-01  1.850e-01   2.239 0.025158 *  
lang_3Other              -7.430e-01  2.926e-01  -2.539 0.011111 *  
mstat_5Unknown           -3.359e-01  7.447e-02  -4.511 6.46e-06 ***
mstat_5Unmarried         -3.934e-01  6.927e-02  -5.680 1.35e-08 ***
mstat_5DivorcedSeparated -1.712e-01  1.636e-01  -1.047 0.295158    
mstat_5Widow             -5.852e-01  2.265e-01  -2.584 0.009761 ** 
relig_affilNo            -2.159e-01  5.431e-02  -3.974 7.06e-05 ***
act_tobYes               -2.062e-01  8.589e-02  -2.400 0.016384 *  
max_ch                    2.801e-02  6.068e-03   4.615 3.92e-06 ***
IC                        6.364e-01  6.040e-02  10.537  < 2e-16 ***
pop_dens                 -2.963e-05  9.136e-06  -3.244 0.001180 ** 
r_pct                    -1.404e-02  1.755e-03  -8.001 1.23e-15 ***
I(RPL_THEMES^2)          -7.787e-01  1.243e-01  -6.264 3.74e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8346.6  on 6386  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8390.6

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_SVI.2 )

broom::tidy(Flu2019_SVI.2 , exponentiate = TRUE)

tbl_regression(Flu2019_SVI.2, label = list(gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
I(age_yrs^2)	1.00	1.00, 1.00	<0.001
ibd_3
CD	—	—
UC	1.06	0.95, 1.18	0.3
Unspecified	0.41	0.16, 0.90	0.038
Gender
Male	—	—
Female	1.19	1.07, 1.32	<0.001
Race
White	—	—
Black	0.55	0.44, 0.69	<0.001
Asian	1.47	1.08, 2.01	0.015
Native	1.07	0.46, 2.48	0.9
Other	0.66	0.48, 0.89	0.007
Ethnicity
NonHispanic	—	—
Hispanic	1.51	1.05, 2.18	0.025
Preferred Language
English	—	—
Other	0.48	0.26, 0.83	0.011
Marital Status
Married	—	—
Unknown	0.71	0.62, 0.83	<0.001
Unmarried	0.67	0.59, 0.77	<0.001
DivorcedSeparated	0.84	0.61, 1.16	0.3
Widow	0.56	0.35, 0.86	0.010
Any Religious Affiliation
Yes	—	—
No	0.81	0.72, 0.90	<0.001
act_tob
No	—	—
Yes	0.81	0.69, 0.96	0.016
Charlson Comorbidity Index	1.03	1.02, 1.04	<0.001
IC	1.89	1.68, 2.13	<0.001
Population Density	1.00	1.00, 1.00	0.001
Percent Republican	0.99	0.98, 0.99	<0.001
I(RPL_THEMES^2)	0.46	0.36, 0.59	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

## Standard Residuals 
Flu2019_SVI.2.std <- rstandard(Flu2019_SVI.2)
plot(Flu2019_SVI.2.std, ylab="Standardized Residuals")

## Margins 
cplot(Flu2019_SVI.2, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

## Model performance graphs
model_performance(Flu2019_SVI.2 )

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
--------------------------------------------------------------------------------------------------
1.787e+05 | 1.789e+05 |     0.059 | 0.479 | 1.143 |    0.651 |      -Inf |       1.561e-04 | 0.540

performance::check_model(Flu2019_SVI.2, panel = TRUE)

Interaction between race and RPL_THEMES

Flu2019_SVI.mod <- glm(flu_2019 ~ I(age_yrs^2) + ibd_3 + age_yrs + gender  + 
                         ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + 
                         max_ch + IC + pop_dens + r_pct + RPL_THEMES +
                         I(RPL_THEMES^2)*race_5,
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_SVI.mod )

Call:
glm(formula = flu_2019 ~ I(age_yrs^2) + ibd_3 + age_yrs + gender + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_THEMES + I(RPL_THEMES^2) * race_5, 
    family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9591  -1.0391  -0.7745   1.1986   2.1004  

Coefficients:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   3.675e-01  2.190e-01   1.678 0.093368 .  
I(age_yrs^2)                 -2.852e-04  7.359e-05  -3.876 0.000106 ***
ibd_3UC                       5.511e-02  5.635e-02   0.978 0.328102    
ibd_3Unspecified             -9.000e-01  4.358e-01  -2.065 0.038885 *  
age_yrs                       1.674e-02  7.194e-03   2.327 0.019948 *  
genderFemale                  1.656e-01  5.304e-02   3.122 0.001798 ** 
ethnic_3Hispanic              4.080e-01  1.861e-01   2.192 0.028384 *  
lang_3Other                  -7.392e-01  2.943e-01  -2.512 0.012007 *  
mstat_5Unknown               -2.903e-01  7.634e-02  -3.802 0.000143 ***
mstat_5Unmarried             -3.377e-01  7.263e-02  -4.650 3.32e-06 ***
mstat_5DivorcedSeparated     -1.638e-01  1.639e-01  -0.999 0.317641    
mstat_5Widow                 -4.862e-01  2.311e-01  -2.104 0.035393 *  
relig_affilNo                -2.216e-01  5.446e-02  -4.069 4.71e-05 ***
act_tobYes                   -2.106e-01  8.619e-02  -2.443 0.014568 *  
max_ch                        2.912e-02  6.093e-03   4.778 1.77e-06 ***
IC                            6.479e-01  6.059e-02  10.693  < 2e-16 ***
pop_dens                     -2.923e-05  9.111e-06  -3.208 0.001338 ** 
r_pct                        -1.386e-02  1.759e-03  -7.880 3.27e-15 ***
RPL_THEMES                   -1.204e+00  3.650e-01  -3.297 0.000976 ***
I(RPL_THEMES^2)               5.044e-01  4.285e-01   1.177 0.239082    
race_5Black                  -5.835e-01  1.791e-01  -3.258 0.001122 ** 
race_5Asian                   3.050e-01  2.025e-01   1.506 0.131959    
race_5Native                  1.276e-01  5.639e-01   0.226 0.820917    
race_5Other                  -4.517e-01  1.894e-01  -2.385 0.017098 *  
I(RPL_THEMES^2):race_5Black  -7.177e-02  3.692e-01  -0.194 0.845893    
I(RPL_THEMES^2):race_5Asian   8.018e-01  1.040e+00   0.771 0.440724    
I(RPL_THEMES^2):race_5Native -2.090e-01  1.922e+00  -0.109 0.913394    
I(RPL_THEMES^2):race_5Other   1.525e-01  5.939e-01   0.257 0.797303    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8329.8  on 6380  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8385.8

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_SVI.mod )

broom::tidy(Flu2019_SVI.mod , exponentiate = TRUE)

tbl_regression(Flu2019_SVI.mod, label = list(gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
I(age_yrs^2)	1.00	1.00, 1.00	<0.001
ibd_3
CD	—	—
UC	1.06	0.95, 1.18	0.3
Unspecified	0.41	0.16, 0.91	0.039
age_yrs	1.02	1.00, 1.03	0.020
Gender
Male	—	—
Female	1.18	1.06, 1.31	0.002
Ethnicity
NonHispanic	—	—
Hispanic	1.50	1.04, 2.17	0.028
Preferred Language
English	—	—
Other	0.48	0.26, 0.84	0.012
Marital Status
Married	—	—
Unknown	0.75	0.64, 0.87	<0.001
Unmarried	0.71	0.62, 0.82	<0.001
DivorcedSeparated	0.85	0.61, 1.17	0.3
Widow	0.61	0.39, 0.96	0.035
Any Religious Affiliation
Yes	—	—
No	0.80	0.72, 0.89	<0.001
act_tob
No	—	—
Yes	0.81	0.68, 0.96	0.015
Charlson Comorbidity Index	1.03	1.02, 1.04	<0.001
IC	1.91	1.70, 2.15	<0.001
Population Density	1.00	1.00, 1.00	0.001
Percent Republican	0.99	0.98, 0.99	<0.001
RPL_THEMES	0.30	0.15, 0.61	<0.001
I(RPL_THEMES^2)	1.66	0.71, 3.83	0.2
Race
White	—	—
Black	0.56	0.39, 0.79	0.001
Asian	1.36	0.91, 2.02	0.13
Native	1.14	0.37, 3.50	0.8
Other	0.64	0.44, 0.92	0.017
I(RPL_THEMES^2) * Race
I(RPL_THEMES^2) * Black	0.93	0.45, 1.92	0.8
I(RPL_THEMES^2) * Asian	2.23	0.29, 17.6	0.4
I(RPL_THEMES^2) * Native	0.81	0.01, 36.7	>0.9
I(RPL_THEMES^2) * Other	1.16	0.35, 3.67	0.8
1 OR = Odds Ratio, CI = Confidence Interval

## Standard Residuals 
Flu2019_SVI.mod.std <- rstandard(Flu2019_SVI.mod)
plot(Flu2019_SVI.mod.std, ylab="Standardized Residuals")

## Margins 
cplot(Flu2019_SVI.mod, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

## Model performance graphs
model_performance(Flu2019_SVI.mod )

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
--------------------------------------------------------------------------------------------------
1.862e+05 | 1.864e+05 |     0.061 | 0.479 | 1.143 |    0.650 |      -Inf |       1.561e-04 | 0.541

performance::check_model(Flu2019_SVI.mod, panel = TRUE)

SVI as a mediator of race -> vaccination (bivariate)

library(mediation)
model.0 <- glm(flu_2019 ~ race_5, family = binomial, data = vax_clean1)
model.m <- lm(RPL_THEMES ~ race_5, data = vax_clean1)
model.y <- glm(flu_2019 ~ race_5 + RPL_THEMES, family = binomial, data = vax_clean1)
results1 <- mediate(model.m, model.y, treat = "race_5", mediator = "RPL_THEMES", 
                    boot = TRUE, sims = 500)
summary(results1)

SVI as a mediator of race -> vaccination

vax_clean1 %>% 
mutate(black = as_factor(race_5),
         black = fct_recode(black, No = "Other",
                  No = "Asian",
                  No = "White", Yes = "Black",
                  No = "Native"),
         black = fct_relevel(black, ref = 'No')) -> race_dich
race_dich$black<-ifelse(race_dich$black=="Yes",1,0)
race_dich %>% 
  mutate(svi_2 = case_when(RPL_THEMES>= 0.5 ~ '1',TRUE ~ "0")) ->  themes_dich
themes_dich$svi_2 = as.numeric(themes_dich$svi_2)
#Step 1: Is there a relationship between X(race) and Y(flu vaccine). Yes - black and other patients were less likely to receive the flu vaccine than Asian and white patients
Flu2019 <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + lang_3 +
                   mstat_5+ relig_affil + max_ch + IC + pop_dens + r_pct + svi_2,
                 family = binomial, data = themes_dich)

# Step 2: Is there a relationship between the mediator (SVI) and race? Yes, patients with higher SVI were more likely to be black 
Flu2019_M <- glm(svi_2 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + lang_3 +
                   mstat_5+ relig_affil + max_ch + IC + pop_dens + r_pct + flu_2019,
                 family = binomial, 
                 data = themes_dich)
summary(Flu2019_M )

Call:
glm(formula = svi_2 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + max_ch + IC + pop_dens + 
    r_pct + flu_2019, family = binomial, data = themes_dich)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5450  -0.7870  -0.6884   0.9929   2.0203  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -1.072e+00  1.740e-01  -6.160 7.29e-10 ***
ibd_3UC                  -3.294e-01  6.191e-02  -5.320 1.04e-07 ***
ibd_3Unspecified          3.496e-01  3.460e-01   1.011  0.31219    
age_yrs                   9.567e-04  2.011e-03   0.476  0.63417    
genderFemale              5.086e-02  5.839e-02   0.871  0.38380    
black                     1.383e+00  1.079e-01  12.815  < 2e-16 ***
ethnic_3Hispanic          4.339e-01  1.805e-01   2.403  0.01625 *  
lang_3Other               9.716e-01  2.513e-01   3.867  0.00011 ***
mstat_5Unknown            3.319e-01  8.241e-02   4.027 5.64e-05 ***
mstat_5Unmarried          2.323e-01  7.901e-02   2.940  0.00328 ** 
mstat_5DivorcedSeparated  7.014e-01  1.648e-01   4.255 2.09e-05 ***
mstat_5Widow              3.793e-01  2.161e-01   1.755  0.07921 .  
relig_affilNo             5.878e-02  5.977e-02   0.984  0.32536    
max_ch                    1.269e-02  6.604e-03   1.921  0.05468 .  
IC                        1.813e-02  6.665e-02   0.272  0.78561    
pop_dens                  2.614e-05  8.088e-06   3.232  0.00123 ** 
r_pct                    -3.789e-03  1.855e-03  -2.043  0.04106 *  
flu_2019                 -3.023e-01  6.039e-02  -5.006 5.55e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7618.1  on 6451  degrees of freedom
Residual deviance: 7223.2  on 6434  degrees of freedom
  (8756 observations deleted due to missingness)
AIC: 7259.2

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_M )

tbl_regression(Flu2019_M, label = list(age_yrs ~ "Age", gender~ "Gender", black ~ "Black or African American", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"))

Characteristic	log(OR)1	95% CI1	p-value
ibd_3
CD	—	—
UC	-0.33	-0.45, -0.21	<0.001
Unspecified	0.35	-0.35, 1.0	0.3
Age	0.00	0.00, 0.00	0.6
Gender
Male	—	—
Female	0.05	-0.06, 0.17	0.4
Black or African American	1.4	1.2, 1.6	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	0.43	0.07, 0.78	0.016
Preferred Language
English	—	—
Other	0.97	0.48, 1.5	<0.001
Marital Status
Married	—	—
Unknown	0.33	0.17, 0.49	<0.001
Unmarried	0.23	0.08, 0.39	0.003
DivorcedSeparated	0.70	0.37, 1.0	<0.001
Widow	0.38	-0.05, 0.80	0.079
Any Religious Affiliation
Yes	—	—
No	0.06	-0.06, 0.18	0.3
max_ch	0.01	0.00, 0.03	0.055
IC	0.02	-0.11, 0.15	0.8
Population Density	0.00	0.00, 0.00	0.001
Percent Republican	0.00	-0.01, 0.00	0.041
flu_2019	-0.30	-0.42, -0.18	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Step 3: Mediation analysis
results <- mediate(Flu2019_M, Flu2019, treat= "black", mediator="svi_2",
                   boot=TRUE, sims=500)

Running nonparametric bootstrap

summary(results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                         Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            -0.0224      -0.0313        -0.01  <2e-16 ***
ACME (treated)            -0.0184      -0.0267        -0.01  <2e-16 ***
ADE (control)             -0.1468      -0.1950        -0.10  <2e-16 ***
ADE (treated)             -0.1429      -0.1898        -0.10  <2e-16 ***
Total Effect              -0.1653      -0.2095        -0.12  <2e-16 ***
Prop. Mediated (control)   0.1354       0.0792         0.22  <2e-16 ***
Prop. Mediated (treated)   0.1115       0.0607         0.19  <2e-16 ***
ACME (average)            -0.0204      -0.0290        -0.01  <2e-16 ***
ADE (average)             -0.1449      -0.1924        -0.10  <2e-16 ***
Prop. Mediated (average)   0.1235       0.0699         0.20  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 6452 


Simulations: 500 

Flu 2019 + SVI Quartile

All Variables

Flu2019_SVI.4 <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_4,
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_SVI.4 )

Call:
glm(formula = flu_2019 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_4, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8713  -1.0406  -0.7844   1.1962   2.0851  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               7.913e-01  1.613e-01   4.907 9.27e-07 ***
ibd_3UC                   5.978e-02  5.619e-02   1.064 0.287360    
ibd_3Unspecified         -9.123e-01  4.345e-01  -2.100 0.035766 *  
age_yrs                  -1.028e-02  1.854e-03  -5.545 2.94e-08 ***
genderFemale              1.825e-01  5.280e-02   3.456 0.000547 ***
race_5Black              -6.178e-01  1.152e-01  -5.364 8.13e-08 ***
race_5Asian               4.062e-01  1.587e-01   2.559 0.010500 *  
race_5Native              8.040e-02  4.279e-01   0.188 0.850970    
race_5Other              -4.356e-01  1.552e-01  -2.807 0.005005 ** 
ethnic_3Hispanic          3.963e-01  1.852e-01   2.140 0.032321 *  
lang_3Other              -7.620e-01  2.924e-01  -2.606 0.009170 ** 
mstat_5Unknown           -3.403e-01  7.540e-02  -4.514 6.37e-06 ***
mstat_5Unmarried         -3.947e-01  7.125e-02  -5.539 3.03e-08 ***
mstat_5DivorcedSeparated -1.653e-01  1.633e-01  -1.012 0.311482    
mstat_5Widow             -6.893e-01  2.243e-01  -3.073 0.002117 ** 
relig_affilNo            -2.032e-01  5.420e-02  -3.749 0.000177 ***
act_tobYes               -1.987e-01  8.584e-02  -2.315 0.020631 *  
max_ch                    2.555e-02  6.001e-03   4.257 2.07e-05 ***
IC                        6.411e-01  6.047e-02  10.602  < 2e-16 ***
pop_dens                 -2.869e-05  9.012e-06  -3.183 0.001456 ** 
r_pct                    -1.380e-02  1.753e-03  -7.871 3.51e-15 ***
RPL_4Second              -2.698e-01  6.236e-02  -4.327 1.51e-05 ***
RPL_4Third               -3.389e-01  7.424e-02  -4.565 4.99e-06 ***
RPL_4Fourth              -5.335e-01  1.016e-01  -5.251 1.51e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8351.5  on 6384  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8399.5

Number of Fisher Scoring iterations: 4

model_performance(Flu2019_SVI.4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8399.529 | 8561.897 |     0.058 | 0.480 | 1.144 |    0.652 |      -Inf |       1.561e-04 | 0.540

broom::glance(Flu2019_SVI.4 )

broom::tidy(Flu2019_SVI.4 , exponentiate = TRUE)

tbl_regression(Flu2019_SVI.4, label = list(ibd_3 ~ "IBD Type", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", RPL_4 ~ "SVI Quartile"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Type
CD	—	—
UC	1.06	0.95, 1.19	0.3
Unspecified	0.40	0.16, 0.89	0.036
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.20	1.08, 1.33	<0.001
race
White	—	—
Black	0.54	0.43, 0.67	<0.001
Asian	1.50	1.10, 2.05	0.010
Native	1.08	0.46, 2.51	0.9
Other	0.65	0.48, 0.87	0.005
Ethnicity
NonHispanic	—	—
Hispanic	1.49	1.03, 2.14	0.032
Preferred Language
English	—	—
Other	0.47	0.26, 0.81	0.009
Marital Status
Married	—	—
Unknown	0.71	0.61, 0.82	<0.001
Unmarried	0.67	0.59, 0.77	<0.001
DivorcedSeparated	0.85	0.61, 1.17	0.3
Widow	0.50	0.32, 0.77	0.002
Any Religious Affiliation
Yes	—	—
No	0.82	0.73, 0.91	<0.001
Active Tobacco Use
No	—	—
Yes	0.82	0.69, 0.97	0.021
Charlson Comorbidity Index	1.03	1.01, 1.04	<0.001
IC	1.90	1.69, 2.14	<0.001
Population density	1.00	1.00, 1.00	0.001
Percent Rebuplican	0.99	0.98, 0.99	<0.001
SVI Quartile
First	—	—
Second	0.76	0.68, 0.86	<0.001
Third	0.71	0.62, 0.82	<0.001
Fourth	0.59	0.48, 0.72	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

## Margins Plot
cplot(Flu2019_SVI.4, "RPL_4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI Quartile")

# Performance 
model_performance(Flu2019_SVI.4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8399.529 | 8561.897 |     0.058 | 0.480 | 1.144 |    0.652 |      -Inf |       1.561e-04 | 0.540

performance::check_model(Flu2019_SVI.4, panel = TRUE)

SVI Quartile as moderator for race -> vaccination

Flu2019_SVI.4race <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_4*race_5,
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_SVI.4race )

Call:
glm(formula = flu_2019 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + RPL_4 * race_5, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6830  -1.0447  -0.7777   1.1921   2.1242  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               7.900e-01  1.621e-01   4.874 1.09e-06 ***
ibd_3UC                   6.005e-02  5.631e-02   1.066 0.286284    
ibd_3Unspecified         -9.029e-01  4.350e-01  -2.076 0.037927 *  
age_yrs                  -1.046e-02  1.859e-03  -5.626 1.84e-08 ***
genderFemale              1.846e-01  5.287e-02   3.492 0.000479 ***
ethnic_3Hispanic          3.894e-01  1.869e-01   2.083 0.037229 *  
lang_3Other              -7.340e-01  2.962e-01  -2.478 0.013200 *  
mstat_5Unknown           -3.386e-01  7.557e-02  -4.481 7.43e-06 ***
mstat_5Unmarried         -3.983e-01  7.140e-02  -5.579 2.42e-08 ***
mstat_5DivorcedSeparated -1.898e-01  1.643e-01  -1.156 0.247857    
mstat_5Widow             -6.771e-01  2.243e-01  -3.019 0.002540 ** 
relig_affilNo            -2.082e-01  5.430e-02  -3.834 0.000126 ***
act_tobYes               -1.953e-01  8.600e-02  -2.271 0.023132 *  
max_ch                    2.569e-02  6.011e-03   4.274 1.92e-05 ***
IC                        6.432e-01  6.061e-02  10.613  < 2e-16 ***
pop_dens                 -2.877e-05  9.085e-06  -3.166 0.001543 ** 
r_pct                    -1.386e-02  1.761e-03  -7.872 3.49e-15 ***
RPL_4Second              -2.399e-01  6.587e-02  -3.642 0.000270 ***
RPL_4Third               -2.875e-01  7.919e-02  -3.630 0.000283 ***
RPL_4Fourth              -6.118e-01  1.164e-01  -5.256 1.47e-07 ***
race_5Black              -2.534e-01  2.481e-01  -1.022 0.307009    
race_5Asian               6.036e-01  2.285e-01   2.642 0.008245 ** 
race_5Native              1.700e-01  6.931e-01   0.245 0.806256    
race_5Other              -5.303e-01  2.136e-01  -2.482 0.013050 *  
RPL_4Second:race_5Black  -5.934e-01  3.366e-01  -1.763 0.077920 .  
RPL_4Third:race_5Black   -7.037e-01  3.325e-01  -2.116 0.034322 *  
RPL_4Fourth:race_5Black  -1.058e-01  3.198e-01  -0.331 0.740666    
RPL_4Second:race_5Asian  -6.452e-01  3.502e-01  -1.843 0.065384 .  
RPL_4Third:race_5Asian   -1.532e-01  4.413e-01  -0.347 0.728495    
RPL_4Fourth:race_5Asian   1.284e+01  1.812e+02   0.071 0.943491    
RPL_4Second:race_5Native -1.250e-01  9.666e-01  -0.129 0.897095    
RPL_4Third:race_5Native   2.028e-01  1.583e+00   0.128 0.898058    
RPL_4Fourth:race_5Native -4.199e-01  1.440e+00  -0.292 0.770570    
RPL_4Second:race_5Other   4.275e-01  3.624e-01   1.180 0.238115    
RPL_4Third:race_5Other   -2.775e-01  4.448e-01  -0.624 0.532728    
RPL_4Fourth:race_5Other   3.338e-01  4.797e-01   0.696 0.486586    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8332.5  on 6372  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8404.5

Number of Fisher Scoring iterations: 11

model_performance(Flu2019_SVI.4race)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8404.530 | 8648.081 |     0.060 | 0.479 | 1.144 |    0.650 |      -Inf |       1.561e-04 | 0.541

broom::glance(Flu2019_SVI.4race )

broom::tidy(Flu2019_SVI.4race , exponentiate = TRUE)

tbl_regression(Flu2019_SVI.4race, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_4 ~ "SVI Quartile", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
ibd_3
CD	—	—
UC	1.06	0.95, 1.19	0.3
Unspecified	0.41	0.16, 0.90	0.038
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.20	1.08, 1.33	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.48	1.02, 2.13	0.037
Preferred Language
English	—	—
Other	0.48	0.26, 0.84	0.013
Marital Status
Married	—	—
Unknown	0.71	0.61, 0.83	<0.001
Unmarried	0.67	0.58, 0.77	<0.001
DivorcedSeparated	0.83	0.60, 1.14	0.2
Widow	0.51	0.32, 0.78	0.003
Any Religious Affiliation
Yes	—	—
No	0.81	0.73, 0.90	<0.001
act_tob
No	—	—
Yes	0.82	0.69, 0.97	0.023
max_ch	1.03	1.01, 1.04	<0.001
IC	1.90	1.69, 2.14	<0.001
Population Density	1.00	1.00, 1.00	0.002
Percent Republican	0.99	0.98, 0.99	<0.001
SVI Quartile
First	—	—
Second	0.79	0.69, 0.89	<0.001
Third	0.75	0.64, 0.88	<0.001
Fourth	0.54	0.43, 0.68	<0.001
Race
White	—	—
Black	0.78	0.47, 1.26	0.3
Asian	1.83	1.18, 2.89	0.008
Native	1.19	0.30, 4.97	0.8
Other	0.59	0.38, 0.89	0.013
SVI Quartile * Race
Second * Black	0.55	0.28, 1.07	0.078
Third * Black	0.49	0.26, 0.95	0.034
Fourth * Black	0.90	0.48, 1.69	0.7
Second * Asian	0.52	0.26, 1.04	0.065
Third * Asian	0.86	0.36, 2.05	0.7
Fourth * Asian	377,697	0.00, NA	>0.9
Second * Native	0.88	0.12, 5.81	0.9
Third * Native	1.22	0.04, 39.0	0.9
Fourth * Native	0.66	0.02, 10.4	0.8
Second * Other	1.53	0.75, 3.12	0.2
Third * Other	0.76	0.30, 1.76	0.5
Fourth * Other	1.40	0.53, 3.51	0.5
1 OR = Odds Ratio, CI = Confidence Interval

# Performance 
performance::check_model(Flu2019_SVI.4race, panel = TRUE)

Flu 2019 + Themes: Multivaraible Model

All Variables

Flu2019_themes <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEME1 + RPL_THEME2 + RPL_THEME3
                      + RPL_THEME4,
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_themes )

Call:
glm(formula = flu_2019 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9017  -1.0388  -0.7782   1.1945   2.0910  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               7.979e-01  1.922e-01   4.151 3.31e-05 ***
ibd_3UC                   5.431e-02  5.623e-02   0.966 0.334151    
ibd_3Unspecified         -8.795e-01  4.347e-01  -2.023 0.043061 *  
age_yrs                  -1.014e-02  1.857e-03  -5.462 4.70e-08 ***
genderFemale              1.820e-01  5.285e-02   3.443 0.000575 ***
race_5Black              -5.785e-01  1.158e-01  -4.997 5.84e-07 ***
race_5Asian               3.517e-01  1.600e-01   2.197 0.027991 *  
race_5Native              9.071e-02  4.278e-01   0.212 0.832077    
race_5Other              -4.396e-01  1.554e-01  -2.829 0.004667 ** 
ethnic_3Hispanic          3.980e-01  1.854e-01   2.147 0.031818 *  
lang_3Other              -7.735e-01  2.926e-01  -2.643 0.008209 ** 
mstat_5Unknown           -3.347e-01  7.552e-02  -4.431 9.37e-06 ***
mstat_5Unmarried         -3.891e-01  7.136e-02  -5.452 4.98e-08 ***
mstat_5DivorcedSeparated -1.513e-01  1.634e-01  -0.926 0.354591    
mstat_5Widow             -6.880e-01  2.245e-01  -3.064 0.002182 ** 
relig_affilNo            -1.976e-01  5.438e-02  -3.633 0.000280 ***
act_tobYes               -1.898e-01  8.608e-02  -2.204 0.027497 *  
max_ch                    2.556e-02  6.008e-03   4.254 2.10e-05 ***
IC                        6.411e-01  6.057e-02  10.585  < 2e-16 ***
pop_dens                 -2.667e-05  9.033e-06  -2.953 0.003150 ** 
r_pct                    -1.162e-02  2.075e-03  -5.599 2.16e-08 ***
RPL_THEME1               -3.882e-01  1.546e-01  -2.512 0.012014 *  
RPL_THEME2               -3.310e-01  1.355e-01  -2.443 0.014571 *  
RPL_THEME3                5.502e-02  1.101e-01   0.500 0.617336    
RPL_THEME4               -2.061e-01  1.149e-01  -1.794 0.072851 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8339.4  on 6383  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8389.4

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_themes )

broom::tidy(Flu2019_themes , exponentiate = TRUE)

tbl_regression(Flu2019_themes, label = list(ibd_3 ~ "IBD Type", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", 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 Type
CD	—	—
UC	1.06	0.95, 1.18	0.3
Unspecified	0.41	0.16, 0.92	0.043
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.20	1.08, 1.33	<0.001
Race
White	—	—
Black	0.56	0.45, 0.70	<0.001
Asian	1.42	1.04, 1.95	0.028
Native	1.09	0.47, 2.54	0.8
Other	0.64	0.47, 0.87	0.005
Ethnicity
NonHispanic	—	—
Hispanic	1.49	1.04, 2.14	0.032
Preferred Language
English	—	—
Other	0.46	0.25, 0.80	0.008
Marital Status
Married	—	—
Unknown	0.72	0.62, 0.83	<0.001
Unmarried	0.68	0.59, 0.78	<0.001
DivorcedSeparated	0.86	0.62, 1.18	0.4
Widow	0.50	0.32, 0.77	0.002
Any Religious Affiliation
Yes	—	—
No	0.82	0.74, 0.91	<0.001
Active Tobacco Use
No	—	—
Yes	0.83	0.70, 0.98	0.027
Charlson Comorbidity Index	1.03	1.01, 1.04	<0.001
IC	1.90	1.69, 2.14	<0.001
Population density	1.00	1.00, 1.00	0.003
Percent Rebuplican	0.99	0.98, 0.99	<0.001
Soceioeconomic Status	0.68	0.50, 0.92	0.012
Household Composition	0.72	0.55, 0.94	0.015
Minority Status and Language	1.06	0.85, 1.31	0.6
Housing and Transportation	0.81	0.65, 1.02	0.073
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(Flu2019_themes)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8389.397 | 8558.529 |     0.060 | 0.479 | 1.143 |    0.651 |      -Inf |       1.561e-04 | 0.541

performance::check_model(Flu2019_themes, panel = TRUE)

# Margins 
cplot(Flu2019_themes, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME1")

cplot(Flu2019_themes, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME2")

cplot(Flu2019_themes, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME3")

cplot(Flu2019_themes, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME4")

Interaction between Theme 1 and Race

Flu2019_themes.mod <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEME1*race_5 + RPL_THEME2 + RPL_THEME3
                      + RPL_THEME4,
              family = "binomial",
              data = vax_clean1)
summary(Flu2019_themes.mod )

Call:
glm(formula = flu_2019 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + RPL_THEME1 * race_5 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9050  -1.0392  -0.7786   1.1943   2.0863  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               8.011e-01  1.927e-01   4.156 3.23e-05 ***
ibd_3UC                   5.496e-02  5.627e-02   0.977 0.328693    
ibd_3Unspecified         -8.793e-01  4.347e-01  -2.023 0.043095 *  
age_yrs                  -1.017e-02  1.858e-03  -5.475 4.38e-08 ***
genderFemale              1.815e-01  5.287e-02   3.433 0.000598 ***
ethnic_3Hispanic          3.975e-01  1.859e-01   2.138 0.032516 *  
lang_3Other              -7.853e-01  2.943e-01  -2.669 0.007613 ** 
mstat_5Unknown           -3.353e-01  7.553e-02  -4.439 9.05e-06 ***
mstat_5Unmarried         -3.891e-01  7.139e-02  -5.451 5.02e-08 ***
mstat_5DivorcedSeparated -1.520e-01  1.635e-01  -0.930 0.352525    
mstat_5Widow             -6.876e-01  2.245e-01  -3.062 0.002196 ** 
relig_affilNo            -1.977e-01  5.439e-02  -3.635 0.000278 ***
act_tobYes               -1.894e-01  8.612e-02  -2.199 0.027882 *  
max_ch                    2.563e-02  6.011e-03   4.264 2.01e-05 ***
IC                        6.413e-01  6.059e-02  10.585  < 2e-16 ***
pop_dens                 -2.666e-05  9.058e-06  -2.944 0.003243 ** 
r_pct                    -1.172e-02  2.099e-03  -5.585 2.34e-08 ***
RPL_THEME1               -3.745e-01  1.626e-01  -2.303 0.021297 *  
race_5Black              -4.985e-01  2.241e-01  -2.224 0.026125 *  
race_5Asian               3.647e-01  2.212e-01   1.649 0.099180 .  
race_5Native              5.863e-01  8.023e-01   0.731 0.464883    
race_5Other              -4.815e-01  2.252e-01  -2.138 0.032503 *  
RPL_THEME2               -3.259e-01  1.366e-01  -2.385 0.017073 *  
RPL_THEME3                5.120e-02  1.105e-01   0.463 0.643241    
RPL_THEME4               -2.100e-01  1.153e-01  -1.821 0.068649 .  
RPL_THEME1:race_5Black   -1.574e-01  3.814e-01  -0.413 0.679899    
RPL_THEME1:race_5Asian   -4.792e-02  6.888e-01  -0.070 0.944535    
RPL_THEME1:race_5Native  -1.383e+00  1.900e+00  -0.728 0.466645    
RPL_THEME1:race_5Other    1.296e-01  5.168e-01   0.251 0.802030    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8733.5  on 6407  degrees of freedom
Residual deviance: 8338.6  on 6379  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 8396.6

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_themes.mod )

broom::tidy(Flu2019_themes.mod , exponentiate = TRUE)

tbl_regression(Flu2019_themes.mod, label = list(age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
ibd_3
CD	—	—
UC	1.06	0.95, 1.18	0.3
Unspecified	0.42	0.16, 0.92	0.043
Age	0.99	0.99, 0.99	<0.001
Gender
Male	—	—
Female	1.20	1.08, 1.33	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.49	1.03, 2.15	0.033
Preferred Language
English	—	—
Other	0.46	0.25, 0.80	0.008
Marital Status
Married	—	—
Unknown	0.72	0.62, 0.83	<0.001
Unmarried	0.68	0.59, 0.78	<0.001
DivorcedSeparated	0.86	0.62, 1.18	0.4
Widow	0.50	0.32, 0.77	0.002
Any Religious Affiliation
Yes	—	—
No	0.82	0.74, 0.91	<0.001
act_tob
No	—	—
Yes	0.83	0.70, 0.98	0.028
Charlson Comorbidity Index	1.03	1.01, 1.04	<0.001
IC	1.90	1.69, 2.14	<0.001
Population density	1.00	1.00, 1.00	0.003
Percent Rebuplican	0.99	0.98, 0.99	<0.001
Soceioeconomic Status	0.69	0.50, 0.95	0.021
Race
White	—	—
Black	0.61	0.39, 0.94	0.026
Asian	1.44	0.94, 2.23	0.10
Native	1.80	0.39, 9.67	0.5
Other	0.62	0.40, 0.96	0.033
Household Composition	0.72	0.55, 0.94	0.017
Minority Status and Language	1.05	0.85, 1.31	0.6
Housing and Transportation	0.81	0.65, 1.02	0.069
Soceioeconomic Status * Race
Soceioeconomic Status * Black	0.85	0.40, 1.80	0.7
Soceioeconomic Status * Asian	0.95	0.24, 3.65	>0.9
Soceioeconomic Status * Native	0.25	0.00, 8.73	0.5
Soceioeconomic Status * Other	1.14	0.41, 3.10	0.8
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(Flu2019_themes.mod)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
8396.591 | 8592.785 |     0.060 | 0.479 | 1.143 |    0.651 |      -Inf |       1.561e-04 | 0.541

performance::check_model(Flu2019_themes.mod, panel = TRUE)

# Margins 
cplot(Flu2019_themes.mod, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME1")

cplot(Flu2019_themes.mod, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME2")

cplot(Flu2019_themes.mod, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME3")

cplot(Flu2019_themes.mod, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME4")

Theme 1 as mediator of race -> vaccination

#Step 1: Is there a relationship between X(race) and Y(flu vaccine). Yes - black and other patients were less likely to receive the flu vaccine than Asian and white patients


race_dich %>% 
  mutate(theme1_2 = case_when(RPL_THEME1>= 0.5 ~ '1',TRUE ~ "0")) ->  themes_dich2
themes_dich2$theme1_2 = as.numeric(themes_dich2$theme1_2)
Flu2019.0.2 <- glm(flu_2019 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + theme1_2, family = binomial, data = themes_dich2)

# Step 2: Is there a relationship between the mediator (SVI) and race? Yes, patients with higher SVI were more likely to be black 
Flu2019_M.2 <- glm(theme1_2 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + flu_2019,
                family = binomial, 
              data = themes_dich2)
summary(Flu2019_M )

Call:
glm(formula = svi_2 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + max_ch + IC + pop_dens + 
    r_pct + flu_2019, family = binomial, data = themes_dich)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5450  -0.7870  -0.6884   0.9929   2.0203  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -1.072e+00  1.740e-01  -6.160 7.29e-10 ***
ibd_3UC                  -3.294e-01  6.191e-02  -5.320 1.04e-07 ***
ibd_3Unspecified          3.496e-01  3.460e-01   1.011  0.31219    
age_yrs                   9.567e-04  2.011e-03   0.476  0.63417    
genderFemale              5.086e-02  5.839e-02   0.871  0.38380    
black                     1.383e+00  1.079e-01  12.815  < 2e-16 ***
ethnic_3Hispanic          4.339e-01  1.805e-01   2.403  0.01625 *  
lang_3Other               9.716e-01  2.513e-01   3.867  0.00011 ***
mstat_5Unknown            3.319e-01  8.241e-02   4.027 5.64e-05 ***
mstat_5Unmarried          2.323e-01  7.901e-02   2.940  0.00328 ** 
mstat_5DivorcedSeparated  7.014e-01  1.648e-01   4.255 2.09e-05 ***
mstat_5Widow              3.793e-01  2.161e-01   1.755  0.07921 .  
relig_affilNo             5.878e-02  5.977e-02   0.984  0.32536    
max_ch                    1.269e-02  6.604e-03   1.921  0.05468 .  
IC                        1.813e-02  6.665e-02   0.272  0.78561    
pop_dens                  2.614e-05  8.088e-06   3.232  0.00123 ** 
r_pct                    -3.789e-03  1.855e-03  -2.043  0.04106 *  
flu_2019                 -3.023e-01  6.039e-02  -5.006 5.55e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7618.1  on 6451  degrees of freedom
Residual deviance: 7223.2  on 6434  degrees of freedom
  (8756 observations deleted due to missingness)
AIC: 7259.2

Number of Fisher Scoring iterations: 4

summary(Flu2019_M.2 )

Call:
glm(formula = theme1_2 ~ ibd_3 + age_yrs + gender + black + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + flu_2019, family = binomial, data = themes_dich2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.7049  -0.7746  -0.6713   1.0204   2.0272  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -1.472e+00  1.827e-01  -8.057 7.80e-16 ***
ibd_3UC                  -2.901e-01  6.269e-02  -4.627 3.71e-06 ***
ibd_3Unspecified          2.641e-02  3.745e-01   0.071 0.943794    
age_yrs                   7.131e-04  2.036e-03   0.350 0.726181    
genderFemale              4.546e-03  5.900e-02   0.077 0.938587    
black                     1.435e+00  1.091e-01  13.149  < 2e-16 ***
ethnic_3Hispanic          4.301e-01  1.833e-01   2.346 0.018953 *  
lang_3Other               1.089e+00  2.538e-01   4.291 1.78e-05 ***
mstat_5Unknown            2.459e-01  8.319e-02   2.956 0.003113 ** 
mstat_5Unmarried          1.218e-01  7.983e-02   1.526 0.127104    
mstat_5DivorcedSeparated  5.944e-01  1.667e-01   3.566 0.000362 ***
mstat_5Widow              7.570e-02  2.277e-01   0.332 0.739552    
relig_affilNo             6.051e-02  6.044e-02   1.001 0.316763    
act_tobYes                5.924e-01  8.564e-02   6.918 4.60e-12 ***
max_ch                    1.356e-02  6.646e-03   2.040 0.041344 *  
IC                        6.700e-02  6.743e-02   0.994 0.320388    
pop_dens                  5.239e-05  1.070e-05   4.898 9.70e-07 ***
r_pct                     3.460e-03  1.990e-03   1.739 0.082038 .  
flu_2019                 -3.593e-01  6.120e-02  -5.871 4.34e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7527.9  on 6407  degrees of freedom
Residual deviance: 7089.2  on 6389  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7127.2

Number of Fisher Scoring iterations: 4

broom::glance(Flu2019_M.2 )

tbl_regression(Flu2019_M.2, label = list(age_yrs ~ "Age", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"))

Characteristic	log(OR)1	95% CI1	p-value
ibd_3
CD	—	—
UC	-0.29	-0.41, -0.17	<0.001
Unspecified	0.03	-0.75, 0.74	>0.9
Age	0.00	0.00, 0.00	0.7
Gender
Male	—	—
Female	0.00	-0.11, 0.12	>0.9
black	1.4	1.2, 1.6	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	0.43	0.06, 0.78	0.019
Preferred Language
English	—	—
Other	1.1	0.59, 1.6	<0.001
Marital Status
Married	—	—
Unknown	0.25	0.08, 0.41	0.003
Unmarried	0.12	-0.03, 0.28	0.13
DivorcedSeparated	0.59	0.26, 0.92	<0.001
Widow	0.08	-0.38, 0.51	0.7
Any Religious Affiliation
Yes	—	—
No	0.06	-0.06, 0.18	0.3
act_tob
No	—	—
Yes	0.59	0.42, 0.76	<0.001
max_ch	0.01	0.00, 0.03	0.041
IC	0.07	-0.06, 0.20	0.3
Population Density	0.00	0.00, 0.00	<0.001
Percent Republican	0.00	0.00, 0.01	0.082
flu_2019	-0.36	-0.48, -0.24	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Step 3: Mediation analysis
results1 <- mediate(Flu2019_M.2, Flu2019.0.2, treat="black", mediator="theme1_2",
                   boot=TRUE, sims=500)

Running nonparametric bootstrap

summary(results1)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                         Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            -0.0261      -0.0363        -0.02  <2e-16 ***
ACME (treated)            -0.0217      -0.0315        -0.01  <2e-16 ***
ADE (control)             -0.1408      -0.1887        -0.09  <2e-16 ***
ADE (treated)             -0.1364      -0.1830        -0.09  <2e-16 ***
Total Effect              -0.1625      -0.2082        -0.12  <2e-16 ***
Prop. Mediated (control)   0.1605       0.0993         0.26  <2e-16 ***
Prop. Mediated (treated)   0.1338       0.0787         0.23  <2e-16 ***
ACME (average)            -0.0239      -0.0335        -0.02  <2e-16 ***
ADE (average)             -0.1386      -0.1859        -0.09  <2e-16 ***
Prop. Mediated (average)   0.1472       0.0883         0.25  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 6408 


Simulations: 500 

Flu 2021 + Total SVI: Multivaraible Model

All variables

Flu2021_SVI <- glm(flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEMES,
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_SVI )

Call:
glm(formula = flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_THEMES, family = "binomial", 
    data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9608  -0.8554  -0.6624   1.1839   2.3224  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.661e-01  1.793e-01   0.927 0.354162    
ibd_3UC                   4.073e-02  6.096e-02   0.668 0.504101    
ibd_3Unspecified         -1.406e+00  6.234e-01  -2.255 0.024121 *  
age_yrs                   6.469e-03  1.985e-03   3.258 0.001121 ** 
genderFemale              2.596e-01  5.775e-02   4.496 6.91e-06 ***
race_5Black              -5.379e-01  1.274e-01  -4.222 2.42e-05 ***
race_5Asian               3.752e-01  1.612e-01   2.327 0.019945 *  
race_5Native             -2.313e-01  4.801e-01  -0.482 0.629985    
race_5Other              -4.068e-01  1.746e-01  -2.329 0.019855 *  
ethnic_3Hispanic          3.722e-01  1.976e-01   1.883 0.059648 .  
lang_3Other              -8.604e-01  3.374e-01  -2.550 0.010776 *  
mstat_5Unknown           -1.564e-01  8.036e-02  -1.946 0.051702 .  
mstat_5Unmarried         -3.710e-01  7.695e-02  -4.822 1.42e-06 ***
mstat_5DivorcedSeparated -2.483e-01  1.783e-01  -1.392 0.163794    
mstat_5Widow             -4.179e-01  2.168e-01  -1.928 0.053899 .  
relig_affilNo            -1.242e-01  5.945e-02  -2.089 0.036697 *  
act_tobYes               -4.813e-01  1.032e-01  -4.665 3.08e-06 ***
max_ch                    4.440e-02  6.216e-03   7.142 9.17e-13 ***
IC                        4.128e-01  6.552e-02   6.300 2.98e-10 ***
pop_dens                 -4.865e-05  1.278e-05  -3.807 0.000141 ***
r_pct                    -2.568e-02  1.986e-03 -12.926  < 2e-16 ***
RPL_THEMES               -9.939e-01  1.214e-01  -8.186 2.71e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7318.2  on 6386  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7362.2

Number of Fisher Scoring iterations: 5

broom::glance(Flu2021_SVI )

broom::tidy(Flu2021_SVI , exponentiate = TRUE)

tbl_regression(Flu2021_SVI, label = list(ibd_3 ~ "IBD Type", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised", RPL_THEMES ~ "Total SVI"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Type
CD	—	—
UC	1.04	0.92, 1.17	0.5
Unspecified	0.25	0.06, 0.72	0.024
Age	1.01	1.00, 1.01	0.001
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Race
White	—	—
Black	0.58	0.45, 0.75	<0.001
Asian	1.46	1.06, 1.99	0.020
Native	0.79	0.29, 1.96	0.6
Other	0.67	0.47, 0.93	0.020
Ethnicity
NonHispanic	—	—
Hispanic	1.45	0.98, 2.13	0.060
Preferred Language
English	—	—
Other	0.42	0.21, 0.79	0.011
Marital Status
Married	—	—
Unknown	0.86	0.73, 1.00	0.052
Unmarried	0.69	0.59, 0.80	<0.001
DivorcedSeparated	0.78	0.55, 1.10	0.2
Widow	0.66	0.43, 1.00	0.054
Any Religious Affiliation
Yes	—	—
No	0.88	0.79, 0.99	0.037
Active Tobacco Use
No	—	—
Yes	0.62	0.50, 0.75	<0.001
Charlson Comorbidity Index	1.05	1.03, 1.06	<0.001
Immunocompromised	1.51	1.33, 1.72	<0.001
Population density	1.00	1.00, 1.00	<0.001
Percent Rebuplican	0.97	0.97, 0.98	<0.001
Total SVI	0.37	0.29, 0.47	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(Flu2021_SVI)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7362.208 | 7511.045 |     0.082 | 0.440 | 1.071 |    0.571 |      -Inf |       1.565e-04 | 0.613

performance::check_model(Flu2021_SVI, panel = TRUE)

# Margins 
cplot(Flu2021_SVI, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

Age and RPL_THEMES as Quadratic

Flu2021_SVI.2 <- glm(flu_2021 ~ I(age_yrs^2) + ibd_3  + gender + race_5 +
                       ethnic_3 + lang_3 + mstat_5 + relig_affil + 
                       act_tob + max_ch + IC + pop_dens + r_pct + I(RPL_THEMES^2),
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_SVI.2 )

Call:
glm(formula = flu_2021 ~ I(age_yrs^2) + ibd_3 + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + I(RPL_THEMES^2), family = "binomial", 
    data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9699  -0.8568  -0.6634   1.1917   2.4034  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.909e-01  1.547e-01   1.234 0.217313    
I(age_yrs^2)              6.237e-05  1.970e-05   3.166 0.001545 ** 
ibd_3UC                   3.852e-02  6.095e-02   0.632 0.527353    
ibd_3Unspecified         -1.402e+00  6.214e-01  -2.256 0.024062 *  
genderFemale              2.614e-01  5.770e-02   4.530 5.90e-06 ***
race_5Black              -5.173e-01  1.281e-01  -4.038 5.39e-05 ***
race_5Asian               3.655e-01  1.611e-01   2.269 0.023276 *  
race_5Native             -2.471e-01  4.800e-01  -0.515 0.606776    
race_5Other              -4.035e-01  1.747e-01  -2.310 0.020881 *  
ethnic_3Hispanic          3.774e-01  1.977e-01   1.909 0.056305 .  
lang_3Other              -8.693e-01  3.377e-01  -2.574 0.010052 *  
mstat_5Unknown           -1.764e-01  7.936e-02  -2.223 0.026245 *  
mstat_5Unmarried         -3.959e-01  7.480e-02  -5.292 1.21e-07 ***
mstat_5DivorcedSeparated -2.600e-01  1.785e-01  -1.457 0.145086    
mstat_5Widow             -4.566e-01  2.193e-01  -2.082 0.037356 *  
relig_affilNo            -1.228e-01  5.951e-02  -2.064 0.039018 *  
act_tobYes               -4.814e-01  1.031e-01  -4.670 3.01e-06 ***
max_ch                    4.382e-02  6.256e-03   7.005 2.47e-12 ***
IC                        4.084e-01  6.547e-02   6.238 4.42e-10 ***
pop_dens                 -4.975e-05  1.284e-05  -3.876 0.000106 ***
r_pct                    -2.576e-02  1.988e-03 -12.960  < 2e-16 ***
I(RPL_THEMES^2)          -1.097e+00  1.426e-01  -7.694 1.42e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7324.4  on 6386  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7368.4

Number of Fisher Scoring iterations: 5

broom::glance(Flu2021_SVI.2 )

broom::tidy(Flu2021_SVI.2 , exponentiate = TRUE)

tbl_regression(Flu2021_SVI.2, label = list(gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
I(age_yrs^2)	1.00	1.00, 1.00	0.002
ibd_3
CD	—	—
UC	1.04	0.92, 1.17	0.5
Unspecified	0.25	0.06, 0.72	0.024
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Race
White	—	—
Black	0.60	0.46, 0.76	<0.001
Asian	1.44	1.05, 1.97	0.023
Native	0.78	0.29, 1.93	0.6
Other	0.67	0.47, 0.93	0.021
Ethnicity
NonHispanic	—	—
Hispanic	1.46	0.98, 2.14	0.056
Preferred Language
English	—	—
Other	0.42	0.21, 0.79	0.010
Marital Status
Married	—	—
Unknown	0.84	0.72, 0.98	0.026
Unmarried	0.67	0.58, 0.78	<0.001
DivorcedSeparated	0.77	0.54, 1.09	0.15
Widow	0.63	0.41, 0.97	0.037
Any Religious Affiliation
Yes	—	—
No	0.88	0.79, 0.99	0.039
act_tob
No	—	—
Yes	0.62	0.50, 0.75	<0.001
Charlson Comorbidity Index	1.04	1.03, 1.06	<0.001
IC	1.50	1.32, 1.71	<0.001
Population Density	1.00	1.00, 1.00	<0.001
Percent Republican	0.97	0.97, 0.98	<0.001
I(RPL_THEMES^2)	0.33	0.25, 0.44	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

## Standard Residuals 
Flu2021_SVI.2.std <- rstandard(Flu2021_SVI.2)
plot(Flu2021_SVI.2.std, ylab="Standardized Residuals")

## Margins 
cplot(Flu2021_SVI.2, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

## Model performance graphs
model_performance(Flu2021_SVI.2 )

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
--------------------------------------------------------------------------------------------------
2.773e+05 | 2.774e+05 |     0.081 | 0.440 | 1.071 |    0.572 |      -Inf |       1.569e-04 | 0.612

performance::check_model(Flu2021_SVI.2, panel = TRUE)

Interaction between race and RPL_THEMES

Flu2021_SVI.mod <- glm(flu_2021 ~ I(age_yrs^2) + ibd_3 + gender + ethnic_3 
                      + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + I(RPL_THEMES^2)*race_5,
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_SVI.mod )

Call:
glm(formula = flu_2021 ~ I(age_yrs^2) + ibd_3 + gender + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + I(RPL_THEMES^2) * race_5, family = "binomial", 
    data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9555  -0.8568  -0.6648   1.1924   2.6438  

Coefficients:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   1.769e-01  1.551e-01   1.140 0.254082    
I(age_yrs^2)                  6.294e-05  1.971e-05   3.193 0.001407 ** 
ibd_3UC                       4.040e-02  6.100e-02   0.662 0.507809    
ibd_3Unspecified             -1.395e+00  6.215e-01  -2.244 0.024817 *  
genderFemale                  2.613e-01  5.772e-02   4.526 6.01e-06 ***
ethnic_3Hispanic              4.019e-01  1.996e-01   2.014 0.044027 *  
lang_3Other                  -8.361e-01  3.404e-01  -2.456 0.014057 *  
mstat_5Unknown               -1.763e-01  7.939e-02  -2.221 0.026372 *  
mstat_5Unmarried             -3.947e-01  7.486e-02  -5.273 1.34e-07 ***
mstat_5DivorcedSeparated     -2.546e-01  1.786e-01  -1.425 0.154076    
mstat_5Widow                 -4.578e-01  2.194e-01  -2.087 0.036895 *  
relig_affilNo                -1.214e-01  5.953e-02  -2.039 0.041438 *  
act_tobYes                   -4.760e-01  1.032e-01  -4.613 3.98e-06 ***
max_ch                        4.378e-02  6.259e-03   6.995 2.66e-12 ***
IC                            4.112e-01  6.551e-02   6.277 3.45e-10 ***
pop_dens                     -4.914e-05  1.282e-05  -3.834 0.000126 ***
r_pct                        -2.560e-02  1.991e-03 -12.861  < 2e-16 ***
I(RPL_THEMES^2)              -1.102e+00  1.582e-01  -6.964 3.30e-12 ***
race_5Black                  -6.011e-01  1.942e-01  -3.096 0.001964 ** 
race_5Asian                   2.816e-01  2.048e-01   1.375 0.169072    
race_5Native                  4.324e-01  7.163e-01   0.604 0.546056    
race_5Other                  -2.716e-01  2.095e-01  -1.296 0.194931    
I(RPL_THEMES^2):race_5Black   2.223e-01  4.058e-01   0.548 0.583715    
I(RPL_THEMES^2):race_5Asian   7.125e-01  1.071e+00   0.666 0.505721    
I(RPL_THEMES^2):race_5Native -6.121e+00  5.972e+00  -1.025 0.305348    
I(RPL_THEMES^2):race_5Other  -8.854e-01  8.281e-01  -1.069 0.284991    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7320.2  on 6382  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7372.2

Number of Fisher Scoring iterations: 5

broom::glance(Flu2021_SVI.mod )

broom::tidy(Flu2021_SVI.mod , exponentiate = TRUE)

tbl_regression(Flu2021_SVI.mod, label = list(gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
I(age_yrs^2)	1.00	1.00, 1.00	0.001
ibd_3
CD	—	—
UC	1.04	0.92, 1.17	0.5
Unspecified	0.25	0.06, 0.72	0.025
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.49	1.01, 2.20	0.044
Preferred Language
English	—	—
Other	0.43	0.21, 0.82	0.014
Marital Status
Married	—	—
Unknown	0.84	0.72, 0.98	0.026
Unmarried	0.67	0.58, 0.78	<0.001
DivorcedSeparated	0.78	0.54, 1.09	0.2
Widow	0.63	0.41, 0.97	0.037
Any Religious Affiliation
Yes	—	—
No	0.89	0.79, 1.00	0.041
act_tob
No	—	—
Yes	0.62	0.51, 0.76	<0.001
Charlson Comorbidity Index	1.04	1.03, 1.06	<0.001
IC	1.51	1.33, 1.72	<0.001
Population Density	1.00	1.00, 1.00	<0.001
Percent Republican	0.97	0.97, 0.98	<0.001
I(RPL_THEMES^2)	0.33	0.24, 0.45	<0.001
Race
White	—	—
Black	0.55	0.37, 0.80	0.002
Asian	1.33	0.88, 1.98	0.2
Native	1.54	0.40, 7.13	0.5
Other	0.76	0.50, 1.14	0.2
I(RPL_THEMES^2) * Race
I(RPL_THEMES^2) * Black	1.25	0.56, 2.75	0.6
I(RPL_THEMES^2) * Asian	2.04	0.23, 16.2	0.5
I(RPL_THEMES^2) * Native	0.00	0.00, 4.03	0.3
I(RPL_THEMES^2) * Other	0.41	0.07, 1.88	0.3
1 OR = Odds Ratio, CI = Confidence Interval

## Standard Residuals 
Flu2021_SVI.mod.std <- rstandard(Flu2021_SVI.mod)
plot(Flu2021_SVI.mod.std, ylab="Standardized Residuals")

## Margins 
cplot(Flu2021_SVI.mod, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI")

## Model performance graphs
model_performance(Flu2021_SVI.mod )

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
--------------------------------------------------------------------------------------------------
2.776e+05 | 2.777e+05 |     0.082 | 0.440 | 1.071 |    0.571 |      -Inf |       1.583e-04 | 0.613

performance::check_model(Flu2021_SVI.mod, panel = TRUE)

SVI as a mediator of race -> vaccination (bivariate)

model.0.2 <- glm(flu_2021 ~ race_5, family = binomial, data = vax_clean1)
model.m.2 <- lm(RPL_THEMES ~ race_5, data = vax_clean1)
model.y.2 <- glm(flu_2021 ~ race_5 + RPL_THEMES, family = binomial, data = vax_clean1)
results2 <- mediate(model.m.2, model.y.2, treat = "race_5", mediator = "RPL_THEMES", 
                    boot = TRUE, sims = 500)
summary(results2)

SVI as a mediator of race -> vaccination

#Step 1: Is there a relationship between X(race) and Y(flu vaccine). Yes - black and other patients were less likely to receive the flu vaccine than Asian and white patients
Flu2021.0 <- glm(flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + svi_2, 
                 family = binomial, 
                 data = themes_dich)

# Step 2: Is there a relationship between the mediator (SVI) and race? Yes, patients with higher SVI were more likely to be black 
Flu2021_M <- glm(svi_2 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + flu_2021,
                family = binomial, 
              data = themes_dich)
summary(Flu2021_M )

Call:
glm(formula = svi_2 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + flu_2021, family = binomial, data = themes_dich)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4316  -0.7845  -0.6761   1.0009   2.1873  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -1.083e+00  1.750e-01  -6.186 6.18e-10 ***
ibd_3UC                  -2.927e-01  6.267e-02  -4.670 3.01e-06 ***
ibd_3Unspecified          2.999e-01  3.579e-01   0.838 0.402091    
age_yrs                   1.586e-03  2.032e-03   0.781 0.434951    
genderFemale              6.913e-02  5.893e-02   1.173 0.240768    
race_5Black               1.353e+00  1.093e-01  12.380  < 2e-16 ***
race_5Asian              -4.999e-01  2.032e-01  -2.460 0.013875 *  
race_5Native             -3.435e-01  5.140e-01  -0.668 0.503943    
race_5Other               4.604e-02  1.644e-01   0.280 0.779372    
ethnic_3Hispanic          4.044e-01  1.928e-01   2.098 0.035912 *  
lang_3Other               1.070e+00  2.609e-01   4.101 4.12e-05 ***
mstat_5Unknown            3.360e-01  8.300e-02   4.048 5.15e-05 ***
mstat_5Unmarried          1.911e-01  7.956e-02   2.402 0.016312 *  
mstat_5DivorcedSeparated  5.985e-01  1.671e-01   3.582 0.000341 ***
mstat_5Widow              3.512e-01  2.194e-01   1.601 0.109375    
relig_affilNo             6.953e-02  6.022e-02   1.155 0.248283    
act_tobYes                4.715e-01  8.654e-02   5.448 5.08e-08 ***
max_ch                    1.527e-02  6.656e-03   2.294 0.021787 *  
IC                        9.783e-03  6.687e-02   0.146 0.883686    
pop_dens                  2.443e-05  7.909e-06   3.089 0.002005 ** 
r_pct                    -5.705e-03  1.898e-03  -3.006 0.002645 ** 
flu_2021                 -3.974e-01  6.805e-02  -5.839 5.24e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7566.5  on 6407  degrees of freedom
Residual deviance: 7117.2  on 6386  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7161.2

Number of Fisher Scoring iterations: 4

broom::glance(Flu2021_M )

tbl_regression(Flu2021_M, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"))

Characteristic	log(OR)1	95% CI1	p-value
ibd_3
CD	—	—
UC	-0.29	-0.42, -0.17	<0.001
Unspecified	0.30	-0.43, 0.99	0.4
Age	0.00	0.00, 0.01	0.4
Gender
Male	—	—
Female	0.07	-0.05, 0.18	0.2
Race
White	—	—
Black	1.4	1.1, 1.6	<0.001
Asian	-0.50	-0.91, -0.11	0.014
Native	-0.34	-1.5, 0.59	0.5
Other	0.05	-0.28, 0.36	0.8
Ethnicity
NonHispanic	—	—
Hispanic	0.40	0.02, 0.78	0.036
Preferred Language
English	—	—
Other	1.1	0.56, 1.6	<0.001
Marital Status
Married	—	—
Unknown	0.34	0.17, 0.50	<0.001
Unmarried	0.19	0.04, 0.35	0.016
DivorcedSeparated	0.60	0.27, 0.92	<0.001
Widow	0.35	-0.09, 0.77	0.11
Any Religious Affiliation
Yes	—	—
No	0.07	-0.05, 0.19	0.2
act_tob
No	—	—
Yes	0.47	0.30, 0.64	<0.001
max_ch	0.02	0.00, 0.03	0.022
IC	0.01	-0.12, 0.14	0.9
Population Density	0.00	0.00, 0.00	0.002
Percent Republican	-0.01	-0.01, 0.00	0.003
flu_2021	-0.40	-0.53, -0.26	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

# Step 3: Mediation analysis
results3 <- mediate(Flu2021_M, Flu2021.0, treat="race_5", mediator="svi_2",
                   boot=TRUE, sims=500)

Warning: treatment and control values do not match factor levels; using White and Black as control and treatment, respectivelyRunning nonparametric bootstrap

summary(results3)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                         Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            -0.0222      -0.0322        -0.02  <2e-16 ***
ACME (treated)            -0.0166      -0.0244        -0.01  <2e-16 ***
ADE (control)             -0.1114      -0.1532        -0.07  <2e-16 ***
ADE (treated)             -0.1058      -0.1445        -0.07  <2e-16 ***
Total Effect              -0.1280      -0.1667        -0.09  <2e-16 ***
Prop. Mediated (control)   0.1734       0.1089         0.27  <2e-16 ***
Prop. Mediated (treated)   0.1300       0.0761         0.23  <2e-16 ***
ACME (average)            -0.0194      -0.0280        -0.01  <2e-16 ***
ADE (average)             -0.1086      -0.1484        -0.07  <2e-16 ***
Prop. Mediated (average)   0.1517       0.0936         0.25  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 6408 


Simulations: 500 

Flu 2021 + SVI Quartile

All Variables

Flu2021_SVI.4 <- glm(flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_4,
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_SVI.4 )

Call:
glm(formula = flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_4, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9280  -0.8541  -0.6670   1.1948   2.3382  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               4.640e-02  1.782e-01   0.260  0.79460    
ibd_3UC                   4.442e-02  6.090e-02   0.729  0.46581    
ibd_3Unspecified         -1.433e+00  6.209e-01  -2.308  0.02101 *  
age_yrs                   6.364e-03  1.983e-03   3.209  0.00133 ** 
genderFemale              2.594e-01  5.768e-02   4.497 6.90e-06 ***
race_5Black              -5.645e-01  1.277e-01  -4.420 9.85e-06 ***
race_5Asian               3.826e-01  1.611e-01   2.375  0.01754 *  
race_5Native             -2.427e-01  4.803e-01  -0.505  0.61331    
race_5Other              -4.083e-01  1.746e-01  -2.338  0.01938 *  
ethnic_3Hispanic          3.641e-01  1.975e-01   1.843  0.06528 .  
lang_3Other              -8.814e-01  3.371e-01  -2.614  0.00894 ** 
mstat_5Unknown           -1.657e-01  8.024e-02  -2.065  0.03894 *  
mstat_5Unmarried         -3.778e-01  7.683e-02  -4.918 8.75e-07 ***
mstat_5DivorcedSeparated -2.675e-01  1.782e-01  -1.501  0.13327    
mstat_5Widow             -4.241e-01  2.165e-01  -1.959  0.05010 .  
relig_affilNo            -1.235e-01  5.938e-02  -2.079  0.03760 *  
act_tobYes               -4.955e-01  1.030e-01  -4.811 1.50e-06 ***
max_ch                    4.414e-02  6.207e-03   7.111 1.15e-12 ***
IC                        4.076e-01  6.542e-02   6.230 4.66e-10 ***
pop_dens                 -5.103e-05  1.295e-05  -3.940 8.15e-05 ***
r_pct                    -2.570e-02  1.996e-03 -12.878  < 2e-16 ***
RPL_4Second              -2.265e-01  6.734e-02  -3.363  0.00077 ***
RPL_4Third               -3.653e-01  8.194e-02  -4.458 8.27e-06 ***
RPL_4Fourth              -7.648e-01  1.188e-01  -6.438 1.21e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7333.0  on 6384  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7381

Number of Fisher Scoring iterations: 5

model_performance(Flu2021_SVI.4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7380.951 | 7543.319 |     0.080 | 0.441 | 1.072 |    0.572 |      -Inf |       1.574e-04 | 0.612

broom::glance(Flu2021_SVI.4 )

broom::tidy(Flu2021_SVI.4 , exponentiate = TRUE)

tbl_regression(Flu2021_SVI.4, label = list(ibd_3 ~ "IBD Type", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised", RPL_4 ~ "SVI Quartile"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
IBD Type
CD	—	—
UC	1.05	0.93, 1.18	0.5
Unspecified	0.24	0.06, 0.69	0.021
Age	1.01	1.00, 1.01	0.001
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Race
White	—	—
Black	0.57	0.44, 0.73	<0.001
Asian	1.47	1.07, 2.01	0.018
Native	0.78	0.29, 1.94	0.6
Other	0.66	0.47, 0.93	0.019
Ethnicity
NonHispanic	—	—
Hispanic	1.44	0.97, 2.11	0.065
Preferred Language
English	—	—
Other	0.41	0.20, 0.78	0.009
Marital Status
Married	—	—
Unknown	0.85	0.72, 0.99	0.039
Unmarried	0.69	0.59, 0.80	<0.001
DivorcedSeparated	0.77	0.54, 1.08	0.13
Widow	0.65	0.42, 0.99	0.050
Any Religious Affiliation
Yes	—	—
No	0.88	0.79, 0.99	0.038
Active Tobacco Use
No	—	—
Yes	0.61	0.50, 0.74	<0.001
Charlson Comorbidity Index	1.05	1.03, 1.06	<0.001
Immunocompromised	1.50	1.32, 1.71	<0.001
Population density	1.00	1.00, 1.00	<0.001
Percent Rebuplican	0.97	0.97, 0.98	<0.001
SVI Quartile
First	—	—
Second	0.80	0.70, 0.91	<0.001
Third	0.69	0.59, 0.81	<0.001
Fourth	0.47	0.37, 0.59	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

## Margins Plot
cplot(Flu2021_SVI.4, "RPL_4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given SVI Quartile")

# Performance 
model_performance(Flu2021_SVI.4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7380.951 | 7543.319 |     0.080 | 0.441 | 1.072 |    0.572 |      -Inf |       1.574e-04 | 0.612

performance::check_model(Flu2021_SVI.4, panel = TRUE)

SVI Quartile as moderator for race -> vaccination

Flu2021_SVI.4race <- glm(flu_2021 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_4*race_5,
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_SVI.4race )

Call:
glm(formula = flu_2021 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + RPL_4 * race_5, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9995  -0.8550  -0.6673   1.1936   2.4968  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               3.392e-02  1.790e-01   0.190  0.84970    
ibd_3UC                   4.714e-02  6.098e-02   0.773  0.43950    
ibd_3Unspecified         -1.430e+00  6.208e-01  -2.304  0.02124 *  
age_yrs                   6.397e-03  1.986e-03   3.221  0.00128 ** 
genderFemale              2.594e-01  5.773e-02   4.494 6.98e-06 ***
ethnic_3Hispanic          3.929e-01  1.993e-01   1.971  0.04872 *  
lang_3Other              -8.549e-01  3.397e-01  -2.517  0.01185 *  
mstat_5Unknown           -1.647e-01  8.035e-02  -2.050  0.04033 *  
mstat_5Unmarried         -3.760e-01  7.692e-02  -4.887 1.02e-06 ***
mstat_5DivorcedSeparated -2.728e-01  1.787e-01  -1.526  0.12689    
mstat_5Widow             -4.298e-01  2.167e-01  -1.984  0.04728 *  
relig_affilNo            -1.241e-01  5.945e-02  -2.088  0.03679 *  
act_tobYes               -4.909e-01  1.031e-01  -4.763 1.91e-06 ***
max_ch                    4.430e-02  6.216e-03   7.126 1.03e-12 ***
IC                        4.128e-01  6.551e-02   6.302 2.94e-10 ***
pop_dens                 -5.179e-05  1.305e-05  -3.968 7.24e-05 ***
r_pct                    -2.582e-02  2.005e-03 -12.881  < 2e-16 ***
RPL_4Second              -2.106e-01  7.130e-02  -2.954  0.00313 ** 
RPL_4Third               -3.252e-01  8.744e-02  -3.719  0.00020 ***
RPL_4Fourth              -7.602e-01  1.359e-01  -5.596 2.20e-08 ***
race_5Black              -3.669e-01  2.665e-01  -1.377  0.16855    
race_5Asian               3.266e-01  2.239e-01   1.458  0.14475    
race_5Native              4.527e-01  7.037e-01   0.643  0.52008    
race_5Other              -2.387e-01  2.269e-01  -1.052  0.29274    
RPL_4Second:race_5Black  -2.305e-01  3.596e-01  -0.641  0.52151    
RPL_4Third:race_5Black   -3.828e-01  3.573e-01  -1.071  0.28395    
RPL_4Fourth:race_5Black  -1.545e-01  3.574e-01  -0.432  0.66561    
RPL_4Second:race_5Asian   2.071e-01  3.502e-01   0.591  0.55437    
RPL_4Third:race_5Asian   -2.730e-01  4.626e-01  -0.590  0.55507    
RPL_4Fourth:race_5Asian   1.763e+00  1.284e+00   1.373  0.16982    
RPL_4Second:race_5Native -9.909e-01  1.096e+00  -0.904  0.36618    
RPL_4Third:race_5Native  -1.210e+01  2.291e+02  -0.053  0.95786    
RPL_4Fourth:race_5Native -1.163e+01  1.816e+02  -0.064  0.94894    
RPL_4Second:race_5Other  -4.882e-01  4.200e-01  -1.162  0.24505    
RPL_4Third:race_5Other   -1.312e-01  4.735e-01  -0.277  0.78177    
RPL_4Fourth:race_5Other  -5.747e-01  6.756e-01  -0.851  0.39502    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7323.6  on 6372  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7395.6

Number of Fisher Scoring iterations: 11

model_performance(Flu2021_SVI.4race)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7395.583 | 7639.134 |     0.081 | 0.440 | 1.072 |    0.571 |      -Inf |       3.770e-04 | 0.612

broom::glance(Flu2021_SVI.4race )

broom::tidy(Flu2021_SVI.4race , exponentiate = TRUE)

tbl_regression(Flu2021_SVI.4race, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_4 ~ "SVI Quartile", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
ibd_3
CD	—	—
UC	1.05	0.93, 1.18	0.4
Unspecified	0.24	0.06, 0.70	0.021
Age	1.01	1.00, 1.01	0.001
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.48	1.00, 2.18	0.049
Preferred Language
English	—	—
Other	0.43	0.21, 0.80	0.012
Marital Status
Married	—	—
Unknown	0.85	0.72, 0.99	0.040
Unmarried	0.69	0.59, 0.80	<0.001
DivorcedSeparated	0.76	0.53, 1.07	0.13
Widow	0.65	0.42, 0.99	0.047
Any Religious Affiliation
Yes	—	—
No	0.88	0.79, 0.99	0.037
act_tob
No	—	—
Yes	0.61	0.50, 0.75	<0.001
max_ch	1.05	1.03, 1.06	<0.001
IC	1.51	1.33, 1.72	<0.001
Population Density	1.00	1.00, 1.00	<0.001
Percent Republican	0.97	0.97, 0.98	<0.001
SVI Quartile
First	—	—
Second	0.81	0.70, 0.93	0.003
Third	0.72	0.61, 0.86	<0.001
Fourth	0.47	0.36, 0.61	<0.001
Race
White	—	—
Black	0.69	0.40, 1.15	0.2
Asian	1.39	0.89, 2.15	0.14
Native	1.57	0.39, 6.74	0.5
Other	0.79	0.50, 1.22	0.3
SVI Quartile * Race
Second * Black	0.79	0.39, 1.61	0.5
Third * Black	0.68	0.34, 1.38	0.3
Fourth * Black	0.86	0.43, 1.74	0.7
Second * Asian	1.23	0.62, 2.44	0.6
Third * Asian	0.76	0.30, 1.85	0.6
Fourth * Asian	5.83	0.50, 134	0.2
Second * Native	0.37	0.04, 2.91	0.4
Third * Native	0.00		>0.9
Fourth * Native	0.00		>0.9
Second * Other	0.61	0.26, 1.37	0.2
Third * Other	0.88	0.33, 2.15	0.8
Fourth * Other	0.56	0.12, 1.88	0.4
1 OR = Odds Ratio, CI = Confidence Interval

# Performance 
performance::check_model(Flu2021_SVI.4race, panel = TRUE)

Flu 2019 + Themes: Multivaraible Model

All Variables

Flu2021_themes <- glm(flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 +
                        RPL_THEME4,
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_themes )

Call:
glm(formula = flu_2021 ~ ibd_3 + age_yrs + gender + race_5 + 
    ethnic_3 + lang_3 + mstat_5 + relig_affil + act_tob + max_ch + 
    IC + pop_dens + r_pct + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9355  -0.8561  -0.6565   1.1821   2.3496  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               6.612e-02  2.128e-01   0.311 0.755988    
ibd_3UC                   3.806e-02  6.109e-02   0.623 0.533222    
ibd_3Unspecified         -1.391e+00  6.225e-01  -2.235 0.025450 *  
age_yrs                   6.660e-03  1.989e-03   3.348 0.000815 ***
genderFemale              2.598e-01  5.783e-02   4.493 7.01e-06 ***
race_5Black              -5.098e-01  1.283e-01  -3.973 7.09e-05 ***
race_5Asian               3.117e-01  1.626e-01   1.917 0.055192 .  
race_5Native             -2.316e-01  4.816e-01  -0.481 0.630507    
race_5Other              -4.250e-01  1.751e-01  -2.427 0.015209 *  
ethnic_3Hispanic          3.543e-01  1.980e-01   1.790 0.073452 .  
lang_3Other              -9.111e-01  3.371e-01  -2.703 0.006871 ** 
mstat_5Unknown           -1.564e-01  8.052e-02  -1.942 0.052125 .  
mstat_5Unmarried         -3.715e-01  7.711e-02  -4.818 1.45e-06 ***
mstat_5DivorcedSeparated -2.432e-01  1.787e-01  -1.361 0.173476    
mstat_5Widow             -4.194e-01  2.170e-01  -1.933 0.053263 .  
relig_affilNo            -1.190e-01  5.970e-02  -1.993 0.046230 *  
act_tobYes               -4.832e-01  1.034e-01  -4.675 2.94e-06 ***
max_ch                    4.454e-02  6.228e-03   7.152 8.53e-13 ***
IC                        4.097e-01  6.566e-02   6.240 4.37e-10 ***
pop_dens                 -4.908e-05  1.341e-05  -3.661 0.000251 ***
r_pct                    -2.187e-02  2.346e-03  -9.326  < 2e-16 ***
RPL_THEME1               -3.862e-01  1.729e-01  -2.234 0.025503 *  
RPL_THEME2               -6.723e-01  1.503e-01  -4.472 7.73e-06 ***
RPL_THEME3                6.391e-02  1.196e-01   0.534 0.593024    
RPL_THEME4               -1.870e-01  1.267e-01  -1.475 0.140132    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7302.3  on 6383  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7352.3

Number of Fisher Scoring iterations: 5

broom::glance(Flu2021_themes )

broom::tidy(Flu2021_themes , exponentiate = TRUE)

tbl_regression(Flu2021_themes, label = list(ibd_3 ~ "IBD Type", age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index", act_tob ~ "Active Tobacco Use", IC ~ "Immunocompromised", 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 Type
CD	—	—
UC	1.04	0.92, 1.17	0.5
Unspecified	0.25	0.06, 0.73	0.025
Age	1.01	1.00, 1.01	<0.001
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Race
White	—	—
Black	0.60	0.47, 0.77	<0.001
Asian	1.37	0.99, 1.88	0.055
Native	0.79	0.29, 1.97	0.6
Other	0.65	0.46, 0.92	0.015
Ethnicity
NonHispanic	—	—
Hispanic	1.43	0.96, 2.09	0.073
Preferred Language
English	—	—
Other	0.40	0.20, 0.75	0.007
Marital Status
Married	—	—
Unknown	0.86	0.73, 1.00	0.052
Unmarried	0.69	0.59, 0.80	<0.001
DivorcedSeparated	0.78	0.55, 1.11	0.2
Widow	0.66	0.43, 1.00	0.053
Any Religious Affiliation
Yes	—	—
No	0.89	0.79, 1.00	0.046
Active Tobacco Use
No	—	—
Yes	0.62	0.50, 0.75	<0.001
Charlson Comorbidity Index	1.05	1.03, 1.06	<0.001
Immunocompromised	1.51	1.33, 1.71	<0.001
Population density	1.00	1.00, 1.00	<0.001
Percent Rebuplican	0.98	0.97, 0.98	<0.001
Soceioeconomic Status	0.68	0.48, 0.95	0.026
Household Composition	0.51	0.38, 0.69	<0.001
Minority Status and Language	1.07	0.84, 1.35	0.6
Housing and Transportation	0.83	0.65, 1.06	0.14
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(Flu2021_themes)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7352.332 | 7521.465 |     0.085 | 0.439 | 1.070 |    0.570 |      -Inf |       1.570e-04 | 0.614

performance::check_model(Flu2021_themes, panel = TRUE)

# Margins 
cplot(Flu2021_themes, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME1")

cplot(Flu2021_themes, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME2")

cplot(Flu2021_themes, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME3")

cplot(Flu2021_themes, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME4")

Interaction between Theme 1 and Race

Flu2021_themes.mod <- glm(flu_2021 ~ ibd_3 + age_yrs + gender  + ethnic_3 + 
                      lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC
                      + pop_dens + r_pct + RPL_THEME1*race_5 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              family = "binomial",
              data = vax_clean1)
summary(Flu2021_themes.mod )

Call:
glm(formula = flu_2021 ~ ibd_3 + age_yrs + gender + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + act_tob + max_ch + IC + 
    pop_dens + r_pct + RPL_THEME1 * race_5 + RPL_THEME2 + RPL_THEME3 + 
    RPL_THEME4, family = "binomial", data = vax_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9491  -0.8568  -0.6563   1.1822   2.4726  

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               5.073e-02  2.133e-01   0.238 0.811999    
ibd_3UC                   4.051e-02  6.112e-02   0.663 0.507433    
ibd_3Unspecified         -1.380e+00  6.224e-01  -2.218 0.026581 *  
age_yrs                   6.714e-03  1.991e-03   3.373 0.000744 ***
genderFemale              2.598e-01  5.786e-02   4.490 7.13e-06 ***
ethnic_3Hispanic          3.756e-01  1.995e-01   1.883 0.059749 .  
lang_3Other              -8.832e-01  3.399e-01  -2.598 0.009374 ** 
mstat_5Unknown           -1.570e-01  8.055e-02  -1.949 0.051254 .  
mstat_5Unmarried         -3.701e-01  7.715e-02  -4.798 1.61e-06 ***
mstat_5DivorcedSeparated -2.437e-01  1.787e-01  -1.363 0.172743    
mstat_5Widow             -4.221e-01  2.170e-01  -1.945 0.051775 .  
relig_affilNo            -1.179e-01  5.972e-02  -1.975 0.048253 *  
act_tobYes               -4.840e-01  1.034e-01  -4.679 2.88e-06 ***
max_ch                    4.442e-02  6.229e-03   7.132 9.91e-13 ***
IC                        4.106e-01  6.569e-02   6.250 4.11e-10 ***
pop_dens                 -4.791e-05  1.337e-05  -3.582 0.000340 ***
r_pct                    -2.198e-02  2.369e-03  -9.278  < 2e-16 ***
RPL_THEME1               -3.155e-01  1.811e-01  -1.742 0.081501 .  
race_5Black              -4.209e-01  2.420e-01  -1.739 0.081996 .  
race_5Asian               4.026e-01  2.232e-01   1.804 0.071228 .  
race_5Native              7.990e-01  9.603e-01   0.832 0.405402    
race_5Other              -2.103e-01  2.474e-01  -0.850 0.395350    
RPL_THEME2               -6.785e-01  1.517e-01  -4.473 7.70e-06 ***
RPL_THEME3                5.248e-02  1.200e-01   0.437 0.661880    
RPL_THEME4               -1.892e-01  1.272e-01  -1.487 0.136965    
RPL_THEME1:race_5Black   -2.000e-01  4.207e-01  -0.475 0.634606    
RPL_THEME1:race_5Asian   -4.183e-01  7.289e-01  -0.574 0.566054    
RPL_THEME1:race_5Native  -3.618e+00  3.146e+00  -1.150 0.250098    
RPL_THEME1:race_5Other   -7.803e-01  6.574e-01  -1.187 0.235272    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7855.1  on 6407  degrees of freedom
Residual deviance: 7298.7  on 6379  degrees of freedom
  (8800 observations deleted due to missingness)
AIC: 7356.7

Number of Fisher Scoring iterations: 5

broom::glance(Flu2021_themes.mod )

broom::tidy(Flu2021_themes.mod , exponentiate = TRUE)

tbl_regression(Flu2021_themes.mod, label = list(age_yrs ~ "Age", pop_dens ~ "Population density", gender~ "Gender", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", mstat_5 ~ "Marital Status", r_pct ~ "Percent Rebuplican", race_5 ~ "Race", max_ch ~ "Charlson Comorbidity Index"), exponentiate = TRUE)

Characteristic	OR1	95% CI1	p-value
ibd_3
CD	—	—
UC	1.04	0.92, 1.17	0.5
Unspecified	0.25	0.06, 0.73	0.027
Age	1.01	1.00, 1.01	<0.001
Gender
Male	—	—
Female	1.30	1.16, 1.45	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	1.46	0.98, 2.14	0.060
Preferred Language
English	—	—
Other	0.41	0.20, 0.78	0.009
Marital Status
Married	—	—
Unknown	0.85	0.73, 1.00	0.051
Unmarried	0.69	0.59, 0.80	<0.001
DivorcedSeparated	0.78	0.55, 1.11	0.2
Widow	0.66	0.43, 1.00	0.052
Any Religious Affiliation
Yes	—	—
No	0.89	0.79, 1.00	0.048
act_tob
No	—	—
Yes	0.62	0.50, 0.75	<0.001
Charlson Comorbidity Index	1.05	1.03, 1.06	<0.001
IC	1.51	1.33, 1.72	<0.001
Population density	1.00	1.00, 1.00	<0.001
Percent Rebuplican	0.98	0.97, 0.98	<0.001
Soceioeconomic Status	0.73	0.51, 1.04	0.082
Race
White	—	—
Black	0.66	0.40, 1.05	0.082
Asian	1.50	0.96, 2.32	0.071
Native	2.22	0.38, 19.2	0.4
Other	0.81	0.49, 1.31	0.4
Household Composition	0.51	0.38, 0.68	<0.001
Minority Status and Language	1.05	0.83, 1.33	0.7
Housing and Transportation	0.83	0.64, 1.06	0.14
Soceioeconomic Status * Race
Soceioeconomic Status * Black	0.82	0.36, 1.87	0.6
Soceioeconomic Status * Asian	0.66	0.15, 2.68	0.6
Soceioeconomic Status * Native	0.03	0.00, 3.71	0.3
Soceioeconomic Status * Other	0.46	0.12, 1.58	0.2
1 OR = Odds Ratio, CI = Confidence Interval

# Model performance 
model_performance(Flu2021_themes.mod)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
7356.665 | 7552.859 |     0.085 | 0.439 | 1.070 |    0.569 |      -Inf |       1.566e-04 | 0.614

performance::check_model(Flu2021_themes.mod, panel = TRUE)

# Margins 
cplot(Flu2021_themes.mod, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME1")

cplot(Flu2021_themes.mod, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME2")

cplot(Flu2021_themes.mod, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME3")

cplot(Flu2021_themes.mod, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Flu Vaccine Given THEME4")

Theme 1 as mediator of race -> vaccination

#Step 1: Is there a relationship between X(race) and Y(flu vaccine). Yes - black and other patients were less likely to receive the flu vaccine than Asian and white patients
Flu2019.0.3 <- glm(flu_2021 ~ age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5
                + relig_affil+ max_ch + pop_dens + r_pct, family = binomial, data = vax_clean1)

# Step 2: Is there a relationship between the mediator (SVI) and race? Yes, patients with higher SVI were more likely to be black 
Flu2019_M.3 <- lm(RPL_THEME1 ~ age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5
                + relig_affil+ max_ch + pop_dens + r_pct + flu_2021 + RPL_THEME2 + RPL_THEME3 + RPL_THEME4,
              data = vax_clean1)
summary(Flu2019_M.3 )
summary(Flu2019_M.3 )
broom::glance(Flu2019_M.3 )
tbl_regression(Flu2019_M.3, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"))

# Step 3: Mediation analysis
results3 <- mediate(Flu2019_M.3, Flu2021_themes, treat="race_5", mediator="RPL_THEME1",
                   boot=TRUE, sims=500)
summary(results3)

Theme 1 as mediator without other themes

Flu2019.0.4 <- glm(flu_2021 ~ age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5
                + relig_affil+ max_ch + pop_dens + r_pct + RPL_THEME1, family = binomial, data = vax_clean1)


Flu2019_M.4 <- lm(RPL_THEME1 ~ age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5
                + relig_affil+ max_ch + pop_dens + r_pct + flu_2021,
              data = vax_clean1)
summary(Flu2019_M.4 )
summary(Flu2019_M.4 )
broom::glance(Flu2019_M.4 )
tbl_regression(Flu2019_M.4, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", pop_dens ~ "Population Density", r_pct ~ "Percent Republican"))

# Step 3: Mediation analysis
results4 <- mediate(Flu2019_M.4, Flu2019.0.4, treat="race_5", mediator="RPL_THEME1",
                   boot=TRUE, sims=500)
summary(results4)

```