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):

  1. Pregnancies: Number of times pregnant
  2. Glucose: Plasma glucose concentration (2-hour test)
  3. BloodPressure: Diastolic blood pressure (mm Hg)
  4. SkinThickness: Triceps skin fold thickness (mm)
  5. Insulin: 2-hour serum insulin (mu U/ml)
  6. BMI: Body mass index (weight in kg/(height in m)^2)
  7. DiabetesPedigreeFunction: Diabetes pedigree function (a function scoring genetic risk)
  8. Age: Age in years
  9. 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.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.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
url <- "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                        0

Question 1: Create and Interpret a Logistic Regression Model - Fit a logistic regression model to predict Outcome using Glucose, BMI, and Age.

#Find number of rows where Glucose OR BMI is NA to see if it is small enough to safely remove
sum(is.na(data$Glucose) | is.na(data$BMI))
## [1] 16
# 16/768 = about 2% of data, so it can be removed safely w/o significantly affecting the analysis
data <- data[!(is.na(data$Glucose) | is.na(data$BMI)),] # deleting the NAs

#Fit a logistic regression model and provide a summary
logistic <- glm(Outcome ~ Glucose + BMI + Age, data=data, family="binomial")
summary(logistic)
## 
## 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
## AIC: 732.96
## 
## Number of Fisher Scoring iterations: 4
#Calculate R^2
1 - (logistic$deviance / logistic$null.deviance)
## [1] 0.25626

The \(R^2\) of about 0.25626 indicates a low goodness-of fit. About 26% of the variance in diabetes outcomes can be explained by the model, and the rest is a matter of random outside factors.

What does the intercept represent (log-odds of diabetes when predictors are zero)?

The log(odds) that a person will have diabetes if their BMI, Glucose, and Age are 0. In reality, this is nonsense, but it serves as an anchor point for the model.

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)?

A one-unit increase in Glucose, BMI, or Age increases the odds of diabetes according to the model. All three predictors are significant because they all have a p-value of less than 0.05.

Question 2: Confusion Matrix and Important Metric

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_probs = data.frame(probability_of_diabetes = logistic$fitted.values, Outcome_num=data_subset$Outcome_num)

# Predicted classes
predicted_probs <- predicted_probs |>
  mutate(Predicted_Outcome_num = ifelse(probability_of_diabetes > 0.5, 1, 0))

# Confusion matrix
confusion_matrix <- 
  data.frame(true_positive = sum(predicted_probs$Predicted_Outcome_num == 1 & predicted_probs$Outcome_num == 1),
             true_negative = sum(predicted_probs$Predicted_Outcome_num == 0 & predicted_probs$Outcome_num == 0),
             false_positive = sum(predicted_probs$Predicted_Outcome_num == 1 & predicted_probs$Outcome_num == 0),
             false_negative = sum(predicted_probs$Predicted_Outcome_num == 0 & predicted_probs$Outcome_num == 1))
#Extract Values:
TN <- confusion_matrix$true_negative
FP <- confusion_matrix$false_positive
FN <- confusion_matrix$false_negative
TP <- confusion_matrix$true_positive

#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.718

Interpret: How well does the model perform? Is it better at detecting diabetes (sensitivity) or non-diabetes (specificity)? Why might this matter for medical diagnosis?

The model performs fairly well. Based on the Accuracy, about 77% of the model’s predictions are correct. The model is much better at predicting non-diabetes than diabetes because it has a higher specificity than sensitivity. About 88% of non-diabetics were correctly predicted, whereas the model only correctly predicted about 57% of diabetics. This would matter for medical diagnosis because the test will often produce false negatives. Doctors should not withhold treatment for diabetes based solely on a prediction from this model.

Question 3: ROC Curve, AUC, and Interpretation

# install.packages("pROC") # if needed
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
# ROC curve
roc_obj <- roc(response = data_subset$Outcome,
               predictor = logistic$fitted.values)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Plot ROC
plot.roc(roc_obj, print.auc = TRUE,
         xlab = "False Positive Rate (1 - Specificity)",
         ylab = "True Positive Rate (Sensitivity)")

# Print AUC value
auc_val <- auc(roc_obj); auc_val
## Area under the curve: 0.828

What does AUC indicate (0.5 = random, 1.0 = perfect)?

There is about an 83% chance that the model gives a random positive outcome a higher probability than a random negative outcome. That means the model performs much better than random chance, but is not close to being 100% accurate.

For diabetes diagnosis, prioritize sensitivity (catching cases) or specificity (avoiding false positives)? Suggest a threshold and explain.

In my opinion, one should prioritize specificity because of the scarcity of diabetes drugs. 0.6 could be a possible threshold that would reduce the number of false positives and ensure that the insulin supply went only to those who certainly had the disease. To me, this seems like the best management of risk because there are extreme negative health consequences for a non-diabetic recieving diabetes treatment and a diabetic not receiving diabetes treatment.