LOGISTIC REGRESSION WITH EQUITY ANALYSIS

WEEK 10 IN-CLASS LAB

Author

Tess Vu

Published

November 10, 2025

Logistic Regression with Equity Analysis

Week 10: MUSA 5080 Public Policy Analytics

This script demonstrates how to:

Build a logistic regression model for binary classification
Evaluate model performance using multiple metrics
Analyze disparate impact across demographic groups
Test different decision thresholds
Make evidence-based policy recommendations

SETUP

Code

# Load required packages

library(tidyverse) # For data manipulation and visualization
library(caret) # For model training and confusion matrices library
library(pROC) # For ROC curves and AUC

# Set random seed for reproducibility

set.seed(2025)

# Configure visualization theme
theme_set(theme_minimal())

PART 1: DATA LOADING AND EXPLORATION

Code

## Load Georgia Department of Corrections recidivism data

# Source: National Institute of Justice Recidivism Challenge

recidivism_data <- read_csv("data/NIJ_s_Recidivism_Challenge_Full_Dataset_20240407.csv") # Examine data structure glimpse(recidivism_data)

# Check outcome variable distribution

table(recidivism_data$Recidivism_Within_3years, useNA = "ifany")


FALSE  TRUE 
10931 14904

Code

# Calculate overall recidivism rate

mean(recidivism_data$Recidivism_Within_3years, na.rm = TRUE)

[1] 0.5768918

Code

# Examine recidivism rates by race

recidivism_data %>%
  filter(!is.na(Race), !is.na(Recidivism_Within_3years)) %>%
  group_by(Race) %>%
  summarise(n = n(),
            recidivism_rate = mean(Recidivism_Within_3years),
            avg_age = mean(Age_at_Release, na.rm = TRUE),
            avg_prior_felonies = mean(Prior_Arrest_Episodes_Felony, na.rm = TRUE)
            ) %>%
  arrange(desc(recidivism_rate))

# A tibble: 2 × 5
  Race      n recidivism_rate avg_age avg_prior_felonies
  <chr> <int>           <dbl>   <dbl>              <dbl>
1 BLACK 14847           0.587      NA                 NA
2 WHITE 10988           0.563      NA                 NA

CRITICAL QUESTION: Why do base rates differ across groups?

Consider both individual factors and systemic factors (differential policing, economic opportunities, program access, etc.)

PART 2: DATA PREPARATION

Code

# Clean and prepare data for modeling

model_data <- recidivism_data %>% # Create binary outcome variable (must be numeric for glm)
  mutate(recidivism = as.integer(Recidivism_Within_3years == TRUE)) %>% # Select relevant features
  select(recidivism, # Outcome
         Age_at_Release, # Demographics
         Gender, Race, Gang_Affiliated, # Risk factors
         Dependents, Prior_Arrest_Episodes_Felony, # Criminal history
         Prior_Arrest_Episodes_Violent, Prior_Conviction_Episodes_Prop, Condition_MH_SA, # Mental health / substance abuse
         Supervision_Risk_Score_First, # Official risk score
         Percent_Days_Employed, # Economic factors
         Education_Level # Education
         ) %>%
  # Remove missing values (in practice, consider imputation)
  na.omit()

# Check final sample characteristics

cat("Final sample size:", nrow(model_data), "individualsn")

Final sample size: 21837 individualsn

Code

cat("Recidivism rate:", round(mean(model_data$recidivism), 3), "nn")

Recidivism rate: 0.604 nn

Code

# Sample size by race

model_data %>%
  count(Race) %>%
  arrange(desc(n))

# A tibble: 2 × 2
  Race      n
  <chr> <int>
1 BLACK 13220
2 WHITE  8617

PART 3: TRAIN-TEST SPLIT

Code

# Create stratified train-test split (maintains outcome distribution)

trainIndex <- createDataPartition(y = model_data$recidivism, p = 0.70, # 70% training, 30% testing
                                  list = FALSE )

train_data <- model_data[trainIndex, ]
test_data <- model_data[-trainIndex, ]

# Verify split preserves outcome distribution

cat("Training set:n")

Training set:n

Code

cat(" N =", nrow(train_data), "n")

 N = 15286 n

Code

cat(" Recidivism rate =", round(mean(train_data$recidivism), 3), "nn")

 Recidivism rate = 0.606 nn

Code

cat("Test set:n")

Test set:n

Code

cat(" N =", nrow(test_data), "n")

 N = 6551 n

Code

cat(" Recidivism rate =", round(mean(test_data$recidivism), 3), "nn")

 Recidivism rate = 0.599 nn

PART 4: MODEL TRAINING

Code

# Fit logistic regression model

# Note: We're excluding Race from predictors to avoid direct discrimination,

# but this doesn't eliminate bias (proxy variables may exist)

logit_model <- glm(recidivism ~
                     Age_at_Release + Dependents + Gang_Affiliated + Prior_Arrest_Episodes_Felony +
                     Prior_Arrest_Episodes_Violent + Prior_Conviction_Episodes_Prop + Percent_Days_Employed +
                     Supervision_Risk_Score_First, data = train_data, family = "binomial") # Specifies logistic regression

# View model summary

summary(logit_model)


Call:
glm(formula = recidivism ~ Age_at_Release + Dependents + Gang_Affiliated + 
    Prior_Arrest_Episodes_Felony + Prior_Arrest_Episodes_Violent + 
    Prior_Conviction_Episodes_Prop + Percent_Days_Employed + 
    Supervision_Risk_Score_First, family = "binomial", data = train_data)

Coefficients:
                                         Estimate Std. Error z value Pr(>|z|)
(Intercept)                              1.536623   0.205913   7.462 8.49e-14
Age_at_Release23-27                     -0.306180   0.078571  -3.897 9.75e-05
Age_at_Release28-32                     -0.626801   0.082161  -7.629 2.37e-14
Age_at_Release33-37                     -0.884262   0.087578 -10.097  < 2e-16
Age_at_Release38-42                     -1.097585   0.095083 -11.543  < 2e-16
Age_at_Release43-47                     -1.328051   0.099313 -13.372  < 2e-16
Age_at_Release48 or older               -1.774217   0.098409 -18.029  < 2e-16
Dependents1                              0.112010   0.051801   2.162 0.030593
Dependents2                              0.118375   0.055002   2.152 0.031382
Dependents3 or more                     -0.022584   0.047740  -0.473 0.636171
Gang_AffiliatedTRUE                      0.760899   0.057000  13.349  < 2e-16
Prior_Arrest_Episodes_Felony1           -1.262064   0.197966  -6.375 1.83e-10
Prior_Arrest_Episodes_Felony10 or more   0.484253   0.197974   2.446 0.014444
Prior_Arrest_Episodes_Felony2           -0.723727   0.196190  -3.689 0.000225
Prior_Arrest_Episodes_Felony3           -0.513387   0.196169  -2.617 0.008869
Prior_Arrest_Episodes_Felony4           -0.320167   0.196595  -1.629 0.103406
Prior_Arrest_Episodes_Felony5           -0.201793   0.197509  -1.022 0.306928
Prior_Arrest_Episodes_Felony6           -0.061271   0.199735  -0.307 0.759025
Prior_Arrest_Episodes_Felony7           -0.180865   0.201464  -0.898 0.369318
Prior_Arrest_Episodes_Felony8            0.046340   0.205038   0.226 0.821197
Prior_Arrest_Episodes_Felony9            0.137792   0.209034   0.659 0.509777
Prior_Arrest_Episodes_Violent1           0.070370   0.044557   1.579 0.114264
Prior_Arrest_Episodes_Violent2           0.109302   0.056685   1.928 0.053827
Prior_Arrest_Episodes_Violent3 or more   0.167071   0.058030   2.879 0.003989
Prior_Conviction_Episodes_Prop1          0.144881   0.047901   3.025 0.002490
Prior_Conviction_Episodes_Prop2          0.230322   0.061184   3.764 0.000167
Prior_Conviction_Episodes_Prop3 or more  0.431137   0.061050   7.062 1.64e-12
Percent_Days_Employed                   -1.124559   0.043608 -25.788  < 2e-16
Supervision_Risk_Score_First             0.022286   0.009194   2.424 0.015352
                                           
(Intercept)                             ***
Age_at_Release23-27                     ***
Age_at_Release28-32                     ***
Age_at_Release33-37                     ***
Age_at_Release38-42                     ***
Age_at_Release43-47                     ***
Age_at_Release48 or older               ***
Dependents1                             *  
Dependents2                             *  
Dependents3 or more                        
Gang_AffiliatedTRUE                     ***
Prior_Arrest_Episodes_Felony1           ***
Prior_Arrest_Episodes_Felony10 or more  *  
Prior_Arrest_Episodes_Felony2           ***
Prior_Arrest_Episodes_Felony3           ** 
Prior_Arrest_Episodes_Felony4              
Prior_Arrest_Episodes_Felony5              
Prior_Arrest_Episodes_Felony6              
Prior_Arrest_Episodes_Felony7              
Prior_Arrest_Episodes_Felony8              
Prior_Arrest_Episodes_Felony9              
Prior_Arrest_Episodes_Violent1             
Prior_Arrest_Episodes_Violent2          .  
Prior_Arrest_Episodes_Violent3 or more  ** 
Prior_Conviction_Episodes_Prop1         ** 
Prior_Conviction_Episodes_Prop2         ***
Prior_Conviction_Episodes_Prop3 or more ***
Percent_Days_Employed                   ***
Supervision_Risk_Score_First            *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 20500  on 15285  degrees of freedom
Residual deviance: 17780  on 15257  degrees of freedom
AIC: 17838

Number of Fisher Scoring iterations: 4

Code

# Interpret coefficients as odds ratios

exp(coef(logit_model))

                            (Intercept)                     Age_at_Release23-27 
                              4.6488629                               0.7362542 
                    Age_at_Release28-32                     Age_at_Release33-37 
                              0.5342983                               0.4130187 
                    Age_at_Release38-42                     Age_at_Release43-47 
                              0.3336758                               0.2649933 
              Age_at_Release48 or older                             Dependents1 
                              0.1696161                               1.1185242 
                            Dependents2                     Dependents3 or more 
                              1.1256661                               0.9776691 
                    Gang_AffiliatedTRUE           Prior_Arrest_Episodes_Felony1 
                              2.1401994                               0.2830692 
 Prior_Arrest_Episodes_Felony10 or more           Prior_Arrest_Episodes_Felony2 
                              1.6229618                               0.4849416 
          Prior_Arrest_Episodes_Felony3           Prior_Arrest_Episodes_Felony4 
                              0.5984653                               0.7260279 
          Prior_Arrest_Episodes_Felony5           Prior_Arrest_Episodes_Felony6 
                              0.8172644                               0.9405683 
          Prior_Arrest_Episodes_Felony7           Prior_Arrest_Episodes_Felony8 
                              0.8345479                               1.0474303 
          Prior_Arrest_Episodes_Felony9          Prior_Arrest_Episodes_Violent1 
                              1.1477367                               1.0729050 
         Prior_Arrest_Episodes_Violent2  Prior_Arrest_Episodes_Violent3 or more 
                              1.1154991                               1.1818388 
        Prior_Conviction_Episodes_Prop1         Prior_Conviction_Episodes_Prop2 
                              1.1559020                               1.2590059 
Prior_Conviction_Episodes_Prop3 or more                   Percent_Days_Employed 
                              1.5390061                               0.3247958 
           Supervision_Risk_Score_First 
                              1.0225361

Code

# INTERPRETATION NOTE:

# Coefficients show log-odds relationship

# Odds ratios > 1 indicate increased risk

# Odds ratios < 1 indicate decreased risk

PART 5: GENERATING PREDICTIONS

Code

# Generate predicted probabilities on test set

test_data <- test_data %>%
  mutate(predicted_prob = predict(logit_model, newdata = test_data, type = "response"))

# Examine distribution of predicted probabilities

summary(test_data$predicted_prob)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.07046 0.45871 0.63226 0.60528 0.76880 0.96374

Code

# Visualize predicted probabilities by actual outcome

ggplot(test_data,
       aes(x = predicted_prob,
           fill = as.factor(recidivism))
       ) +
  geom_histogram(bins = 50, alpha = 0.7, position = "identity"
                 ) +
  scale_fill_manual(values = c("steelblue", "coral"),
                    labels = c("No Recidivism", "Recidivism")
                    ) +
  labs(title = "Distribution of Predicted Probabilities",
       x = "Predicted Probability of Recidivism",
       y = "Count",
       fill = "Actual Outcome")

PART 6: MODEL EVALUATION - OVERALL PERFORMANCE

6.1: Confusion Matrix at Default Threshold (0.5)

Code

test_data <- test_data %>%
  mutate(predicted_class_50 = ifelse(predicted_prob > 0.5, 1, 0))

# Create confusion matrix

cm_50 <- confusionMatrix(as.factor(test_data$predicted_class_50),
  as.factor(test_data$recidivism), positive = "1") # "1" is the positive class (recidivism)

print(cm_50)

Confusion Matrix and Statistics

          Reference
Prediction    0    1
         0 1294  719
         1 1334 3204
                                          
               Accuracy : 0.6866          
                 95% CI : (0.6752, 0.6978)
    No Information Rate : 0.5988          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.3215          
                                          
 Mcnemar's Test P-Value : < 2.2e-16       
                                          
            Sensitivity : 0.8167          
            Specificity : 0.4924          
         Pos Pred Value : 0.7060          
         Neg Pred Value : 0.6428          
             Prevalence : 0.5988          
         Detection Rate : 0.4891          
   Detection Prevalence : 0.6927          
      Balanced Accuracy : 0.6546          
                                          
       'Positive' Class : 1

Code

# Extract key metrics

cat("nKey Metrics at Threshold = 0.5:n")

nKey Metrics at Threshold = 0.5:n

Code

cat("Sensitivity (Recall):", round(cm_50$byClass["Sensitivity"], 3),
    "- Proportion of actual recidivists correctly identifiedn")

Sensitivity (Recall): 0.817 - Proportion of actual recidivists correctly identifiedn

Code

cat("Specificity:", round(cm_50$byClass["Specificity"], 3), "- Proportion of non-recidivists correctly identifiedn")

Specificity: 0.492 - Proportion of non-recidivists correctly identifiedn

Code

cat("Precision (PPV):", round(cm_50$byClass["Precision"], 3),
    "- Proportion of predicted recidivists who actually recidivatedn")

Precision (PPV): 0.706 - Proportion of predicted recidivists who actually recidivatedn

Code

cat("Accuracy:", round(cm_50$overall["Accuracy"], 3), "nn")

Accuracy: 0.687 nn

6.2: ROC Curve and AUC

Code

# Generate ROC curve

roc_obj <- roc(response = test_data$recidivism,
  predictor = test_data$predicted_prob )

# Plot ROC curve

ggroc(roc_obj,
      color = "steelblue",
      size = 1.5
      ) +
  geom_abline(slope = 1,
              intercept = 1,
              linetype = "dashed",
              color = "gray50"
              ) +
  labs(title = "ROC Curve: Recidivism Prediction Model",
       subtitle = paste0("AUC = ", round(auc(roc_obj), 3)),
       x = "1 - Specificity (False Positive Rate)",
       y = "Sensitivity (True Positive Rate)"
       ) +
  annotate("text",
           x = 0.5,
           y = 0.3,
           label = "Random Classifiern(AUC = 0.5)",
           color = "gray50",
           size = 3
           ) +
  coord_fixed()

Code

# Print AUC

cat("Area Under the Curve (AUC):", round(auc(roc_obj), 3), "n")

Area Under the Curve (AUC): 0.732 n

Code

cat("Interpretation: AUC of", round(auc(roc_obj), 2), "indicates", ifelse(auc(roc_obj) > 0.8, "good", "acceptable"), "discriminationnn")

Interpretation: AUC of 0.73 indicates acceptable discriminationnn

6.3: Performance Across Multiple Thresholds

Code

# Test three different thresholds

thresholds_to_test <- c(0.3, 0.5, 0.7)

threshold_results <- map_df(thresholds_to_test, function(thresh) { # Generate predictions at this threshold
  preds <- ifelse(test_data$predicted_prob > thresh, 1, 0)

# Calculate confusion matrix
  cm <- confusionMatrix(as.factor(preds),
                        as.factor(test_data$recidivism),
                        positive = "1")

# Extract metrics
  data.frame(threshold = thresh,
             accuracy = cm$overall["Accuracy"],
             sensitivity = cm$byClass["Sensitivity"],
             specificity = cm$byClass["Specificity"],
             precision = cm$byClass["Precision"],
             f1_score = cm$byClass["F1"],
             # Calculate false positive and false negative rates
             fpr = 1 - cm$byClass["Specificity"],
             fnr = 1 - cm$byClass["Sensitivity"] ) })

# Display results

print(threshold_results)

             threshold  accuracy sensitivity specificity precision  f1_score
Accuracy...1       0.3 0.6476874   0.9686464   0.1685693 0.6349206 0.7670569
Accuracy...2       0.5 0.6866127   0.8167219   0.4923896 0.7060379 0.7573573
Accuracy...3       0.7 0.6209739   0.5034412   0.7964231 0.7868526 0.6140215
                   fpr        fnr
Accuracy...1 0.8314307 0.03135356
Accuracy...2 0.5076104 0.18327810
Accuracy...3 0.2035769 0.49655876

Code

# Visualize threshold trade-offs

threshold_results %>%
  select(threshold, sensitivity, specificity, precision) %>%
  pivot_longer(cols = c(sensitivity, specificity, precision),
               names_to = "metric",
               values_to = "value") %>%
  ggplot(aes(x = threshold,
             y = value,
             color = metric,
             group = metric)
         ) +
  geom_line(size = 1.2
            ) +
  geom_point(size = 3
             ) +
  scale_color_brewer(palette = "Set1",
                     labels = c("Precision", "Sensitivity", "Specificity")
                     ) +
  scale_x_continuous(breaks = thresholds_to_test
                     ) +
  labs(title = "Performance Metrics Across Thresholds",
       subtitle = "The threshold-performance trade-off",
       x = "Probability Threshold",
       y = "Metric Value",
       color = "Metric"
       ) +
  theme(legend.position = "bottom")

CRITICAL INSIGHT:

As threshold increases:

Fewer people flagged as high-risk (more conservative)
Sensitivity decreases (miss more actual recidivists)
Specificity increases (fewer false accusations)
Precision usually increases (predictions more accurate when made)

PART 7: EQUITY ANALYSIS - GROUP-WISE PERFORMANCE

7.1: Calculate Performance Metrics by Race

Code

# First, calculate overall metrics for comparison

overall_metrics <- test_data %>%
  summarise(Group = "Overall",
            N = n(),
            Base_Rate = mean(recidivism),
            Sensitivity = sum(predicted_class_50 == 1 & recidivism == 1) / sum(recidivism == 1),
            Specificity = sum(predicted_class_50 == 0 & recidivism == 0) / sum(recidivism == 0),
            Precision = sum(predicted_class_50 == 1 & recidivism == 1) / sum(predicted_class_50 == 1),
            FPR = sum(predicted_class_50 == 1 & recidivism == 0) / sum(recidivism == 0),
            FNR = sum(predicted_class_50 == 0 & recidivism == 1) / sum(recidivism == 1))

# Calculate metrics by race

race_metrics <- test_data %>%
  group_by(Race) %>%
  summarise(N = n(),
            Base_Rate = mean(recidivism),
            Sensitivity = sum(predicted_class_50 == 1 & recidivism == 1) / sum(recidivism == 1),
            Specificity = sum(predicted_class_50 == 0 & recidivism == 0) / sum(recidivism == 0),
            Precision = sum(predicted_class_50 == 1 & recidivism == 1) / sum(predicted_class_50 == 1),
            FPR = sum(predicted_class_50 == 1 & recidivism == 0) / sum(recidivism == 0),
            FNR = sum(predicted_class_50 == 0 & recidivism == 1) / sum(recidivism == 1) ) %>%
  rename(Group = Race)

# Combine overall and by-group metrics

equity_analysis <- bind_rows(overall_metrics, race_metrics)

# Display results

print("=" %>% rep(80) %>% paste0(collapse = ""))

[1] "================================================================================"

Code

print("EQUITY ANALYSIS: Model Performance by Race (Threshold = 0.5)")

[1] "EQUITY ANALYSIS: Model Performance by Race (Threshold = 0.5)"

Code

print("=" %>% rep(80) %>% paste0(collapse = ""))

[1] "================================================================================"

Code

print(equity_analysis, width = Inf)

# A tibble: 3 × 8
  Group       N Base_Rate Sensitivity Specificity Precision   FPR   FNR
  <chr>   <int>     <dbl>       <dbl>       <dbl>     <dbl> <dbl> <dbl>
1 Overall  6551     0.599       0.817       0.492     0.706 0.508 0.183
2 BLACK    3931     0.598       0.846       0.438     0.691 0.562 0.154
3 WHITE    2620     0.600       0.773       0.575     0.732 0.425 0.227

Code

print("=" %>% rep(80) %>% paste0(collapse = ""))

[1] "================================================================================"

7.2: Visualize Disparate Impact

Code

# Create bar chart comparing key metrics across groups

race_metrics %>%
  select(Group, FPR, FNR, Sensitivity, Specificity) %>%
  pivot_longer(cols = c(FPR, FNR, Sensitivity, Specificity),
               names_to = "Metric", values_to = "Rate") %>%
  ggplot(aes(x = Group, y = Rate, fill = Metric)
         ) +
  geom_col(position = "dodge",
           alpha = 0.8
           ) +
  scale_fill_brewer(palette = "Set2"
                    ) +
  labs(title = "Model Performance Disparities Across Racial Groups",
       subtitle = "Using threshold = 0.5",
       x = "Race",
       y = "Rate",
       fill = "Metric",
       caption = "FPR = False Positive Rate, FNR = False Negative Rate"
       ) +
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1))

7.3: Testing Threshold Adjustments for Equity

Code

# Can we achieve more equitable outcomes by adjusting thresholds?

# Function to calculate key rates at a given threshold

calc_rates_by_threshold <- function(data, threshold) {
  data %>%
    mutate(pred = ifelse(predicted_prob > threshold, 1, 0)) %>%
    summarise( threshold = threshold,
               FPR = sum(pred == 1 & recidivism == 0) / sum(recidivism == 0),
               FNR = sum(pred == 0 & recidivism == 1) / sum(recidivism == 1),
               Sensitivity = sum(pred == 1 & recidivism == 1) / sum(recidivism == 1),
               Specificity = sum(pred == 0 & recidivism == 0) / sum(recidivism == 0) ) }

# Test range of thresholds for each racial group

threshold_range <- seq(0.3, 0.7, by = 0.05)

threshold_by_race <- test_data %>%
  nest_by(Race) %>%
  reframe(map_df(threshold_range,
                 ~calc_rates_by_threshold(data, .x))) %>%
  ungroup()

# Visualize FPR across thresholds by race

ggplot(threshold_by_race,
       aes(x = threshold, y = FPR,color = Race)
       ) +
  geom_line(size = 1.2
            ) +
  geom_point(size = 2
            ) +
  geom_vline(xintercept = 0.5,
             linetype = "dashed",
             alpha = 0.5
             ) +
  labs(title = "False Positive Rate by Threshold and Race",
       subtitle = "Could different thresholds equalize FPR?",
       x = "Probability Threshold",
       y = "False Positive Rate",
       color = "Race")

Code

# Visualize FNR across thresholds by race

ggplot(threshold_by_race,
       aes(x = threshold, y = FNR, color = Race)
       ) +
  geom_line(size = 1.2
            ) +
  geom_point(size = 2
             ) +
  geom_vline(xintercept = 0.5,
             linetype = "dashed",
             alpha = 0.5
             ) +
  labs(title = "False Negative Rate by Threshold and Race",
       subtitle = "The trade-off: equalizing FPR may worsen FNR disparities",
       x = "Probability Threshold",
       y = "False Negative Rate",
       color = "Race")

PART 8: SYNTHESIS AND RECOMMENDATIONS

Code

# Fill in the sections below.

cat("\n")

Code

cat("\nKEY FINDINGS & POLICY IMPLICATIONSn")


KEY FINDINGS & POLICY IMPLICATIONSn

Code

cat("\n\n")

Code

cat("1. MODEL PERFORMANCE:\n")

1. MODEL PERFORMANCE:

Code

cat(" - AUC:", round(auc(roc_obj), 3), "\n")

 - AUC: 0.732

Code

cat(" - Overall accuracy at threshold = 0.5:", round(cm_50$overall["Accuracy"], 3), "\n")

 - Overall accuracy at threshold = 0.5: 0.687

Code

cat(" - Interpretation: [Your assessment here]\n\n")

 - Interpretation: [Your assessment here]

Code

cat("2. THRESHOLD TRADE-OFFS:\n")

2. THRESHOLD TRADE-OFFS:

Code

cat(" - Lower threshold (0.3): Higher sensitivity, more false positives\n")

 - Lower threshold (0.3): Higher sensitivity, more false positives

Code

cat(" - Higher threshold (0.7): Higher specificity, more false negatives\n")

 - Higher threshold (0.7): Higher specificity, more false negatives

Code

cat(" - Decision depends on: relative costs of each error type\n\n")

 - Decision depends on: relative costs of each error type

Code

cat("3. EQUITY CONCERNS:\n")

3. EQUITY CONCERNS:

Code

cat(" - [Identify which groups face higher FPR or FNR]\n")

 - [Identify which groups face higher FPR or FNR]

Code

cat(" - [Discuss potential causes: differential base rates, proxy variables]\n")

 - [Discuss potential causes: differential base rates, proxy variables]

Code

cat(" - [Consider impossibility of perfect fairness when base rates differ]\n\n")

 - [Consider impossibility of perfect fairness when base rates differ]

Code

cat("4. DEPLOYMENT RECOMMENDATION:\n")

4. DEPLOYMENT RECOMMENDATION:

Code

cat(" Should this model be used for parole decisions?\n")

 Should this model be used for parole decisions?

Code

cat(" - Technical readiness: [Your assessment]\n")

 - Technical readiness: [Your assessment]

Code

cat(" - Ethical considerations: [Your analysis]\n")

 - Ethical considerations: [Your analysis]

Code

cat(" - Required safeguards: [Your recommendations]\n")

 - Required safeguards: [Your recommendations]

Code

cat(" - Alternatives: [Your suggestions]\n\n")

 - Alternatives: [Your suggestions]

ADDITIONAL ANALYSES (OPTIONAL EXTENSIONS)

Extension 1: Calibration Analysis

Code

# Check if predicted probabilities match actual outcomes

test_data %>%
  mutate(prob_bin = cut(predicted_prob,
                        breaks = seq(0, 1, by = 0.1))
         ) %>%
  group_by(prob_bin) %>%
  summarise(mean_predicted = mean(predicted_prob),
            mean_observed = mean(recidivism),
            n = n() ) %>%
  filter(n > 10) %>%
  ggplot(aes(x = mean_predicted, y = mean_observed, size = n)
         ) +
  geom_point(alpha = 0.6, color = "steelblue"
             ) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red"
              ) +
  labs(title = "Calibration Plot",
       subtitle = "Are predicted probabilities well-calibrated?",
       x = "Mean Predicted Probability",
       y = "Mean Observed Recidivism Rate",
       size = "N in bin"
       ) +
  coord_fixed()

Extension 2: Feature Importance

Code

# Which variables matter most?

coef_df <- summary(logit_model)$coefficients %>%
  as.data.frame() %>% rownames_to_column("variable") %>%
  mutate(odds_ratio = exp(Estimate)) %>%
  filter(variable != "(Intercept)") %>%
  arrange(desc(abs(Estimate)))

ggplot(coef_df,
       aes(x = reorder(variable, Estimate),
           y = odds_ratio)
       ) +
  geom_point(size = 3, color = "steelblue"
             ) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red"
             ) +
  coord_flip() +
  scale_y_continuous(trans = "log10") +
  labs(title = "Feature Importance (Odds Ratios)",
       subtitle = "OR > 1 increases recidivism risk; OR < 1 decreases risk",
       x = "Variable",
       y = "Odds Ratio (log scale)")