Submit an Rmd notebook as well as the PDF generated from the notebook on Gradescope. Pages on Gradescope should be assigned to questions.
In this part of the project, you will analyze the COMPAS recidivism risk scores released by ProPublica. COMPAS is a commercial tool that assigns each criminal defendant a risk score (a decile score from 1 to 10) that is meant to predict whether the defendant will re-offend. Scores like these are used in the United States to inform decisions about bail, sentencing, and parole.
The dataset contains records for over 10,000 defendants from Broward County, Florida, including each defendant’s COMPAS score and whether they actually went on to re-offend (“recidivate”) within two years. This lets us ask two intertwined questions:
These two questions cannot be answered in isolation. A central lesson of this part of the project is that there are several reasonable, mutually incompatible definitions of “fair,” and a score can satisfy one while violating another. You will compute the relevant quantities yourself and form your own view.
The risk-prediction system you analyze here is a real one that has been used to make decisions about real people. Keep that in mind as you interpret your numbers.
The dataset is hosted by ProPublica. You can read it directly from the web:
url <- "https://github.com/propublica/compas-analysis/raw/master/compas-scores-two-years.csv"
raw <- read.csv(url, stringsAsFactors = FALSE)Following ProPublica’s analysis, we keep only rows for which the COMPAS screening happened within 30 days of the arrest, for which a recidivism outcome is recorded, and for which the charge degree is a felony or misdemeanour (not an ordinary traffic offence). Apply this filter exactly as below:
library(dplyr)
compas <- raw %>%
filter(days_b_screening_arrest <= 30,
days_b_screening_arrest >= -30,
is_recid != -1,
c_charge_degree != "O",
score_text != "N/A")This leaves 6172 defendants. The first few rows of the columns you will use look like this:
sex age race priors_count c_charge_degree decile_score score_text two_year_recid
1 Male 69 Other 0 F 1 Low 0
2 Male 34 African-American 0 F 3 Low 1
3 Male 24 African-American 4 F 4 Low 1
4 Male 44 Other 0 M 1 Low 0
5 Male 41 Caucasian 14 F 6 Medium 1
6 Male 43 Other 3 F 4 Low 0The columns you will need:
decile_score — the COMPAS risk score, an integer from 1 (lowest risk) to 10 (highest risk).score_text — COMPAS’s own coarse label: Low, Medium, or High.two_year_recid — the ground truth: 1 if the defendant was arrested for a new offence within two years, 0 otherwise. (Overall, about 45.5% of defendants in the filtered data recidivated.)race — self-reported race. The two largest groups are African-American (3175 defendants) and Caucasian (2103); the remaining groups (Hispanic, Other, Asian, Native American) are much smaller. For Parts 1–3 you will compare the two largest groups.age, sex, priors_count (number of prior offences), and c_charge_degree (F = felony, M = misdemeanour) are candidate predictors.In this part of the project, we treat a defendant as “predicted to re-offend” when their decile score exceeds a threshold that you choose — this is exactly the threshold-on-a-score situation from the ICU project, except the score now comes from COMPAS rather than from a model you fit yourself.
Your code should be general enough that if the dataset were changed (to one with the same column names), the code would still run. In particular, do not hard-code the number 6172 or the names of the racial groups where a group_by or a computed value would do.
You will be graded on the correctness of your code as well as on the professionalism of the presentation of your report. The text should be clear; figures should visualize the data well and have appropriate labels and legends. Use ggplot for plots and dplyr/tidyverse for data manipulation.
When you are asked to compute something, use code (rather than computing things outside of R and typing in the answer), and include the code in the report.
Restrict attention to defendants whose race is African-American or Caucasian. For each of these two groups, plot a histogram of decile_score (the scores 1–10). Put the two histograms side by side, or overlay them, so that the shapes can be compared directly.
Describe in words how the two distributions differ. (You are not asked to explain why they differ here — just to read the plot.)
We will say a defendant is classified as high-risk when their decile_score is greater than 4 (so that Low scores are “low-risk” and Medium/High scores are “high-risk”). Treat two_year_recid as the ground truth.
Write a function that, given a subset of the data, returns:
Report these three numbers separately for African-American and for Caucasian defendants. Comment on what you find: in which direction does each error rate differ between the groups, and why might that be troubling if the score is used to set bail?
The choice of threshold 4 in Part 2 was arbitrary. Sweep the threshold over the values 0.5, 1.5, 2.5, ..., 9.5 (i.e. every place a defendant could move from “low-risk” to “high-risk”). For each threshold, compute the False Positive Rate and the True Positive Rate separately for each of the two groups.
Produce a plot of True Positive Rate vs. False Positive Rate (an ROC curve) with one curve per group. Does COMPAS appear to rank defendants equally well within each group? Is the ranking quality (the shape of each ROC curve) the same thing as the fairness you examined in Part 2? Explain in two or three sentences.
Now build your own score instead of using COMPAS’s. First, split the filtered data into training, validation, and test sets (roughly 60/20/20) by sampling row indices.
Using the training set, fit a logistic regression that predicts two_year_recid from age and priors_count:
model <- glm(two_year_recid ~ age + priors_count, family = binomial, data = train)Report the fitted coefficients and interpret the sign and rough magnitude of each one in plain language (what does the age coefficient say about who is predicted to re-offend?). Then compare the model’s correct-classification rate on the training set to its rate on the validation set, using a threshold of 0.5 on the predicted probability. Is there evidence of over-fitting?
Try adding at least four other variables (individually or in combination) to the model from Part 4 — for example sex, c_charge_degree, or interactions. Use the validation set (not the test set) to decide whether each addition helps. Report what you tried and which model you would choose, and justify the choice using the validation results. Explain why it would be a mistake to make this decision using the test set.
In Part 2 you (probably) found that a single threshold produces different false positive rates for the two groups. One notion of fairness asks that the false positive rates be equal across groups.
Using your chosen model from Part 5, find group-specific thresholds — one for African-American defendants and one for Caucasian defendants — such that the two groups have (approximately) equal false positive rates, while making the overall correct-classification rate as high as you can subject to that constraint. Describe your search procedure and report the thresholds you chose and the resulting error rates. Use the validation set to choose the thresholds, and report the final numbers on the test set.
In two or three sentences, state the cost of this fairness criterion: what do you give up, and for whom, in order to equalize false positive rates?
Choose two variables in the dataset (other than the pair you were told to use in Part 4) and make a single figure that shows how they relate to each other and to the COMPAS score or to recidivism. Explain what the figure shows and why it is relevant to the fairness questions above.
So far your only model has been logistic regression, which is just a neural network with no hidden layer. In this part you will fit a network with one hidden layer of 10 units using the nnet package, and output its weights.
For the neural network, throw in the kitchen sink: use every predictor that is known at the time of the COMPAS screening. There is exactly one rule — you must not use any column that records a future arrest, because those columns leak the answer. They are recorded after the screening and encode the very outcome you are trying to predict, so a model that uses them would look brilliant in your report and be useless in practice. The future-arrest columns to exclude are the target two_year_recid itself, is_recid, is_violent_recid, violent_recid, every r_* and vr_* column (the new charge), and start, end, event. (It also makes sense to drop COMPAS’s own score columns — decile_score, score_text, the v_* scores — and free-text/identifier/date columns, since the point is to predict the outcome from raw inputs.)
A reasonable kitchen-sink predictor set is:
predictors <- c("sex", "age", "age_cat", "race",
"juv_fel_count", "juv_misd_count", "juv_other_count",
"priors_count", "c_charge_degree")
form <- as.formula(paste("two_year_recid ~", paste(predictors, collapse = " + ")))Fit the network on the training set:
library(nnet)
set.seed(1)
nn <- nnet(form, data = train, size = 10, maxit = 300, decay = 0.01, trace = FALSE)size = 10 is the number of hidden units; decay is a small amount of weight regularization (which keeps the optimization well-behaved); and because nnet is sensitive to the starting point, you must set.seed so your report is reproducible.
Output the weights with
summary(nn) # weights, grouped by their destination unit
nn$wts # the same weights as a single numeric vector
nn$coefnames # which predictor each input i1, i2, ... corresponds tosummary(nn) reports the architecture and every weight. Multi-level factors are expanded into 0/1 dummy columns (so race, with six levels, contributes five inputs, and so on); with this predictor set the predictors expand to 14 inputs, and with 10 hidden units and 1 output the network is 14-10-1 and has 161 weights. Your output should begin like this:
a 14-10-1 network with 161 weights
options were - decay=0.01
b->h1 i1->h1 i2->h1 i3->h1 i4->h1 i5->h1 i6->h1 i7->h1 i8->h1 i9->h1
-1.71 2.22 -0.40 -0.01 -2.03 -0.15 0.40 2.15 -0.01 7.09
i10->h1 i11->h1 i12->h1 i13->h1 i14->h1
-0.03 0.59 5.19 0.13 2.68
b->h2 i1->h2 ...
...
b->o h1->o h2->o ... h10->oHere iN->hj is the weight from input N to hidden unit j, b->hj is the bias into hidden unit j, and hj->o is the weight from hidden unit j to the output. Use nn$coefnames to state in your report which predictor each input number (i1…i14) corresponds to.
A single train/validation/test split gives a noisy estimate of how good each model is — you might just have gotten a lucky split. To compare the neural network with logistic regression fairly, resample 10 times: repeat the whole split-and-fit procedure 10 times with a different random partition each time, and average the results.
Concretely, write a function that does one resample — it draws a fresh 60/20/20 split, fits both the logistic-regression model and the 10-unit neural network on the training set, and returns the correct-classification rate of each model on the validation set and on the test set. To make this a fair, apples-to-apples comparison of the two model families, give both models the same kitchen-sink predictors (form) from Part 8. Then call it 10 times and average:
acc <- function(prob, truth) mean((prob > 0.5) == (truth == 1))
compare_once <- function() {
n <- nrow(compas)
idx <- sample(n)
tr <- idx[1:(0.6*n)]; va <- idx[(0.6*n+1):(0.8*n)]; te <- idx[(0.8*n+1):n]
train <- compas[tr, ]; val <- compas[va, ]; test <- compas[te, ]
g <- glm(form, family = binomial, data = train)
nn <- nnet(form, data = train, size = 10, maxit = 300, decay = 0.01, trace = FALSE)
c(glm_val = acc(predict(g, val, type = "response"), val$two_year_recid),
nn_val = acc(as.vector(predict(nn, val)), val$two_year_recid),
glm_test = acc(predict(g, test, type = "response"), test$two_year_recid),
nn_test = acc(as.vector(predict(nn, test)), test$two_year_recid))
}
set.seed(1)
results <- replicate(10, compare_once())
rowMeans(results)A typical run gives numbers close to
glm_val nn_val glm_test nn_test
0.672 0.670 0.675 0.673Report your averaged validation and test accuracies for both models. Then answer: is the difference between the neural network and logistic regression larger or smaller than the spread across the 10 resamples? (Look at apply(results, 1, sd) — here it is about 0.01, several times larger than the gap between the two models.) What does that tell you about whether the extra flexibility of the hidden layer is buying you anything on this dataset?
The output of summary(nn) lists, for every hidden unit, the weights coming into it from each input (plus a bias), and the single weight carrying that unit’s activation to the output. This lets you give each unit a rough “job description.”
For each of the 10 hidden units, read off its incoming weights (use nn$coefnames to translate i1…i14 into predictor names) and report which input(s) it responds to most strongly (largest absolute weights) and the sign of its connection to the output (hj->o). Then, in one sentence per unit, interpret what that unit appears to detect. For example, a unit with a large positive weight on priors_count and a positive h->o weight is acting as a “many-priors ⇒ higher predicted risk” detector, whereas a unit whose weights are all near zero is doing little.
Finally, comment on the limits of this exercise: hidden units in a network are not guaranteed to be individually meaningful (several can be redundant, or split one concept between them), which is part of why a neural network is harder to explain to a defendant than the logistic regression of Part 4. Tie this back to the fairness discussion: why does interpretability matter when a model is used to make decisions about people?
Finally, fit a model that is at the opposite end of the interpretability spectrum: a decision tree, using the rpart package. As with the neural network, throw in the kitchen sink — use the same form (every predictor known at screening, no future-arrest columns). A tree handles many predictors gracefully: it will simply ignore the ones it does not find useful, so you can see which features it chooses to split on. Fit it on the training set and visualize it:
library(rpart)
set.seed(1)
tree <- rpart(update(form, factor(two_year_recid) ~ .),
data = train, method = "class")
plot(tree, margin = 0.1) # draw the branches
text(tree, use.n = TRUE) # label the splits and leaves(If you have the rpart.plot package, rpart.plot(tree) produces a nicer figure — either is acceptable.)
Include the tree in your report and describe it: which variable does the tree split on first, and what does that say about which feature is most predictive of recidivism? Follow one path from the root to a leaf and state, in plain language, the rule it represents (e.g. “a defendant with more than X prior offences and age below Y is predicted to re-offend”).
Compare the three model families you have now fit — logistic regression, the neural network, and the decision tree — in one paragraph: which is easiest to explain to the person being scored, which performed best in Part 9, and how you would weigh accuracy against explainability if this model were actually used to set bail.