Exemple lorsque l'utilisation de la précision comme mesure de résultat entraînera une conclusion erronée

J'examine différentes mesures de performance pour les modèles prédictifs. Beaucoup a été écrit sur les problèmes d'utilisation de la précision, au lieu de quelque chose de plus continu pour évaluer les performances du modèle. Frank Harrell http://www.fharrell.com/post/class-damage/ fournit un exemple lorsque l'ajout d'une variable informative à un modèle entraînera une baisse de la précision, une conclusion clairement contre-intuitive et erronée.

Cependant, dans ce cas, cela semble être dû au fait d'avoir des classes déséquilibrées, et donc cela peut être résolu simplement en utilisant une précision équilibrée à la place ((sens + spec) / 2). Existe-t-il un exemple où l'utilisation de la précision sur un ensemble de données équilibré conduira à des conclusions clairement fausses ou contre-intuitives?

Éditer

Je recherche quelque chose où la précision chutera même lorsque le modèle est clairement meilleur, ou que l'utilisation de la précision entraînera une sélection faussement positive de certaines fonctionnalités. Il est facile de faire des exemples de faux négatifs, où la précision est la même pour deux modèles où l'un est clairement meilleur en utilisant d'autres critères.

Je vais tricher.

Plus précisément, j'ai souvent soutenu (par exemple ici ) que la partie statistique de la modélisation et de la prédiction s'étend uniquement à la réalisation de prédictions probabilistes pour les appartenances aux classes (ou à la fourniture de densités prédictives, dans le cas de la prévision numérique). Traiter une instance spécifique comme si elle appartenait à une classe spécifique (ou des prédictions ponctuelles dans le cas numérique), n'est plus correctement des statistiques. Cela fait partie de l' aspect théorique de la décision .

Et les décisions ne doivent pas seulement être fondées sur la prédiction probabiliste, mais aussi sur les coûts des erreurs de classification et sur une multitude d' autres actions possibles . Par exemple, même si vous n'avez que deux classes possibles, "malade" contre "sain", vous pourriez avoir un large éventail d'actions possibles selon la probabilité qu'un patient souffre de la maladie, de le renvoyer chez lui parce qu'il est presque certainement en bonne santé, de lui donner deux aspirines, d'exécuter des tests supplémentaires, d'appeler immédiatement une ambulance et de le mettre sous assistance respiratoire.

L'évaluation de la précision suppose une telle décision. La précision en tant que métrique d'évaluation pour la classification est une erreur de catégorie .

Donc, pour répondre à votre question, je vais suivre le chemin d'une telle erreur de catégorie. Nous considérerons un scénario simple avec des classes équilibrées où une classification sans tenir compte des coûts d'une mauvaise classification nous induira en effet mal en erreur.

---

Supposons qu'une épidémie de Gutrot malin sévisse dans la population. Heureusement, nous pouvons filtrer tout le monde facilement pour un trait ( ), et nous savons que la probabilité de développer MG dépend linéairement de , pour un paramètre ( ). Le trait est uniformément distribué dans la population.$t$$0\leq t \leq 1$$t$$p=\gamma t$$\gamma$$0\leq \gamma \leq 1$$t$

Heureusement, il existe un vaccin. Malheureusement, c'est cher et les effets secondaires sont très inconfortables. (Je vais laisser votre imagination vous fournir les détails.) Cependant, ils valent mieux que de souffrir de MG.

Dans l'intérêt de l'abstraction, je postule qu'il n'y a en effet que deux modes d'action possibles pour un patient donné, compte tenu de leur valeur de trait : soit vacciner, soit ne pas vacciner.$t$

Ainsi, la question est: comment devrions-nous décider qui vacciner et qui ne pas le faire, étant donné ? Nous serons utilitaires à ce sujet et visons à avoir les coûts totaux escomptés les plus bas. Il est évident que cela revient à choisir un seuil et à vacciner tout le monde avec le .$t$$\theta$$t\geq\theta$

---

Le modèle et la décision 1 sont axés sur la précision. Montez un modèle. Heureusement, nous connaissons déjà le modèle. Choisissez le seuil qui maximise la précision lors de la classification des patients et vaccinez tout le monde avec le . Nous voyons facilement que est le nombre magique - tout le monde avec a plus de chances de contracter MG que non, et vice versa, donc ce seuil de probabilité de classification maximisera précision. En supposant des classes équilibrées, , nous vaccinerons la moitié de la population. Curieusement, si , nous ne vaccinerons personne$\theta$$t\geq \theta$$\theta=\frac{1}{2\gamma}$$t\geq \theta$$\gamma=1$$\gamma<\frac{1}{2}$. (Nous sommes principalement intéressés par des classes équilibrées, alors ne tenons pas compte du fait que nous laissons juste une partie de la population mourir d'une mort douloureuse horrible.)

Il va sans dire que cela ne prend pas en compte les coûts différentiels de classification erronée.

---

Le modèle et la décision 2 tirent parti à la fois de notre prédiction probabiliste («étant donné votre trait , votre probabilité de contracter MG est ») et de la structure des coûts.$t$$\gamma t$

D'abord, voici un petit graphique. L'axe horizontal donne le trait, l'axe vertical la probabilité MG. Le triangle ombré donne la proportion de la population qui contractera la MG. La ligne verticale donne une particulière . La ligne horizontale en pointillés à rendra les calculs ci-dessous un peu plus simples à suivre. Nous supposons , juste pour vous faciliter la vie.$\theta$$\gamma\theta$$\gamma>\frac{1}{2}$

Donnons nos noms de coûts et calculons leurs contributions aux coûts totaux attendus, étant donné et (et le fait que le trait est uniformément distribué dans la population).$\theta$$\gamma$

• Soit le coût pour un patient vacciné et qui aurait contracté la MG. Étant donné , la proportion de la population qui engage ce coût est le trapèze ombré en bas à droite avec l'aire $c^+_+$$\theta$ $$(1-\theta)\gamma\theta + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta).$$
• Soit le coût pour un patient vacciné qui n'aurait pas contracté de MG. Étant donné , la proportion de la population qui engage ce coût est le trapèze non ombré en haut à droite avec l'aire $c^-_+$$\theta$ $$(1-\theta)(1-\gamma) + \frac{1}{2}(1-\theta)(\gamma-\gamma\theta).$$
• Soit le coût pour un patient qui n'est pas vacciné et n'aurait pas contracté de MG. Étant donné , la proportion de la population qui engage ce coût est le trapèze non ombré en haut à gauche avec la zone $c^-_-$$\theta$ $$\theta(1-\gamma\theta) + \frac{1}{2}\theta\gamma\theta.$$
• Soit le coût pour un patient non vacciné et qui aurait contracté la MG. Étant donné , la proportion de la population qui engage ce coût est le triangle ombré en bas à gauche avec l'aire $c^+_-$$\theta$ $$\frac{1}{2}\theta\gamma\theta.$$

(Dans chaque trapèze, je calcule d'abord l'aire du rectangle, puis j'ajoute l'aire du triangle.)

Les coûts totaux attendus sont

$$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.$$

En différenciant et en mettant la dérivée à zéro, nous obtenons que les coûts attendus sont minimisés par

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

Ceci est seulement égal à la valeur maximisant la précision de pour une structure de coût très spécifique, à savoir si et seulement si ou $\theta$

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

Par exemple, supposons que pour les classes équilibrées et que les coûts soient Ensuite, la précision maximisant des coûts attendus de , tandis que la minimisation des coûts des coûts attendus de .$\gamma=1$

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

Dans cet exemple, baser nos décisions sur des classifications non probabilistes qui maximisent la précision conduit à plus de vaccinations et à des coûts plus élevés que d'utiliser une règle de décision qui utilise explicitement les structures de coûts différentielles dans le contexte d'une prédiction probabiliste.

---

Conclusion: la précision n'est un critère de décision valable que si

• il y a une relation un à un entre les classes et les actions possibles
• et les coûts des actions appliquées aux classes suivent une structure très spécifique.

Dans le cas général, l'évaluation de la précision pose une mauvaise question et la maximisation de la précision est une erreur dite de type III: fournir la bonne réponse à la mauvaise question.

---

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

Il pourrait être utile d'ajouter un autre exemple, peut-être plus simple, à l'excellente réponse de Stephen.

Considérons un test médical, dont le résultat est normalement distribué, à la fois chez les personnes malades et chez les personnes en bonne santé, avec des paramètres différents bien sûr (mais pour simplifier, supposons une homoscédasticité, c'est-à-dire que la variance est la même):Notons la prévalence de la maladie avec (c'est-à-dire ), donc ceci, avec ce qui précède, qui sont essentiellement des distributions conditionnelles, spécifie entièrement la distribution conjointe.

$$\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)$

Ainsi, la matrice de confusion avec le seuil (c'est-à-dire que ceux dont les résultats des tests sont supérieurs à sont classés comme malades) est$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}.$$

---

Approche basée sur la précision

La précision est

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

nous prenons sa dérivée wrt , la mettons égale à 0, multiplions avec et réorganisons un peu: Le premier terme ne peut pas être zéro, donc la seule façon dont le produit peut être zéro est si le deuxième terme est zéro:$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*}$$ La solution est donc $$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}.$$

Notez que cela - bien sûr - ne dépend pas des coûts.

Si les classes sont équilibrées, l'optimum est la moyenne des valeurs de test moyennes chez les personnes malades et en bonne santé, sinon il est déplacé en fonction du déséquilibre.

---

Approche basée sur les coûts

En utilisant la notation de Stephen, le coût global attendu estPrenez sa dérivée par rapport à et fixez-la à zéro:

$$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*}$$ en utilisant la notation que j'ai introduite dans mes commentaires ci-dessous la réponse de Stephen, c.-à-d. et .$c_d^{+}=c_{-}^{+}-c_{+}^{+}$$c_d^{-}=c_{+}^{-}-c_{-}^{-}$

Le seuil optimal est donc donné par la solution de l'équationDeux choses sont à noter ici:

$$\boxed{\frac{\varphi_{+}\left(b\right)}{\varphi_{-}\left(b\right)}=\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}}.$$

1. Ce résultat est totalement générique et fonctionne pour n'importe quelle distribution des résultats de test, pas seulement normal. ( dans ce cas signifie bien sûr la fonction de densité de probabilité de la distribution, pas la densité normale.)$\varphi$
2. Quelle que soit la solution de , c'est sûrement une fonction de . (C'est-à-dire que nous voyons immédiatement à quel point les coûts comptent - en plus du déséquilibre de classe!)$b$$\frac{\left(1-p\right)c_d^{-}}{pc_d^{+}}$

Je serais vraiment intéressé de voir si cette équation a une solution générique pour (paramétrée par les s), mais je serais surpris.$b$$\varphi$

Néanmoins, nous pouvons le faire pour normal! s annuler sur le côté gauche, nous avons donc donc la solution est$\sqrt{2\pi\sigma^2}$

$$\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*}$$ $$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^{+}}.$$

(Comparez-le au résultat précédent! Nous voyons qu'ils sont égaux si et seulement si , c'est-à-dire que les différences de coût de mauvaise classification par rapport au coût d'une classification correcte sont les mêmes chez les malades et les sains. gens.)$c_d^{-}=c_d^{+}$

---

Une courte démonstration

Disons que (c'est assez naturel médicalement), et que (on peut toujours l'obtenir en divisant les coûts par , c'est-à-dire en mesurant chaque coût en unités ). Disons que la prévalence est . Supposons également que , et .$c_{-}^{-}=0$$c_{+}^{+}=1$$c_{+}^{+}$$c_{+}^{+}$$p=0.2$$\mu_{-}=9.5$$\mu_{+}=10.5$$\sigma=1$

Dans ce cas:

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

Le résultat est (les points représentent le coût minimum, et la ligne verticale montre le seuil optimal avec l'approche basée sur la précision):

Nous pouvons très bien voir comment l'optimum basé sur les coûts peut être différent de l'optimum basé sur la précision. Il est instructif de réfléchir à la raison: s'il est plus coûteux de classer une personne malade en mauvaise santé que l'inverse ( est élevé, est faible) que le le seuil descend, car nous préférons classer plus facilement dans la catégorie malade, en revanche, s'il est plus coûteux de classer une personne saine malade par erreur que l'inverse ( est faible, est élevé) que le seuil monte, car nous préférons classer plus facilement dans la catégorie sain. (Vérifiez-les sur la figure!)$c_{-}^{+}$$c_{+}^{-}$$c_{-}^{+}$$c_{+}^{-}$

---

Un exemple concret

Jetons un regard sur un exemple empirique, au lieu d'une dérivation théorique. Cet exemple sera fondamentalement différent de deux aspects:

• Au lieu de supposer la normalité, nous utiliserons simplement les données empiriques sans une telle hypothèse.
• Au lieu d'utiliser un seul test, et ses résultats dans ses propres unités, nous utiliserons plusieurs tests (et les combinerons avec une régression logistique). Un seuil sera donné à la probabilité finale prédite. C'est en fait l'approche préférée, voir le chapitre 19 - Diagnostic - dans le BBR de Frank Harrell .

L' ensemble de données ( acathde l'emballage Hmisc) provient de la banque de données sur les maladies cardiovasculaires de l'Université Duke et indique si le patient souffrait d'une maladie coronarienne importante, telle qu'évaluée par cathétérisme cardiaque, ce sera notre étalon-or, à savoir le véritable état de la maladie et le "test" "sera la combinaison de l'âge, du sexe, du taux de cholestérol et de la durée des symptômes du sujet:

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 )

Il vaut la peine de tracer les risques prévus à l'échelle logit, pour voir à quel point ils sont normaux (essentiellement, c'est ce que nous avons supposé précédemment, avec un seul test!):

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

Eh bien, ce n'est pas normal ...

Continuons et calculons le coût global attendu:

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 )

Et tracons-la pour tous les coûts possibles (une note de calcul: nous n'avons pas besoin d'itérer sans réfléchir par des nombres de 0 à 1, nous pouvons parfaitement reconstruire la courbe en la calculant pour toutes les valeurs uniques des probabilités prédites):

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

Nous pouvons très bien voir où nous devons mettre le seuil pour optimiser le coût global attendu (sans utiliser de sensibilité, de spécificité ou de valeurs prédictives n'importe où!). C'est la bonne approche.

Il est particulièrement instructif de comparer ces mesures:

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 )

Nous pouvons maintenant analyser les métriques qui sont parfois spécifiquement annoncées comme étant capables de proposer un seuil optimal sans frais, et le comparer avec notre approche basée sur les coûts! Utilisons les trois métriques les plus utilisées:

• Précision (maximiser la précision)
• Règle de Youden (maximiser )$Sens+Spec-1$
• Règle du peuple (minimiser )$\left(1-Sens\right)^2+\left(1-Spec\right)^2$

(Par souci de simplicité, nous soustraireons les valeurs ci-dessus de 1 pour le Youden et la règle de précision afin d'avoir un problème de minimisation partout.)

Voyons les résultats:

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 )

Cela concerne bien sûr une structure de coûts spécifique, , , , (cela n'a évidemment d'importance que pour la décision de coût optimal). Pour étudier l'effet de la structure des coûts, choisissons simplement le seuil optimal (au lieu de tracer toute la courbe), mais tracer-le en fonction des coûts. Plus précisément, comme nous l'avons déjà vu, le seuil optimal ne dépend des quatre coûts que par le biais du rapport , nous allons donc tracer le seuil optimal en fonction de cela, ainsi que le type généralement utilisé mesures qui n'utilisent pas de coûts:$c_{-}^{-}=0$$c_{+}^{+}=1$$c_{+}^{-}=2$$c_{-}^{+}=4$$c_d^{-}/c_d^{+}$

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

Les lignes horizontales indiquent les approches qui n'utilisent pas de coûts (et sont donc constantes).

Encore une fois, nous voyons bien qu'à mesure que le coût supplémentaire des erreurs de classification dans le groupe sain augmente par rapport à celui du groupe malade, le seuil optimal augmente: si nous ne voulons vraiment pas que les personnes en bonne santé soient classées comme malades, nous utiliserons un seuil plus élevé (et inversement, bien sûr!).

Et, enfin, nous voyons encore une fois pourquoi ces méthodes qui n'utilisent pas de coûts ne sont pas ( et ne peuvent pas! ) Être toujours optimales.