HW 10
Javier Mantilla
Introduction:
In this homework, you will apply logistic regression to a real-world dataset: the Pima Indians Diabetes Database. This dataset contains medical records from 768 women of Pima Indian heritage, aged 21 or older, and is used to predict the onset of diabetes (binary outcome: 0 = no diabetes, 1 = diabetes) based on physiological measurements.
The data is publicly available from the UCI Machine Learning Repository and can be imported directly.
Dataset URL: https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.data.csv
Columns (no header in the CSV, so we need to assign them manually):
- Pregnancies: Number of times pregnant
- Glucose: Plasma glucose concentration (2-hour test)
- BloodPressure: Diastolic blood pressure (mm Hg)
- SkinThickness: Triceps skin fold thickness (mm)
- Insulin: 2-hour serum insulin (mu U/ml)
- BMI: Body mass index (weight in kg/(height in m)^2)
- DiabetesPedigreeFunction: Diabetes pedigree function (a function scoring genetic risk)
- Age: Age in years
- Outcome: Class variable (0 = no diabetes, 1 = diabetes)
Task Overview: You will load the data, build a logistic regression model to predict diabetes onset using a subset of predictors (Glucose, BMI, Age), interpret the model, evaluate it with a confusion matrix and metrics, and analyze the ROC curve and AUC.
Cleaning the dataset Don’t change the following code
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.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 errorsurl <- "https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.data.csv"
data <- read.csv(url, header = FALSE)
colnames(data) <- c("Pregnancies", "Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI", "DiabetesPedigreeFunction", "Age", "Outcome")
data$Outcome <- as.factor(data$Outcome)
# Handle missing values (replace 0s with NA because 0 makes no sense here)
data$Glucose[data$Glucose == 0] <- NA
data$BloodPressure[data$BloodPressure == 0] <- NA
data$BMI[data$BMI == 0] <- NA
colSums(is.na(data))## Pregnancies Glucose BloodPressure
## 0 5 35
## SkinThickness Insulin BMI
## 0 0 11
## DiabetesPedigreeFunction Age Outcome
## 0 0 0Question 1: Create and Interpret a Logistic Regression Model - Fit a logistic regression model to predict Outcome using Glucose, BMI, and Age.
Provide the model summary.
Calculate and interpret R²: 1 - (model\(deviance / model\)null.deviance). What does it indicate about the model’s explanatory power?
## Enter your code here
logistic_regression <- glm(Outcome ~ Glucose + BMI + Age, data = data, family = "binomial")
summary(logistic_regression)##
## Call:
## glm(formula = Outcome ~ Glucose + BMI + Age, family = "binomial",
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.032377 0.711037 -12.703 < 2e-16 ***
## Glucose 0.035548 0.003481 10.212 < 2e-16 ***
## BMI 0.089753 0.014377 6.243 4.3e-10 ***
## Age 0.028699 0.007809 3.675 0.000238 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 974.75 on 751 degrees of freedom
## Residual deviance: 724.96 on 748 degrees of freedom
## (16 observations deleted due to missingness)
## AIC: 732.96
##
## Number of Fisher Scoring iterations: 4R-Squared
r_squared <- 1 - (logistic_regression$deviance/logistic_regression$null.deviance)
r_squared## [1] 0.25626For this model, we got an R-squared of 0.25626, or 25.626%. This means that 25.626% of the variation in whether or not one get’s diabetes (outcome variable) is predicted by age, glucose and BMI.
What does the intercept represent (log-odds of diabetes when predictors are zero)?
The intercept of -9.032377 shows the log-odds of diabetes when all of the predictors, (Age, BMI and Glucose) are zero.
For each predictor (Glucose, BMI, Age), does a one-unit increase raise or lower the odds of diabetes? Are they significant (p-value < 0.05)?
For each predictor, a one-unit increase raises the odds of diabetes because all of their coefficients are positive. All of our variables are significant as their p-values are all less than 0.05.
Glucose: For glucose, our coefficient is 0.035548, which is greater than 0 meaning that a one-unit increase raises the odds of diabetes. Glucose is statistically significant as it’s p-value is < 2e-16 , which is less than our alpha of 0.05.
BMI: For BMI, our coefficient is 0.089753, which is greater than 0 meaning that a one-unit increase raises the odds of diabetes. BMI is statistically significant as it’s p-value is 4.3e-10 , which is less than our alpha of 0.05.
Age: For age, our coefficient is 0.028699, which is greater than 0 meaning that a one-unit increase raises the odds of diabetes. Age is statistically significant as it’s p-value is 0.000238, which is less than our alpha of 0.05.
Question 2: Confusion Matrix and Important Metric
Predict probabilities using the fitted model.
Create predicted classes with a 0.5 threshold (1 if probability > 0.5, else 0).
Build a confusion matrix (Predicted vs. Actual Outcome).
Calculate and report the metrics:
Accuracy: (TP + TN) / Total Sensitivity (Recall): TP / (TP + FN) Specificity: TN / (TN + FP) Precision: TP / (TP + FP)
Use the following starter code
# Keep only rows with no missing values in Glucose, BMI, or Age
data_subset <- data[complete.cases(data[, c("Glucose", "BMI", "Age")]), ]
#Create a numeric version of the outcome (0 = no diabetes, 1 = diabetes).This is required for calculating confusion matrices.
data_subset$Outcome_num <- ifelse(data_subset$Outcome == "1", 1, 0)
# Predicted probabilities
predicted.data <- data.frame(
probability.of.outcome=logistic_regression$fitted.values,
outcome=data_subset$Outcome)
predicted.data <- predicted.data[
order(predicted.data$probability.of.outcome, decreasing=FALSE),]
predicted.data$rank <- 1:nrow(predicted.data)
head(predicted.data)## probability.of.outcome outcome rank
## 681 0.01422697 0 1
## 63 0.01490161 0 2
## 618 0.01550448 0 3
## 98 0.01718930 0 4
## 91 0.02040067 0 5
## 590 0.02132939 0 6# Predicted classes
predicted.probs <- logistic_regression$fitted.values
predicted.classes <- ifelse(predicted.probs > 0.5, 1, 0)
# Confusion matrix
confusion_matrix <- table(
Predicted = factor(predicted.classes, levels = c(0, 1)),
Actual = factor(data_subset$Outcome_num, levels = c(0, 1))
)
confusion_matrix## Actual
## Predicted 0 1
## 0 429 114
## 1 59 150#Extract Values:
TN <- 429
FP <- 59
FN <- 114
TP <- 150
#Metrics
accuracy <- (TP + TN) / (TP +TN + FP + FN)
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
precision <- TP / (TP + FP)
cat("Accuracy:", round(accuracy, 3), "\nSensitivity:", round(sensitivity, 3), "\nSpecificity:", round(specificity, 3), "\nPrecision:", round(precision, 3))## Accuracy: 0.77
## Sensitivity: 0.568
## Specificity: 0.879
## Precision: 0.718Interpret: How well does the model perform? Is it better at detecting diabetes (sensitivity) or non-diabetes (specificity)? Why might this matter for medical diagnosis?
Accuracy: This model performs with 77% accuracy, which is pretty decent in general. However, in context this model is used for health reasons and having a 77% accuracy in predicting diabetes is not good at all in the medical field as it needs to be over 90% accuracy at least.
Sensitivity: This model predicts sensitivity at a 56.8% rate, meaning that it correctly identifies only 56.8% of people with diabetes.
Specificity: This model predicts specificity at a 87.9% rate, meaning that it correctly identifies 87.9% of people without diabetes.
This model is better at detecting non-diabetes (specificity) rather than diabetes (sensitivity) because it predicts non-diabetes at 87.9% rate in comparison to diabetes at 56.8% rate.
This may matter for medical diagnosis because the low sensitivity will lead to false negatives. This will matter because this model’s purpose is to predict diabetes, so in order to make sure the patients are properly diagnosed they will need to raise the sensitivity to make sure they properly diagnose patients.
Question 3: ROC Curve, AUC, and Interpretation
Plot the ROC curve, use the “data_subset” from Q2.
Calculate AUC.
#Enter your code here
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var# ROC curve & AUC on full data
roc_obj <- roc(response = data_subset$Outcome_num,
predictor = logistic_regression$fitted.values,
levels = c(0, 1),
direction = "<") # smaller prob = Healthy
# Print AUC value
auc_val <- auc(roc_obj); auc_val## Area under the curve: 0.828# Plot ROC with AUC displayed
plot.roc(roc_obj, print.auc = TRUE, legacy.axes = TRUE,
xlab = "False Positive Rate (1 - Specificity)",
ylab = "True Positive Rate (Sensitivity)")
What does AUC indicate (0.5 = random, 1.0 = perfect)?
- AUC: The area under the curve for this model is 0.828, or 82.8%. This is much better than random because it is greater than 50%, but not as close to 1.0 or perfect as we would like it. Still good though.
For diabetes diagnosis, prioritize sensitivity (catching cases) or specificity (avoiding false positives)? Suggest a threshold and explain.
For diabetes cases, I think that it is best to prioritize sensitivity (catching cases) rather than specificity (avoiding false positives) because it is important to catch as many true diabetes cases as possible than to avoid false positives because if missing a diagnosis will delay treatment and cause serious health issues because diabetes was not diagnosed. A threshold for sensitivity should be 97%, which is higher than the normal 95% because I think that this extra 2% will allow more people to live and get diagnosed properly, which will save lives in the long run.