library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(base)
#library(tigerstats)
library(mosaic)
library(pROC) # for ROC/AUC
library(caret)Predicting Incidence of Heart Disease as a Result of Lifestyle Habits
Lifestyle Habits
A mock data set similar to The Behavioral Risk Factor Surveillance System is being used in this project to create a model that will reliably predict the incidence of heart disease based off a number of ‘lifestyle habits’. Lifestyle habits are defined for the purpose of the study as : smoking, drinking, and the different food consumption.
The Behavioral Risk Factor Surveillance System (BRFSS) is a state-based national system of telephone surveys in the United States that collects data on health-related risk behaviors, chronic health conditions, and the use of preventive services for adults aged 18 or older.
Data Import
#data import and manipulation
cvd <- read_csv("D:/python stuff/CVD_cleaned.csv")
cvd <-as.data.frame(cvd)
cvd <- cvd %>%
rename(
height = `Height_(cm)`,
weight = `Weight_(kg)`,
age = Age_Category,
Fruit = Fruit_Consumption,
Fries = FriedPotato_Consumption,
Alcohol = Alcohol_Consumption,
Smoking = Smoking_History,
Greens = Green_Vegetables_Consumption) %>% mutate(
Heart_Disease = factor(ifelse(Heart_Disease == "Yes", "1", "0"), levels = c("0","1")),
Smoking = factor(ifelse(Smoking == "Yes", 1, 0)),
Exercise = factor(Exercise))
#subset dataframe
lifestyle_df <- cvd[,c(3,4,15:19)]The dataset has 19 variables with a sample size of 308854. There is strong class imbalance between the Heart Disease variable with ‘0’ accounting for ~92% of the dataset.
## Data Quality Check + +`{r message=FALSE, warning=FALSE} +# Missingness by variable +
sum(is.na(cvd))[1] 0lifestyle_df |> count(Heart_Disease) |> mutate(Proportion = round(n / sum(n), 3)) Heart_Disease n Proportion
1 0 283883 0.919
2 1 24971 0.081ggplot(lifestyle_df, aes(x = Heart_Disease, fill = Heart_Disease)) + geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5) +
labs(title = "Heart Disease Class Distribution", x = "Heart Disease", y = "Count") +
scale_y_continuous(labels = scales::label_comma())+
theme_bw() +
theme(legend.position = "none") 
Data Exploration and Visualization
#data exploration
head(cvd) General_Health Checkup Exercise Heart_Disease Skin_Cancer
1 Poor Within the past 2 years No 0 No
2 Very Good Within the past year No 1 No
3 Very Good Within the past year Yes 0 No
4 Poor Within the past year Yes 1 No
5 Good Within the past year No 0 No
6 Good Within the past year No 0 No
Other_Cancer Depression Diabetes Arthritis Sex age height weight BMI
1 No No No Yes Female 70-74 150 32.66 14.54
2 No No Yes No Female 70-74 165 77.11 28.29
3 No No Yes No Female 60-64 163 88.45 33.47
4 No No Yes No Male 75-79 180 93.44 28.73
5 No No No No Male 80+ 191 88.45 24.37
6 No Yes No Yes Male 60-64 183 154.22 46.11
Smoking Alcohol Fruit Greens Fries
1 1 0 30 16 12
2 0 0 30 0 4
3 0 4 12 3 16
4 0 0 30 30 8
5 1 0 8 4 0
6 0 0 12 12 12str(cvd)'data.frame': 308854 obs. of 19 variables:
$ General_Health: chr "Poor" "Very Good" "Very Good" "Poor" ...
$ Checkup : chr "Within the past 2 years" "Within the past year" "Within the past year" "Within the past year" ...
$ Exercise : Factor w/ 2 levels "No","Yes": 1 1 2 2 1 1 2 2 1 1 ...
$ Heart_Disease : Factor w/ 2 levels "0","1": 1 2 1 2 1 1 2 1 1 1 ...
$ Skin_Cancer : chr "No" "No" "No" "No" ...
$ Other_Cancer : chr "No" "No" "No" "No" ...
$ Depression : chr "No" "No" "No" "No" ...
$ Diabetes : chr "No" "Yes" "Yes" "Yes" ...
$ Arthritis : chr "Yes" "No" "No" "No" ...
$ Sex : chr "Female" "Female" "Female" "Male" ...
$ age : chr "70-74" "70-74" "60-64" "75-79" ...
$ height : num 150 165 163 180 191 183 175 165 163 163 ...
$ weight : num 32.7 77.1 88.5 93.4 88.5 ...
$ BMI : num 14.5 28.3 33.5 28.7 24.4 ...
$ Smoking : Factor w/ 2 levels "0","1": 2 1 1 1 2 1 2 2 2 1 ...
$ Alcohol : num 0 0 4 0 0 0 0 3 0 0 ...
$ Fruit : num 30 30 12 30 8 12 16 30 12 12 ...
$ Greens : num 16 0 3 30 4 12 8 8 12 12 ...
$ Fries : num 12 4 16 8 0 12 0 8 4 1 ...life_long <- lifestyle_df %>%
pivot_longer(
cols = c(Greens, Fruit,Fries,Alcohol),
names_to = "Food",
values_to = "Value")
ggplot(life_long,aes(Heart_Disease,y=Value, fill = Food=='Alcohol'))+
geom_col(position = "dodge")+
theme_minimal()+
scale_fill_manual(
values = c(`0` = "#440154",
`1` = "#FDE725"))Warning: No shared levels found between `names(values)` of the manual scale and the
data's fill values.
No shared levels found between `names(values)` of the manual scale and the
data's fill values.
life_long %>% select(Heart_Disease,Smoking,Food,Value) %>% group_by(Heart_Disease,Food) %>% summarise(n = sum(Value)) %>%
ggplot(aes(Food,n, fill=Heart_Disease))+
geom_col( position = 'dodge')+
scale_y_continuous(labels = scales::label_comma())+
theme_bw()`summarise()` has grouped output by 'Heart_Disease'. You can override using the
`.groups` argument.
ggplot(life_long) +
aes(x = Food, y = Value, fill = Heart_Disease) +
geom_col() +
scale_fill_hue(direction = 1) +
theme_minimal() +
facet_wrap(vars(Heart_Disease), scales = "free_y")
#life_long %>% select(Heart_Disease,Smoking) %>% filter(Heart_Disease=='1') %>%
# ggplot(life_long,aes(Heart_Disease,fill=Smoking))+
# geom_bar(position = 'dodge')+facet_grid(~Heart_Disease)# Smoking x Heart Disease
ggplot(life_long) +
aes(x = Heart_Disease, fill = Smoking) +
geom_bar(aes(y = after_stat(count / sum(count))), position = "dodge2") +
scale_fill_manual(values = c(`0` = "#440154", `1` = "#FDE725")) +
labs(title = "Smoking x Heart Disease", y = "Proportion") +
theme_bw() +
facet_wrap(vars(Heart_Disease), scales = "free")
# Exercise x Heart Disease
ggplot(life_long) +
aes(x = Heart_Disease, fill = Exercise) +
geom_bar(aes(y = after_stat(count / sum(count))), position = "dodge2") +
scale_fill_manual(values = c(No = "#440154", Yes = "#FDE725")) +
labs(title = "Exercise x Heart Disease", y = "Proportion") +
theme_bw() +
facet_wrap(vars(Heart_Disease), scales = "free")
From the graphs above, we can see that smoking is much higher represented in the heart disease group. Inversely, given no heart disease ‘no smoking’ represents much higher proportion of the population.
A different trend is observed for exercise such that while exercise is in higher proportion in both Heart Disease populations, it is much more likely for those who do not have heart disease to indicate that they do exercise.
A logistic regression to determine what variables are significant predictors of heart disease diagnosis
# compute class weights to correct for 92/8 imbalance
class_counts <- table(lifestyle_df$Heart_Disease)
class_weights <- 1 / class_counts[lifestyle_df$Heart_Disease]
#class_weights
#main effects
lifestyle_df$Heart_Disease<-relevel(lifestyle_df$Heart_Disease, ref = "0")
life_lm <- glm(Heart_Disease ~ ., data = lifestyle_df,
family = binomial());summary(life_lm)
Call:
glm(formula = Heart_Disease ~ ., family = binomial(), data = lifestyle_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.2029434 0.0173526 -126.952 < 2e-16 ***
ExerciseYes -0.6219205 0.0144355 -43.083 < 2e-16 ***
Smoking1 0.7679578 0.0136314 56.337 < 2e-16 ***
Alcohol -0.0187079 0.0009083 -20.596 < 2e-16 ***
Fruit 0.0002526 0.0002857 0.884 0.377
Greens -0.0022513 0.0004931 -4.566 4.98e-06 ***
Fries -0.0070181 0.0008540 -8.218 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 173478 on 308853 degrees of freedom
Residual deviance: 167313 on 308847 degrees of freedom
AIC: 167327
Number of Fisher Scoring iterations: 5life_lm_reduced <- glm(Heart_Disease ~ . - Fruit, data = lifestyle_df,
family = binomial())
anova(life_lm_reduced, life_lm, test = "Chisq") # no leading +Analysis of Deviance Table
Model 1: Heart_Disease ~ (Exercise + Smoking + Alcohol + Fruit + Greens +
Fries) - Fruit
Model 2: Heart_Disease ~ Exercise + Smoking + Alcohol + Fruit + Greens +
Fries
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 308848 167313
2 308847 167313 1 0.78005 0.3771cbind(minus_fruit = AIC(life_lm_reduced), with_fruit = AIC(life_lm)) minus_fruit with_fruit
[1,] 167325.4 167326.6round(exp(cbind(OR = coef(life_lm_reduced),
confint.default(life_lm_reduced))), 3) OR 2.5 % 97.5 %
(Intercept) 0.111 0.108 0.115
ExerciseYes 0.538 0.523 0.553
Smoking1 2.153 2.097 2.212
Alcohol 0.981 0.980 0.983
Greens 0.998 0.997 0.999
Fries 0.993 0.991 0.995The main effects model indicates that only fruit consumption is not a statistical significant co-variate, and the model could likely be improved by removing it. Smoking has a high odds of heart disease with (~2.15 OR) . Surprisingly, Fruit, Greens and Fries Consumption shows effectively little to no difference in odds of a positive heart disease (1.0, .99, .99,1.00) and the confidence intervals are very close to 1. Alcohol consumption indicates a 2.2% decrease in odds for every 1 extra alcoholic drink taken. Only exercise seems to have a protective effect with a 0.53 (0.980, 0.983) odds of a heart attack compared to those that don’t exercise. The seeming protective effects of fries consumption on heart-disease outcome should be investigated as there may be confounding variable or it may indicate that given a heart disease diagnoses may reduce alcohol and fries intake.
Model Creation and Comparisons
library(broom)
or_df <- tidy(life_lm_reduced, exponentiate = TRUE, conf.int = TRUE) %>%
filter(term != "(Intercept)")
ggplot(or_df, aes(x = estimate, y = reorder(term, estimate))) +
geom_point(size = 3, colour = "#440154") +
geom_errorbar(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
labs(title = "Odds ratios (95% CI)", x = "Odds ratio", y = NULL) +
theme_bw()`height` was translated to `width`.
library(broom)
aug_df <- augment(life_lm_reduced, type.predict = "response", type.residuals = "deviance")
ggplot(aug_df, aes(x = .fitted, y = .resid)) +
geom_point(alpha = 0.1) +
geom_smooth(se = FALSE, colour = "#440154") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
labs(title = "Deviance residuals vs fitted values",
x = "Fitted probability", y = "Deviance residual") +
theme_bw()`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
life_lm_reduced <- glm(Heart_Disease ~ . - Fruit, data = lifestyle_df, family = binomial())
+ anova(life_lm_reduced, life_lm, test = 'Chisq') Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 308848 167313.4 NA NA NA
2 308847 167312.6 1 0.780054 0.3771246#comparing AIC
cbind(minus_fruit=AIC(life_lm_reduced),with_fruit=AIC(life_lm)) minus_fruit with_fruit
[1,] 167325.4 167326.6Removing fruit from the main model improved the model’s AIC score and comparing the models with a chi-square test indicates that they are not statistically different.
life_lm_interact <- glm(Heart_Disease ~ .^2, data = lifestyle_df,
family = binomial())
summary(life_lm_interact)
Call:
glm(formula = Heart_Disease ~ .^2, family = binomial(), data = lifestyle_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.216e+00 2.971e-02 -74.588 < 2e-16 ***
ExerciseYes -6.155e-01 3.215e-02 -19.146 < 2e-16 ***
Smoking1 7.343e-01 3.100e-02 23.689 < 2e-16 ***
Alcohol -1.821e-02 2.647e-03 -6.881 5.96e-12 ***
Fruit 1.847e-03 6.766e-04 2.730 0.006326 **
Greens 2.094e-03 1.248e-03 1.679 0.093216 .
Fries -1.388e-02 2.085e-03 -6.661 2.72e-11 ***
ExerciseYes:Smoking1 1.225e-01 2.950e-02 4.153 3.28e-05 ***
ExerciseYes:Alcohol -2.903e-04 2.031e-03 -0.143 0.886324
ExerciseYes:Fruit -1.660e-03 6.225e-04 -2.667 0.007653 **
ExerciseYes:Greens -2.718e-03 1.073e-03 -2.534 0.011277 *
ExerciseYes:Fries 2.169e-05 1.733e-03 0.013 0.990012
Smoking1:Alcohol -4.690e-04 1.979e-03 -0.237 0.812627
Smoking1:Fruit -3.857e-04 5.771e-04 -0.668 0.503975
Smoking1:Greens -4.616e-03 9.946e-04 -4.641 3.47e-06 ***
Smoking1:Fries 5.149e-03 1.769e-03 2.910 0.003614 **
Alcohol:Fruit -1.142e-05 4.043e-05 -0.283 0.777489
Alcohol:Greens 4.452e-05 6.715e-05 0.663 0.507308
Alcohol:Fries -3.320e-05 1.126e-04 -0.295 0.768028
Fruit:Greens -2.883e-05 1.639e-05 -1.759 0.078629 .
Fruit:Fries 4.453e-05 3.087e-05 1.443 0.149151
Greens:Fries 1.191e-04 3.229e-05 3.687 0.000227 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 173478 on 308853 degrees of freedom
Residual deviance: 167229 on 308832 degrees of freedom
AIC: 167273
Number of Fisher Scoring iterations: 5anova(life_lm_reduced, life_lm_interact, test = "Chisq")Analysis of Deviance Table
Model 1: Heart_Disease ~ (Exercise + Smoking + Alcohol + Fruit + Greens +
Fries) - Fruit
Model 2: Heart_Disease ~ (Exercise + Smoking + Alcohol + Fruit + Greens +
Fries)^2
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 308848 167313
2 308832 167229 16 84.1 3.01e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cbind(AIC(life_lm_reduced), AIC(life_lm_interact)) [,1] [,2]
[1,] 167325.4 167273.3exp(cbind(OR = coef(life_lm), confint.default(life_lm))) OR 2.5 % 97.5 %
(Intercept) 0.1104775 0.1067833 0.1142995
ExerciseYes 0.5369123 0.5219343 0.5523202
Smoking1 2.1553601 2.0985377 2.2137211
Alcohol 0.9814660 0.9797203 0.9832148
Fruit 1.0002526 0.9996927 1.0008128
Greens 0.9977512 0.9967874 0.9987159
Fries 0.9930065 0.9913457 0.9946700ggplot(lifestyle_df, aes(sample = Fruit)) +
stat_qq() +
stat_qq_line(colour = "steelblue") +
labs(title = "Q-Q plot: Fruit consumption") +
theme_bw()
var.test(lifestyle_df$Fruit ~ lifestyle_df$Heart_Disease)
F test to compare two variances
data: lifestyle_df$Fruit by lifestyle_df$Heart_Disease
F = 1.0387, num df = 283882, denom df = 24970, p-value = 5.209e-05
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
1.019821 1.057831
sample estimates:
ratio of variances
1.038726 t.test(lifestyle_df$Fruit ~ lifestyle_df$Heart_Disease,
mu = 0, var.equal = FALSE, alternative = "greater")
Welch Two Sample t-test
data: lifestyle_df$Fruit by lifestyle_df$Heart_Disease
t = 11.326, df = 29720, p-value < 2.2e-16
alternative hypothesis: true difference in means between group 0 and group 1 is greater than 0
95 percent confidence interval:
1.564253 Inf
sample estimates:
mean in group 0 mean in group 1
29.98316 28.15314 var.test(lifestyle_df$Fries ~ lifestyle_df$Heart_Disease)
F test to compare two variances
data: lifestyle_df$Fries by lifestyle_df$Heart_Disease
F = 1.0195, num df = 283882, denom df = 24970, p-value = 0.0397
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
1.000910 1.038215
sample estimates:
ratio of variances
1.019464 t.test(lifestyle_df$Fries ~ lifestyle_df$Heart_Disease,
mu = 0, var.equal = FALSE, alternative = "two.sided")
Welch Two Sample t-test
data: lifestyle_df$Fries by lifestyle_df$Heart_Disease
t = 5.1697, df = 29628, p-value = 2.361e-07
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
0.1803644 0.4006544
sample estimates:
mean in group 0 mean in group 1
6.320104 6.029594 Model Testing
#model test
# train/test split
train_idx <- createDataPartition(lifestyle_df$Heart_Disease, p = 0.8, list = FALSE)
train_df <- lifestyle_df[ train_idx, ]
test_df <- lifestyle_df[-train_idx, ]
train_weights <- class_weights[train_idx]
final_model <- glm(Heart_Disease ~ . - Fruit, data = train_df,
family = binomial())
# predictions on test set
pred_prob <- predict(final_model, newdata = test_df, type = "response")
pred_class <- factor(ifelse(pred_prob > 0.5, "1", "0"), levels = c("0", "1"))
# confusion matrix
confusionMatrix(pred_class, test_df$Heart_Disease, positive = "1")Confusion Matrix and Statistics
Reference
Prediction 0 1
0 56776 4994
1 0 0
Accuracy : 0.9192
95% CI : (0.917, 0.9213)
No Information Rate : 0.9192
P-Value [Acc > NIR] : 0.5038
Kappa : 0
Mcnemar's Test P-Value : <2e-16
Sensitivity : 0.00000
Specificity : 1.00000
Pos Pred Value : NaN
Neg Pred Value : 0.91915
Prevalence : 0.08085
Detection Rate : 0.00000
Detection Prevalence : 0.00000
Balanced Accuracy : 0.50000
'Positive' Class : 1
# ROC / AUC
roc_obj <- roc(test_df$Heart_Disease, pred_prob)Setting levels: control = 0, case = 1Setting direction: controls < casesauc(roc_obj)Area under the curve: 0.6558ggroc(roc_obj, colour = "#440154") +
geom_abline(slope = 1, intercept = 1, linetype = "dashed", colour = "grey60") +
labs(title = paste0("ROC curve (AUC = ", round(auc(roc_obj), 3), ")"),
x = "Specificity", y = "Sensitivity") +
theme_bw()
The model selected produced an AUC value of 0.653 indicating that the model does not do a great job of correctly classifying those with heart disease. The model itself could be much improved through the use of investigating additional , possible confounding, variables.The result is not surpr ising considering that only Exercise and Smoking were the only two covariates with significant odds that were much different from 1.0 in the logistic regression.