Sensibilité, Spécificité et ROC pour un test diagnostic

Pour les utilisateurs de R, on commence par charger un certain nombre de modules intéressants.

# le package moderne standard pour manipuler ses données
library(tidyverse)
# pour faire des graphes multiples
library(cowplot)
# pour afficher joliment les tableaux dans un document rmarkdown
library(knitr)

# quelques couleurs manuelles
bleufonce <- "#3d5468"
bleuclair <- "#5b7c98"
rose <- "#ff5555"

Introduction

Test diagnostic

Etant donné un patient \(X\), caractérisé par un ensemble de variables, un test diagnostic est simplement un classificateur, qui classe en deux catégories distinctes :

• Diagnostiqué positif, notés \(D_X = 1\).
• Diagnostiqué négatif, notés \(D_X = 0\).

On ne s’intéresse donc, ni plus, ni moins, qu’à une simple variable de Bernoulli. Tout l’objet de la section suivante consiste à étudier la distribution de cette variable parmi la population de malades et de non-malades.

Threshold sur mesure continue

Un cas classique auquel on s’intéressera ensuite consiste à classer un individu \(X\) selon une certaine mesure continue \(T(X)\), en définissant un threshold \(t\) et en fixant :

\[ D_X = \delta_{T(X) > t} \]

Les individus diagnostiqués positifs sont alors tous ceux dont la mesure \(T(X)\) est supérieure au threshold \(t\).

Reste ensuite à savoir si cette valeur sépare bien nos individus malades de nos individus non-malades. Ce sera l’objet de la troisième section.

Comparaison à un gold standard

Le gold standard donne la vérité

On dispose d’un test diagnostic appelé gold standard qui classe nos individus en deux premières catégories qui vont constituer notre référence :

• Malades, notés \(M=1\).
• Non-Malades, notés \(M=0\).

Une équipe de recherche biomédicale met au point un nouveau test diagnostic, et souhaite connaître les caractéristiques de ce test. Ce test classe nos individus en deux catégories :

• Diagnostiqués positifs, notés \(D=1\).
• Diagnostiqués négatifs, notés \(D=0\).

Questions et définitions classiques

On peut se poser tout un tas de questions à propos de la performance de notre test diagnostic, Notamment par exemple les plus courantes :

Quelle est la probabilité qu’un patient soit diagnostiqué positif \(D=1\) sachant qu’il est malade \(M=1\) ?

C’est ce qu’on appelle la sensibilité du test : \(\mathbb{P}(D = 1 | M=1)\).

Quelle est la probabilité qu’un patient soit diagnostiqué négatif \(D=0\) sachant qu’il est non-malade \(M=0\) ?

C’est ce qu’on appelle la spécificité du test : \(\mathbb{P}(D = 0 | M=0)\).

Etant donné qu’un patient est diagnostiqué positif \(D = 1\), quelle est la probabilité qu’il soit malade ?

C’est ce qu’on appelle la valeur prédictive positive : \(\mathbb{P}(M=1 | D=1)\).

Etant donné qu’un patient est diagnostiqué négatif \(D = 0\), quelle est la probabilité qu’il ne soit pas malade ?

C’est ce qu’on appelle la valeur prédictive négative : \(\mathbb{P}(M=0 | D=0)\).

Quelle est la probabilité qu’un patient soit classé correctement ?

C’est l’exactitude du test (on parle en général en anglais en utilisant le terme accuracy).

Comparaison en population générale

Pour nous fixer les idées sur un scénario particulier, nous allons considérer le scénario suivant :

• un échantillon de taille totale \(n = 1e3\) patients.
• chaque patient a une probabilité d’être malade \(\mathbb{P}(M=1) = 0.1\), c’est ce qu’on appelle la prévalence.
• sachant qu’il est malade, il a une probabilité d’être diagnostiqué positif de \(\mathbb{P}(D = 1 | M = 1) = 0.8\) (sensibilité).
• sachant qu’il est non-malade, il a une probabilité d’être diagnostiqué positif de \(\mathbb{P}(D = 1 | M = 0) = 0.1\) (1 - spécificité).

On peut simuler une réalisation de ce scénario de comparaison de test diagnostic :

n <- 1e3

p_M <- 0.1
preval_th <- p_M

p_D_given_M <- 0.8
sens_th <- p_D_given_M

p_D_given_NM <- 0.1
spec_th <- 1 - p_D_given_NM

# Simulation des malades et non-malades
n_M <- rbinom(1, n, p_M)
n_NM <- n - n_M

# Puis des diagnostics positifs et négatifs parmi les malades et non malades
n_D_M <- rbinom(1, n_M, p_D_given_M)
n_ND_M <- n_M - n_D_M
n_D_NM <- rbinom(1, n_NM, p_D_given_NM)
n_ND_NM <- n_NM - n_D_NM

# Organisation dans un tableau de données
lM <- "Malade\n(M=1)"
lNM <- "Non-Malade\n(M=0)"
lD <- "Diagnostiqué positif\n(D=1)"
lND <- "Diagnostiqué négatif\n(D=0)"

d <- data.frame(maladie = factor(c(lM, lM, lNM, lNM), levels=c(lNM, lM)),
                diagnostic = as.factor(c(lD, lND, lD, lND)),
                n = c(n_D_M, n_ND_M, n_D_NM, n_ND_NM))

# Plot des valeurs
d %>% ggplot(aes(x = maladie, y = diagnostic)) +
    geom_raster(aes(fill=n)) +
    geom_text(aes(label=n)) +
    scale_x_discrete(position = "top") +
    theme_bw()

L’avantage de faire notre comparaison en population générale est que le nombre de malades et non-malades est directement représentatif de la prévalence.

Toutes les quantités définies ci-dessus s’estiment donc de façon très naturelle en prenant en compte les bonnes cases du tableau :

sens_estim <- n_D_M / (n_D_M + n_ND_M)
spec_estim <- n_ND_NM / (n_ND_NM + n_D_NM)
vpp_estim <- n_D_M / (n_D_M + n_D_NM)
vpn_estim <- n_ND_NM / (n_ND_NM + n_ND_M)
acc_estim <- (n_D_M + n_ND_NM) / n

Sans surprise, les estimations ponctuelles de sensibilité et spécificité sont proches des valeurs fixées dans la simulation, avec respectivement 0.8278689 et 0.9077449, à comparer à 0.8 et 0.9.

La valeur prédictive positive (VPP) est estimée à 0.5549451, la valeur prédictive négative (VPN) est estimée à 0.9743276, ce qu’on peut comparer aux valeurs théoriques qu’on peut obtenir après un rapide calcul. Pour la VPP théorique, on a :

\[\begin{align*} VPP &= \mathbb{P}(M = 1 | D = 1) \\ &= \frac{ \mathbb{P}(M = 1, D = 1) }{ \mathbb{P}(D = 1) }\\ &= \frac{ \mathbb{P}(D = 1 | M = 1) \times \mathbb{P}(M = 1) }{ \mathbb{P}(D = 1 | M = 1) \times \mathbb{P}(M=1) + \mathbb{P}(D = 1 | M = 0) \times \mathbb{P}(M=0) }\\ &= \frac{ \text{sensibilité} \times \text{prévalence} }{ \text{sensibilité} \times \text{prévalence} + (1 - \text{spécificité}) \times (1 - \text{prévalence}) } \end{align*}\]

Numériquement, on obtient 0.4705882.

Pour la VPN théorique, on a :

\[\begin{align*} VPN &= \mathbb{P}(M = 0 | D = 0) \\ &= \frac{ \mathbb{P}(M = 0, D = 0) }{ \mathbb{P}(D = 0) }\\ &= \frac{ \mathbb{P}(D = 0 | M = 0) \times \mathbb{P}(M = 0) }{ \mathbb{P}(D = 0 | M = 0) \times \mathbb{P}(M=0) + \mathbb{P}(D = 0 | M = 1) \times \mathbb{P}(M=1) }\\ &= \frac{ \text{spécificité} \times (1 - \text{prévalence}) }{ \text{spécificité} \times (1 - \text{prévalence}) + (1 - \text{sensibilité}) \times \text{prévalence} } \end{align*}\]

Numériquement, avec les valeurs prises dans les simulations, on obtient 0.9759036.

Enfin, l’accuracy est estimée à 0.898. La valeur théorique, quant à elle, peut être calculée à partir des valeurs choisies dans la simulation :

\[\begin{align*} ACC &= \mathbb{P}((D = 1, M = 1) \cup (D = 0, M = 0))\\ &= \mathbb{P}(D = 1, M = 1) + \mathbb{P}(D = 0, M = 0)\\ &= \mathbb{P}(D = 1 | M = 1) \times \mathbb{P}(M = 1) + \mathbb{P}(D = 0 | M = 0) * \mathbb{P}(M = 0) \\ &= \text{sensibilité} \times \text{prévalence} + \text{spécificité} \times (1 - \text{prévalence}) \end{align*}\]

Numériquement, on obtient 0.89.

Comparaison dans une population controlant la maladie

Lorsque la prévalence d’une maladie est très faible, il peut être pénible et/ou trop coûteux de faire la comparaison sur un échantillon représentatif de la population générale.

Plus précisément, on peut toujours le faire, mais le nombre de malades dans notre échantillon risque d’être très faible, et la variance associée à notre estimateur de sensibilité sera très grande.

Dans ce cas, le design de l’étude peut être modifié, pour effectuer la comparaison du test diagnostic vs. gold standard sur un nombre fixé de patients malades et non-malades.

On peut à nouveau se fixer les idées sur un scénario simulé :

• un échantillon de 60 malades et 80 non-malades.
• sachant qu’il est malade, un patient a une probabilité d’être diagnostiqué positif de \(\mathbb{P}(D = 1 | M = 1) = 0.8\) (sensibilité).
• sachant qu’il est non-malade, un patient a une probabilité d’être diagnostiqué positif de \(\mathbb{P}(D = 1 | M = 0) = 0.1\) (1 - spécificité).

On peut simuler une réalisation de ce scénario de comparaison de test diagnostic :

p_D_given_M <- 0.8
sens_th <- p_D_given_M

p_D_given_NM <- 0.1
spec_th <- 1 - p_D_given_NM

# Simulation des malades et non-malades
n_M <- 60
n_NM <- 80

# Puis des diagnostics positifs et négatifs parmi les malades et non malades
n_D_M <- rbinom(1, n_M, p_D_given_M)
n_ND_M <- n_M - n_D_M
n_D_NM <- rbinom(1, n_NM, p_D_given_NM)
n_ND_NM <- n_NM - n_D_NM

# Organisation dans un tableau de données
lM <- "Malade\n(M=1)"
lNM <- "Non-Malade\n(M=0)"
lD <- "Diagnostiqué positif\n(D=1)"
lND <- "Diagnostiqué négatif\n(D=0)"

d <- data.frame(maladie = factor(c(lM, lM, lNM, lNM), levels=c(lNM, lM)),
                diagnostic = as.factor(c(lD, lND, lD, lND)),
                n = c(n_D_M, n_ND_M, n_D_NM, n_ND_NM))

# Plot des valeurs
d %>% ggplot(aes(x = maladie, y = diagnostic)) +
    geom_raster(aes(fill=n)) +
    geom_text(aes(label=n)) +
    scale_x_discrete(position = "top") +
    theme_bw()

A présent, que peut-on faire pour analyser de telles données ?

Comme précédemment, on peut toujours s’intéresser à la sensibilité et à la spécificité, puisqu’il s’agit de quantités qui peuvent directement être estimées parmi nos malades et nos non-malades.

sens_estim <- n_D_M / (n_D_M + n_ND_M)
spec_estim <- n_ND_NM / (n_ND_NM + n_D_NM)

Sans surprise, les estimations ponctuelles de sensibilité et spécificité sont proches des valeurs fixées dans la simulation, avec respectivement 0.75 et 0.9125, à comparer à 0.8 et 0.9.

En revanche, impossible de calculer les autres quantités d’intérêt, telles VPP et VPN, qui nécessitent de contrôler le nombre de patients diagnostiqués positifs et négatifs. Si on souhaite calculer ces quantités, on est contraint de reposer sur une information externe à propos de la prévalence de la maladie, de façon à pouvoir utiliser les formules dérivées précédemment.

De même, le calcul de l’accuracy doit reposer sur une information de prévalence externe.

Comparaison dans une population contrôlant les diagnostiqués

Le cas symétrique consiste à effectuer la même étude, mais cette fois sur un échantillon constitué d’un nombre fixé de patients diagnostiqués positifs \(D=1\) et de patients diagnostiqués négatifs \(D=0\).

Un tel cas design d’étude semble toutefois bien moins réaliste que les deux précédents. Il pourrait se rencontrer si le nouveau test diagnostic est réalisé en première intention, et que le gold standard, plus coûteux et long à réaliser, n’est réalisé que plus tard.

Dans un tel cas, on ne pourrait estimer que VPP et VPN à partir des données. Il nous faudrait alors une information de prévalence externe pour pouvoir en déduire quelque chose sur la sensibilité et la spécificité du test diagnostic.

Threshold sur mesure continue

Contexte

Intéressons-nous à présent au cas d’un test diagnostic reposant sur une mesure continue \(T(X)\) sur nos individus \(X\), et sur l’établissement d’un threshold \(t\) permettant de classer en \(D=1\) ou \(D=0\).

Le test est alors simplement \(D = \delta_{T(X) > t}\).

Supposons à présent qu’on possède les données de \(T(X)\) pour un ensemble d’individus, pour lesquels on a également accès au classement en \(M=1\) ou \(M=0\) donné par un gold standard.

LA question centrale d’établissement d’un nouveau test adéquat est donc :

Quelle valeur de threshold optimale \(t\) faut-il prendre ?

Et la question subsidiaire, extrêmement importante :

Que signifie l’ “optimalité”, dont il est question dans la question précédente ?

Scenario simulé

Ici encore, on va se fixer les idées en simulant un jeu de données fictif selon un scénario bien choisi :

• la prévalence, \(\mathbb{P}(M = 1)\) est de 0.1.
• Sachant qu’un individu est malade, \(T(X) \sim \mathcal{N}(\mu_M, \sigma_M)\).
• Sachant qu’un individu est non-malade, \(T(X) \sim \mathcal{N}(\mu_{\bar{M}}, \sigma_\bar{M})\).

n <- 1e3
preval <- 0.1
mu_NM <- 4
mu_M <- 6
sigma_M <- 1
sigma_NM <- 2

d <- data.frame(id = seq(1, n),
                 maladie = rbinom(n, size=1, prob=preval)) %>%
    mutate(mesure = rnorm(n, 
                          mean=ifelse(maladie==1, mu_M, mu_NM),
                          sd=ifelse(maladie==1, sigma_M, sigma_NM) ))

d %>% ggplot(aes(x = mesure)) +
    geom_histogram() +
    theme_bw()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Sensibilité, spécificité, accuracy

Mettons qu’on fixe une valeur de threshold \(t\). On dispose alors d’un test diagnostic, qu’on peut comparer à notre gold standard comme décrit dans la section précédente. Par exemple, pour \(t = 5\), on obtient :

get_contingency_table <- function(t, data){
    table <- data %>% mutate(diagnostic = ifelse(mesure > t, 1, 0)) %>%
        group_by(maladie, diagnostic) %>%
        summarise(nb = n()) %>% ungroup()

    return (table)
}

# Plot des valeurs
table <- get_contingency_table(t=5, data=d)

## `summarise()` has grouped output by 'maladie'. You can override using the
## `.groups` argument.

table %>% mutate(maladie = ifelse(maladie==1, "malade", "non-malade"),
             diagnostic = ifelse(diagnostic==1, "diagnostic positif", "diagnostic négatif")) %>%
    ggplot(aes(x = maladie, y = diagnostic)) +
    geom_raster(aes(fill=nb)) +
    geom_text(aes(label=nb)) +
    scale_x_discrete(position = "top") +
    theme_bw()

On peut alors en déduire la sensibilité et la spécificité du test :

get_sens_spec_acc <- function(table){
    n_M_D <- table %>% filter(maladie==1, diagnostic==1) %>% pull(nb)
    n_M_D <- ifelse(length(n_M_D) == 0, 0, n_M_D)

    n_M_ND <- table %>% filter(maladie==1, diagnostic==0) %>% pull(nb)
    n_M_ND <- ifelse(length(n_M_ND) == 0, 0, n_M_ND)

    n_NM_ND <- table %>% filter(maladie==0, diagnostic==0) %>% pull(nb)
    n_NM_ND <- ifelse(length(n_NM_ND) == 0, 0, n_NM_ND)

    n_NM_D <- table %>% filter(maladie==0, diagnostic==1) %>% pull(nb)
    n_NM_D <- ifelse(length(n_NM_D) == 0, 0, n_NM_D)

    sens_estim <- n_M_D / (n_M_D + n_M_ND)
    spec_estim <- n_NM_ND / (n_NM_ND + n_NM_D)
    acc_estim <- (n_M_D + n_NM_ND) / (n_M_D + n_M_ND + n_NM_ND + n_NM_D)

    return (list(sens = sens_estim, spec = spec_estim, acc = acc_estim))
}
estim <- get_sens_spec_acc(table)

On obtient ici des valeurs numériques respectivement à 0.8658537 pour la sensibilité, 0.6797386 pour la spécificité, et 0.695 pour l’accuracy.

Choix du seuil

On peut ensuite aller plus loin, et calculer ces valeurs pour toutes les valeurs de threshold qui pourraient nous intéresser, par exemple toutes les valeurs entre -2 et 12, dans l’espoir de guider le choix de \(t\).

t_list <- seq(from=-2, to=12, by=0.2)
sens_list <- c()
spec_list <- c()
acc_list <- c()

for (t in t_list){
    table <- get_contingency_table(t, d)
    estim <- get_sens_spec_acc(table)
    
    sens_list <- c(sens_list, estim$sens)
    spec_list <- c(spec_list, estim$spec)
    acc_list <- c(acc_list, estim$acc)
}

d_charac <- data.frame(t = rep(t_list, 3),
                          val = c(sens_list, spec_list, acc_list),
                          type = c(rep("sensibilité", length(t_list)),
                                   rep("spécificité", length(t_list)),
                                   rep("accuracy", length(t_list)) ) )

d_charac %>% ggplot(aes(x=t, y=val, col=type)) +
    theme_bw() +
    geom_point() +
    geom_line()

On pourrait donc penser à choisir :

• le seuil qui maximise la sensibilité,
• le seuil qui maximise la spécificité,
• le seuil qui maximise l’accuracy.

Selon ce qui fait le plus sens dans une situation donnée.

Courbe ROC

Une autre représentation graphique très appréciée consiste à représenter cette fois-ci la sensibilité en fonction de (1 - spécificité).

Cette courbe s’appelle la courbe ROC, pour Receiver Operating Characteristic, ou plus simplement, dans le domaine des stats médicales, courbe sensibilité/spécificité, et la voici dans le cas qui nous intéresse :

d_roc <- data.frame(sensibility = sens_list,
                    specificity = spec_list,
                    one_minus_specificity = 1 - spec_list,
                    threshold = t_list)

d_roc %>% ggplot(aes(y = sensibility, x = one_minus_specificity)) +
    geom_point() +
    geom_line() +
    xlab("1 - spécificité") +
    ylab("sensibilité") +
    geom_abline(intercept = 0, slope = 1, col="gray") +
    theme_bw()

Chaque point de cette courbe correspond à une valeur de seuil différente. Elle permet de se faire une bonne idée du trade-off qui existe entre la sensbilité et la spécificité :

• en haut à droite on se place au point qui maximise la sensibilité, au prix d’une spécificité faible,
• en bas à gauche, on se place au point qui maximise la spécificité au prix d’une sensibilité faible,
• entre les deux, sur la courbe, on peut choisir le point qui nous semble être le meilleur compromis entre les deux,
• en haut à gauche, au point de coordonnées (0,1), il y aurait un test diagnostic parfait, qui ne ferait aucun faux positif ni aucun faux négatif.

Notons que la droite grise superposée au graphe correspond à ce qu’on obtient si on diagnostique nos individus au hasard, en jouant à pile ou face. La distance à cette courbe permet d’apprécier à quel point la valeur sur laquelle on se base pour établir le diagnostic est distribuée différemment entre les malades et les non-malades.

A présent qu’on a bien mis les mains dans le cambouis pour comprendre comment ça fonctionne, place à la magie de R et des multiples packages qui existent ! Evidemment, toutes ces opérations pénibles à faire peuvent être réalisées en trois ligne de code à l’aide du package ROCR (par exemple).

library(ROCR)
# on créé un objet de "prédiction" à partir de notre mesure et du gold standard
pred <- prediction(d$mesure, d$maladie)
# on calcule les indicateurs de "tpr" : true positive rate = sensibilité
# et "fpr" : false positive rate = 1 - specificité,
# pour toutes les valeurs de threshold possibles
perf <- performance(pred, "tpr", "fpr")
plot(perf)

Notons que, si on souhaite tout de même utiliser ggplot pour améliorer ce graphique, la meilleure option consiste à regarder le contenu de l’object perf et à en utiliser des morceaux bien choisis. Par exemple :

str(perf)

## Formal class 'performance' [package "ROCR"] with 6 slots
##   ..@ x.name      : chr "False positive rate"
##   ..@ y.name      : chr "True positive rate"
##   ..@ alpha.name  : chr "Cutoff"
##   ..@ x.values    :List of 1
##   .. ..$ : num [1:1001] 0 0.00109 0.00218 0.00327 0.00436 ...
##   ..@ y.values    :List of 1
##   .. ..$ : num [1:1001] 0 0 0 0 0 0 0 0 0 0 ...
##   ..@ alpha.values:List of 1
##   .. ..$ : num [1:1001] Inf 10.5 10.29 10.15 9.75 ...

data.frame(sensibility = perf@y.values[[1]], one_minus_specificity = perf@x.values[[1]]) %>%
    ggplot(aes(y = sensibility, x = one_minus_specificity)) +
    geom_line() +
    xlab("1 - spécificité") +
    ylab("sensibilité") +
    geom_abline(intercept = 0, slope = 1, col="gray") +
    theme_bw()

Un autre package alternatif, pROC, est également disponible. Voici quelques lignes pour débuter avec :

library(pROC)
roc(maladie ~ mesure, d) %>% plot.roc()

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

De noombreuses options sont disponibles pour modifier le comportement de la fonction roc et de la fonction d’affichage plot.roc.

ROC AUC

Imaginons qu’on souhaite, à présent, comparer deux (ou plus) tests diagnostics. Si on sait déjà qu’on souhaite prioriser, par exemple, l’accuracy, la question est relativement simple à résoudre :

• on prend le seuil maximisant l’accuracy dans le premier test,
• on prend le seuil maximisant l’accuracy dans le second test,
• et on prend le test (et le seuil) maximisant l’accuracy entre les deux tests.

En revanche, lorsqu’on ne sait pas bien choisir si on souhaite prioriser sensibilité, spécificité ou accuracy, on peut s’appuyer sur une nouvelle quantité intéressante basée sur la courbe ROC : la ROC AUC, pour ROC Area Under The Curve.

L’idée est très simple : plus l’aire sous la courbe est importante, meilleur sera considéré notre test diagnostic. Simulons un deuxième jeu de données pour visualiser ce que ça donne, avec une autre mesure permettant de diagnostiquer notre maladie :

n <- 1e3
preval <- 0.1

mu1_NM <- 4
mu1_M <- 6
sigma1_M <- 1
sigma1_NM <- 2

mu2_NM <- 3
mu2_M <- 6
sigma2_M <- 3
sigma2_NM <- 3

d <- data.frame(id = seq(1, n),
                 maladie = rbinom(n, size=1, prob=preval)) %>%
    mutate(mesure1 = rnorm(n, 
                          mean=ifelse(maladie==1, mu1_M, mu1_NM),
                          sd=ifelse(maladie==1, sigma1_M, sigma1_NM)),
           mesure2 = rnorm(n, 
                          mean=ifelse(maladie==1, mu2_M, mu2_NM),
                          sd=ifelse(maladie==1, sigma2_M, sigma2_NM) ))

On va utiliser le package pROC pour nous calculer la courbe ROC, ainsi que l’aire sous la courbe :

r1 <- roc(maladie ~ mesure1, d)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

r2 <- roc(maladie ~ mesure2, d)

## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

auc1 <- r1$auc
auc2 <- r2$auc

data.frame(sens = c(r1$sensitivities, r2$sensitivities),
           one_minus_spec = c(1-r1$specificities, 1-r2$specificities),
           test = c(rep("1", n+1), rep("2", n+1)) ) %>%
    ggplot(aes(x=one_minus_spec, y=sens, col=test)) +
    geom_line() +
    theme_bw() +
    ylab("sensibilité") +
    xlab("1 - spécificité")

L’aire sous la courbe vaut 0.8060096 pour le test diagnostic numéro 1, et 0.8002912 pour le test diagnostic numéro 2. Sur la base de la ROC AUC, on choisirait donc préférentiellement le test diagnostic numéro 1. Notons qu’on choisirait également le test diagnostic numéro 1 si on souhaite obtenir un test diagnostic ayant une forte sensibilité. En revanche, on choisirait préférentiellement le test diagnostic numéro 2 si on souhaitait obtenir un test diagnostic à plus forte spécificité.

Conclusion

On peut identifier un certain nombre de take-home messages sur ce sujet :

• la sensibilité est la probabilité d’être diagnostiqué positif sachant qu’on est malade. \[ \mathbb{P}(D = 1 | M = 1) \]
• la spécificité est la probabilité d’être diagnostiqué négatif sachant qu’on est non-malade. \[ \mathbb{P}(D = 0 | M = 0) \]
• la valeur prédictive positive (VPP) est la probabilité d’être malade sachant qu’on est diagnostiqué positif. \[ \mathbb{P}(M = 1 | D = 1) \]
• la valeur prédictive négative (VPN) est la probabilité d’être non-malade sachant qu’on est diagnostiqué négatif. \[ \mathbb{P}(M = 0 | D = 0) \]
• dans un échantillon représentatif de la population générale, la comparaison d’un test diagnostic au gold standard permet d’estimer la sensibilité, la spécificité, la VPP et la VPN à partir de la table de contingence.
• dans un échantillon contenant un nombre fixé de malades et non-malades, on ne peut estimer que la sensibilité et la spécificité. Le calcul de la VPP et VPN requièrent d’avoir accès à la prévalence de la maladie.
• si on souhaite créer un test diagnostic à partir d’une mesure continue, on peut calculer la sensibilité et spécificité pour toutes les valeurs de seuil possibles, et construire une courbe ROC. On peut alors choisir le seuil que l’on souhaite, maximisant plutôt la sensibilité ou la spécificité.
• l’aire sous la courbe ROC (ROC AUC) permet de comparer des tests diagnostics lorsqu’on ne sait pas si on doit prioriser une bonne sensibilité ou une bonne spécificité.