{"word_count": 1178, "character_count": 10713}
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, entreet sa version contrainteOn cherche alors à tester l'hypothèsecontre l'hypothèse alternativePour ç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),&#62; election=read.table(+ "http://freakonometrics.blog.free.fr/public/data/205cf26f15a43974fdb7d3b25397613d.csv",+ sep=";",header=TRUE,dec=",")&#62; head(election)Code.du.dpt Inscrits Votants VoixH X.VoixExpH VoixS X.VoixExpS1           1      592     506   195      41.05   280      58.952           1      215     175    82      52.23    75      47.773           1     8208    6447  3009      50.41  2960      49.594           1     1152     949   340      38.16   551      61.845           1      105      83    35      44.30    44      55.706           1     1702    1461   630      46.46   726      53.54On va nettoyer un peu la base&#62; 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 !&#62; sbelection$Dpt=(paste("D",&#62; sbelection$Code.du.département,sep=""))&#62; I=as.numeric(as.character(sbelection$Code.du.département))&#60;10&#62; sbelection$Dpt[I]=(paste("D0",+ sbelection$Code.du.département,sep=""))&#62; 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.&#62; 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).&#62; COEF=coefficients(summary(reg))La sortie ressemble à&#62; COEFEstimate Std. Error   t value Pr(&#62;|t|)DptD01 42.77423  0.4397101  97.27826        0DptD02 52.39423  0.4572016 114.59765        0DptD03 56.89434  0.5468760 104.03518        0DptD04 51.06544  0.7807748  65.40355        0DptD05 50.70525  0.8403979  60.33481        0DptD06 35.69903  0.3246129 109.97414        0DptD07 53.44185  0.5564489  96.04089        0DptD08 51.86319  0.6360901  81.53434        0DptD09 64.67228  0.8013389  80.70528        0DptD10 42.61811  0.6182921  68.92877        0[...]DptD89 46.87270  0.5642231  83.07476        0DptD90 50.47561  0.9063735  55.68964        0DptD91 53.41748  0.3172310 168.38670        0DptD92 49.54305  0.2850597 173.79888        0DptD93 65.29062  0.3338954 195.54217        0DptD94 56.48855  0.3196643 176.71211        0DptD95 53.89382  0.3344373 161.14778        0On 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)&#62; chg=which(diff(as.numeric(as.character(+ sbelection$Code.du.département)))!=0)&#62; newelection=sbelection[1+c(0,chg),]&#62; PredDept=predict(reg,newdata=newelection)On peut ensuite essayer de faire un petite dessin pour visualiser cette prédiction,&#62; chg=which(diff(as.numeric(as.character(+ sbelection$Code.du.département)))!=0)&#62; newelection=sbelection[1+c(0,chg),]&#62; PredDept=predict(reg,newdata=newelection)&#62; nomdpt=read.table(+ "http://freakonometrics.blog.free.fr/public/data/departements.txt",sep=";")&#62; nomdpt$nom=tolower(nomdpt$V2)&#62; dpt=as.character(c(paste("0",+ 1:9,sep=""),c(10:19,21:95)))&#62; noms=as.character(nomdpt[nomdpt$V3%in%dpt,2])&#62; couleurliste=rev(rainbow(n=100))&#62; library(maps)&#62; francemap&#60;-map(database="france")&#62; dpt=noms&#62; 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&#62; match=match.map(francemap,dpt)&#62; francemap$names[which(is.na(match.map(+ francemap,dpt[35]))==FALSE)][1] "Indre-et-Loire" "Indre"car en l'occurrence, seule l'Indre aurait du être obtenue. On va faire la correspondance à la main,&#62; match=rep(NA,length(francemap$names))&#62; for(i in 1:length(francemap$names)){+ if(sum(dpt==francemap$names[i])&#62;0)+ {match[i]=which(dpt==francemap$names[i])}+ }Une fois les départements identifiés, on colorie,&#62; color=couleur[match]&#62; 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),&#62; k=34On 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)&#62; alpha=.7&#62; CI.COEF=confint(reg,level=alpha)On peut visualiser ces intervalles, avec une bande horizontale pour Rennes,&#62; plot(c(1:19,21:95),COEF[,1],ylim=c(35,+ 60),pch=3,cex=.5)&#62; segments(c(1:19,21:95),CI.COEF[,1],+ c(1:19,21:95),CI.COEF[,2])&#62; lk=which((CI.COEF[,1]&#60;CI.COEF[k,2])&#38;+ (CI.COEF[,2]&#62;CI.COEF[k,1]))&#62; 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.&#62; sbelection$Dpt=relevel(sbelection$Dpt,"D35")&#62; reg35=lm(X.VoixExpH~Dpt,+ weights=Votants,data=sbelection)&#62; COEF35=coefficients(summary(reg35))&#62; COEF35[c(1,lk[lk&#60;k]+1,lk[lk&#62;k]),]Estimate Std. Error     t value  Pr(&#62;|t|)(Intercept) 55.66226653  0.3275407 169.9399801 0.0000000DptD11       0.59043789  0.6272106   0.9413710 0.3465210DptD32       0.96845567  0.7855154   1.2328920 0.2176241DptD36      -0.03007518  0.7414655  -0.0405618 0.9676455DptD62       0.56238714  0.4247897   1.3239189 0.1855384DptD75      -0.06226653  0.4098013  -0.1519432 0.8792326DptD81      -0.11970183  0.6067040  -0.1972986 0.8435950On 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)&#62; sbelection$DptGroupe=sbelection$Dpt&#62; sbelection$DptGroupe[sbelection$Dpt%+ in%substr(rownames(COEF[lk,]),4,6)]=paste("D",k+1,sep="")&#62; 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)&#62; 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)&#62; SCR0=sum(sbelection$Votants*residuals(model0)^2)&#62; SCR1=sum(sbelection$Votants*residuals(model1)^2)&#62; dl0=model0$df.residual&#62; dl1=model1$df.residual&#62; (F=(SCR0-SCR1)/(dl0-dl1)/SCR1*dl1)[1] 0.917397On obtient la statistique de Fisher, et le plus simple est alors de calculer la -value&#62; 1-pf(F,dl0-dl1,dl1)[1] 0.4809411On a fait (sans le savoir) une analyse de la variance, qui peut se faire plus simplement (?) en utilisant la fonction suivante&#62; anova(model0,model1)Analysis of Variance TableModel 1: X.VoixExpH ~ DptGroupeModel 2: X.VoixExpH ~ DptRes.Df        RSS Df Sum of Sq      F Pr(&#62;F)1  36117 22805386122  36111 2280191044  6    347568 0.9174 0.4809Et 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,&#62; library(maps)&#62; francemap&#60;-map(database="france")&#62; dpt=noms&#62; couleur="red"&#62; match=match.map(francemap,dpt)&#62; color=couleur[match]&#62; 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).