Tag Archives: test

Introduction aux tests nonparamétriques

Aujourd’hui, on va commencer la partie du cours sur les tests, avec une introduction par les tests non-paramétriques. On parlera (dans la continuité de ce qu’on a fait jusqu’à présent) de tests d’ajustement de lois, mais aussi de tests plus divers, comme des tests d’indépendance (entre observations successives).

Le test le plus classique est le test des signes. L’idée est de regarder puis de compter le nombre de fois où https://latex.codecogs.com/gif.latex?X_{i-1}\leq%20X_i (différence positive) et le nombre de fois où https://latex.codecogs.com/gif.latex?X_{i-1}\leq%20X_i (différence négative).

> (X=rnorm(20))
 [1]  0.81352495  0.78076709 -0.58654359  0.52805014 -0.31101282 -0.65022034 -0.93619729 -0.16483729 -1.97396753 -0.22506985  1.14959851 -0.38895448 -0.09298421
> (S=c("N","P")[1+1*(diff(X)>0)])
 [1] "N" "N" "P" "N" "N" "N" "P" "N" "P" "P" "N" "P" "P" "N" "P" "N" "P" "P" "N"

On verra ce test en cours. On peut par exemple compter le nombre de fois où la différence était positive,

> sim=function(a=0){
+ n=100
+ epsilon=X=rnorm(n)
+ for(t in 2:n) X[t]=a*X[t-1]+epsilon[t]
+ S=c("N","P")[1+1*(diff(X)>0)]
+ return(sum(S=="P"))}

(je fais ici une fonction un peu plus compliqué pour voir ce qui se passe si les données ne sont pas indépendantes, justement). Si on génère plein d’échantillons, on peut voir que le nombre de variations positives suit une loi normale, si les observations sont effectivement indépendantes,

> library(ggplot2)
> set.seed(1)
> vectS=unlist(lapply(rep(0,1000),sim))
> ggplot(data=data.frame(vectS = vectS), 
+ aes(x = vectS)) +
+ geom_histogram(aes(y = ..density..),
+ binwidth = 1, fill = "light blue", 
+ col = "black") +
+ stat_function(fun=function(x) dnorm(x,
+ (n-1)/2+.5,sqrt((n+1)/12)),
+ col = "red", size = 1) + ylab("Densite") +
+ ggtitle("Histogramme")

Continue reading Introduction aux tests nonparamétriques

Multiple Tests, an Introduction

Last week, a student asked me about multiple tests. More precisely, she ran an experience over – say – 20 weeks, with the same cohort of – say – 100 patients. An we observe some https://latex.codecogs.com/gif.latex?X_{i,t}

size=100
nb=20
set.seed(1)
X=matrix(rnorm(size*nb),size,nb)

(here, I just generate some fake data). I can visualize some trajectories, over the 20 weeks,

library(RColorBrewer)
cl1=brewer.pal(12,"Set3")[12]
cl2=brewer.pal(8,"Set2")[2:7]
cl=c(cl1,cl2)
boxplot(X)
for(i in 1:length(cl)){
lines(1:20,X[i,],type="b",col=cl[i])  }

http://freakonometrics.hypotheses.org/files/2014/09/Selection_585.png

Continue reading Multiple Tests, an Introduction

Temperatures Series as Random Walks

Last year, I did mention in a post that unit-root tests are dangerous, because they might lead us to strange models. For instance, in a post, I did obtain that the temperature observed in January 2013, in Montréal, might be considered as a random walk process (or at leat an integrated process). The code to extract the data has changed (since the website has been updated), so here, we use

library(RCurl)
library(XML)
options(RCurlOptions = list(useragent = "R"))
HEURE=0:23
extracttemp=function(Y,M,D){
url=paste(
"http://climate.weather.gc.ca/climateData/hourlydata_e.html?timeframe=1&Prov=QC&StationID=5415&Year=",Y,"&Month=",
M,"&Day=",D,sep="")
wp <- getURLContent(url)
doc <- htmlParse(wp, asText = TRUE) 
docName(doc) <- url
tmp <- readHTMLTable(doc)
basejour=data.frame(Year=Y,Month=M,Day=D,
Hour=HEURE,Temp=as.numeric(as.character(data.frame(tmp[2])[,2]))[2:25])
return(basejour)}
B=NULL
for(y in 1955:2013){
for(d in 1:31){
B=rbind(B,extracttemp(y,1,d))}}

Here are all the temperatures observed, and 2013,

plot(B$X,B$Temp,cex=.5,col="light blue",xlab="January, in Montreal",ylab="Temperature (Celsius)")
I=which(B$Year==2013)
lines(B$X[I],B$Temp[I],col="red")

In the previous post, one test only was used, and one year was considered. I was wondering if this behavior was observed only with temperature of 2013 (or not), and how the other tests (mentioned in a previous post too) were performing.

I might need a function, because those tests cannot be used if there is a missing value, even only one. So I did use the value observed one hour before (just to make sure that the tests can be done)

correcty=function(Y){
I=which(is.na(Y))
	if(length(I)==0){Yc=Y}
	if(length(I)>0){Yc=Y;for(i in I) Yc[i]=Yc[i-1]}	
return(Yc)
}

Now, we can compute the p-values, for all the years, and the three different three (keeping in mind that two test if the series is non-stationary, and one if the series is stationary)

DF=matrix(NA,2013-1954,3)
library(urca)
for(y in 1955:2013){
Z=B$Temp[which(B$Year==y)]
	Zc=correcty(Z)
	DF[y-1954,2]=as.numeric(pp.test(Zc)$p.value)	
	DF[y-1954,1]=as.numeric(kpss.test(Zc)$p.value)	
	DF[y-1954,3]=as.numeric(adf.test(Zc)$p.value)
}

Visually, if red means stationary, and blue means non-stationary, we get

DFP=DF
DFP[,1]=DF[,1]<.05
DFP[,2:3]=DF[,2:3]>.05
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")
plot(0:1,0:1,xlim=c(1950,2015),ylim=c(0,3),axes=FALSE,xlab="",ylab="")
axis(1)
text(1952,.5,"KPSS")
text(1952,1.5,"PP")
text(1952,2.5,"ADF")
for(y in 1955:2013){
for(i in 1:3){
polygon(y+c(-1,-1,1,1)/2.2,i-.5+c(-1,1,1,-1)/2.2,col=CL[1+(DFP[y-1954,i]==1)*5],border=NA)}}

Quite frequently, we conclude that the temperature is a random walk. Which does not make sense (from a physical point of view). But again, it might come from the fact that temperature are stationary, but with some fractional behavior (as suggested in the previous post).

Unit Root Tests

This week, in the MAT8181 Time Series course, we’ve discussed unit root tests. According to Wold’s theorem, if https://latex.codecogs.com/gif.latex?(Y_t) is  (weakly) stationnary then

https://latex.codecogs.com/gif.latex?%20%20%20%20%20Y_{t}=\sum%20_{{j=0}}^{\infty%20}\psi_{j}\varepsilon%20_{{t-j}}+\xi%20_{t}

where https://latex.codecogs.com/gif.latex?%20(\varepsilon%20_{{t}}) is the innovation process, and where https://latex.codecogs.com/gif.latex?%20(\xi%20_{{t}}) is some deterministic series (just to get a result as general as possible). Observe that

https://latex.codecogs.com/gif.latex?\sum%20_{{j=0}}^{{\infty%20}}|\psi_{{j}}|^{2}%20%3C%20\infty

as discussed in a previous post. To go one step further, there is also the Beveridge-Nelson decomposition : an integrated of order one process, defined as

https://latex.codecogs.com/gif.latex?%20%20%20%20%20\Delta%20Y_{t}=(1-L)%20Y_t=\sum%20_{{j=0}}^{\infty%20}\psi_{j}\varepsilon%20_{{t-j}}+\xi=\Psi(L)\varepsilon%20_{{t}}+\xican be represented as

a linear trend https://latex.codecogs.com/gif.latex?+ a random walk https://latex.codecogs.com/gif.latex?+ a stationary remaining term

i.e.

https://latex.codecogs.com/gif.latex?%20Y_{t}=\underbrace{Y_0%20+%20\xi%20t%20%C2%A0}+\underbrace{\Psi(1)\sum_{j=1}^t\varepsilon%20_{{i}}}+\underbrace{\tilde\Psi(L)\varepilon_0-\tilde\Psi(L)\varepsilon_t}

where https://latex.codecogs.com/gif.latex?%20\tilde\Psi(\cdot) is the polynomial with terms https://latex.codecogs.com/gif.latex?%20\tilde\psi_j, where

https://latex.codecogs.com/gif.latex?%20\tilde\psi_j%20=\sum_{i=j+1}^\infty\psi_i

For unit-root tests, we will use various representation of the process. In order to illustrate the implementation of those tests, consider the following series

> E=rnorm(240)
> X=cumsum(E)
> plot(X,type="l")
  • Dickey Fuller (standard)

Here, for the simple version of the Dickey-Fuller test, we assume that

https://latex.codecogs.com/gif.latex?%20Y_t=\alpha+\beta%20t+\varphi%20Y_{t-1}+\varepsilon_t

and we would like to test if https://latex.codecogs.com/gif.latex?%20\varphi=1 (or not). We can write the previous representation as

https://latex.codecogs.com/gif.latex?%20\Delta%20Y_t=\alpha+\beta%20t+[\varphi-1]%20Y_{t-1}+\varepsilon_t

so we simply have to test if the regression coefficient in the linear regression is – or not – null. Which can be done with Student’s test. If we consider the previous model without the linear drift, we have to consider the following regression

> lags=0
> z=diff(X)
> n=length(z)
> z.diff=embed(z, lags+1)[,1]
> z.lag.1=X[(lags+1):n]
> summary(lm(z.diff~0+z.lag.1 ))

Call:
lm(formula = z.diff ~ 0 + z.lag.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.84466 -0.55723 -0.00494  0.63816  2.54352 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.005609   0.007319  -0.766    0.444

Residual standard error: 0.963 on 238 degrees of freedom
Multiple R-squared:  0.002461,	Adjusted R-squared:  -0.00173 
F-statistic: 0.5873 on 1 and 238 DF,  p-value: 0.4442

Our testing procedure will be based on the Student’s t value,

> summary(lm(z.diff~0+z.lag.1 ))$coefficients[1,3]
[1] -0.7663308

which is exactly the value computed using

> library(urca)
> df=ur.df(X,type="none",lags=0)
> df

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: -0.7663

The interpretation of this value can be done using critical values (99%, 95%, 90%)

> qnorm(c(.01,.05,.1)/2)
[1] -2.575829 -1.959964 -1.644854

If the statistics exceeds those values, then the series is not stationnary, since we cannot reject the assumption that https://latex.codecogs.com/gif.latex?%20\varphi-1=0. So we might conclude that there is a unit root. Actually, those critical values are obtained using

> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 

Call:
lm(formula = z.diff ~ z.lag.1 - 1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.84466 -0.55723 -0.00494  0.63816  2.54352 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.005609   0.007319  -0.766    0.444

Residual standard error: 0.963 on 238 degrees of freedom
Multiple R-squared:  0.002461,	Adjusted R-squared:  -0.00173 
F-statistic: 0.5873 on 1 and 238 DF,  p-value: 0.4442

Value of test-statistic is: -0.7663 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

The problem with R is that there are several packages that can be used for unit root tests. Just to mention another one,

> library(tseries)
> adf.test(X,k=0)

	Augmented Dickey-Fuller Test

data:  X
Dickey-Fuller = -2.0433, Lag order = 0, p-value = 0.5576
alternative hypothesis: stationary

We do have here also a test where the null hypothesis is that there is a unit root. But the p-value is quite different. What is odd is that we have

> 1-adf.test(X,k=0)$p.value
[1] 0.4423705
> df@testreg$coefficients[4]
[1] 0.4442389

(but I think it is a coincidence).

  • Augmented Dickey Fuller

It is possible to had some lags in the regression. For instance, we can consider

https://latex.codecogs.com/gif.latex?%20\Delta%20Y_t=\alpha+\beta%20t+[\varphi-1]%20Y_{t-1}+\psi%20\Delta%20Y_{t-1}+\varepsilon_t

Again, we have to check if one coefficient is null, or not. And this can be done using Student’s t test.

> lags=1
> z=diff(X)
> n=length(z)
> z.diff=embed(z, lags+1)[,1]
> z.lag.1=X[(lags+1):n]
> k=lags+1
> z.diff.lag = embed(z, lags+1)[, 2:k]
> summary(lm(z.diff~0+z.lag.1+z.diff.lag ))

Call:
lm(formula = z.diff ~ 0 + z.lag.1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87492 -0.53977 -0.00688  0.64481  2.47556 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657

Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared:  0.003292,	Adjusted R-squared:  -0.005155 
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777

> summary(lm(z.diff~0+z.lag.1+z.diff.lag ))$coefficients[1,3]
[1] -0.7328138

This value is the one obtained using

> df=ur.df(X,type="none",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 

Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87492 -0.53977 -0.00688  0.64481  2.47556 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657

Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared:  0.003292,	Adjusted R-squared:  -0.005155 
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777

Value of test-statistic is: -0.7328 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

And again, other pckages can be used:

> adf.test(X,k=1)

	Augmented Dickey-Fuller Test

data:  X
Dickey-Fuller = -1.9828, Lag order = 1, p-value = 0.5831
alternative hypothesis: stationary

Hopefully, the conclusion is the same (we should reject the assumption that the series is stationary, but I am not sure about the computation of the p-value).

  • Augmented Dickey Fuller with trend and drift

So far, we have not included the drift in our model. But this is simple to do (this will be called the augmented version of the previous procedure): we just have to include a constant in the regression,

> summary(lm(z.diff~1+z.lag.1+z.diff.lag ))

Call:
lm(formula = z.diff ~ 1 + z.lag.1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.91930 -0.56731 -0.00548  0.62932  2.45178 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared:  0.02313,	Adjusted R-squared:  0.01482 
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393

The statistics of interest are obtained here considering some analysis of variance outputs, where this model is compared with the one without the integrated part, and the drift,

> summary(lm(z.diff~1+z.lag.1+z.diff.lag ))$coefficients[2,3]
[1] -2.303948
> anova(lm(z.diff ~ z.lag.1 + 1 + z.diff.lag),lm(z.diff ~ 0 + z.diff.lag))$F[2]
[1] 2.732912

Those two values are the ones obtained also with

> df=ur.df(X,type="drift",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.91930 -0.56731 -0.00548  0.62932  2.45178 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared:  0.02313,	Adjusted R-squared:  0.01482 
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393

Value of test-statistic is: -2.3039 2.7329 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

And we can also include a linear trend,

> temps=(lags+1):n
> summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))

Call:
lm(formula = z.diff ~ 1 + temps + z.lag.1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87727 -0.58802 -0.00175  0.60359  2.47789 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
temps       -0.0004194  0.0009767  -0.429   0.6680  
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared:  0.0239,	Adjusted R-squared:  0.01139 
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287

> summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))$coefficients[3,3]
[1] -1.98282
> anova(lm(z.diff ~ z.lag.1 + 1 + temps+ z.diff.lag),lm(z.diff ~ 1+ z.diff.lag))$F[2]
[1] 2.737086

while R function returns

> df=ur.df(X,type="trend",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.87727 -0.58802 -0.00175  0.60359  2.47789 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
tt          -0.0004194  0.0009767  -0.429   0.6680  
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared:  0.0239,	Adjusted R-squared:  0.01139 
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287

Value of test-statistic is: -1.9828 1.8771 2.7371 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47
  • KPSS test

Here, in the KPSS testing procedure, two models can be considerd : with a drift, or with a linear trend. Here, the null hypothesis is that the series is stationnary.

With a drift, the code is

> summary(ur.kpss(X,type="mu"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 4 lags. 

Value of test-statistic is: 0.972 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.73

while it will be, in the case there is a trend

> summary(ur.kpss(X,type="tau"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: tau with 4 lags. 

Value of test-statistic is: 0.5057 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

One more time, it is possible to use another package to get the same test (but again, a different output)

> kpss.test(X,"Level")

	KPSS Test for Level Stationarity

data:  X
KPSS Level = 1.1997, Truncation lag parameter = 3, p-value = 0.01

> kpss.test(X,"Trend")

	KPSS Test for Trend Stationarity

data:  X
KPSS Trend = 0.6234, Truncation lag parameter = 3, p-value = 0.01

At least, there is some kind of consistency, since we keep rejecting the stationnary assumption, for that series.

  • Philipps-Perron test

The Philipps-Perron test is based on the ADF procedure. The code is here

> PP.test(X)

	Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller = -2.0116, Truncation lag parameter = 4, p-value = 0.571

with again, a possible alternative with the other package

> pp.test(X)

	Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller Z(alpha) = -7.7345, Truncation lag parameter = 4, p-value
= 0.6757
alternative hypothesis: stationary
  •  Comparison

I will not spend more time comparing the different codes, in R, to run those tests. Let us spend some additional time on a quick comparison of those three procedure. Let us generate some autoregressive processes, with more or less autocorrelation, as well as some random walk, and let us see how those tests perform :

> n=100
> AR=seq(1,.7,by=-.01)
> P=matrix(NA,3,31)
> M1=matrix(NA,1000,length(AR))
> M2=matrix(NA,1000,length(AR))
> M3=matrix(NA,1000,length(AR))

> for(i in 1:(length(AR)+1)){
+ for(s in 1:1000){
+ if(i==1) X=cumsum(rnorm(n))
+ if(i!=1) X=arima.sim(n=n,list(ar=AR[i]))
+ library(urca)
+ M2[s,i]=as.numeric(pp.test(X)$p.value)
+ M1[s,i]=as.numeric(kpss.test(X)$p.value)
+ M3[s,i]=as.numeric(adf.test(X)$p.value)
+ }}

Here, we would like to count how many times the p-value of our tests exceed 5%,

> prop05=function(x) mean(x>.05)
+ P[1,]=1-apply(M1,2,prop05)
+ P[2,]=apply(M2,2,prop05)
+ P[3,]=apply(M3,2,prop05)
+ }
> plot(AR,P[1,],type="l",col="red",ylim=c(0,1),ylab="proportion of non-stationnary 
+ series",xlab="autocorrelation coefficient")
> lines(AR,P[2,],type="l",col="blue")
> lines(AR,P[3,],type="l",col="green")
> legend(.7,1,c("ADF","KPSS","PP"),col=c("green","red","blue"),lty=1,lwd=1)

 

We can see here how poorly Dickey-Fuller test behave, since a 50% (at least) of our autoregressive processes are considered as non-stationnary.

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.

Playing cards in Vegas?

In a previous post, a few weeks ago, I mentioned that I will be in Las Vegas by the end of July. And I took the opportunity to write a post on roulette(s). Since some colleagues told me I should take some time to play poker there, I guess I have to understand how to play poker… so I went back to basics on cards, and shuffling techniques.

Now, I have to confess that I have been surprised, while I was looking for mathematical models for shuffling, to find so many deterministic techniques (and results related to algebra, and cycles).

On http://mathworld.wolfram.com/ for instance, one can find nice articles on so-called in-shuffle or out-shuffle techniques. There is also a great article, Golomb (1961), but mainly on algebraic properties of permutations by cutting and shuffling, as well as Diaconis, Kantor and Graham’s The Mathematics of Perfect Shuffle And if you look at Monge’s shuffle, you can find a deterministic recursive relationship. As a statistician (or applied probabilist), I should confess that I did not find answer to the question I wanted to ask : how long should we shuffle before getting cards randomly sorted in ours hands ?

  • Randomness (from a statistician perspective)

First, I need to define (as properly as possible) a notion of “cards randomly sorted“. Consider a game with 32 cards. Why 32 ? Mathematicians will tell you that 32 is a great number, since it is a power of 2, so there might be interesting (algebraic) properties when shuffling. From a computational point of view, 32 is smaller than 52, so my random generations will run faster. This is basically why I used 32. 10 would have been better, but not realistic with cards.

So, our 32 cards can be seen as a vector, or a list, of 32 items, say

In order to assess if my cards are randomly sorted, let us get back to number properties (real valued numbers). If there were 10 cards, the list can be seen as an element of the following set

(or to be more specific, a subset of that set, since numbers have to be different – it has to be a permutation – we cannot have duplicates, we’ll get back to that point in a few seconds). Let us see this list as a decimal number, with 10 digits. More precisely,

Now,  it is natural to say that cards are randomly sorted is this number is uniformly distributed on the unit interval, isn’t it ? (if we use the same shuffle many times, with the same starting point)

Well, if we think about it twice, uniform on the unit interval is probably not the proper distribution, since (as mentioned above) all digits have to be different. For instance, the smallest number would be  and the largest  . But as we will see, it this uniform assumption might not be too strong, actually.

And if we want to get back to our initial problem, with 32 cards, we simply have to use a decomposition in the 32-basis.

So if we have an algorithm to shuffle cards, we just have to run it several times (with the same starting value) and see when  starts to be uniformly distributed. We start with a Dirac distribution, we have some kind of transition matrix, we expect our limiting distribution to be uniform and we wonder when the limiting distribution is reached… And from a statistical point of view, that should not be that difficult to assess, since we do have several goodness of fit tests that can be used.

Actually, it is possible to check if our technique passes the test of a uniform distribution, when digit are randomly generated (without replacement). The code to generate  is

> j = 32
> X3 = (0:(j-1))[sample(1:j)] 
> x3 = sum(j^(-(1:j))*X3)

If we run it a few times, and check if the assumption of a uniform distribution is valid (on samples with, say, 500 observations),

> P3=NULL
> for(i in 1:10000){
+   U3=NULL
+   for(s in 1:500){
+     X3 =(0:(j-1))[sample(1:j)] 
+     x3 =sum(j^(-(1:j))*X3)
+     U3 =c(U3,x3)}
+   P3 =c(P3,ks.test(U3,punif)$p.value)
+ }

in 95% of the scenarios, the -value exceeds 5%

> mean(P3>.05)
[1] 0.9529

(which is something we should have under the null), More precisely, we can check that the -value is uniformly distributed on the unit interval.

> hist(P3,freq=FALSE)

So assuming that our number is uniform on the unit interval might be a good notion for “cards are randomly sorted“.

What we need now is some shuffling algorithms. Or to be more specific, some feasible shuffling algorithm. I mean here that I just start playing with cards, so it should be some techniques that I should be able to perform, to understand how it works…. So you will have to wait a few weeks before I start talking about the riffle or dovetail shuffle (you know the kind of shuffle in which half of the deck is held in each hand, and then cards are released by the thumbs so that they fall to the table interleaved… like in the movies) !

  • Top in at random shuffle, and related (simple) algorithm

My first algorithm is simple: the top-in at random shuffle. We start with the following ordering

    N=1:m

There are https://latex.codecogs.com/gif.latex%20?m cards, and n denote the place where the card on top will go.

    n=sample(2:m,size=1)
    if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])  
    if(n==m) N=c(N[2:n],N[1])

Then, we repeat that transfer of the card on top several times.

schuffle1=function(m,ns=10){
  N=1:m
  for(i in 1:ns)
    {
    n=sample(2:m,size=1)
    if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])  
    if(n==m) N=c(N[2:n],N[1])
    }
return(N)}

Now, it is also possible to consider a bottom-in at random shuffle. The idea is the same, the only difference it that you start from the card at the bottom of the deck. But that would be the same as the one before (in terms of time before reaching randomness)

    n=sample(1:(m-1),size=1)
    if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])  
    if(n==1) N=c(N[m],N[1:(m-1)])

Why not mixing ? Randomly. Call it randomly mixed top-bottom in at random shuffle. You start either with the card on top, or at bottom (with identical probability), of the deck and then move the card somewhere,

     card=sample(c("top","bottom"),size=1)
     if(card=="top"){
       n=sample(2:m,size=1)
       if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])  
       if(n==m) N=c(N[2:n],N[1])}
     if(card=="bottom"){
       n=sample(1:(m-1),size=1)
       if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])  
       if(n==1) N=c(N[m],N[1:(m-1)])}

All those codes can be together (within the same function),

schuffle1=function(m,ns=10,which="top"){
  N=1:m
if(which=="top"){
  for(i in 1:ns)
    {
    n=sample(2:m,size=1)
    if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])  
    if(n==m) N=c(N[2:n],N[1])
    }}
if(which=="bottom"){
  for(i in 1:ns)
    {
    n=sample(1:(m-1),size=1)
    if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])  
    if(n==1) N=c(N[m],N[1:(m-1)])
    }}
  if(which=="mixed"){
    for(i in 1:ns)
    {card=sample(c("top","bottom"),size=1)
     if(card=="top"){
       n=sample(2:m,size=1)
       if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])  
       if(n==m) N=c(N[2:n],N[1])
       }
     if(card=="bottom"){
       n=sample(1:(m-1),size=1)
       if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])  
       if(n==1) N=c(N[m],N[1:(m-1)])
       }
    }}  
  return(N)}

But why do we take only one card ? It won’t be more complex to take 2. Or 3. Or more.

  • Tops in at random shuffle, and related (mixed) algorithm

Yes, I used tops to say that we would take several cards on top of the deck. Say a random number of cards. And then, the strategy is the same, so the previous code is (slightly) adapted, as follows

     k=sample(1:(m-1),size=1)
     n=sample((k+1):m,size=1); if(k==m-1) n=m
     if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])  
     if(n==m) N=c(N[(k+1):n],N[1:k])

The idea is the following, here

As earlier, it is possible to take cards at the bottom of the deck, or, one more time, to use a mixed strategy. The codes would be

     card=sample(c("top","bottom"),size=1)
     if(card=="top"){
		 k=sample(1:(m-1),size=1)
		 n=sample((k+1):m,size=1); if(k==m-1) n=m
		 if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])  
		 if(n==m) N=c(N[(k+1):n],N[1:k])}
     if(card=="bottom"){
		 k=sample(2:m,size=1)
		 n=sample(1:(k-1),size=1); if(k==1) n=1
		 if(n>1)  N=c(N[1:(n-1)],N[k:m],N[n:(k-1)])  
		 if(n==1) N=c(N[k:m],N[n:(k-1)])}

Again, it is possible to have all those codes in the same function,

schuffle2=function(m,ns=10,which="top"){
  N=1:m
  if(which=="top"){
    for(i in 1:ns)
    {
      k=sample(1:(m-1),size=1)
      n=sample((k+1):m,size=1); if(k==m-1) n=m
      if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])  
      if(n==m) N=c(N[(k+1):n],N[1:k])
    }}
  if(which=="bottom"){
    for(i in 1:ns)
    {
      k=sample(2:m,size=1)
      n=sample(1:(k-1),size=1); if(k==1) n=1
      if(n>1)  N=c(N[1:(n-1)],N[k:m],N[n:(k-1)])  
      if(n==1) N=c(N[k:m],N[n:(k-1)])
    }}
  if(which=="mixed"){
    for(i in 1:ns)
    {card=sample(c("top","bottom"),size=1)
     if(card=="top"){
		 k=sample(1:(m-1),size=1)
		 n=sample((k+1):m,size=1); if(k==m-1) n=m
		 if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])  
		 if(n==m) N=c(N[(k+1):n],N[1:k])
     }
     if(card=="bottom"){
		 k=sample(2:m,size=1)
		 n=sample(1:(k-1),size=1); if(k==1) n=1
		 if(n>1)  N=c(N[1:(n-1)],N[k:m],N[n:(k-1)])  
		 if(n==1) N=c(N[k:m],N[n:(k-1)])
     }
    }}  
  return(N)}
  • How long should we shuffle before having cards randomly sorted ?

With the codes mentioned above, it is possible to run generations of shuffles,

distu=function(k=100,j=32){
	U1B=U1T=U1M=U2B=U2T=U2M=U3=NULL
	for(s in 1:100){
		X1T=(0:(j-1))[schuffle1(j,k,"top")] 
		X1B=(0:(j-1))[schuffle1(j,k,"bottom")] 
		X1M=(0:(j-1))[schuffle1(j,k,"mixed")] 
		X2T=(0:(j-1))[schuffle2(j,k,"top")] 
		X2B=(0:(j-1))[schuffle2(j,k,"bottom")] 
		X2M=(0:(j-1))[schuffle2(j,k,"mixed")]
		X3 =(0:(j-1))[sample(1:j)] 

		x1T=sum(j^(-(1:j))*X1T)
		x1B=sum(j^(-(1:j))*X1B)
		x1M=sum(j^(-(1:j))*X1M)
		x2T=sum(j^(-(1:j))*X2T)
		x2B=sum(j^(-(1:j))*X2B)
		x2M=sum(j^(-(1:j))*X2M)
		x3 =sum(j^(-(1:j))*X3)

		U1T=c(U1T,x1T)
		U1B=c(U1B,x1B)
		U1M=c(U1M,x1M)
		U2T=c(U2T,x2T)
		U2B=c(U2B,x2B)
		U2M=c(U2M,x2M)
		U3 =c(U3,x3)
    }
	B=list(U1T=U1T,...)
}

and then, we run tests to see if the samples can be assumed to be uniformly distributed on the unit interval, e.g. for the very first kind first shuffle describe above, it would be

ks.test(B$U1T,punif)$p.value

More precisely, we use the following function, to estimate to proportion of scenarios where the -value exceeds 5%,

PV=function(k){
	P1B=P1T=P1M=P2B=P2T=P2M=P3=NULL
	for(i in 1:10000){
        B=dist(k,j=32)
		P1T=c(P1T,ks.test(B$U1T,punif)$p.value)
		P1M=c(P1M,ks.test(B$U1M,punif)$p.value)
		P2T=c(P2T,ks.test(B$U2T,punif)$p.value)
		P2M=c(P2M,ks.test(B$U2M,punif)$p.value)
		P3 =c(P3,ks.test(B$U3,punif)$p.value)}
	return(list(
		p1T=mean(P1T>.05),
		p1M=mean(P1M>.05),
		p2T=mean(P2T>.05),
		p2M=mean(P2M>.05),				
		p3=mean(P3>.05)))}

If we plot the results, we have

K=1:100
MP=Vectorize(PV)(K)
plot(K,MP[1,],col="red",type="b",ylim=0:1,pch=1)
lines(K,MP[2,],type="b",pch=19,col="red")
lines(K,MP[3,],col="blue",type="b",pch=1)
lines(K,MP[4,],type="b",pch=19,col="blue")
lines(K,MP[5,],type="b",pch=3,col="black")

Here, we look at the proportion of -values that exceed 5%. We can pretend that we have a uniform distribution if that proportion is close to 95%. So basically, we just have to see when we reached for the first time the 95% region.If we zoom in the upper part of the graph, we get

With 32 cards,

  • with a top in at random, we have to shuffle about 70 or 80 cards before having a randomly sorted set of cards. Which is large, but which is quite intuitive. One can imagine that it might take a while before getting the cards at bottom much higher in the pack,
  • with a randomly mixed top in at random strategy, it is faster, slightly (we do not have that problem with cards at bottom that stay at bottom), since it takes about 60 or 70 cards.
  • with a tops in a random, it goes again faster, with about 35 rounds,
  • with a randomly mixed tops-bottoms in at random, it takes about 10 to 15 rounds.

Those results were obtained on tests on samples of size 100. The same code ran on a server during the week-end, with samples of size 500. Note that the output is rather close,

Note that those algorithm were mentioned because they were feasible, not only from a computational point of view, but when playing with real cards, in paper. Like with kids. I can actually ask my kids to perform those shuffle techniques next time we play with cards. The good thing is that randomly mixed tops-bottoms in at random shuffle technique: kids can do it 10 times, and cards should be randomly ordered in the deck…

Now, for those willing to see more algorithms, there are the so-called Fisher-Yates (also Knuth) shuffle. But may I keep that for another post ?

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.

Tests on tail index for extremes

Since several students got the intuition that natural catastrophes might be non-insurable (underlying distributions with infinite mean), I will post some comments on testing procedures for extreme value models.

A natural idea is to use a likelihood ratio test (for composite hypotheses). Let http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif denote the parameter (of our parametric model, e.g. the tail index), and we would like to know whether http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif is smaller or larger than http://freakonometrics.blog.free.fr/public/perso5/lrtest22.gif (where in the context of finite versus infinite mean http://freakonometrics.blog.free.fr/public/perso5/lrtest23.gif). I.e. either http://freakonometrics.blog.free.fr/public/perso5/lrtest21.gif belongs to the set http://freakonometrics.blog.free.fr/public/perso5/lrtest-10.gif or to its complementary http://freakonometrics.blog.free.fr/public/perso5/lrtest-11.gif. Consider the maximum likelihood estimator http://freakonometrics.blog.free.fr/public/perso5/lrtest24.gif, i.e.

http://freakonometrics.blog.free.fr/public/perso5/lrtest-9.gif

Let http://freakonometrics.blog.free.fr/public/perso5/lrtest25.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-3.gif denote the constrained maximum likelihood estimators on http://freakonometrics.blog.free.fr/public/perso5/lrtest26.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest27.gif respectively,

http://freakonometrics.blog.free.fr/public/perso5/lrtest-12.gif

http://freakonometrics.blog.free.fr/public/perso5/lrtest-2.gif

Either http://freakonometrics.blog.free.fr/public/perso5/lrtest-13.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-6.gif (on the left), or http://freakonometrics.blog.free.fr/public/perso5/lrtest-14.gif and http://freakonometrics.blog.free.fr/public/perso5/lrtest-7.gif (on the right)

So likelihood ratios

http://freakonometrics.blog.free.fr/public/perso5/lrtest-15.gif      http://freakonometrics.blog.free.fr/public/perso5/lrtest-16.gif

 are either equal to

http://freakonometrics.blog.free.fr/public/perso5/lrtest-19.gif      http://freakonometrics.blog.free.fr/public/perso5/lrtest-18.gif

or

http://freakonometrics.blog.free.fr/public/perso5/lrtest-20.gif        http://freakonometrics.blog.free.fr/public/perso5/lrtest-17.gif

If we use the code mentioned in the post on profile likelihood, it is easy to derive that ratio. The following graph is the evolution of that ratio, based on a GPD assumption, for different thresholds,

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> library(evir)
> X=base1$Loss.in.DKM
> U=seq(2,10,by=.2)
> LR=P=ES=SES=rep(NA,length(U))
> for(j in 1:length(U)){
+ u=U[j]
+ Y=X[X>u]-u
+ loglikelihood=function(xi,beta){
+ sum(log(dgpd(Y,xi,mu=0,beta))) }
+ XIV=(1:300)/100;L=rep(NA,300)
+ for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value }
+ plot(XIV,L,type="l")
+ PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)}
+ (L0=(OPT=optimize(f=PL,interval=c(0,10)))$objective)
+ profilelikelihood=function(beta){
+ -loglikelihood(1,beta) }
+ (L1=optim(par=1,fn=profilelikelihood)$value)
+ LR[j]=L1-L0
+ P[j]=1-pchisq(L1-L0,df=1)
+ G=gpd(X,u)
+ ES[j]=G$par.ests[1]
+ SES[j]=G$par.ses[1]
+ }
>
> plot(U,LR,type="b",ylim=range(c(0,LR)))
> abline(h=qchisq(.95,1),lty=2)

with on top the values of the ratio (the dotted line is the quantile of a chi-square distribution with one degree of freedom) and below the associated p-value

> plot(U,P,type="b",ylim=range(c(0,P)))
> abline(h=.05,lty=2)

In order to compare, it is also possible to look at confidence interval for the tail index of the GPD fit,

> plot(U,ES,type="b",ylim=c(0,1))
> lines(U,ES+1.96*SES,type="h",col="red")
> abline(h=1,lty=2)

To go further, see Falk (1995), Dietrich, de Haan & Hüsler (2002), Hüsler & Li (2006) with the following table, or Neves & Fraga Alves (2008). See also here or there (for the latex based version) for an old paper I wrote on that topic.

Régression sur des variables catégorielles

Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression

> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, 
+ family=poisson(link="log"))
> summary(reg0)

Call:
glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), 
    data = base, offset = log(exposition))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5717  -0.3968  -0.2996  -0.1547  12.6722  

Coefficients:
      Estimate Std. Error z value Pr(>|z|)    
zoneB -2.54187    0.06287  -40.43   <2e-16 ***
zoneA -2.54912    0.05285  -48.23   <2e-16 ***
zoneC -2.38525    0.03753  -63.56   <2e-16 ***
zoneD -2.13454    0.03878  -55.05   <2e-16 ***
zoneE -2.00204    0.03965  -50.49   <2e-16 ***
zoneF -2.06932    0.11547  -17.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 50966  on 50000  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(reg0,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
        1         2         3         4         5 
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,

> base$zone=relevel(base$zone,"B")
> regB=glm(nbre~zone,offset=log(exposition),data=base,
+ family=poisson(link="log"))
> summary(regB)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54187    0.06287 -40.431  < 2e-16 ***
zoneA       -0.00725    0.08213  -0.088 0.929661
zoneC        0.15662    0.07322   2.139 0.032432 *
zoneD        0.40733    0.07387   5.514 3.50e-08 ***
zoneE        0.53983    0.07433   7.263 3.80e-13 ***
zoneF        0.47255    0.13148   3.594 0.000325 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(regB,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
1         2         3         4         5
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,

> base$zone=relevel(base$zone,"A")
> reg=glm(nbre~zone,offset=log(exposition),
> data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54912    0.05285 -48.232  < 2e-16 ***
zoneB        0.00725    0.08213   0.088 0.929661
zoneC        0.16387    0.06482   2.528 0.011471 *
zoneD        0.41458    0.06555   6.324 2.54e-10 ***
zoneE        0.54708    0.06607   8.280  < 2e-16 ***
zoneF        0.47980    0.12699   3.778 0.000158 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.

> base$zonesimple=base$zone
> base$zonesimple[base$zone%in%c("A","B")]="A"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54612    0.04046 -62.937  < 2e-16 ***
zonesimpleC  0.16087    0.05518   2.915  0.00355 **
zonesimpleD  0.41158    0.05604   7.345 2.06e-13 ***
zonesimpleE  0.54408    0.05665   9.605  < 2e-16 ***
zonesimpleF  0.47681    0.12235   3.897 9.74e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence

> base$zonesimple=relevel(base$zonesimple,"C")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.38525    0.03753 -63.557  < 2e-16 ***
zonesimpleA -0.16087    0.05518  -2.915  0.00355 **
zonesimpleD  0.25071    0.05396   4.646 3.39e-06 ***
zonesimpleE  0.38321    0.05460   7.019 2.24e-12 ***
zonesimpleF  0.31593    0.12142   2.602  0.00927 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.06932    0.11547 -17.921  < 2e-16 ***
zonesimpleC -0.31593    0.12142  -2.602  0.00927 **
zonesimpleA -0.47681    0.12235  -3.897 9.74e-05 ***
zonesimpleD -0.06522    0.12181  -0.535  0.59232
zonesimpleE  0.06727    0.12209   0.551  0.58161
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles

> base$zonesimple[base$zone%in%c("D","E","F")]="F"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5660  -0.3959  -0.3004  -0.1547  12.5929

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.07182    0.02696 -76.853  < 2e-16 ***
zonesimpleC -0.31344    0.04621  -6.783 1.18e-11 ***
zonesimpleA -0.47431    0.04861  -9.757  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15698  on 49997  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> library(car)
> linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0"))
Linear hypothesis test

Hypothesis:
zonesimpleD = 0
zonesimpleE = 0

Model 1: restricted model
Model 2: nbre ~ zonesimple

Res.Df Df  Chisq Pr(>Chisq)
1  49997
2  49995  2 5.7073    0.05763 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression

> library(tree)
> arbre=tree(nbre~zone+offset(log(exposition)),
+ data=base,split="gini")
> plot(arbre)
> text(arbre)

De la créativité des gangsters

Pendant mon séjour récent en Nouvelle Angleterre, j’ai survolé le livre de Leonard Mlodinow, the drunkard’s walk. Et au hasard de mes lectures, je suis tombé sur la petite histoire suivante

Autrement dit, en utilisant les cinq derniers chiffres d’une quantité économique comme le us treasury balance, on aurait un générateur de nombre aléatoire…
Par contre la suite est un peu plus surprenante,


http://freakonometrics.hypotheses.org/files/2015/12/sammy.gifLa loi de Benford (que j’avais pu évoquer ici ou ) parle des premiers chiffres, mais cette fois on parle deslast five digits. Donc visiblement l’évoquation n’est pas pertinente ici. Mais qui sait ? Ca reste malgré tout une histoire intéressante. Considérons – histoire de tester cette légende – les deux sources de données suivantes,http://treasurydirect.gov/ et http://economagic.com/. Sur ce dernier site, un petit travail de mise en forme des données est nécessaire.

b1=read.table("http://freakonometrics.free.fr/debtus1.txt",
   header=TRUE,sep="\t")
b2=read.table("http://freakonometrics.free.fr/debtus2.txt",
   header=TRUE)
X1=as.character(b1$Dollar.Amount)
n1=nchar(X1)
Y1=substr(X1,n1-8,n1-3)
X1=as.numeric(substr(Y1,1,2))*1000+as.numeric(substr(Y1,4,6))
x=X1/100000
X2=b2$DEBT
Y2=trunc(as.numeric(X2))
X2=as.character(Y2)
n2=nchar(X2)
Y2=substr(X2,n2-4,n2)
y=as.numeric(Y2)/100000
y=y[y<1]

Pour rappel, un générateur aléatoire (standard) vérifie deux propriétés importantes

  • les nombres doivent être tirés suivant une loi uniforme sur [0,1], i.e. ici, si on divise les nombres à 5 chiffres par 10000,
  • les tirages doivent être indépendants entre eux.

La première propriété semble assez naturelle, et correspond à l’histoire racontée dans un commentaire posté ici (expliquant comment un casino avait été au bord de la faillite car une roulette faisait sortir certains chiffres trop souvent, et j’essayais de comprendre comment utiliser l’information qu’un chiffre sort plus souvent). La seconde est probablement encore plus importante.

  • Visualisation de la distribution

La première idée est de visualiser la densité de nos séries de chiffres. Pour éviter les problèmes de bord (et comme c’est juste en introduction) on va utiliser un histogramme, et pas une estimation à noyau.

hist(x,col="red")
hist(y,col="blue")

On obtient pour la première série la courbe rouge, et pour la seconde la courbebleue,

On note qu’a priori, pour la première série, l’hypothèse d’uniformité n’est peut être pas la plus réaliste…

  • Test de Kolmogorov-Smirnov

On peut aussi mettre en œuvre un test de Kolmogorov-Smirnov afin de tester si la loi uniforme est adaptée.:

> ks.test(x,"punif")
 
	One-sample Kolmogorov-Smirnov test
 
data:  x 
D = 0.1047, p-value = 0.01645
alternative hypothesis: two-sided 
 
> ks.test(y,"punif")
 
	One-sample Kolmogorov-Smirnov test
 
data:  y 
D = 0.0456, p-value = 0.3581
alternative hypothesis: two-sided

On retrouve ici confirmée l’intuition précédante: la loi uniforme est pertinente pour la seconde série, pas la première.

  • Les autocorrélations de la série

Travaillons uniquement sur la seconde série. On peut étudier l’autocorrélation de notre série de nombres, ou peut-être un peu plus malin, sur les quantiles gaussiens associés (les autocorrélations étant intéressantes pour les séries gaussiennes),:

plot(acf(y))
plot(acf(qnorm(y)))

ie. pour la série brute

et pour la série normalisée,

Bref, on pourrait être tenté de valider l’hypothèse d’indépendance entre les tirages.
  • Run test (de Wald–Wolfowitz)

L’idée est de comparer une série de chiffres à la médiane, s’ils sont plus grands, on note + (ou A) et sinon – (ou B). On crée alors une série du genre “+++−−++−−++++++−−−” et on compte les séries de + et les séries de -, les runs,

library(lawstat)
runs.test(y,plot=TRUE)

Runs Test - Two sided

data:  y
Standardized Runs Statistic = -0.2462, p-value = 0.8055

Bref, la légende me semble à prendre avec des pincettes (car fonction de la source considérée), même si l’idée est intéressante (si l’on met de côté les aspects d’aléa moral). Et l’analyse sur la loi de Benford ne semble pas valide: les derniers chiffres sur les grands nombres ne se comportent pas du tout comme les premiers.

STT2700, estimation, tests et coupes du monde

Mercredi, pour le dernier cours, nous allons revenir sur l’estimation, les tests, et plus généralement sur la modélisation statistique. Pour cela, j’avais pensé travailler sur les nombres de buts marqués, par match, lors de différentes coupes du monde de soccer (1982, 1998 et 2010). Je ne mets pas l’intégralité du code aujourd’hui, l’idée est pour l’instant de mettre en ligne des données qui serviront à répondre aux questions qui seront posées mercredi. Le code (accompagné – éventuellement – d’explications théoriques) sera posté par la suite.

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)

Les boxplot associés à ces trois échantillons sont les suivants,

On va se poser des questions autour de ces données, par exemple voir s’il est vraisemblance (ou pas) que le nombre moyen de but dans un match (avant prolongation, s’il y en a eu). On peut commencé par essayer de se demander quel modèle utiliser. Classiquement, la loi de Poisson est la plus utilisée (en plus, c’est la seule loi qui est autorisée lorsqu’on publie un billet le 1er avril). Les histogrammes sont les suivants

boxplot(S82,S98,S10,col=c("red","yellow","blue"),
label=("1982","1998","2010"))
hist(S82,breaks=0:11,col="red")
hist(S98,breaks=0:11,col="yellow")
hist(S10,breaks=0:11,col="blue")

Si on compare les fonctions de répartition empiriques à celles de lois de Poisson ajustées par maximum de vraisemblance, on obtient, pour la coupe du monde de 1982

et pour celle de 2010,

Visuellement, l’ajustement semble relativement bon, surtout en 2010. On peut aussi faire un test du chi-deux,

> library(vcd)
> (GF=goodfit(S10,type="poisson"))

Observed and fitted values for poisson distribution
with parameters estimated by ML

 count observed     fitted
     0        7  6.6409703
     1       17 15.0459484
     2       13 17.0442384
     3       14 12.8719508
     4        7  7.2907534
     5        5  3.3036226
     6        0  1.2474617
     7        1  0.4037543

> summary(GF)

	 Goodness-of-fit test for poisson distribution

                      X^2 df  P(> X^2)
Likelihood Ratio 5.586765  5 0.3485255

On voit que l’on accepte l’ajustement par une loi de Poisson. Pour ceux qui veulent une visualisation, sur la figure ci-dessous, on a la densité d’une loi du chi-deux. Le premier trait vertical est la valeur observée, et l’aire jaune est alors la p-value (qui excède largement 5%). En rouge on a 5%, donc le second trait vertical est la borne de la région critique associé au test pour une erreur de première espèce valant 5%,

On peut ensuite faire plein de tests.  On suppose que . On va pouvoir tester

http://freakonometrics.free.fr/test-soccer-04.gif

contre une hypothèse alternative

http://freakonometrics.free.fr/test-<br /><br /> soccer-06.gif

Comme on a une hypothèse sur la loi des observations qui semble robuste, on peut utiliser un test de type rapport de vraisemblance.
On peut aussi faire un test de la forme

http://freakonometrics.free.fr/test-soccer-09.gif

contre

http://freakonometrics.free.fr/test-soccer-10.gif

(histoire de tester des hypothèses simples – qui ont une interprétation). Sinon, comme ce qui nous intéresse, c’est de savoir si on a plus de trois buts par matchs, on peut définir la variable binomiale http://freakonometrics.free.fr/test-soccer-03.gif, en notant que

http://freakonometrics.free.fr/test-soccer-02.gif

est une proportion – donc facile à tester – qui nous intéresse ici compte tenu du problème que l’on cherchera à résoudre. On pourra alors tester, par exemple

http://freakonometrics.free.fr/test-soccer-08.gif

contre

http://freakonometrics.free.fr/test-soccer-07.gif

Ces derniers tests sont alors facile à mettre en œuvre,

> Z=(S10>=3)*1
> prop.test(sum(Z),length(Z),p=1/2,alternative="less")

	1-sample proportions test with continuity correction

data:  sum(Z) out of length(Z), null probability 1/2 
X-squared = 1.2656, df = 1, p-value = 0.1303
alternative hypothesis: true p is less than 0.5 
95 percent confidence interval:
 0.0000000 0.5322764 
sample estimates:
       p 
0.421875

On peut aussi faire des tests de moyenne sur

http://freakonometrics.free.fr/test-soccer-11.gif

Un test de l’hypothèse

http://freakonometrics.free.fr/test-soccer-04.gif

contre une hypothèse alternative

http://freakonometrics.free.fr/test-soccer-06.gif

s’écrit alors
> t.test(S10,mu=3,alternative ="less")

	One Sample t-test

data:  S10 
t = -3.7763, df = 63, p-value = 0.0001775
alternative hypothesis: true mean is less than 3 
95 percent confidence interval:
     -Inf 2.590273 
sample estimates:
mean of x 
 2.265625
Mais on l’aura l’occasion de revoir tous les points du cours, y compris aller peut être un peu plus loin, par exemple sur la comparaison de moyenne entre échantillons,
> t.test(S82,S98,var.equal=FALSE)

	Welch Two Sample t-test

data:  S82 and S98 
t = 0.427, df = 85.266, p-value = 0.6704
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.5503669  0.8514658 
sample estimates:
mean of x mean of y 
 2.807692  2.657143

Test d’ajustement et lâcher de bombes

Avant les applications demain en cours, un petit billet (presque d’actualité) sur une application évoquée vendredi dernier sur le test du chi-deux: l’ajustement de lois.

Le problème est le suivant: un pays se fait bombarder. Et les dirigeants doivent se demander si certaines cibles sont visées (auquel cas il peut être légitime de les déplacer) ou au contraire si les tirs sont aléatoires. Quand je dis que des cibles sont visées, j’entends aussi par là que les pilotes savent viser. Car les pilotes ont (je pense, ou j’espère) un carnet de route à suivre, avec des cibles à viser

Ce que l’on va se demander c’est plutôt si les pilotes savent – ou arrive à – viser. Bref, le problème peut se modéliser en supposant que le lancer de bombes se fait (globalement) selon un processus de Poisson. Localement, si les tirs sont aléatoires, on devrait observer des tirages de lois de Poisson. C’est tout du moins la théorie que R. D. Clarke (alors actuaire chez Prudential Assurance Company) a utilisé pendant la seconde guerre mondiale, sur les bombardements à Londres (en ligne ici). Cette histoire (vraie) est reprise par Thomas Pynchon dans Gravity’s Rainbow (en ligne ici)

(l’exemple est même repris dans Feller). Bref, le problème n’est pas de savoir quel serait le paramètre de la loi de Poisson, mais si la loi de Poisson est adaptée, ou pas. C’est ce qu’on appelle un problème d’ajustement de loi.

  • Test du chi-deux et ajustement de lois

Jusqu’à présent, on avait supposé que les observations suivaient une certaine loi, e.g. une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/chi2-16.gif, et on cherchait à tester une hypothèse de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-13.gif

versus

http://freakonometrics.hypotheses.org/files/2015/12/chi2-14.gif

Ici on va cherche à utiliser un test sur une loi multinomiale, de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-01.gif

versus

http://freakonometrics.hypotheses.org/files/2015/12/chi2-02.gif

L’hypothèse nulle est ici une égalité vectorielle,

http://freakonometrics.hypotheses.org/files/2015/12/chi2-03.gif

ou encore

http://freakonometrics.hypotheses.org/files/2015/12/chi2-56.gif

Dans le cas d’un test d’ajustement de lois de Poisson, si on suppose que les observations suivent une loi http://freakonometrics.hypotheses.org/files/2015/12/chi2-09.gif, on va utiliser un test sur une loi multinomiale, avec

http://freakonometrics.hypotheses.org/files/2015/12/chi2-08.gif

Le soucis est que, puisque l’on se limite à un vecteur de taille finie, le vecteur ne sera pas dans le simplexe (cf ici). Donc classiquement, pour la dernière valeur, on retient une hypothèse de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-10.gif

Le test est alors basé sur la statistique de Pearson

http://freakonometrics.hypotheses.org/files/2015/12/chi2-11.gif

qui va suivre (si effectivement les observations suivent une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/chi2-09.gif) une loi http://freakonometrics.hypotheses.org/files/2015/12/chi2-12.gif.
Rappelons que cette statistique peut aussi s’écrire

http://freakonometrics.hypotheses.org/files/2015/12/chi2-21.gif

 

  • Application pour étudier la précision d’un lancer de bombes

Pendant la second guerre mondiale, R.D. Clarke étudia les impacts de bombes V1 et V2 tombées dans une zone de 144 km2 dans le sud de Londres (l’article original, publié après guerre en 1946, est en ligne ici). Il divisa cette zone en 576 zones de 0,25 km2 et compta le nombre d’impact dans chacune des zones. Il obtint plus précisément la distribution suivante

nbre impacts par zone 0 1 2 3 4 5 et plus
fréquence (nbre zones) 229 211 93 35 7 1

(en fait, on sait que le “5
et plus” correspond à 7, car sait qu’il y a eu 537 bombes sur 576 zones). Avant de se lancer tête baissée, réfléchissons un peu au type de loi que l’on pourrait utiliser. Pour cela, notons http://freakonometrics.hypotheses.org/files/2015/12/Nb.gif le nombre de points tombés dans un ensemble http://freakonometrics.hypotheses.org/files/2015/12/asubb.gif (ou http://freakonometrics.hypotheses.org/files/2015/12/calA.gif désigne la région globale). Si on suppose qu’un nombre aléatoire http://freakonometrics.hypotheses.org/files/2015/12/Npois.gif de points sont lancés aléatoirement dans http://freakonometrics.hypotheses.org/files/2015/12/calA.gif, et que http://freakonometrics.hypotheses.org/files/2015/12/Npois.gif suit une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/lambda.gif, alors http://freakonometrics.hypotheses.org/files/2015/12/Nb.gif suit une loi de Poisson de paramètre

http://freakonometrics.hypotheses.org/files/2015/12/ratio-aire.gif

Bref, si on regarde plusieurs plusieurs régions http://freakonometrics.hypotheses.org/files/2015/12/calB.gif (de même taille, éventuellement – afin de garder toutes les observations – formant une partition de http://freakonometrics.hypotheses.org/files/2015/12/calA.gif) et si on observe une loi de Poisson, c’est que dans la région http://freakonometrics.hypotheses.org/files/2015/12/calA.gif, les tirs sont faits au hasard.

> (n=c(229,211,93,35,7,0,0,1))
[1] 229 211  93  35   7   0   0   1
> y=0:7
> (lambda=sum(n*y)/sum(n))
[1] 0.9322917
> prob=dpois(y,lambda)
> freq.theo=sum(n)*prob
> freq.emp =n
> cbind(y,trunc(freq.theo),freq.emp)
     y     freq.emp
[1,] 0 226      229
[2,] 1 211      211
[3,] 2  98       93
[4,] 3  30       35
[5,] 4   7        7
[6,] 5   1        0
[7,] 6   0        0
[8,] 7   0        1

On peut commencer par regarder la log-vraisemblance

> logvrais=function(L){sum(log(dpois(y,L))*n)}
> param=seq(.5,1.5,by=.025)
> LV=sapply(param,logvrais)
> plot(param,LV,type="b",col="blue")
http://freakonometrics.hypotheses.org/files/2015/12/logvrais1-bombes.png

La statistique du chi-deux, elle, ressemble à ça

> chi2=function(L)sum(n)*sum(((n/sum(n)-dpois(y,L))^2)/dpois(y,L))
> C2=sapply(param,chi2)
> plot(param,C2,type="b",col="red")xxx.
http://freakonometrics.hypotheses.org/files/2015/12/chi2-bombes.png

Bon, le soucis c’est qu’en tronquant le vecteur (on suppose que le nombre maximum d’impact est 7), la somme des probabilités ne fait pas exactement un. On peut regrouper dans une même classe les fréquences élevées (comme le fait Clarke dans le papier initial d’ailleurs).

> (n=c(229,211,93,35,7,1))
[1] 229 211  93  35   7   1
> y=0:4
> prob=c(dpois(y,lambda),1-ppois(4,lambda))
> freq.theo=sum(n)*prob
> freq.emp =n
> (Q=sum(1
[1] 1.169155
> 1-pchisq(Q,length(y)-1)
[1] 0.8831505

On retrouve les quantités évoquées dans l’article de Clarke. La  statistique de test vaut 1.16 et la p-value associée est de l’ordre de 88%. On va donc accepter l’hypothèse de loi de Poisson,

> chi2=function(L){
+ sum(n)*sum(2)^2)/
+ c(dpois(y,L),1-ppois(4,L)) )}
> C2=sapply(param,chi2)
> plot(param,C2,type="b",col="red")

La valeur de la statistique du chi-deux en fonction du paramètre de la loi de Poisson est représentée ci-dessous,

http://freakonometrics.hypotheses.org/files/2015/12/chi2-bombes-v2.png

et le trait horizontal est la valeur seuil de la région critique (en dessous, on accepte l’ajustement d’une loi de Poisson),

> abline(h=qchisq(.95,length(y)-1),lty=2)

Il existe plusieurs fonctions sous R permettant de faire des choses semblables,

> library(vcd)
> (n=c(229,211,93,35,7,0,0,1))
[1] 229 211  93  35   7   0   0   1
> nsim=c(rep(y[0],n[0]),rep(y[1],n[1]),
       rep(y[2],n[2]),rep(y[3],n[3]),
       rep(y[4],n[4]),rep(y[5],n[5]),
       rep(y[6],n[6]),rep(y[7],n[7]),
       rep(y[8],n[8]))
> gf=goodfit(nsim,type="poisson",method="ML")
> summary(gf)
 
 Goodness-of-fit test for poisson distribution
 
 X^2 df P(> X^2)
Likelihood Ratio 4.007784 3 0.2606249
 
> gf=goodfit(nsim,type="poisson",method="MinChisq")
> summary(gf)
 
 Goodness-of-fit test for poisson distribution
 
 X^2 df P(> X^2)
Pearson 1.275499 3 0.7349592

Bref, quelle que soit la méthode utilisée, on notera que l’on accepte toujours l’hypothèse d’une loi de Poisson. Autrement dit les bombes étaient envoyés un peu au hasard dans Londres….

Et si on y réfléchit un peu, d’un point de vue de la théorie des jeux, c’est effectivement une stratégie optimale…

  1. freq.theo-freq.emp)^2)/freq.theo []
  2. n/sum(n)-c(dpois(y,L),1-ppois(4,L []

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 !