Ejemplo en el que el uso de la precisión como medida de resultado llevará a una conclusión errónea

Question

Estoy estudiando diversas medidas de rendimiento para modelos predictivos. Se ha escrito mucho sobre los problemas de utilizar la precisión, en lugar de algo más continuo para evaluar el rendimiento de los modelos. Frank Harrell ofrece una ejemplo cuando añadir una variable informativa a un modelo conlleva un descenso de la precisión, conclusión claramente contraintuitiva y errónea.

Sin embargo, en este caso, esto parece ser causado por tener clases desequilibradas, y por lo tanto se puede resolver simplemente mediante el uso de precisión equilibrada en su lugar ((sens+spec)/2). ¿Hay algún ejemplo en el que el uso de la precisión en un conjunto de datos equilibrado lleve a conclusiones claramente erróneas o contrarias a la intuición?

Editar

Busco algo en lo que la precisión disminuya incluso cuando el modelo sea claramente mejor, o que el uso de la precisión conduzca a una selección de falsos positivos de algunas características. Es fácil hacer ejemplos falsos negativos, donde la precisión es la misma para dos modelos donde uno es claramente mejor usando otros criterios.

Answer 1

Haré trampa.

En concreto, he argumentado a menudo (p. ej, aquí ) que el estadística parte del modelado y la predicción se extiende sólo a hacer probabilístico predicciones de pertenencia a clases (o dar densidades de predicción, en el caso de la predicción numérica). Tratamiento de una instancia específica como si perteneciera a una clase específica (o punto predicciones en el caso numérico), ya no es propiamente estadística. Forma parte de la decisión aspecto teórico.

Y las decisiones no sólo deben basarse en la predicción probabilística, sino también en los costes de las clasificaciones erróneas y en una serie de factores. otras acciones posibles . Por ejemplo, incluso si sólo tiene dos clases posibles, "enfermo" frente a "sano", podría tener una amplia gama de posibles acciones dependiendo de cómo probablemente es que un paciente padezca la enfermedad, desde enviarle a casa porque es casi seguro que está sano, hasta darle dos aspirinas, hacerle pruebas adicionales, llamar inmediatamente a una ambulancia y ponerle respiración asistida.

Evaluar la precisión presupone una decisión de este tipo. _La precisión como métrica de evaluación de la clasificación es una error de categoría ._

Así que, para responder a su pregunta, recorreré el camino de un error de categoría de este tipo. Consideraremos un escenario sencillo con clases equilibradas en el que clasificar sin tener en cuenta los costes de una clasificación errónea nos llevará a un error grave.

---

Supongamos que una epidemia de Gutrot Maligna hace estragos entre la población. Afortunadamente, podemos examinar a todos fácilmente por algún rasgo $t$ ( $0\leq t \leq 1$ ), y sabemos que la probabilidad de desarrollar MG depende linealmente de $t$ , $p=\gamma t$ para algún parámetro $\gamma$ ( $0\leq \gamma \leq 1$ ). El rasgo $t$ se distribuye uniformemente en la población.

Afortunadamente, existe una vacuna. Desgraciadamente, es cara, y los efectos secundarios son muy incómodo. (Dejaré que su imaginación aporte los detalles.) Sin embargo, es mejor que sufran de MG.

En aras de la abstracción, propongo que sólo hay dos posibles cursos de acción para cualquier paciente, dado su valor de rasgo $t$ O se vacuna, o no se vacuna.

Así pues, la pregunta es: ¿cómo debemos decidir a quién vacunar y a quién no, teniendo en cuenta $t$ ? Estaremos utilitario sobre esto y aspirar a tener los costes totales previstos más bajos. Es obvio que esto se reduce a elegir un umbral $\theta$ y vacunar a todos con $t\geq\theta$ .

---

El modelo y la decisión 1 se basan en la precisión. Ajustar un modelo. Afortunadamente, ya conozca el modelo. Elegir el umbral $\theta$ que maximiza la precisión al clasificar a los pacientes y vacunar a todos con $t\geq \theta$ . Vemos fácilmente que $\theta=\frac{1}{2\gamma}$ es el número mágico - todo el mundo con $t\geq \theta$ tiene más posibilidades de contraer MG que no, y viceversa, por lo que esta umbral de probabilidad de clasificación maximizará la precisión. Suponiendo clases equilibradas, $\gamma=1$ vacunaremos a la mitad de la población. Curiosamente, si $\gamma<\frac{1}{2}$ vacunaremos nadie . (Nos interesan sobre todo las clases equilibradas, así que no tengamos en cuenta que acabamos de dejar que parte de la población muera de una Horrible Muerte Dolorosa).

Huelga decir que esto no tiene en cuenta los costes diferenciales de la clasificación errónea.

---

El modelo y la decisión 2 aprovechan tanto nuestra predicción probabilística ("dado su rasgo $t$ su probabilidad de contraer MG es $\gamma t$ ") y la estructura de costes.

En primer lugar, he aquí un pequeño gráfico. El eje horizontal indica el rasgo, el eje vertical la probabilidad de MG. El triángulo sombreado da la proporción de la población que contraerá MG. La línea vertical indica $\theta$ . La línea horizontal discontinua en $\gamma\theta$ hará que los cálculos a continuación sean un poco más sencillos de seguir. Suponemos que $\gamma>\frac{1}{2}$ sólo para hacer la vida más fácil.

Pongamos nombre a nuestros costes y calculemos su contribución a los costes totales previstos, dado que $\theta$ y $\gamma$ (y el hecho de que el rasgo esté uniformemente distribuido en la población).

• Sea $c^+_+$ denotar el coste para un paciente vacunado que hubiera contraído MG. Dado $\theta$ la proporción de la población que incurre en este coste es el trapecio sombreado de la parte inferior derecha con área $$ (1-\theta)\gamma\theta + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta). $$
• Sea $c^-_+$ denote el coste para un paciente vacunado y que no han contraído MG. Dado $\theta$ la proporción de la población que incurre en este coste es el trapezoide no sombreado de la parte superior derecha con área $$ (1-\theta)(1-\gamma) + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta). $$
• Sea $c^-_-$ denotan el coste para un paciente que es no vacunados y no han contraído MG. Dado $\theta$ la proporción de la población que incurre en este coste es el trapezoide sin sombrear de la parte superior izquierda con área $$ \theta(1-\gamma\theta) + \frac{1}{2}\theta\gamma\theta. $$
• Sea $c^+_-$ denotan el coste para un paciente que es no vacunado y habría contraído la MG. Dado $\theta$ la proporción de la población que incurre en este coste es el triángulo sombreado de la parte inferior izquierda con el área $$ \frac{1}{2}\theta\gamma\theta. $$

(En cada trapecio, calculo primero el área del rectángulo y luego añado el área del triángulo).

Los costes totales previstos son $$ 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. $$

Diferenciando y fijando la derivada a cero, obtenemos que los costes esperados se minimizan mediante $$ \theta^\ast = \frac{c^-_+-c^-_-}{\gamma(c^+_-+c^-_+-c^+_+-c^-_-)}.$$

Esto sólo es igual al valor de maximización de la precisión de $\theta$ para una estructura de costes muy específica, a saber, si y sólo si $$ \frac{1}{2\gamma} = \frac{c^-_+-c^-_-}{\gamma(c^+_-+c^-_+-c^+_+-c^-_-)},$$ o $$ \frac{1}{2} = \frac{c^-_+-c^-_-}{c^+_-+c^-_+-c^+_+-c^-_-}.$$

Por ejemplo, supongamos que $\gamma=1$ para clases equilibradas y que los costes $$ c^+_+ = 1, \quad c^-_+=2, \quad c^+_-=10, \quad c^-_-=0.$$ A continuación, la maximización de la precisión $\theta=\frac{1}{2}$ supondrá unos costes previstos de $1.875$ mientras que la minimización de costes $\theta=\frac{2}{11}$ supondrá unos costes previstos de $1.318$ .

En este ejemplo, basar nuestras decisiones en clasificaciones no probabilísticas que maximizaban la precisión condujo a más vacunaciones y costes más elevados que utilizar una regla de decisión que utilizara explícitamente las estructuras de costes diferenciales en el contexto de una predicción probabilística.

---

En resumen: la precisión sólo es un criterio de decisión válido si

• existe una relación de uno a uno entre clases y posible acciones
• y los costes de las acciones aplicadas a las clases siguen una estructura muy específica.

En el caso general, evaluar la precisión es una pregunta equivocada, y maximizar la precisión es un error del tipo III: dar la respuesta correcta a una pregunta equivocada.

---

Código R:

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)

Answer 2

Quizá valga la pena añadir otro ejemplo, tal vez más directo, a la excelente respuesta de Stephen.

Consideremos una prueba médica, cuyo resultado se distribuye normalmente, tanto en enfermos como en sanos, por supuesto con parámetros diferentes (pero para simplificar, supongamos homocedasticidad, es decir, que la varianza es la misma): $$\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*}.$$ Denotemos la prevalencia de la enfermedad con $p$ (es decir $D\oplus\sim Bern\left(p\right)$ ), por lo que ésta, junto con las anteriores, que son esencialmente distribuciones condicionales, especifica completamente la distribución conjunta.

Así pues, la matriz de confusión con umbral $b$ (es decir, aquellos con resultados superiores a $b$ se clasifican como enfermos) es $$\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}.$$

---

Enfoque basado en la precisión

La precisión es $$p\left(1-\Phi_{+}\left(b\right)\right)+\left(1-p\right)\Phi_{-}\left(b\right),$$

tomamos su derivada respecto a $b$ igual a 0, multiplicar por $\sqrt{1\pi\sigma^2}$ y reorganizar un poco: $$\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*}$$ El primer término no puede ser cero, por lo que la única forma de que el producto sea cero es que el segundo término sea cero: $$\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*}$$ Así que la solución es $$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}.$$

Tenga en cuenta que esto -por supuesto- no depende de los costes.

Si las clases están equilibradas, el óptimo es la media de los valores medios de las pruebas en enfermos y sanos; en caso contrario, se desplaza en función del desequilibrio.

---

Enfoque basado en los costes

Utilizando la notación de Stephen, el coste global esperado es $$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).$$ Tomar su derivada respecto a $b$ y ponerlo a cero: $$\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*}$$ utilizando la notación que introduje en mis comentarios debajo de la respuesta de Stephen, es decir $c_d^{+}=c_{-}^{+}-c_{+}^{+}$ y $c_d^{-}=c_{+}^{-}-c_{-}^{-}$ .

Por tanto, el umbral óptimo viene dado por la solución de la ecuación $$\boxed{\frac{\varphi_{+}\left(b\right)}{\varphi_{-}\left(b\right)}=\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}}.$$ Aquí hay que señalar dos cosas:

1. Este resultado es totalmente genérico y trabaja para cualquier distribución de los resultados de las pruebas, no sólo normal. ( $\varphi$ en ese caso, por supuesto, significa la función de densidad de probabilidad de la distribución, no la densidad normal).
2. Cualquiera que sea la solución para $b$ es, seguramente es una función de $\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}$ . (Es decir, vemos inmediatamente cómo importan los costes, además del desequilibrio de clases).

Me interesaría mucho saber si esta ecuación tiene una solución genérica para $b$ (parametrizado por el $\varphi$ s), pero me sorprendería.

No obstante, ¡podemos solucionarlo con normalidad! $\sqrt{2\pi\sigma^2}$ s se cancelan en el lado izquierdo, por lo que tenemos $$\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*}$$ por lo tanto la solución es $$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^{+}}.$$

(¡Compárelo con el resultado anterior! Vemos que son iguales si y sólo si $c_d^{-}=c_d^{+}$ es decir, las diferencias en el coste de la clasificación errónea en comparación con el coste de la clasificación correcta es el mismo en personas enfermas y sanas).

---

Breve demostración

Digamos que $c_{-}^{-}=0$ (es bastante natural desde el punto de vista médico), y que $c_{+}^{+}=1$ (siempre podemos obtenerlo dividiendo los costes por $c_{+}^{+}$ es decir, midiendo cada coste en $c_{+}^{+}$ unidades). Digamos que la prevalencia es $p=0.2$ . Además, digamos que $\mu_{-}=9.5$ , $\mu_{+}=10.5$ y $\sigma=1$ .

En este caso:

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 )
        } )

El resultado es (los puntos representan el coste mínimo, y la línea vertical muestra el umbral óptimo con el enfoque basado en la precisión):

Podemos ver muy bien cómo el óptimo basado en los costes puede ser diferente del óptimo basado en la precisión. Es instructivo pensar por qué: si es más costoso clasificar a un enfermo erróneamente como sano que al revés ( $c_{-}^{+}$ es alta, $c_{+}^{-}$ es bajo) que el umbral baja, ya que preferimos clasificar más fácilmente en la categoría de enfermo, por otro lado, si es más costoso clasificar erróneamente como enfermo a una persona sana que al revés ( $c_{-}^{+}$ es baja, $c_{+}^{-}$ es alto) que el umbral sube, ya que preferimos clasificar más fácilmente en la categoría sano. (¡Compruébelo en la figura!)

---

Un ejemplo de la vida real

Veamos un ejemplo empírico, en lugar de una derivación teórica. Este ejemplo será diferente básicamente desde dos aspectos:

• En lugar de suponer la normalidad, utilizaremos simplemente los datos empíricos sin ninguna suposición de este tipo.
• En lugar de utilizar una única prueba, y sus resultados en sus propias unidades, utilizaremos varias pruebas (y las combinaremos con una regresión logística). El umbral se dará a la probabilidad predicha final. Este es en realidad el enfoque preferido, véase el capítulo 19 - Diagnóstico - en BBR de Frank Harrell .

En conjunto de datos ( acath del paquete Hmisc ) procede del Banco de Datos de Enfermedades Cardiovasculares de la Universidad de Duke, y contiene si el paciente tenía una enfermedad coronaria significativa, evaluada mediante cateterismo cardíaco, éste será nuestro patrón oro, es decir, el verdadero estado de la enfermedad, y la "prueba" será la combinación de la edad, el sexo, el nivel de colesterol y la duración de los síntomas del sujeto:

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 )

Merece la pena trazar los riesgos predichos en escala logit, para ver hasta qué punto son normales (en esencia, eso es lo que suponíamos anteriormente, ¡con una sola prueba!):

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

Bueno, no son normales...

Sigamos y calculemos el coste global previsto:

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 )

Y vamos a trazarla para todos los costes posibles (una nota computacional: no necesitamos iterar sin sentido a través de números de 0 a 1, podemos reconstruir perfectamente la curva calculándola para todos los valores únicos de probabilidades previstas):

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 )
        } )

Podemos ver muy bien dónde debemos poner el umbral para optimizar el coste global previsto (¡sin utilizar la sensibilidad, la especificidad o los valores predictivos en ningún sitio!). Este es el enfoque correcto.

Resulta especialmente instructivo contrastar estas métricas:

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 )

Ahora podemos analizar esas métricas que a veces se anuncian específicamente como capaces de llegar a un corte óptimo sin costes, ¡y contrastarlo con nuestro enfoque basado en los costes! Utilicemos las tres métricas más utilizadas:

• Precisión (maximizar la precisión)
• Regla de Youden (maximizar $Sens+Spec-1$ )
• Regla de Topleft (minimizar $\left(1-Sens\right)^2+\left(1-Spec\right)^2$ )

(Para simplificar, restaremos los valores anteriores de 1 para la regla de Youden y la de Precisión, de modo que tengamos un problema de minimización en todas partes).

Veamos los resultados:

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 )

Por supuesto, esto se refiere a una estructura de costes específica, $c_{-}^{-}=0$ , $c_{+}^{+}=1$ , $c_{+}^{-}=2$ , $c_{-}^{+}=4$ (obviamente, esto sólo importa para la decisión de coste óptimo). Para investigar el efecto de la estructura de costes, elijamos sólo el umbral óptimo (en lugar de trazar toda la curva) y representémoslo en función de los costes. Más concretamente, como ya hemos visto, el umbral óptimo depende de los cuatro costes sólo a través de la variable $c_d^{-}/c_d^{+}$ por lo que vamos a trazar el límite óptimo en función de esto, junto con las métricas típicamente utilizadas que no utilizan costes:

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 )
        } )

Las líneas horizontales indican los enfoques que no utilizan costes (y que, por tanto, son constantes).

Una vez más, vemos que a medida que aumenta el coste adicional de la clasificación errónea en el grupo sano en comparación con el grupo enfermo, aumenta el umbral óptimo: si realmente no queremos que las personas sanas se clasifiquen como enfermas, utilizaremos un umbral más alto (¡y al revés, por supuesto!).

Y, por último, volvemos a ver por qué los métodos que no utilizan costes no son ( ¡y no puedo! ) sea siempre óptima.