Obtaining predicted values (Y=1 or 0) from a logistic regression model fit 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.
 A: 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. 

# 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"))

