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
> 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 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 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 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 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
. Pour interpréter ce 5%, on peut générer des échantillons, sous
(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 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
> 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 .
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…
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (November 4, 2016). Interpréter un test de Wald. Freakonometrics. Retrieved October 3, 2024 from https://doi.org/10.58079/ov5d
One thought on “Interpréter un test de Wald”