Tag Archives: élection

Report des voix entre les deux tours aux élections présidentielles en France

Allez, on va profiter des élections présidentielles, à deux tours, en France, pour jouer un peu, maintenant que les données commencent à être disponibles. On peut avoir les données par département (une centaine) ou par commune (36000)

# département
loc_dpt_1 = "https://www.data.gouv.fr/fr/datasets/r/312b99d0-27a7-429e-be8a-b930975545e9"
loc_dpt_2 = "https://www.data.gouv.fr/fr/datasets/r/c7b72e78-5e78-4976-995f-30533f0a35d3"
# commune
loc_com_1 = "https://www.data.gouv.fr/fr/datasets/r/54782507-e795-4f9d-aa70-ed06feba22e3"
loc_com_2 = "https://www.data.gouv.fr/fr/datasets/r/1601b027-48d1-42b1-96f1-b0fb49f55aa3"
download.file(loc_dpt_1,"dpt1.csv")
download.file(loc_dpt_2,"dpt2.csv")
d1 = read.csv("dpt1.csv")
d2 = read.csv("dpt2.csv")
download.file(loc_com_1,"com1.csv")
download.file(loc_com_2,"com2.csv")
c1 = read.csv("com1.csv")
c2 = read.csv("com2.csv")

Commençons par les départements. On peut récupérer les pourcentages de voix obtenus par les différents candidats, au premier tour, et regarder les pourcentages obtenus par Marine Le Pen. Par exemple, on peut regarder, département par département, le score de Marine Le Pen au premier et au second tour, pour vérifier que le report des voix s’opère bien

D1 = d1[d1$cand_nom=="LE PEN",c("dep_code","cand_rapport_inscrits","cand_rapport_exprim")]
D2 = d2[,c("CodeDépartement","LE.PEN.ins","LE.PEN.exp")]
names(D1) = c("dep_code","ins1","exp1")
names(D2) = c("dep_code","ins2","exp2")
D = merge(D1,D2)
par(mfrow=c(1,2))
plot(D$ins1,D$ins2,xlab="inscrits (%), premier tour",ylab="inscrits (%), second tour",col=colr2[1])
plot(D$exp1,D$exp2,xlab="exprimés (%), premier tour",ylab="exprimés (%), second tour",col=colr2[2],main="Le Pen / Le Pen",cex.main=1.1)

Comme on pouvait s’y attendre, on a une belle droite de régression, avec les nombres de voix obtenus, relativement aux inscrits, en bleu, et aux exprimés, à droite. Pour Emmanuel Macron au premier et Marine Le Pen au second tour, on a une belle corrélation négative, mais moins prononcée qu’auparavant. Autrement dit, un département qui votait peu Macron au premier tour a beaucoup voté Le Pen au second

Mais là, on parlait des deux candidats qui se sont maintenus. Si on regarde l’autre candidat de l’extrême droite, les reports sont assez étranges… après, les taux étaient “relativement” faibles au premier tour pour Erik Zemmour.

Pour le candidat en tête à gauche, on a une corrélation négative, mais sur les votes exprimées, on observera des outliers en haut à droite, avec manifestement un vote d’opposition au président sortant dans certains départements

Pour le candidat écologiste, la corrélation est elle encore plus clair

autrement dit, quand on a voté écologiste au premier tour, on avait peu tendance à voter Le Pen (ou disons les départements où on votait fortement écologiste n’ont pas voté Le Pen au second #ecologicalfallacy, etc).

Mais on peut faire le même exercice par commune. Mon seul soucis a été que les codes communes dans les deux bases étaient différents, donc j’ai utilisé les noms des communes, ce que je n’aime pas faire (avec les accents, les tirets, etc, la fusion peut être assez mauvaises… mais là ça marche bien). Je me suis ensuite limité aux grosses communes, de plus de 10000 inscrits,

C1 = c1[c1$cand_nom=="LE PEN",c("commune_name","cand_rapport_inscrits","cand_rapport_exprim","votants_nb")]
C2 = c2[,c("Commune","LE.PEN.ins","LE.PEN.exp")]
names(C1) = c("com_code","ins1","exp1","taille")
names(C2) = c("com_code","ins2","exp2")
D = merge(C1,C2)
D = D[D$taille>5000,]
par(mfrow=c(1,2))
plot(D$ins1,D$ins2,xlab="inscrits (%), premier tour",ylab="inscrits (%), second tour",col=colr2[1])
plot(D$exp1,D$exp2,xlab="exprimés (%), premier tour",ylab="exprimés (%), second tour",col=colr2[2],main="Le Pen / Le Pen",cex.main=1.1)

On peut alors faire le même exercice, avec Le Pen au premier et second tour

ou Macron au premier, Le Pen au second

On observe des beaux nuages, sur la gauche, respectivement en haut et en bas… Beaucoup de communes qui avaient peu voté pour Le Pen semblent avoir eu peur, et on massivement voté au second pour elle. Alors que ce n’est pas trop le cas pour celles où Le Pen avait obtenu beaucoup de voix. Et pour Macron, c’est en bas à gauche, autrement dit, beaucoup de communes qui n’ont pas voté Macron au premier tour n’ont pas pour autant voté Le Pen au second. Pour ceux que ça amuse, les codes sont disponibles pour regarder qui sont ces communes dans les coins.

Pour Zemmour, c’est assez confus, j’avoue, le personnage est déroutant (ou dégoutant ? comme le suggérait le correcteur automatique)

ou alors Jean Luc Mélanchon

ou la relativement franche corrélation des écologistes, qui semblent résolument aller en sens inverse avec le vote d’extrême droite…

 

Les sondages et le bruit

Après les sondages pour les présidentielles américaines l’été passé, on commence à crouler sous les sondages, en France, pour la prochaine élection présentielle. Avec les questions que ça induit. Et loin de moi l’idée de faire un billet sur les sondeurs sont des gros nuls, parce que c’est quand même plus compliqué que ça. Et il est normal de faire des erreurs, comme le rappelait  Ranganathan (2014). Et comme on m’a convié dans deux semaines à une table ronde sur les sondages, je me suis dit qu’il serait temps que je travaille un peu.

Même si ça semble être des banalités, les élections américaines et les élections en France sont très différentes.

Premier point, les 538 grands électeurs, contre un vote direct en France. Quand on estime la probabilité de gagner on obtient des choses très différentes.

On peut faire des simulations pour s’en rendre compte

simu=function(i){
n=1+2000
vote=sample(c("A","B"),size=n,replace=TRUE,prob=c(.51,.49))
region=sample(1:50,size=n,replace=TRUE)
france=mean(vote=="A")
us=aggregate(vote,by=list(region),function(x) mean(x=="A"))$x
prop=c(france,mean(us>1/2))
names(prop)=c("FR","US")
return(prop)
}
S=Vectorize(simu)(1:1e4)
plot(density(S[1,]))

En France, en interrogeant 2000 personnes, on a une très faible incertitude sur probabilités réelles, et dans 82% des scénarios, le candidats qui a 51% d’intention de vote est déclaré vainqueur.

> mean(S[1,]<.5)
[1] 0.1841

Aux Etats-Unis (ou disons une simulation d’une élection avec un collège électoral, avec ici 50 grands électeurs), c’est une autre histoire.

plot(density(S[2,]))

Dans un tiers des cas, on se trompe sur le vainqueur.

> mean(S[2,]<.5)
[1] 0.3319

Et la distribution de la couleur des grands électeurs peut être très volatile… allant de 30% à 70% pour A, alors la moitié du peuple a voté pour lui.

Deviner qui va gagner est un exercice manifestement compliqué aux Etats-Unis, à cause de ces grands électeurs. Et je ne parle pas du fait qu’estimer une probabilité très finement ne présage en rien du résultat final, comme  je l’avais dit dans un précédant billet.

Deuxième point, le nombre de candidats dans la course. Un simplifiant un peu, ils sont deux aux Etats-Unis, avec (quasiment) un vote en un tour, alors qu’en France, ils sont plus nombreux, avec un vote en deux tours, et une règle qui dit que seuls les deux premiers accèdent au second tour.

Commençons par le cas simple, avec quatre candidats, et une élection à un tour (à la majorité)

simu=function(i){
n=1+2000
vote=sample(c("A","B","C","D"),size=n,replace=TRUE,prob=c(.27,.26,.24,.23))
tb=table(vote)
return(names(which.max(tb)))
}
S=Vectorize(simu)(1:1e4)
barplot(table(S))

Pour passer au second tour, on imagine qu’il faut comprendre comment se passe le report des voix. On va supposer qu’il existe une espèce d’ordre, entre les candidats, A, B, C et D. Si son candidat ne passe pas le prochain tour, on vote pour un voisin. Autrement dit, A ne votera jamais pour C, ni pour D.  Mais si B est candidat, il votera pour lui. Par contre, si B a le choix entre A et C, on va supposer que l’électorat se partage en deux, entre A et C. Par contre, il ne votera jamais pour D.

simu=function(i){
n=1+2000
probabilite=c(.27,.26,.24,.23)
vote=sample(c("A","B","C","D"),size=n,replace=TRUE,prob=probabilite)
tb=table(vote)+rnorm(4)/1000
perdants=LETTERS[1:4][(rank(tb)%in%1:2)]
if(("A"%in%perdants)&("B"%in%perdants)) prob2=c(0,0,probabilite[2]+probabilite[3],probabilite[4])
if(("A"%in%perdants)&("C"%in%perdants)) prob2=c(0,probabilite[2]+probabilite[1]+probabilite[3]/2,0,probabilite[4]+probabilite[3]/2)
if(("A"%in%perdants)&("D"%in%perdants)) prob2=c(0,probabilite[2]+probabilite[1],probabilite[4]+probabilite[3],0)
if(("B"%in%perdants)&("C"%in%perdants)) prob2=c(probabilite[2]+probabilite[1],0,0,probabilite[4]+probabilite[3])
if(("B"%in%perdants)&("D"%in%perdants)) prob2=c(probabilite[1]+probabilite[2]/2,0,probabilite[3]+probabilite[2]/2+probabilite[4],0)
if(("C"%in%perdants)&("D"%in%perdants)) prob2=c(probabilite[1],probabilite[2]+probabilite[3],0,0)
vote=sample(c("A","B","C","D"),size=n,replace=TRUE,prob=prob2)
tb=table(vote)
return(LETTERS[which.max(tb)])
}
S=Vectorize(simu)(1:1e4)
barplot(table(S))

Le résultat est vraiment différent. Et la clé est le report des voix, entre les deux tours. Aussi, deviner qui va gagner est un exercice manifestement compliqué en France, à cause de ces deux tours.

Au delà des systèmes de votes, il existe quand même des similitudes, quand on parle d’incertitude dans les sondages. Classiquement, on pense au biais d’estimation. Parce qu’il correspond à ce que toute personne ayant fait un peu de statistiques (ce qui représentera bientôt la majorité des français, si j’en crois les réformes du programme du lycée) doit savoir (je n’ai pas dit compris). C’est le fameux « le sondage IFOP publié jeudi dans France-Soir a, quant à lui, été mené auprès de 1 046 personnes, soit une marge d’erreur d’environ 2,5 % » . Il s’agit de la marge d’erreur.

Mais ce n’est pas la seule source d’erreur, loin de là ! La base des sondages, c’est d’avoir un échantillon aléatoire de votants (au sens de personnes qui iront voter), et d’extrapoler ensuite sur la population entière, éventuellement en pondérant les réponses. Les poids sont là pour corriger de différences observables entre la population entière, et la population interrogée. Classiquement, on va utiliser des variables démographiques (âge, genre, profession, région, etc) dans la population des électeurs. Par exemple, si sur les listes électorales on a 50% d’hommes et 50% de femmes, mais que dans notre échantillon, il y a deux fois plus d’hommes que de femmes, on va généralement attribuer un poids de 2 pour les femmes, et de 1 pour les hommes afin de corriger. Malgré tout, on peut observer une différence entre les votes estimés, et la réalité le jour du vote. Une des raisons est que la population des votants peut être différente de celle des populations inscrites sur listes électorales. Si les hommes s’abstiennent beaucoup plus que les femmes, la correction utilisée ne suffira pas. C’est ce qu’on pourrait appeler une erreur de framing. Cette dernière correspond à l’erreur que l’on fait quand la population réelle et la population de l’échantillon diffère. L’exemple classique était lié à l’utilisation du téléphone pour constituer des listes de sondés. Avec cette méthode, les personnes sans téléphone (ou sur liste rouge) étaient systématiquement exclues.

Le second problème est le traitement des non-réponses. Si la non-réponse est un phénomène purement aléatoire, imprévisible, dans ce cas, ça va juste augmenter la variance. C’est ce qu’on observe sur l’exemple ci-dessous. Si tout le monde répond, on a

simu=function(i){
n=2000
vote=sample(c("A","B"),size=n,replace=TRUE,prob=c(.51,.49))
return(mean(vote=="A"))
}
S=Vectorize(simu)(1:1e4)

Supposons maintenant que 50% des gens ne répondent pas

simu=function(i){
n=2000
vote=sample(c("A","B",NA),size=n,replace=TRUE,prob=c(.51/2,.49/2,.5))
return(mean(vote[!is.na(vote)]=="A"))
}
S=Vectorize(simu)(1:1e4)

On a davantage de dispersion, mais on n’a pas de biais. Autrement dit, le nombre de cas où A est effectivement en tête sur le sondage est du même ordre de grandeur, avec ou sans non-réponse.

Par contre, si les électeurs de A ont davantage tendance à ne pas répondre, c’est une toute autre histoire

simu=function(i){
n=2000
vote=sample(c("A",NA,"B",NA),size=n,replace=TRUE,prob=c(.51*.8,.51*.2,.49*.9,.49*.1))
return(mean(vote[!is.na(vote)]=="A"))
}
S=Vectorize(simu)(1:1e4)

Autrement dit, si la non-réponse est liée au vote, on a un biais causé par la non-réponse. C’est le cas lorsque les personnes votant pour un des candidats n’osent pas révéler. C’est ce que montre très bien Fuchs, Bossert & Stukowski (2013).

Un autre exemple connu en sondage est lié à l’ordre des questions, comme le soulignait McFarland (1981), ou même les mots utilisés, pour reprendre Smith (1987). On peut rapprocher ce biais au biais dit de spécification, lorsqu’il y a une mauvaise compréhension de la question.

Enfin (et je m’arrêterais là), il reste un biais important, qui est un biais d’indécision, pour personnes indécises ou qui pourraient juste changer d’avis entre deux dates.

Compte tenu de tous ces biais, on voit qu’il est compliqué d’estimer l’erreur associée à un sondage. En tous cas, on se doute que le fameux « Le sondage IFOP publié jeudi dans France-Soir a, quant à lui, été mené auprès de 1 046 personnes, soit une marge d’erreur d’environ 2,5 % » implique une (forte) sous-estimation de l’incertitude. En fait, si on en croit Shirani-Mehr et al. (2016), la marge d’erreur pourrait être deux fois supérieure, avec un RMSE de ±7% sur l’intervalle de confiance et non pas ±3%. J’aurais donc tendance à croire que les sondages rajoutent essentiellement une information très bruitée dans le climat pré-élections.

Votes Populaires et Grands Électeurs

La semaine dernière, j’avais mis un billet en ligne où je m’interrogeais sur la pertinence de faire commencer l’axe des ordonnées en 0, quand on regardait les élections américaines, en représentant le nombre de votes obtenus. Mon point était que 0 n’avait pas vraiment de sens, quand on regardait républicains et démocrates (en 2016 en tous cas) car ces deux parties obtiendront toujours plusieurs dizaines de millions de voix, quoi qu’il se passe.

Mais les votes populaires importent finalement peu, compte tenu du mode de scrutin. Le plus importants est le nombre de grands électeurs du collège électoral.

> electoral=as.matrix(data.frame(
+ red=c(306,206,173),
+ blue=c(232,332,365)))
> electoral=electoral[3:1,]
> rownames(electoral)=c(2008,2012,2016)
> barplot(t(electoral), col=c("red","darkblue")

Faire commencer le graphique en 0 a ici davantage de sens, car il serait possible d’imaginer qu’un candidat obtienne très peu de grands électeurs. Et force est de constater que ce graphique est  parlant. On peut même tenter d’aller un peu plus loin, en superposant les deux… Le premier soucis est qu’on ne peut pas superposer deux graphiques en barres, donc il faut revenir aux lignes (ou aux points, mais comme c’est une évolution dans le temps, on peut les relier). Le second soucis est que les deux n’ont pas la même échelle (pour les grands électeurs, on en a en tout 539, alors que la population américaine se compte en dizaine de millions d’habitants). On pourrait par exemple tracer

ou encore (de manière toute aussi valide)

Bref, les doubles axes posent soucis (mais ce n’est pas nouveau). La solution serait peut être de travailler en pourcentage. Pour les grands électeurs, on peut mettre le pourcentage obtenu parmi les 539 électeurs. Pour les nombres de votes populaires, on peut se contenter (je pense) de comparer les deux (i.e. pourcentage de voix obtenues par le rouge par rapport au nombre voix obtenues par le rouge et le bleu)

C’est plus parlant, non ?

Analyse des Données et Cartes

Vendredi, on continuera en cours la classification (non supervisée) et en particulier les méthodes hiérarchiques, et les arbres. On va utiliser la base suivante, avec les résultats des élections présidentielles (premier tour) de 2012,

elections2012 <- read.csv("http://komodo.regardscitoyens.org/public/presidentielles2012/resultats/resultats_departements_final_T1.csv",
sep=";",header=TRUE,dec=",")

La base qu’on garde contient juste les voix exprimées, en pourcentage

voix <- which(substr(names(elections2012),1,12)=="X..Voix.Exp.")
X <- as.matrix(elections2012[,voix])
colnames(X) <- substr(names(elections2012)[voix],13,nchar(names(elections2012)[voix]))
rownames(X) <- elections2012[,3]

Parfois, il convient de normaliser les données. Ici, sur les département, ça n’est pas nécessaire. Par contre, sur les candidats, ça le serait (comme pour l’ACP en fait).  Pour visualiser les distances (entre départements), on peut utiliser

heatmap(X)

On peut faire une CAH, sur notre matrice de distance entre rangs

cah <- hclust(dist(X))
plot(cah,cex=.6)

et si on souhaite garder 5 groupes (par exemple) on utilise

rect.hclust(cah,k=5)
groups.5 <- cutree(cah,5)

On peut aussi visualiser les classes avec la fonction

library(dendroextras)
plot(colour_clusters(cah,k=5))

Pour comprendre qui se trouve dans nos groupes, on peut utiliser

aggregate(X,list(groups.5),mean)

qui va nous renvoyer les “votes  moyens” pour chaque candidat, dans chaque groupe. Maintenant, pour en finir avec la visualisation, on peut utiliser le code suivant (un peu long, certes) qui permettra de visualiser sur un carte les différents groupes

carte_classe <- function(groupes){
library(stringr)
elections2012$dep <- elections2012$Libellé.du.département
elections2012$dep <- tolower(elections2012$dep)
elections2012$dep <- str_replace_all(elections2012$dep, pattern = " |-|'|/", replacement = "")
library(maps)
france<-map(database="france")
france$dep <- france$names
france$dep <- tolower(france$dep)
france$dep <- str_replace_all(france$dep, pattern = " |-|'|/", replacement = "")
corresp_noms <- elections2012[, c("Libellé.du.département", "dep")]
corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud"
col2001<-groupes+1
names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms$Libellé.du.département)]
color <- col2001[match(france$dep, names(col2001))]
map(database="france", fill=TRUE, col=color)
}

carte_classe(groups.5)

Avec ces fonctions, on devrait pouvoir tester différentes méthodes pour constituer des groupes.

Voting Twice in France

On the Monkey Cage blog, Baptiste Coulmont (a.k.a. @coulmont) recently uploaded a post entitled “You can vote twice ! The many political appeals of proxy votes in France“, coauthored with Joël Gombin (a.k.a. @joelgombin), and myself. The study was initially written in French as mentioned in a previous post. Baptiste posted additional information on his blog (http://coulmont.com/blog/…) and I also wanted to post some lines of code, to mention a model that was not used in that study (more complex to analyze, but more realistic, and with the same conclusions). The econometric study is based on aggregated voted, with a possible ecological misinterpretation.

  • Regression Model: Possible Explanatory Variables

The first idea was to model proxies using a binomial regression, per pooling station  where  denote the number of proxy vote, per station , and  denotes the number of voters. Proportion  can be a function of possible explanatory variables (on Baptiste’s blog there are additional information about the datasets, obtained from insee.fr and opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

For instance, assume that  is a function of the proportion of owner of the place people live in, denoted  in the neighborhood of the pooling station,

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference]

We can consider a logistic regression

or a logistic regression with splines, if we do not want to assume a linear model

With cubic splines, the code is

> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])

while a standard logistic regression would be

> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

Here it is (the confidence region is for the spline regression) with on blue the first round of the Presidential election, and in red, the second round (in France, it’s a two-round system)

(the legend of the y axis is not correct). We can consider as explanatory variable the rate of H.L.M., low-cost housing or council housing,

If I like the graph, unfortunately, the interpretation of coefficient  might be complicated

> summary(regt1l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset1, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.9549   -1.5722    0.0319    1.6292   13.1303  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.70811    0.01516  -244.6   <2e-16 ***
taux         1.49666    0.04012    37.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12507  on 836  degrees of freedom
Residual deviance: 11065  on 835  degrees of freedom
AIC: 15699

Number of Fisher Scoring iterations: 4

> summary(regt2l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset2, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-15.4872   -1.7817   -0.1615    1.6035   12.5596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24272    0.01230 -263.61   <2e-16 ***
taux         1.45816    0.03266   44.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9424.7  on 836  degrees of freedom
Residual deviance: 7362.3  on 835  degrees of freedom
AIC: 12531

Number of Fisher Scoring iterations: 4

So we did consider a standard linear regression model, for the proxy rate, per station,

(again, either a model with splines, or a standard linear model). The code is

> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

Here is again the evolution as a function of the rate of owner of their homes,

The graph is rather close to the one before, and here, the interpretation of the summary table is more conventional,

> summary(regt1l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.9994 -0.2926  0.0011  0.3173  3.2072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.021268   0.001739   12.23   <2e-16 ***
taux        0.054371   0.004812   11.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared:  0.1326,	Adjusted R-squared:  0.1316 
F-statistic: 127.7 on 1 and 835 DF,  p-value: < 2.2e-16

> summary(regt2l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.9029 -0.4148 -0.0338  0.4029  3.4907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.033909   0.001866   18.17   <2e-16 ***
taux        0.079749   0.005165   15.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared:  0.2221,	Adjusted R-squared:  0.2212 
F-statistic: 238.4 on 1 and 835 DF,  p-value: < 2.2e-16

We have used those codes to produce the graphs mentioned in the post. But before mentioning the residuals of the multiple model we considered, I wanted to share some awesome code that produce maps (I can say that those codes are awesome since Baptiste wrote most of them).

  • Visualization of Residuals on a Map of Paris

To plot the neighborhood of the pooling stations, one more time the post on Baptiste’s blog, explains how the shapefile was obtained from cartelec.net

> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")

To visualize the proxy rate (the average of round one and round two), here is the code

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

To get nice colors, function of the rates, we use

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
>(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

and finally

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

If we consider a model with three explanatory variable, to explain the proxy rate,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

we can plot the residuals using

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
>(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

It might not be a pure random spatial noise… But we could not get better with our small set of covariates.

Le Vote par Procuration en France

La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“.

Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition.

  • Modèle de régression, et recherche de variables explications

Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que  où  est le nombre de procurations dans le bureau de vote , et où  est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que , la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

Si on suppose que  est fonction  le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote),

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference]

il est légitime de tenter une régression logistique,

voire un lissage par splines, si on pense que le lien peut ne pas être linéaire,

Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques)

> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])

et, pour la régression logistique standard (linéaire)

> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration).

On peut faire la même chose sur le taux de logement HLM, dans le quartier,

Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de  laisse à désirer (toute suggestion est la bienvenue !).

> summary(regt1l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset1, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.9549   -1.5722    0.0319    1.6292   13.1303  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.70811    0.01516  -244.6   <2e-16 ***
taux         1.49666    0.04012    37.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12507  on 836  degrees of freedom
Residual deviance: 11065  on 835  degrees of freedom
AIC: 15699

Number of Fisher Scoring iterations: 4

> summary(regt2l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset2, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-15.4872   -1.7817   -0.1615    1.6035   12.5596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24272    0.01230 -263.61   <2e-16 ***
taux         1.45816    0.03266   44.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9424.7  on 836  degrees of freedom
Residual deviance: 7362.3  on 835  degrees of freedom
AIC: 12531

Number of Fisher Scoring iterations: 4

 

On a alors voulu comparer avec un modèle qui me semble moins juste, mais qui est plus simple à interpréter, où on suppose que le taux de procuration (par bureau) est expliqué par un modèle linéaire

(que l’on peut aussi lisser, pour vérifier que le lien est effectivement linéaire). Le code est ici

> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre)

Cette fois, on a une sortie de régression plus classique,

> summary(regt1l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.9994 -0.2926  0.0011  0.3173  3.2072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.021268   0.001739   12.23   <2e-16 ***
taux        0.054371   0.004812   11.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared:  0.1326,	Adjusted R-squared:  0.1316 
F-statistic: 127.7 on 1 and 835 DF,  p-value: < 2.2e-16

> summary(regt2l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.9029 -0.4148 -0.0338  0.4029  3.4907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.033909   0.001866   18.17   <2e-16 ***
taux        0.079749   0.005165   15.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared:  0.2221,	Adjusted R-squared:  0.2212 
F-statistic: 238.4 on 1 and 835 DF,  p-value: < 2.2e-16

On note plusieurs choses de ces graphiques. (i) les deux types de régression donnent des modèles pour les taux de procuration très très proches. Donc autant prendre le plus simple à interpréter. (ii) le lissage n’apporte rien, et le modèle linéaire semble pertinent. On a ainsi regardé plusieurs variables, et on en a retenu un certain nombre, pour tenter un modèle multiple. Avant de parler des résidus de notre modèle, je devrais peut être prendre quelques lignes pour parler un peu de cartographie (une nouvelle fois, Baptiste m’a fait découvrir de belles fonctions).

  • Visualiser les bureaux de vote à Paris

Pour récupérer le fond de carte, avec les bureaux de vote, on utilise (là encore, je renvoie au blog de Baptiste, qui explique l’utilisation de données de cartelec.net)

> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")

Si on veut visualiser, par exemple, le taux de procuration (disons la moyenne entre les deux tours), on utilise les données suivantes

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur,

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

Reste à conclure, en faisant une visualisation graphique de nos données

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

Dans ce cas, la visualisation des résidus donne

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux.

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

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-2.png

et sa version contrainte

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-1.png

On cherche alors à tester l’hypothèse

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-3.png

contre l’hypothèse alternative

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-4.png

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é

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-9.png

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.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-10.png

Fisher a montré que sous l’hypothèse H0,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-11.png

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…

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)
[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)

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).

Ex-æquo aux élections ?

La saison 1 de Mad Men se termine avec l’élection de 1960, Nixon contre Kennedy, qui s’est jouée dans un mouchoir de poche.

L’autre jour, Benoit Rittaud évoquait sur http://images.math.cnrs.fr/ la possibilité que des candidats soient ex-æquo dans une élection, et les aspects légaux (car visiblement la constitution française n’a pas vraiment prévu ce cas). Le plus intéressant aurait été – à mon avis (mais mon avis est fortement biaisé) – de calculer la probabilité qu’un tel évènement survienne. Par exemple, qu’à l’issu du premier tour, le deuxième et le troisième candidat soient ex-æquo, de telle sorte qu’on ne puisse déterminer qui ira au second tour. Justement, un des commentaire affirme “la probabilité qu’un tel événement arrive est environ l’inverse du nombre typique de voix séparant le deuxième homme (ou femme) du troisième”. N’ayant réussi à démontrer ce résultat, je me suis dit qu’on pourrait se lancer dans des simulations pour quantifier cette probabilité.

Commençons par supposer qu’il y a quatre candidats (et la possibilité de s’abstenir). On suppose que le taux d’abstention est de 20%, et que les intentions de vote pour chacun des candidats soit de 20%.

EXE=function(N=100,ns=1000000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("A","B","C","D","E"),size=N,prob=rep(.2,5),
replace=TRUE)
tb=table(M)
Ms=rev(sort(tb[-(names(tb)=="A")]))
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Le soucis est que faire plusieurs centaines de millions de tirages de plus de 40 millions de lois multinomiales pourrait prendre du temps (genre vraiment beaucoup). Une stratégie intuitive peut être de faire des simulations pour une population plus faible, puis d’extrapoler.

E=function(N) EXE(N)
taille=c(20,50,100,200,500,1000,2000)
proba=Vectorize(E)(taille)
plot(taille,proba,type="b",log="xy")

base=data.frame(y=proba,x=taille)
reg=lm(log(y)~I(log(x)),data=base)
s2=summary(reg)$sigma
m=exp(predict(reg,newdata=
data.frame(x=41194689)))*exp(s2^2/2)

L’extrapolation linéaire n’étant peut-être pas parfaite, on peut tenter un ajustement quadratique (comme cela a été fait sur la figure ci-dessus)

reg2=lm(log(y)~I(log(x))+I(log(x)^2),data=base)
s2=summary(reg2)$sigma
m2=exp(predict(reg2,newdata=data.frame(x=41194689)))*exp(s2^2/2)

La différence de prédiction avec les deux modèles est significative comme on le voit sur le graphique. Toutefois, avoir des ex-æquo avec une très très grande population est tellement rare qu’il faut réellement beaucoup de simulations pour avoir un ordre de grandeur valide. L’extrapolation semble être une alternative valide ici (mais je suis ouvert à toute critique, ou discussion).

Maintenant que l’on a une méthodologie (ou presque), utilisons des ordres de grandeurs réalistes. Par exemple en se basant sur l’élection présidentielle de 2002. 41.19 millions d’électeurs, 28.4% d’abstentions, et pour les votes exprimés, un candidat en tête avec 19.9% des votes, un second à 16.8% et le troisième à 16.2% (et les autres plus loin). Les intentions de votes pour les principaux candidats figurent dans la loi multinomiale ci-dessous

EXE=function(N=100,ns=1000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("abs","JC","JMLP","LJ","FB","AL","JPC"
,"NM","autres"),size=N,prob=c(.284,.1988*.716,
.1686*.716,.1618*.716,.0684*.716,.0572*.716,
.0533*.716,.0525*.716,.1714104),replace=TRUE)
tb=table(M);
Ms=rev(sort(tb[-(names(tb)%in%c("abs","autres"))]));
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Si on en croit l’affirmation du commentaire, la probabilité d’avoir eu des ex-æquo (entre le second et le troisième candidat) serait de l’ordre de

> 1/(41194689*(.1686*.716-.1618*.716))
[1] 4.985823e-06

Je dois avouer que les calculs n’ont rien donné (on verra tout à l’heure ce qui se passe si on change la définition d’ex-æquo). Numériquement, les calculs tournent encore (après 4 jours) sur le serveur… Par contre, si la population ne contenait que 10000 électeurs, la probabilité d’avoir des ex-æquo (entre le 2ème et le 3ème candidat) serait de l’ordre de 3 sur 10000 (avec un intervalle de confiance que je n’ai toujours pas calculé)

> EXE(10000,100000000)
[1] 3.245562e-05

A titre d’illustration, on peut regarder ce qui s’est passé, commune par commune en France (les données sont en ligne sur le site du Ministère de l’Intérieur). L’exercice est biaisé parce qu’il y a de l’hétérogénéité spatiale… mais ça permet toujours d’illustrer… Les nombres de voix obtenues, pour le second et le troisième candidat, par commune est le suivant

> election=read.table(
+ "/Users/UQAM/Documents/Pres2002Tour1.csv",sep=";")
+ election=election[-1,]
> VOIX=election[,seq(14,74,by=4)]
> DEUXIEME=TROISIEME=TOTAL=rep(NA,nrow(VOIX))
> for(i in 1:nrow(VOIX)){
+ u=rev(sort(as.numeric(as.character(
+ election[i,seq(14,74,by=4)]))))
+ TOTAL[i]=sum(u)
+ DEUXIEME[i]=u[2]
+ TROISIEME[i]=u[3]
+ }

Sur 37000 communes,

> nrow(election)
[1] 37513

ou plutôt sur les 24000 communes qui ont eu plus de 5000 votes exprimés

> sum((TOTAL>5000))
[1] 24230

il y a eu 56 cas avec de deuxième et troisième candidats ex-aequo,

> base=data.frame(TOTAL,DEUXIEME,TROISIEME,COMMMUNE=election[,3])
> selection=(TOTAL>5000)
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 56

On a même le détail,

> sousbase[indice,]
TOTAL DEUXIEME TROISIEME                     COMMMUNE
1134   5819      699       699                        Crouy
1507   5386      717       717          Louroux-Bourbonnais
2603   6235     1151      1151                  Saint-Peray
2916   6346      935       935                        Coucy
2970   7529     1413      1413                   Haraucourt
5089   5567      717       717                      Vignats

(etc). On peut alors regarder, en fonction de la taille des communes (ici le nombre de votes exprimés) la probabilité d’avoir un ex-æquo, dans une commune donnée. On va subdivisionner en 50 classes (pour commencer), et compter à chaque fois le nombre d’ex-æquo,

> n=50
> niveau=seq(1/(2*n),1-1/(2*n),by=1/n)
> base$CLASSE=cut(TOTAL,quantile(TOTAL,seq(0,1,by=1/n)),
+ labels=niveau)
> base$EXAEQUO=(base$DEUXIEME==base$TROISIEME)*1
> m=by(base$EXAEQUO, base$CLASSE, sum)

Par exemple pour les 2% des plus grosses villes, on a eu 5 cas de communes avec des ex-aequos

------------------------------------------------------------
base$CLASSE: 0.99
[1] 5

Si on regarde plus en détails, ce sont les communes avec plus de 7500 votes exprimés,

> quantile(TOTAL,.98)
98%
7427.76
> selection=(TOTAL>quantile(TOTAL,.98))
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 5
> indice=which(sousbase$DEUXIEME==sousbase$TROISIEME)
> sousbase[indice,]
TOTAL DEUXIEME TROISIEME       COMMMUNE CLASSE EXAEQUO
2970   7529     1413      1413     Haraucourt   0.99       1
14359  7496     1392      1392          Renac   0.99       1
25028  8104     1381      1381         Heloup   0.99       1
30014  8087     1340      1340         Barnay   0.99       1
36770  7435      764       764 Pont-sur-Yonne   0.99       1

On peut représenter ces nombres de cas d’ex-æquo et faire une régression de Poisson (car on fait du comptage ici),

> taille=length(TOTAL)/n
> plot(quantile(TOTAL,niveau),m/taille)
> sb=data.frame(niveau,exaequo=as.vector(m))
> library(splines)
> reg=glm(exaequo~bs(niveau),family=poisson,data=sb)
> u=seq(0,1,by=.001)
> lines(quantile(TOTAL,u),predict(reg,
+ newdata=data.frame(niveau=u),
+ type="response")/taille,lwd=2,col="red")

pour des facilités de lecture, on a en ordonnées des probabilités et en abscisse les tailles moyennes des villes

On peut faire la même chose avec 200 classes,

Bref, pour une commune de taille moyenne (disons moyenne supérieure, avec 5000 votes exprimés), on a 2 chances sur 1000 d’avoir un second et un troisième candidat ex-æquo. Ce qui laisserait entendre que pour 40 millions, la probabilité devrait être faible, voire très faible…

Bon, au delà du problème conceptuel des ex-æquo (ou légal, car je doute qu’il soit possible d’avoir une égalité parfaite, le Ministère de l’Intérieur ferait un recompte et trouverait forcément un résultat différent), il semble que la probabilité semble beaucoup plus faible que celle obtenue par la règle ad-hoc….

Mais peut-être serait-il réaliste de dire que l’élection n’a pas su dégager une majorité (c’est le but d’une élection il me semble) si le nombre de voix séparant deux candidats est trop faible. Par exemple, si moins de 100 voix séparent les candidats, on refait un vote ! Pour revenir sur mon exemple précédant, avec 10000 votants

> EXE100(10000,1000000)
[1] 0.012641

soit une chance sur 100 (je me suis autorisé à faire moins de simulations car l’évènement me semble moins improbable… donc il y a besoin de moins de simulations pour l’observer à l’occasion). Si on va plus loin, et que l’on fait des simulations que l’on projette, comme auparavant, pour estimer la probabilité d’avoir des ex-æquo non pas pour une ville, mais au niveau de la France entière, on est plutôt aux alentours de

> (m=exp(predict(reg,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
1.505422e-12

pour un modèle linéaire, et

> (m=exp(predict(reg2,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
4.074588e-98

pour un modèle quadratique (avec la confiance – relativement faible – que l’on peut accorder à une telle projection). Ce qui donne une différence énorme… Si on visualise tout ça,

lines(c(1,100000000),exp(predict(reg,newdata=
data.frame(x=c(1,100000000))))*exp
(s2^2/2),lty=2,col="light blue")

que l’on représente en bleu, en vert un ajustement quadratique

lines(10^(seq(1,8,by=.1)),exp(predict(reg,newdata=
data.frame(x=10^(seq(1,8,by=.1)))))*exp
(s2^2/2),lty=2,col="dark green")

alors que le calcul suggéré dans le commentaire donnerait la courbe rouge,

lines(c(1,100000000),1/(c(1,100000000)*
(.1686*.716-.1618*.716)),col="red",lty=2)

Peut-on parler d’évènement improbable ? A titre de comparaison, quand on parle de quelque chose d’improbable, ou qui arrive avec une probabilité infinitésimale, on pense au loto. A l’ancien loto (français), il fallait choisir 6 numéros parmi 49 (sans remise). La probabilité de gagner le gros lot était de

Pour l’Euromillions, il faut choisir 5 numéros parmi 50 numéros auxquels il faut ajouter les combinaisons possibles avec les Étoiles, soit 2 numéros à choisir parmi 9. La probabilité de gagner le gros lot est alors de

Bref, tomber sur des ex-æquo lors d’une élection nationale, en France, c’est a priori plus rare que gagner le gros lot en jouant une fois au loto… On peut donc dire qu’il ne devrait pas y avoir d’ex-æquo à l’issu du premier tour !

et c’est parti pour les élections…

Ce soir, @imparibus reprenait  une information du Point (ici). “Un sondage de l’institut Sensus diffusé jeudi crédite Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7% […]  la marge d’erreur du sondage est de 2,2%, a expliqué un reponsable de l’enquête [… ce qui donne une ] grande probabilité que les 2 candidats soient à égalité”. @imparibus affirmait que c’était faux,@adelaigue me demandant mon avis…
Bon, je suis flâté qu’Alexandre m’appelle à l’aide, mais je suis plutôt nul en sondage… Bref, pour faire une réponse intelligente et détaillée, il me faudra un peu plus de place que ce que ne me propose Twitter…
Soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial01.png la probabilité de voter pour R, et https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial02.png la probabilité de voter pour S. On note https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial03.png la probabilité de voter pour un autre candidat, https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial04.png.
Si https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial05.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial06.png sont les fréquences empiriques obtenues à partir de https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png personnes sondées, rappelons que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial07.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial08.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial09.png

Aussi,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial11.png

Si on veut tester https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial12.png (hypothèse https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png) car c’est finalement de ça qu’on parle ici, l’idée la plus naturelle serait d’introduire la statistique de test d’égalité des proportions,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial14.png

qui suit – à la louche – une loi normale sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png.
C’est bien joli tout ça, mais je ne connais pas https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png ! En fait, on va essayer de s’aider de la petite phrase “la marge d’erreur du sondage est de 2,2%”
Alors si j’en crois le site d’ipsos (ici),  cette marge d’erreur correspond à 2000 personnes interrogées…
Or ici https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial05.png =46% et et https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial06.png = 42.7%, avec

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial16.png

soit,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial18.png

Si on compare à une loi normale, on a 2.6% de chances d’atteindre une telle valeur…
On peut aussi faire un test du chi-deux. En cas d’égalité des proportions, on devrait avoir https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial19.png=44.75%
et si on utilise un test du chi-deux pour mesurer la distance entre les deux lois multinomiale (sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png et sous H1), on a comme distance du chi-deux

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial21.png

soit ici https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png=3.71 avec 2000 personnes interrogées. Sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png devrait suivre une loi du chi-deux à 1 degré de liberté, ce qui donne une p-value de 5.26%.
On peut aussi tenter un test de rapport de vraisemblance. Sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png, on a

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial24.png

alors que sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial23.png

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial25.png

(en utilisant la forme des lois multinomiales). Le rapport des vraisemblances s’écrit alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial26.png

soit ici https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial27.png= 10.436. La statistique de test est icihttps://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial27.png qui doit suivre, sous https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial13.png une loi du chi-deux à un degré de liberté. La p-value associée est alors de l’ordre de 0,1%….
Bon, il semble qu’on puisse raffiner un peu les tests, en reprenant les idées de Nass (1959) et de Hoel (1938) – mais qui sont intéressant si on a peu de données (ce qui n’est pas vraiment le cas ici).
Bref, j’aurais tendance à être plutôt d’accord avec @imparibus. Maintenant, je suis nul en sondage, donc si quelqu’un a d’autres idées, je suis preneur ! Surtout que j’ai peut-être fait quelques erreurs de calculs ici ou là….

Paradoxe, vote et démocratie

Le paradoxe le plus classique sur les élections est le paradoxe de Condorcet, correspondant à une intransitivité des préférences. Pour rappels, dans le cours d’économie de l’incertain, j’avais postulé que les relations de préférences étaient une relation d’ordre (et donc étaient transitives) sans que cela ne soulève la moindre objection. Prenons un cas pratique: un sondage LH2 pour RMC-BFM TV-20 minutes, diffusé le 12 mars 2007 sur les élections présidentiels, donnait les intentions de vote suivant,

  • au premier tour, Sarkozy (28) > Royal (26) > Bayrou (22)
  • au second tour, Sarkozy (52) > Royal (48) mais  Bayrou (55) > Sarkozy (45)

Autrement dit, on peut A>B, et B>C, mais sans pour autant avoir A>C.

Pour Condorcet, il n’existait pas de système simple assurant cette cohérence, à savoir une fonction de choix social indiscutable, permettant d’agréger des préférences individuelles en préférences sociales. Arrow démontre, sous réserve d’acceptation de quelques hypothèses, qu’il n’existerait pas de système du tout assurant la cohérence, hormis celui où la fonction de choix social coïnciderait avec les choix d’un seul individu, parfois surnommé dictateur, indépendemment du reste de la population.

Le théorème d’Arrow peut être énoncé sous la forme suivante.

Pour au moins 3 options de choix et deux individus, il n’existe pas de fonction de choix social satisfaisant les propriétés suivantes :

  1. Universalité : la fonction de choix social doit être définie pour tout profil de préférences logiquement possible. Cela signifie qu’on ne peut éliminer le choix de certains individus du groupe des individus qui déterminent le choix collectif.
  2. Non-dictature : il n’existe aucun individu pour lequel la liste de ses choix personnels coïncide avec la fonction de choix social, indépendamment des préférences des autres ;
  3. Unanimité : lorsque tous les individus ont les mêmes préférences, la fonction de choix social doit associer ces mêmes préférences à la société.
  4. Indifférence des Options Non-Pertinentes : le classement relatif de deux options ne doit dépendre que de leur position relative pour les individus et non du classement d’options tierces ; si l’on ne considère qu’un sous-ensemble d’options, la fonction ne doit pas aboutir à un autre classement de ce sous-ensemble.

Sur les paradoxes du vote, je renvoie à une très bonne émission qui avait été diffusée dans Archimède sur Arte en janvier 2001. L’intégralité du texte se trouve ici,

Analyse des données, contrôle continu

Pour le contrôle continu d’analyse des données (qui ne concerne que les élèves du Master SLEB, qui concerne apparemment tout le monde *), je vous demande de travailler, au choix

  • sur l’analyse par département des votes lors de l’élection présidentielle de 2007. Vous disposez de deux feuilles xls, avec le résultat du 1er tour, et le résultat du 2nd tour. Je mets aussi à votre disposition plein de variables supplémentaires, via des données du ministère de l’intérieur (ici), en particulier les tableaux 19 et 20.
  • sur l’analyse des conditions de vie des ménages. La base est au format sas ou bien au format xls (après avoir enfin trouvé une licence SAS pour faire la conversion). Le descriptif de la base figure ci-dessous (ou ici)

En vous mettant par groupe de 2 ou 3, je vous demande un petit rapport de 15 pages (au plus), commentant ce que vous avez pu noter.

* Voilà d’ailleurs les informations dont je dispose,