Tag Archives: puissance

Niveau et puissance de test (paramétrique)

Lors du dernier cours de statistique, nous avons vu comme se construisent des tests, et comment on compare leur puissance. Considérons pour cela des observations \{x_1,\cdots,x_n\} tirées suivant une loi uniforme \mathcal{U}([0,\theta]). Et on se pose la question de savoir si l’hypothèse H_0:\theta\leq\theta_0 est rejetée, ou acceptée, avec \theta_0=1 par exemple. L’hypothèse alternative sera ici le complémentaire, autrement dit H_1:\theta>\theta_0.

Comme on l’a vu, la construction d’un test se fait en utilisant une statistique de test (ou plus simplement un estimateur du paramètre inconnu), avec \psi(\mathbf{x})=\mathbf{1}(\widehat{\theta}(\mathbf{x}))\leq c. On prend la décision D_1 (consistant à rejeter H_0) si \psi(\mathbf{x})=1, et la décision D_0 (consistant à ‘accepter’ H_0) si \psi(\mathbf{x})=0.

Le risque d’erreur de première espèce est basé sur \mathbb{E}[\psi(\mathbf{X})]. Si l’hypothèse H_0 est une hypothèse simple H_0:\theta=\theta_0, le risque d’erreur de première espèce est https://latex.codecogs.com/gif.latex?\alpha=\mathbb{E}[\psi(\mathbf{X})\vert%20X_1,\cdots,X_n\sim%20F_{\theta_0}], soit \alpha=\mathbb{P}[\widehat{\theta}(\mathbf{X})>c \vert X_1,\cdots,X_n\sim F_{\theta_0}] pour un seuil c. On va déterminer le seuil c de manière à ce que le risque de première espèce soit 5%.

Pour cela, il nous fait un estimateur \widehat{\theta} dont on connaîtrait la loi. Commençons par l’estimateur du maximum de vraisemblance. Rappelons que l’estimateur du maximum de vraisemblance est \widehat{\theta}=\max\{x_i\}. Nous avons vu – dans la section sur les statistiques d’ordre – que \mathbb{P}[\max\{X_i\}]=\left(\frac{x}{\theta}\right)^n

pour x\leq\theta. Aussi, pour un seuil c, \mathbb{P}[\widehat{\theta}>c]=1-\left(\frac{c}{\theta}\right)^n

pour c\leq\theta. Le risque de première espèce \alpha=\underset{\theta\leq\theta_0}{\max}\left\lbrace\mathbb{P}\big[\widehat{\theta}\geq c\vert X\sim F_{\theta}\big]\right\rbrace

soit \alpha=\underset{\theta\leq\theta_0}{\max}\left\lbrace1-\left(\frac{c}{\theta}\right)^n\right\rbrace

qui atteindra son maximum en \theta_0. Aussi \alpha=1-\left(\frac{c}{\theta_0}\right)^n

soit c=\theta_0\cdot[1-\alpha]^{1/n}. Les régions d’acceptation et de rejet sont les suivantes,

Continue reading Niveau et puissance de test (paramétrique)

La normalité ? avant tout une histoire de taille

Un billet rapide pour préciser le fond de ma pensée lorsque j’explique aux étudiants (qui travaillent sur leur devoir) que sur les bases de données avec beaucoup d’observations, il est très dur d’accepter l’hypothèse de normalité.

Pour ça, simulons des jeux de données: on a une seule variable explicative, un bruit gaussien, homoscédastique. Mais on va mal spécifier le modèle. En particulier, le vrai modèle ne sera pas linéaire. Voilà ce qu’on simule, et que l’on va chercher à ajuster.

n=50
X1=runif(n)
epsilon=rnorm(n)
Y=-2+8*X1-4*(X1<.2)+epsilon
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)

Autrement dit, les résidus (du modèle linéaire) ne devraient pas être gaussiens, car on n’ajuste pas le bon modèle. Maintenant, mon point est le suivant: avec peu de données, la plupart des tests de normalité ne vont pas rejeter l’hypothèse de normalité. Alors qu’avec beaucoup de données, l’hypothèse de normalité sera presque toujours rejetée.
Prenons par exemple le test Lillie, évoqué dans un vieux billet, qui est une adaptation du test de Kolmogorov-Smirnov quand on compare la distribution empirique avec la loi normale dont les paramètres ont été obtenues par maximum de vraisemblance. Sur le graphique ci-dessous, à gauche, on a en abscisse la taille de l’échantillon (qui va de 10 à 1000), et en ordonnée la p-value du test Lillie. Ou plus précisément, la courbe bleue est la région où s’est trouvée la p-value dans 90% des cas (la courbe supérieure étant le quantile à 95% de la p-value sur 10,000 échantillons simulés). La courbe rouge est la valeur moyenne de la p-value. A droite, on a – toujours en fonction de la taille de l’échantillon – la proportion des échantillons simulés pour lesquels la p-value était supérieure à 5%.

Autrement dit, plus la taille de l’échantillon est grande, plus on a de chance de rejeter (fort justement d’ailleurs) l’hypothèse de normalité. Et ce, uniquement sur un effet taille. Ou pour être plus rigoureux, sur les petits échantillons, on accepte beaucoup trop souvent (à tort cette fois) l’hypothèse de normalité.
On peut aussi regarder le test de Shapiro-Wilk

ou d’Anderson-Darling,

ou encore de Cramér-von Mises

Le code R pour faire ces tests est relativement simple,

X1=runif(n)
Y=-2+8*X1-4*(X1<.2)+rnorm(n)
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)
AD[n,i]=ad.test(residuals(reg))$p.value
SW[n,i]=shapiro.test(residuals(reg))$p.value
CV[n,i]=cvm.test(residuals(reg))$p.value
LL[n,i]=lillie.test(residuals(reg))$p.value

Pour ceux qui veulent aller plus loin, se cache derrière la notion de puissance du test…

Mal calculer les probabilités tue les probabilités

Suite à un tweet d’Olivier, j’ai découvert avec stupéfaction l’article paru dans Libé sur le risque nucléaire (ici).

On apprend que la probabilité d’accident majeur est d’environ 0.0003 0,0003 par an pour chaque réacteur. Or l’Europe possède 143 réacteurs et donc, sur 30 ans, la probabilité d’un incident majeur est

http://freakonometrics.hypotheses.org/files/2015/12/nucle-05.gif

et les auteurs de l’article concluent que “la probabilité d’occurrence d’un accident majeur sur ces parcs serait donc […] de plus de 100% pour l’Union européenne”

Etienne Ghys a déjà – brillamment – discuté ce chiffre avant hier (ici). Le soucis est que si http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif est la proba annuelle d’incident par an et par réacteur, en supposant les centrales indépendantes, la probabilité d’avoir un incident sur http://freakonometrics.hypotheses.org/files/2015/12/nucle-3.gif réacteurs sur une durée de http://freakonometrics.hypotheses.org/files/2015/12/nucle-04.gif années est

http://freakonometrics.hypotheses.org/files/2015/12/nucle-01.gif

Comme http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif est supposé petit, un développement limité donnerait

http://freakonometrics.hypotheses.org/files/2015/12/nucle-06.gif

Malheureusement, dans les développements limités de ce genre, s’il faut que http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif soit petit, encore faut-il que la puissance ne soit pas trop grande….
En fait, si on retient que pour un réacteur nucléaire à eau pressurisée (REP) tels ceux exploités en Europe de l’Ouest, le risque de fusion du cœur est estimé à 0.00005 par centrale et par an (en reprenant les estimateurs mentionnés sur wikipedia, ici), la vraie probabilité d’avoir un incident majeur en Europe serait plutôt de l’ordre de

http://freakonometrics.hypotheses.org/files/2015/12/nuclea-07.gif

(avec une probabilité 10 fois supérieure, Etienne arrivait à 72%). C’est marrant car cette erreur commise en croyant qu’un développement limité serait légitime quand on utilise des probabilité annuelle sur des horizons plus longs était évoqué ici il y a un an jour pour jour, sur des problèmes d’actualisation….

Tests statistiques et p value

Mercredi, nous avons abordé le concept de http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value en cours (avant de parler tout à l’heure de tests asymptotiques). Intuitivement, la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est la probabilité que la statistique de test prenne une valeur au moins aussi extrême que celle qui a été observée, si l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif est vraie.

  • Fisher vs. Neyman (et Pearson)

Dans l’approche de Fisher, on suppose que http://freakonometrics.blog.free.fr/public/perso2/testp-01.gif et on cherche à tester http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif. On dispose pour cela d’un échantillon http://freakonometrics.blog.free.fr/public/perso2/testp-02.gif.
On choisit une statistique http://freakonometrics.blog.free.fr/public/perso2/testp-03.gif telle que plus http://freakonometrics.blog.free.fr/public/perso2/testic-5.gif est grande, plus on s’éloigne de http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif (on calcule une distance à l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif, quelque chose comme ça). On calcule alors la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value comme

http://freakonometrics.blog.free.fr/public/perso2/testp-07.gif

et on rejette http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif si http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif est trop petite. La notion de “plus petite” est liée à la comparaison avec un seuil prédéfini, souvent 5% (cf plus loin pour une discussion).
C’est un peu heuristique, mais c’est comme cela que ça fonctionne.
Dans l’approche de Neyman-Pearson, on cherche à testerhttp://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif contre une hypothèse alternative, du genre http://freakonometrics.blog.free.fr/public/perso2/testp-11.gif. Si http://freakonometrics.blog.free.fr/public/perso2/testp-12.gif, on va construire une région de rejet de la forme

http://freakonometrics.blog.free.fr/public/perso2/testp-13.gif

où http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est choisi en fonction d’une probabilité d’erreur de première espèce souhaitée a priori, i.e.

http://freakonometrics.blog.free.fr/public/perso2/testp-15.gif

Par exemple, considérons un échantillon suivant une loi http://freakonometrics.blog.free.fr/public/perso2/testp-16.gif. On dispose de http://freakonometrics.blog.free.fr/public/perso2/testp-17.gif observations pour tester http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif. L’échantillon est le suivant

> set.seed(4)
> (X=rnorm(10))
 [1]  0.2167549 -0.5424926  0.8911446  0.5959806  1.6356180  0.6892754
 [7] -1.2812466 -0.2131445  1.8965399  1.7768632

La statistique proposée par Fisher (je ne discuterais pas le choix optimal de la statistique – si quelqu’un a des références, les commentaires sont ouvertes) est la moyenne normalisée,

http://freakonometrics.blog.free.fr/public/perso2/testp-21.gif

On sait que sous H_0, cette statistique suit une loi normale centrée réduite.  La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value se calcule alors simplement

> (T=sqrt(10)*mean(X))
[1] 1.791523
 
> 2*(1-pnorm(T,mean=0,sd=0)))
[1] 0.07320942

Graphiquement, cela se représente de la manière suivante.

Cette probabilité étant faible, on aurait tendance vouloir rejeter http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif.
Dans l’optique de Neyman Pearson, si on teste http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif contre http://freakonometrics.blog.free.fr/public/perso2/testp-22.gif, la région critique s’obtient à l’aide du rapport de vraisemblance

http://freakonometrics.blog.free.fr/public/perso2/testp-25.gif

ou encore

http://freakonometrics.blog.free.fr/public/perso2/testp-26.gif

ce qui se traduit par

http://freakonometrics.blog.free.fr/public/perso2/testp-27.gif

où http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est tel que

http://freakonometrics.blog.free.fr/public/perso2/testp-29.gif

Or sous H_0, la somme suit une loi normale centrée de variance http://freakonometrics.blog.free.fr/public/perso2/testp-30.gif. Ici, http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est donné par

> qnorm(1-alpha,mean=0,sd=sqrt(10))
[1] 5.201484

Or pour rappel, la somme vaut ici

> sum(X)
[1] 5.665293

On est dans la région de rejet: on rejette alors http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif.

  • p-value pour les tests usuels de R

La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est intéressante en pratique car elle permet à l’utilisateur de ne pas rentrer trop dans les détails de la construction du test, et de la loi de la statistique considérée.
Par exemple (comme celui abordé ici), pour un test d’égalité de proportion, sous R, la commande est

> prop.test(nx,n,.5)
 
1-sample proportions test with continuity correction
 
data:  nx out of n, null probability 0.5
X-squared = 0.05, df = 1, p-value = 0.823
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.3204804 0.7617145
sample estimates:
p
0.55

Et de manière générale, tous les tests sous R renvoient une http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value.

  • p ne signifie pas puissance

(mais probabilité). Dans un test simple,http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif contrehttp://freakonometrics.blog.free.fr/public/perso2/testh13.gif, la puissance est la probabilité de rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif sous l’hypothèse où http://freakonometrics.blog.free.fr/public/perso2/H1.gif est vraie (rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif avec raison si on veut). C’est donc une fonction de  http://freakonometrics.blog.free.fr/public/perso2/testh14.gif.
Pour un test multiple, avec une hypothèse alternative de la forme http://freakonometrics.blog.free.fr/public/perso2/testp-36.gif, on définit alors une fonction puissance pour toute valeur  http://freakonometrics.blog.free.fr/public/perso2/testh14.gif.  On a alors des graphiques de fonction puissance, comme ceux évoqués ici. Bref, la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est une valeur, alors que la puissance sera une fonction, définie sur http://freakonometrics.blog.free.fr/public/perso2/testp-40.giflorsque H_1 est une hypothèse de la forme http://freakonometrics.blog.free.fr/public/perso2/testp-41.gif.

  • L’interprétation de p

La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est la probabilité d’obtenir (au moins) la statistique obtenue pour l’échantillon si http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie. Classiquement, on rejette http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif si la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est inférieure au seuil http://freakonometrics.blog.free.fr/public/perso2/testh04.gif .
Certains interprètent la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value comme la probabilité que http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie, http://freakonometrics.blog.free.fr/public/perso2/testp-55.gif. Ce qui n’a pas de sens si l’on s’en tient au formalisme énoncé en cours: http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie, ou pas.
Ce n’est pas non plus la probabilité de rejeter, à tort, http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif , http://freakonometrics.blog.free.fr/public/perso2/testp-57.gif, correspondant à la puissance du test (cf la discussion dans le paragraphe précédant).

  • Pourquoi 5% ?

En pratique, on prend toujours comme valeur critique 5%, avec la règle décision simple suivante:

  • http://freakonometrics.blog.free.fr/public/perso2/testp-65.gif: si http://freakonometrics.blog.free.fr/public/perso2/testp-61.gif on accepte http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • http://freakonometrics.blog.free.fr/public/perso2/testp-66.gif: si http://freakonometrics.blog.free.fr/public/perso2/testp-60.gif on rejette http://freakonometrics.blog.free.fr/public/perso2/testp-33.gif

C’est Fisher qui semble être à l’origine de ce 5%, en 1925 “The value for which P=0.05, or 1 in 20, is 1.96 or nearly 2; it is convenient to take this point as a limit in judging whether a deviation ought to be considered significant or not. Deviations exceeding twice the standard deviation are thus formally regarded as significant. Using this criterion we should be led to follow up a false indication only once in 22 trials, even if the statistics were the only guide available.” On retrouve ce “un sur vingt” un an plus tard, “Either there is something in the treatment, or a coincidence has occurred such as does not occur more than once in twenty trials”.
Voilà ce que l’on retient souvent de la lecture des conseils de Fisher.
Maintenant, il lui arrive aussi, à l’occasion, de penser que 4% ce n’est pas significatif, e.g. “P is between .02 and .05. The result must be judged significant, though barely so; in view of the data we cannot ignore the possibility that on this field, and in conjunction with the other manures used, nitrate of soda has conserved the fertility better than sulphate of ammonia; the data do not, however, demonstrate this point beyond the possibility of doubt.” Ou parfois pense que 8% est significatif “P=.089. Thus a larger value of 2 would be obtained by chance only 8.9 times in a hundred, from a series of values in random order. There is thus some reason to suspect that the distribution of rainfall in successive years is not wholly fortuitous, but that some slowly changing cause is liable to affect in the same direction the rainfall of a number of consecutive years.”