Let's say that I have an object of class glm (corresponding to a logistic regression model) and I'd like to turn the predicted probabilities given by predict.glm using the argument type="response" into binary responses, i.e. $Y=1$ or $Y=0$. What's the quickest & most canonical way to do this in R?

While, again, I'm aware of predict.glm, I don't know where exactly the cutoff value $P(Y_i=1|\hat X_{i})$ lives -- and I guess this is my main stumbling block here.


Once you have the predicted probabilities, it is up to you what threshold you would like to use. You may choose the threshold to optimize sensitivity, specificity or whatever measure it most important in the context of the application (some additional info would be helpful here for a more specific answer). You may want to look at ROC curves and other measures related to optimal classification.

Edit: To clarify this answer somewhat I'm going to give an example. The real answer is that the optimal cutoff depends on what properties of the classifier are important in the context of the application. Let $Y_{i}$ be the true value for observation $i$, and $\hat{Y}_{i}$ be the predicted class. Some common measures of performance are

(1) Sensitivity: $P(\hat{Y}_i=1 | Y_i=1)$ - the proportion of '1's that are correctly identified as so.

(2) Specificity: $P(\hat{Y}_i=0 | Y_i=0)$ - the proportion of '0's that are correctly identified as so

(3) (Correct) Classification Rate: $P(Y_i = \hat{Y}_i)$ - the proportion of predictions that were correct.

(1) is also called True Positive Rate, (2) is also called True Negative Rate.

For example, if your classifier were aiming to evaluate a diagnostic test for a serious disease that has a relatively safe cure, the sensitivity is far more important that the specificity. In another case, if the disease were relatively minor and the treatment were risky, specificity would be more important to control. For general classification problems, it is considered "good" to jointly optimize the sensitivity and specification - for example, you may use the classifier that minimizes their Euclidean distance from the point $(1,1)$:

$$ \delta = \sqrt{ [P(Y_i=1 | \hat{Y}_i=1)-1]^2 + [P(Y_i=0 | \hat{Y}_i=0)-1]^2 }$$

$\delta$ could be weighted or modified in another way to reflect a more reasonable measure of distance from $(1,1)$ in the context of the application - euclidean distance from (1,1) was chosen here arbitrarily for illustrative purposes. In any case, all of these four measures could be most appropriate, depending on the application.

Below is a simulated example using prediction from a logistic regression model to classify. The cutoff is varied to see what cutoff gives the "best" classifier under each of these three measures. In this example the data comes from a logistic regression model with three predictors (see R code below plot). As you can see from this example, the "optimal" cutoff depends on which of these measures is most important - this is entirely application dependent.

Edit 2: $P(Y_i = 1 | \hat{Y}_i = 1)$ and $P(Y_i = 0 | \hat{Y}_i = 0)$, the Positive Predictive Value and Negative Predictive Value (note these are NOT the same as sensitivity and specificity) may also be useful measures of performance.

enter image description here

# data y simulated from a logistic regression model 
# with with three predictors, n=10000
x = matrix(rnorm(30000),10000,3)
lp = 0 + x[,1] - 1.42*x[2] + .67*x[,3] + 1.1*x[,1]*x[,2] - 1.5*x[,1]*x[,3] +2.2*x[,2]*x[,3] + x[,1]*x[,2]*x[,3]
p = 1/(1+exp(-lp))
y = runif(10000)<p

# fit a logistic regression model
mod = glm(y~x[,1]*x[,2]*x[,3],family="binomial")

# using a cutoff of cut, calculate sensitivity, specificity, and classification rate
perf = function(cut, mod, y)
   yhat = (mod$fit>cut)
   w = which(y==1)
   sensitivity = mean( yhat[w] == 1 ) 
   specificity = mean( yhat[-w] == 0 ) 
   c.rate = mean( y==yhat ) 
   d = cbind(sensitivity,specificity)-c(1,1)
   d = sqrt( d[1]^2 + d[2]^2 ) 
   out = t(as.matrix(c(sensitivity, specificity, c.rate,d)))
   colnames(out) = c("sensitivity", "specificity", "c.rate", "distance")

s = seq(.01,.99,length=1000)
OUT = matrix(0,1000,4)
for(i in 1:1000) OUT[i,]=perf(s[i],mod,y)
legend(0,.25,col=c(2,"darkgreen",4,"darkred"),lwd=c(2,2,2,2),c("Sensitivity","Specificity","Classification Rate","Distance"))
  • 2
    $\begingroup$ (+1) Very nice answer. I like the example. Is there a ready interpretation that you know of to motivate the use of the Euclidean distance that you've given? I also think it might be interesting to point out in this context that the ROC curve is essentially obtained by making a post hoc modification of the intercept estimate of the logistic model. $\endgroup$ – cardinal Apr 1 '12 at 1:59
  • $\begingroup$ @Cardinal, I know that thresholds for binary classification are often chosen based on which point on the ROC curve is closest to (1,1) - euclidean distance was arbitrarily the default definition of "distance" in my example $\endgroup$ – Macro Apr 1 '12 at 2:12
  • $\begingroup$ I see. I thought there might be an intuitive interpretation of this quantity in terms of underlying model that I was not seeing. (Maybe there is[?]) $\endgroup$ – cardinal Apr 1 '12 at 2:20
  • $\begingroup$ Maybe there is because I've found the point that the sensitivity and specificity curves intersect is the same place that $\delta$ is minimized... $\endgroup$ – Macro Apr 1 '12 at 2:45
  • $\begingroup$ I don't have an $R$ terminal to use at the moment, but I'm wondering if the effect you're seeing is possibly due in part to the fact that your simulation should be generating a roughly equal number of 0 and 1 responses. If you vary the proportion of ones, does the relationship still hold? $\endgroup$ – cardinal Apr 1 '12 at 3:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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