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[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,
> 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…
> 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) [1] 0.917397
On obtient la statistique de Fisher, et le plus simple est alors de calculer la -value
> 1-pf(F,dl0-dl1,dl1) [1] 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)
OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (October 27, 2012). Une application du test de Fisher, partie 1. Freakonometrics. Retrieved December 5, 2024 from https://doi.org/10.58079/oun7