Pour la fin du cours de statistique, on joue enfin avec des données. Après les taux de réussite au brevet des collèges, on a travaillé un peu sur la survenance des ouragans aux Etats-Unis. La base utilisée est la suivante

> base=read.table(
"http://freakonometrics.free.fr/ouragan.csv",
sep=";",header=TRUE,dec=",")
> base=base[1:207,1:7]

Le nombre annuel d’ouragan est la série suivante

>  TB <- table(base$Year)
>  years <- as.numeric(names(TB))
>  counts <- as.numeric(TB)
>  years0=(1900:2005)[which(!(1900:2005)%in%years)]
>  db <- data.frame(years=c(years,years0),
         counts=c(counts,rep(0,length(years0))))
>  X=db$counts
>  plot(db,type="h")

On peut comparer les fréquences observées avec ce que donnerait une loi de Poisson (paramètre obtenu par maximum de vraisemblance)

>  lambda=mean(X)
>  n=length(X)
>  cbind(X,dpois(0:6,lambda)*n)
  [,1]      [,2]
0   16 15.038430
1   34 29.367500
2   27 28.674870
3   13 18.665717
4    5  9.112744
5    6  3.559128
6    5  1.158396

On peut alors faire un test d’ajustement,

>  library(vcd)
>  fit=goodfit(X,"poisson")
>  plot(fit)

>  summary(fit)
 
	 Goodness-of-fit test for poisson distribution
 
                      X^2 df   P(> X^2)
Likelihood Ratio 14.17686  5 0.01452415

On voit que le test d’ajustement d’une loi Poisson n’est pas très concluant ici. Notons que ni la loi binomiale

>  fit=goodfit(X,"binomial")
>  plot(fit)
>  summary(fit)
 
	 Goodness-of-fit test for binomial distribution
 
                      X^2 df     P(> X^2)
Likelihood Ratio 46.21196  5 8.222959e-09

ni la loi binomiale négative ne semble marcher.

>  fit=goodfit(X,"nbinomial")
>  plot(fit)
>  summary(fit)
 
	 Goodness-of-fit test for nbinomial distribution
 
                      X^2 df   P(> X^2)
Likelihood Ratio 10.91507  4 0.02753523

On peut faire un test de stabilité temporelle du paramètre d’une loi de Poisson, à l’aide du test de Przyborowski & Wilenski (1940).  Un ancien package permettait de mettre en place ce test. Mais on peut heureusement utiliser une autre fonction, avec une version différente du test (on pourra aussi consulter Krishnamoorthy & Thomson (2004) pour aller plus loin). L’idée du premier test est assez simple: notons que si sur la première période (de durée ), on a un processus de Poisson de paramètre  et sur la seconde période (de durée ), alors

et

Un test de  s’écrit  Notons que la loi
conditionnelle de  sachant  suit une loi binomiale, donc le paramètre de probabilité est fonction de ce ratio. Plus précisément, conditionnellement à,  suit un loi binomiale  avec

A partir de là, il est facile de mettre en place une stratégie de test. Une autre stratégie est d’utiliser

>  seuil=1960
>  X1=db[db$years<seuil,"counts"]
>  X2=db[db$years>=seuil,"counts"]
>  N1=length(db$years<seuil)
>  N2=length(db$years>=seuil)
>  poisson.test(c(sum(X1),sum(X2)),c(N1,N2))
 
	Comparison of Poisson rates
 
data:  c(sum(X1), sum(X2)) time base: c(N1, N2)
count1 = 85, expected count1 = 103.5, p-value = 0.01216
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.5218662 0.9265666
sample estimates:
rate ratio 
 0.6967213

 

Bref, on peut accepter que le nombre annuel d’ouragan n’est pas un processus de Poisson homogène dans le temps…