Load Libraries

Functions

##Print
myprint=function(x, lab="") {print(x)%>%kbl(caption=lab)%>%kable_classic(html_font="Cambria")}


## VIF and Effect Sizes


mydef=function(x){ 
  pseudo=round(1-x$deviance/x$null.deviance,3)
  vif=round(max(vif(x)),3)
  VIF=vif
  R2=pseudo
  return(myprint(noquote(rbind(VIF, R2))))
  }

Read CSV, Pre-processed in Python

y2019=read.csv('D:/MI/MI.CSV', encoding="UTF-8")

Train / Test

set.seed(1234)
mys=sample(seq(1, nrow(y2019)),.2*nrow(y2019),replace=TRUE)
mytrain=y2019[-mys,]
mytest=y2019[mys,]
print(c(nrow(mytrain),nrow(mytest)))

## [1] 285856  69835

Apply Survey Weights

mydata=y2019
options(survey.lonely.psu = "adjust")
svy2019 = svydesign(id=~1,strata = ~Stratum,weights = ~Weights,data = mydata)

train2019=svydesign(id=~1,strata = ~Stratum,weights = ~Weights,data = mytrain)
test2019=svydesign(id=~1,strata = ~Stratum,weights = ~Weights,data = mytest)

Unweighted

options(scipen=9999)
library(psych)
myprint(as.data.frame(apply(y2019,2,function(x)round(mean(x),3))))

##                            apply(y2019, 2, function(x) round(mean(x), 3))
## MI                                                                  0.069
## Male                                                                0.441
## Veteran                                                             0.141
## Rent_Home                                                           0.184
## Poor_Health                                                         0.060
## Smoker                                                              0.428
## Chew_Snuff                                                          0.029
## Percent_Drink                                                       0.154
## High_Cholesterol                                                    0.398
## High_BP                                                             0.463
## High_BMI                                                            0.643
## Poor_Health_Percent                                                 0.159
## Poor_Mental_Health_Percent                                          0.114
## Depression                                                          0.183
## Stroke                                                              0.053
## Asthma                                                              0.134
## Skin_Cancer                                                         0.118
## Cancer                                                              0.120
## COPD                                                                0.097
## Kidney                                                              0.045
## Arthritis                                                           0.388
## No_Health_Plan                                                      0.068
## No_Doctor                                                           0.130
## Cost                                                                0.093
## No_Checkup                                                          0.169
## Metropolitan                                                        0.683
## Weights                                                           508.551
## Stratum                                                        297642.443
## PSU                                                        2019004649.714
## Age_65                                                              0.443
## Age_55_to_64                                                        0.241
## Age_45_to_54                                                        0.175
## Black                                                               0.074
## Hispanic                                                            0.073
## Other_Race                                                          0.064
## Previous_Marriage                                                   0.320
## Never_Married                                                       0.095
## Income_LT25K                                                        0.201
## Income_LT75K                                                        0.335
## Income_DKR                                                          0.175
## Pre_High_School                                                     0.073
## High_School                                                         0.261
## Post_High_School                                                    0.273
## Retired_Unable                                                      0.449
## Out_of_Work                                                         0.032
## Other_Not_Working                                                   0.058
## Fair_Health                                                         0.153
## Good_Health                                                         0.324
## Poor_Exercise                                                       0.279
## Percent_Drink_2                                                     0.101
## Diabetic                                                            0.161
## Prediabetic                                                         0.025
## Division_D1                                                         0.117
## Division_D2                                                         0.049
## Division_D3                                                         0.105
## Division_D4                                                         0.171
## Division_D5                                                         0.182
## Division_D6                                                         0.063
## Division_D7                                                         0.069
## Division_D8                                                         0.131
## Division_D9                                                         0.095

Descriptive Charts, Weighted Data

myd=matrix(rep(0,61*4), 61,4)

for (i in 1:length(svy2019$variables)){
myd[i,1]=colnames(svy2019$variables[i])
myd[i,2]=round(svymean(~svy2019$variables[,i], svy2019),3)
myd[i,3]=round(svyvar(~svy2019$variables[,i], svy2019)^.5,3)
myd[i,4]=round(svytotal(~svy2019$variables[,i], svy2019)/1000000,3)
}
colnames(myd)=c('Variable', 'Mean', 'SD','Total in Millions')
myd=noquote(myd)
myd%>%kbl()%>%kable_classic(html_font='Cambria')

Model 1: Demographics (D)

myformula1=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Black+Hispanic+Other_Race+Previous_Marriage+Never_Married+
            Veteran"))

m1 = survey::svyglm(myformula1,design=train2019,
                    family=quasibinomial,maxit = 100)

pm1=plot_model(m1, main="1. Demographics", show.values=TRUE, show.p=TRUE, value.offset=.4)
mydef(m1)

##      [,1]
## VIF 7.158
## R2  0.083

pm1

Model 2: D + Socioeconomics (S)

myformula2=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Previous_Marriage+Never_Married+Veteran+
            Income_LT25K+Income_LT75K+Income_DKR+Rent_Home+
            Retired_Unable+Out_of_Work+Other_Not_Working+
            Pre_High_School+High_School+Post_High_School"))

m2 = survey::svyglm(myformula2,design=train2019,
                    family=quasibinomial,maxit = 100)

mydef(m2)

##      [,1]
## VIF 8.332
## R2  0.117

pm2=plot_model(m2, main="2. Socioeconomics", show.values=TRUE, show.p=TRUE, value.offset=.4)
pm2

Model 3: D + S + Geography (G)

myformula3=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Previous_Marriage+Never_Married+Veteran+
            Income_LT25K+Income_LT75K+Income_DKR+Rent_Home+
            Retired_Unable+Out_of_Work+Other_Not_Working+
            Pre_High_School+High_School+Post_High_School+
            Division_D1+Division_D2+Division_D3+Division_D4+Division_D5+
            Division_D6+Division_D7+Division_D8+Division_D9+Metropolitan"))

m3 = survey::svyglm(myformula3,design=train2019,
                    family=quasibinomial,maxit = 100)
mydef(m3)

##      [,1]
## VIF 8.417
## R2  0.118

pm3=plot_model(m3, main="3. Geography", show.values=TRUE, show.p=TRUE, value.offset=.4)
pm3

Model 4: D + S + G + Behavioral / Risk (B)

myformula4=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Previous_Marriage+Never_Married+Veteran+
            Income_LT25K+Income_LT75K+Income_DKR+Rent_Home+
            Retired_Unable+
            Pre_High_School+High_School+Post_High_School+
            Division_D3+Division_D4+
            Division_D6+Division_D7+Metropolitan+
            Poor_Health+Fair_Health+Good_Health+
            Smoker+Chew_Snuff+Poor_Exercise+Percent_Drink+Percent_Drink_2+
            High_Cholesterol+High_BP+High_BMI+
            Diabetic+Prediabetic+Poor_Health_Percent+
            Poor_Mental_Health_Percent+Depression+
            Stroke+Asthma+Skin_Cancer+Cancer+COPD+Kidney+Arthritis"))

m4 = survey::svyglm(myformula4,design=train2019,
                    family=quasibinomial,maxit = 100)
mydef(m4)

##       [,1]
## VIF 18.044
## R2   0.209

pm4=plot_model(m4, main="4. Behavior", show.values=TRUE, show.p=TRUE, value.offset=.4)
pm4

Model 5: D + S + G + B + Access (A)

myformula5=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Previous_Marriage+Never_Married+Veteran+
            Income_LT25K+Income_LT75K+Income_DKR+Rent_Home+
            Retired_Unable+
            Pre_High_School+High_School+Post_High_School+
            Division_D4+
            Poor_Health+Fair_Health+Good_Health+
            Smoker+Percent_Drink+Percent_Drink_2+
            High_Cholesterol+High_BP+
            Diabetic+Poor_Health_Percent+
            Stroke+COPD+Kidney+Arthritis+
            No_Health_Plan+No_Doctor+Cost+No_Checkup"))

m5 = survey::svyglm(myformula5,design=train2019,
                    family=quasibinomial,maxit = 100)
mydef(m5)

##       [,1]
## VIF 17.605
## R2   0.210

pm5=plot_model(m5, main="5. Access",show.values=TRUE, show.p=TRUE, value.offset=.4)
pm5

Model 6: Significant Variables Only

myformula6=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Previous_Marriage+Never_Married+Veteran+
            Income_LT25K+Income_LT75K+Income_DKR+Rent_Home+
            Retired_Unable+
            Pre_High_School+High_School+Post_High_School+
            Division_D4+
            Poor_Health+Fair_Health+Good_Health+
            Smoker+Percent_Drink+Percent_Drink_2+
            High_Cholesterol+High_BP+
            Diabetic+Poor_Health_Percent+
            Stroke+COPD+Kidney+Arthritis+
            No_Doctor+Cost+No_Checkup"))

m6 = survey::svyglm(myformula6,design=train2019,
                    family=quasibinomial,maxit = 500)
mydef(m6)

##       [,1]
## VIF 17.602
## R2   0.210

pm6=plot_model(m6, main="6. Total, Only Significant",show.values=TRUE, show.p=TRUE,)
pm6

Model 7: All Variables

myformula7=as.formula(paste("MI","~Age_45_to_54+Age_55_to_64+Age_65+Male+
            Black+Hispanic+Other_Race+Previous_Marriage+Never_Married+
            Veteran+
            Income_LT25K+Income_LT75K+Income_DKR+Rent_Home+
            Retired_Unable+Out_of_Work+Other_Not_Working+
            Pre_High_School+High_School+Post_High_School+
            Division_D1+Division_D2+Division_D3+Division_D4+Division_D5+
            Division_D6+Division_D7+Division_D8+Division_D9+Metropolitan+
            Poor_Health+Fair_Health+Good_Health+
            Smoker+Chew_Snuff+Poor_Exercise+Percent_Drink+Percent_Drink_2+
            High_Cholesterol+High_BP+High_BMI+
            Diabetic+Prediabetic+Poor_Health_Percent+
            Poor_Mental_Health_Percent+Depression+
            Stroke+Asthma+Skin_Cancer+Cancer+COPD+Kidney+Arthritis+
            No_Health_Plan+No_Doctor+Cost+No_Checkup"))

m7 = survey::svyglm(myformula7,design=train2019,
                    family=quasibinomial,maxit = 500)
mydef(m7)

##       [,1]
## VIF 18.156
## R2   0.212

pm7=plot_model(m7, main="7. Total, All Variables", show.values=TRUE, show.p=TRUE,)
pm7

Optimal Cut Points for Quasibinomial and Confusion Matrix on Unseen Data

cutpoints=predict(m6, mytest, type='response')
for (i in 1:length(cutpoints)) {
    if (cutpoints[i]>.17){cutpoints[i]=1}else{cutpoints[i]=0}} #.17 m6


print(noquote(c('F1 Score, No MI: ',F1_Score(cutpoints, mytest$MI, positive='0'))))

## [1] F1 Score, No MI:  0.939864247385317

print(noquote(c('Precision, No MI: ',Precision(cutpoints, mytest$MI, positive='0'))))

## [1] Precision, No MI:  0.925941886234917

print(noquote(c('Recall, No MI: ',Recall(cutpoints, mytest$MI, positive='0'))))

## [1] Recall, No MI:    0.954211669918003

print(noquote(c('F1 Score, MI: ',F1_Score(cutpoints, mytest$MI, positive='1'))))

## [1] F1 Score, MI:     0.338744309885768

print(noquote(c('Precision, MI: ',Precision(cutpoints, mytest$MI, positive='1'))))

## [1] Precision, MI:    0.405844824037868

print(noquote(c('Recall, MI: ',Recall(cutpoints, mytest$MI, positive='1'))))

## [1] Recall, MI:       0.290683962264151

print(noquote(c('Accuracy: ',Accuracy(cutpoints, mytest$MI))))

## [1] Accuracy:         0.889754421135534

print(noquote(c('Pseudo-R2: ', 1-m6$devianc/m6$null.deviance)))

## [1] Pseudo-R2:        0.210038168165197

Group Data for Plotting

t1=pm1$data[,c(1,2,5,6)]
t1$Group=rep("1. Demography", nrow(t1))
t2=pm2$data[,c(1,2,5,6)]
t2$Group=rep("2. SES", nrow(t2))
t3=pm3$data[,c(1,2,5,6)]
t3$Group=rep("3. Geog.", nrow(t3))
t4=pm4$data[,c(1,2,5,6)]
t4$Group=rep("4. Behavior", nrow(t4))
t5=pm5$data[,c(1,2,5,6)]
t5$Group=rep("5. Access", nrow(t5))
t6=pm6$data[,c(1,2,5,6)]
t6$Group=rep("6. Significant Only", nrow(t6))
t7=pm7$data[,c(1,2,5,6)]
t7$Group=rep("7. All Variables", nrow(t7))

ttot=rbind(t1,t2,t3,t4,t5,t6,t7)
ttot$Group=as.factor(ttot$Group)

Plot

ggplot(data=ttot,
    aes(x = term,y = estimate, ymin = .5, ymax = 2.0 ))+
    geom_point(aes(col=Group))+
    geom_hline(aes(fill=Group),yintercept=1, linetype=2)+
    xlab('')+ ylab("Odds Ratio (95% Confidence Interval)")+
    geom_errorbar(aes(ymin=conf.low,
                      ymax=conf.high,col=Group),width=0.5,cex=1)+ 
  facet_grid(~Group)+
        theme(plot.title=element_text(size=16,face="bold"),
        axis.text.y=element_text(size=8),
        axis.text.x=element_text(size=8,face="bold", angle=90),
        axis.title=element_text(size=8,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))+
        guides(colour=FALSE)+
        coord_flip()

## Warning: geom_hline(): Ignoring `mapping` because `yintercept` was provided.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

Table of Predictors for Complete Model

a1=model_parameters(m6, df_method='wald')
a1$CI_high=round(exp(a1$CI_high),3)
a1$CI_low=round(exp(a1$CI_low),3)
a1$Coefficient=round(exp(a1$Coefficient),3)
a1$p=round(a1$p,3)
a1$SE=NULL
a1$t=NULL
a1$df_error=NULL
a1$CI=NULL
myprint(a1)

## Parameter           | Log-Odds |       95% CI |      p
## ------------------------------------------------------
## (Intercept)         | 3.00e-03 | [0.00, 0.00] | < .001
## Age_45_to_54        |     1.74 | [1.45, 2.10] | < .001
## Age_55_to_64        |     2.25 | [1.89, 2.68] | < .001
## Age_65              |     3.30 | [2.76, 3.95] | < .001
## Male                |     2.24 | [2.08, 2.41] | < .001
## Previous_Marriage   |     1.13 | [1.05, 1.22] | 0.001 
## Never_Married       |     0.79 | [0.70, 0.90] | < .001
## Veteran             |     1.12 | [1.03, 1.23] | 0.008 
## Income_LT25K        |     1.30 | [1.15, 1.46] | < .001
## Income_LT75K        |     1.20 | [1.08, 1.33] | < .001
## Income_DKR          |     1.14 | [1.01, 1.27] | 0.028 
## Rent_Home           |     1.10 | [1.01, 1.19] | 0.037 
## Retired_Unable      |     1.29 | [1.19, 1.40] | < .001
## Pre_High_School     |     1.19 | [1.06, 1.34] | 0.004 
## High_School         |     1.10 | [1.01, 1.20] | 0.030 
## Post_High_School    |     1.11 | [1.02, 1.21] | 0.020 
## Division_D4         |     1.11 | [1.04, 1.20] | 0.004 
## Poor_Health         |     2.77 | [2.40, 3.18] | < .001
## Fair_Health         |     2.15 | [1.93, 2.40] | < .001
## Good_Health         |     1.53 | [1.40, 1.67] | < .001
## Smoker              |     1.50 | [1.41, 1.61] | < .001
## Percent_Drink       |     0.37 | [0.22, 0.61] | < .001
## Percent_Drink_2     |     2.19 | [1.28, 3.75] | 0.004 
## High_Cholesterol    |     1.65 | [1.54, 1.76] | < .001
## High_BP             |     1.73 | [1.60, 1.87] | < .001
## Diabetic            |     1.31 | [1.22, 1.42] | < .001
## Poor_Health_Percent |     1.12 | [1.01, 1.24] | 0.032 
## Stroke              |     2.92 | [2.66, 3.20] | < .001
## COPD                |     1.59 | [1.47, 1.73] | < .001
## Kidney              |     1.38 | [1.24, 1.53] | < .001
## Arthritis           |     1.16 | [1.09, 1.25] | < .001
## No_Doctor           |     0.84 | [0.74, 0.95] | 0.005 
## Cost                |     1.23 | [1.11, 1.36] | < .001
## No_Checkup          |     0.78 | [0.69, 0.87] | < .001

Forest Plots of Training Set Estimates and Complete Set Estimates

pm8=plot_model(m6, main="1. Training Set, Significant Variables")
pm9=plot_model(m8, main="2. Full Set, Significant Variables")

t8=pm8$data[,c(1,2,5,6)]
t8$Group=rep("1. Training Set, Significant Variables", nrow(t8))
t9=pm9$data[,c(1,2,5,6)]
t9$Group=rep("2. Full Set, Significant Variables", nrow(t9))


ttot2=rbind(t8,t9)
ttot2$Group=as.factor(ttot2$Group)

ggplot(data=ttot2,
    aes(x = term,y = estimate, ymin = .5, ymax = 2.0 ))+
    geom_point(aes(col=Group))+
    geom_hline(aes(fill=Group),yintercept =1, linetype=2)+
    xlab('')+ ylab("Odds Ratio (95% Confidence Interval)")+
    geom_errorbar(aes(ymin=conf.low,
                      ymax=conf.high,col=Group),width=0.5,cex=1)+ 
  facet_grid(~Group)+
        theme(plot.title=element_text(size=16,face="bold"),
        axis.text.y=element_text(size=10),
        axis.text.x=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))+
        guides(colour=FALSE)+
        coord_flip()

## Warning: geom_hline(): Ignoring `mapping` because `yintercept` was provided.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

Coefficient Comparison, Training versus Full

a1$Interval=noquote(paste0(format(a1$Coefficient,nsmall=3),' (',format(a1$CI_low,nsmall=3),',', format(a1$CI_high,nsmall=3), ')'))
a2$Interval=noquote(paste0(format(a2$Coefficient,nsmall=3),' (',format(a2$CI_low,nsmall=3),',', format(a2$CI_high,nsmall=3), ')'))
newdata=cbind(a1$Parameter,a1$Interval,a2$Interval)
colnames(newdata)=c("Variable","Training Data Odds Ratio (95% CI)", "Full Data Odds Ratio (95% CI)")
myprint(newdata)

##       Variable              Training Data Odds Ratio (95% CI)
##  [1,] "(Intercept)"         "0.003 (0.002,0.003)"            
##  [2,] "Age_45_to_54"        "1.744 (1.451,2.095)"            
##  [3,] "Age_55_to_64"        "2.249 (1.886,2.682)"            
##  [4,] "Age_65"              "3.302 (2.760,3.950)"            
##  [5,] "Male"                "2.241 (2.083,2.410)"            
##  [6,] "Previous_Marriage"   "1.134 (1.052,1.223)"            
##  [7,] "Never_Married"       "0.790 (0.697,0.896)"            
##  [8,] "Veteran"             "1.125 (1.032,1.226)"            
##  [9,] "Income_LT25K"        "1.296 (1.149,1.461)"            
## [10,] "Income_LT75K"        "1.199 (1.084,1.327)"            
## [11,] "Income_DKR"          "1.137 (1.014,1.274)"            
## [12,] "Rent_Home"           "1.096 (1.006,1.194)"            
## [13,] "Retired_Unable"      "1.288 (1.186,1.398)"            
## [14,] "Pre_High_School"     "1.188 (1.056,1.336)"            
## [15,] "High_School"         "1.100 (1.009,1.198)"            
## [16,] "Post_High_School"    "1.108 (1.016,1.208)"            
## [17,] "Division_D4"         "1.114 (1.036,1.197)"            
## [18,] "Poor_Health"         "2.765 (2.405,3.179)"            
## [19,] "Fair_Health"         "2.153 (1.934,2.397)"            
## [20,] "Good_Health"         "1.526 (1.395,1.670)"            
## [21,] "Smoker"              "1.504 (1.408,1.607)"            
## [22,] "Percent_Drink"       "0.368 (0.221,0.613)"            
## [23,] "Percent_Drink_2"     "2.188 (1.277,3.746)"            
## [24,] "High_Cholesterol"    "1.648 (1.539,1.764)"            
## [25,] "High_BP"             "1.730 (1.603,1.867)"            
## [26,] "Diabetic"            "1.315 (1.220,1.416)"            
## [27,] "Poor_Health_Percent" "1.117 (1.009,1.235)"            
## [28,] "Stroke"              "2.915 (2.657,3.197)"            
## [29,] "COPD"                "1.591 (1.466,1.726)"            
## [30,] "Kidney"              "1.377 (1.238,1.532)"            
## [31,] "Arthritis"           "1.164 (1.088,1.245)"            
## [32,] "No_Doctor"           "0.837 (0.738,0.949)"            
## [33,] "Cost"                "1.225 (1.107,1.356)"            
## [34,] "No_Checkup"          "0.777 (0.690,0.874)"            
##       Full Data Odds Ratio (95% CI)
##  [1,] "0.003 (0.002,0.003)"        
##  [2,] "1.597 (1.349,1.891)"        
##  [3,] "2.102 (1.791,2.467)"        
##  [4,] "3.066 (2.606,3.607)"        
##  [5,] "2.171 (2.032,2.320)"        
##  [6,] "1.117 (1.043,1.197)"        
##  [7,] "0.796 (0.707,0.897)"        
##  [8,] "1.159 (1.073,1.253)"        
##  [9,] "1.297 (1.163,1.446)"        
## [10,] "1.218 (1.112,1.335)"        
## [11,] "1.155 (1.038,1.285)"        
## [12,] "1.078 (0.996,1.167)"        
## [13,] "1.297 (1.202,1.400)"        
## [14,] "1.232 (1.106,1.372)"        
## [15,] "1.118 (1.032,1.210)"        
## [16,] "1.091 (1.008,1.181)"        
## [17,] "1.127 (1.055,1.204)"        
## [18,] "2.882 (2.542,3.267)"        
## [19,] "2.189 (1.986,2.412)"        
## [20,] "1.548 (1.426,1.681)"        
## [21,] "1.483 (1.396,1.576)"        
## [22,] "0.382 (0.238,0.612)"        
## [23,] "2.138 (1.299,3.517)"        
## [24,] "1.653 (1.553,1.759)"        
## [25,] "1.683 (1.569,1.805)"        
## [26,] "1.331 (1.244,1.425)"        
## [27,] "1.124 (1.021,1.237)"        
## [28,] "2.950 (2.712,3.209)"        
## [29,] "1.558 (1.446,1.677)"        
## [30,] "1.362 (1.235,1.502)"        
## [31,] "1.168 (1.098,1.242)"        
## [32,] "0.853 (0.759,0.959)"        
## [33,] "1.234 (1.125,1.353)"        
## [34,] "0.807 (0.724,0.900)"