Category Archives: Tests

Puissance et tests statistiques (simples)

Pour comprendre la notion de puissance dans un test statistique, reprenons Statistical power analysis for the behavioral sciences de Jacob Cohen. Comme il le rappelle Fisher (dans The Design of Experiments) ne souhaite pas “prouver” qu’une hypothèse (appelée H_0) est valide, mais il espère, en réalité, rejeter l’hypothèse en question (c’est ce que rappellent d’ailleurs Hubbard et Bayarri en 2003). Autrement dit, supposons que l’on cherche à faire un test médical, et à détecter si une personne est malade (ou pas). Si une personne est généralement malade dès lors qu’un taux de quelque chose est faible (noté https://latex.codecogs.com/gif.latex?\theta), on va chercher à tester https://latex.codecogs.com/gif.latex?H_0:\theta_\star%3E\theta_0. Dans un test statistique, on espère rejeter l’hypothèse nulle https://latex.codecogs.com/gif.latex?H_0.

The power of a statistical test of a null hypothesis is the probability that it will lead to the rejection of the null hypothesis.

Continue reading Puissance et tests statistiques (simples)

Tests basés sur la vraisemblance – score

Une autre grandeur intéressant est le score, qui est la dérivée de la vraisemblance. Intuitivement (c’est l’idée de la condition du premier ordre),   et  seront proches si les dérivées en ces points sont proches. En  la dérivée est nulle, donc on va se demander ici, tout simplement, si la dérivée en  est proche de 0. Ou pas.

Continue reading Tests basés sur la vraisemblance – score

Tests basés sur la vraisemblance – Rapport de Vraisemblance

Chose promise, chose due. J’avais dit qu’on parlerait du test de rapport de vraisemblance. L’idée – visuelle – est d’avoir une lecture dans l’autre sens : au lieu de se demander si  et  sont proches, on va se demander si   et  sont proches.

Si la fonction de vraisemblance est suffisamment régulière, on se pose la même question.

Continue reading Tests basés sur la vraisemblance – Rapport de Vraisemblance

Les tests et la logique, modus tollens

Quand on apprend la logique, on apprend la notion de modus tollens, correspondant à une notion de contraposition. Si on a une proposition du genre , alors la proposition contraposée est .  Et on apprend que les deux propositions sont équivalentes (je fais de la logique classique). Par exemple si  correspond à “feu” et  à “fumée“,  signifie que tout feu fait de la fumée. Si cette affirmation est vrai, alors il n’y a pas de fumée sans feu, et donc , autrement dit, s’il n’y a pas de fumée, c’est qu’il n’y a pas de feu.

Continue reading Les tests et la logique, modus tollens

Loi multinomiale et loi du chi-deux

La semaine passée, en cours, j’avais rappelé que quand https://latex.codecogs.com/gif.latex?\boldsymbol{N}=(N_{1},\cdots,N_{k}) décrivait le compte de https://latex.codecogs.com/gif.latex?n variable multinomiales prenant https://latex.codecogs.com/gif.latex?k modalités, la variable

https://latex.codecogs.com/gif.latex?\tilde{X}_i%20=%20\frac{N_i-np_i}{\sqrt{np_i\color{red}{(1-p_i)}}}

suit asymptotiquement une loi https://latex.codecogs.com/gif.latex?\mathcal{N}(0,1). Et plus généralement, on peut montrer que https://latex.codecogs.com/gif.latex?\tilde{\boldsymbol{X}}\sim\mathcal{N}(\boldsymbol{0},\star). Le soucis est que la matrice de variance covariance n’est pas la matrice identité. Pire que ça, elle n’est pas diagonale. Encore pire, elle n’est pas inversible.

On ne peut alors pas utiliser le joli résultat qui nous dit que si https://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}), alors https://latex.codecogs.com/gif.latex?Y=(\boldsymbol{X}-\boldsymbol{\mu})^{\text{\sffamily%20T}}\boldsymbol{\Sigma}^{-1}(\boldsymbol{X}-\boldsymbol{\mu}) suit une loi du chi-deux à https://latex.codecogs.com/gif.latex?k degrés de liberté (dimension du vecteur https://latex.codecogs.com/gif.latex?\boldsymbol{X}). Car il faut que https://latex.codecogs.com/gif.latex?\boldsymbol{\Sigma} soit inversible. Or ici, je pense qu’on peut montrer que

https://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=\prod_{i=1}^k%20p_i%20\cdot\left(\sum_{i=1}^k%20p_i%20-1\right)

(à un facteur https://latex.codecogs.com/gif.latex?n près). Et comme on a une loi de probabilité, le terme de droite est nul. Bref, https://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=0 et son inverse n’existe alors pas. Cela dit, https://latex.codecogs.com/gif.latex?0 n’est valeur propre qu’une fois, d’où la dimension https://latex.codecogs.com/gif.latex?k-1 qui apparaît sans cesse (à condition que les https://latex.codecogs.com/gif.latex?p_i soient non-nuls, ce qui était une hypothèse dans la définition de la loi multinomiale : toutes les modalités existent vraiment).

L’idée est alors de poser

https://latex.codecogs.com/gif.latex?X_i%20=%20\frac{N_i-np_i}{\sqrt{np_i}}

suit asymptotiquement une loi https://latex.codecogs.com/gif.latex?\mathcal{N}(0,\star). Et plus généralement, là encore, https://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{0},\star).

Ce qui est amusant, c’est que dans le cas binomiale (j’avais dit qu’il était intéressant de voir la loi multinomiale comme une généralisation de la loi binomiale),

https://latex.codecogs.com/gif.latex?\tilde{X}_1^2=X_1^2+X_2^2

puisque

https://latex.codecogs.com/gif.latex?\frac{(N_1-np_1)^2}{np_1(1-p_1)}=\frac{(N_1-np_1)^2}{np_1}+\frac{(N_2-np_2)^2}{np_2}

Et en fait, un peut de calculs rapides permettent de noter que

https://latex.codecogs.com/gif.latex?(N_1-np_1)^2=(N_2-np_2)^2

de telle sorte que https://latex.codecogs.com/gif.latex?\tilde{X}_1^2=\tilde{X}_2^2 alors que https://latex.codecogs.com/gif.latex?{X}_1^2\neq%20{X}_2^2. Bref, c’est clairement sur cette dernière quantité qu’il faut travailler !

Et là on fait des calculs.

L’astuce est d’utiliser la formule de Stirling

https://latex.codecogs.com/gif.latex?n!\sim%20n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}

alors pour https://latex.codecogs.com/gif.latex?n_1+\cdots+n_k=n

https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=\frac{n!}{n_1%20!%20\cdots%20n_k%20!}p_1^{n_1}\cdots%20p_k^{n_k}soit

https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim\frac{n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}}{\prod_{j=1}^k%20[n_j^{n_j+\frac{1}{2}}e^{-n_j}\sqrt{2\pi}]}p_1^{n_1}\cdots%20p_k^{n_k}

et donc

https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim%20\kappa%20\prod_{j=1}^k%20\left(\frac{np_j}{n_j}\right)^{n_j+\frac{1}{2}}

de telle sorte que

https://latex.codecogs.com/gif.latex?\log%20\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=k+\sum_{j=1}^k%20\left(n_j+\frac{1}{2}\right)%20\log%20%20\left(\frac{np_j}{n_j}\right)On peut alors faire des développements limités, en faisant intervenir

https://latex.codecogs.com/gif.latex?x_i=\frac{n_i-np_i}{\sqrt{np_i}}

et non pas la version avec https://latex.codecogs.com/gif.latex?\sqrt{np_i(1-p_i)} au dénominateur. On peut alors réécrire la somme en faisant intervenir https://latex.codecogs.com/gif.latex?k-1 termes (et non pas https://latex.codecogs.com/gif.latex?k), exactement comme dans les calculs faits dans le premier cours, avec les calculs de variances. Bref, un peu de calcul supplémentaire permet de montrer que

https://latex.codecogs.com/gif.latex?\sum_{i=1}^k%20X_i^2=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i}\overset{\mathcal{L}}{\rightarrow}%20\chi^2(k-1)

C’est donc

https://latex.codecogs.com/gif.latex?Q=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i}

et non pas

https://latex.codecogs.com/gif.latex?\tilde{Q}=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i\color{red}{(1-p_i)}}

qui sera utilisée comme statistique de test par la suite.

Pour utiliser le test du chi-deux comme test d’ajustement de loi, je pourrais reprendre l’exemple des buts lors des coupes du monde de Football

soccer1982=read.table("http://freakonometrics.free.fr/soccer1982")
S82=(soccer1982$V1+soccer1982$V2)
soccer1998=read.table("http://freakonometrics.free.fr/soccer1998")
S98=(soccer1998$V1+soccer1998$V2)
soccer2010=read.table("http://freakonometrics.free.fr/soccer2010")
S10=(soccer2010$V1+soccer2010$V3)

Identification of ARMA processes

Last week (in the MAT8181 course) in order to identify the orders of an ARMA process, we’ve seen the eacf method, and I mentioned the scan method, introduced in Tsay and Tiao (1985). The code below – to produce the output of the scan procedure – has been adapted from an old code by Steve Chen (where I included a visualization of the p-values, with the following colors)

The procedure was described in the course, last Thursday,

arma.scan=function(z,ar.max=15,ma.max=15,alpha=0.01)
{
  ym=function(z,t,m){return(z[t:(t-m)])}
  n=length(z)
  z=z - mean(z)
  cmax=ma.max + 1
  rmax=ar.max + 1
  corref=matrix(0,nrow=rmax,ncol=cmax)
  cmj.table=matrix(0,nrow=rmax,ncol=cmax)
  pv=matrix(0,nrow=rmax,ncol=cmax)
  mark=matrix(rep("X",(rmax)*(cmax)),nrow=rmax,ncol=cmax)
  Rnames=paste("AR",0:(ar.max),sep="-")
  Cnames=paste("MA",0:(ma.max),sep="-")
  rownames(corref)=Rnames
  colnames(corref)=Cnames
  rownames(cmj.table)=Rnames
  colnames(cmj.table)=Cnames
  rownames(pv)=Rnames
  colnames(pv)=Cnames
  rownames(mark)=Rnames
  colnames(mark)=Cnames
  for (m in 0:ar.max)
  {
   m1=m+1
   for (j in 0:ma.max)
   {
   j1=j+1 
   if (m == 0 && j != 0)  
   {
      racf=acf(z,plot=FALSE)$acf[1:(j+1)]    
      lamb=racf[j+1]^2    
      corref[m1,j]=round(lamb,4)
      dmj=1 + 2*sum(racf[1:j]^2)
      cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)
      pvalue =pchisq(cmj,1,lower.tail=FALSE)
      pv[m1,j]=round(pvalue,4)
      cmj.table[m1,j]=round(cmj,4)
      mark[m1,j]=ifelse(pvalue > alpha,"O","X")    
    } 
    else if (m != 0 && j == 0) 
    {
      racf=pacf(z,plot=FALSE)$acf[1:(m+1)]
      lamb=racf[m+1]^2
      corref[m1,j1]=round(lamb,4)
      dmj = 1
      cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)    
      pvalue =pchisq(cmj,1,lower.tail=FALSE)
      pv[m1,j1]=round(pvalue,4)
      cmj.table[m1,j1]=round(cmj,4)    
      mark[m1,j1]=ifelse(pvalue > alpha,"O","X")
    } 
    else
    {        
      mat1=matrix(0,nrow=m1,ncol=m1)
      mat2=matrix(0,nrow=m1,ncol=m1) 
      mat3=matrix(0,nrow=m1,ncol=m1)
      mat4=matrix(0,nrow=m1,ncol=m1)     
      for (t in (j+m+2):n)
      {
         tj1=t-j-1
         ym1=ym(z,tj1,m)
         ym2=ym(z,t,m)    

         mat1=mat1 + as.matrix(ym1)%*%ym1    
         mat2=mat2 + as.matrix(ym1)%*%ym2    
         mat3=mat3 + as.matrix(ym2)%*%ym2    
         mat4=mat4 + as.matrix(ym2)%*%ym1                
      }  
      b1=solve(mat1)%*%mat2
      b2=solve(mat3)%*%mat4
      A=b2%*%b1
      eig <-eigen(A)
      eig.val <-eig$values
      eig.val=Re(eig.val)
      eig.len=length(eig.val)
      eig.vector=eig$vectors
      lamb=min(eig.val)
      eig.vector0=eig.vector[,which.min(eig.val)]
      eig.vector0 = eig.vector0/eig.vector0[1]
      resid=(1:n)*0 
      for (t in (j+m+1):n)
      {
        z0=z[seq(t,t-m,-1)]      
        resid[t]=sum(z0 * eig.vector0)
      } 
      jm1=j + m + 1
      rx=Re(resid[jm1:n])
      racf=acf(rx,plot=FALSE)$acf[1:j]
      dmj=1 + 2*sum(racf^2)
      cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)     
      pvalue =pchisq(cmj,df=1,lower.tail=FALSE)
      corref[m1,j1]=round(lamb,4)     
      pv[m1,j1]=round(pvalue,4)
      cmj.table[m1,j1]=round(cmj,4)    
      mark[m1,j1]=ifelse(pvalue > alpha,"O","X")
    }
   } 
  } 

  cat("\n\nSCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)\n\n")
  cat("Estimates of Squared Canonical Correlation \n\n")
  print(corref)
  cat("\n\nC(m,j)\n\n")
  print(cmj.table)
  cat("\n\nChi-Square(1) Test p-value\n\n")
  print(pv)
  cat("\nSCAN Matrix \n\n")
  print(mark)

plot(0:1,0:1,col="white",xlim=c(0,nrow(pv)-1),ylim=c(0,ncol(pv)-1),axes=FALSE,xlab="AR",ylab="MA")
axis(1); axis(2)
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")[c(1,2,3,5)]
cpv=matrix(as.numeric(cut(as.vector(pv),c(-1,.01,.05,.1,2))),nrow(pv),ncol(pv))
for(i in 1:nrow(pv)){
for(j in 1:ncol(pv)){
 polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
 col=CL[cpv[i,j]])
}}
}

Consider the following simulated time series,

> s=arima.sim(n=200,model=list(ar=c(0,0,0,.4,0,0,0,.5),ma=c(0,0,1))) 
> plot(s,type="l")

The output is here

> arma.scan(s,6,6)

SCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)

Estimates of Squared Canonical Correlation 

       MA-0   MA-1   MA-2   MA-3   MA-4   MA-5   MA-6
AR-0 0.0614 0.0104 0.1862 0.3516 0.0971 0.0128 0.0000
AR-1 0.0302 0.0294 0.1501 0.0943 0.0855 0.0127 0.0385
AR-2 0.3070 0.2781 0.2140 0.0006 0.1589 0.1884 0.2243
AR-3 0.1627 0.0037 0.1927 0.2311 0.1379 0.0207 0.0376
AR-4 0.2087 0.3947 0.3653 0.3075 0.1502 0.1364 0.1013
AR-5 0.1677 0.1219 0.0110 0.0263 0.0332 0.0350 0.0044
AR-6 0.0114 0.0485 0.0561 0.0427 0.0009 0.0089 0.0308

C(m,j)

        MA-0    MA-1    MA-2    MA-3   MA-4   MA-5    MA-6
AR-0  4.1161  0.6585 12.0315 20.6512 4.5388 0.5620  0.0000
AR-1  6.1127  1.9499  9.9356  4.9145 4.7219 0.4642  1.9015
AR-2 72.6011 19.1679 14.3512  0.0337 7.9668 9.6479 11.4573
AR-3 34.9724  0.2386 10.1620 13.4082 6.7875 0.8725  1.4071
AR-4 45.8691 27.5070 19.1422 20.2835 7.3339 5.5374  3.5874
AR-5 35.7981  8.0498  0.6280  1.3543 1.8470 1.7930  0.2338
AR-6  2.2147  3.1466  3.5990  1.9904 0.0511 0.4816  1.6440

Chi-Square(1) Test p-value

       MA-0   MA-1   MA-2   MA-3   MA-4   MA-5   MA-6
AR-0 0.0425 0.4171 0.0005 0.0000 0.0331 0.4534 0.0000
AR-1 0.0134 0.1626 0.0016 0.0266 0.0298 0.4957 0.1679
AR-2 0.0000 0.0000 0.0002 0.8543 0.0048 0.0019 0.0007
AR-3 0.0000 0.6252 0.0014 0.0003 0.0092 0.3503 0.2355
AR-4 0.0000 0.0000 0.0000 0.0000 0.0068 0.0186 0.0582
AR-5 0.0000 0.0046 0.4281 0.2445 0.1741 0.1806 0.6287
AR-6 0.1367 0.0761 0.0578 0.1583 0.8212 0.4877 0.1998

SCAN Matrix 

     MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6
AR-0 "O"  "O"  "X"  "X"  "O"  "O"  "X" 
AR-1 "O"  "O"  "X"  "O"  "O"  "O"  "O" 
AR-2 "X"  "X"  "X"  "O"  "X"  "X"  "X" 
AR-3 "X"  "O"  "X"  "X"  "X"  "O"  "O" 
AR-4 "X"  "X"  "X"  "X"  "X"  "O"  "O" 
AR-5 "X"  "X"  "O"  "O"  "O"  "O"  "O" 
AR-6 "O"  "O"  "O"  "O"  "O"  "O"  "O"

with the following graph

Of course, it is possible to ask for larger values,

> arma.scan(s,12,12)

The graph is now

Maximum Likelihood versus Goodness of Fit

Thursday, I got an interesting question from a colleague of mine (JP). I mean, the way I understood the question turned out to be a nice puzzle (but I have to confess I might have misunderstood). The question is the following : consider a i.i.d. sample https://latex.codecogs.com/gif.latex?\{X_1,\cdots,X_n\} of continuous variables. We would like to choose between two (parametric) families for the distribution, https://latex.codecogs.com/gif.latex?\mathcal{F}=\{F_{\boldsymbol%20\alpha};\boldsymbol%20\alpha\in\mathcal{A}\} and https://latex.codecogs.com/gif.latex?\mathcal{G}=\{G_{\boldsymbol%20\beta};\boldsymbol%20\beta\in\mathcal{B}\}. If we use maximum likelihood techniques, we get two estimators, one for each family, https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol%20\alpha} and https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol%20\beta}. Clearly, https://latex.codecogs.com/gif.latex?F_{\widehat{\boldsymbol%20\alpha}}(\cdot) is a much better than https://latex.codecogs.com/gif.latex?G_{\widehat{\boldsymbol%20\beta}}(\cdot), in the sense of a standard goodness of fit test (e.g. Kolmogorov-Smirnov since the sample is assumed to be obtained from a continuous variable). Does that mean that family is https://latex.codecogs.com/gif.latex?\mathcal{F} (somehow) better than family https://latex.codecogs.com/gif.latex?\mathcal{G} ?

This is my interpretation of the question, and I found it amusing. So I will try to show (using simulated samples) that some odd situations can easily be obtained.

Consider a sample from a mixture of log-normal distributions,

>  set.seed(228)
>  X=exp(c(rnorm(50,1,1),rnorm(50,2,1.2)))

Consider two standard families for positive random variables: a Gamma distribution and a lognormal distribution.

> library(MASS)
> ab=fitdistr(X,"gamma")
> ms=fitdistr(X,"lognormal")

If we want to visualized those two distributions, let us use

> vab=pgamma(u,ab$estimate[1],ab$estimate[2])
> vms=plnorm(u,ms$estimate[1],ms$estimate[2])
> plot(ecdf(X))
> lines(u,vab,col="red")
> lines(u,vms,col="blue")

Here, we get

What else can we say ? actually, we can also compute Kolmogorov-Smirnov statistic,

https://latex.codecogs.com/gif.latex?D_n=\sup_x%20|\widehat%20F_n(x)-F_\star(x)|where

https://latex.codecogs.com/gif.latex?\widehat%20F_n(x)={1%20\over%20n}\sum_{i=1}^n%20\boldsymbol{1}_{X_i\leq%20x}

This can be done using

> ks.test(X,"plnorm",ms$estimate[1],ms$estimate[2])

One-sample Kolmogorov-Smirnov test

data:  X
D = 0.0693, p-value = 0.7231
alternative hypothesis: two-sided

> ks.test(X,"pgamma",ab$estimate[1],ab$estimate[2])

One-sample Kolmogorov-Smirnov test

data:  X
D = 0.148, p-value = 0.02507
alternative hypothesis: two-sided

From a theoretical point of view, we should not look at the p-values, since the null-distribution is based on a fixed distribution, not a fitted one (see the Lilliefors tests for normal samples). But still. The Gamma distribution seems to be very far away from the true distribution. The statistics is twice the one we have with our lognormal distribution. And one p-value is 72%, while the other one is 2.5%. Here, we should prefer this lognormal distribution to that Gamma one. But here, we did consider only one distribution in each family. Does that mean that we cannot find one Gamma distribution that will be better than all possible lognormal distributions ? Better, for instance, according to Kolmogorov-Smirnov statistics…

Well, it is possible to use another strategy to find appropriate parameters. We can minimize this statistic actually. Define

> ks1=function(ms) {m=ms[1];s=ms[2];ks.test(X,"plnorm",m,s)$statistic}
> ks2=function(ab) {a=ab[1];b=ab[2];ks.test(X,"pgamma",a,b)$statistic}

and compute

> n1=nlm(ks1,c(ms$estimate[1],ms$estimate[2]))
> n1
$minimum
[1] 0.05252692

$estimate
[1] 1.547437 1.121864
> n2=nlm(ks2,c(ab$estimate[1],ab$estimate[2]))
> n2
$minimum
[1] 0.04737725

$estimate
[1] 1.1449041 0.167041

So here, it is possible to find a distribution much closer to the empirical sample, within the Gamma family actually.

>  vab=pgamma(u,n2$estimate[1],n2$estimate[2])
>  vms=plnorm(u,n1$estimate[1],n1$estimate[2])
>  lines(u,vab,col="red",lwd=2)
>  lines(u,vms,col="blue",lwd=2)

What would be the point here? Maybe just the idea that the maximum likelihood estimator is only one estimator among a lot of them. And if it has interesting asymptotic properties, on small samples, it might not be the best estimator to consider…

And to be completely honest, I’ve been cheating here… I mean, not really cheating (not more than any researcher using a statistical test to validate the findings). But here, I did fix the seed of the random number generator. Actually, such example does not occur that frequently. Here, out of 1000 samples, I got this odd conclusion almost 15 times. And the smaller the sample, the more likely we can observe that story, where the maximum likelihood estimator can be far away from the best fit. Here is the proportion of opposite conclusions, as a function of the sample size,

> SIM=function(ns=1000,n=100){
+ t=0
+ for(s in 1:ns){
+  set.seed(s)
+  X=exp(c(rnorm(n/2,1,1),rnorm(n/2,2,1.2)))
+  ks1=function(ms) {m=ms[1];s=ms[2];ks.test(X,"plnorm",m,s)$statistic}
+  ks2=function(ab) {a=ab[1];b=ab[2];ks.test(X,"pgamma",a,b)$statistic}
+  library(MASS)
+  ab=fitdistr(X,"gamma")
+  ms=fitdistr(X,"lognormal")
+  n1=nlm(ks1,c(ms$estimate[1],ms$estimate[2]))
+  n2=nlm(ks2,c(ab$estimate[1],ab$estimate[2]))
+  if( (ks.test(X,"plnorm",ms$estimate[1],ms$estimate[2])$statistic-
+  ks.test(X,"pgamma",ab$estimate[1],ab$estimate[2])$statistic)
+ *(n1$minimum-n2$minimum)<=0 ) t=t+1
+ }
+ return(t/ns)}

> VM=rep(NA,20)
> VS=seq(10,200,by=10)
> for(i in 1:20){VM[i]=SIM(n=VS[i],ns=1000)}
> plot(VS,VM,type="p")

So to provide a more complete answer to JP’s question, with a very large sample, I guess that your intuition should be valid…. but clearly not on a small sample.

De la significativité (statistique), suite

Suite à mon rapide billet sur l’article donnant les conclusions d’une étude des effets de la e-cigarette (et la comparaison avec le patch pour arrêter de fumer), j’avais entendu cette phrase “ça veut dire quoi ‘statistiquement non significatif’ ?“. Comme c’est une très bonne question, qui plus est relativement complexe, j’ai pris un peu de temps pour y répondre… Il ne s’agit pas de complexité mathématique, en s’entend car tous les cours de statistique abordent cette notion, mais plus d’une notion complexe à vulgariser. On va revenir un peu sur l’histoire des tests de significativité et sur les travaux de Ronald Fisher (pour comprendre l’origine des 5% utilisé comme seuil lors de l’utilisation de la p-value), mais il est possible d’aller voir directement les exemples.

  • Un peu d’histoire

Si on veut commencer par la préhistoire, les travaux les plus anciens qui posent la question de la significativité sont ceux de John Arbuthnot, en 1710. En prenant 82 années de données de natalité, à Londres, il avait observé que chacune des années, il y avait eu plus de naissances de garçons que de filles (évoqué dans Hacking (1975)).

Il s’était demandé, avec la terminologie d’il y a plus de 300 ans, si cette différence pouvait être jugée comme ‘statistiquement significative. Je pense que l’on peut utiliser le terme ‘statistiquement‘ car John Arbuthnot pose la question en termes probabilistes, très clairement,

Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un ‘test de significativité‘, là encore avec notre terminologie actuelle. Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste (que l’on verra formalisé bien plus tard), en se demandant s’il est ‘extrêmement probable‘ qu’il y ait une différence entre les deux mesures. Il avait alors fait un test de comparaison des valeurs moyennes, et noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important. A la fin du XIXème siècle, Francis Edgeworth reprendra cette idée (sur la comparaison de moyennes) lorsqu’il écrira en 1885,

où on retrouve les premières idées qui serviront de base aux travaux ultérieurs. Car pour bien comprendre l’origine de la notion de significativité en statistique, il faut relire les travaux de Ronald Fisher des années 20, mais aussi de Student (William Gosset) et Karl Pearson. Ronald Fisher s’interroge, dans son papier de 1921, sur la valeur d’un coefficient de corrélation (c’est dans ce papier qu’il introduit formelle sa transformation avec une tangente hyperbolique z = {1 \over 2}\ln\left({1+r \over 1-r}\right) = \text{arctanh}(r) pour faire un test de nullité d’un coefficient de corrélation), où il note

Aussi, dans le papier, une valeur (ici une corrélation) peut être jugée comme ‘significative‘ (c’est à dire ‘significativement non nulle‘). Encore une fois, l’idée n’est pas nouvelle, Edgeworth aurait dit ‘non accidentelle‘ mais Fisher pose ici la terminologie, en parlant de ‘significativité‘. L’année suivante, en 1922, il parle de l’idée de tester la significativité, dans le cadre d’un modele de régression (on est toujours dans le contexte de la corrélation) ou en 1924, on retrouve l’idée de tester la significativité, Mais c’est surtout dans son livre Statistical Methods for Research Workers que l’on va trouver une formalisation claire (et très illustrée) de cette idée, et surtout la connecter à la notion de p-value (on ne peut pas parler de significatité sans parler de p-value). Pour l’anecdote, dans l’ouvrage de Ronald Fisher, on y parle de p-value (oui, avec une majuscule). Cela dit, la paternité de la p-value est à attribuer davantage à Pearson, par exemple dans son papier datant de 1900, ou encore, sur un autre exemple avec en plus l’interprétation usuelle en terme de significativité (sans toutefois utiliser ce terme), puisque Pearson dit qu’on a 1 chance sur 8 d’avoir une valeur aussi improbable. Plus loin, il dit aussi qu’arriver avec 1 chance sur 70 peut être qualifier de très improbable (un long débat suivra sur le seuil à partir duquel, effectivement, un événement peut être qualifié d’improbable) Mais revenons au livre de Ronald Fisher, et plus spécifiquement au chapitre 4, qui propose précisément d’étudier la ‘significativité de la différence entre deux moyennes

On y retrouve le test de Pearson basé sur la statistique du chi-deux, comme test d’ajustement (goodness of fit) et d’indépendance. et la comparaison de la p-value avec (ce qui va correspondre aux mythiques) 5%. Il parle plus loin d’une ‘ligne conventionnelle‘ que l’on tracera à 5%, Il faut d’ailleurs un lien entre une p-value à 5% (correspondant à un événement qui survient avec 1 chance sur 20) et le fait de s’éloigner de 2 écart-types de la moyenne pour une loi normale (qui survient avec 1 chances sur 22) Tout au long du livre, on retrouve évoqué le chiffre mythique 5% qui restera à la postérité. Le plus intéressant pour apprendre comme juger la significativité d’une valeur, on peut regarder les exemples, comme l’exemple 11 l’exemple 12, l’exemple 27, ou encore l’exemple 32 (pour en prendre juste quelques uns) (on retrouve l’utilisation de 2 écart-types pour une loi normale). La discussion amorcée dans ce livre sera poursuivie dans son papier de 1926, qui commence avec Le papier est passionnant, car il pose les bases de ce que font tous les statisticiens depuis 80 ans (il n’est pas surprenant de voir cet article reproduit dans le second volume des Breakthroughs in Statistics), où le seuil de 5% est encore utilisé. Il ne cesse tout au long de l’article de justifier cette règle de ‘un sur vingt‘, tout en autorisant d’autres valeurs, comme 2%, ou 10% (Fisher semble moins orthodoxe que certains utilisateurs de la notion de significativité) Il faudra ensuite attendre les travaux de Jerzy Neyman et Egon Pearson. qui vont axiomatiser et formaliser les tests statistiques, mais ça dépasse (de loin) l’objectif de mon billet d’aujourd’hui. Donc on va repousser à plus tard la présentation de leurs apports. Pour illustrer tout ce qu’on vient de voir, considérons deux exemples simples.

  • Petit exemple, tester l’indépendance

Comme Ronald Fisher, commençons par ce qui me semble le plus simple car on n’a pas besoin de parler modèle probabiliste ici. On se demande si deux variables sont indépendantes, ou pas. Deux cas peuvent se présenter:

  • on a des variables continues, et on veut voir si elles sont non-corrélées (oui, pour commencer on va simplifier, et dire que si elles sont non-corrélées… ça va se traduire par une notion que l’on peut assimiler – rapidement – à une notion d’indépendance. Promis, d’ici 6 mois il y aura des dizaines de billets sur le sujet sur le blog puisque je redonne mon cours sur la dépendance cet hiver). Je laisse les plus motivés jouer avec les données qui traînent sur le blog pour voir si la taille et le poids sont corrélées, ou utiliser un exemple que j’avais l’habitude de présenter dans mes cours de modèles de régression, où on se demande si – chez les femmes – le tour de poitrine et le tour de taille sont corrélés.
  • on a des variables catégorielles, par exemple la couleur des yeux, et la couleur des cheveux. Et on se demande si les deux sont couleurs sont indépendantes, ou pas.

Regardons ce dernier exemple, avec un petit jeu de 25 observations.

> X=c(rep("Brun",8+5),rep("Blond",9+3))
> Y=c(rep("Noisette",8),rep("Bleu",5),rep("Noisette",3),rep("Bleu",9))

On peut résumer l’information dans un tableau de contingence,

> T=as.matrix(table(X,Y))
> T
       Y
X       Bleu Noisette
  Blond    9        3
  Brun     5        8

On va essayer de comparer ce tableau avec ce qu’on devait avoir si les deux variables étaient indépendantes (car c’est ce que l’on souhaite tester : on veut savoir si le tableau que l’on a est significativement différent de ce qu’on aurait avec des variables indépendantes). Il faut ici ressortir la définition de l’indépendance entre deux variables aléatoires, \mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B) = \mathbb{P}(X\in\mathcal A) \times\mathbb{P}(Y\in\mathcal B) pour tout \mathcal A et \mathcal B. Par exemple \mathcal A peut être avoir des cheveux blonds, et \mathcal B avoir des yeux bleus. Pour rappel, on avait en tout 25 personnes ayant ces deux caractéristiques en même temps, donc empiriquement, on avait observé une probabilité de 9/25 soit 36% (c’est le terme de gauche). Pour le terme de droite, on notera que \mathbb{P}(X\in\mathcal A), c’est la probabilité d’avoir les cheveux blonds. Empiriquement, on en avait 9+3 sur 25, soit 48%. Et \mathbb{P}(Y\in\mathcal B) c’est la probabilité d’avoir les yeux bleus. Or on avait 9+5 personnes avec ces caractéristiques, soit 56%. Si les deux variables étaient indépendantes, avec la formule précédente, la probabilité d’avoir à la fois les cheveux blonds et les yeux bleus serait 48% x 56%, soit 26.88%. A comparer aux 36% observés. En fait, au lieu de comparer des pourcentages, on peut aussi comparer des nombres: comme il y a 25 personnes en tout, 26.88% signifie que dans notre groupe, 6.72 personnes devraient avoir les cheveux blonds et les yeux bleus. On peut s’amuser à faire ces calculs pour tous les \mathcal A et \mathcal B.

> Tind=chisq.test(T)$expected
> Tind
       Y
X       Bleu Noisette
  Blond 6.72     5.28
  Brun  7.28     5.72

Visuellement, on peut mettre ça en forme,

(sur le graphique de droite, les traits noirs correspondent aux nombres sous hypothèse d’indépendance, et en couleur, les valeurs observées). Bon, maintenant, il faut réfléchir un peu… On se demande si nos deux tableaux sont proches ou pas. Ou, pour utiliser notre mot de la journée, significativement proches. Alors là, ça va devenir technique le temps d’un rapide paragraphe. Dans chaque case du tableau, on compte des gens. Sur nos n personnes, on peut s’attendre à en observer un nombre aléatoire, noté N_{\mathcal{A},\mathcal{B}}, qui suit une loi binomiale N_{\mathcal{A},\mathcal{B}}\sim\mathcal{B}(n,\mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B)) (c’est la définition de la loi binomiale). On notera que quand n est grand, la loi normale est alors proche de la loi Gaussienne. Dans chaque case, on a des lois normales (non indépendantes, car les nombres sont contraints, par exemple au total, il faut n personnes), et en prenant les carrés des lois normales (centrées et réduites), on va obtenir une loi du chi-deux (je renvoie vers un précédant billet pour les aspects techniques). Prendre le carré est intéressant, et naturel, car ça fait penser à une distance Euclidienne (usuelle). On va alors poser
Q=\sum_{\mathcal{A},\mathcal{B}} \frac{[n_{\mathcal{A},\mathcal{B}}-n_{\mathcal{A},\mathcal{B}}^\perp]^2}{n_{\mathcal{A},\mathcal{B}}^\perp}
n_{\mathcal{A},\mathcal{B}} est le nombre de personnes observées, de caractéristiques (jointes) \mathcal A et \mathcal B. et n^\perp_{\mathcal{A},\mathcal{B}} est ce qu’on aurait si les variables étaient indépendantes (calculé auparavant), soit

n^\perp_{\mathcal{A},\mathcal{B}}=n\times \frac{n_{\mathcal{A}}}{n}\times \frac{n_{\mathcal{B}}}{n}

avec les notations naturelles. Ici,

> sum((T-Tind)^2/Tind)

[1] 3.380994

Bon… on en fait quoi maintenant de ce 3.38 ? C’est là qu’on utilise le petit résultat évoqué auparavant, qui nous dit que, si n est grand, et que les variables sont effectivement indépendantes, alors Q va suivre une loi du chi-deux (en l’occurrence à 1 degré de liberté, pour tenir compte des diverses contraintes). On va alors se demander s’il est possible, ou vraisemblable, d’avoir 3.38 avec une loi du chi-deux.

> pchisq(3.380994,df=1)
[1] 0.9340477
> 1-pchisq(3.380994,df=1)
[1] 0.06595227

Aussi, une loi du chi-deux à un degré de liberté a 93.4% chances d’être inférieure à la valeur obtenue, ou encore

\mathbb{P}(Q>3.38) \approx 6.59\% avec Q\sim \chi^2(1)

C’est cette grandeur que l’on appelle p-value et que l’on va essayer de comparer au seuil désormais mythique de 5%,

  • si cette probabilité est inférieure à 5%, on va rejeter l’hypothèse d’indépendance
  • si cette probabilité est supérieure à 5%, on va retenir l’hypothèse d’indépendance (pour faire simple)
> chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 3.381, df = 1, p-value = 0.06595

Les plus attentifs devraient me faire remarquer que c’est vaseux mon histoire : tout ce que je raconte est valide si n est grand, et dire que 25 c’est grand… c’est limite ! Tout d’abord, je dirais que c’est ce qui est utilisé dans un paquet d’études en médecine (peut-être un peu plus, mais guère plus). Ensuite je dirais que, malheureusement, la plupart des résultats que l’on voit dans les cours de statistiques sont de ce genre, à savoir des résultats asymptotiques, valides seulement si n est grand (c’est pareil en économie, ou des tests asymptotiques sont parfois évoqués avec 25 années d’observations). Mais heureusement, en s’accrochant un peu, on peut s’affranchir de cette hypothèse. A condition de faire un peu de bootstrap. En fait, l’idée est très très simple. On peut utiliser des générateurs de nombres aléatoires pour générer des tableaux de contingences pour lesquels les variables sont indépendantes. Car générer des variables aléatoires indépendantes, c’est ce qui est le plus simple. En fait, ici, on va se contenter de mélanger les variables,

> s=1
> set.seed(s)
> Xs=sample(X)
> Ys=sample(Y)

En faisant une permutation, je m’assure d’avoir autant de personnes dont la couleur des cheveux est \mathcal A que dans la population initiale (je ne fais que permuter) et de personnes dont la couleur des yeux est \mathcal B. Par contre, comme je mélange indépendamment, je suis certain que les variables sont indépendantes,

> Ts=as.matrix(table(Xs,Ys))
> Ts
       Ys
Xs      Bleu Noisette
  Blond    7        5
  Brun     7        6

Générer un tableau de contingence, c’est juste pour illustrer, car l’idée est de générer mille (ou plus encore) tableaux de contingence obtenus en simulant des variables indépendantes, et de voir ce que cette statistique du chi-deux donne

> Q=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ Xs=sample(X)
+ Ys=sample(Y)
+ Ts=as.matrix(table(Xs,Ys))
+ Q[s]=sum((Ts-Tind)^2/Tind)}

Si on regarde les valeurs obtenues

>  table(Q) 
Q
0.0509906759906761  0.337162837162837   1.06560106560107   1.92411754911755 
               298                262                200                133 
  3.38099400599401   4.81185481185481    6.9971694971695   9.00037462537462 
                60                 33                  9                  3 
  11.9141275391275    14.489676989677 
                 1                  1

autrement dit, dans 29.8% des scénarios, on avait une distance (un statistique) valant 0.05, dans 26.2% des scénarios, on avait 0.33, et dans 20% des scénarios, on avait 1.06. Etc. Si on cumule on obtient

>  cumsum(table(Q))/10
0.0509906759906761  0.33162837162837   1.06560106560107    1.92411754911755 
3.38099400599401    4.81185481185481    6.9971694971695     9.00037462537462 
11.9141275391275     14.489676989677 

Pour rappel, on avait eu

>  sum((T-Tind)^2/Tind) 
[1] 3.380994

qui peut effectivement être obtenu avec des variables indépendantes. Mais seulement dans 6% des scénarios. En fait, dans 95.3% des scénarios, on a eu une valeur inférieure ou égale à celle obtenue sur notre échantillon. Et dans 4.7% des scénarios, on a strictement dépassé cette valeur. On pourrait parler de p-value obtenue par bootstrap, ou rééchantillonage. On voit que l’on flirte ici avec le seuil des 5%… Visusellement, on a

> plot(c(0,as.numeric(names(table(Q)))),c(0,
+ cumsum(table(Q))/10),type="s",xlab="",ylab="")
> abline(v=sum((T-Tind)^2/Tind),lty=2)
> lines(seq(0,15,by=.01),pchisq(seq(0,15,by=.01),1)*100,col="red")

Ici, on est bien embêté pour conclure quoi que ce soit… Ce qui est amusant, c’est qu’avec un échantillon plus grand (trois fois plus grand par exemple), et les mêmes pourcentages, il n’y a plus aucune ambiguïté,

> T=T*3
> T
       Y
X       Bleu Noisette
  Blond   27        9
  Brun    15       24
>  Tind=chisq.test(T)$expected
> Tind Y X Bleu Noisette Blond 20.16 15.84 Brun 21.84 17.16

que l’on peut encore visualiser ci-dessous,

Ici, on obtient une p-value de l’ordre de 0.1%,

>  chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 10.143, df = 1, p-value = 0.001449

On peut aussi faire du bootstrap si on trouve que 75 observations, c’est trop faible, mais la conclusion est la même,

  • Petit exemple, sur la comparaison de moyenne de deux groupes

Cette fois, on a deux populations, et on se demande si les deux groupes sont ‘comparables‘. On peut penser à comparer le nombre de personnes guéries au sein de deux populations, une qui a eu un placebo, et l’autre un médicament. Cela dit, cet exemple simple est en fait compliqué car ‘guérie’ n’est peut-être pas une notion claire. Prenons plutôt des quantités que l’on peut mesurer. Disons la taille des gens, et demandons nous si dans un groupe, les gens sont significativement plus grands que dans l’autre. Je vais tricher un peu… je vais prendre deux populations que je sais être différentes, a priori. Un groupes d’hommes, et un groupes de femmes. On va noter \{x_1,\cdots,x_{n_X}\} et \{y_1,\cdots,y_{n_Y}\} ces deux échantillons, car il va falloir formaliser un minimum.

> 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)] 
> Davis=Davis[order(Davis$height),]
> attach(Davis)
> set.seed(1)
> X=sample(height[sex=="F"],size=6)
> Y=sample(height[sex=="M"],size=5)

On a ici 5 personnes dans le premier groupe et 6 personnes dans le second. Et on se demande si les tailles sont significativement différentes entre les deux groupes,

> (mx=mean(X))
[1] 165.5
> (sx=sd(X))
[1] 4.679744
> (my=mean(Y))
[1] 178
> (sy=sd(Y))
[1] 7.968689

On veut savoir si l’écart entre 165,5 cm et 178 cm est significativement différent de zéro (on pourrait aussi se demander s’il est significativement positif, mais ça serait une autre histoire). On va faire une hypothèse ici : on va supposer que nos échantillons sont obtenus comme des tirages de variables aléatoires, indépendantes, x_i=X_i(\omega) et y_j=Y_j(\omega). Le modèle qu’on va supposer est que la taille des gens suit une loi normale, conditionnellement à leur sexe,

\left\{\begin{array}{l}X_i\sim\mathcal{N}(\mu_X,\sigma_X^2) \\Y_j\sim\mathcal{N}(\mu_Y,\sigma_Y^2)\end{array}\right.

Si on veut visualiser nos données, on peu utiliser

> u=seq(150,200,by=.25)
> v1=dnorm(u,mx,sx)
> v2=dnorm(u,my,sy)
> plot(u,v1,col="red",type="l",ylim=c(0,.09),axes=FALSE)
> axis(1)
> lines(u,v2,col="blue")
> points(X,rep(.075,length(X)),col="red")
> points(Y,rep(.075,length(Y)),col="blue")
> abline(v=mx,lty=2,col="red")
> abline(v=my,lty=2,col="blue")

 

Le plus simple aurait été de supposer que les deux lois ont la même variance, mais quand on regarde les densités sur le graphique ci-dessus, il faut rester crédible. Bref, on va utiliser le test de Welch. L’idée est d’utiliser un petit exercice de probabilité. Soient \{X_1,\cdots,X_{n_X}\} et \{Y_1,\cdots,Y_{n_Y}\} des variables aléatoires, Gaussiennes (comme supposé juste avant de faire un dessin) et indépendantes. Posons alors
\overline{X}=\frac{1}{n_X}\sum_{i=1}^{n_X} X_i
et
s_X=\frac{1}{n_X-1}\sum_{i=1}^{n_X} [X_i-\overline{X}]^2
(avec les quantités similaires pour le second échantillon) alors si \mu_X=\mu_Y, et si on note
s = \sqrt{{s_X^2 \over n_X} + {s_Y^2 \over n_Y}}.
la variable aléatoire
T = {\overline{X} - \overline{Y} \over s[latex]<br /> va suivre une loi de Student (connue depuis les travaux de <a href="http://freakonometrics.hypotheses.org/1968">William Gosset</a>), avec</p> [latex display="true"]\frac{(s_X^2/n_X + s_Y^2/n_Y)^2}{(s_X^2/n_X)^2/(n_X-1) + (s_Y^2/n_Y)^2/(n_Y-1)}

degrés de liberté. On notera qu'il ne s'agit pas ici d'un résultat asymptotique. Par contre, il repose intégralement sur une hypothèse - forte - de normalité des tailles des individus. L'idée est alors la même qu'auparavant : on va calculer cette statistique sur notre échantillon,

> nx=length(X)
> ny=length(Y)
> s=sqrt(sx^2/nx+sy^2/ny)
> T=(mx-my)/s
> T
[1] -3.091371

et se demander s'il est vraisemblable d'avoir une telle valeur avec une loi de Student. Le nombre de degrés de liberté est ici

> d=(sx^2/nx+sy^2/ny)^2/((sx^2/nx)^2/(nx-1)+(sy^2/ny)^2/(ny-1))
> d
[1] 6.218682

On pourrait alors calculer
\mathbb{P}(T>t)T\sim\mathcal{S}td(d)

mais ça serait un test pour savoir si la différence est significativement positive, ce qui n'est pas ce qu'on avait dit que l'on cherchait : on veut ici significativement non-nulle.

\mathbb{P}(\vert T\vert> \vert t\vert)T\sim\mathcal{S}td(d)

Si on veut visualiser un peu tout ça, on peut utiliser

> u=seq(-4,4,by=.01)
> v=dt(u,df=d)
> plot(u,v,type="l")
> abline(v=T,col="red")
> u=seq(-4,T,length=100)
> polygon(c(u,rev(u)),c(dt(u,df=d),rep(0,100)),col="red",border=NA)

Pour calculer ces probabilités, on utilise

> pt(T,d)+(1-pt(-T,d))
[1] 0.02038577

et comme la loi de Student est symétrique par rapport à l’origine,

> pt(T,d)*2           
[1] 0.02038577

Aussi, la p-value est de l'ordre de 2%, ce qui est suffisamment petit pour dire que, oui, nos moyennes sont différentes. C'est d'ailleurs ce que renvoie la commande

> t.test(X,Y)

	Welch Two Sample t-test

data:  X and Y
t = -3.0914, df = 6.219, p-value = 0.02039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.31039  -2.68961
sample estimates:
mean of x mean of y 
    165.5     178.0

On avait vu tout à l'heure que l'idée de la valeur de 5% était liée au scénario qui se répète avec 1 chance sur 20 : si on pouvait générer d'autres échantillons, indépendants, alors 1 fois sur 20 on devrait atteindre une statistique de test de l'ordre de celle observée. Pour générer des échantillons vraiment indépendants, on avait utilisé du bootstrap, en rééchantillonnant dans notre population de test. Ici, on a bien plus que nos 11 observations, on en a en réalité 200. On peut se demande ce qui se serait passé si on avait constitué d'autres échantillons,

> Diff=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ X=sample(height[sex=="F"],size=6)
+ Y=sample(height[sex=="M"],size=5) 
+ Diff[s]=(mean(X)-mean(Y))/(sqrt(sd(X)^2/nx+sd(Y)^2/ny))}
> hist(Diff,col="light blue",border="white")

Sur aucun échantillon, la moyenne du premier groupe n'a dépassé la moyenne du second. Alors qu'on aurait pu penser en avoir une vingtaine (2% des mille scénarios générés). Ceci vient peut-être du fait que l'hypothèse de normalité pour les tailles peut être trop restrictive...

  • Pour aller plus loin

Le billet aura été un peu long, mais on aura eu le temps de voir en détails d'où vient cette notion de significativité statistique, et ses liens avec la p-value. On aura même eu le temps de voir sur deux exemples l'importance des hypothèses sous-jacentes, en particulier si le test est asymptotique, ou s'il repose sur des hypothèses fortes de lois. Mais ce n'était qu'une courte mise en bouche. Pour aller plus loin, je suggère la lecture de Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle, M.V. (2008) Sound and fury: McCloskey and significance testing in economics, Liberman, M. (2013). "Significance", in 1885 and today et le passionnant Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Et sur les p-value, Goodman, S. (2002) A Dirty Dozen:Twelve P-Value Misconceptions. Je devrais aussi citer Ziliak, S. & McCloskey, D. (2008). The cult of statistical significance, mais je ne l'ai pas encore lu. Par contre, Xi'an en parlait sur son blog. Et si j'étais courageux, je tenterais une lecture bayésienne de ce billet.... promis, dès que je trouve un peu de temps !

Logistic regression and categorical covariates

A short post to get back – for my nonlife insurance course – on the interpretation of the output of a regression when there is a categorical covariate. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)
> tail(db)
     Y       X1       X2 X3
995  1 4.801836 20.82947  A
996  1 9.867854 24.39920  C
997  1 5.390730 21.25119  D
998  1 6.556160 20.79811  D
999  1 4.710276 21.15373  A
1000 1 6.631786 19.38083  A

Let us run a logistic regression on that dataset

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.45885    1.04646  -4.261 2.04e-05 ***
X1           0.51664    0.11178   4.622 3.80e-06 ***
X2           0.21008    0.07247   2.899 0.003745 ** 
X3B          1.74496    0.49952   3.493 0.000477 ***
X3C         -0.03470    0.35691  -0.097 0.922543    
X3D          0.08004    0.34916   0.229 0.818672    
X3E          2.21966    0.56475   3.930 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 397.69  on 993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

Here, the reference is modality . Which means that for someone with characteristics , we predict the following probability

where  denotes the cumulative distribution function of the logistic distribution

For someone with characteristics , we predict the following probability

For someone with characteristics , we predict the following probability

(etc.) Here, if we accept  (against ), it means that modality  cannot be considerd as different from .

A natural idea can be to change the reference modality, and to look at the -values. If we consider the following loop, we get

> M = matrix(NA,5,5)
> rownames(M)=colnames(M)=LETTERS[1:5]
> for(k in 1:5){
+ db$X3 = relevel(X3,LETTERS[k])
+ reg = glm(Y~X1+X2+X3,family=binomial,data=db)
+ M[levels(db$X3)[-1],k] = summary(reg)$coefficients[4:7,4]
+ } 
> M
             A            B            C            D            E
A           NA 0.0004771853 9.225428e-01 0.8186723647 8.482647e-05
B 4.771853e-04           NA 4.841204e-04 0.0009474491 4.743636e-01
C 9.225428e-01 0.0004841204           NA 0.7506242347 9.194193e-05
D 8.186724e-01 0.0009474491 7.506242e-01           NA 1.730589e-04
E 8.482647e-05 0.4743636442 9.194193e-05 0.0001730589           NA

and if we simply want to know if the -value exceeds – or not – 5%, we get the following,

> M.TF = M>.05
> M.TF
      A     B     C     D     E
A    NA FALSE  TRUE  TRUE FALSE
B FALSE    NA FALSE FALSE  TRUE
C  TRUE FALSE    NA  TRUE FALSE
D  TRUE FALSE  TRUE    NA FALSE
E FALSE  TRUE FALSE FALSE    NA

The first column is obtained when  is the reference, and then, we see which parameter should be considered as null. The interpretation is the following:

  •  and  are not different from 
  •  is not different from 
  •  and  are not different from 
  •  and  are not different from 
  •  is not different from 

Note that we only have, here, some kind of intuition. So, let us run a more formal test. Let us consider the following regression (we remove the intercept to get a model easier to understand)

> library(car)
> db$X3=relevel(X3,"A")
> reg=glm(Y~0+X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
    Estimate Std. Error z value Pr(>|z|)    
X1   0.51664    0.11178   4.622 3.80e-06 ***
X2   0.21008    0.07247   2.899  0.00374 ** 
X3A -4.45885    1.04646  -4.261 2.04e-05 ***
X3E -2.23919    1.06666  -2.099  0.03580 *  
X3D -4.37881    1.04887  -4.175 2.98e-05 ***
X3C -4.49355    1.06266  -4.229 2.35e-05 ***
X3B -2.71389    1.07274  -2.530  0.01141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1386.29  on 1000  degrees of freedom
Residual deviance:  397.69  on  993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

It is possible to use Fisher test to test if some coefficients are equal, or not (more generally if some linear constraints are satisfied)

> linearHypothesis(reg,c("X3A=X3C","X3A=X3D","X3B=X3E"))
Linear hypothesis test

Hypothesis:
X3A - X3C = 0
X3A - X3D = 0
- X3E  + X3B = 0

Model 1: restricted model
Model 2: Y ~ 0 + X1 + X2 + X3

  Res.Df Df  Chisq Pr(>Chisq)
1    996                     
2    993  3 0.6191      0.892

Here, we clearly accept the assumption that the first three factors are equal, as well as the last two. What is the next step? Well, if we believe that there are mainly two categories,  and , let us create that factor,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.39439    1.02791  -4.275 1.91e-05 ***
X1           0.51378    0.11138   4.613 3.97e-06 ***
X2           0.20807    0.07234   2.876  0.00402 ** 
X3bisBE      1.94905    0.36852   5.289 1.23e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 398.31  on 996  degrees of freedom
AIC: 406.31

Number of Fisher Scoring iterations: 7

Here, all the categories are significant. So we do have a proper model.

Comparing quantiles for two samples

Recently, for a research paper, I got some samples, and I wanted to compare them. Not to compare their means (by construction, all of them were centered) but there dispersion. And not their variance, but more their quantiles. Consider the following boxplot type function, where everything here is quantile related (which is not the case for standard boxplot, see http://freakonometrics.hypotheses.org/4138, in French)

> boxplotqbased=function(x){
+ q=quantile(x[is.na(x)==FALSE],c(.05,.25,.5,.75,.95))
+ plot(1,1,col="white",axes=FALSE,xlab="",ylab="",
+ xlim=range(X),ylim=c(1-.6,1+.6))
+ polygon(c(q[2],q[2],q[4],q[4]),1+c(-.4,.4,.4,-.4))
+ segments(q[1],1-.4,q[1],1+.4)
+ segments(q[5],1,q[4],1)
+ segments(q[5],1-.4,q[5],1+.4)
+ segments(q[1],1,q[2],1)
+ segments(q[3],1-.4,q[3],1+.4,lwd=2)
+ xt=x[(x<q[1])|(x>q[5])]
+ points(xt,rep(1,length(xt)))
+ axis(1)
+ }

(one can easily adapt the code for lists, e.g.). Consider for instance temperature, when the (linear) trend is removed (see http://freakonometrics.hypotheses.org/1016 for a discussion on that series, in Paris),

from January 1st till December 31st. Let us remove now the seasonal cycle, i.e. we do have here the difference with the average seasonal temperature (with here upper and lower quantiles),

Seasonal boxplots are here (with Autumn on top, then Summer, Spring and Winter, below),

If we zoom in, we do have (where upper and lower segments are 95% and 5% quantiles, while classically, boxes are related to the 75% and the 25% quantiles)

Is there a (standard) test to compare quantiles – some of them perhaps ? Can we compare easily quantiles when have two (or more) samples ?

Note that this example on temperature could be related to other old posts (see e.g. http://freakonometrics.hypotheses.org/2190), but the research paper was on a very different topic.

Consider two (i.i.d.) samples https://latex.codecogs.com/gif.latex?\{x_1,\cdots,x_m\} and https://latex.codecogs.com/gif.latex?\{y_1,\cdots,y_n\}, considered as realizations of random variables https://latex.codecogs.com/gif.latex?X and https://latex.codecogs.com/gif.latex?Y. In all statistical courses, tests on the average are always considered, i.e.

https://latex.codecogs.com/gif.latex?H_0:\mathbb{E}(X)=\mathbb{E}(Y)

against

https://latex.codecogs.com/gif.latex?H_1:\mathbb{E}(X)\neq\mathbb{E}(Y)

Usually, the idea in courses is to start with a one sample test, and to test something like

https://latex.codecogs.com/gif.latex?H_0:\mathbb{E}(X)=\mu_\star

against

https://latex.codecogs.com/gif.latex?H_1:\mathbb{E}(X)\neq\mu_\star

The idea is to assume that samples are from Gaussian variables,

https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T%20=%20\frac{\overline{x}%20-%20\mu_\star}{\widehat{\sigma}/\sqrt{n}}
Under https://latex.codecogs.com/gif.latex?H_0https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T has a Student t distribution. All that can be found in any Statistics 101 course. We can derive https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20pvalue, computing probabilities that https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T exceeds the observed values (for two sided tests, the probability that the absolute value of https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T exceed the absolute value of the observed statistics). This test is closely related to the construction of confidence intervals for https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20\mu. If https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20\mu_\star belongs to the confidence interval, then it might be a suitable value. The graphical representation of this test is related to the following graph

Here the observed value was 1,96, i.e. the https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20pvalue (the area in red above) is exactly 5%.

To compare means, the standard test is based on

https://latex.codecogs.com/gif.latex?T%20=%20{\overline{x}%20-%20\overline{y}%20\over%20%20\displaystyle\sqrt{{s_x^2%20\over%20m}%20+%20{s_y^2%20\over%20n}}%20}

which has – under https://latex.codecogs.com/gif.latex?H_0 – a Student-t distribution, with https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20\nu degrees of freedom, where

https://latex.codecogs.com/gif.latex?\nu%20=%20\frac{(s_x^2/m%20+%20s_y^2/n)^2}{(s_x^2/m)^2/(m-1)%20+%20(s_y^2/n)^2/(n-1)}.

Here, the graphical representation is the following,

But tests on quantiles are rarely considered in statistical courses. In a general setting,define quantiles as

https://latex.codecogs.com/gif.latex?Q_X(p)=\inf\left\{%20x\in%20\mathbb%20R%20:%20p%20\le%20\mathbb%20P(X\leq%20x)%20\right\}

one might be interested to test
https://latex.codecogs.com/gif.latex?H_0:Q_X(p)=Q_Y(p)
against
https://latex.codecogs.com/gif.latex?H_1:Q_X(p)\neq%20Q_Y(p)
for some https://latex.codecogs.com/gif.latex?p\in(0,1). Note that we might be interested also to test if

https://latex.codecogs.com/gif.latex?H_0:Q_X(p_k)=%20Q_Y(p_k)
for all https://latex.codecogs.com/gif.latex%20?%20%20k, for some vector of probabilities https://latex.codecogs.com/gif.latex?\boldsymbol{p}=(p_1,\cdots,p_d)\in(0,1)^d.
One can imagine that this multiple test will be more complex. But more interesting, e.g. a test on boxplots (are the four quantiles equal ?).  Let us start with something a bit more simple: a test on quantiles for one sameple, and the derivation of a confidence interval for quantiles.

  • Quantiles for one sample

The important idea here is that it should be extremely simple to get https://latex.codecogs.com/gif.latex?pvalues. Consider the following sample, and let us run a test to assess if the median can be zero.

> set.seed(1)
> X=rnorm(20)
> sort(X)
[1] -2.21469989 -0.83562861 -0.82046838 -0.62645381 -0.62124058 -0.30538839
[7] -0.04493361 -0.01619026  0.18364332  0.32950777  0.38984324  0.48742905
[13]  0.57578135  0.59390132  0.73832471  0.82122120  0.94383621  1.12493092
[19]  1.51178117  1.59528080
> sum(X<=0)
[1] 8

Here, 8 observations (out of 20, i.e. 40%) were below zero. But we do know the distribution of https://latex.codecogs.com/gif.latex%20?%20%20N the number of observation below the target

https://latex.codecogs.com/gif.latex?N=\sum_{i=1}^n%20\boldsymbol{1}(X_i\leq%20x_\star)

It is a binomial distribution. Under https://latex.codecogs.com/gif.latex?H_0, it is a binomial distribution https://latex.codecogs.com/gif.latex?\mathcal{B}(n,p_\star) where https://latex.codecogs.com/gif.latex?p_\star is the probability target (here 50% since the test is on the median). Thus, one can easily compute the https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20p-value,

> plot(n,dbinom(n,size=20,prob=0.50),type="s",xlab="",ylab="",col="white")
> abline(v=sum(X<=0),col="red")
> for(i in 1:sum(X<=0)){
+ polygon(c(n[i],n[i],n[i+1],n[i+1]),
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ polygon(21-c(n[i],n[i],n[i+1],n[i+1]),
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ }
> lines(n,dbinom(n,size=20,prob=0.50),type="s")

which yields

Here, the https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20p-value is

> 2*pbinom(sum(X<=0),20,.5)
[1] 0.5034447

Here the probability is easy to compute. But one can observe that there is some kind of disymmetry here. Actually, if the observed value was not 8, but 12, some minor changes should be done (to keep some symmetry),

> plot(n,dbinom(n,size=20,prob=0.50),type="s",xlab="",ylab="",col="grey")
> abline(v=20-sum(X<=0),col="red")
> for(i in 1:sum(X<=0)){
+ polygon(c(n[i],n[i],n[i+1],n[i+1])-1,
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ polygon(21-c(n[i],n[i],n[i+1],n[i+1])-1,
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ }
> lines(n-1,dbinom(n,size=20,prob=0.50),type="s")

Based on those observations, one can easily write a code to test if the https://latex.codecogs.com/gif.latex?p_\star-quantile of a sample is https://latex.codecogs.com/gif.latex?x_\star. Or not. For a two sided test, consider

> quantile.test=function(x,xstar=0,pstar=.5){
+ n=length(x)
+ T1=sum(x<=xstar)
+ T2=sum(x< xstar)
+ p.value=2*min(1-pbinom(T2-1,n,pstar),pbinom(T1,n,pstar))
+ return(p.value)}

Here, we have

> quantile.test(X)
[1] 0.5034447

Now, based on that idea, due to the duality between confidence intervals and tests, one can easily write a function that computes confidence interval for quantiles,

> quantile.interval=function(x,pstar=.5,conf.level=.95){
+ n=length(x)
+ alpha=1-conf.level
+ r=qbinom(alpha/2,n,pstar)
+ alpha1=pbinom(r-1,n,pstar)
+ s=qbinom(1-alpha/2,n,pstar)+1
+ alpha2=1-pbinom(s-1,n,pstar)
+ c.lower=sort(x)[r]
+ c.upper=sort(x)[s]
+ conf.level=1-alpha1-alpha2
+ return(list(interval=c(c.lower,c.upper),confidence=conf.level))}
> quantile.interval(X,.50,.95)
$interval
[1] -0.3053884  0.7383247

$confidence
[1] 0.9586105

Because of the use of non-asymptotic distributions, we can not get exactly a 95% confidence interval. But it is not that bad, here.

  • Comparing quantiles for two samples

Now, to compare quantiles for two samples… it is more complicated. Exact tests are discussed in Kosorok (1999) (see http://bios.unc.edu/~kosorok/…) or in Li, Tiwari and Wells (1996) (see http://jstor.org/…). For the computational aspects, as mentioned in a post published almost one year ago on http://nicebread.de/… there is a function to compare quantiles for two samples.

> install.packages("WRS")
> library("WRS")

Some multiple tests on quantiles can be performed here. For instance, on the temperature, if we compare quantiles for Winter and Summer (on only 1,000 observations since it can be long to run that function), i.e. 5%, 25%, 75% and 95%,

> qcomhd(Z1[1:1000],Z2[1:1000],q=c(.05,.25,.75,.95))
q   n1   n2      est.1      est.2 est.1_minus_est.2     ci.low     ci.up     p_crit p.value signif
1 0.05 1000 1000 -6.9414084 -6.3312131       -0.61019530 -1.6061097 0.3599339 0.01250000   0.220     NO
2 0.25 1000 1000 -3.3893867 -3.1629541       -0.22643261 -0.6123292 0.2085305 0.01666667   0.322     NO
3 0.75 1000 1000  0.5832394  0.7324498       -0.14921041 -0.4606231 0.1689775 0.02500000   0.338     NO
4 0.95 1000 1000  3.7026388  3.6669997        0.03563914 -0.5078507 0.6067754 0.05000000   0.881     NO

or if we compare quantiles for Winter and Summer

> qcomhd(Z1[1:1000],Z3[1:1000],q=c(.05,.25,.75,.95))
q   n1  n2      est.1     est.2 est.1_minus_est.2     ci.low       ci.up     p_crit p.value signif
1 0.05 1000 984 -6.9414084 -6.438318        -0.5030906 -1.3748624  0.39391035 0.02500000   0.278     NO
2 0.25 1000 984 -3.3893867 -3.073818        -0.3155683 -0.7359727  0.06766466 0.01666667   0.103     NO
3 0.75 1000 984  0.5832394  1.010454        -0.4272150 -0.7222362 -0.11997409 0.01250000   0.012    YES
4 0.95 1000 984  3.7026388  3.873347        -0.1707078 -0.7726564  0.37160846 0.05000000   0.539     NO

(the following graphs are then plotted)

Those tests are based on the procedure proposed in Wilcox, Erceg-Hurn,  Clark and Carlson (2013), online on http://tandfonline.com/…. They rely on the use of bootstrap samples. The idea is quite simple actually (even if, in the paper, they use Harrell–Davis estimator to estimate quantiles, i.e. a weighted sum of ordered statistics – as described in http://freakonometrics.hypotheses.org/1755 – but the idea can be understood with any estimator): we generate several bootstrap samples, and compute the median for all of them (since our interest was initially on the median)

>  Q=rep(NA,10000)
>  for(b in 1:10000){
+  Q[b]=quantile(sample(X,size=20,replace=TRUE),.50)
+  }

Then, to derive a confidence interval (with, say, 95% confidence), we compute quantiles of those median estimates,

> quantile(Q,c(.025,.975))
     2.5%     97.5% 
-0.175161  0.666113

We can actually visualize the distribution of that bootstrap median,

> hist(Q)

Now, if we want to compare medians from two independent samples, the strategy is rather similar: we bootstrap the two samples – independently – then compute the median, and keep in mind the difference. Then, we will look if the difference is significantly different from 0. E.g.

> set.seed(2)
> Y=rnorm(50,.6)
> QX=QY=D=rep(NA,10000)
> for(b in 1:10000){
+ QX[b]=quantile(sample(X,size=length(X),replace=TRUE),.50)
+ QY[b]=quantile(sample(Y,size=length(Y),replace=TRUE),.50)
+ D[b]=QY[b]-QX[b]
+ }

The 95% confidence interval obtained from the bootstrap difference is

> quantile(D,c(.025,.975))
      2.5%      97.5% 
-0.2248471  0.9204888

which is rather close to was can be obtained with the R function

> qcomhd(X,Y,q=.5)
    q n1 n2    est.1     est.2 est.1_minus_est.2    ci.low     ci.up p_crit p.value signif
1 0.5 20 50 0.318022 0.5958735        -0.2778515 -0.923871 0.1843839   0.05    0.27     NO

(where the difference is here the oppositive of mine). And when testing for 2 (or more) quantiles, Bonferroni method can be used to take into account that those tests cannot be considered as independent.

Test, valeur critique et p-value

Un petit complément suite au cours de mercredi dernier, pour insister sur l’importance de la p-value dans la lecture de la sortie d’un test.

  • Les erreurs dans un test statistique

Mais avant, rappelons qu’un test est une prise de décision: accepter ou rejeter une hypothèse. Et qu’on peut commettre une erreur. Ou pour être plus précis, on peut commettre deux types d’erreur,
• accepter l’hypothèse alors que cette dernière est fausse
• rejeter l’hypothèse alors que cette dernière était vraie
Pour reprendre une terminologie plus médicale, un test de grossesse peut dire à une femme qu’elle n’est pas enceinte, alors qu’elle l’est; ou dire qu’elle l’est, alors qu’elle ne l’est pas (voir tous les exemples dans les exercices de probabilités de l’examen P de la SOA, ou le cours ACT2121).
Formellement, on a deux probabilités,
• la probabilité d’accepter à tort notre hypothèse (on parlera d’erreur de second espèce), \beta
• la probabilité de rejeter à tort notre hypothèse (on parlera d’erreur de première espèce) \alpha
Dans un monde idéal on voudrait que les deux probabilités soient aussi petites que possibles… Mais c’est impossible, et le plus souvent, baisser une des probabilités se fait en augmentant l’autre. Les cas extrêmes étant
• avoir un test de grossesse qui déclare tout le monde enceinte: on ne rejette alors jamais à tort (on ne rejette jamais tout court en fait), mais on a un fort taux d’acceptation à tort,
• avoir un test de grossesse qui ne déclare personne enceinte: on n’accepte jamais à tort (car on n’accepte jamais) mais on a un fort taux de rejet à tort.
Bref, on a un arbitrage à faire entre deux types d’erreurs. Souvent, en pratique on va demander à contrôler l’erreur de première espèce (i.e. \alpha de l’ordre de 5%), et on chercher a un test qui, à \alpha donné, possède la plus faible erreur de première espèce. Voilà en gros pour la théorie: on se donne un seuil de significativité \alpha, qui correspond à la probabilité d’erreur de premier type. Et on va chercher à tester si une hypothèse H_0 est vraie, l’alternative étant une hypothèse H_1.

H_0 vraie H_1 vraie
accepter H_0 OK erreur
type 2
rejeter H_0 erreur
type 1
OK
  • La “valeur critique”

La notion de valeur critique a été introduite dans Neyman & Pearson (1928). Cette valeur dépend de la forme de l’hypothèse alternative, en particulier savoir si le test est bilatéral, unilatéral à gauche, ou unilatéral à droite. Pour un test donné, la valeur critique peut-être vue comme la valeur limite a partir de laquelle on pourra rejeter H_0 avec un seuil de significativité donné.

  • La p-value

La p-value a été introduite dans Gibbons & Pratt (1975), meme si on peut retrouve l’idée beaucoup plus tôt, comme Pearson (1900), qui propose de calculer “the probability that the observed value of the chi-square statistic would be exceeded under the null hypothesis“. La p-value est la probabilité, sous H_0, d’obtenir une statistique aussi extrême (pour ne pas dire aussi grande) que la valeur observée sur l’échantillon. Aussi, pour un seuil de significativité \alpha donné, on compare p et \alpha, afin d’accepter, ou de rejeter H_0,
• si p\leq\alpha, on va rejeter l’hypothèse H_0 (en faveur de H_1)
• si p>\alpha, on va rejeter H_1 (en faveur de H_0).
On peut alors interpréter la p-value comme le plus petit seuil de significativité pour lequel l’hypothèse nulle est acceptée. Gibbons & Pratt (1975) reviennent longuement sur les interprétations, et surtout les mauvaises interprétations, de cette p-value.

  • Valeur critique versus p-value

Si on formalise un peu, on peut vouloir tester H_0:\theta=\theta_0 contre H_1:\theta>theta_0 (par exemple). De manière très générale, on dispose d’une statistique de test T qui a pour loi, sous H_0, F_{\theta_0}(\cdot) (que l’on supposera continue). Notons qu’on peut considérer une hypothèse alternative de la forme H_1:\theta\neq\theta_0, c’est juste plus pénible parce qu’il faut travailler sur \vert T\vert, et calculer des probabilités à gauche, ou à droite. Donc pour notre exemple, on va prendre un test unilatéral.
Dans l’approche classique (telle que présentée dans tous les cours de statistiques), on se donne un seul d’acceptation \alpha petit (disons 5%), et on cherche une valeur critique T_{1-alpha} telle que

Pour ceux qui se souviennent de leur cours de stats, cela peut faire penser à la puissance du test, définie par

\pi(\theta\vert \alpha)=\mathbb{P}(T\geq T_{1-\alpha}\vert \theta)=1-F_{\theta}(T_{1-\alpha})

Formellement, la p-value associée au test T est la variable aléatoire P définie par
P=1-F_{\theta_0}(T).
Donc effectivement, la p-value et la puissance sont liées, puisque

\mathbb{P}(P\leq \alpha\vert \theta)=\pi(\theta\vert \alpha)

autrement dit, la puissance peut-être vue comme la fonction de répartition de la p-value.

  • Intérêt computationnel de la p-value

D’un point de vue computationnel, la p-value est l’outil le plus important pour interpréter la sortie d’un test. Commençons par un test simple, comme une comparaison de moyennes. On cherche ici à tester H_0:\mu_X=\mu_Y contre H_1:\mu_X>\mu_Y pour des moyennes calculées sur deux groupes. Pour reprendre l’exemple abordé dans un précédant billet, on a les notes obtenues en ACT6420 par deux groupes différents. Et on veut savoir s’ils sont vraiment différents (ci-dessous le nombre de bonnes réponses, sur 40 questions, on travaillera ensuite sur la note sur 100)

image manquante

La statistique de test est ici

T = \frac{\overline{X} - \overline{Y}}{\displaystyle{ \sqrt{ {s_X^2 \over n_X} + {s_Y^2 \over n_Y} }}}

et sous H_0, T va suivre une loi de Student à \nu degrés de liberté, où \nu est donné par la relation de Welch–Satterthwaite (d’après Satterwaite (1946) et Welch (1947)),

\nu = {{\left( {s_X^2 \over n_X} + {s_Y^2 \over n_Y}\right)^2 } \over {{s_X^4 \over n_X^2 \cdot \left({n_X-1}\right)}+{s_Y^4 \over n_Y^2 \cdot \left({n_Y-1}\right)}}}

Numériquement, on a ici

> Xbar=mean(X)
> Ybar=mean(Y)
> Sx2=var(X)
> Sy2=var(Y)
> nX=length(X)
> nY=length(Y)
> (T=(Xbar-Ybar)/sqrt(Sx2/nX+Sy2/nY))
[1] -2.155754

et pour les degrés de liberté

> (nu=(Sx2/nX+Sy2/nY)^2/(Sx2^2/nX^2/(nX-1)+
+ Sy2^2/nY^2/(nY-1)))
[1] 36.35279

La valeur critique est obtenue en lisant dans les tables,

(car ici on a des probabilité pour un test bilatéral dans la table) comme on apprenait dans les cours de statistique au siècle passé. D’un point de vue informatique, on cherche à savoir si on est à gauche, ou à droite de la valeur critique

> qt(.05,df=nu)
[1] -1.687865

image manquante

On peut aussi calculer la p-value,

> pt(T,df=nu)
[1] 0.01889768

Si on regarde, sous R, il existe des fonctions de tests, pour comparer des moyennes. Et dans ce cas, la sortie est

> t.test(X,Y,alternative = "less")

Welch Two Sample t-test

data:  X and Y
t = -2.1558, df = 36.353, p-value = 0.0189
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf -1.772507
sample estimates:
mean of x mean of y
48.75000  56.91667

Autrement dit, on a automatiquement la p-value, et qui permet rapidement d’interpréter le test. Moralité, si on sait interpréter une p-value (et que l’on vérifié au préalable les conditions d’application d’un test), on peut faire tous les tests que l’on veut !
Si on veut faire un peu plus compliqué, on peut regarder la distribution des notes, et se demander si une loi \mathcal{N}(60,15^2) serait possible (par exemple, ça sera notre hypothèse H_0, l’hypothèse alternative étant que ce n’est pas cette loi). Pour faire ce test, il existe le test de Kolmogorov-Smirnov. La statistique de test est ici

T=\sup\{\vert \widehat{F}_n(x)-F_0(x)\vert ,x\in\mathbb{R}\}

F_0(\cdot) est la fonction de répartition de la loi \mathcal{N}(60,15^2), et \widehat{F}_n(\cdot) est la fonction de répartition empirique

\widehat{F}_n(x)=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(x_i\leq x)

La loi de T n’est pas simple, ou moins simple qu’une loi de Student (cf Marsaglia, Tsang & Wang (2003) par exemple). En revanche, on a les p-values automatiquement,

> ks.test(Y, "pnorm", 60, 15)

One-sample Kolmogorov-Smirnov test

data:  Y
D = 0.1421, p-value = 0.5796
alternative hypothesis: two-sided

Aussi, on peut accepter ici l’hypothèse nulle. On peut d’ailleurs faire un petit dessin pour s’en convaincre,

> Femp=function(x) mean(Y<=x)
> plot(0:100,Vectorize(Femp)(0:100),type="s")
> lines(0:100,pnorm(0:100,60,15),col="red")

image manquante


Et ça va nous servir dans ce cours ? A priori oui… parce qu’on parlera du test de Student (pour tester si une variable dans une régression est significative), du test de Fisher (pour tester si plusieurs variables dans une régression sont significatives, ou plus généralement si une contrainte – linéaire – sur les coefficients peut être acceptée), du test de Chow (pour tester des ruptures dans un modèle linéaire, mais c’est un test de Fisher un peu déguisé), du test d’Anderson-Darling (pour tester si des résidus sont Gaussiens), du test de Breuch-Pagan voire le test de White (pour tester si les résidus peuvent être considérés de variance constante), du test de Durbin-Watson (pour tester s’il n’y a pas d’auto-corrélation dans la série des résidus), du test de Dickey-Fuller (pour tester si une série temporelle est – ou n’est pas – stationnaire), des tests de Franses (pour tester si une série peut être considérée comme saisonnière, ou pas), du test de Ljung-Box (pour tester si un bruit est un bruit blanc)… Et j’en oublie un paquet. Donc quand il est dit (dans le plan de cours) que le cours de statistique est un prérequis, il ne s’agit pas de l’avoir suivi, mais bel et bien de l’avoir compris, car on passera notre temps à utiliser des notions entrevues dans ce cours.

We keep breaking records ? so what ?… Get statistical perspective….

This summer, we have been told that some financial series broke some records (here, in French)

For instance, the French CAC40 had negative return for 11 consecutive days (which has never been seen, so far).

> library(tseries)
> x<-get.hist.quote("^FCHI")
> Y=x$Close
> Z=diff(log(Y))
> RUN=rle(as.character(Z>=0))$lengths
> n=length(RUN)
> LOSS=RUN[seq(2,n,by=2)]
> GAIN=RUN[seq(1,n,by=2)]
> TG=sort(table(GIN))
> TG[as.character(1:13)]
GAIN
   1    2    3    4    5    6    7    8    9 <NA> <NA> <NA>   13
 645  336  170   72   63   21    7    3    4   NA   NA   NA    1
 
> TL=sort(table(LOSS))
> TL[as.character(1:15)]
LOSS
   1    2    3    4    5    6    7    8    9 <NA>   11 <NA> <NA> 
 664  337  186   68   42   14    5    3    1   NA    1   NA   NA 
 
> TR=sort(table(RUN))
> TR[as.character(1:15)]
RUN
   1    2    3    4    5    6    7    8    9 <NA>   11 <NA>   13 
1309  673  356  140  105   35   12    6    5   NA    1   NA    1

Indeed 11 consecutive days of negative returns is a record. But one should keep in mind the fact that the real records for runs is 13 consecutive days with positivereturns…
But what does that mean ? Can we still assume time independence of log-returns (since today, a lot of financial models are still based on that assumption) ?
Actually. if financial series were time-independence, such a probability, indeed, should be rather small. At least on 11 or 10 runs. Something like

http://freakonometrics.blog.free.fr/public/perso3/cacindp03.gif

(assuming that each day, the probability to observe a negative return is 50%). But maybe not over 25 years (6250 trading days): the probability to observe a sub-sequence of 10 consecutive negative value (with daily probability of one half) over 6250 observations will be much larger. My guess is that is would be

http://freakonometrics.blog.free.fr/public/perso3/cacindp02.gif

where at the numerator we have the number of favourable cases over the total number of cases. At the numerator, the first number the number of cases where the first 10 (at least) are negative, then for the second one, we count the number of cases where the first is positive, then the next 10 (at least) are negative (and then the second is positive and then the next 10 are negative, the third is positive etc). For those interested by more details (and a more general formula on runs), an answer can be found here.
But note that the probability is quite large… So it is not that unlikely to observe such a sequence over 25 years.
A classical idea when looking at time series is to look at the autocorrelation function of the returns,

which might suggest that there is no correlation with past returns. But it should be possible to do more advanced tests.
On the CAC40 series, we can run an independence run test on the latest 100 consecutive days, and look at the p-value,

> library(lawstat)
> u=as.vector(Z[(n-100):n])
> runs.test(u,plot=TRUE)
 
	Runs Test - Two sided
 
data:  u 
Standardized Runs Statistic = -0.4991, p-value = 0.6177

The B’s here are returns lower than the median (almost null, so they might be considered as negative returns). With such a high p-value, we accept the null hypothesis, i.e. time independence.
If we consider a moving-time window

we can see that we accept the assumption of independence most the the time.
Actually, here, the time window is 100 days (+/- 50 days). But it is possible to consider 200 days,

or even 400 days,

So, except if we focus on 2006, it looks like we should reject the idea of time dependence in financial markets.
It is also possible to look more carefully at the distribution of runs, and to compare it with the case of independent samples (here we consider monte carlo generation of sequences having the same size),

> m=length(Z)
> ns=100000
> HIST=matrix(NA,ns,15)
> for(j in 1:ns){
+ XX=sample(c("A","B"),size=m,replace=TRUE)
+ RUNX=rle(as.character(XX))$lengths
+ S=sort(table(RUNX))
+ HIST[j,]=S[as.character(1:15)]
+ }
> meana=function(x){sum(x[is.na(x)==FALSE])/length(x)}
> cbind(TR[as.character(1:15)],apply(HIST,2,meana),
+       round(m/(2^(1+1:15))))
 
     [,1]       [,2] [,3]
1    1309 1305.12304 1305
2     673  652.46513  652
3     356  326.21119  326
4     140  163.05101  163
5     105   81.52366   82
6      35   40.74539   41
7      12   20.38198   20
8       6   10.16383   10
9       5    5.09871    5
10     NA    2.56239    3
11      1    1.26939    1
12     NA    0.63731    1
13      1    0.31815    0
14     NA    0.15812    0
15     NA    0.08013    0

The first column above is the empirical frequency of runs of length 1,2,3, etc. The second one is the average frequencies obtained on random simulation of independent sample. The third one is the theoretical frequency based on a (geometric distribution with mean 1).

Here again, it looks like our time series behave like an independent sample. Here is also a nice paper by Mark Schilling on the longest run of heads.
So it is not that odd to observe such a series of losses on financial markets….

Cochrane, Pearson et le test du chi-deux

En cours, nous avons poursuivi sur la loi multinomiale, et le test du chi-deux. Je voulais mettre un petit billet pour récapituler les différents points, et montrer une application numérique (nous reviendrons en détails mercredi sur des applications des outils vus jusqu’à présent).

  • Inférence avec la loi multinomiale

On suppose qu’une variable http://freakonometrics.blog.free.fr/public/maths/coch-01.gif peut prendre http://freakonometrics.blog.free.fr/public/maths/coch-02.gif modalités, notées http://freakonometrics.blog.free.fr/public/maths/coch-03.gif, avec http://freakonometrics.blog.free.fr/public/maths/coch-04.gif. On posera

http://freakonometrics.blog.free.fr/public/maths/coch-05.gif

en notant que http://freakonometrics.blog.free.fr/public/maths/coch-06.gif appartient au simplexe de http://freakonometrics.blog.free.fr/public/maths/coch-07.gif au sens où

http://freakonometrics.blog.free.fr/public/maths/coch-08.gif

On a vu que l’estimateur du maximum de vraisemblance s’obtenait en faisant un peu d’optimisation sous contrainte, et que

http://freakonometrics.blog.free.fr/public/maths/coch-10.gif

(en reprenant les notations du cours). On avait montré que

http://freakonometrics.blog.free.fr/public/maths/coch-11.gif

et on a vu

http://freakonometrics.blog.free.fr/public/maths/coch-13.gif

(ce qui peut se retrouver en introduisant la variable binomiale http://freakonometrics.blog.free.fr/public/maths/coch-16.gif). Mais plus généralement, on finira les calculs permettant d’établir que, pour http://freakonometrics.blog.free.fr/public/maths/coch-17.gif

http://freakonometrics.blog.free.fr/public/maths/coch-18.gif

(ce qui permet d’obtenir la matrice de variance covariance de http://freakonometrics.blog.free.fr/public/maths/coch-20.gif).
En utilisant le théorème central limite on peut montrer que

http://freakonometrics.blog.free.fr/public/maths/coch-23.gif

Sous une forme multivariée, on écrira

http://freakonometrics.blog.free.fr/public/maths/coch-25.gif
http://freakonometrics.blog.free.fr/public/maths/coch-26.gif avec http://freakonometrics.blog.free.fr/public/maths/coch-27.gif et pour http://freakonometrics.blog.free.fr/public/maths/coch-17.gif, http://freakonometrics.blog.free.fr/public/maths/coch-28.gif.

On notera que la somme de la ième colonne de http://freakonometrics.blog.free.fr/public/maths/coch-29.gif est http://freakonometrics.blog.free.fr/public/maths/coch-30.gif.
Aussi, http://freakonometrics.blog.free.fr/public/maths/coch-29.gif n’est pas inversible (c’est le fait que notre paramètre appartient au simplexe).
Pour s’en sortir, la première idée est d’utiliser un peu d’algèbre linéaire. Une matrice http://freakonometrics.blog.free.fr/public/maths/coch-31.gif est une matrice de projection si elle est idempotente, i.e. http://freakonometrics.blog.free.fr/public/maths/coch-32.gif. Ses valeurs propres sont alors 0 ou 1, et si http://freakonometrics.blog.free.fr/public/maths/coch-35.gif est le nombre de fois où 1 est valeur propre, et si http://freakonometrics.blog.free.fr/public/maths/coch-36.gif, alors http://freakonometrics.blog.free.fr/public/maths/coch-37.gif.
Posons http://freakonometrics.blog.free.fr/public/maths/coch-38.gif. Alors

http://freakonometrics.blog.free.fr/public/maths/coch-39.gif

Or compte tenu de la forme de http://freakonometrics.blog.free.fr/public/maths/coch-29.gif,

http://freakonometrics.blog.free.fr/public/maths/coch-40.gif

qui est une matrice de projection dont la trace est http://freakonometrics.blog.free.fr/public/maths/coch-41.gif (qui est aussi le nombre de fois où 1 est valeur propre). Donc

http://freakonometrics.blog.free.fr/public/maths/coch-42.gif

Le test du chi-deux sera basé sur le fait qu’asymmptotiquement

http://freakonometrics.blog.free.fr/public/maths/coch-44.gif

Une autre idée consiste à construire http://freakonometrics.blog.free.fr/public/maths/coch-41.gif variables aléatoires qui seront indépendantes. Mais on peut plutôt regarder les applications de ce test, en particulier comme test d’indépendance.
Pour information, Frank Yates a proposé un correction “pour continuité“, ici. La statistique considérée est alors

http://fre<br /><br /><br /><br />
akonometrics.blog.free.fr/public/maths/coch-46.gif
  • Retour sur le théorème de Cochrane

Soit http://freakonometrics.blog.free.fr/public/maths/coch-50.gif de dimension http://freakonometrics.blog.free.fr/public/maths/coch-51.gif. Posons http://freakonometrics.blog.free.fr/public/maths/coch-59.gif, pour http://freakonometrics.blog.free.fr/public/maths/coch-04.gif, où on notera http://freakonometrics.blog.free.fr/public/maths/coch-60.gif le rang de http://freakonometrics.blog.free.fr/public/maths/coch-62.gif, en supposant que les http://freakonometrics.blog.free.fr/public/maths/coch-62.gif sont positives semidefinies, alors on a équivalence entre

  •  http://freakonometrics.blog.free.fr/public/maths/coch-63.gif
  • http://freakonometrics.blog.free.fr/public/maths/coch-64.gif pour http://freakonometrics.blog.free.fr/public/maths/coch-04.gif,
  • les http://freakonometrics.blog.free.fr/public/maths/coch-65.gif sont des variables indépendantes.

Les http://freakonometrics.blog.free.fr/public/maths/coch-65.gif s’interprètent comme des longueurs (euclidienne) de projections d’un vecteur Gaussien sur des sous-espaces orthogonaux (de dimension respective http://freakonometrics.blog.free.fr/public/maths/coch-60.gif). Si ces vecteurs sont indépendants, et suivent des lois du chi-deux à http://freakonometrics.blog.free.fr/public/maths/coch-60.gif degrés de libertés, avec http://freakonometrics.blog.free.fr/public/maths/coch-63.gif, alors les sous-espaces sont orthogonaux, et supplémentaires. On peut y voir une espèce d’extension du théorème de Pythagore, en remplaçant la notion de vecteurs orthogonaux par des variables indépendantes suivant des lois du chi-deux, et la somme des carrés des longueurs devient la somme des degrés de liberté.

  • Application comme test d’indépendance

Considérons deux variables http://freakonometrics.blog.free.fr/public/maths/coch-66.gif pouvant prendre toutes les deux deux modalités (disons deux lois binomiales). On est alors face a une loi multinomiale à 4 modalités

  • http://freakonometrics.blog.free.fr/public/maths/coch-79.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-70.gif
  • http://freakonometrics.blog.free.fr/public/maths/coch-78.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-73.gif
  • http://freakonometrics.blog.free.fr/public/maths/10gif.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-74.gif
  • http://freakonometrics.blog.free.fr/public/maths/coch-77.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-75.gif

Un test d’indépendance revient à tester si la loi multinomiale peut s’écrire

http://freakonometrics.blog.free.fr/public/maths/chi2-ab.gif
http://freakonometrics.blog.free.fr/public/maths/chi_ab2.gif
http://freakonometrics.blog.free.fr/public/maths/chi_ab3.gif
http://freakonometrics.blog.free.fr/public/maths/chi_ab4.gif

pour des vecteurs http://freakonometrics.blog.free.fr/public/maths/chi-a.gif et http://freakonometrics.blog.free.fr/public/maths/chi-b.gif. tels que http://freakonometrics.blog.free.fr/public/maths/chi-a2.gif et http://freakonometrics.blog.free.fr/public/maths/chi-b2.gif. On a alors http://freakonometrics.blog.free.fr/public/maths/chi212121.gif contraintes sur les paramètres. Ces deux contraintes additionnelles font que la statistique de test s’écrit

http://freakonometrics.blog.free.fr/public/maths/CHI-INDEP.gif

qui va suivre asymptotiquement une loi http://freakonometrics.blog.free.fr/public/maths/CHI1.gif i.e. http://freakonometrics.blog.free.fr/public/maths/CHI12.gif d’après le théorème de Cochrane.

  • Peine de mort pour les condamnées pour meurtre en Floride 1976-1987

en fonction de la “race” du meurtrier et de la victime,

  • meurtrier de “race blanche” et victime de “race blanche“: 53 condamnés à mort, 414 non condamnés à mort
  • meurtrier de “race blanche” et victime de “race noire“: 0 condamné à mort, 16 non condamnés à mort
  • meurtrier de “race noire“et victime de “race blanche“: 11 condamnés à mort, 37 non condamnés à mort
  • meurtrier de “race noire“et victime de “race noire“: 4 condamnés à mort, 139 non condamnés à mort

On peut alors faire des tests d’indépendance, entre la “race” du meurtrier et le verdict par exemple.

MEURTRIER=matrix(c(53+0,11+4,414+16,139+37),2,2)
VICTIME  =matrix(c(53+11,0+4,414+37,139+16),2,2)
n=sum(MEURTRIER)
(PROBMEUTR=MEURTRIER/n)
           [,1]      [,2]
[1,] 0.07863501 0.6379822
[2,] 0.02225519 0.2611276

SL=rowSums(PROBMEUTR)
SC=colSums(PROBMEUTR)
(MEUTRINDEP=outer(SL, SC, "*"))
           [,1]      [,2]
[1,] 0.07229966 0.6443176
[2,] 0.02859055 0.2547922

(Q=n*sum((PROBMEUTR - MEUTRINDEP)^2/MEUTRINDEP))
[1] 1.468519

(Qcorrec=n*sum((abs(PROBMEUTR - MEUTRINDEP)-.5/n)^2/MEUTRINDEP))
[1] 1.144741

pchisq(Qcorrec, (2-1)*(2-1), lower.tail
 = FALSE)
[1] 0.2846528

qchisq(.95, (2-1)*(2-1))
[1] 3.841459

chisq.test(MEURTRIER)

Pearson's Chi-squared test with Yates' continuity correction

data:  MEURTRIER 
X-squared = 1.1447, df = 1, p-value = 0.2847

On rejette donc l’hypothèse d’indépendance.

Want to say one thing and the exact oppositive with strong confidence ?

No need to do politics. Just take a statistical course. And I do not talk about misinterpretation of statistics, but I talk about the mathematical foundations of statistical tests.
Consider the following parametric test, with a one-dimensional parameter: http://freakonometrics.blog.free.fr/public/perso2/test-lies-01.gif versus http://freakonometrics.blog.free.fr/public/perso2/test-lies-02.gif, for some (fixed) http://freakonometrics.blog.free.fr/public/perso2/test-lies-03.gif. A standard way of doing such a test is to consider an rejection region http://freakonometrics.blog.free.fr/public/perso2/test-lies-05.gif. The test works as follows: consider a sample http://freakonometrics.blog.free.fr/public/perso2/test-lies-06.gif,

  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-07.gif, then we accept http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-09.gif, the we reject http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif

For instance, consider the case of a Bernoulli sample, with probability http://freakonometrics.blog.free.fr/public/perso2/test-lies-62.gif. The standard idea is to define

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

The rejection region is then based on statistic http://freakonometrics.blog.free.fr/public/perso2/test-lies-210.gif,

  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-25.gif, then we accept http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-22.gif, the we reject http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif

where threshold http://freakonometrics.blog.free.fr/public/perso2/test-lies-26.gif is taken so that the probability to make a first type error is http://freakonometrics.blog.free.fr/public/perso2/test-lies-28.gif(say 5%) using the Gaussian approximation for z. Here

http://freakonometrics.blog.free.fr/public/perso2/test-lies-30.gif

Thus, the acceptation region is then the green area below, while the rejection region is the red one, for http://freakonometrics.blog.free.fr/public/perso2/test-lies-210.gif.

Consider now the exact opposite test (with the same http://freakonometrics.blog.free.fr/public/perso2/test-lies-03.gif), http://freakonometrics.blog.free.fr/public/perso2/test-lies-51.gifversus http://freakonometrics.blog.free.fr/public/perso2/test-lies-52.gif. Here, we use the same statistics, and the test is

  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-22.gif, then we accept http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-25.gif, the we reject http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif

where now

http://freakonometrics.blog.free.fr/public/perso2/test-lies-50.gif

Thus, now, the acceptation region is then the green area below, while the rejection region is the red one.

So if we summarize what we just said,

  • in the region on the left below, both test agree that http://freakonometrics.blog.free.fr/public/perso2/test-lies-55.gif
  • in the region on the right below, both test agree that http://freakonometrics.blog.free.fr/public/perso2/test-lies-57.gif
  • and in the region in blue, in the middle, the two tests disagree (one claims that http://freakonometrics.blog.free.fr/public/perso2/test-lies-55.gif, and the other one that http://freakonometrics.blog.free.fr/public/perso2/test-lies-57.gif)

Here is the evolution of the region as a function of http://freakonometrics.blog.free.fr/public/perso2/test-lies-56.gif (the size of the sample) when the sample frequency is 20%. With a small sample size, we can hardly say anything.

n=seq(1,100)
p=0.2
x1=p+qnorm(.95)*sqrt(p*(1-p)/n)
x2=p+qnorm(.05)*sqrt(p*(1-p)/n)
plot(n,x1,type="l",ylim=c(0,1))
polygon(c(n,rev(n)),c(x1,rev(x2)),col="light blue",border=NA)
lines(n,x1,lwd=2,col="red")
lines(n,x2,lwd=2,col="red")

One might say that those bounds are based on a Gaussian approximation which is not correct when http://freakonometrics.blog.free.fr/public/perso2/test-lies-56.gif is too small. So we can compute exact bounds,
y1=qbinom(.95,size=n,prob=p)/n
y2=qbinom(.05,size=n,prob=p)/n
polygon(c(n,rev(n)),c(y1,rev(y2)),col="blue",border=NA)
lines(n,y1,lwd=2,col="red")
lines(n,y2,lwd=2,col="red")

and we get

This is what we can observe if we use R statistical procedures, either the asymptotic one,

> prop.test(2,10,.5,alternative="less")
 
1-sample proportions test with continuity correction
 
data:  2 out of 10, null probability 0.5
X-squared = 2.5, df = 1, p-value = 0.05692
alternative hypothesis: true p is less than 0.5
95 percent confidence interval:
0.0000000 0.5100219
sample estimates:
p
0.2
 
> prop.test(2,10,.5,alternative="greater")
 
1-sample proportions test with continuity correction
 
data:  2 out of 10, null probability 0.5
X-squared = 2.5, df = 1, p-value = 0.943
alternative hypothesis: true p is greater than 0.5
95 percent confidence interval:
0.04368507 1.00000000
sample estimates:
p
0.2

or a more accurate one

> binom.test(2,10,.5,alternative="less")
 
Exact binomial test
 
data:  2 and 10
number of successes = 2, number of trials = 10, p-value = 0.05469
alternative hypothesis: true probability of success is less than 0.5
95 percent confidence interval:
0.0000000 0.5069013
sample estimates:
probability of success
0.2
 
> binom.test(2,10,.5,alternative="greater")
 
Exact binomial test
 
data:  2 and 10
number of successes = 2, number of trials = 10, p-value = 0.9893
alternative hypothesis: true probability of success is greater than 0.5
95 percent confidence interval:
0.03677144 1.00000000
sample estimates:
probability of success
0.2

Here, when the sample frequency is 20% and http://freakonometrics.blog.free.fr/public/perso2/test-lies-56.gif is equal to 10, we accept at the same time that theta is higher than 50% and lower than 50%.
And obviously it is not only a theoretical problem: it has obviously some strong implications. This morning, a good friend mentioned a post published some months ago, online here, about discrimination, and the lack of women with academic positions in mathematics, in France. As claimed by the author of the post“A Paris VI, meilleure université française selon son président, sur 11 postes de maitres de conférences, 5 filles classées premières. Il y a donc des filles excellentes ? A Toulouse, sur 4 postes, 2 filles premières. Parité parfaite. Mais à côté de cela, Bordeaux, 4 postes, 0 fille première. Littoral, 3 postes, 0 fille, Nice, 5 postes, 0 fille, Rennes, 7 postes, 0 fille…”.
Consider the latter one: in Rennes, out of 7 people hired last year, no woman. So in some sense, it looks obvious that there is some kind of discrimination ! Zero out of seven ! Well, if we consider the fact that around 30% of PhD thesis in mathematics were defended by women those years, we can also try to see is there if no “positive discrimination“, i.e. test http://freakonometrics.blog.free.fr/public/perso2/test-lies-60.gif where theta is the probability to hire a woman (just to be a little bit provocative).

> prop.test(0,7,.3,alternative="less")
 
1-sample proportions test with continuity correction
 
data:  0 out of 7, null probability 0.3
X-squared = 1.7415, df = 1, p-value = 0.09347
alternative hypothesis: true p is less than 0.3
95 percent confidence interval:
0.0000000 0.3719021
sample estimates:
p
0
 
Warning message:
In prop.test(0, 7, 0.3, alternative = "less") :
Chi-squared approximation may be incorrect
> binom.test(0,7,.3,alternative="less")
 
Exact binomial test
 
data:  0 and 7
number of successes = 0, number of trials = 7, p-value = 0.08235
alternative hypothesis: true probability of success is less than 0.3
95 percent confidence interval:
0.0000000 0.3481637
sample estimates:
probability of success
0

With no woman hired that year, we can still pretend that there was some kind of “positive discrimination“. An note that we do accept – with more confidence – the assumption of “positive discrimination” if we look at all universities together,

> prop.test(5+2,11+4+4+3+5+7,.3,alternative="less")
 
1-sample proportions test with continuity correction
 
data:  5 + 2 out of 11 + 4 + 4 + 3 + 5 + 7, null probability 0.3
X-squared = 1.021, df = 1, p-value = 0.1561
alternative hypothesis: true p is less than 0.3
95 percent confidence interval:
0.0000000 0.3556254
sample estimates:
p
0.2058824
 
> binom.test(5+2,11+4+4+3+5+7,.3,alternative="less")
 
Exact binomial test
 
data:  5 + 2 and 11 + 4 + 4 + 3 + 5 + 7
number of successes = 7, number of trials = 34, p-value = 0.1558
alternative hypothesis: true probability of success is less than 0.3
95 percent confidence interval:
0.0000000 0.3521612
sample estimates:
probability of success
0.2058824

So obviously, with small sample, almost anything can be claimed !

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.