Erhalten von vorhergesagten Werten (Y = 1 oder 0) aus einer logistischen Regressionsmodellanpassung

50

Nehmen wir an, ich habe ein Klassenobjekt glm(das einem logistischen Regressionsmodell entspricht) und möchte die predict.glmmit dem Argument angegebenen vorhergesagten Wahrscheinlichkeiten type="response"in binäre Antworten umwandeln, dh oder . Was ist der schnellste und kanonischste Weg, dies in R zu tun?Y=1Y=0

Auch predict.glmwenn mir bewusst ist , dass ich nicht genau weiß, wo der - und ich denke, dies ist mein Hauptstolperstein hier.P(Y.ich=1|X^ich)

Tetragrammaton
quelle

Antworten:

77

Sobald Sie die vorhergesagten Wahrscheinlichkeiten haben, liegt es an Ihnen, welchen Schwellenwert Sie verwenden möchten. Sie können den Schwellenwert auswählen, um die Empfindlichkeit, die Spezifität oder das Maß zu optimieren, das im Kontext der Anwendung am wichtigsten ist (einige zusätzliche Informationen wären hier für eine spezifischere Antwort hilfreich). Möglicherweise möchten Sie sich die ROC-Kurven und andere Kennzahlen ansehen, die sich auf die optimale Klassifizierung beziehen.

Bearbeiten: Um diese Antwort etwas zu klären, werde ich ein Beispiel geben. Die eigentliche Antwort lautet, dass der optimale Cutoff davon abhängt, welche Eigenschaften des Klassifikators im Kontext der Anwendung wichtig sind. Sei der wahre Wert für Beobachtung und die vorhergesagte Klasse. Einige gebräuchliche Leistungsmaßstäbe sind i Y iY.ichichY.^ich

(1) Empfindlichkeit: - der Anteil der , die korrekt als identifiziert werden.P(Y.^ich=1|Y.ich=1)

(2) Spezifität: - der Anteil der , die korrekt als identifiziert werdenP(Y.^ich=0|Y.ich=0)

(3) (Richtige) Klassifizierungsrate: - der Anteil der Vorhersagen, die richtig waren.P(Y.ich=Y.^ich)

(1) wird auch als True Positive Rate bezeichnet, (2) wird auch als True Negative Rate bezeichnet.

Wenn Ihr Klassifikator beispielsweise einen diagnostischen Test für eine schwere Krankheit mit relativ sicherer Heilung durchführen wollte, ist die Empfindlichkeit weitaus wichtiger als die Spezifität. In einem anderen Fall, wenn die Krankheit relativ gering und die Behandlung riskant wäre, wäre die Kontrolle der Spezifität wichtiger. Bei allgemeinen Klassifizierungsproblemen wird es als "gut" angesehen, die Empfindlichkeit und Spezifikation gemeinsam zu optimieren. Beispielsweise können Sie den Klassifizierer verwenden, der den euklidischen Abstand von Punkt minimiert :(1,1)

δ=[P(Y.ich=1|Y.^ich=1)-1]2+[P(Y.ich=0|Y.^ich=0)-1]2

δ könnte gewichtet oder auf andere Weise modifiziert werden, um ein vernünftigeres Maß für den Abstand von im Kontext der Anwendung wiederzugeben - der euklidische Abstand von (1,1) wurde hier zu Veranschaulichungszwecken willkürlich gewählt. In jedem Fall könnten alle diese vier Maßnahmen je nach Anwendung am geeignetsten sein.(1,1)

Nachfolgend finden Sie ein simuliertes Beispiel, das die Vorhersage eines logistischen Regressionsmodells zur Klassifizierung verwendet. Der Cutoff wird variiert, um zu sehen, welcher Cutoff den "besten" Klassifikator unter jedem dieser drei Takte ergibt. In diesem Beispiel stammen die Daten aus einem logistischen Regressionsmodell mit drei Prädiktoren (siehe R-Code unten). Wie Sie diesem Beispiel entnehmen können, hängt die "optimale" Abschaltung davon ab, welche dieser Maßnahmen am wichtigsten ist - dies ist vollständig anwendungsabhängig.

P(Y.ich=1|Y.^ich=1)P(Y.ich=0|Y.^ich=0)

Bildbeschreibung hier eingeben

# 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"))
Makro
quelle
2
(+1) Sehr nette Antwort. Ich mag das Beispiel. Gibt es eine fertige Interpretation, die Sie kennen, um die Verwendung der von Ihnen angegebenen euklidischen Distanz zu motivieren? Ich halte es auch für interessant, in diesem Zusammenhang darauf hinzuweisen, dass die ROC-Kurve im Wesentlichen durch eine post-hoc-Modifikation der Abschnittsschätzung des logistischen Modells erhalten wird.
Kardinal
@Cardinal, ich weiß, dass Schwellenwerte für die binäre Klassifizierung häufig anhand des Punkts auf der ROC-Kurve ausgewählt werden, der (1,1) am nächsten liegt. Der euklidische Abstand war in meinem Beispiel willkürlich die Standarddefinition für "Abstand"
Makro
Aha. Ich dachte, es könnte eine intuitive Interpretation dieser Größe in Bezug auf das zugrunde liegende Modell geben, die ich nicht gesehen habe. (Vielleicht gibt es [?])
Kardinal
δ
R