Tests non-paramétriques et simulations

Lors du dernier cours de statistique, nous avons présenter les tests d’ajustment de lois. Nous avions illustré ces tests à partir de la taille d’individus (déjà utilisé pour présenter l’ajustement de lois, et l’estimation de densité), correspondant à https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_n\}.

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)]
> X=Davis$height
> n=length(X)

On notera https://latex.codecogs.com/gif.latex?%20(x_{i:n}) la statistique d’ordre, au sens où

https://latex.codecogs.com/gif.latex?%20x_{1:n}\leq%20x_{2:n}\leq\cdots\leq%20x_{n-1:n}\leq%20x_{n:n}

Parmi les outils graphiques, nous avons vu le PP plot (graphique probabilité-probabilité) et le QQ plot (graphique quantile). Le code pour créer un PP plot peut être le suivant

> PP=function(Y,F=pnorm){
+   n=length(Y)
+   x=F(sort(Y))
+   y=seq(1/n/2,1-1/n/2,by=1/n)
+   return(data.frame(x=x,y=y))
+ }

qui représente (à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?\left\{F_0(x_{i:n}),\frac{i}{n}\right\}

et celui pour le QQ plot serait

> QQ=function(Y,Q=qnorm){
+   n=length(Y)
+   x=Q(seq(1/n/2,1-1/n/2,by=1/n))
+   y=sort(Y)
+   return(data.frame(x=x,y=y))
+ }

qui représente (toujours à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?%20\left\{F_0^{-1}\left(\frac{i}{n}\right),x_{i:n}\right\}

où https://latex.codecogs.com/gif.latex?%20F_0 est la loi que l’on cherche à tester, au sens où on a https://latex.codecogs.com/gif.latex?%20H_0:F=F_0 avec comme hypothèse alternative https://latex.codecogs.com/gif.latex?%20H_1:F\neq%20F_0.

Rappelons que la fonction de réparition empirique est

https://latex.codecogs.com/gif.latex?%20\widehat{F}(x)=\frac{1}{n}\sum_{i=1}^n\boldsymbol{1}(x_i\leq%20x)

Nous avions vu que https://latex.codecogs.com/gif.latex?\widehat{F}(x_{i:n})=i/n.

En cours, nous avions utilisé une séquence https://latex.codecogs.com/gif.latex?\{1/n,2/n,\cdots,(n-1)/n,1\} alors qu’ici je centre ces valeurs, en prenant https://latex.codecogs.com/gif.latex?\{1/2n,3/2n,\cdots,(2n-3)/2n,(2n-1)/2n\}. La première raison est que le décallage à droite n’est pas joli, et surtout https://latex.codecogs.com/gif.latex?F^{-1}(1) n’a peut pas de sens (si la loi n’est pas bornée).

On va ici chercher à tester une loi normale https://latex.codecogs.com/gif.latex?\mathcal{N}(\mu_0,\sigma_0^2), notée https://latex.codecogs.com/gif.latex?%20F_0

> m0=170
> s0=9
> F0=function(x) pnorm(x,m0,s0)
> Q0=function(x) qnorm(x,m0,s0)

Le graphique probabilité-probabilité est ici

> plot(PP(X,F=F0),xlab="Theoretical Probabilities",
+      ylab="Empirical Probabilities")
> abline(a=0,b=1,col=rgb(1,0,0,.4))

alors que le graphique quantile-quantile est

> plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",
+      ylab="Empirical Quantiles")
> abline(a=0,b=1,col=rgb(1,0,0,.4))

Si https://latex.codecogs.com/gif.latex?%20F_0 est la vraie loi, les PP plot et QQ plot sont des nuages de points qui doivent être alignés suivant la première bissectrice. Comme on fait des statistiques, “alignés” doit se comprendre comme “presque alignés”, ou “pas trop éloignés de la première bissectrice”. Pour savoir si on n’est pas trop loin, on peut utiliser des simulations. Il suffit de simuler des échantillons de même taille que \boldsymbol{x} suivant la loi https://latex.codecogs.com/gif.latex?%20F_0, et de tracer les PP et QQ plots

> MVPP=MVQQ=matrix(NA,500,length(X))
> for(s in 1:500) MVPP[s,]=PP(Q0(runif(length(X))),F0)$x
> for(s in 1:500) MVQQ[s,]=QQ(Q0(runif(length(X))),Q0)$y

Commençons par le graphique quantile-quantile

> plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",
+      ylab="Empirical Quantiles",col="white")
> abscisse=Q0(seq(1/n/2,1-1/n/2,by=1/n))
> for(s in 1:500) lines(abscisse,MVQQ[s,],col="light blue",lwd=.5)
> points(QQ(X,Q=Q0))
> abline(a=0,b=1,col=rgb(1,0,0,.4))

On peut calculer les quantiles empiriques, vus verticalement

> Q1=apply(MVQQ,2,function(x) quantile(x,.025))
> Q2=apply(MVQQ,2,function(x) quantile(x,.975))
> lines(abscisse,Q1,col="red")
> lines(abscisse,Q2,col="red")

Ce sont les courbes rouges ci-dessus. On peut aussi représenter les points qui sortes de la région de confiance, ponctuellement

> plot(QQ(X,Q=Q0),xlab="Theoretical Quantiles",
+      ylab="Empirical Quantiles",col="white")
> polygon(c(abscisse,rev(abscisse)),c(Q1,rev(Q2)),border=NA,col="yellow")
> id=((QQ(X,Q=Q0)$y>=Q1)&(QQ(X,Q=Q0)$y<=Q2))+1
> points(QQ(X,Q=Q0),col=c('red','blue')[id])

Notons qu’ici 5 points, sur 200, sortent des régions de confiance.

Dans le cas du graphique probabilité-probabilité, on utilise

> plot(PP(X,F=F0),xlab="Theoretical Probabilities",
+      ylab="Empirical Probabilities",col="white")
> ordonnee=seq(1/n/2,1-1/n/2,by=1/n)
> for(s in 1:500) lines(MVPP[s,],ordonnee,col="light blue",lwd=.5)
> points(PP(X,F=F0))
> abline(a=0,b=1,col=rgb(1,0,0,.4))

Cette fois, commes les invariants dans les graphiques sont sur l’axe des ordonnées, la recherche de quantile se fait horizontalement (ce qui est peu commun)

> Q1=apply(MVPP,2,function(x) quantile(x,.025))
> Q2=apply(MVPP,2,function(x) quantile(x,.975))
> lines(Q1,ordonnee,col="red")
> lines(Q2,ordonnee,col="red")

On peut visualiser des histoires de points hors de la région d’acceptation sur le graphique ci-dessous (en insistant sur la lecture horizontale des quantiles)

> plot(PP(X,F=F0))
> polygon(c(Q1,rev(Q2)),c(ordonnee,rev(ordonnee)), border=NA,col="yellow")
> id=((PP(X,F=F0)$x>=Q1)&(PP(X,F=F0)$x<=Q2))+1
> points(PP(X,F=F0),col=c('red','blue')[id],pch=c(19,1)[id])

avec le fait que les points hors de la région d’acceptation se voient horizontalement ici.

Un autre outils, plus formel, est le test de Kolmogorov-Smirnov. Cet dernier est basé sur la fonction de répartition empirique. On peut tracer cette dernière à l’aide de la fonction suivante

> Femp = function(Y){
+   n=length(Y)
+   x=sort(Y)
+   y=seq(1/n/2,1-1/n/2,by=1/n)
+   return(data.frame(x=x,y=y))
+ }
> plot(Femp(X),type="s",xlim=c(150,200))

Une alternative est d’uilitser la fonction suivante

> F=Vectorize(function(x) mean((X<=x)))
> x=seq(150,200,by=.5)
> plot(x,F(x),type="s",ylab=" Probabilities")

Notons que ce graphique est très proche des PP et QQ plot

> plot(Femp(X)$x,Q0(Femp(X)$y),type="s")

En transformant par l’inverse de https://latex.codecogs.com/gif.latex?%20F_0 en ordonnées, on retrouve le QQ plot, et en transformant par https://latex.codecogs.com/gif.latex?%20F_0 en abscisse, on retrouve le PP plot,

> plot(F0(Femp(X)$x),Femp(X)$y,type="s")

Là encore, on peut faire des simulations, pour voir où se trouverait la fonction de répartition, si on simulait des échantillons suivant https://latex.codecogs.com/gif.latex?%20F_0

> MVF=matrix(NA,500,length(x))
> for(s in 1:500){
+   Xs=Q0(runif(length(X)))
+   Fs=Vectorize(function(x) mean((Xs<=x)))
+   MVF[s,]=Fs(x)
+ }

On peeut superposer les courbes

> plot(x,F(x),type="s",xlab="Values",ylab="Cumulated Probabilities",col="white")
> for(s in 1:500) lines(x,MVF[s,],col="light blue",lwd=.5,type="s")
> lines(x,F(x),type="s")
> lines(x,F0(x),type="l",col=rgb(1,0,0,.5))

La distance utilisée pour le test est ici

https://latex.codecogs.com/gif.latex?\Vert%20\widehat{F}-F_0\Vert_{\infty}=\sup_{x\in\mathbb{R}}\{\vert\widehat{F}(x)-F_0(x)\vert\}

Comme nous l’avions vu, cette distance peut être vu comme un maximum, en cherchant la valeur maximale des

https://latex.codecogs.com/gif.latex?\vert\widehat{F}(x_i^-)-F_0(x_i^-)\vert,\vert\widehat{F}(x_i^+)-F_0(x_i^+)\vert

> D=function(Y,F0){
+   F=Vectorize(function(x) mean((Y<=x)))
+   n=length(Y)
+   x=sort(Y)
+   d1=abs(F(x+1e-6)-F0(x+1e-6))
+   d2=abs(F(x-1e-6)-F0(x-1e-6))
+   return(max(c(d1,d2)))
+ }

Ici, la distance entre la fonction de répartition et la loi normale que l’on cherche à tester est

> D(X,F0)
[1] 0.06296857

On retrouve cette statistique à l’aide de la fonction

> ks.test(X,F0)

One-sample Kolmogorov-Smirnov test

data:  X
D = 0.062969, p-value = 0.406
alternative hypothesis: two-sided

Pour comprendre d’où sort la p-value, là encore, on peut utiliser nos simulations

> VD1=rep(NA,500)
> for(s in 1:500){
+   Xs=rnorm(length(X),m0,s0)
+   VD1[s]=D(Xs,function(x) pnorm(x,m0,s0))
+ }

Dans 38% des scénarios, la distance dépasse la distance observée (0.0629)

> mean(D(X,F0)<=VD1)
[1] 0.38

ce qui est comparable avec les 40% donnés ici. On peut visualiser cette p-value ci-dessous

> dst=density(VD1)
> plot(dst,col="darkblue")
> id=which(dst$x>D(X,F0))
> polygon(c(dst$x[id],rev(dst$x[id])),c(dst$y[id],   + rep(0,length(id))),border=NA,col=rgb(0,0,1,.4))
> abline(v=D(X,F0),lty=2)


One thought on “Tests non-paramétriques et simulations”

Leave a Reply

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