The data set that I will be working with is a data set that displays information about patients based on variables that affect heart failure. The main dependent variable is DEATH_EVENT which tells whether a patient died of heart failure or not. The independent variables in this data set are age, anemia, high blood pressure, creatinine phosphokinase, diabetes, ejection fraction, sec, platelets, serum creatinine, serum sodium, smoking, time, and death event.

In this data set, high blood pressure is a boolean variable. 1 means they have high blood pressure and 0 means they do not. Additionally, age is age of the patient at the time the data was taken. Lastly, time is the follow-up period measured in days.

I intend to explore find the best equation to predict the probability of someone dying of heart failure. I also want to see which age_range died the most. Additionally I want to see how quickly people died after the experiment

Load Libraries

#Loading libraries
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── 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(vcd)
## Loading required package: grid
library(vcdExtra)
## Loading required package: gnm
## 
## Attaching package: 'vcdExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     summarise
library(car)
## Loading required package: carData
## 
## Attaching package: 'carData'
## 
## The following object is masked from 'package:vcdExtra':
## 
##     Burt
## 
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:vcdExtra':
## 
##     summarise
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

Load Heart Failure Dataset

#Loading dataset
heart <- read_csv("/Users/nnokigbo/Downloads/heart_failure_clinical_records_dataset.csv")
## Rows: 299 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (13): age, anaemia, creatinine_phosphokinase, diabetes, ejection_fractio...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Viewing all the variables in the dataset
names(heart)
##  [1] "age"                      "anaemia"                 
##  [3] "creatinine_phosphokinase" "diabetes"                
##  [5] "ejection_fraction"        "high_blood_pressure"     
##  [7] "platelets"                "serum_creatinine"        
##  [9] "serum_sodium"             "sex"                     
## [11] "smoking"                  "time"                    
## [13] "DEATH_EVENT"
head(heart)
## # A tibble: 6 × 13
##     age anaemia creatinine_phosphokinase diabetes ejection_fraction
##   <dbl>   <dbl>                    <dbl>    <dbl>             <dbl>
## 1    75       0                      582        0                20
## 2    55       0                     7861        0                38
## 3    65       0                      146        0                20
## 4    50       1                      111        0                20
## 5    65       1                      160        1                20
## 6    90       1                       47        0                40
## # ℹ 8 more variables: high_blood_pressure <dbl>, platelets <dbl>,
## #   serum_creatinine <dbl>, serum_sodium <dbl>, sex <dbl>, smoking <dbl>,
## #   time <dbl>, DEATH_EVENT <dbl>

Finding the best logistic model

Let’s look at a logistic full model and assess which variables are significant to death by heart failure and may lower our AIC by removing them.

# Making sure our response variable is binary
unique(heart$DEATH_EVENT)
## [1] 1 0
# Full logistic regression model
full_model <- glm(data = heart, DEATH_EVENT ~ high_blood_pressure + diabetes + platelets + age + anaemia + creatinine_phosphokinase + ejection_fraction + serum_creatinine + serum_sodium + sex + smoking + time, family = binomial())
summary(full_model)
## 
## Call:
## glm(formula = DEATH_EVENT ~ high_blood_pressure + diabetes + 
##     platelets + age + anaemia + creatinine_phosphokinase + ejection_fraction + 
##     serum_creatinine + serum_sodium + sex + smoking + time, family = binomial(), 
##     data = heart)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.018e+01  5.657e+00   1.801 0.071774 .  
## high_blood_pressure      -1.027e-01  3.587e-01  -0.286 0.774688    
## diabetes                  1.451e-01  3.512e-01   0.413 0.679380    
## platelets                -1.200e-06  1.889e-06  -0.635 0.525404    
## age                       4.742e-02  1.580e-02   3.001 0.002690 ** 
## anaemia                  -7.470e-03  3.605e-01  -0.021 0.983467    
## creatinine_phosphokinase  2.222e-04  1.779e-04   1.249 0.211684    
## ejection_fraction        -7.666e-02  1.633e-02  -4.695 2.67e-06 ***
## serum_creatinine          6.661e-01  1.815e-01   3.670 0.000242 ***
## serum_sodium             -6.698e-02  3.974e-02  -1.686 0.091855 .  
## sex                      -5.337e-01  4.139e-01  -1.289 0.197299    
## smoking                  -1.349e-02  4.126e-01  -0.033 0.973915    
## time                     -2.104e-02  3.014e-03  -6.981 2.92e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 375.35  on 298  degrees of freedom
## Residual deviance: 219.55  on 286  degrees of freedom
## AIC: 245.55
## 
## Number of Fisher Scoring iterations: 6

We have an AIC of 245.55. Some variables with large p-values such as smoking, anemia, high blood pressure, diabetes, and platelets. Let’s remove them.

#Making second model
model2 <- glm(data = heart, DEATH_EVENT ~ age + creatinine_phosphokinase + ejection_fraction + serum_creatinine + serum_sodium + sex + time, family = binomial())
summary(model2)
## 
## Call:
## glm(formula = DEATH_EVENT ~ age + creatinine_phosphokinase + 
##     ejection_fraction + serum_creatinine + serum_sodium + sex + 
##     time, family = binomial(), data = heart)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              10.3375665  5.5190911   1.873 0.061061 .  
## age                       0.0464158  0.0155516   2.985 0.002839 ** 
## creatinine_phosphokinase  0.0002216  0.0001761   1.258 0.208212    
## ejection_fraction        -0.0761854  0.0162053  -4.701 2.59e-06 ***
## serum_creatinine          0.6688599  0.1782683   3.752 0.000175 ***
## serum_sodium             -0.0702160  0.0389939  -1.801 0.071751 .  
## sex                      -0.5126188  0.3669166  -1.397 0.162384    
## time                     -0.0208448  0.0029326  -7.108 1.18e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 375.35  on 298  degrees of freedom
## Residual deviance: 220.22  on 291  degrees of freedom
## AIC: 236.22
## 
## Number of Fisher Scoring iterations: 6

We successfully achieved a lower AIC of 236.22. Let’s see what else we can take off. The largest p-values in this set include sex and creatinine phosphokinase. Let’s remove them as see how our AIC is affected.

# Making third model
model3 <- glm(data = heart, DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + serum_sodium + time, family = binomial())
summary(model3)
## 
## Call:
## glm(formula = DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + 
##     serum_sodium + time, family = binomial(), data = heart)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        9.493034   5.405768   1.756  0.07907 .  
## age                0.042466   0.015030   2.825  0.00472 ** 
## ejection_fraction -0.073430   0.015785  -4.652 3.29e-06 ***
## serum_creatinine   0.685990   0.174044   3.941 8.10e-05 ***
## serum_sodium      -0.064557   0.038377  -1.682  0.09254 .  
## time              -0.020895   0.002916  -7.166 7.74e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 375.35  on 298  degrees of freedom
## Residual deviance: 223.49  on 293  degrees of freedom
## AIC: 235.49
## 
## Number of Fisher Scoring iterations: 6

Our AIC lowered again. Let’s remove serum_sodium and see if we can lower our AIC even more.

# Making another model
model4 <- glm(data = heart, DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + time, family = binomial())
summary(model4)
## 
## Call:
## glm(formula = DEATH_EVENT ~ age + ejection_fraction + serum_creatinine + 
##     time, family = binomial(), data = heart)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.604473   1.036111   0.583  0.55962    
## age                0.043326   0.014872   2.913  0.00358 ** 
## ejection_fraction -0.074804   0.015555  -4.809 1.52e-06 ***
## serum_creatinine   0.719785   0.174597   4.123 3.75e-05 ***
## time              -0.020611   0.002881  -7.153 8.48e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 375.35  on 298  degrees of freedom
## Residual deviance: 226.30  on 294  degrees of freedom
## AIC: 236.3
## 
## Number of Fisher Scoring iterations: 5

Our AIC increased so I’ve decided it’s best to stay with model3 as our best model.

Now that we’ve found a successful model, let’s check for co-linearity

#Checking for co-linearity
vif(model3)
##               age ejection_fraction  serum_creatinine      serum_sodium 
##          1.053111          1.133484          1.079122          1.028355 
##              time 
##          1.096862

All the variables in my model are below 4 so it is safe to assume multi-colinearity is not an issue. Yay!

Predicintg Heart Failure

It’s safe to say that model3 is a good model. The equation for model3 is: *log_odds(DEATH_EVENT) = -0.042(age) - 0.073(ejection_fraction) + 0.685(serum_creatinine) - 0.064(serum_sodium) - 0.0208(time)

Let’s calculate the probability of dying of heart failure given a patients age is 43, they have an ejection fraction of 32, their creatinine serum dose was 1.50, their sodium serum dose was 139 and their follow up time was 32 days.

# Calulating log odds of dying of heart failure with given data
logodds <- (-0.042*42) - (0.073*32) + (0.685*1.50) - (0.064*139) - (0.0208*32)
logodds
## [1] -12.6341
# Calulating probability of dying of heart failure with given data
exp(logodds)/(1+exp(logodds))
## [1] 3.258957e-06

The probability of dying of heart failure is 0.000325%

Death Event and Age

Plotting Age and DEATH_EVENT

heart1 <- heart |>
  mutate(sex = ifelse(sex == 1, "Men", "Women"))
log_plot <- ggplot(data = heart1, aes(y = DEATH_EVENT, x = age, shape = sex)) + 
  geom_jitter(width = 0, height = 0.05, alpha = 0.5) +
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial"))
log_plot
## `geom_smooth()` using formula = 'y ~ x'

I added a sigmoid curve, and we can see that as age increases, the possibility for a death via heart failure increases.

Adding a New Variable called age_range

Now seeing that age was important in factors contributing to death by heart failure, I want to see how age groups varied when it came to dying by heart failure.

# Seeing the minimum and maximum for age
summary(heart$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   40.00   51.00   60.00   60.83   70.00   95.00
# Making age a ctegorical variable by separating them into 6 groups of 10
heart_mod <- heart |>
  filter(!is.na(age)) |> #Getting rid of Na's
  mutate(age_range = case_when((age >= 40 & age <= 50) ~ "40-50",
                              (age > 50 & age <= 60) ~ "51-60",
                              (age > 60 & age <= 70) ~ "61-70",
                              (age > 70 & age <= 80) ~ "71-80",
                              (age > 80 & age <= 90) ~ "81-90",
                              (age > 90 & age <= 100) ~ "91-100",)) |>
  mutate(diedOrNot = ifelse(DEATH_EVENT == 1, "Died", "Did Not Die")) # Changing variable to more readable valiues

Plotting death event with age

# Plotting death with age in a bar graph
ageplot <- ggplot(heart_mod, aes(x = diedOrNot, fill = age_range)) +
  geom_bar() +
  scale_fill_brewer(palette = "PuRd", name = "Age Range") +
  labs(
    title = "Counts of Deaths vs Non-Deaths based on Age",
    y = "Amount of Deaths and No Deaths",
    x = "Death Event",
    caption = "Collected at the Faisalabad Institute of Cardiology\n and at the Allied Hospital in Faisalabad\n (Punjab, Pakistan)"
  ) + theme_classic()

ageplot

# Adding interactivity
ageplot <- ggplotly(ageplot)
ageplot

Death Event and Time

Plotting Time and DEATH_EVENT

heart2 <- heart |>
  mutate(sex = ifelse(sex == 1, "Men", "Women"))
log_plot2 <- ggplot(data = heart2, aes(y = DEATH_EVENT, x = time, shape = sex)) + 
  geom_jitter(width = 0, height = 0.05, alpha = 0.5) +
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial"))
log_plot2
## `geom_smooth()` using formula = 'y ~ x'

This graph shows that as time goes on, the possibility of a death event decreases. Most deaths caused by heart failure happen in the beginning.

Adding a New Variable called time_range

Now I’ll be separating the variable time by 30 days and defining them as months

summary(heart$time)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     4.0    73.0   115.0   130.3   203.0   285.0
heart_time <- heart |>
  filter(!is.na(time)) |> #Getting rid of Na's
  #Separating the variable by months
  mutate(time_range = case_when((time >= 1 & time <= 30) ~ "First Month",
                              (time > 30 & time <= 60) ~ "Second Month",
                              (time > 60 & time <= 90) ~ "Third Month",
                              (time > 90 & time <= 120) ~ "Fourth Month",
                              (time > 120 & time <= 150) ~ "Fifth Month",
                              (time > 150 & time <= 180) ~ "Sixth Month",
                              (time > 180 & time <= 210) ~ "Seventh Month",
                              (time > 210 & time <= 240) ~ "Eighth Month",
                              (time > 240 & time <= 270) ~ "Ninth Month",
                              (time > 270 & time <= 300) ~ "Tenth Month",)) |>
  mutate(diedOrNot = ifelse(DEATH_EVENT == 1, "Died", "Did Not Die")) # Changing variable to more readable valiues

Plotting death event and time range in a bar graph.

# Plotting death with time in a bar graph
timeplot <- ggplot(heart_time, aes(x = diedOrNot, fill = time_range)) +
  geom_bar() +
  scale_fill_brewer(palette = "PRGn", name = "Time Range", breaks=c('First Month', 'Second Month', 'Third Month', 'Fourth Month', 'Fifth Month', 'Sixth Month', 'Seventh Month', 'Eighth Month', 'Ninth Month', 'Tenth Month')) + #Reorganizing Legend
  labs(
    title = "Counts of Deaths vs Non-Deaths based on Time",
    y = "Amount of Deaths and No Deaths",
    x = "Death Event",
    caption = "Collected at the Faisalabad Institute of Cardiology\n and at the Allied Hospital in Faisalabad\n (Punjab, Pakistan)"
  ) + theme_classic()

timeplot

timeplot <- ggplotly(timeplot) #Adding interactivity
timeplot

End Essay

I did not do much of anything to cleant he data up because the data was already clean. However, I did add two know variables time_range and age_range to make age and time a categorical variable. My first visualization comparing age_range and death event represents the amount of people who died and did not die based on age range. One thing I found interesting that the younger people in the dataset died the most of heart failure. However this could be because younger people took up the majority of the data points in the data set. Another thing I found that was interesting that most of the people who died, died very early on after the experiment, typically a month or two after. I wish that I had got to include the amount of people who died per age range and time range per capita. The data is a bit misleading because as people in the 40-60 age range died more, they also were the largest group of people in the data set.