Ein Beispiel für die Verwendung der Genauigkeit als Ergebnismaß führt zu einer falschen Schlussfolgerung

Ich untersuche verschiedene Leistungsmessungen für Vorhersagemodelle. Es wurde viel über Probleme bei der Verwendung von Genauigkeit geschrieben, anstatt über etwas Kontinuierlicheres, um die Modellleistung zu bewerten. Frank Harrell http://www.fharrell.com/post/class-damage/ liefert ein Beispiel, wenn das Hinzufügen einer informativen Variablen zu einem Modell zu einem Rückgang der Genauigkeit, einer eindeutig eingängigen und falschen Schlussfolgerung führt.

In diesem Fall scheint dies jedoch auf unausgeglichene Klassen zurückzuführen zu sein, und daher kann es nur durch Verwendung einer ausgeglichenen Genauigkeit ((sens + spec) / 2) gelöst werden. Gibt es ein Beispiel, bei dem die Verwendung der Genauigkeit eines ausgewogenen Datensatzes zu eindeutig falschen oder nicht intuitiven Schlussfolgerungen führt?

Bearbeiten

Ich bin auf der Suche nach etwas, bei dem die Genauigkeit sinkt, selbst wenn das Modell deutlich besser ist, oder bei dem die Verwendung der Genauigkeit zu einer falsch positiven Auswahl einiger Merkmale führt. Es ist einfach, falsch negative Beispiele zu erstellen, bei denen die Genauigkeit für zwei Modelle gleich ist, bei denen eines anhand anderer Kriterien eindeutig besser ist.

Ich werde betrügen.

Insbesondere habe ich oft (z. B. hier ) argumentiert, dass sich der statistische Teil der Modellierung und Vorhersage nur auf probabilistische Vorhersagen für Klassenmitgliedschaften erstreckt (oder bei numerischen Vorhersagen Vorhersagedichten angibt). Die Behandlung eine bestimmte Instanz , als ob es zu einer bestimmten Klasse (oder gehörte Punkt Vorhersagen in dem numerischen Fall) ist, Statistiken nicht mehr richtig. Es ist Teil des entscheidungstheoretischen Aspekts.

Und Entscheidungen sollten nicht nur auf der Wahrscheinlichkeitsvorhersage beruhen, sondern auch auf den Kosten von Fehlklassifizierungen und auf einer Vielzahl anderer möglicher Maßnahmen . Selbst wenn Sie beispielsweise nur zwei mögliche Klassen haben, "krank" oder "gesund", können Sie eine Vielzahl möglicher Maßnahmen ergreifen, je nachdem, wie wahrscheinlich es ist, dass ein Patient an der Krankheit leidet und ihn nach Hause schickt, weil er ist mit ziemlicher Sicherheit gesund, ihm zwei Aspirin zu geben, zusätzliche Tests durchzuführen, sofort einen Krankenwagen zu rufen und ihn lebenserhaltend zu behandeln.

Die Beurteilung der Genauigkeit setzt eine solche Entscheidung voraus. Die Genauigkeit als Bewertungsmetrik für die Klassifizierung ist ein Kategoriefehler .

Um Ihre Frage zu beantworten, werde ich den Weg eines solchen Kategoriefehlers beschreiten. Wir werden ein einfaches Szenario mit ausgewogenen Klassen betrachten, in dem eine Klassifizierung ohne Berücksichtigung der Kosten einer Fehlklassifizierung uns tatsächlich stark irreführt.

---

Angenommen, eine Epidemie von bösartigem Gutrot ist in der Bevölkerung weit verbreitet. Glücklicherweise können wir jeden leicht auf ein Merkmal ( ) untersuchen, und wir wissen, dass die Wahrscheinlichkeit der Entwicklung von MG linear von abhängt , p = γ t für einen Parameter γ ( 0 ≤ γ ≤ 1 ) . Das Merkmal t ist gleichmäßig in der Bevölkerung verteilt.$t$$0\leq t \leq 1$$t$$p=\gamma t$$\gamma$$0\leq \gamma \leq 1$$t$

Zum Glück gibt es einen Impfstoff. Leider ist es teuer und die Nebenwirkungen sind sehr unangenehm. (Ich lasse Ihre Fantasie die Details liefern.) Sie sind jedoch besser, als unter MG zu leiden.

Im Interesse der Abstraktion gehe ich davon aus, dass es für einen bestimmten Patienten aufgrund seines Merkmalswerts $t$ tatsächlich nur zwei mögliche Vorgehensweisen gibt : entweder impfen oder nicht impfen.

Somit ist die Frage: Wie sollen wir entscheiden , wer zu impfen und die nicht zu, da $t$ ? Wir werden diesbezüglich nützlich sein und die niedrigsten erwarteten Gesamtkosten anstreben. Es ist offensichtlich, dass es darauf ankommt, einen Schwellenwert $\theta$ zu wählen und jeden mit $t\geq\theta$ zu impfen .

---

Modell und Entscheidung 1 sind genauigkeitsabhängig. Passen Sie ein Modell an. Glücklicherweise haben wir bereits kennen das Modell. Wählen Sie den Schwellenwert $\theta$ , der die Genauigkeit bei der Klassifizierung von Patienten maximiert , und impfen Sie alle mit $t\geq \theta$ . Wir sehen leicht, dass $\theta=\frac{1}{2\gamma}$ ist die magische Zahl - jeder mit$t\geq \theta$hat eine höhere Wahrscheinlichkeit, MG zu kontrahieren als nicht, und umgekehrt, sodass dieseKlassifizierungswahrscheinlichkeitsschwelledie Genauigkeit maximiert. Unter der Annahme ausgewogener Klassen,$\gamma=1$, werden wir die Hälfte der Bevölkerung impfen. Lustigerweise, wenn$\gamma<\frac{1}{2}$ , wir werdenniemandenimpfen. (Wir sind hauptsächlich an ausgeglichenen Klassen interessiert. Lassen Sie uns also außer Acht, dass wir nur einen Teil der Bevölkerung eines schrecklichen schmerzhaften Todes sterben lassen.)

Dies berücksichtigt natürlich nicht die unterschiedlichen Kosten einer Fehlklassifizierung.

---

Modell und Entscheidung 2 nutzen sowohl unsere probabilistische Vorhersage ("angesichts Ihres Merkmals $t$ beträgt Ihre Wahrscheinlichkeit, MG zu bekommen, $\gamma t$ ") als auch die Kostenstruktur.

Hier ist zunächst eine kleine Grafik. Die horizontale Achse gibt das Merkmal an, die vertikale Achse die MG-Wahrscheinlichkeit. Das schattierte Dreieck gibt den Anteil der Bevölkerung an, die MG bekommen wird. Die vertikale Linie gibt ein bestimmtes $\theta$ . Die horizontale gestrichelte Linie bei $\gamma\theta$ vereinfacht die folgenden Berechnungen etwas. Wir nehmen $\gamma>\frac{1}{2}$ , nur um das Leben leichter zu machen.

Geben wir unsere Kostennamen an und berechnen ihre Beiträge zu den erwarteten Gesamtkosten unter Berücksichtigung von $\theta$ und $\gamma$ (und der Tatsache, dass das Merkmal gleichmäßig in der Bevölkerung verteilt ist).

• Lassen Sie $c^+_+$ bezeichnen die Kosten für einen Patienten, geimpft und hätte MG zusammengezogen. Bei $\theta$ ist der Anteil der Bevölkerung, der diese Kosten verursacht, das schattierte Trapez unten rechts mit der Fläche $$(1-\theta)\gamma\theta + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta).$$
• Lassen Sie $c^-_+$ bezeichnen die Kosten für einen Patienten, geimpft und würde nicht haben MG zusammengezogen. Bei $\theta$ ist der Anteil der Bevölkerung, der diese Kosten verursacht, das nicht schattierte Trapez oben rechts mit der Fläche $$(1-\theta)(1-\gamma) + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta).$$
• Lassen Sie $c^-_-$ die Kosten für einen Patienten bezeichnen , wer nicht geimpft und würde nicht haben MG zusammengezogen. Bei $\theta$ ist der Anteil der Bevölkerung, der diese Kosten verursacht, das nicht schattierte Trapez oben links mit der Fläche $$\theta(1-\gamma\theta) + \frac{1}{2}\theta\gamma\theta.$$
• Lassen Sie $c^+_-$ bezeichnen die Kosten für einen Patienten, ist nicht geimpft und hätte MG zusammengezogen. Bei $\theta$ ist der Anteil der Bevölkerung, der diese Kosten verursacht, das schattierte Dreieck unten links mit Fläche $$\frac{1}{2}\theta\gamma\theta.$$

(In jedem Trapez berechne ich zuerst die Fläche des Rechtecks ​​und füge dann die Fläche des Dreiecks hinzu.)

Die erwarteten Gesamtkosten betragen

$$c^+_+\bigg((1-\theta)\gamma\theta + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta)\bigg) + c^-_+\bigg((1-\theta)(1-\gamma) + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta)\bigg) + c^-_-\bigg(\theta(1-\gamma\theta) + \frac{1}{2}\theta\gamma\theta\bigg) + c^+_-\frac{1}{2}\theta\gamma\theta.$$

Wenn wir die Ableitung differenzieren und auf Null setzen, erhalten wir, dass die erwarteten Kosten durch θ ∗ = c - + - c - - minimiert werden.

$$\theta^\ast = \frac{c^-_+-c^-_-}{\gamma(c^+_-+c^-_+-c^+_+-c^-_-)}.$$

Dies ist nur gleich dem Genauigkeitsmaximierungswert von $\theta$ für eine sehr spezifische Kostenstruktur, nämlich genau dann, wenn

$$\frac{1}{2\gamma} = \frac{c^-_+-c^-_-}{\gamma(c^+_-+c^-_+-c^+_+-c^-_-)},$$ oder $$\frac{1}{2} = \frac{c^-_+-c^-_-}{c^+_-+c^-_+-c^+_+-c^-_-}.$$

Nehmen wir als Beispiel an, dass $\gamma=1$ für ausgeglichene Klassen ist und dass die Kosten

$$c^+_+ = 1, \quad c^-_+=2, \quad c^+_-=10, \quad c^-_-=0.$$ Dann ist die Genauigkeitsmaximierung$\theta=\frac{1}{2}$$1.875$$\theta=\frac{2}{11}$$1.318$

In diesem Beispiel führte die Grundlage unserer Entscheidungen für nicht-probabilistische Klassifikationen, die die Genauigkeit maximierten, zu mehr Impfungen und höheren Kosten als die Verwendung einer Entscheidungsregel, bei der die unterschiedlichen Kostenstrukturen explizit im Rahmen einer probabilistischen Vorhersage verwendet wurden.

---

Fazit: Genauigkeit ist nur dann ein gültiges Entscheidungskriterium, wenn

• Es gibt eine Eins-zu-Eins-Beziehung zwischen Klassen und möglichen Aktionen
• und die Kosten für Aktionen, die auf Klassen angewendet werden, folgen einer sehr spezifischen Struktur.

Im allgemeinen Fall stellt die Bewertung der Genauigkeit eine falsche Frage, und die Maximierung der Genauigkeit ist ein sogenannter Typ-III-Fehler: die richtige Antwort auf die falsche Frage.

---

R-Code:

rm(list=ls())
gamma <- 0.7

cost_treated_positive <- 1          # cost of treatment, side effects unimportant
cost_treated_negative <- 2          # cost of treatment, side effects unnecessary
cost_untreated_positive <- 10       # horrible, painful death
cost_untreated_negative <- 0        # nothing

expected_cost <- function ( theta ) {
    cost_treated_positive * ( (1-theta)*theta*gamma + (1-theta)*(gamma-gamma*theta)/2 ) +
    cost_treated_negative * ( (1-theta)*(1-gamma) + (1-theta)*(gamma-gamma*theta)/2 ) +
    cost_untreated_negative *( theta*(1-gamma*theta) + theta*gamma*theta/2 ) +
    cost_untreated_positive * theta*gamma*theta/2
}

(theta <- optim(par=0.5,fn=expected_cost,lower=0,upper=1,method="L-BFGS-B")$par)
(cost_treated_negative-cost_untreated_negative)/
    (gamma*(cost_treated_negative+cost_untreated_positive-cost_treated_positive-cost_untreated_negative))

plot(c(0,1),c(0,1),type="n",bty="n",xaxt="n",xlab="Trait t",yaxt="n",ylab="MG probability")
rect(0,0,1,1)
axis(1,c(0,theta,1),c(0,"theta",1),lty=0,line=-1)
axis(2,c(0,1),lty=0,line=-1,las=1)
axis(4,c(0,gamma,1),c(0,"gamma",1),lty=0,line=-1.8,las=1)
polygon(c(0,1,1),c(0,0,gamma),col="lightgray")
abline(v=theta,col="red",lwd=2)
abline(h=gamma*theta,lty=2,col="red",lwd=2)

expected_cost(1/(2*gamma))
expected_cost(theta)

Es könnte sich lohnen, Stephens hervorragender Antwort ein weiteres, vielleicht einfacheres Beispiel hinzuzufügen.

Betrachten wir einen medizinischen Test, dessen Ergebnis normalerweise sowohl bei kranken als auch bei gesunden Menschen mit unterschiedlichen Parametern verteilt ist (aber der Einfachheit halber nehmen wir Homoskedastizität an, dh die Varianz ist dieselbe):

$$\begin{gather*}T \mid D \ominus \sim \mathcal{N}\left(\mu_{-},\sigma^2\right)\\T \mid D \oplus \sim \mathcal{N}\left(\mu_{+},\sigma^2\right)\end{gather*}.$$ $p$$D\oplus\sim Bern\left(p\right)$

$b$$b$

$$\begin{pmatrix} & D\oplus & D\ominus\\ T\oplus & p\left(1-\Phi_{+}\left(b\right)\right) & \left(1-p\right)\left(1-\Phi_{-}\left(b\right)\right)\\ T\ominus & p\Phi_{+}\left(b\right) & \left(1-p\right)\Phi_{-}\left(b\right)\\ \end{pmatrix}.$$

---

Genauigkeitsbasierter Ansatz

$$p\left(1-\Phi_{+}\left(b\right)\right)+\left(1-p\right)\Phi_{-}\left(b\right),$$

$b$$\sqrt{1\pi\sigma^2}$

$$\begin{gather*} -p\varphi_{+}\left(b\right)+\varphi_{-}\left(b\right)-p\varphi_{-}\left(b\right)=0\\ e^{-\frac{\left(b-\mu_{-}\right)^2}{2\sigma^2}}\left[\left(1-p\right)-pe^{-\frac{2b\left(\mu_{-}-\mu_{+}\right)+\left(\mu_{+}^2-\mu_{-}^2\right)}{2\sigma^2}}\right]=0\end{gather*}$$ $$\begin{gather*}\left(1-p\right)-pe^{-\frac{2b\left(\mu_{-}-\mu_{+}\right)+\left(\mu_{+}^2-\mu_{-}^2\right)}{2\sigma^2}}=0\\-\frac{2b\left(\mu_{-}-\mu_{+}\right)+\left(\mu_{+}^2-\mu_{-}^2\right)}{2\sigma^2}=\log\frac{1-p}{p}\\ 2b\left(\mu_{+}-\mu_{-}\right)+\left(\mu_{-}^2-\mu_{+}^2\right)=2\sigma^2\log\frac{1-p}{p}\\ \end{gather*}$$ $$b^{\ast}=\frac{\left(\mu_{+}^2-\mu_{-}^2\right)+2\sigma^2\log\frac{1-p}{p}}{2\left(\mu_{+}-\mu_{-}\right)}=\frac{\mu_{+}+\mu_{-}}{2}+\frac{\sigma^2}{\mu_{+}-\mu_{-}}\log\frac{1-p}{p}.$$

Beachten Sie, dass dies natürlich nicht von den Kosten abhängt.

Wenn die Klassen ausgeglichen sind, ist das Optimum der Durchschnitt der mittleren Testwerte bei kranken und gesunden Menschen, andernfalls wird es aufgrund des Ungleichgewichts verschoben.

---

Kostenbasierter Ansatz

$$c_{+}^{+}p\left(1-\Phi_{+}\left(b\right)\right) + c_{+}^{-}\left(1-p\right)\left(1-\Phi_{-}\left(b\right)\right) + c_{-}^{+} p\Phi_{+}\left(b\right) + c_{-}^{-} \left(1-p\right)\Phi_{-}\left(b\right).$$ $b$ $$\begin{gather*} -c_{+}^{+} p \varphi_{+}\left(b\right)-c_{+}^{-}\left(1-p\right)\varphi_{-}\left(b\right)+c_{-}^{+}p\varphi_{+}\left(b\right)+c_{-}^{-}\left(1-p\right)\varphi_{-}\left(b\right)=\\ =\varphi_{+}\left(b\right)p\left(c_{-}^{+}-c_{+}^{+}\right)+\varphi_{-}\left(b\right)\left(1-p\right)\left(c_{-}^{-}-c_{+}^{-}\right)=\\ = \varphi_{+}\left(b\right)pc_d^{+}-\varphi_{-}\left(b\right)\left(1-p\right)c_d^{-}= 0,\end{gather*}$$ unter Verwendung der Notation, die ich in meinen Kommentaren unter Stephens Antwort eingeführt habe, dh $c_d^{+}=c_{-}^{+}-c_{+}^{+}$ und $c_d^{-}=c_{+}^{-}-c_{-}^{-}$.

Die optimale Schwelle ergibt sich daher aus der Lösung der Gleichung

$$\boxed{\frac{\varphi_{+}\left(b\right)}{\varphi_{-}\left(b\right)}=\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}}.$$ Zwei Dinge sollten hier beachtet werden:

1. Dieses Ergebnis ist völlig allgemein gehalten und funktioniert für jede Verteilung der Testergebnisse, nicht nur normal. ($\varphi$ in diesem Fall bedeutet natürlich die Wahrscheinlichkeitsdichtefunktion der Verteilung, nicht die normale Dichte.)
2. Was auch immer die Lösung für $b$ ist, es ist sicherlich eine Funktion von $\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}$. (Dh wir sehen sofort, wie wichtig Kosten sind - zusätzlich zum Klassenungleichgewicht!)

Es würde mich wirklich interessieren, ob diese Gleichung eine generische Lösung für hat $b$ (parametrisiert durch die $\varphi$s), aber ich wäre überrascht.

Trotzdem können wir es normal ausarbeiten! $\sqrt{2\pi\sigma^2}$s Abbrechen auf der linken Seite, also haben wir

$$\begin{gather*} e^{-\frac{1}{2}\left(\frac{\left(b-\mu_{+}\right)^2}{\sigma^2}-\frac{\left(b-\mu_{-}\right)^2}{\sigma^2}\right)}=\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}} \\ \left(b-\mu_{-}\right)^2-\left(b-\mu_{+}\right)^2 =2\sigma^2 \log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}} \\ 2b\left(\mu_{+}-\mu_{-}\right)+\left(\mu_{-}^2-\mu_{+}^2\right) =2\sigma^2 \log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}\end{gather*}$$ daher ist die Lösung $$b^{\ast}=\frac{\left(\mu_{+}^2-\mu_{-}^2\right)+2\sigma^2 \log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}}{2\left(\mu_{+}-\mu_{-}\right)}=\frac{\mu_{+}+\mu_{-}}{2}+\frac{\sigma^2}{\mu_{+}-\mu_{-}}\log \frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}.$$

(Vergleichen Sie das vorherige Ergebnis! Wir sehen, dass sie genau dann gleich sind, wenn $c_d^{-}=c_d^{+}$dh die Unterschiede zwischen den Kosten für die Fehlklassifizierung und den Kosten für die korrekte Klassifizierung sind bei kranken und gesunden Menschen gleich.)

---

Eine kurze Demonstration

Sagen wir $c_{-}^{-}=0$ (es ist medizinisch ganz natürlich), und das $c_{+}^{+}=1$ (Wir können es immer erhalten, indem wir die Kosten mit teilen $c_{+}^{+}$dh durch Messen aller Kosten in $c_{+}^{+}$Einheiten). Nehmen wir an, die Prävalenz ist$p=0.2$. Sagen wir das auch$\mu_{-}=9.5$, $\mu_{+}=10.5$ und $\sigma=1$.

In diesem Fall:

library( data.table )
library( lattice )

cminusminus <- 0
cplusplus <- 1
p <- 0.2
muminus <- 9.5
muplus <- 10.5
sigma <- 1

res <- data.table( expand.grid( b = seq( 6, 17, 0.1 ),
                                cplusminus = c( 1, 5, 10, 50, 100 ),
                                cminusplus = c( 2, 5, 10, 50, 100 ) ) )
res$cost <- cplusplus*p*( 1-pnorm( res$b, muplus, sigma ) ) +
  res$cplusminus*(1-p)*(1-pnorm( res$b, muminus, sigma ) ) +
  res$cminusplus*p*pnorm( res$b, muplus, sigma ) +
  cminusminus*(1-p)*pnorm( res$b, muminus, sigma )

xyplot( cost ~ b | factor( cminusplus ), groups = cplusminus, ylim = c( -1, 22 ),
        data = res, type = "l", xlab = "Threshold",
        ylab = "Expected overall cost", as.table = TRUE,
        abline = list( v = (muplus+muminus)/2+
                         sigma^2/(muplus-muminus)*log((1-p)/p) ),
        strip = strip.custom( var.name = expression( {"c"^{"+"}}["-"] ),
                              strip.names = c( TRUE, TRUE ) ),
        auto.key = list( space = "right", points = FALSE, lines = TRUE,
                         title = expression( {"c"^{"-"}}["+"] ) ),
        panel = panel.superpose, panel.groups = function( x, y, col.line, ... ) {
          panel.xyplot( x, y, col.line = col.line, ... )
          panel.points( x[ which.min( y ) ], min( y ), pch = 19, col = col.line )
        } )

Das Ergebnis ist (Punkte stellen die minimalen Kosten dar und die vertikale Linie zeigt den optimalen Schwellenwert mit dem auf Genauigkeit basierenden Ansatz):

Wir können sehr gut sehen, wie sich das kostenbasierte Optimum vom genauigkeitsbasierten Optimum unterscheiden kann. Es ist lehrreich darüber nachzudenken, warum: wenn es teurer ist, kranke Menschen fälschlicherweise gesund zu klassifizieren als umgekehrt ($c_{-}^{+}$ ist hoch, $c_{+}^{-}$ ist niedrig) als der Schwellenwert sinkt, da wir es vorziehen, leichter in die Kategorie krank einzustufen, wenn es andererseits teurer ist, ein gesundes Volk zu fälschen, das fälschlicherweise krank ist, als umgekehrt ($c_{-}^{+}$ ist niedrig, $c_{+}^{-}$ist hoch) als der Schwellenwert steigt, da wir es vorziehen, leichter in die Kategorie gesund zu klassifizieren. (Überprüfen Sie diese auf der Abbildung!)

---

Ein reales Beispiel

Schauen wir uns ein empirisches Beispiel anstelle einer theoretischen Ableitung an. Dieses Beispiel unterscheidet sich grundlegend von zwei Aspekten:

• Anstatt Normalität anzunehmen, werden wir einfach die empirischen Daten ohne eine solche Annahme verwenden.
• Anstatt einen einzelnen Test und seine Ergebnisse in eigenen Einheiten zu verwenden, werden wir mehrere Tests verwenden (und diese mit einer logistischen Regression kombinieren). Der Schwellenwert wird für die endgültige vorhergesagte Wahrscheinlichkeit angegeben. Dies ist eigentlich der bevorzugte Ansatz, siehe Kapitel 19 - Diagnose - in Frank Harrells BBR .

Der Datensatz ( acathaus dem PaketHmisc ) stammt aus der Datenbank für kardiovaskuläre Erkrankungen der Duke University und enthält, ob der Patient eine signifikante Koronarerkrankung hatte, wie durch Herzkatheterisierung festgestellt. Dies ist unser Goldstandard, dh der wahre Krankheitsstatus und der "Test" "wird die Kombination aus Alter, Geschlecht, Cholesterinspiegel und Dauer der Symptome des Probanden sein:

library( rms )
library( lattice )
library( latticeExtra )
library( data.table )

getHdata( "acath" )
acath <- acath[ !is.na( acath$choleste ), ]
dd <- datadist( acath )
options( datadist = "dd" )

fit <- lrm( sigdz ~ rcs( age )*sex + rcs( choleste ) + cad.dur, data = acath )

Es lohnt sich, die vorhergesagten Risiken im Logit-Maßstab aufzuzeichnen, um zu sehen, wie normal sie sind (im Wesentlichen haben wir dies zuvor mit einem einzigen Test angenommen!):

densityplot( ~predict( fit ), groups = acath$sigdz, plot.points = FALSE, ref = TRUE,
             auto.key = list( columns = 2 ) )

Nun, sie sind kaum normal ...

Lassen Sie uns fortfahren und die erwarteten Gesamtkosten berechnen:

ExpectedOverallCost <- function( b, p, y, cplusminus, cminusplus,
                                 cplusplus = 1, cminusminus = 0 ) {
  sum( table( factor( p>b, levels = c( FALSE, TRUE ) ), y )*matrix(
    c( cminusminus, cplusminus, cminusplus, cplusplus ), nc = 2 ) )
}

table( predict( fit, type = "fitted" )>0.5, acath$sigdz )

ExpectedOverallCost( 0.5, predict( fit, type = "fitted" ), acath$sigdz, 2, 4 )

Und zeichnen wir es für alle möglichen Kosten auf (ein rechnerischer Hinweis: Wir müssen nicht gedankenlos durch Zahlen von 0 bis 1 iterieren, wir können die Kurve perfekt rekonstruieren, indem wir sie für alle eindeutigen Werte der vorhergesagten Wahrscheinlichkeiten berechnen):

ps <- sort( unique( c( 0, 1, predict( fit, type = "fitted" ) ) ) )

xyplot( sapply( ps, ExpectedOverallCost,
                p = predict( fit, type = "fitted" ), y = acath$sigdz,
                cplusminus = 2, cminusplus = 4 ) ~ ps, type = "l", xlab = "Threshold",
        ylab = "Expected overall cost", panel = function( x, y, ... ) {
          panel.xyplot( x, y, ... )
          panel.points( x[ which.min( y ) ], min( y ), pch = 19, cex = 1.1 )
          panel.text( x[ which.min( y ) ], min( y ), round( x[ which.min( y ) ], 3 ),
                      pos = 3 )
        } )

Wir können sehr gut sehen, wo wir den Schwellenwert setzen sollten, um die erwarteten Gesamtkosten zu optimieren (ohne irgendwo Sensitivität, Spezifität oder Vorhersagewerte zu verwenden!). Dies ist der richtige Ansatz.

Es ist besonders lehrreich, diese Metriken gegenüberzustellen:

ExpectedOverallCost2 <- function( b, p, y, cplusminus, cminusplus,
                                  cplusplus = 1, cminusminus = 0 ) {
  tab <- table( factor( p>b, levels = c( FALSE, TRUE ) ), y )
  sens <- tab[ 2, 2 ] / sum( tab[ , 2 ] )
  spec <- tab[ 1, 1 ] / sum( tab[ , 1 ] )
  c( `Expected overall cost` = sum( tab*matrix( c( cminusminus, cplusminus, cminusplus,
                                                   cplusplus ), nc = 2 ) ),
     Sensitivity = sens,
     Specificity = spec,
     PPV = tab[ 2, 2 ] / sum( tab[ 2, ] ),
     NPV = tab[ 1, 1 ] / sum( tab[ 1, ] ),
     Accuracy = 1 - ( tab[ 1, 1 ] + tab[ 2, 2 ] )/sum( tab ),
     Youden = 1 - ( sens + spec - 1 ),
     Topleft = ( 1-sens )^2 + ( 1-spec )^2
  )
}

ExpectedOverallCost2( 0.5, predict( fit, type = "fitted" ), acath$sigdz, 2, 4 )

res <- melt( data.table( ps, t( sapply( ps, ExpectedOverallCost2,
                                        p = predict( fit, type = "fitted" ),
                                        y = acath$sigdz,
                                        cplusminus = 2, cminusplus = 4 ) ) ),
             id.vars = "ps" )

p1 <- xyplot( value ~ ps, data = res, subset = variable=="Expected overall cost",
              type = "l", xlab = "Threshold", ylab = "Expected overall cost",
              panel=function( x, y, ... ) {
                panel.xyplot( x, y,  ... )
                panel.abline( v = x[ which.min( y ) ],
                              col = trellis.par.get()$plot.line$col )
                panel.points( x[ which.min( y ) ], min( y ), pch = 19 )
              }  )
p2 <- xyplot( value ~ ps, groups = variable,
              data = droplevels( res[ variable%in%c( "Expected overall cost",
                                                     "Sensitivity",
                                                     "Specificity", "PPV", "NPV" ) ] ),
              subset = variable%in%c( "Sensitivity", "Specificity", "PPV", "NPV" ),
              type = "l", xlab = "Threshold", ylab = "Sensitivity/Specificity/PPV/NPV",
              auto.key = list( columns = 3, points = FALSE, lines = TRUE ) )
doubleYScale( p1, p2, use.style = FALSE, add.ylab2 = TRUE )

Wir können jetzt die Metriken analysieren, für die manchmal speziell geworben wird, um einen optimalen Cutoff ohne Kosten zu erzielen, und dies unserem kostenbasierten Ansatz gegenüberstellen! Verwenden wir die drei am häufigsten verwendeten Metriken:

• Genauigkeit (maximale Genauigkeit)
• Youden-Regel (maximieren $Sens+Spec-1$)
• Topleft-Regel (minimieren $\left(1-Sens\right)^2+\left(1-Spec\right)^2$)

(Der Einfachheit halber subtrahieren wir die obigen Werte von 1 für die Youden- und die Genauigkeitsregel, damit wir überall ein Minimierungsproblem haben.)

Lassen Sie uns die Ergebnisse sehen:

p3 <- xyplot( value ~ ps, groups = variable,
              data = droplevels( res[ variable%in%c( "Expected overall cost", "Accuracy",
                                                     "Youden", "Topleft"  ) ] ),
              subset = variable%in%c( "Accuracy", "Youden", "Topleft"  ),
              type = "l", xlab = "Threshold", ylab = "Accuracy/Youden/Topleft",
              auto.key = list( columns = 3, points = FALSE, lines = TRUE ),
              panel = panel.superpose, panel.groups = function( x, y, col.line, ... ) {
                panel.xyplot( x, y, col.line = col.line, ... )
                panel.abline( v = x[ which.min( y ) ], col = col.line )
                panel.points( x[ which.min( y ) ], min( y ), pch = 19, col = col.line )
              } )
doubleYScale( p1, p3, use.style = FALSE, add.ylab2 = TRUE )

Dies betrifft natürlich eine bestimmte Kostenstruktur, $c_{-}^{-}=0$, $c_{+}^{+}=1$, $c_{+}^{-}=2$, $c_{-}^{+}=4$(Dies ist natürlich nur für die optimale Kostenentscheidung von Bedeutung). Um den Effekt der Kostenstruktur zu untersuchen, wählen wir nur den optimalen Schwellenwert (anstatt die gesamte Kurve zu verfolgen), zeichnen ihn jedoch als Funktion der Kosten auf. Genauer gesagt, wie wir bereits gesehen haben, hängt die optimale Schwelle von den vier Kosten nur durch die ab$c_d^{-}/c_d^{+}$ Verhältnis, also zeichnen wir den optimalen Cutoff als Funktion davon zusammen mit den normalerweise verwendeten Metriken, die keine Kosten verwenden:

res2 <- data.frame( rat = 10^( seq( log10( 0.02 ), log10( 50 ), length.out = 500 ) ) )
res2$OptThreshold <- sapply( res2$rat,
                             function( rat ) ps[ which.min(
                               sapply( ps, Vectorize( ExpectedOverallCost, "b" ),
                                       p = predict( fit, type = "fitted" ),
                                       y = acath$sigdz,
                                       cplusminus = rat,
                                       cminusplus = 1,
                                       cplusplus = 0 ) ) ] )

xyplot( OptThreshold ~ rat, data = res2, type = "l", ylim = c( -0.1, 1.1 ),
        xlab = expression( {"c"^{"-"}}["d"]/{"c"^{"+"}}["d"] ), ylab = "Optimal threshold",
        scales = list( x = list( log = 10, at = c( 0.02, 0.05, 0.1, 0.2, 0.5, 1,
                                                   2, 5, 10, 20, 50 ) ) ),
        panel = function( x, y, resin = res[ ,.( ps[ which.min( value ) ] ),
                                             .( variable ) ], ... ) {
          panel.xyplot( x, y, ... )
          panel.abline( h = resin[variable=="Youden"] )
          panel.text( log10( 0.02 ), resin[variable=="Youden"], "Y", pos = 3 )
          panel.abline( h = resin[variable=="Accuracy"] )
          panel.text( log10( 0.02 ), resin[variable=="Accuracy"], "A", pos = 3 )
          panel.abline( h = resin[variable=="Topleft"] )
          panel.text( log10( 0.02 ), resin[variable=="Topleft"], "TL", pos = 1 )
        } )

Horizontale Linien geben die Ansätze an, die keine Kosten verwenden (und daher konstant sind).

Wir sehen wieder gut, dass mit steigenden zusätzlichen Kosten für eine Fehlklassifizierung in der gesunden Gruppe im Vergleich zu denen der erkrankten Gruppe die optimale Schwelle steigt: Wenn wir wirklich nicht wollen, dass gesunde Menschen als krank eingestuft werden, werden wir einen höheren Grenzwert verwenden (und umgekehrt natürlich!).

Und schließlich sehen wir noch einmal, warum Methoden, die keine Kosten verwenden, nicht immer optimal sind ( und auch nicht können! ).