Comment réaliser l'équivalent du t-test en bayésien ? La méthode permet d'estimer la moyenne de l'échantillon, et fournit également une estimation de la normalité de la série. Les outliers sont automatiquement pris en compte.
aperçu de la séquenceUn médecin fait l'hypothèse qu'il voit en moyenne ses patients 5 fois par an.
Pour évaluer la validité de cette hypothèse, il sélectionne 10 patients au hasard et note le nombre de leurs visites.
nvisites = c(9,10,8,4,8,3,0,10,15,9)
Le t-test suppose que les données sont distribuées selon une loi normale, avec $\sigma$ inconnu.
# pas d'évidence de non normalité
# mais il serait surprenant de rejeter H0 avec si peu d'observations
shapiro.test(nvisites)
##
## Shapiro-Wilk normality test
##
## data: nvisites
## W = 0.94032, p-value = 0.5565
# remarquer au passage que des outliers sont détectés
shapiro.test(c(nvisites, 30))
##
## Shapiro-Wilk normality test
##
## data: c(nvisites, 30)
## W = 0.81781, p-value = 0.01614
#test non directionnel
t.test(nvisites,mu=5)
##
## One Sample t-test
##
## data: nvisites
## t = 1.9355, df = 9, p-value = 0.08492
## alternative hypothesis: true mean is not equal to 5
## 95 percent confidence interval:
## 4.561253 10.638747
## sample estimates:
## mean of x
## 7.6
On accepte $H_0$ au seuil de 5 pourcent pour l'alternative $H_1 \colon \theta \neq 5$
# test unilateral à droite
t.test(nvisites,mu=5, alternative = "greater")
##
## One Sample t-test
##
## data: nvisites
## t = 1.9355, df = 9, p-value = 0.04246
## alternative hypothesis: true mean is greater than 5
## 95 percent confidence interval:
## 5.137587 Inf
## sample estimates:
## mean of x
## 7.6
On rejette $H_0$ au seuil de 5 pourcent pour l'alternative $H_1 \colon \theta > 5$
Le test unidirectionnel revient à encoder une information a priori (comme en bayésien) : le médecin pense que ça ne peut pas être inférieur à 5. Et dans ce cas le test est plus puissant (disposant de davantage d'informations), et détecte un effet.
La test bayésien ne nécessite pas que les données soient distribuées selon une loi normale
La procédure implémentées dans BEST suppose que les données sont distribuées selon un loi de student, afin de modéliser la présence éventuelle d'outliers
Dans la modélisation bayésienne, on modélise également ce paramètre $\nu$, ce qui donne une information sur le caractère plus ou moins éloigné de la normalité de la série de données.
Les distributions a priori par défaut des paramètres de la student à estimer par la fonction BESTmcmc sont :
La student qui modélise les données est donc proche d'une normale a priori, mais $\nu$ pourra a posteriori se positionner plus bas que 30 si la série contient des outliers.
library(BEST)
## Loading required package: HDInterval
BESTout = BESTmcmc(nvisites, parallel = TRUE)
## Waiting for parallel processing to complete...
## done.
# indiquer la région opérationelle d'équivalence pour représentation sur le graphique
plot(BESTout, compVal=5, ROPE=c(4.5,5.5))
summary(BESTout)
## mean median mode HDI% HDIlo HDIup compVal %>compVal
## mu 7.67 7.68 7.79 95 4.473 10.90 0 100
## sigma 4.71 4.45 4.11 95 2.248 7.77
## nu 31.38 22.60 6.82 95 1.000 89.39
## log10nu 1.31 1.35 1.52 95 0.428 2.08
## effSz 1.79 1.72 1.63 95 0.559 3.04 0 100
print(BESTout)
## MCMC fit results for BEST analysis:
## 100002 simulations saved.
## mean sd median HDIlo HDIup Rhat n.eff
## mu 7.668 1.615 7.684 4.473 10.897 1.000 56010
## nu 31.382 28.877 22.605 1.000 89.388 1.000 19127
## sigma 4.712 1.498 4.450 2.248 7.774 1.001 25892
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## 'Rhat' is the potential scale reduction factor (at convergence, Rhat=1).
## 'n.eff' is a crude measure of effective sample size.
Lecture des résultats :
On accepte l'hypothèse $H_0$ pour un ROPE qui va de 4.5 à 5,5, car le ROPE est entièrement contenu dans la région hautement probable HDI à 95 pourcent.
La taille de l'effet (effsz) est en moyenne de 1.80 (sous-entendu de plus que 5). Ce paramètre s'interprète surtout lorsqu'on compare la moyenne de 2 populations.
On remarque que $\nu$ est resté autour de 30 en moyenne, même si son HDI est très large. C'est compatible avec le test de shapiro. Le faible nombre d'observation explique le HDI très large.
library(BEST)
BESTout = BESTmcmc(c(nvisites,30), parallel = TRUE)
## Waiting for parallel processing to complete...done.
# indiquer la région opérationelle d'équivalence pour représentation sur le graphique
plot(BESTout, compVal=5, ROPE=c(4.5,5.5))
summary(BESTout)
## mean median mode HDI% HDIlo HDIup compVal %>compVal
## mu 8.92 8.81 8.72 95 4.133 14.08 0 99.9
## sigma 7.28 7.02 6.82 95 2.153 12.55
## nu 22.31 12.80 3.17 95 1.000 74.80
## log10nu 1.08 1.11 1.30 95 0.117 1.95
## effSz 1.41 1.31 1.24 95 0.250 2.74 0 99.9
print(BESTout)
## MCMC fit results for BEST analysis:
## 100002 simulations saved.
## mean sd median HDIlo HDIup Rhat n.eff
## mu 8.916 2.469 8.807 4.133 14.08 1 46130
## nu 22.311 25.706 12.805 1.000 74.80 1 10542
## sigma 7.285 2.750 7.018 2.153 12.55 1 16704
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## 'Rhat' is the potential scale reduction factor (at convergence, Rhat=1).
## 'n.eff' is a crude measure of effective sample size.
On remarque que $\nu$ est trombé à 22 en moyenne, même si son HDI est très large. C'est encore compatible avec le test de shapiro. Le faible nombre d'observation explique le HDI très large.
Le médecin a toujours constaté que c'était environ 5 par le passé.
On peut encoder cet a priori dans la modélisation bayésienne.
# on suppose a priori que le nombre de visites moyen est 5
BESTout = BESTmcmc(nvisites, parallel = TRUE, priors=list(muM=5))
## Waiting for parallel processing to complete...done.
plot(BESTout, compVal=5, ROPE=c(4.5,5.5))
summary(BESTout)
## mean median mode HDI% HDIlo HDIup compVal %>compVal
## mu 7.66 7.68 7.69 95 4.440 10.87 0 100
## sigma 4.69 4.45 3.99 95 2.283 7.66
## nu 32.64 23.50 7.12 95 0.538 93.43
## log10nu 1.33 1.37 1.48 95 0.443 2.11
## effSz 1.79 1.73 1.62 95 0.596 3.06 0 100
print(BESTout)
## MCMC fit results for BEST analysis:
## 100002 simulations saved.
## mean sd median HDIlo HDIup Rhat n.eff
## mu 7.660 1.616 7.678 4.4402 10.874 1 52978
## nu 32.642 30.186 23.503 0.5379 93.434 1 18362
## sigma 4.692 1.451 4.449 2.2832 7.663 1 27232
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## 'Rhat' is the potential scale reduction factor (at convergence, Rhat=1).
## 'n.eff' is a crude measure of effective sample size.
Résultat, la moyenne a posteriori est plus basse