Exemple t-test

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

Méthode classique : t-test

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.

Méthode bayésienne

• 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

  • En effet, les queues de distribution d'une loi de student sont plus épaisses que celle d'une normale, et s'accommodent donc mieux de la présence éventuelle d'outliers
  • Une distribution de student est paramétrée par $\nu$ le degré de liberté qui varie de 1 à l'infini.
  • si $\nu=\infty$, la student devient une loi normale exacte.
  • si $\nu > 30$, on est déjà très proche d'une normale

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 :

• pour $\mu$ une normale $N(\bar y, 1000 * sd(y))$
• pour $\sigma$ une loi uniforme $U(\frac{sd(y)}{1000},sd(y) \times 1000)$
• pour $\nu$, moyenne 30 et écart type 30, $\nu >1$ par contrainte

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.

exemple 1 sans outliers dixit shapiro

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 :

• Rhat statistique de Gelman rubin (scale reduction) converge bien vers 1
• n.eff est bien supérieur à 10000
• on a donc une bonne convergence des chaînes MCMC

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.

Exemple 2 avec outliers

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.

Exemple 3 si on a un a priori sur la moyenne

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

• le ROPE n'est pas entièrement contenu dans le HDI (ni complètement extérieur), on doit conclure à la nécessité de relever davantage d'observations. Ou bien revoir la ROPE, qui est fixée arbitrairement.