library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.2     ✔ tidyr     1.3.2
## ✔ readr     2.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
library(MASS)
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
setwd("~/My Stuff/Masters/MPA/PAD 6833/Data")
tab1<-read.delim("PUDF_base1_1q2020_tab.txt")
diabetes_subset <- subset(
  tab1,
  grepl("^E11", PRINC_DIAG_CODE)
)
diabetes_subset$PAT_AGE <- as.numeric(diabetes_subset$PAT_AGE)
## Warning: NAs introduced by coercion
unique(diabetes_subset$PAT_AGE)
##  [1] 12 17 15 16 13 19 18 24 14 11 20  9 21 25 10  6  7 26  4 23  2  5  8  3 NA
## [26] 22
c_diabetes <- diabetes_subset %>%
  filter(
    !is.na(PAT_AGE),
    !PAT_AGE %in% 22:26
  )

#22-26 were filtered out because they were their own data subset including patients that had HIV or were on drugs/alcohol

c_diabetes$PAT_AGE <- as.numeric(c_diabetes$PAT_AGE)
unique(c_diabetes$PAT_AGE)
##  [1] 12 17 15 16 13 19 18 14 11 20  9 21 10  6  7  4  2  5  8  3
prop.table(table(c_diabetes$PAT_AGE))*100
## 
##           2           3           4           5           6           7 
##  0.03213368  0.02410026  0.26510283  0.34543702  0.28920308  0.92384319 
##           8           9          10          11          12          13 
##  1.41388175  2.21722365  4.08901028  5.63142674  8.40295630 10.10604113 
##          14          15          16          17          18          19 
## 12.61246787 13.39170951 12.96593830 10.33097686  7.39877892  5.16548843 
##          20          21 
##  3.09286632  1.30141388

#Because these numbers are associated with an age range (ie 06 corresponds to 18-19) this table is indicating that the biggest proportion of the population is 15 (60-64) 16 which is ages 65-69. 14 is 55-59

min(c_diabetes$PAT_AGE)
## [1] 2
max(c_diabetes$PAT_AGE)
## [1] 21

#Youngest age range: Code 02 = 1-4 years Oldest age range: Code 21 = 90+ years (remember for future analysis that any correlation is not based on getting older in a year but being in an older age group)

c_diabetes <- c_diabetes %>%
filter(!is.na(FIRST_PAYMENT_SRC),
!FIRST_PAYMENT_SRC %in% c("`", "09", "ZZ", "LI", "LM", "AM", "WC")) %>%
mutate(PRIVATE_INS = ifelse(FIRST_PAYMENT_SRC %in% c("MA", "MB", "MC", "VA", "CH", "OF", "TV"),0,1))
table(c_diabetes$PRIVATE_INS)
## 
##    0    1 
## 5397 5083

A value of 1 indicates private insurance coverage, while a value of 0 indicates government-sponsored insurance coverage.

0 = Government insurance

Medicare Part A (MA) Medicare Part B (MB) Medicaid (MC) Veteran Administration (VA) CHAMPUS (CH) Other Federal Program (OF) Title V (TV)

1 = Private/Non-government insurance

PPO (12) POS (13) EPO (14) Indemnity (15) HMO Medicare Risk (16) HMO (HM) Commercial Insurance (CI) Blue Cross/Blue Shield (BL) Disability Insurance (DS) Other Non-federal Programs (11)

table(c_diabetes$PRIVATE_INS)
## 
##    0    1 
## 5397 5083
prop.table(table(c_diabetes$PRIVATE_INS))*100
## 
##        0        1 
## 51.49809 48.50191
c_diabetes <- c_diabetes %>%
filter(!SEX_CODE %in% c("U", "'")) %>%
mutate(SEX_M = ifelse(SEX_CODE == "M", 1, 0))

#Male =1

table(c_diabetes$SEX_M)
## 
##    0    1 
## 4564 5875
prop.table(table(c_diabetes$SEX_M))*100
## 
##        0        1 
## 43.72066 56.27934
c_diabetes <- c_diabetes %>% filter (!RACE %in% c("'", NA)) %>% mutate(WHITE= ifelse(RACE == "4", 1, 0))
table(c_diabetes$WHITE)
## 
##    0    1 
## 3547 6892
prop.table(table(c_diabetes$WHITE))*100
## 
##        0        1 
## 33.97835 66.02165
c_diabetes <- c_diabetes %>% filter (!ETHNICITY %in% c("'", NA)) %>% mutate(HISPANIC= ifelse(ETHNICITY == "1", 1, 0))
table(c_diabetes$HISPANIC)
## 
##    0    1 
## 7279 3160
prop.table(table(c_diabetes$HISPANIC))*100
## 
##       0       1 
## 69.7289 30.2711
is.numeric(c_diabetes$ILLNESS_SEVERITY)
## [1] TRUE
is.numeric(c_diabetes$RISK_MORTALITY)
## [1] TRUE
c_diabetes <- c_diabetes %>%
  filter(ILLNESS_SEVERITY != 0)
illness_model <- lm(ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE + PRIVATE_INS, data=c_diabetes)

summary(illness_model)
## 
## Call:
## lm(formula = ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + 
##     WHITE + PRIVATE_INS, data = c_diabetes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7921 -0.4739 -0.2520  0.5457  1.9961 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.044712   0.042436  48.183  < 2e-16 ***
## PAT_AGE      0.028904   0.002539  11.382  < 2e-16 ***
## SEX_M        0.074716   0.015004   4.980 6.47e-07 ***
## HISPANIC     0.065687   0.016276   4.036 5.48e-05 ***
## WHITE       -0.069512   0.015812  -4.396 1.11e-05 ***
## PRIVATE_INS -0.115845   0.015022  -7.712 1.35e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7594 on 10433 degrees of freedom
## Multiple R-squared:  0.02366,    Adjusted R-squared:  0.0232 
## F-statistic: 50.57 on 5 and 10433 DF,  p-value: < 2.2e-16

#Illness severity indicates the extent of physiologic decompesnation Assignment of a severity of illness score from the All Patient Refined (APR) Diagnosis Related Group (DRG) from the 3Mâ„¢ APR-DRG Grouper.

plot(illness_model,which=1)

plot(illness_model,which=3)

bptest(illness_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  illness_model
## BP = 64.257, df = 5, p-value = 1.598e-12

bptest p-value = 1.819e-10, which is less than 0.05, we can conclude the model is likely heteroscedastic.

plot(illness_model,which=2)

shapiro.test(illness_model$residuals) this cant be done because the model is too big

durbinWatsonTest(illness_model)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.05214514      1.895649       0
##  Alternative hypothesis: rho != 0

#durbinWatson P-value is 0 which means that it likely contains correlation of residuals

illness_model_changed<-rlm(ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE + PRIVATE_INS, data=c_diabetes)
durbinWatsonTest(illness_model_changed)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.05181661      1.896307       0
##  Alternative hypothesis: rho != 0

#To try to mitigate the independent errors violation the model was run as a robust linear model instead of a linear model. This did not work and the p-value is still 0.

qqnorm(resid(illness_model))
qqline(resid(illness_model))

vif(illness_model)
##     PAT_AGE       SEX_M    HISPANIC       WHITE PRIVATE_INS 
##    1.031501    1.002735    1.012184    1.015236    1.020160

#As none of the numbers are >10 there are no variables that are strongly correlated with one another

raintest(illness_model)
## 
##  Rainbow test
## 
## data:  illness_model
## Rain = 1.0081, df1 = 5220, df2 = 5213, p-value = 0.3854

#raintest p-value is 0.3884 which is greater than 0.05 indicating the model is likely linear

plot(illness_model)

mortality_model <- lm(RISK_MORTALITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE + PRIVATE_INS, data=c_diabetes)

summary(mortality_model)
## 
## Call:
## lm(formula = RISK_MORTALITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE + 
##     PRIVATE_INS, data = c_diabetes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8347 -0.4692 -0.1187  0.6203  2.5130 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.425855   0.044756  31.858  < 2e-16 ***
## PAT_AGE      0.064425   0.002678  24.055  < 2e-16 ***
## SEX_M        0.051846   0.015825   3.276  0.00106 ** 
## HISPANIC     0.004120   0.017166   0.240  0.81034    
## WHITE       -0.120382   0.016676  -7.219 5.61e-13 ***
## PRIVATE_INS -0.140617   0.015843  -8.876  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8009 on 10433 degrees of freedom
## Multiple R-squared:  0.06795,    Adjusted R-squared:  0.06751 
## F-statistic: 152.1 on 5 and 10433 DF,  p-value: < 2.2e-16

#Risk Mortality Indicates the likelihood of dying.Assignment of a risk of mortality score from the All Patient Refined (APR) Diagnosis Related Group (DRG) from the 3Mâ„¢ APR-DRG Grouper.

plot(mortality_model,which=1)

plot(mortality_model,which=3)

bptest(mortality_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  mortality_model
## BP = 21.499, df = 5, p-value = 0.0006518
plot(illness_model,which=2)

durbinWatsonTest(illness_model)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.05214514      1.895649       0
##  Alternative hypothesis: rho != 0
raintest(mortality_model)
## 
##  Rainbow test
## 
## data:  mortality_model
## Rain = 0.98044, df1 = 5220, df2 = 5213, p-value = 0.7622
durbinWatsonTest(mortality_model)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.02768931      1.944509   0.008
##  Alternative hypothesis: rho != 0
plot(mortality_model, which=3)

bptest(mortality_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  mortality_model
## BP = 21.499, df = 5, p-value = 0.0006518

#bptest p-value = 0.0006518, which is less than 0.05, we can conclude the model is likely heteroscedastic.

plot(mortality_model, which=2) 

durbinWatsonTest(mortality_model)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.02768931      1.944509   0.008
##  Alternative hypothesis: rho != 0

#durbinWatson P-value is 0.1 which means that it likely contains correlation of residuals

qqnorm(resid(mortality_model))
qqline(resid(mortality_model))

vif(mortality_model)
##     PAT_AGE       SEX_M    HISPANIC       WHITE PRIVATE_INS 
##    1.031501    1.002735    1.012184    1.015236    1.020160

#As none of the numbers are >10 there are no variables that are strongly correlated with one another

raintest(mortality_model)
## 
##  Rainbow test
## 
## data:  mortality_model
## Rain = 0.98044, df1 = 5220, df2 = 5213, p-value = 0.7622

#rainest p-value is 0.7622 which is greater than 0.05 indicating the model is likely linear

plot(mortality_model)

table(c_diabetes$ILLNESS_SEVERITY)
## 
##    1    2    3    4 
## 1015 4611 4032  781
prop.table(table(c_diabetes$ILLNESS_SEVERITY))*100
## 
##         1         2         3         4 
##  9.723154 44.170898 38.624389  7.481560
table(c_diabetes$RISK_MORTALITY)
## 
##    1    2    3    4 
## 1971 4191 3701  576
prop.table(table(c_diabetes$RISK_MORTALITY))*100
## 
##        1        2        3        4 
## 18.88112 40.14752 35.45359  5.51777
ill_ordinal_model <- polr(
  factor(ILLNESS_SEVERITY) ~
    PAT_AGE +
    SEX_M +
    HISPANIC +
    WHITE +
    PRIVATE_INS,
  data = c_diabetes,
  Hess = TRUE
)

summary(ill_ordinal_model)
## Call:
## polr(formula = factor(ILLNESS_SEVERITY) ~ PAT_AGE + SEX_M + HISPANIC + 
##     WHITE + PRIVATE_INS, data = c_diabetes, Hess = TRUE)
## 
## Coefficients:
##                Value Std. Error t value
## PAT_AGE      0.07462   0.006313  11.820
## SEX_M        0.17943   0.037215   4.822
## HISPANIC     0.18049   0.040249   4.484
## WHITE       -0.19218   0.039241  -4.897
## PRIVATE_INS -0.31108   0.037324  -8.335
## 
## Intercepts:
##     Value    Std. Error t value 
## 1|2  -1.2754   0.1064   -11.9844
## 2|3   1.1545   0.1054    10.9580
## 3|4   3.5530   0.1117    31.8156
## 
## Residual Deviance: 23706.36 
## AIC: 23722.36
summary(illness_model)
## 
## Call:
## lm(formula = ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + 
##     WHITE + PRIVATE_INS, data = c_diabetes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7921 -0.4739 -0.2520  0.5457  1.9961 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.044712   0.042436  48.183  < 2e-16 ***
## PAT_AGE      0.028904   0.002539  11.382  < 2e-16 ***
## SEX_M        0.074716   0.015004   4.980 6.47e-07 ***
## HISPANIC     0.065687   0.016276   4.036 5.48e-05 ***
## WHITE       -0.069512   0.015812  -4.396 1.11e-05 ***
## PRIVATE_INS -0.115845   0.015022  -7.712 1.35e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7594 on 10433 degrees of freedom
## Multiple R-squared:  0.02366,    Adjusted R-squared:  0.0232 
## F-statistic: 50.57 on 5 and 10433 DF,  p-value: < 2.2e-16
mort_ordinal_model <- polr(
  factor(RISK_MORTALITY) ~
    PAT_AGE +
    SEX_M +
    HISPANIC +
    WHITE +
    PRIVATE_INS,
  data = c_diabetes,
  Hess = TRUE
)

summary(mort_ordinal_model)
## Call:
## polr(formula = factor(RISK_MORTALITY) ~ PAT_AGE + SEX_M + HISPANIC + 
##     WHITE + PRIVATE_INS, data = c_diabetes, Hess = TRUE)
## 
## Coefficients:
##                Value Std. Error t value
## PAT_AGE      0.15380   0.006426 23.9328
## SEX_M        0.12390   0.036795  3.3673
## HISPANIC     0.01779   0.039868  0.4462
## WHITE       -0.28510   0.038819 -7.3443
## PRIVATE_INS -0.33725   0.036878 -9.1451
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2  0.4988  0.1049     4.7566
## 2|3  2.4316  0.1077    22.5709
## 3|4  4.9836  0.1168    42.6748
## 
## Residual Deviance: 24467.40 
## AIC: 24483.40
summary(mortality_model)
## 
## Call:
## lm(formula = RISK_MORTALITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE + 
##     PRIVATE_INS, data = c_diabetes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8347 -0.4692 -0.1187  0.6203  2.5130 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.425855   0.044756  31.858  < 2e-16 ***
## PAT_AGE      0.064425   0.002678  24.055  < 2e-16 ***
## SEX_M        0.051846   0.015825   3.276  0.00106 ** 
## HISPANIC     0.004120   0.017166   0.240  0.81034    
## WHITE       -0.120382   0.016676  -7.219 5.61e-13 ***
## PRIVATE_INS -0.140617   0.015843  -8.876  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8009 on 10433 degrees of freedom
## Multiple R-squared:  0.06795,    Adjusted R-squared:  0.06751 
## F-statistic: 152.1 on 5 and 10433 DF,  p-value: < 2.2e-16

#Visual inspection of the Q-Q plot suggested departures from normality. Because the dependent variable is measured on a four-level ordinal scale, some deviation from normality was expected. Given the large sample size. Consistent across modeling approaches when using ordinal

chisq.test(
  table(c_diabetes$PRIVATE_INS,
        c_diabetes$ILLNESS_SEVERITY)
)
## 
##  Pearson's Chi-squared test
## 
## data:  table(c_diabetes$PRIVATE_INS, c_diabetes$ILLNESS_SEVERITY)
## X-squared = 117.21, df = 3, p-value < 2.2e-16
prop.table(
  table(c_diabetes$PRIVATE_INS,
        c_diabetes$ILLNESS_SEVERITY),
  1
) * 100
##    
##             1         2         3         4
##   0  8.883107 39.918231 43.151831  8.046831
##   1 10.616845 48.695136 33.807829  6.880190
chisq.test(
  table(c_diabetes$PRIVATE_INS,
        c_diabetes$RISK_MORTALITY)
)
## 
##  Pearson's Chi-squared test
## 
## data:  table(c_diabetes$PRIVATE_INS, c_diabetes$RISK_MORTALITY)
## X-squared = 161.52, df = 3, p-value < 2.2e-16
prop.table(
  table(c_diabetes$PRIVATE_INS,
        c_diabetes$RISK_MORTALITY),
  1
) * 100
##    
##             1         2         3         4
##   0 15.090132 38.914700 40.234157  5.761011
##   1 22.914195 41.459075 30.367734  5.258996