setwd("~/My Stuff/Masters/MPA/PAD 6833/Data")
tab1<-read.delim("PUDF_base1_1q2020_tab.txt")
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
diabetes_subset <- subset(
  tab1,
  grepl("^E11", PRINC_DIAG_CODE)
)
c_diabetes<- diabetes_subset |> drop_na(PAT_AGE)
c_diabetes$PAT_AGE <- as.numeric(c_diabetes$PAT_AGE)
## Warning: NAs introduced by coercion
c_diabetes <- c_diabetes %>%
  filter(!FIRST_PAYMENT_SRC %in% c("`", "", "ZZ", NA)) %>%
  mutate(
    INSURED = ifelse(FIRST_PAYMENT_SRC == "09", 0, 1)
  )
c_diabetes <- c_diabetes %>% filter (!SEX_CODE %in% c("U", "'", NA)) %>% mutate(SEX_M = ifelse(SEX_CODE == "M", 0, 1))
c_diabetes <- c_diabetes %>% filter (!RACE %in% c("'","5", NA)) %>% mutate(WHITE= ifelse(RACE == "4", 1, 0))
c_diabetes <- c_diabetes %>% filter (!ETHNICITY %in% c("'", NA)) %>% mutate(HISPANIC= ifelse(ETHNICITY == "1", 1, 0))
illness_model1 <- lm(ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE, data=c_diabetes)
summary(illness_model1)
## 
## Call:
## lm(formula = ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + 
##     WHITE, data = c_diabetes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6376 -0.4634 -0.3124  0.5584  1.8826 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.170059   0.036700  59.129  < 2e-16 ***
## PAT_AGE      0.021776   0.002199   9.904  < 2e-16 ***
## SEX_M       -0.086135   0.015779  -5.459 4.91e-08 ***
## HISPANIC     0.084696   0.018891   4.483 7.43e-06 ***
## WHITE       -0.075423   0.019157  -3.937 8.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7585 on 9632 degrees of freedom
## Multiple R-squared:  0.013,  Adjusted R-squared:  0.01259 
## F-statistic: 31.71 on 4 and 9632 DF,  p-value: < 2.2e-16
plot(illness_model1,which=1)

raintest(illness_model1)
## 
##  Rainbow test
## 
## data:  illness_model1
## Rain = 1.0104, df1 = 4819, df2 = 4813, p-value = 0.3595

#Rain test P-value = 0.3595 which is greater than .05 which means it is likely linear.

durbinWatsonTest(illness_model1)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.03919167      1.921537       0
##  Alternative hypothesis: rho != 0

#the p-value is less than 0.05 meaning that i likely contains correlation of residuals

plot(illness_model1, which=3)

bptest(illness_model1)
## 
##  studentized Breusch-Pagan test
## 
## data:  illness_model1
## BP = 25.941, df = 4, p-value = 3.253e-05

#The p-value is below .05 so it is likely heteroscedastic

plot(illness_model1, which=2)

vif(illness_model1)
##  PAT_AGE    SEX_M HISPANIC    WHITE 
## 1.046545 1.039599 1.115261 1.116768

#all are below 5 and 10 so there is probably no multicolinearity

#The only failure was that of the Durbin-watson test so there are correlations between residuals.

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
illness_model_robust <- rlm(ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + WHITE, data=c_diabetes)
summary(illness_model_robust)
## 
## Call: rlm(formula = ILLNESS_SEVERITY ~ PAT_AGE + SEX_M + HISPANIC + 
##     WHITE, data = c_diabetes)
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6471 -0.4600 -0.3057  0.5625  1.8960 
## 
## Coefficients:
##             Value   Std. Error t value
## (Intercept)  2.1680  0.0359    60.3176
## PAT_AGE      0.0224  0.0022    10.4217
## SEX_M       -0.0867  0.0155    -5.6083
## HISPANIC     0.0942  0.0185     5.0921
## WHITE       -0.0895  0.0188    -4.7705
## 
## Residual standard error: 0.7608 on 9632 degrees of freedom
durbinWatsonTest(illness_model_robust)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.03896337      1.921995       0
##  Alternative hypothesis: rho != 0

#This is still 0 so running it as a rlm did not fix it. I am uncertain what would!