Obtención de valores pronosticados (Y = 1 o 0) a partir de un ajuste del modelo de regresión logística

50

Digamos que tengo un objeto de clase glm(correspondiente a un modelo de regresión logística) y me gustaría convertir las probabilidades predichas dadas predict.glmusando el argumento type="response"en respuestas binarias, es decir, o . ¿Cuál es la forma más rápida y canónica de hacer esto en R?Y = 0Y=1Y=0

Si bien, una vez más, predict.glmsé que no sé exactamente dónde vive el valor de corte , y supongo que este es mi principal obstáculo aquí.P(Yi=1|X^i)

tetragrammaton
fuente

Respuestas:

77

Una vez que tenga las probabilidades predichas, depende de usted qué umbral le gustaría usar. Puede elegir el umbral para optimizar la sensibilidad, especificidad o cualquier medida que sea más importante en el contexto de la aplicación (alguna información adicional sería útil aquí para obtener una respuesta más específica). Es posible que desee ver las curvas ROC y otras medidas relacionadas con la clasificación óptima.

Editar: Para aclarar un poco esta respuesta, voy a dar un ejemplo. La respuesta real es que el límite óptimo depende de qué propiedades del clasificador son importantes en el contexto de la aplicación. Sea el verdadero valor para la observación , y sea ​​la clase predicha. Algunas medidas comunes de rendimiento son i Y iYiiY^i

(1) Sensibilidad: - la proporción de '1' que se identifica correctamente como tal.P(Y^i=1|Yi=1)

(2) Especificidad: - la proporción de '0 que se identifica correctamente como talP(Y^i=0|Yi=0)

(3) Tasa de clasificación (correcta): : la proporción de predicciones correctas.P(Yi=Y^i)

(1) también se llama tasa positiva verdadera, (2) también se llama tasa negativa verdadera.

Por ejemplo, si su clasificador tenía como objetivo evaluar una prueba de diagnóstico para una enfermedad grave que tiene una cura relativamente segura, la sensibilidad es mucho más importante que la especificidad. En otro caso, si la enfermedad fuera relativamente menor y el tratamiento fuera riesgoso, la especificidad sería más importante de controlar. Para problemas generales de clasificación, se considera "bueno" optimizar conjuntamente la sensibilidad y la especificación; por ejemplo, puede usar el clasificador que minimiza su distancia euclidiana desde el punto :(1,1)

δ=[P(Yi=1|Y^i=1)1]2+[P(Yi=0|Y^i=0)1]2

( 1 , 1 )δ podría ser ponderado o modificado de otra manera para reflejar una medida más razonable de la distancia desde en el contexto de la aplicación: la distancia euclidiana desde (1,1) se eligió aquí arbitrariamente con fines ilustrativos. En cualquier caso, estas cuatro medidas podrían ser las más apropiadas, dependiendo de la aplicación.(1,1)

A continuación se muestra un ejemplo simulado que utiliza la predicción de un modelo de regresión logística para clasificar. El límite se varía para ver qué límite da el "mejor" clasificador en cada una de estas tres medidas. En este ejemplo, los datos provienen de un modelo de regresión logística con tres predictores (ver el código R debajo de la gráfica). Como puede ver en este ejemplo, el límite "óptimo" depende de cuál de estas medidas es más importante: esto depende completamente de la aplicación.

Edición 2: y , el valor predictivo positivo y el valor predictivo negativo (tenga en cuenta que estos NO son lo mismo como sensibilidad y especificidad) también pueden ser medidas útiles de rendimiento.P(Yi=1|Y^i=1)P(Yi=0|Y^i=0)

ingrese la descripción de la imagen aquí

# 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"))
Macro
fuente
2
(+1) Muy buena respuesta. Me gusta el ejemplo ¿Hay alguna interpretación que conozcas para motivar el uso de la distancia euclidiana que has dado? También creo que podría ser interesante señalar en este contexto que la curva ROC se obtiene esencialmente haciendo una modificación post hoc de la estimación de intercepción del modelo logístico.
cardenal
@Cardinal, sé que los umbrales para la clasificación binaria a menudo se eligen en función de qué punto de la curva ROC está más cerca de (1,1): la distancia euclidiana era arbitrariamente la definición predeterminada de "distancia" en mi ejemplo
Macro
Veo. Pensé que podría haber una interpretación intuitiva de esta cantidad en términos del modelo subyacente que no estaba viendo. (Tal vez hay [?])
cardenal
Tal vez sí porque descubrí que las curvas de sensibilidad y especificidad se cruzan en el mismo lugar donde se minimiza ...δ
Macro
No tengo una terminal para usar en este momento, pero me pregunto si el efecto que está viendo se debe en parte al hecho de que su simulación debería generar un número aproximadamente igual de respuestas 0 y 1. Si varía la proporción de unos, ¿la relación aún se mantiene? R
cardenal