Load Packages

library(tidyverse)
library(codebookr)
library(summarytools)
library(broom) 
library(performance)
library(gt)
library(gtsummary)
library(janitor)
library(forcats)
library(margins)
library(ggplot2)
library(expss)
library(glmtoolbox)
library(DescTools)

Import data

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

Data Cleaning

Include only immunocompromised cohort and people from Michigan

pneumo_MI <- SDOH_data_01.08.2023[ which(SDOH_data_01.08.2023$PATIENT_STATE_CD=='MI'), ]

pneumo_MI_IC <- pneumo_MI[ which(pneumo_MI$IC==1), ]

Include only patients with at least 1 year of data

pneumo_MI_IC_year <- pneumo_MI_IC[ which(pneumo_MI_IC$pyears > 1), ]

Make variables for prevnar, pvax, and total pneumonia vaccine

pneumo_MI_IC_year %>% 
mutate(pvax_2 = case_when(pvax>= 1 ~ '1',TRUE ~ "0")) %>% 
  mutate(prevnar_2 = case_when(prevnar>= 1 ~ '1',TRUE ~ "0")) -> pneumo_1

pneumo_1$pvax_2 = as.numeric(pneumo_1$pvax_2)
pneumo_1$prevnar_2 = as.numeric(pneumo_1$prevnar_2)

pneumo_1 %>% 
  mutate(pneumo_count = pvax_2 + prevnar_2) -> pneumo_clean

pneumo_clean %>% 
  mutate(pneumo_2 = case_when(pneumo_count>= 2 ~ '1',TRUE ~ "0")) -> pneumo_clean1

pneumo_clean1$pneumo_2 = as.numeric(pneumo_clean1$pneumo_2)

Baseline characteristics

pneumo_clean1 %>% 
  dplyr::select(ibd_3, age_yrs, gender, race_5, ethnic_3, lang_3, relig_affil, mstat_5, act_tob, max_ch, insurance, IC, insurance, pop_dens,r_pct, RPL_THEMES, RPL_4, RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4, pvax_2, prevnar_2, pneumo_2) -> pneumo_baseline
pneumo_baseline %>% tbl_summary(
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)")
Characteristic N = 5,7751
IBD Diagnosis
    CD 3,630 (63%)
    UC 2,133 (37%)
    IC 0 (0%)
    Unknown 12 (0.2%)
Age 43 (19)
Gender
    Male 2,777 (48%)
    Female 2,998 (52%)
Race
    White 4,989 (86%)
    Black 410 (7.1%)
    Asian or Pacific Islander 130 (2.3%)
    American Indian or Alaska Native 22 (0.4%)
    Other 224 (3.9%)
Ethnicity
    Hispanic 117 (2.1%)
    Non-Hispanic 5,504 (98%)
    (Missing) 154
Preferred Language
    English 5,734 (99%)
    Other 41 (0.7%)
Any Religious Affiliation 3,040 (55%)
    (Missing) 262
Marital Status
    Married 2,063 (45%)
    Single 2,298 (50%)
    Divorced/Separated 150 (3.3%)
    Widowed 83 (1.8%)
    (Missing) 1,181
Active Tobacco Use 709 (12%)
    (Missing) 20
Charlson Comorbidity Index 2.8 (4.5)
Insurance Type
    Private Insurance 3,958 (69%)
    Medicaid 965 (17%)
    Medicare 793 (14%)
    Other Governmental 34 (0.6%)
    Other 24 (0.4%)
    (Missing) 1
Immunocompromised 5,775 (100%)
Population Density of Census Tract 2,203 (2,351)
    (Missing) 123
Percentage Republican Voters in Census Tract 45 (18)
    (Missing) 807
Social Vulnerability Index 0.37 (0.26)
    (Missing) 18
SVI by Quartile
    First 2,322 (40%)
    Second 1,733 (30%)
    Third 1,113 (19%)
    Fourth 589 (10%)
    (Missing) 18
Socioeconomic Status 0.35 (0.25)
    (Missing) 38
Household Composition 0.39 (0.26)
    (Missing) 18
Minority and Language Status 0.48 (0.29)
    (Missing) 16
Housing and Transportation 0.43 (0.29)
    (Missing) 23
pvax_2 2,500 (43%)
prevnar_2 2,801 (49%)
pneumo_2 1,898 (33%)
1 n (%); Mean (SD)

Bivariate Analysis

tbl_pneumo_biv <-
  tbl_uvregression(
    pneumo_clean1[c("pneumo_2", "ibd_3", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "mstat_5", "relig_affil", "act_tob", "max_ch", "IC", "insurance", "pop_dens", "RPL_THEMES", "RPL_THEME1", "RPL_THEME2", "RPL_THEME3", "RPL_THEME4")],
    method = glm,
    y = pneumo_2,
    method.args = list(family = binomial),
    exponentiate = TRUE)
print(tbl_pneumo_biv, method = render)
`...` must be empty.
✖ Problematic argument:
• method = render
Characteristic N OR1 95% CI1 p-value
IBD Diagnosis 5,775
    CD — —
    UC 0.99 0.88, 1.11 0.9
    Unknown 0.18 0.01, 0.95 0.11
Age 5,775 1.01 1.01, 1.01 <0.001
Gender 5,775
    Male — —
    Female 1.00 0.89, 1.11 >0.9
Race 5,775
    White — —
    Black 0.63 0.50, 0.80 <0.001
    Asian or Pacific Islander 1.08 0.74, 1.54 0.7
    American Indian or Alaska Native 1.12 0.45, 2.63 0.8
    Other 0.74 0.54, 0.99 0.044
Ethnicity 5,621
    Hispanic — —
    Non-Hispanic 1.10 0.74, 1.65 0.7
Preferred Language 5,775
    English — —
    Other 0.84 0.41, 1.62 0.6
Marital Status 4,594
    Married — —
    Single 0.66 0.58, 0.75 <0.001
    Divorced/Separated 0.88 0.62, 1.23 0.4
    Widowed 0.49 0.29, 0.81 0.007
Any Religious Affiliation 5,513
    Yes — —
    No 0.84 0.75, 0.94 0.002
Active Tobacco Use 5,755
    No — —
    Yes 0.78 0.65, 0.92 0.004
Charlson Comorbidity Index 5,775 1.04 1.02, 1.05 <0.001
Immunocompromised 5,775
Insurance Type 5,774
    Private Insurance — —
    Medicaid 0.64 0.54, 0.74 <0.001
    Medicare 1.06 0.90, 1.24 0.5
    Other Governmental 0.59 0.25, 1.24 0.2
    Other 0.50 0.17, 1.25 0.2
Population Density of Census Tract 5,652 1.00 1.00, 1.00 0.4
Social Vulnerability Index 5,757 0.62 0.50, 0.77 <0.001
Socioeconomic Status 5,737 0.64 0.51, 0.79 <0.001
Household Composition 5,757 0.58 0.47, 0.71 <0.001
Minority and Language Status 5,759 1.11 0.91, 1.34 0.3
Housing and Transportation 5,752 0.71 0.59, 0.87 <0.001
1 OR = Odds Ratio, CI = Confidence Interval 
NULL

Multivariable Models: Full Pneumonia Vaccine

Complete Vaccination + SVI Continuous

pneumo_SVI <- glm(pneumo_2 ~ ibd_3 + age_yrs + gender + race_5 + 
                    ethnic_3 + mstat_5 + relig_affil
                      + lang_3 + act_tob + max_ch + IC + insurance
                      + pop_dens + RPL_THEMES,
              family = "binomial", 
              data = pneumo_clean1)
summary(pneumo_SVI )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2820  -0.9373  -0.8075   1.3394   2.0080  

Coefficients: (1 not defined because of singularities)
                                         Estimate Std. Error z value
(Intercept)                            -8.041e-02  2.918e-01  -0.276
ibd_3UC                                -1.255e-01  6.920e-02  -1.813
ibd_3Unknown                           -1.294e+01  2.139e+02  -0.061
age_yrs                                 4.778e-03  2.690e-03   1.776
genderFemale                           -1.051e-02  6.674e-02  -0.158
race_5Black                            -3.792e-01  1.454e-01  -2.608
race_5Asian or Pacific Islander         2.235e-01  2.186e-01   1.022
race_5American Indian or Alaska Native  3.576e-01  5.510e-01   0.649
race_5Other                            -2.656e-01  2.143e-01  -1.239
ethnic_3Non-Hispanic                   -4.803e-01  2.486e-01  -1.932
mstat_5Single                          -2.672e-01  8.602e-02  -3.107
mstat_5Divorced/Separated              -7.455e-02  1.882e-01  -0.396
mstat_5Widowed                         -7.898e-01  2.772e-01  -2.849
relig_affilNo                          -9.169e-02  6.831e-02  -1.342
lang_3Other                             1.341e-01  4.692e-01   0.286
act_tobYes                             -1.342e-01  1.075e-01  -1.248
max_ch                                  2.197e-02  8.015e-03   2.741
IC                                             NA         NA      NA
insuranceMedicaid                      -1.912e-01  1.027e-01  -1.862
insuranceMedicare                      -1.552e-01  1.133e-01  -1.370
insuranceOther Governmental            -8.643e-01  5.590e-01  -1.546
insuranceOther                         -5.881e-01  5.741e-01  -1.024
pop_dens                                4.061e-05  1.453e-05   2.794
RPL_THEMES                             -3.858e-01  1.401e-01  -2.754
                                       Pr(>|z|)   
(Intercept)                             0.78285   
ibd_3UC                                 0.06978 . 
ibd_3Unknown                            0.95175   
age_yrs                                 0.07567 . 
genderFemale                            0.87485   
race_5Black                             0.00910 **
race_5Asian or Pacific Islander         0.30659   
race_5American Indian or Alaska Native  0.51630   
race_5Other                             0.21528   
ethnic_3Non-Hispanic                    0.05340 . 
mstat_5Single                           0.00189 **
mstat_5Divorced/Separated               0.69203   
mstat_5Widowed                          0.00439 **
relig_affilNo                           0.17953   
lang_3Other                             0.77498   
act_tobYes                              0.21200   
max_ch                                  0.00613 **
IC                                           NA   
insuranceMedicaid                       0.06259 . 
insuranceMedicare                       0.17081   
insuranceOther Governmental             0.12205   
insuranceOther                          0.30566   
pop_dens                                0.00520 **
RPL_THEMES                              0.00589 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5389.6  on 4197  degrees of freedom
Residual deviance: 5277.4  on 4175  degrees of freedom
  (1577 observations deleted due to missingness)
AIC: 5323.4

Number of Fisher Scoring iterations: 12
broom::glance(pneumo_SVI )
broom::tidy(pneumo_SVI , exponentiate = TRUE)
tbl_regression(pneumo_SVI, exponentiate = TRUE)
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Characteristic OR1 95% CI1 p-value
IBD Diagnosis
    CD — —
    UC 0.88 0.77, 1.01 0.070
    Unknown 0.00 >0.9
Age 1.00 1.00, 1.01 0.076
Gender
    Male — —
    Female 0.99 0.87, 1.13 0.9
Race
    White — —
    Black 0.68 0.51, 0.91 0.009
    Asian or Pacific Islander 1.25 0.81, 1.91 0.3
    American Indian or Alaska Native 1.43 0.46, 4.20 0.5
    Other 0.77 0.50, 1.16 0.2
Ethnicity
    Hispanic — —
    Non-Hispanic 0.62 0.38, 1.01 0.053
Marital Status
    Married — —
    Single 0.77 0.65, 0.91 0.002
    Divorced/Separated 0.93 0.64, 1.34 0.7
    Widowed 0.45 0.26, 0.77 0.004
Any Religious Affiliation
    Yes — —
    No 0.91 0.80, 1.04 0.2
Preferred Language
    English — —
    Other 1.14 0.44, 2.82 0.8
Active Tobacco Use
    No — —
    Yes 0.87 0.71, 1.08 0.2
Charlson Comorbidity Index 1.02 1.01, 1.04 0.006
Immunocompromised
Insurance Type
    Private Insurance — —
    Medicaid 0.83 0.67, 1.01 0.063
    Medicare 0.86 0.69, 1.07 0.2
    Other Governmental 0.42 0.12, 1.15 0.12
    Other 0.56 0.16, 1.57 0.3
Population Density of Census Tract 1.00 1.00, 1.00 0.005
Social Vulnerability Index 0.68 0.52, 0.89 0.006
1 OR = Odds Ratio, CI = Confidence Interval 
# Hosmer-Lemeshow Goodness-of-Fit Test
hltest(pneumo_SVI)

   The Hosmer-Lemeshow goodness-of-fit test

         Statistic =  5.35251 
degrees of freedom =  8 
           p-value =  0.71932 
# C-Statistic/AUROC 
Cstat(pneumo_SVI)
[1] 0.5964078
# Model performance 
model_performance(pneumo_SVI)
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading
# Indices of model performance

AIC      |     AICc |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
-----------------------------------------------------------------------------------------------------------
5323.446 | 5323.711 | 5469.320 |     0.026 | 0.468 | 1.124 |    0.629 |      -Inf |       5.835e-04 | 0.562
performance::check_model(pneumo_SVI)
Variable `Component` is not in your data frame :/

# Margins 

cplot(pneumo_SVI, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Pneumococcal Vaccine Given SVI")
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading

Complete Vaccination + SVI Quartiles

pneumo_SVI4 <- glm(pneumo_2 ~ ibd_3 + age_yrs + gender + race_5 + 
                    ethnic_3 + mstat_5 + relig_affil
                      + lang_3 + act_tob + max_ch + IC + insurance
                      + pop_dens + RPL_4,
              family = "binomial", 
              data = pneumo_clean1)
summary(pneumo_SVI4 )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3107  -0.9349  -0.8034   1.3392   1.9840  

Coefficients: (1 not defined because of singularities)
                                         Estimate Std. Error z value
(Intercept)                            -8.492e-02  2.905e-01  -0.292
ibd_3UC                                -1.244e-01  6.925e-02  -1.797
ibd_3Unknown                           -1.294e+01  2.134e+02  -0.061
age_yrs                                 4.749e-03  2.691e-03   1.765
genderFemale                           -1.081e-02  6.677e-02  -0.162
race_5Black                            -3.859e-01  1.459e-01  -2.644
race_5Asian or Pacific Islander         2.294e-01  2.188e-01   1.048
race_5American Indian or Alaska Native  3.475e-01  5.527e-01   0.629
race_5Other                            -2.716e-01  2.143e-01  -1.267
ethnic_3Non-Hispanic                   -4.795e-01  2.488e-01  -1.928
mstat_5Single                          -2.696e-01  8.608e-02  -3.132
mstat_5Divorced/Separated              -7.485e-02  1.883e-01  -0.398
mstat_5Widowed                         -7.840e-01  2.774e-01  -2.826
relig_affilNo                          -9.025e-02  6.835e-02  -1.320
lang_3Other                             1.186e-01  4.691e-01   0.253
act_tobYes                             -1.322e-01  1.076e-01  -1.229
max_ch                                  2.224e-02  8.019e-03   2.774
IC                                             NA         NA      NA
insuranceMedicaid                      -1.909e-01  1.027e-01  -1.859
insuranceMedicare                      -1.559e-01  1.132e-01  -1.378
insuranceOther Governmental            -8.529e-01  5.590e-01  -1.526
insuranceOther                         -5.661e-01  5.744e-01  -0.985
pop_dens                                4.000e-05  1.453e-05   2.753
RPL_4Second                            -2.016e-01  7.883e-02  -2.557
RPL_4Third                             -2.309e-01  9.462e-02  -2.440
RPL_4Fourth                            -3.001e-01  1.299e-01  -2.310
                                       Pr(>|z|)   
(Intercept)                             0.77003   
ibd_3UC                                 0.07232 . 
ibd_3Unknown                            0.95163   
age_yrs                                 0.07762 . 
genderFemale                            0.87144   
race_5Black                             0.00819 **
race_5Asian or Pacific Islander         0.29443   
race_5American Indian or Alaska Native  0.52955   
race_5Other                             0.20513   
ethnic_3Non-Hispanic                    0.05390 . 
mstat_5Single                           0.00173 **
mstat_5Divorced/Separated               0.69095   
mstat_5Widowed                          0.00471 **
relig_affilNo                           0.18671   
lang_3Other                             0.80047   
act_tobYes                              0.21912   
max_ch                                  0.00554 **
IC                                           NA   
insuranceMedicaid                       0.06302 . 
insuranceMedicare                       0.16832   
insuranceOther Governmental             0.12705   
insuranceOther                          0.32440   
pop_dens                                0.00591 **
RPL_4Second                             0.01057 * 
RPL_4Third                              0.01467 * 
RPL_4Fourth                             0.02088 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5389.6  on 4197  degrees of freedom
Residual deviance: 5273.9  on 4173  degrees of freedom
  (1577 observations deleted due to missingness)
AIC: 5323.9

Number of Fisher Scoring iterations: 12
broom::glance(pneumo_SVI4 )
broom::tidy(pneumo_SVI4 , exponentiate = TRUE)
tbl_regression(pneumo_SVI4, exponentiate = TRUE)
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Characteristic OR1 95% CI1 p-value
IBD Diagnosis
    CD — —
    UC 0.88 0.77, 1.01 0.072
    Unknown 0.00 >0.9
Age 1.00 1.00, 1.01 0.078
Gender
    Male — —
    Female 0.99 0.87, 1.13 0.9
Race
    White — —
    Black 0.68 0.51, 0.90 0.008
    Asian or Pacific Islander 1.26 0.81, 1.92 0.3
    American Indian or Alaska Native 1.42 0.46, 4.17 0.5
    Other 0.76 0.50, 1.15 0.2
Ethnicity
    Hispanic — —
    Non-Hispanic 0.62 0.38, 1.01 0.054
Marital Status
    Married — —
    Single 0.76 0.64, 0.90 0.002
    Divorced/Separated 0.93 0.64, 1.34 0.7
    Widowed 0.46 0.26, 0.77 0.005
Any Religious Affiliation
    Yes — —
    No 0.91 0.80, 1.04 0.2
Preferred Language
    English — —
    Other 1.13 0.43, 2.78 0.8
Active Tobacco Use
    No — —
    Yes 0.88 0.71, 1.08 0.2
Charlson Comorbidity Index 1.02 1.01, 1.04 0.006
Immunocompromised
Insurance Type
    Private Insurance — —
    Medicaid 0.83 0.67, 1.01 0.063
    Medicare 0.86 0.68, 1.07 0.2
    Other Governmental 0.43 0.12, 1.16 0.13
    Other 0.57 0.16, 1.61 0.3
Population Density of Census Tract 1.00 1.00, 1.00 0.006
SVI by Quartile
    First — —
    Second 0.82 0.70, 0.95 0.011
    Third 0.79 0.66, 0.95 0.015
    Fourth 0.74 0.57, 0.95 0.021
1 OR = Odds Ratio, CI = Confidence Interval 
# Hosmer-Lemeshow Goodness-of-Fit Test
hltest(pneumo_SVI4)

   The Hosmer-Lemeshow goodness-of-fit test

         Statistic =  2.42343 
degrees of freedom =  8 
           p-value =  0.96521 
# C-Statistic/AUROC 
Cstat(pneumo_SVI4)
[1] 0.5976039
# Model performance 
model_performance(pneumo_SVI4)
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading
# Indices of model performance

AIC      |     AICc |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
-----------------------------------------------------------------------------------------------------------
5323.870 | 5324.181 | 5482.429 |     0.027 | 0.468 | 1.124 |    0.628 |      -Inf |       5.835e-04 | 0.562
performance::check_model(pneumo_SVI4)
Variable `Component` is not in your data frame :/

# Margins 
cplot(pneumo_SVI4, "RPL_4", what = "prediction", main = "Predicted Likelihood of Pneumococcal Vaccine Given SVI")
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading

Complete Vaccination + All Themes

pneumo_SVI_themes <- glm(pneumo_2 ~ ibd_3 + age_yrs + gender + race_5 + 
                    ethnic_3 + mstat_5 + relig_affil
                      + lang_3 + act_tob + max_ch + IC + insurance
                      + pop_dens + RPL_THEME1 + RPL_THEME2 + 
                     RPL_THEME3 + RPL_THEME4, 
              family = "binomial", 
              data = pneumo_clean1)
summary(pneumo_SVI_themes )

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3177  -0.9351  -0.8003   1.3286   2.0604  

Coefficients: (1 not defined because of singularities)
                                         Estimate Std. Error z value
(Intercept)                            -2.867e-02  3.030e-01  -0.095
ibd_3UC                                -1.186e-01  6.950e-02  -1.706
ibd_3Unknown                           -1.291e+01  2.125e+02  -0.061
age_yrs                                 4.707e-03  2.698e-03   1.744
genderFemale                           -1.389e-02  6.701e-02  -0.207
race_5Black                            -4.262e-01  1.468e-01  -2.904
race_5Asian or Pacific Islander         1.462e-01  2.215e-01   0.660
race_5American Indian or Alaska Native  3.460e-01  5.513e-01   0.628
race_5Other                            -2.774e-01  2.146e-01  -1.293
ethnic_3Non-Hispanic                   -4.827e-01  2.500e-01  -1.931
mstat_5Single                          -2.756e-01  8.640e-02  -3.189
mstat_5Divorced/Separated              -7.033e-02  1.885e-01  -0.373
mstat_5Widowed                         -7.787e-01  2.776e-01  -2.805
relig_affilNo                          -8.785e-02  6.871e-02  -1.279
lang_3Other                             7.105e-02  4.692e-01   0.151
act_tobYes                             -1.311e-01  1.078e-01  -1.216
max_ch                                  2.185e-02  8.063e-03   2.711
IC                                             NA         NA      NA
insuranceMedicaid                      -1.838e-01  1.035e-01  -1.777
insuranceMedicare                      -1.481e-01  1.138e-01  -1.301
insuranceOther Governmental            -8.750e-01  5.609e-01  -1.560
insuranceOther                         -5.788e-01  5.743e-01  -1.008
pop_dens                                2.197e-05  1.622e-05   1.354
RPL_THEME1                              1.224e-01  2.059e-01   0.594
RPL_THEME2                             -4.378e-01  1.787e-01  -2.450
RPL_THEME3                              1.714e-01  1.271e-01   1.348
RPL_THEME4                             -2.288e-01  1.416e-01  -1.615
                                       Pr(>|z|)   
(Intercept)                             0.92461   
ibd_3UC                                 0.08802 . 
ibd_3Unknown                            0.95154   
age_yrs                                 0.08109 . 
genderFemale                            0.83577   
race_5Black                             0.00369 **
race_5Asian or Pacific Islander         0.50929   
race_5American Indian or Alaska Native  0.53023   
race_5Other                             0.19615   
ethnic_3Non-Hispanic                    0.05348 . 
mstat_5Single                           0.00143 **
mstat_5Divorced/Separated               0.70909   
mstat_5Widowed                          0.00503 **
relig_affilNo                           0.20102   
lang_3Other                             0.87962   
act_tobYes                              0.22397   
max_ch                                  0.00672 **
IC                                           NA   
insuranceMedicaid                       0.07564 . 
insuranceMedicare                       0.19329   
insuranceOther Governmental             0.11875   
insuranceOther                          0.31357   
pop_dens                                0.17573   
RPL_THEME1                              0.55231   
RPL_THEME2                              0.01428 * 
RPL_THEME3                              0.17772   
RPL_THEME4                              0.10628   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5359.9  on 4174  degrees of freedom
Residual deviance: 5238.9  on 4149  degrees of freedom
  (1600 observations deleted due to missingness)
AIC: 5290.9

Number of Fisher Scoring iterations: 12
broom::glance(pneumo_SVI_themes )
broom::tidy(pneumo_SVI_themes , exponentiate = TRUE)
tbl_regression(pneumo_SVI_themes, exponentiate = TRUE)
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredWarning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Characteristic OR1 95% CI1 p-value
IBD Diagnosis
    CD — —
    UC 0.89 0.77, 1.02 0.088
    Unknown 0.00 >0.9
Age 1.00 1.00, 1.01 0.081
Gender
    Male — —
    Female 0.99 0.86, 1.12 0.8
Race
    White — —
    Black 0.65 0.49, 0.87 0.004
    Asian or Pacific Islander 1.16 0.75, 1.78 0.5
    American Indian or Alaska Native 1.41 0.46, 4.16 0.5
    Other 0.76 0.49, 1.14 0.2
Ethnicity
    Hispanic — —
    Non-Hispanic 0.62 0.38, 1.01 0.053
Marital Status
    Married — —
    Single 0.76 0.64, 0.90 0.001
    Divorced/Separated 0.93 0.64, 1.34 0.7
    Widowed 0.46 0.26, 0.78 0.005
Any Religious Affiliation
    Yes — —
    No 0.92 0.80, 1.05 0.2
Preferred Language
    English — —
    Other 1.07 0.41, 2.65 0.9
Active Tobacco Use
    No — —
    Yes 0.88 0.71, 1.08 0.2
Charlson Comorbidity Index 1.02 1.01, 1.04 0.007
Immunocompromised
Insurance Type
    Private Insurance — —
    Medicaid 0.83 0.68, 1.02 0.076
    Medicare 0.86 0.69, 1.08 0.2
    Other Governmental 0.42 0.12, 1.14 0.12
    Other 0.56 0.16, 1.59 0.3
Population Density of Census Tract 1.00 1.00, 1.00 0.2
Socioeconomic Status 1.13 0.75, 1.69 0.6
Household Composition 0.65 0.45, 0.92 0.014
Minority and Language Status 1.19 0.93, 1.52 0.2
Housing and Transportation 0.80 0.60, 1.05 0.11
1 OR = Odds Ratio, CI = Confidence Interval 
# Hosmer-Lemeshow Goodness-of-Fit Test
hltest(pneumo_SVI_themes)

   The Hosmer-Lemeshow goodness-of-fit test

         Statistic =  5.74428 
degrees of freedom =  8 
           p-value =  0.67585 
# C-Statistic/AUROC 
Cstat(pneumo_SVI_themes)
[1] 0.6015217
# Model performance 
model_performance(pneumo_SVI_themes)
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading
# Indices of model performance

AIC      |     AICc |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
-----------------------------------------------------------------------------------------------------------
5290.864 | 5291.203 | 5455.623 |     0.028 | 0.467 | 1.124 |    0.627 |      -Inf |       5.867e-04 | 0.563
performance::check_model(pneumo_SVI_themes)
Variable `Component` is not in your data frame :/

# Margins 
cplot(pneumo_SVI_themes, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Pneumococcal Vaccine Given Theme1")
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading

cplot(pneumo_SVI_themes, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Pneumococcal Vaccine Given Theme2")
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading

cplot(pneumo_SVI_themes, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Pneumococcal Vaccine Given Theme3")
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading

cplot(pneumo_SVI_themes, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Pneumococcal Vaccine Given Theme4")
Warning: prediction from a rank-deficient fit may be misleadingWarning: prediction from a rank-deficient fit may be misleading

LS0tCnRpdGxlOiAiUG5ldW1vY29jY2FsIE1vZGVscyIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICB0aGVtZXM6IHBhcGVyCiAgICB0b2M6IHllcwogICAgdG9jX2Zsb2F0OiB5ZXMKZWRpdG9yX29wdGlvbnM6CiAgY2h1bmtfb3V0cHV0X3R5cGU6IGlubGluZQpkYXRlOiAiMjAyMy0wMS0wOSIKLS0tCgojIExvYWQgUGFja2FnZXMgCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShjb2RlYm9va3IpCmxpYnJhcnkoc3VtbWFyeXRvb2xzKQpsaWJyYXJ5KGJyb29tKSAKbGlicmFyeShwZXJmb3JtYW5jZSkKbGlicmFyeShndCkKbGlicmFyeShndHN1bW1hcnkpCmxpYnJhcnkoamFuaXRvcikKbGlicmFyeShmb3JjYXRzKQpsaWJyYXJ5KG1hcmdpbnMpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShleHBzcykKbGlicmFyeShnbG10b29sYm94KQpsaWJyYXJ5KERlc2NUb29scykKYGBgCgojIEltcG9ydCBkYXRhIApgYGB7cn0KbG9hZCgifi9EZXNrdG9wL1ItQ29kZS9TRE9IX0FMTC9TRE9IX2RhdGFfMDEuMDguMjAyMy5yZGEiKQpgYGAKCiMgRGF0YSBDbGVhbmluZyB7LnRhYnNldH0KCiMjIEluY2x1ZGUgb25seSBpbW11bm9jb21wcm9taXNlZCBjb2hvcnQgYW5kIHBlb3BsZSBmcm9tIE1pY2hpZ2FuIApgYGB7cn0KcG5ldW1vX01JIDwtIFNET0hfZGF0YV8wMS4wOC4yMDIzWyB3aGljaChTRE9IX2RhdGFfMDEuMDguMjAyMyRQQVRJRU5UX1NUQVRFX0NEPT0nTUknKSwgXQoKcG5ldW1vX01JX0lDIDwtIHBuZXVtb19NSVsgd2hpY2gocG5ldW1vX01JJElDPT0xKSwgXQpgYGAKCiMjIEluY2x1ZGUgb25seSBwYXRpZW50cyB3aXRoIGF0IGxlYXN0IDEgeWVhciBvZiBkYXRhIApgYGB7cn0KcG5ldW1vX01JX0lDX3llYXIgPC0gcG5ldW1vX01JX0lDWyB3aGljaChwbmV1bW9fTUlfSUMkcHllYXJzID4gMSksIF0KYGBgCgojIyBNYWtlIHZhcmlhYmxlcyBmb3IgcHJldm5hciwgcHZheCwgYW5kIHRvdGFsIHBuZXVtb25pYSB2YWNjaW5lIApgYGB7cn0KcG5ldW1vX01JX0lDX3llYXIgJT4lIAptdXRhdGUocHZheF8yID0gY2FzZV93aGVuKHB2YXg+PSAxIH4gJzEnLFRSVUUgfiAiMCIpKSAlPiUgCiAgbXV0YXRlKHByZXZuYXJfMiA9IGNhc2Vfd2hlbihwcmV2bmFyPj0gMSB+ICcxJyxUUlVFIH4gIjAiKSkgLT4gcG5ldW1vXzEKCnBuZXVtb18xJHB2YXhfMiA9IGFzLm51bWVyaWMocG5ldW1vXzEkcHZheF8yKQpwbmV1bW9fMSRwcmV2bmFyXzIgPSBhcy5udW1lcmljKHBuZXVtb18xJHByZXZuYXJfMikKCnBuZXVtb18xICU+JSAKICBtdXRhdGUocG5ldW1vX2NvdW50ID0gcHZheF8yICsgcHJldm5hcl8yKSAtPiBwbmV1bW9fY2xlYW4KCnBuZXVtb19jbGVhbiAlPiUgCiAgbXV0YXRlKHBuZXVtb18yID0gY2FzZV93aGVuKHBuZXVtb19jb3VudD49IDIgfiAnMScsVFJVRSB+ICIwIikpIC0+IHBuZXVtb19jbGVhbjEKCnBuZXVtb19jbGVhbjEkcG5ldW1vXzIgPSBhcy5udW1lcmljKHBuZXVtb19jbGVhbjEkcG5ldW1vXzIpCmBgYAoKIyBCYXNlbGluZSBjaGFyYWN0ZXJpc3RpY3MgCmBgYHtyfQpwbmV1bW9fY2xlYW4xICU+JSAKICBkcGx5cjo6c2VsZWN0KGliZF8zLCBhZ2VfeXJzLCBnZW5kZXIsIHJhY2VfNSwgZXRobmljXzMsIGxhbmdfMywgcmVsaWdfYWZmaWwsIG1zdGF0XzUsIGFjdF90b2IsIG1heF9jaCwgaW5zdXJhbmNlLCBJQywgaW5zdXJhbmNlLCBwb3BfZGVucyxyX3BjdCwgUlBMX1RIRU1FUywgUlBMXzQsIFJQTF9USEVNRTEsIFJQTF9USEVNRTIsIFJQTF9USEVNRTMsIFJQTF9USEVNRTQsIHB2YXhfMiwgcHJldm5hcl8yLCBwbmV1bW9fMikgLT4gcG5ldW1vX2Jhc2VsaW5lCnBuZXVtb19iYXNlbGluZSAlPiUgdGJsX3N1bW1hcnkoCiAgICAgICAgc3RhdGlzdGljID0gbGlzdChhbGxfY29udGludW91cygpIH4gInttZWFufSAoe3NkfSkiKSwKICAgICAgICBtaXNzaW5nX3RleHQgPSAiKE1pc3NpbmcpIikKYGBgCgojIEJpdmFyaWF0ZSBBbmFseXNpcyAKYGBge3J9CnRibF9wbmV1bW9fYml2IDwtCiAgdGJsX3V2cmVncmVzc2lvbigKICAgIHBuZXVtb19jbGVhbjFbYygicG5ldW1vXzIiLCAiaWJkXzMiLCAiYWdlX3lycyIsICJnZW5kZXIiLCAicmFjZV81IiwgImV0aG5pY18zIiwgImxhbmdfMyIsICJtc3RhdF81IiwgInJlbGlnX2FmZmlsIiwgImFjdF90b2IiLCAibWF4X2NoIiwgIklDIiwgImluc3VyYW5jZSIsICJwb3BfZGVucyIsICJSUExfVEhFTUVTIiwgIlJQTF9USEVNRTEiLCAiUlBMX1RIRU1FMiIsICJSUExfVEhFTUUzIiwgIlJQTF9USEVNRTQiKV0sCiAgICBtZXRob2QgPSBnbG0sCiAgICB5ID0gcG5ldW1vXzIsCiAgICBtZXRob2QuYXJncyA9IGxpc3QoZmFtaWx5ID0gYmlub21pYWwpLAogICAgZXhwb25lbnRpYXRlID0gVFJVRSkKcHJpbnQodGJsX3BuZXVtb19iaXYsIG1ldGhvZCA9IHJlbmRlcikKYGBgCgojIE11bHRpdmFyaWFibGUgTW9kZWxzOiBGdWxsIFBuZXVtb25pYSBWYWNjaW5lIAojIyBDb21wbGV0ZSBWYWNjaW5hdGlvbiArIFNWSSBDb250aW51b3VzIApgYGB7cn0KcG5ldW1vX1NWSSA8LSBnbG0ocG5ldW1vXzIgfiBpYmRfMyArIGFnZV95cnMgKyBnZW5kZXIgKyByYWNlXzUgKyAKICAgICAgICAgICAgICAgICAgICBldGhuaWNfMyArIG1zdGF0XzUgKyByZWxpZ19hZmZpbAogICAgICAgICAgICAgICAgICAgICAgKyBsYW5nXzMgKyBhY3RfdG9iICsgbWF4X2NoICsgSUMgKyBpbnN1cmFuY2UKICAgICAgICAgICAgICAgICAgICAgICsgcG9wX2RlbnMgKyBSUExfVEhFTUVTLAogICAgICAgICAgICAgIGZhbWlseSA9ICJiaW5vbWlhbCIsIAogICAgICAgICAgICAgIGRhdGEgPSBwbmV1bW9fY2xlYW4xKQpzdW1tYXJ5KHBuZXVtb19TVkkgKQpicm9vbTo6Z2xhbmNlKHBuZXVtb19TVkkgKQpicm9vbTo6dGlkeShwbmV1bW9fU1ZJICwgZXhwb25lbnRpYXRlID0gVFJVRSkKdGJsX3JlZ3Jlc3Npb24ocG5ldW1vX1NWSSwgZXhwb25lbnRpYXRlID0gVFJVRSkKCiMgSG9zbWVyLUxlbWVzaG93IEdvb2RuZXNzLW9mLUZpdCBUZXN0CmhsdGVzdChwbmV1bW9fU1ZJKQoKIyBDLVN0YXRpc3RpYy9BVVJPQyAKQ3N0YXQocG5ldW1vX1NWSSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKHBuZXVtb19TVkkpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChwbmV1bW9fU1ZJKQoKIyBNYXJnaW5zIAoKY3Bsb3QocG5ldW1vX1NWSSwgIlJQTF9USEVNRVMiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBMaWtlbGlob29kIG9mIFBuZXVtb2NvY2NhbCBWYWNjaW5lIEdpdmVuIFNWSSIpCmBgYAoKIyBDb21wbGV0ZSBWYWNjaW5hdGlvbiArIFNWSSBRdWFydGlsZXMgCmBgYHtyfQpwbmV1bW9fU1ZJNCA8LSBnbG0ocG5ldW1vXzIgfiBpYmRfMyArIGFnZV95cnMgKyBnZW5kZXIgKyByYWNlXzUgKyAKICAgICAgICAgICAgICAgICAgICBldGhuaWNfMyArIG1zdGF0XzUgKyByZWxpZ19hZmZpbAogICAgICAgICAgICAgICAgICAgICAgKyBsYW5nXzMgKyBhY3RfdG9iICsgbWF4X2NoICsgSUMgKyBpbnN1cmFuY2UKICAgICAgICAgICAgICAgICAgICAgICsgcG9wX2RlbnMgKyBSUExfNCwKICAgICAgICAgICAgICBmYW1pbHkgPSAiYmlub21pYWwiLCAKICAgICAgICAgICAgICBkYXRhID0gcG5ldW1vX2NsZWFuMSkKc3VtbWFyeShwbmV1bW9fU1ZJNCApCmJyb29tOjpnbGFuY2UocG5ldW1vX1NWSTQgKQpicm9vbTo6dGlkeShwbmV1bW9fU1ZJNCAsIGV4cG9uZW50aWF0ZSA9IFRSVUUpCnRibF9yZWdyZXNzaW9uKHBuZXVtb19TVkk0LCBleHBvbmVudGlhdGUgPSBUUlVFKQoKIyBIb3NtZXItTGVtZXNob3cgR29vZG5lc3Mtb2YtRml0IFRlc3QKaGx0ZXN0KHBuZXVtb19TVkk0KQoKIyBDLVN0YXRpc3RpYy9BVVJPQyAKQ3N0YXQocG5ldW1vX1NWSTQpCgojIE1vZGVsIHBlcmZvcm1hbmNlIAptb2RlbF9wZXJmb3JtYW5jZShwbmV1bW9fU1ZJNCkKcGVyZm9ybWFuY2U6OmNoZWNrX21vZGVsKHBuZXVtb19TVkk0KQoKIyBNYXJnaW5zIApjcGxvdChwbmV1bW9fU1ZJNCwgIlJQTF80Iiwgd2hhdCA9ICJwcmVkaWN0aW9uIiwgbWFpbiA9ICJQcmVkaWN0ZWQgTGlrZWxpaG9vZCBvZiBQbmV1bW9jb2NjYWwgVmFjY2luZSBHaXZlbiBTVkkiKQpgYGAKCiMjIENvbXBsZXRlIFZhY2NpbmF0aW9uICsgQWxsIFRoZW1lcyAKYGBge3J9CnBuZXVtb19TVklfdGhlbWVzIDwtIGdsbShwbmV1bW9fMiB+IGliZF8zICsgYWdlX3lycyArIGdlbmRlciArIHJhY2VfNSArIAogICAgICAgICAgICAgICAgICAgIGV0aG5pY18zICsgbXN0YXRfNSArIHJlbGlnX2FmZmlsCiAgICAgICAgICAgICAgICAgICAgICArIGxhbmdfMyArIGFjdF90b2IgKyBtYXhfY2ggKyBJQyArIGluc3VyYW5jZQogICAgICAgICAgICAgICAgICAgICAgKyBwb3BfZGVucyArIFJQTF9USEVNRTEgKyBSUExfVEhFTUUyICsgCiAgICAgICAgICAgICAgICAgICAgIFJQTF9USEVNRTMgKyBSUExfVEhFTUU0LCAKICAgICAgICAgICAgICBmYW1pbHkgPSAiYmlub21pYWwiLCAKICAgICAgICAgICAgICBkYXRhID0gcG5ldW1vX2NsZWFuMSkKc3VtbWFyeShwbmV1bW9fU1ZJX3RoZW1lcyApCmJyb29tOjpnbGFuY2UocG5ldW1vX1NWSV90aGVtZXMgKQpicm9vbTo6dGlkeShwbmV1bW9fU1ZJX3RoZW1lcyAsIGV4cG9uZW50aWF0ZSA9IFRSVUUpCnRibF9yZWdyZXNzaW9uKHBuZXVtb19TVklfdGhlbWVzLCBleHBvbmVudGlhdGUgPSBUUlVFKQoKIyBIb3NtZXItTGVtZXNob3cgR29vZG5lc3Mtb2YtRml0IFRlc3QKaGx0ZXN0KHBuZXVtb19TVklfdGhlbWVzKQoKIyBDLVN0YXRpc3RpYy9BVVJPQyAKQ3N0YXQocG5ldW1vX1NWSV90aGVtZXMpCgojIE1vZGVsIHBlcmZvcm1hbmNlIAptb2RlbF9wZXJmb3JtYW5jZShwbmV1bW9fU1ZJX3RoZW1lcykKcGVyZm9ybWFuY2U6OmNoZWNrX21vZGVsKHBuZXVtb19TVklfdGhlbWVzKQoKIyBNYXJnaW5zIApjcGxvdChwbmV1bW9fU1ZJX3RoZW1lcywgIlJQTF9USEVNRTEiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBMaWtlbGlob29kIG9mIFBuZXVtb2NvY2NhbCBWYWNjaW5lIEdpdmVuIFRoZW1lMSIpCmNwbG90KHBuZXVtb19TVklfdGhlbWVzLCAiUlBMX1RIRU1FMiIsIHdoYXQgPSAicHJlZGljdGlvbiIsIG1haW4gPSAiUHJlZGljdGVkIExpa2VsaWhvb2Qgb2YgUG5ldW1vY29jY2FsIFZhY2NpbmUgR2l2ZW4gVGhlbWUyIikKY3Bsb3QocG5ldW1vX1NWSV90aGVtZXMsICJSUExfVEhFTUUzIiwgd2hhdCA9ICJwcmVkaWN0aW9uIiwgbWFpbiA9ICJQcmVkaWN0ZWQgTGlrZWxpaG9vZCBvZiBQbmV1bW9jb2NjYWwgVmFjY2luZSBHaXZlbiBUaGVtZTMiKQpjcGxvdChwbmV1bW9fU1ZJX3RoZW1lcywgIlJQTF9USEVNRTQiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBMaWtlbGlob29kIG9mIFBuZXVtb2NvY2NhbCBWYWNjaW5lIEdpdmVuIFRoZW1lNCIpCmBgYAoK