Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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.

share|improve this question

1 Answer 1

up vote 34 down vote accepted

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.

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 true-positive and true-negative rates (note these are NOT the same as sensitivity and specificity) may also be useful measures of performance. For example, if you were trying to design a diagnostic for when a disease outbreak would occur (in the future), a high true-positive rate would be very desirable, since that would mean that if an outbreak is going to occur, you're very likely to predict that one will occur (and you can implement some intervention). The code could be modified to calculate these instead - I'll leave that to you :)

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")
   return(out)
}

s = seq(.01,.99,length=1000)
OUT = matrix(0,1000,4)
for(i in 1:1000) OUT[i,]=perf(s[i],mod,y)
plot(s,OUT[,1],xlab="Cutoff",ylab="Value",cex.lab=1.5,cex.axis=1.5,ylim=c(0,1),type="l",lwd=2,axes=FALSE,col=2)
axis(1,seq(0,1,length=5),seq(0,1,length=5),cex.lab=1.5)
axis(2,seq(0,1,length=5),seq(0,1,length=5),cex.lab=1.5)
lines(s,OUT[,2],col="darkgreen",lwd=2)
lines(s,OUT[,3],col=4,lwd=2)
lines(s,OUT[,4],col="darkred",lwd=2)
box()
legend(0,.25,col=c(2,"darkgreen",4,"darkred"),lwd=c(2,2,2,2),c("Sensitivity","Specificity","Classification Rate","Distance"))
share|improve this answer
    
+1 for mentioning ROC curve (which is the second question) - care to donate an example? –  Roman Luštrik Mar 28 '12 at 7:07
    
+1, the 'edit' makes this answer truly exceptional, bravo! –  gung Mar 31 '12 at 16:09
    
(+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. –  cardinal Apr 1 '12 at 1:59
    
@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 –  Macro Apr 1 '12 at 2:12
    
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[?]) –  cardinal Apr 1 '12 at 2:20

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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