Interpréter un test de Wald

Jeudi, dans le cadre du cours de statistique mathématiques, nous avions passé un peu de temps à mettre en oeuvre le test de Wald, et discuter la construction et l’interprétation d’un test. On reste ici sur un test paramètrique. La vraie valeur du paramètre, on va la fixer, pour générer des données

> p=.5

En l’occurence, on va considérer un modèle de Bernoulli, i.e. une succession de tirages de pièces (pile ou face)

> n=20
> set.seed(1)
> echantillon=sample(c("pile","face"),
+    size=n,replace=TRUE,
+    prob=c(p,1-p))
> echantillon
[1] "face" "face" "pile" "pile" "face" "pile" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face" "pile" "pile" "face" "pile"

Classiquement, on va recoder en https://latex.codecogs.com/gif.latex?\{0,1\}

> X=(echantillon=="pile")*1
> X
[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1

La vraisemblance est ici donnée par

> L=function(p) prod(dbinom(X,1,p))
> logL=function(p) sum(log(dbinom(X,1,p)))
> xp=seq(0,1,by=.01)
> yL=Vectorize(L)(xp)
> plot(xp,yL,type="l")

mais on va plutôt travailler sur la log-vraisemblance

> ylogL=Vectorize(logL)(xp)
> plot(xp,ylogL,type="l")

Le maximum de cette fonction est obtenu, numériquement, avec

> optim(par=.2,fn=function(p) -logL(p))
$par
[1] 0.55

$value
[1] 13.76278

Le maximum est ici atteint en 0,55, ce qui correspond d’ailleurs à l’estiimateur de la méthode des moments

> mean(X)
[1] 0.55

Passons maintenant à la phase de test. Notre hypothèse https://latex.codecogs.com/gif.latex?H_0 est que la pièce n’est pas truquée,

> p0 = 1/2     #  Hypoth?se H0
> phat = mean(X)

Commençons par calculer l’information de Fisher

> IF=function(p) (p*(1-p))^(-1)
> IF(phat)
[1] 4.040404

On a utilisé l’expression vue en cours, car on a fait le calcul, au tableau. Mais on peut aussi faire un peu de calcul formel (et utiliser le faire que la dérivée seconde est linéaire, de telle sorte que l’espérance de la dérivée seconde de la log vraisemblance est la valeur de la dérivée seconde en l’espérance)

> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> df = D(f,"p")
> d2f = D(df,"p")
> p = phat
> -eval(d2f)
[1] 80.80808

ce qui correspond bien à l’information de Fisher pour l’échantillon

> n*IF(phat)
[1] 80.80808

On peut alors calculer la statistique de Wald. La procédure de test est basée sur le quantile d’une loi du chi-deux à un degré de liberté (qui correspond à la loi – asymtptotique – sous https://latex.codecogs.com/gif.latex?H_0 de la statistique de Wald)

> T=n*IF(phat)*(phat-p0)^2
> seuil=qchisq(p=.95,df=1)
> T<=seuil # TRUE : accepter H0
[1] TRUE

On accepte l’hypothèse https://latex.codecogs.com/gif.latex?H_0  si la statistique est inférieure au seuil. Une autre solution est de calculer la p value

> pvalue=1-pchisq(T,df=1)
> pvalue>0.05 # TRUE : accepter H0
[1] TRUE

On rejette l’hypothèse https://latex.codecogs.com/gif.latex?H_0  si la p-value est inférieure à 5%. 5% est ici le seuil de significativité, correspondant à la probabilité de commettre une erreur de première espèce, i.e. rejeter à tort l’hypothèse https://latex.codecogs.com/gif.latex?H_0 . Pour interpréter ce 5%, on peut générer des échantillons, sous https://latex.codecogs.com/gif.latex?H_0  (i.e. avec la probabilité 1/2), et compter la proportion des scénarios où la statistique de test dépasse le seuil

> simu_Wald=function(graine){
+   set.seed(graine)
+   X=sample(0:1,size=n,replace=TRUE)
+   phat=mean(X)
+   T=n*IF(phat)*(phat-1/2)^2
+   decision=T>qchisq(p=.95,df=1)
+   pvalue=1-pchisq(T,df=1)
+   return(list(
+     echantillon=X,
+     param=phat,
+     stat=T,
+     decision=decision*1,
+     pvalue=pvalue))
+ }
>
> D=Vectorize(function(s){
+   simu_Wald(s)$decision})(1:1e5)
> mean(D)
[1] 0.04094

On rejette à tort https://latex.codecogs.com/gif.latex?H_0  dans environ 5% des cas. Ce qui correspond effectivement au choix du seuil que nous avions fait.

On peut aussi visualier les valeurs des moyennes obtenues dans les 1000 premiers échantillons, et les représenter en rouge quand on a pris la décision de rejeter https://latex.codecogs.com/gif.latex?H_0

> P=Vectorize(function(s){
+   simu_Wald(s)$param})(1:1e5)
> plot(P[1:1000],1:1000,col=D[1:1000]+1)

On s’aperçoit ici que le test de Wald est un test sur la moyenne. Ici, si la proportion de ‘pile’ est comprise entre 0,3 et 0,7, on va accepter l’hypothèse https://latex.codecogs.com/gif.latex?H_0 .

On peut aussi regarder la distribution de la p-value

> PV=Vectorize(function(s){
+   simu_Wald(s)$pvalue})(1:1e5)
> hist(PV)

On est loin de la loi uniforme ici, mais ça vient du fait qu’on est sur une loi discrète (avec peu d’observations). Cela dit, la p-value excède 5% dans 95% des cas…

> mean(PV>.05)
[1] 0.95906

On peut ensuite regarder la distribution de la statistique de test

> St=Vectorize(function(s){
+   simu_Wald(s)$stat})(1:1e5)
> hist(St,prob=TRUE, breaks =seq(0,50,by=.2),
+      xlim=c(0,10))
> x=seq(0,50,by=.2)
> lines(x,dchisq(x,df=1),lwd=2,col="red")

Là encore, on est loin de la loi du chi-deux, mais c’est normal: c’est asymptotiquement que la statistique de Wald suit une loi du chi-deux. Notons juste que le quantile empirique, à 95% est ici

> quantile(St,.95)
95%
3.809524

qui n’est pas loi du quantile théorique de la loi du chi-deux

> seuil
[1] 3.841459

Finallement, on peut regarder la (fonction) puissance de ce test

> simu_puissance=function(graine,p1){
+   set.seed(graine)
+   X=sample(0:1,size=n,replace=TRUE,
+            prob=c(1-p1,p1))
+   phat=mean(X)
+   T=n*IF(phat)*(phat-1/2)^2
+   decision=T>qchisq(p=.95,df=1)
+   return(decision)
+ }
> puissance=function(p1){
+   p=Vectorize(function(s){
+     simu_puissance(s,p1)})(1:1e5)
+   mean(p)
+ }
> P1=seq(.05,.95,by=.05)
> PUISSANCE=Vectorize(puissance)(P1)
> plot(P1,PUISSANCE,type="b")

où on retrouve l’allure classique de la puissance d’un test bilatéral.

Pour conclure, notons que la statistique de Wald est (numériquement) proche de celle du test du rapport de vraisemblance

> T=2*(logL(phat)-logL(p0))
> T
[1] 0.2003347

et du test du score

> sc=function(x,p) x/p-(1-x)/(1-p)
> Sc=Vectorize(function(x) sc(x,p0))(X)
> T=n*(mean(Sc))^2/IF(phat)
> T
[1] 0.198

ce qui est rassurant…


One thought on “Interpréter un test de Wald”

Leave a Reply

Your email address will not be published. Required fields are marked *