# Somewhere else, part 18

A nice post (found a couple of days ago) on an insteresting topic, I should perhaps discuss more often on this blog

Also, this week, a series of posts and articles about L’Aquila trial,

Still, a lot of very interesting posts found here, and there,   Also several documentaries, found, online

Et comme toujours, quelques articles en francais, Did I miss something ?

# Une application du test de Fisher, partie 1

Un petit (?) billet pour revenir sur un point évoqué mardi dernier en cours, sur l’utilisation du test de Fisher (et des techniques d’analyse de variance). Je vais juste évoquer rapidement un vieux billet sur l’histoire du 5% utilisé lorsque l’on calcule la -value (ce qui complètera la réponse faite en classe mardi lorsque la question a été posée). Mardi dernier, j’avais évoqué rapidement le test d’une contrainte linéaire sur les paramètres. On se demande ici si les deux modèles suivants sont vraiment différents, entre et sa version contrainte On cherche alors à tester l’hypothèse contre l’hypothèse alternative Pour ça, on peut utiliser ce qu’a proposé Ronald Fisher, à savoir d’utiliser la somme des carrés des résidus (qui correspondent – à un facteur multiplicatif près – à l’estimateur de la variance du bruit), et il a posé Je parlais tout à l’heure de choix de modèle… on notera que l’on en est pas très loin car on peut réécrire cette expression en faisant apparaître le  des régressions, contraintes et noncontraintes. Fisher a montré que sous l’hypothèse H0, i.e. la loi de la statistique de Fisher est précisément une loi de Fisher, Pour aller plus loin sur l’analyse de la variance (ANOVA, analysis of variance), je renvois àFaraway (2005).

J’avais dit oralement qu’une application intéressant du test de nullité de plusieurs coefficient dans une régression était son utilisation pour des variables catégorielles, et pour faire du regroupement de modalité.

Considérons un cas d’école. On dispose des résultats des votes dans plusieurs villes, en France, lors des dernières élections présidentielles (ces bases sont en ligne sur http://www.data.gouv.fr/), opposant François Hollande (avec un H) et Nicolas Sarkozy (avec un S),

```> election=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/
205cf26f15a43974fdb7d3b25397613d.csv",
+ sep=";",header=TRUE,dec=",")
> head(election)
Code.du.dpt Inscrits Votants VoixH X.VoixExpH VoixS X.VoixExpS
1           1      592     506   195      41.05   280      58.95
2           1      215     175    82      52.23    75      47.77
3           1     8208    6447  3009      50.41  2960      49.59
4           1     1152     949   340      38.16   551      61.84
5           1      105      83    35      44.30    44      55.70
6           1     1702    1461   630      46.46   726      53.54```

On va nettoyer un peu la base

```> sbelection=election[election\$Code.du.département
+ %in%as.character(1:95),]```

car il y a des départements d’outre mers et la Corse, qui a une codification particulière. Bref, on retient les 95 départements métropolitain, moins le 20ème qui est la Corse. On va commencer par créer une variable de département qui soit un facteur qualitatif. Pour l’instant, le département est un numéro, entre 1 et 95 mais ce n’est aucunement un nombre !

```> sbelection\$Dpt=(paste("D",
> sbelection\$Code.du.département,sep=""))
> I=as.numeric(as.character(sbelection\$Code.du.département))<10
> sbelection\$Dpt[I]=(paste("D0",
+ sbelection\$Code.du.département,sep=""))
> sbelection\$Dpt=factor(sbelection\$Dpt)```

On peut alors faire une régression, où on peut essayer de modéliser le taux de vote pour François Hollande (par exemple) en fonction de variables explicatives, et en particulier, on veut voir s’il n’y a pas un effet spatial. Bref, on va régresser sur le département. Et d’ailleurs que sur le département, inutile de mettre d’autre variable pour illustrer.

```> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)```

Par contre, je fais une régression pondérée, car les villes (et villages) ne font pas du tout la même taille. C’est exactement la même idée que de faire une moyenne pondérée pour avoir le taux moyen dans une région (on rediscutera des régressions pondérées en fin de cours).

`> COEF=coefficients(summary(reg))`

La sortie ressemble à

```> COEF
Estimate Std. Error   t value Pr(>|t|)
DptD01 42.77423  0.4397101  97.27826        0
DptD02 52.39423  0.4572016 114.59765        0
DptD03 56.89434  0.5468760 104.03518        0
DptD04 51.06544  0.7807748  65.40355        0
DptD05 50.70525  0.8403979  60.33481        0
DptD06 35.69903  0.3246129 109.97414        0
DptD07 53.44185  0.5564489  96.04089        0
DptD08 51.86319  0.6360901  81.53434        0
DptD09 64.67228  0.8013389  80.70528        0
DptD10 42.61811  0.6182921  68.92877        0```

[…]

```DptD89 46.87270  0.5642231  83.07476        0
DptD90 50.47561  0.9063735  55.68964        0
DptD91 53.41748  0.3172310 168.38670        0
DptD92 49.54305  0.2850597 173.79888        0
DptD93 65.29062  0.3338954 195.54217        0
DptD94 56.48855  0.3196643 176.71211        0
DptD95 53.89382  0.3344373 161.14778        0```

On peut d’ailleurs faire une prédiction, par département (en prenant comme prédiction celle obtenue pour la première ville du département)

```> chg=which(diff(as.numeric(as.character(
+ sbelection\$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)```

On peut ensuite essayer de faire un petite dessin pour visualiser cette prédiction,

```> chg=which(diff(as.numeric(as.character(
+ sbelection\$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)
> nomdpt=read.table(
+ "http://freakonometrics.blog.free.fr/
public/data/departements.txt",sep=";")
> nomdpt\$nom=tolower(nomdpt\$V2)
> dpt=as.character(c(paste("0",
+ 1:9,sep=""),c(10:19,21:95)))
> noms=as.character(nomdpt[nomdpt\$V3%in%dpt,2])
> couleurliste=rev(rainbow(n=100))
> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur=couleurliste[trunc((PredDept-35)*3)]```

C’est un peu long (et compliqué) car il faut identifier le bon département sur la carte R (qui est codée en lettres). D’ailleurs la fonction suivante ne marche pas

```> match=match.map(francemap,dpt)
> francemap\$names[which(is.na(match.map(
+ francemap,dpt))==FALSE)]
 "Indre-et-Loire" "Indre"```

car en l’occurrence, seule l’Indre aurait du être obtenue. On va faire la correspondance à la main,

```> match=rep(NA,length(francemap\$names))
> for(i in 1:length(francemap\$names)){
+ if(sum(dpt==francemap\$names[i])>0)
+ {match[i]=which(dpt==francemap\$names[i])}
+ }```

Une fois les départements identifiés, on colorie,

```> color=couleur[match]
> map(database="france", fill=TRUE, col=color)```

C’est joli, mais comme toujours c’est un peu compliqué de faire des cartes… A titre d’information, les couleurs correspondent aux valeurs suivantes (avec en abscisse, le taux de vote obtenu par François Hollande) Bref, on a beaucoup de modalité car il y a 95 département. Mais on devrait pouvoir regrouper. Si on ne sait pas trop par où commencer, prenons un département au hasard, disons le 35 (i.e. Rennes),

`> k=34`

On peut ensuite regarder les estimations des coefficients, ou mieux, les intervalles de confiances des différents estimateurs, et voir ceux qui chevauche celui de Rennes (on va prendre des intervalles de confiance à 70% histoire d’avoir des chances que les deux coefficients soient vraiment proches)

```> alpha=.7
> CI.COEF=confint(reg,level=alpha)```

On peut visualiser ces intervalles, avec une bande horizontale pour Rennes,

```> plot(c(1:19,21:95),COEF[,1],ylim=c(35,
+ 60),pch=3,cex=.5)
> segments(c(1:19,21:95),CI.COEF[,1],
+ c(1:19,21:95),CI.COEF[,2])
> lk=which((CI.COEF[,1]<CI.COEF[k,2])&
+ (CI.COEF[,2]>CI.COEF[k,1]))
> segments(c(1:19,21:95)[lk],CI.COEF[lk,1],
+ c(1:19,21:95)[lk],CI.COEF[lk,2],col="red",lwd=2)```

et en rouge, les intervalles qui la chevauchent Ça y est, on a une liste potentielle de six départements à regrouper… On va prendre comme modalité de référence le département qui nous intéresse, et voir si les coefficients des autres départements sont (ou pas) significatifs.

```> sbelection\$Dpt=relevel(sbelection\$Dpt,"D35")
> reg35=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)
> COEF35=coefficients(summary(reg35))
> COEF35[c(1,lk[lk<k]+1,lk[lk>k]),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 55.66226653  0.3275407 169.9399801 0.0000000
DptD11       0.59043789  0.6272106   0.9413710 0.3465210
DptD32       0.96845567  0.7855154   1.2328920 0.2176241
DptD36      -0.03007518  0.7414655  -0.0405618 0.9676455
DptD62       0.56238714  0.4247897   1.3239189 0.1855384
DptD75      -0.06226653  0.4098013  -0.1519432 0.8792326
DptD81      -0.11970183  0.6067040  -0.1972986 0.8435950```

On a confirmation que ces six départements sont non significativement (individuellement) différents de celui de référence (et donc on accepte l’hypothèse de nullité du coefficient). Mais peut-on pour autant les regrouper? C’est là que l’on peut faire un test de Fisher. Pour cela on fait deux régressions. La première est celle avec un regroupement de départements en une modalité unique (les sept sont désormais ensemble)

```> sbelection\$DptGroupe=sbelection\$Dpt
> sbelection\$DptGroupe[sbelection\$Dpt%
+ in%substr(rownames(COEF[lk,]),4,6)]=paste("D",k+1,sep="")
> model0=lm(X.VoixExpH~DptGroupe,
+ weights=Votants,data=sbelection)```

et la seconde est celle que l’ont vient de faire (qui sera notre hypothèse alternative, i.e. il ne faut pas regrouper les départements ensemble car ils sont différents – ou au moins l’un est différent)

```> model1=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)```

On peut alors faire le test de Fisher (en continuant à corriger de l’effet taille de la population, et en faisant des sommes pondérées)

```> SCR0=sum(sbelection\$Votants*residuals(model0)^2)
> SCR1=sum(sbelection\$Votants*residuals(model1)^2)
> dl0=model0\$df.residual
> dl1=model1\$df.residual
> (F=(SCR0-SCR1)/(dl0-dl1)/SCR1*dl1)
 0.917397```

On obtient la statistique de Fisher, et le plus simple est alors de calculer la -value

```> 1-pf(F,dl0-dl1,dl1)
 0.4809411```

On a fait (sans le savoir) une analyse de la variance, qui peut se faire plus simplement (?) en utilisant la fonction suivante

```> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36117 2280538612
2  36111 2280191044  6    347568 0.9174 0.4809```

Et on accepte ici l’hypothèse que ces sept départements peuvent être regroupé en un seul.

Enfin, si statistiquement, il y a du sens à regrouper les modalités, on peut se demander si ce regroupement de régions géographique à du sens,

```> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur="red"
> match=match.map(francemap,dpt)
> color=couleur[match]
> map(dataLa première  estbase="france", fill=TRUE, col=color)``` Les régions que l’on souhaite associer à Rennes sont fort éloignées, et les regrouper n’a peut-être pas de sens… (à suivre donc).