Le padbolome est un cancer rare touchant 0,1% de la population. Donald fait un test de dépistage qui se révèle positif (= indique “malade”). Pas de chance ! Mais au fond, quelle est la probabilité que Donald soit réellement malade ?
Fiabilité du test:
#données du problème
p_malade=0.001 #prior
p_fauxpos=0.03
p_fauxneg=0.05
On simule une population de N individus, pour lesquels l'occurrence de la maladie est tirée aléatoirement avec la probabilité p_malade.
Ensuite on simule le résultats d'une campagne de dépistage pour cette population. Pour chaque individu, on tire aléatoirement le résultat du test en fonction de si la personne est malade ou non, et des probabilités de faux positifs et faux négatifs données.
# generation des malades et non-malades
# 0 code le malade, 1 code l'individu sain
genere_population=function(n=10000) {
pop<-ifelse(runif(N,0,1)<p_malade,0,1)
}
# on fait les tests de depistage :
# 0 si test positif, 1 si test negatif (conforme au codage maladie)
genere_test=function(pop){
#generer d'abord des random numbers
#utilisés après pour choisir les résultats des tests
testgen=runif(length(pop),0,1)
# génération des résultats de test
test<-ifelse((pop==1&testgen>p_fauxpos) # individu sain et test negatif
| (pop==0&testgen<p_fauxneg),1,0) # individu malade et test negatif
}
Maintenant on utilise ces fonctions de simulation pour évaluer la probabilité que Donald soit malade sachant que son test est positif.
risque_malade_si_positif=function(n){
pop=genere_population(n)
test=genere_test(pop)
nbpositifs=sum(test==0)
nbpositifmalade=sum(!(test|pop)) # éléments tels que test=0 et pop=0
nbpositifmalade/nbpositifs
}
Maintenant, on fait l'expérience de simulation (= on fait tourner le code précédent). Évidemment, la précision obtenue sur l'estimation de cette probabilité dépend de la taille de la population simulée. Si l'on ne simule que 100 personnes, il y a fort à parier que l'on ne verra aucun malade…
N=100000 #taille de la population testée
risque_malade_si_positif(N)
## [1] 0.03375248
Bon, 1 estimation cela ne donne pas trop d'indications sur la précision. Si on répétait l'expérience?
estimator=function(size=10000,nbrep=10){
estimates<-c()
for (i in 1:nbrep){
estimates[i]<-risque_malade_si_positif(size)
}
estimates
}
On peut maintenant lancer une “campagne d'expériences”
nb_replications=100
estimates=estimator(N,nb_replications)
plot(estimates,type="p",ylab="Estimations de la proba d'être malade si le test positif")
summary(estimates)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02156 0.02811 0.03038 0.03025 0.03231 0.04036
Question subsidiaire : est-il équivalent de simuler 100 x 10 000 personnes et 10 x 100 000 personnes?
N=1000000
nb_replications=10
estimates=estimator(N,nb_replications)
plot(estimates,type="p",ylab="Estimations de la proba d'être malade si le test positif")
summary(estimates)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02924 0.02995 0.03036 0.03049 0.03112 0.03163