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 patients from Michigan, > 50 years, with at least 6 months of data

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

zoster_MI_50 <- zoster_MI[ which(zoster_MI$age_yrs>50), ]

zoster_MI_50_data <- zoster_MI_50[ which(zoster_MI_50$pyears>0.5), ]

Make variable for fully vaccinated

zoster_MI_50_data %>% 
mutate(full_shingrix = case_when(total_shingrix>= 2 ~ '1',TRUE ~ "0")) -> shingrix_clean1

shingrix_clean1$full_shingrix = as.numeric(shingrix_clean1$full_shingrix)

Baseline characteristics

shingrix_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, full_shingrix) -> shingrix_baseline
shingrix_baseline %>% tbl_summary(
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)")
Characteristic N = 5,8611
IBD Diagnosis
    CD 2,591 (44%)
    UC 3,221 (55%)
    IC 8 (0.1%)
    Unknown 41 (0.7%)
Age 66 (10)
Gender
    Male 2,555 (44%)
    Female 3,306 (56%)
Race
    White 5,263 (90%)
    Black 286 (4.9%)
    Asian or Pacific Islander 100 (1.7%)
    American Indian or Alaska Native 22 (0.4%)
    Other 190 (3.2%)
Ethnicity
    Hispanic 76 (1.3%)
    Non-Hispanic 5,575 (99%)
    (Missing) 210
Preferred Language
    English 5,798 (99%)
    Other 63 (1.1%)
Any Religious Affiliation 3,816 (68%)
    (Missing) 261
Marital Status
    Married 3,403 (70%)
    Single 956 (20%)
    Divorced/Separated 296 (6.1%)
    Widowed 237 (4.8%)
    (Missing) 969
Active Tobacco Use 671 (12%)
    (Missing) 56
Charlson Comorbidity Index 5.6 (5.8)
Insurance Type
    Private Insurance 3,082 (53%)
    Medicaid 538 (9.2%)
    Medicare 2,164 (37%)
    Other Governmental 46 (0.8%)
    Other 23 (0.4%)
    (Missing) 8
Immunocompromised 2,272 (49%)
    (Missing) 1,235
Population Density of Census Tract 1,933 (2,093)
    (Missing) 125
Percentage Republican Voters in Census Tract 46 (17)
    (Missing) 842
Social Vulnerability Index 0.35 (0.25)
    (Missing) 12
SVI by Quartile
    First 2,419 (41%)
    Second 1,777 (30%)
    Third 1,150 (20%)
    Fourth 503 (8.6%)
    (Missing) 12
Socioeconomic Status 0.33 (0.25)
    (Missing) 35
Household Composition 0.39 (0.26)
    (Missing) 12
Minority and Language Status 0.46 (0.29)
    (Missing) 10
Housing and Transportation 0.42 (0.28)
    (Missing) 23
full_shingrix 1,124 (19%)
1 n (%); Mean (SD)

Bivariate Analysis

tbl_zoster_biv <-
  tbl_uvregression(
    shingrix_baseline[c("full_shingrix", "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 = full_shingrix,
    method.args = list(family = binomial),
    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
print(tbl_zoster_biv, method = render)
`...` must be empty.
✖ Problematic argument:
• method = render
Characteristic N OR1 95% CI1 p-value
IBD Diagnosis 5,861
    CD — —
    UC 1.24 1.08, 1.41 0.002
    IC 0.00 >0.9
    Unknown 0.66 0.22, 1.54 0.4
Age 5,861 1.01 1.00, 1.01 0.007
Gender 5,861
    Male — —
    Female 1.05 0.92, 1.20 0.5
Race 5,861
    White — —
    Black 0.55 0.37, 0.78 0.001
    Asian or Pacific Islander 2.25 1.47, 3.39 <0.001
    American Indian or Alaska Native 0.93 0.27, 2.50 0.9
    Other 1.01 0.69, 1.44 >0.9
Ethnicity 5,651
    Hispanic — —
    Non-Hispanic 0.88 0.52, 1.58 0.7
Preferred Language 5,861
    English — —
    Other 0.61 0.27, 1.21 0.2
Marital Status 4,892
    Married — —
    Single 0.83 0.69, 1.00 0.056
    Divorced/Separated 0.74 0.53, 1.01 0.065
    Widowed 0.99 0.71, 1.37 >0.9
Any Religious Affiliation 5,600
    Yes — —
    No 1.03 0.89, 1.18 0.7
Active Tobacco Use 5,805
    No — —
    Yes 0.41 0.31, 0.52 <0.001
Charlson Comorbidity Index 5,861 1.02 1.01, 1.03 <0.001
Immunocompromised 4,626 1.07 0.93, 1.23 0.3
Insurance Type 5,853
    Private Insurance — —
    Medicaid 0.36 0.26, 0.48 <0.001
    Medicare 0.69 0.60, 0.79 <0.001
    Other Governmental 0.42 0.14, 0.96 0.066
    Other 1.21 0.43, 2.91 0.7
Population Density of Census Tract 5,736 1.00 1.00, 1.00 0.006
Social Vulnerability Index 5,849 0.22 0.16, 0.29 <0.001
Socioeconomic Status 5,826 0.16 0.12, 0.21 <0.001
Household Composition 5,849 0.18 0.14, 0.24 <0.001
Minority and Language Status 5,851 1.65 1.31, 2.07 <0.001
Housing and Transportation 5,838 0.47 0.37, 0.60 <0.001
1 OR = Odds Ratio, CI = Confidence Interval 
NULL

Multivariable Models

Zoster + SVI Continuous

zoster_SVI <- glm(full_shingrix ~ 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 = shingrix_clean1)
summary(zoster_SVI )

Call:
glm(formula = full_shingrix ~ 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 = shingrix_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2564  -0.7564  -0.6145  -0.3633   2.3525  

Coefficients:
                                         Estimate Std. Error z value
(Intercept)                            -1.855e+00  5.198e-01  -3.569
ibd_3UC                                 9.908e-02  8.942e-02   1.108
ibd_3IC                                -9.704e+00  1.970e+02  -0.049
ibd_3Unknown                           -7.853e-01  7.692e-01  -1.021
age_yrs                                 2.105e-02  5.392e-03   3.905
genderFemale                            1.094e-01  8.663e-02   1.262
race_5Black                            -5.430e-01  2.518e-01  -2.157
race_5Asian or Pacific Islander         7.187e-01  2.804e-01   2.563
race_5American Indian or Alaska Native  7.056e-01  6.293e-01   1.121
race_5Other                            -7.128e-02  2.937e-01  -0.243
ethnic_3Non-Hispanic                   -6.285e-01  3.732e-01  -1.684
mstat_5Single                           4.507e-02  1.188e-01   0.380
mstat_5Divorced/Separated               6.587e-02  1.945e-01   0.339
mstat_5Widowed                         -9.118e-02  2.026e-01  -0.450
relig_affilNo                           5.444e-02  9.189e-02   0.592
lang_3Other                            -2.488e-01  5.030e-01  -0.495
act_tobYes                             -6.125e-01  1.730e-01  -3.540
max_ch                                  1.792e-02  7.335e-03   2.443
IC                                      2.484e-01  8.892e-02   2.794
insuranceMedicaid                      -5.101e-01  2.008e-01  -2.540
insuranceMedicare                      -5.631e-01  1.112e-01  -5.064
insuranceOther Governmental            -4.594e-01  5.552e-01  -0.827
insuranceOther                          5.371e-01  5.636e-01   0.953
pop_dens                                8.206e-05  1.982e-05   4.141
RPL_THEMES                             -1.341e+00  1.948e-01  -6.887
                                       Pr(>|z|)    
(Intercept)                            0.000359 ***
ibd_3UC                                0.267830    
ibd_3IC                                0.960709    
ibd_3Unknown                           0.307324    
age_yrs                                9.43e-05 ***
genderFemale                           0.206826    
race_5Black                            0.031020 *  
race_5Asian or Pacific Islander        0.010381 *  
race_5American Indian or Alaska Native 0.262151    
race_5Other                            0.808248    
ethnic_3Non-Hispanic                   0.092139 .  
mstat_5Single                          0.704289    
mstat_5Divorced/Separated              0.734896    
mstat_5Widowed                         0.652631    
relig_affilNo                          0.553528    
lang_3Other                            0.620828    
act_tobYes                             0.000400 ***
max_ch                                 0.014577 *  
IC                                     0.005213 ** 
insuranceMedicaid                      0.011091 *  
insuranceMedicare                      4.11e-07 ***
insuranceOther Governmental            0.407982    
insuranceOther                         0.340552    
pop_dens                               3.46e-05 ***
RPL_THEMES                             5.71e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3687.4  on 3487  degrees of freedom
Residual deviance: 3508.0  on 3463  degrees of freedom
  (2373 observations deleted due to missingness)
AIC: 3558

Number of Fisher Scoring iterations: 10
broom::glance(zoster_SVI )
broom::tidy(zoster_SVI , exponentiate = TRUE)
tbl_regression(zoster_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 occurredWarning: collapsing to unique 'x' values
Characteristic OR1 95% CI1 p-value
IBD Diagnosis
    CD — —
    UC 1.10 0.93, 1.32 0.3
    IC 0.00 >0.9
    Unknown 0.46 0.07, 1.68 0.3
Age 1.02 1.01, 1.03 <0.001
Gender
    Male — —
    Female 1.12 0.94, 1.32 0.2
Race
    White — —
    Black 0.58 0.35, 0.93 0.031
    Asian or Pacific Islander 2.05 1.17, 3.54 0.010
    American Indian or Alaska Native 2.03 0.53, 6.67 0.3
    Other 0.93 0.51, 1.62 0.8
Ethnicity
    Hispanic — —
    Non-Hispanic 0.53 0.26, 1.14 0.092
Marital Status
    Married — —
    Single 1.05 0.83, 1.32 0.7
    Divorced/Separated 1.07 0.72, 1.55 0.7
    Widowed 0.91 0.61, 1.35 0.7
Any Religious Affiliation
    Yes — —
    No 1.06 0.88, 1.26 0.6
Preferred Language
    English — —
    Other 0.78 0.27, 1.98 0.6
Active Tobacco Use
    No — —
    Yes 0.54 0.38, 0.75 <0.001
Charlson Comorbidity Index 1.02 1.00, 1.03 0.015
Immunocompromised 1.28 1.08, 1.53 0.005
Insurance Type
    Private Insurance — —
    Medicaid 0.60 0.40, 0.88 0.011
    Medicare 0.57 0.46, 0.71 <0.001
    Other Governmental 0.63 0.18, 1.70 0.4
    Other 1.71 0.52, 4.98 0.3
Population Density of Census Tract 1.00 1.00, 1.00 <0.001
Social Vulnerability Index 0.26 0.18, 0.38 <0.001
1 OR = Odds Ratio, CI = Confidence Interval 
# Hosmer-Lemeshow Goodness-of-Fit Test
hltest(zoster_SVI)

   The Hosmer-Lemeshow goodness-of-fit test

         Statistic =  12.48954 
degrees of freedom =  8 
           p-value =  0.13066 
# C-Statistic/AUROC 
Cstat(zoster_SVI)
[1] 0.6489151
# Model performance 
model_performance(zoster_SVI)
# Indices of model performance

AIC      |     AICc |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
-----------------------------------------------------------------------------------------------------------
3558.026 | 3558.402 | 3711.953 |     0.049 | 0.405 | 1.006 |    0.503 |  -197.690 |       2.867e-04 | 0.672
performance::check_model(zoster_SVI)
Variable `Component` is not in your data frame :/

# Margins 

cplot(zoster_SVI, "RPL_THEMES", what = "prediction", main = "Predicted Likelihood of Zoster Vaccine Given SVI")

Full Shingrix + SVI by Quartile

zoster_SVI4 <- glm(full_shingrix ~ 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 = shingrix_clean1)
summary(zoster_SVI4 )

Call:
glm(formula = full_shingrix ~ 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 = shingrix_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2686  -0.7603  -0.6210  -0.3633   2.3721  

Coefficients:
                                         Estimate Std. Error z value
(Intercept)                            -2.046e+00  5.183e-01  -3.949
ibd_3UC                                 1.037e-01  8.930e-02   1.162
ibd_3IC                                -9.881e+00  1.970e+02  -0.050
ibd_3Unknown                           -7.552e-01  7.682e-01  -0.983
age_yrs                                 2.126e-02  5.388e-03   3.945
genderFemale                            1.091e-01  8.654e-02   1.260
race_5Black                            -5.825e-01  2.522e-01  -2.309
race_5Asian or Pacific Islander         7.557e-01  2.801e-01   2.698
race_5American Indian or Alaska Native  6.401e-01  6.284e-01   1.019
race_5Other                            -8.421e-02  2.938e-01  -0.287
ethnic_3Non-Hispanic                   -6.362e-01  3.727e-01  -1.707
mstat_5Single                           3.379e-02  1.185e-01   0.285
mstat_5Divorced/Separated               4.716e-02  1.943e-01   0.243
mstat_5Widowed                         -9.487e-02  2.024e-01  -0.469
relig_affilNo                           5.748e-02  9.176e-02   0.626
lang_3Other                            -2.532e-01  5.022e-01  -0.504
act_tobYes                             -6.222e-01  1.729e-01  -3.598
max_ch                                  1.779e-02  7.324e-03   2.429
IC                                      2.456e-01  8.885e-02   2.765
insuranceMedicaid                      -5.229e-01  2.007e-01  -2.605
insuranceMedicare                      -5.784e-01  1.109e-01  -5.218
insuranceOther Governmental            -4.656e-01  5.559e-01  -0.838
insuranceOther                          5.188e-01  5.625e-01   0.922
pop_dens                                7.870e-05  1.977e-05   3.981
RPL_4Second                            -1.894e-01  9.622e-02  -1.969
RPL_4Third                             -6.672e-01  1.278e-01  -5.220
RPL_4Fourth                            -8.868e-01  2.164e-01  -4.098
                                       Pr(>|z|)    
(Intercept)                            7.86e-05 ***
ibd_3UC                                0.245351    
ibd_3IC                                0.959989    
ibd_3Unknown                           0.325531    
age_yrs                                7.97e-05 ***
genderFemale                           0.207515    
race_5Black                            0.020919 *  
race_5Asian or Pacific Islander        0.006983 ** 
race_5American Indian or Alaska Native 0.308352    
race_5Other                            0.774379    
ethnic_3Non-Hispanic                   0.087846 .  
mstat_5Single                          0.775484    
mstat_5Divorced/Separated              0.808243    
mstat_5Widowed                         0.639313    
relig_affilNo                          0.531020    
lang_3Other                            0.614158    
act_tobYes                             0.000321 ***
max_ch                                 0.015134 *  
IC                                     0.005699 ** 
insuranceMedicaid                      0.009191 ** 
insuranceMedicare                      1.81e-07 ***
insuranceOther Governmental            0.402294    
insuranceOther                         0.356338    
pop_dens                               6.87e-05 ***
RPL_4Second                            0.048998 *  
RPL_4Third                             1.79e-07 ***
RPL_4Fourth                            4.16e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3687.4  on 3487  degrees of freedom
Residual deviance: 3517.2  on 3461  degrees of freedom
  (2373 observations deleted due to missingness)
AIC: 3571.2

Number of Fisher Scoring iterations: 10
broom::glance(zoster_SVI4 )
broom::tidy(zoster_SVI4 , exponentiate = TRUE)
tbl_regression(zoster_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 1.11 0.93, 1.32 0.2
    IC 0.00 >0.9
    Unknown 0.47 0.07, 1.72 0.3
Age 1.02 1.01, 1.03 <0.001
Gender
    Male — —
    Female 1.12 0.94, 1.32 0.2
Race
    White — —
    Black 0.56 0.33, 0.90 0.021
    Asian or Pacific Islander 2.13 1.22, 3.67 0.007
    American Indian or Alaska Native 1.90 0.49, 6.23 0.3
    Other 0.92 0.50, 1.60 0.8
Ethnicity
    Hispanic — —
    Non-Hispanic 0.53 0.26, 1.13 0.088
Marital Status
    Married — —
    Single 1.03 0.82, 1.30 0.8
    Divorced/Separated 1.05 0.71, 1.52 0.8
    Widowed 0.91 0.61, 1.34 0.6
Any Religious Affiliation
    Yes — —
    No 1.06 0.88, 1.27 0.5
Preferred Language
    English — —
    Other 0.78 0.27, 1.97 0.6
Active Tobacco Use
    No — —
    Yes 0.54 0.38, 0.75 <0.001
Charlson Comorbidity Index 1.02 1.00, 1.03 0.015
Immunocompromised 1.28 1.07, 1.52 0.006
Insurance Type
    Private Insurance — —
    Medicaid 0.59 0.39, 0.87 0.009
    Medicare 0.56 0.45, 0.70 <0.001
    Other Governmental 0.63 0.18, 1.69 0.4
    Other 1.68 0.51, 4.88 0.4
Population Density of Census Tract 1.00 1.00, 1.00 <0.001
SVI by Quartile
    First — —
    Second 0.83 0.68, 1.00 0.049
    Third 0.51 0.40, 0.66 <0.001
    Fourth 0.41 0.26, 0.62 <0.001
1 OR = Odds Ratio, CI = Confidence Interval 
# Hosmer-Lemeshow Goodness-of-Fit Test
hltest(zoster_SVI4)

   The Hosmer-Lemeshow goodness-of-fit test

         Statistic =  9.36475 
degrees of freedom =  8 
           p-value =  0.31247 
# C-Statistic/AUROC 
Cstat(zoster_SVI4)
[1] 0.6439193
# Model performance 
model_performance(zoster_SVI4)
# Indices of model performance

AIC      |     AICc |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
-----------------------------------------------------------------------------------------------------------
3571.153 | 3571.590 | 3737.394 |     0.046 | 0.406 | 1.008 |    0.504 |  -197.317 |       2.867e-04 | 0.671
performance::check_model(zoster_SVI4)
Variable `Component` is not in your data frame :/

# Margins 
cplot(zoster_SVI4, "RPL_4", what = "prediction", main = "Predicted Likelihood of Zoster Vaccine Given SVI Quartile")

Fully vaccinated + All themes

zoster_SVI_themes <- glm(full_shingrix ~ 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 = shingrix_clean1)
summary(zoster_SVI_themes )

Call:
glm(formula = full_shingrix ~ 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 = shingrix_clean1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2408  -0.7618  -0.5816  -0.3324   2.6161  

Coefficients:
                                         Estimate Std. Error z value
(Intercept)                            -1.623e+00  5.339e-01  -3.039
ibd_3UC                                 8.084e-02  9.055e-02   0.893
ibd_3IC                                -9.886e+00  1.970e+02  -0.050
ibd_3Unknown                           -7.231e-01  7.702e-01  -0.939
age_yrs                                 1.814e-02  5.469e-03   3.316
genderFemale                            8.503e-02  8.788e-02   0.968
race_5Black                            -6.329e-01  2.550e-01  -2.482
race_5Asian or Pacific Islander         4.240e-01  2.850e-01   1.488
race_5American Indian or Alaska Native  6.437e-01  6.426e-01   1.002
race_5Other                            -1.004e-01  2.973e-01  -0.338
ethnic_3Non-Hispanic                   -6.230e-01  3.785e-01  -1.646
mstat_5Single                           1.712e-02  1.202e-01   0.142
mstat_5Divorced/Separated               5.664e-02  1.985e-01   0.285
mstat_5Widowed                         -5.102e-02  2.040e-01  -0.250
relig_affilNo                           4.352e-02  9.339e-02   0.466
lang_3Other                            -2.631e-01  5.073e-01  -0.519
act_tobYes                             -5.747e-01  1.743e-01  -3.298
max_ch                                  1.704e-02  7.460e-03   2.284
IC                                      2.798e-01  9.028e-02   3.099
insuranceMedicaid                      -4.352e-01  2.023e-01  -2.151
insuranceMedicare                      -5.100e-01  1.131e-01  -4.509
insuranceOther Governmental            -3.680e-01  5.596e-01  -0.658
insuranceOther                          6.303e-01  5.717e-01   1.102
pop_dens                                4.226e-05  2.243e-05   1.884
RPL_THEME1                             -1.041e+00  2.765e-01  -3.763
RPL_THEME2                             -1.125e+00  2.402e-01  -4.682
RPL_THEME3                              4.156e-01  1.637e-01   2.539
RPL_THEME4                              2.567e-01  1.884e-01   1.363
                                       Pr(>|z|)    
(Intercept)                            0.002371 ** 
ibd_3UC                                0.372031    
ibd_3IC                                0.959971    
ibd_3Unknown                           0.347819    
age_yrs                                0.000913 ***
genderFemale                           0.333261    
race_5Black                            0.013074 *  
race_5Asian or Pacific Islander        0.136768    
race_5American Indian or Alaska Native 0.316461    
race_5Other                            0.735681    
ethnic_3Non-Hispanic                   0.099794 .  
mstat_5Single                          0.886760    
mstat_5Divorced/Separated              0.775452    
mstat_5Widowed                         0.802563    
relig_affilNo                          0.641188    
lang_3Other                            0.604069    
act_tobYes                             0.000975 ***
max_ch                                 0.022352 *  
IC                                     0.001943 ** 
insuranceMedicaid                      0.031443 *  
insuranceMedicare                      6.51e-06 ***
insuranceOther Governmental            0.510726    
insuranceOther                         0.270275    
pop_dens                               0.059594 .  
RPL_THEME1                             0.000168 ***
RPL_THEME2                             2.84e-06 ***
RPL_THEME3                             0.011102 *  
RPL_THEME4                             0.173024    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3647.2  on 3462  degrees of freedom
Residual deviance: 3418.7  on 3435  degrees of freedom
  (2398 observations deleted due to missingness)
AIC: 3474.7

Number of Fisher Scoring iterations: 10
broom::glance(zoster_SVI_themes )
broom::tidy(zoster_SVI_themes , exponentiate = TRUE)
tbl_regression(zoster_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 1.08 0.91, 1.30 0.4
    IC 0.00 >0.9
    Unknown 0.49 0.07, 1.79 0.3
Age 1.02 1.01, 1.03 <0.001
Gender
    Male — —
    Female 1.09 0.92, 1.29 0.3
Race
    White — —
    Black 0.53 0.31, 0.86 0.013
    Asian or Pacific Islander 1.53 0.87, 2.66 0.14
    American Indian or Alaska Native 1.90 0.48, 6.44 0.3
    Other 0.90 0.49, 1.59 0.7
Ethnicity
    Hispanic — —
    Non-Hispanic 0.54 0.26, 1.15 0.10
Marital Status
    Married — —
    Single 1.02 0.80, 1.28 0.9
    Divorced/Separated 1.06 0.71, 1.55 0.8
    Widowed 0.95 0.63, 1.41 0.8
Any Religious Affiliation
    Yes — —
    No 1.04 0.87, 1.25 0.6
Preferred Language
    English — —
    Other 0.77 0.26, 1.97 0.6
Active Tobacco Use
    No — —
    Yes 0.56 0.40, 0.78 <0.001
Charlson Comorbidity Index 1.02 1.00, 1.03 0.022
Immunocompromised 1.32 1.11, 1.58 0.002
Insurance Type
    Private Insurance — —
    Medicaid 0.65 0.43, 0.95 0.031
    Medicare 0.60 0.48, 0.75 <0.001
    Other Governmental 0.69 0.20, 1.88 0.5
    Other 1.88 0.56, 5.57 0.3
Population Density of Census Tract 1.00 1.00, 1.00 0.060
Socioeconomic Status 0.35 0.20, 0.61 <0.001
Household Composition 0.32 0.20, 0.52 <0.001
Minority and Language Status 1.52 1.10, 2.09 0.011
Housing and Transportation 1.29 0.89, 1.87 0.2
1 OR = Odds Ratio, CI = Confidence Interval 
# Hosmer-Lemeshow Goodness-of-Fit Test
hltest(zoster_SVI_themes)

   The Hosmer-Lemeshow goodness-of-fit test

         Statistic =  12.77206 
degrees of freedom =  9 
           p-value =  0.17319 
# C-Statistic/AUROC 
Cstat(zoster_SVI_themes)
[1] 0.6731749
# Model performance 
model_performance(zoster_SVI_themes)
# Indices of model performance

AIC      |     AICc |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
-----------------------------------------------------------------------------------------------------------
3474.724 | 3475.197 | 3646.921 |     0.064 | 0.401 | 0.998 |    0.494 |  -194.897 |       2.888e-04 | 0.679
performance::check_model(zoster_SVI_themes)
Variable `Component` is not in your data frame :/

# Margins 

cplot(zoster_SVI_themes, "RPL_THEME1", what = "prediction", main = "Predicted Likelihood of Zoster Vaccine Given Theme1")

cplot(zoster_SVI_themes, "RPL_THEME2", what = "prediction", main = "Predicted Likelihood of Zoster Vaccine Given Theme2")

cplot(zoster_SVI_themes, "RPL_THEME3", what = "prediction", main = "Predicted Likelihood of Zoster Vaccine Given Theme3")

cplot(zoster_SVI_themes, "RPL_THEME4", what = "prediction", main = "Predicted Likelihood of Zoster Vaccine Given Theme4")

LS0tCnRpdGxlOiAiWm9zdGVyIE1vZGVscyIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICB0aGVtZXM6IHBhcGVyCiAgICB0b2M6IHllcwogICAgdG9jX2Zsb2F0OiB5ZXMKZWRpdG9yX29wdGlvbnM6CiAgY2h1bmtfb3V0cHV0X3R5cGU6IGlubGluZQpkYXRlOiAiMjAyMy0wMS0wOSIKLS0tCgojIExvYWQgUGFja2FnZXMgCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShjb2RlYm9va3IpCmxpYnJhcnkoc3VtbWFyeXRvb2xzKQpsaWJyYXJ5KGJyb29tKSAKbGlicmFyeShwZXJmb3JtYW5jZSkKbGlicmFyeShndCkKbGlicmFyeShndHN1bW1hcnkpCmxpYnJhcnkoamFuaXRvcikKbGlicmFyeShmb3JjYXRzKQpsaWJyYXJ5KG1hcmdpbnMpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShleHBzcykKbGlicmFyeShnbG10b29sYm94KQpsaWJyYXJ5KERlc2NUb29scykKYGBgCgojIEltcG9ydCBkYXRhIApgYGB7cn0KbG9hZCgifi9EZXNrdG9wL1ItQ29kZS9TRE9IX0FMTC9TRE9IX2RhdGFfMDEuMDguMjAyMy5yZGEiKQpgYGAKCiMgRGF0YSBDbGVhbmluZyB7LnRhYnNldH0KCiMjIEluY2x1ZGUgb25seSBwYXRpZW50cyBmcm9tIE1pY2hpZ2FuLCA+IDUwIHllYXJzLCB3aXRoIGF0IGxlYXN0IDYgbW9udGhzIG9mIGRhdGEgCmBgYHtyfQp6b3N0ZXJfTUkgPC0gU0RPSF9kYXRhXzAxLjA4LjIwMjNbIHdoaWNoKFNET0hfZGF0YV8wMS4wOC4yMDIzJFBBVElFTlRfU1RBVEVfQ0Q9PSdNSScpLCBdCgp6b3N0ZXJfTUlfNTAgPC0gem9zdGVyX01JWyB3aGljaCh6b3N0ZXJfTUkkYWdlX3lycz41MCksIF0KCnpvc3Rlcl9NSV81MF9kYXRhIDwtIHpvc3Rlcl9NSV81MFsgd2hpY2goem9zdGVyX01JXzUwJHB5ZWFycz4wLjUpLCBdCgpgYGAKCiMjIE1ha2UgdmFyaWFibGUgZm9yIGZ1bGx5IHZhY2NpbmF0ZWQgCmBgYHtyfQp6b3N0ZXJfTUlfNTBfZGF0YSAlPiUgCm11dGF0ZShmdWxsX3NoaW5ncml4ID0gY2FzZV93aGVuKHRvdGFsX3NoaW5ncml4Pj0gMiB+ICcxJyxUUlVFIH4gIjAiKSkgLT4gc2hpbmdyaXhfY2xlYW4xCgpzaGluZ3JpeF9jbGVhbjEkZnVsbF9zaGluZ3JpeCA9IGFzLm51bWVyaWMoc2hpbmdyaXhfY2xlYW4xJGZ1bGxfc2hpbmdyaXgpCmBgYAoKIyBCYXNlbGluZSBjaGFyYWN0ZXJpc3RpY3MgCmBgYHtyfQpzaGluZ3JpeF9jbGVhbjEgJT4lIAogIGRwbHlyOjpzZWxlY3QoaWJkXzMsIGFnZV95cnMsIGdlbmRlciwgcmFjZV81LCBldGhuaWNfMywgbGFuZ18zLCByZWxpZ19hZmZpbCwgbXN0YXRfNSwgYWN0X3RvYiwgbWF4X2NoLCBpbnN1cmFuY2UsIElDLCBpbnN1cmFuY2UsIHBvcF9kZW5zLHJfcGN0LCBSUExfVEhFTUVTLCBSUExfNCwgUlBMX1RIRU1FMSwgUlBMX1RIRU1FMiwgUlBMX1RIRU1FMywgUlBMX1RIRU1FNCwgZnVsbF9zaGluZ3JpeCkgLT4gc2hpbmdyaXhfYmFzZWxpbmUKc2hpbmdyaXhfYmFzZWxpbmUgJT4lIHRibF9zdW1tYXJ5KAogICAgICAgIHN0YXRpc3RpYyA9IGxpc3QoYWxsX2NvbnRpbnVvdXMoKSB+ICJ7bWVhbn0gKHtzZH0pIiksCiAgICAgICAgbWlzc2luZ190ZXh0ID0gIihNaXNzaW5nKSIpCmBgYAoKIyBCaXZhcmlhdGUgQW5hbHlzaXMgCmBgYHtyfQp0Ymxfem9zdGVyX2JpdiA8LQogIHRibF91dnJlZ3Jlc3Npb24oCiAgICBzaGluZ3JpeF9iYXNlbGluZVtjKCJmdWxsX3NoaW5ncml4IiwgImliZF8zIiwgImFnZV95cnMiLCAiZ2VuZGVyIiwgInJhY2VfNSIsICJldGhuaWNfMyIsICJsYW5nXzMiLCAibXN0YXRfNSIsICJyZWxpZ19hZmZpbCIsICJhY3RfdG9iIiwgIm1heF9jaCIsICJJQyIsICJpbnN1cmFuY2UiLCAicG9wX2RlbnMiLCAiUlBMX1RIRU1FUyIsICJSUExfVEhFTUUxIiwgIlJQTF9USEVNRTIiLCAiUlBMX1RIRU1FMyIsICJSUExfVEhFTUU0IildLAogICAgbWV0aG9kID0gZ2xtLAogICAgeSA9IGZ1bGxfc2hpbmdyaXgsCiAgICBtZXRob2QuYXJncyA9IGxpc3QoZmFtaWx5ID0gYmlub21pYWwpLAogICAgZXhwb25lbnRpYXRlID0gVFJVRSkKcHJpbnQodGJsX3pvc3Rlcl9iaXYsIG1ldGhvZCA9IHJlbmRlcikKYGBgCgojIE11bHRpdmFyaWFibGUgTW9kZWxzIHsudGFic2V0fQoKIyMgWm9zdGVyICsgU1ZJIENvbnRpbnVvdXMgCmBgYHtyfQp6b3N0ZXJfU1ZJIDwtIGdsbShmdWxsX3NoaW5ncml4IH4gaWJkXzMgKyBhZ2VfeXJzICsgZ2VuZGVyICsgcmFjZV81ICsgCiAgICAgICAgICAgICAgICAgICAgZXRobmljXzMgKyBtc3RhdF81ICsgcmVsaWdfYWZmaWwKICAgICAgICAgICAgICAgICAgICAgICsgbGFuZ18zICsgYWN0X3RvYiArIG1heF9jaCArIElDICsgaW5zdXJhbmNlCiAgICAgICAgICAgICAgICAgICAgICArIHBvcF9kZW5zICsgUlBMX1RIRU1FUywKICAgICAgICAgICAgICBmYW1pbHkgPSAiYmlub21pYWwiLCAKICAgICAgICAgICAgICBkYXRhID0gc2hpbmdyaXhfY2xlYW4xKQpzdW1tYXJ5KHpvc3Rlcl9TVkkgKQpicm9vbTo6Z2xhbmNlKHpvc3Rlcl9TVkkgKQpicm9vbTo6dGlkeSh6b3N0ZXJfU1ZJICwgZXhwb25lbnRpYXRlID0gVFJVRSkKdGJsX3JlZ3Jlc3Npb24oem9zdGVyX1NWSSwgZXhwb25lbnRpYXRlID0gVFJVRSkKCiMgSG9zbWVyLUxlbWVzaG93IEdvb2RuZXNzLW9mLUZpdCBUZXN0CmhsdGVzdCh6b3N0ZXJfU1ZJKQoKIyBDLVN0YXRpc3RpYy9BVVJPQyAKQ3N0YXQoem9zdGVyX1NWSSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKHpvc3Rlcl9TVkkpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbCh6b3N0ZXJfU1ZJKQoKIyBNYXJnaW5zIAoKY3Bsb3Qoem9zdGVyX1NWSSwgIlJQTF9USEVNRVMiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBMaWtlbGlob29kIG9mIFpvc3RlciBWYWNjaW5lIEdpdmVuIFNWSSIpCmBgYAoKIyMgRnVsbCBTaGluZ3JpeCArIFNWSSBieSBRdWFydGlsZSAKYGBge3J9Cnpvc3Rlcl9TVkk0IDwtIGdsbShmdWxsX3NoaW5ncml4IH4gaWJkXzMgKyBhZ2VfeXJzICsgZ2VuZGVyICsgcmFjZV81ICsgCiAgICAgICAgICAgICAgICAgICAgZXRobmljXzMgKyBtc3RhdF81ICsgcmVsaWdfYWZmaWwKICAgICAgICAgICAgICAgICAgICAgICsgbGFuZ18zICsgYWN0X3RvYiArIG1heF9jaCArIElDICsgaW5zdXJhbmNlCiAgICAgICAgICAgICAgICAgICAgICArIHBvcF9kZW5zICsgUlBMXzQsCiAgICAgICAgICAgICAgZmFtaWx5ID0gImJpbm9taWFsIiwgCiAgICAgICAgICAgICAgZGF0YSA9IHNoaW5ncml4X2NsZWFuMSkKc3VtbWFyeSh6b3N0ZXJfU1ZJNCApCmJyb29tOjpnbGFuY2Uoem9zdGVyX1NWSTQgKQpicm9vbTo6dGlkeSh6b3N0ZXJfU1ZJNCAsIGV4cG9uZW50aWF0ZSA9IFRSVUUpCnRibF9yZWdyZXNzaW9uKHpvc3Rlcl9TVkk0LCBleHBvbmVudGlhdGUgPSBUUlVFKQoKIyBIb3NtZXItTGVtZXNob3cgR29vZG5lc3Mtb2YtRml0IFRlc3QKaGx0ZXN0KHpvc3Rlcl9TVkk0KQoKIyBDLVN0YXRpc3RpYy9BVVJPQyAKQ3N0YXQoem9zdGVyX1NWSTQpCgojIE1vZGVsIHBlcmZvcm1hbmNlIAptb2RlbF9wZXJmb3JtYW5jZSh6b3N0ZXJfU1ZJNCkKcGVyZm9ybWFuY2U6OmNoZWNrX21vZGVsKHpvc3Rlcl9TVkk0KQoKIyBNYXJnaW5zIApjcGxvdCh6b3N0ZXJfU1ZJNCwgIlJQTF80Iiwgd2hhdCA9ICJwcmVkaWN0aW9uIiwgbWFpbiA9ICJQcmVkaWN0ZWQgTGlrZWxpaG9vZCBvZiBab3N0ZXIgVmFjY2luZSBHaXZlbiBTVkkgUXVhcnRpbGUiKQpgYGAKCiMjIEZ1bGx5IHZhY2NpbmF0ZWQgKyBBbGwgdGhlbWVzIApgYGB7cn0Kem9zdGVyX1NWSV90aGVtZXMgPC0gZ2xtKGZ1bGxfc2hpbmdyaXggfiBpYmRfMyArIGFnZV95cnMgKyBnZW5kZXIgKyByYWNlXzUgKyAKICAgICAgICAgICAgICAgICAgICBldGhuaWNfMyArIG1zdGF0XzUgKyByZWxpZ19hZmZpbAogICAgICAgICAgICAgICAgICAgICAgKyBsYW5nXzMgKyBhY3RfdG9iICsgbWF4X2NoICsgSUMgKyBpbnN1cmFuY2UKICAgICAgICAgICAgICAgICAgICAgICsgcG9wX2RlbnMgKyBSUExfVEhFTUUxICsgUlBMX1RIRU1FMiArCiAgICAgICAgICAgICAgICAgICAgUlBMX1RIRU1FMyArIFJQTF9USEVNRTQsCiAgICAgICAgICAgICAgZmFtaWx5ID0gImJpbm9taWFsIiwgCiAgICAgICAgICAgICAgZGF0YSA9IHNoaW5ncml4X2NsZWFuMSkKc3VtbWFyeSh6b3N0ZXJfU1ZJX3RoZW1lcyApCmJyb29tOjpnbGFuY2Uoem9zdGVyX1NWSV90aGVtZXMgKQpicm9vbTo6dGlkeSh6b3N0ZXJfU1ZJX3RoZW1lcyAsIGV4cG9uZW50aWF0ZSA9IFRSVUUpCnRibF9yZWdyZXNzaW9uKHpvc3Rlcl9TVklfdGhlbWVzLCBleHBvbmVudGlhdGUgPSBUUlVFKQoKIyBIb3NtZXItTGVtZXNob3cgR29vZG5lc3Mtb2YtRml0IFRlc3QKaGx0ZXN0KHpvc3Rlcl9TVklfdGhlbWVzKQoKIyBDLVN0YXRpc3RpYy9BVVJPQyAKQ3N0YXQoem9zdGVyX1NWSV90aGVtZXMpCgojIE1vZGVsIHBlcmZvcm1hbmNlIAptb2RlbF9wZXJmb3JtYW5jZSh6b3N0ZXJfU1ZJX3RoZW1lcykKcGVyZm9ybWFuY2U6OmNoZWNrX21vZGVsKHpvc3Rlcl9TVklfdGhlbWVzKQoKIyBNYXJnaW5zIAoKY3Bsb3Qoem9zdGVyX1NWSV90aGVtZXMsICJSUExfVEhFTUUxIiwgd2hhdCA9ICJwcmVkaWN0aW9uIiwgbWFpbiA9ICJQcmVkaWN0ZWQgTGlrZWxpaG9vZCBvZiBab3N0ZXIgVmFjY2luZSBHaXZlbiBUaGVtZTEiKQpjcGxvdCh6b3N0ZXJfU1ZJX3RoZW1lcywgIlJQTF9USEVNRTIiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBMaWtlbGlob29kIG9mIFpvc3RlciBWYWNjaW5lIEdpdmVuIFRoZW1lMiIpCmNwbG90KHpvc3Rlcl9TVklfdGhlbWVzLCAiUlBMX1RIRU1FMyIsIHdoYXQgPSAicHJlZGljdGlvbiIsIG1haW4gPSAiUHJlZGljdGVkIExpa2VsaWhvb2Qgb2YgWm9zdGVyIFZhY2NpbmUgR2l2ZW4gVGhlbWUzIikKY3Bsb3Qoem9zdGVyX1NWSV90aGVtZXMsICJSUExfVEhFTUU0Iiwgd2hhdCA9ICJwcmVkaWN0aW9uIiwgbWFpbiA9ICJQcmVkaWN0ZWQgTGlrZWxpaG9vZCBvZiBab3N0ZXIgVmFjY2luZSBHaXZlbiBUaGVtZTQiKQpgYGAKCgo=