I should find the optimal threshold to minimize both the false positive rate and false negative rate. An equal weight between these two rates should be assumed. I write the following code:

data=read.csv( url("https://raw.githubusercontent.com/propublica/compas-analysis/master/compas-scores-two-years.csv"), sep=",")
pred=prediction(data$decile_score/10, data$two_year_recid)
perf=performance(pred, measure="fnr",x.measure="fpr")

opt.cut = function(perf, pred)
    cut.ind = mapply(FUN=function(x, y, p){
        d = (x - 0)^2 + (y-1)^2
        ind = which(d == min(d))
        c(False_negative_rate = 1-y[[ind]], False_positive_rate = x[[ind]], 
            cutoff = p[[ind]])
    }, [email protected], [email protected], pred@cutoffs)

print(opt.cut(perf, pred))

It throws out this result:

False_negative_rate    0
False_positive_rate    0
cutoff               Inf

However, I think there is something wrong with my code.

  • $\begingroup$ I like to use pacman, and p_load, because if it isn't installed it gets installed. $\endgroup$ Oct 26, 2022 at 13:30

1 Answer 1


Here are the libraries I used:


I had to run the pull several times to get the csv to download.


Your data has 7214 rows and 53 columns, but you are only looking at 2. They are decile_score and two_year_recid.

Here is how I check the type, so I know they can plot well together.

data %>% select(decile_score, two_year_recid) %>% str()

This is the result:

'data.frame':   7214 obs. of  2 variables:
 $ decile_score  : int  1 3 4 8 1 1 6 4 1 3 ...
 $ two_year_recid: int  0 1 1 0 0 0 1 0 0 1 ...

I looks like two numbers. Looks like the bottom is binomial, and the top is integers 1 to 10.

Here is a plot of the elements of prediction, to start to make sure that whatever "prediction" gives isn't junk.

ggplot(data, aes(x=decile_score/10, y=two_year_recid)) + 
  geom_point() +
  geom_smooth(method="glm", formula = 'y ~ x ')

Here is the plot it yields: enter image description here

Yes, nearly everything is wrong with the plot, but it still helps us ask questions.

In a very general sense, as decile score goes up, rate of recidivism goes up. As decile score goes down, rate of recidivism goes down. The place where recidivism chance is 50/50 looks to be aproximately the 50th percentile, or maybe a little above it.

Instead of using your pred, I'm going to use a proper linear model, and I'm going to turn recidivism into a proper data type. It is currently a float, but it should be a boolean, and the fit and prediction functions can respond to that.

df <- data.frame(x=data$decile_score/10,
                 y=data$two_year_recid %>% as.factor)

est_glm <- glm(y~x, data=df, family = binomial()) 

This yields: enter image description here

The formula being fit in logistic regression is this:
enter image description here

The place of 50/50 odds occurs when:
$\beta_0 + \beta_1 \cdot x = 0$

solving that for $x$ gives:
$ -1.408958 + 2.650818 \cdot x = 0 $
$ 2.650818 \cdot x = 1.408958 $
$ x = \frac {1.408958} {2.650818} $

which evaluates to 53.15183%

This suggests that whatever decile is, it slightly under-predicts recidivism.

Some of the diagnostic plots from the glm fit can be interesting.

This suggests that the linear model in the exponent might not be the best. (I would test a cubic, and compare AIC values) enter image description here

This one says there may be some outliers, a few points that corrupt the whole analysis.
enter image description here

Here is a violin plot that might speak to the nature of the population. enter image description here

Why are there so many more non-recidivists at lower decile than high decile when the population of recidivists is nearly constant across deciles? There is a great reason for it, but without reading the paper, or understanding the sampling this can look like a surprising result.

So now for the tnr/tpr parts.

Here is the code I used:

cut_list <- seq(from=min(Preds), to=max(Preds), length.out=101)
store <- matrix(NA, nrow=length(cut_list), ncol=3) %>% as.data.frame()
names(store) <- c("cutoff", "tnr", "tpr")

for(i in 1:length(cut_list)){

  preds2 <- ifelse(Preds > cut_list[i],1,0) %>% as.vector() %>% as.factor()
  tp_num <- length(which(preds2 ==1 & df$y == 1))
  tp_den <- length(which(preds2 ==1 & df$y == 1)) + length(which(preds2 ==0 & df$y == 1))
  tpr <- tp_num/tp_den
  tn_num <- length(which(preds2 ==0 & df$y == 0))
  tn_den <- length(which(preds2 ==0 & df$y == 0)) + length(which(preds2 ==1 & df$y == 0))
  tnr <- tn_num/tn_den
  store[i,1] <- cut_list[i]
  store[i,2] <- tnr
  store[i,3] <- tpr

I first plot each like this:

store2 <- melt(store, id.vars=1)

ggplot(store2, aes(x=cutoff, y=value, color=variable)) + 
  geom_point() + geom_path()

which gives this: enter image description here

There are many ways to consider these values. Optimal means you have a measure of best, but I don't see one clearly shown here. If you consider them equal, you can get something like this.

Here is the code:

store3 <- store %>% mutate(both=tnr+tpr)

ggplot(store3, aes(x=cutoff, y=both)) + 
  geom_point() + geom_path() +
  geom_vline(xintercept = 0.415, color="Red") +
  geom_text(data = NULL, x = 0.415-0.025, y = 1.05, label = "x=0.415") +
  geom_vline(xintercept = 0.478, color="Green") +
  geom_text(data = NULL, x = 0.478+0.025, y = 1.05, label = "x=0.478") +
  geom_vline(xintercept = 0.5*0.478+0.5*0.415, color="Black") 

Here is the plot:

enter image description here

While you can get the same (big-chunk) result between 41.5% and 47.8%, you are likely to get the best result around 44.65%, which suggests that a likelihood above decile score above 44.65% is more likely than not to have recidivism within 2 years.

Like I said, I think the linear model has problems, and the data has a few points doing more of the work than they should. I don't trust someone else's model if I can't understand it, because that's a really fast way to make a big error, and I don't know where decile_score came from in detail.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.