Deciding threshold for glm logistic regression model in R

Question

I have some data with predictors and a binary target. Eg:

df <- data.frame(a=sort(sample(1:100,30)), b= sort(sample(1:100,30)), 
                 target=c(rep(0,11),rep(1,4),rep(0,4),rep(1,11)))

I trained a logistic regresion model using glm()

model1 <- glm(formula= target ~ a + b, data=df, family=binomial)

Now I'm trying to predict the output (for the example, the same data should suffice)

predict(model1, newdata=df, type="response")

This generates a vector of probability numbers. But I want to predict the actual class. I could use round() on the probablity numbers, but this assumes that anything below 0.5 is class '0', and anything above is class '1'. Is this a correct assumption? Even when the population of each class may not be equal (or close to equal)? Or is there a way to estimate this threshold?

Beeson

there are different criteria, one for example is the point where the sum of sensitivity and specificity is maximal, see for example this question: #23132397 – Doormat

@Doormat Thanks! But it would be certainly incorrect to use the threshold as the population fraction, right? That is, if in the population, 30% cases are '0's, and 70% '1's, a naive estimate would be to use 0.3 as the threshold. But this would not be a logical way to approach this? – Beeson

You can find a great tutorial on the subject here: hopstat.wordpress.com/2014/12/19/… – Franconian

Answer 1

The best threshold (or cutoff) point to be used in glm models is the point which maximises the specificity and the sensitivity. This threshold point might not give the highest prediction in your model, but it wouldn't be biased towards positives or negatives. The ROCR package contain functions that can help you do this. check the performance() function in this package. It is going to get you what you're looking for. Here's a picture of what you are expecting to get:

After finding the cutoff point, I normally write a function myself to find the number of datapoints that has their prediction value above the cutoff, and match it with the group they belong to.

Anguine

could you provide a more specific code that would generate the above graph? Also, how can cutoff values be between 0 and 14 for probabilities that take values between 0 and 1? – Stephenstephenie

I added baseR/ggplot approaches below! – Sinfonietta

Answer 2

The gold standard for determining good model parameters, including "what threshold should I set" for logistic regression, is cross-validation.

The general idea is to hold out one or more parts of your training set and choose the threshold that maximizes the number of correct classifications on this held-out set, but Wikipedia can give you many more details.

Cower

Since we'd be tuning the threshold parameter on the cross validation data, ostensibly, that would require a third held out set for evaluation to report an unbiased expected error? – Beeson

@user2175594, Yes that is correct. Traditionally you would have at least three separate partitions of your data: training, validation, and test (evaluation). However, if you are doing something like k-fold cross-validation, then training and validation are essentially the same set re-partitioned in multiple ways. – Cower

The tuning threshold doesn't have to be determined by accuracy (# correct). It can be a metric like AUC or F1 score. – Stenosis

Answer 3

Tooled around trying to replicate first graph. Given a predictions <- prediction(pred,labels) object, then:

baseR approach

plot(unlist(performance(predictions, "sens")@x.values), unlist(performance(predictions, "sens")@y.values), 
     type="l", lwd=2, ylab="Specificity", xlab="Cutoff")
par(new=TRUE)
plot(unlist(performance(predictions, "spec")@x.values), unlist(performance(predictions, "spec")@y.values), 
     type="l", lwd=2, col='red', ylab="", xlab="")
axis(4, at=seq(0,1,0.2),labels=z)
mtext("Specificity",side=4, padj=-2, col='red')

ggplot2 approach

sens <- data.frame(x=unlist(performance(predictions, "sens")@x.values), 
                   y=unlist(performance(predictions, "sens")@y.values))
spec <- data.frame(x=unlist(performance(predictions, "spec")@x.values), 
                   y=unlist(performance(predictions, "spec")@y.values))

sens %>% ggplot(aes(x,y)) + 
  geom_line() + 
  geom_line(data=spec, aes(x,y,col="red")) +
  scale_y_continuous(sec.axis = sec_axis(~., name = "Specificity")) +
  labs(x='Cutoff', y="Sensitivity") +
  theme(axis.title.y.right = element_text(colour = "red"), legend.position="none") 

Sinfonietta

Answer 4

To get the threshold in the data with the closest sensitivity and specificity values (ie. the cross over in the graphs above) programmatically you can you this code which gets pretty close:

predictions = prediction(PREDS, LABELS)

sens = cbind(unlist(performance(predictions, "sens")@x.values), unlist(performance(predictions, "sens")@y.values))
spec = cbind(unlist(performance(predictions, "spec")@x.values), unlist(performance(predictions, "spec")@y.values))
sens[which.min(apply(sens, 1, function(x) min(colSums(abs(t(spec) - x))))), 1]

Bynum

Answer 5

There are 12 methods implemented in the function PresenceAbsence::optimal.thresholds of PresenceAbsence package.

This is also covered in Freeman, E. A., & Moisen, G. G. (2008). A comparison of the performance of threshold criteria for binary classification in terms of predicted prevalence and kappa. Ecological Modelling, 217(1-2), 48-58.

Instar

Answer 6

predict.glm only has type = c("link", "response", "terms"), and "response" is predicted class probabilities, almost always the one you want. It's standard for logistic regression to use 0.5 as a threshold, as argmax conditional class probability, argmax_y P(Y_train=y | X_train=x) is the Bayes (optimal) classifier. Parsnip's predict() for logistic_reg() does this by default. The conditional probability is not the same as the proportion of your response variable in the input data, which is marginal P(Y_train = y), not dependent on X_train.

However, you don't have to use 0.5 as a threshold, since you can tune the threshold based on any performance metric from the class predictions, such as maximizing F1 score or balancing precision/recall.

Stenosis

Answer 7

You can try the below:

perfspec <- performance(prediction.obj = pred, measure="spec", x.measure="cutoff")

plot(perfspec)

par(new=TRUE)

perfsens <- performance(prediction.obj = pred, measure="sens", x.measure="cutoff")

plot(perfsens)

Donaugh