Tests statistiques et intervalles de confiance

Comme je l’avais énoncé rapidement à la fin du chapitre sur les intervalles de confiance, il existe une forme de dualité entre tests et intervalles de confiance. Par exemple, si on considère un test sur la moyenne, de la forme http://freakonometrics.blog.free.fr/public/perso2/testic-01.gifcontre http://freakonometrics.blog.free.fr/public/perso2/testic-02.gif dans un modèle Gaussien, http://freakonometrics.blog.free.fr/public/perso2/test-ic03.gif, la région d’acceptation (de niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif) sera de la forme

http://freakonometrics.blog.free.fr/public/perso2/testic-05.gif
c’est à dire si

http://freakonometrics.blog.free.fr/public/perso2/testic-06.gif
Aussi, pour qu’une valeur http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif soit acceptée avec un niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif (erreur de première espèce associée à un test bilatéral), une condition nécessaire et suffisante est que http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif appartienne à l’intervalle de confiance (symétrique) centré sur la moyenne empirique.
Il y a ainsi dualité entre deux concepts

  • prendre une valeur acceptée pour un test bilatéral de niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif
  • prendre une valeur appartenant à l’intervalle de confiance symétrique de niveau 1-http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif

Considérons par exemple le test asymptotique du rapport de vraisemblance (évoqué ici), dont la région d’acceptation sera de la forme

http://freakonometrics.blog.free.fr/public/perso2/tesp100.gif
Étant donné un échantillon http://freakonometrics.blog.free.fr/public/perso2/testh03.gif, l’intervalle de confiance pour  est

http://freakonometrics.blog.free.fr/public/perso2/tesp101.gif
Numériquement, reprenons l’exemple du jeu de pile ou face abordé ici.

> n=20; alpha=.05
> X=sample(0:1,size=n,replace=TRUE)
> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=0.5,method="BFGS")$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
>
> TLR=function(p0){2*(logL(pml)-logL(p0))<qchisq(1-alpha,df=1)}
> (IC=sapply(p,TLR))
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
[37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[61]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[73]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[97] FALSE FALSE FALSE FALSE FALSE
> c(p[min(which(IC==TRUE))],p[max(which(IC==TRUE))])
[1] 0.34 0.75


Graphiquement, on visualise l’intervalle de confiance (à 95%) entre les deux traits verticaux.
De manière duale, à partir d’un intervalle de confiance http://freakonometrics.blog.free.fr/public/perso2/testp-103.gif pour , on peut construire un test. Si http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif appartient à http://freakonometrics.blog.free.fr/public/perso2/testp-103.gif, on accepte H_0,  sinon on rejette cette hypothèse. Par construction,

http://freakonometrics.blog.free.fr/public/perso2/testp-102.gif
Par exemple l’intervalle de confiance usuel – à 95% – pour une proportion est de la forme suivante

http://freakonometrics.blog.free.fr/public/perso2/testp107.gif
Si http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif appartient à cet intervalle, on accepte http://freakonometrics.blog.free.fr/public/perso2/testp109.gif (contre l’hypothèse bilatérale http://freakonometrics.blog.free.fr/public/perso2/testp101.gif).

> a=mean(X)-1.96/sqrt(n)*sqrt(mean(X)*(1-mean(X)))
> b=mean(X)+1.96/sqrt(n)*sqrt(mean(X)*(1-mean(X)))
> c(a,b)
[1] 0.3319638 0.7680362


La région entre les deux traits verticaux est la région d’acceptation.

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.”

Wald, score et rapport de vraisemblance

Vendredi, nous devrions aborder un peu la problématique des tests asymptotiques. Avant, rappelons que l’estimateur du maximum de vraisemblance de  vérifie la propriété asymptotique suivante.

Or on a une propriété intéressante sur les convergences de vecteurs Gaussiens (pour rester assez général). Soit  une suite de vecteurs aléatoires telle que  où . Alors

où  (pour un paramètre univarié, ça sera alors ). Mais ce résultat n’est intéressant que si  est connue. En fait, si on a la suite des variances vérifie , alors

Traduit sur notre estimateur du maximum de vraisemblance, cela signifie que

ou encore

Trois tests peut alors être mis en œuvre (et tant qu’à faire, on peut essayer de l’illustrer sur un exemple: encore et toujours du pile/face).

> set.seed(1)
> n=20
> X=sample(0:1,size=n,replace=TRUE)
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
> LL=sapply(p,logL)
> plot(p,LL,type="l",col="red",lwd=2)
> p0=.5
> points(p0,logL(p0),pch=3,cex=1.5,lwd=2)
> abline(v=p0,lty=2)

On a alors 20 tirages de pile ou face, on a obtenu 11 piles, et on se demande si la pièce est “juste“.Mais avant toute chose, commençons par calculer l’estimateur du maximum de vraisemblance, ainsi que le score et l’information de Fisher. Dans ce modèle de variables de Bernoulli, on connaît des formes explicites de ces quantités. Mais ici, on va utiliser la méthode la plus générale qui soit, i.e. on va maximiser la log-vraisemblance, puis on va calculer le score et l’information de Fisher à la valeur de  que l’on souhaite tester.

> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=p0,method="BFGS")$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> Df = D(f, "p")
> Df2 = D(Df, "p")
> p=p0
> score=eval(Df)
> (IF=-eval(Df2))
[1] 80
> 1/(p0*(1-p0)/n)
[1] 80

On note que l’information de Fisher calculée par double différenciation de la log-vraisemblance est identique à la formule analytique. Pour rappels, on a

Trois tests (équivalent asymptotiquement comme on le verra en cours) sont alors possibles.

Tout d’abord le test de Wald propose de travailler sur la différence entre l’estimateur du maximum de vraisemblance, et la valeur que l’on cherche à tester (comme le montre le dessin ci-dessous). Cette différence doit être “petite” si  est la vraie valeur.On peut alors utiliser la statistique suivante

Asymptotiquement, en utilisant le résultat évoqué au début, cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> alpha=0.05
> (T=(pml-p0)^2*IF)
[1] 0.1999970
> T<qchisq(1-alpha,df=1)
[1] TRUE

Ici, on accepte l’hypothèse que la pièce n’est pas pipée.

Ensuite le test du rapport de vraisemblance propose de travailler sur les valeurs de la log-vraisemblance. Là encore, cette différence doit être “petite” si  est la vraie valeur. Graphiquement, on mesure une distance normalisée

La statistique est celle d’un test bilatéral,

On pose alors

(le 2 sera détaillé un peu en cours, c’est mis là de manière à avoir trois tests équivalents). Asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=2*(logL(pml)-logL(p0)))
[1] 0.2003347
> T<qchisq(1-alpha,df=1)
[1] TRUE

Là encore, on accepte l’hypothèse que la pièce n’est pas pipée.

Enfin, le test du score propose de travailler sur la pente de la log-vraisemblance en  . Cette pente doit être “petite” si  est la vraie valeur.

La pente correspond au score,

La statistique de test est alors ici

Et là encore, asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=slope^2/IF)
[1] 0.2
> T<qchisq(1-alpha,df=1)
[1] TRUE

Et le test accepte ici encore l’hypothèse que la pièce n’est pas pipée.