Category Archives: Sondages

Les marchés prédictifs comme technique de prévision

Les dernières élections présidentielles ont remis en avant l’importance des soudages, utilisés comme outils de prévision, même si les sondeurs s’en défendent. Comme l’ont dit Niels Bohr et Pierre Dac, « la prévision est difficile surtout lorsqu’elle concerne l’avenir », et des solutions alternatives ont été envisagées. Comme c’est la mode, le « big data » a été mentionné (deviner les intentions de vote à partir des tweets ou des informations publiées sur une page Facebook), mais les marchés prédictifs (et les sites de paris) ont aussi été largement utilisés. Mais est-il si simple d’utiliser les prix (ou les cotes) sur ces marchés prédictifs pour en déduire des probabilités ? Et peut-on vraiment les utiliser ?

Sur quoi parie-t-on ?

Comme le rappelle Rhode & Strumpf (2008) utiliser les paris pour connaître les croyances et donc les probabilités qu’un évènement se réalise n’est pas nouveau. Par exemple lors des élections papales au Vatican. En 1549, Matteo Dandolo c(ambassadeur de Vénitie) notait (raconté dans Baumgartner (2003)) que « it is therefore more than clear that the merchants are very well informed about the state of the poll, and that the Cardinals’ attendants in Conclave (i conclavisti) go partners with them in the wagers, which thus causes many tens of thousands of crowns to change hands ». Et les marchés de paris lors des élections ont été populaires aux Etats-Unis, jusqu’à la seconde guerre mondiale. Rhode & Strumpf (2008) avance plusieurs raisons pour la perte d’intérêt au cours de la seconde moitié du XXème siècle : les améliorations des techniques de sondages et la légalisation des paris sur chevaux (avec pour conséquence la disparition du marché illégal des paris).

Plus récemment, avec le développement des sites et donc d’un marché des paris, on a vu apparaître des paris (contrats vendus sur le marché en ligne intrade.com, liquidé en 2015) sur la survenance de tremblement de terre de magnitude dévastatrice (plus de 9 sur l’échelle de Richter), le vainqueur de la cérémonie des Oscars, la diminution de la couverture de la glace en Arctique entre deux années, ou encore l’observation du boson de Higgs. Polgreen et al. (2007) évoquent des paris sur la survenance d’épidémies de grippe. Un autre exemple – que l’on va détailler un peu – est celui du passage d’une loi sur l’énergie (et l’instauration de mesures contraignantes de limitation des émissions de gaz à effet de serre) aux Etats-Unis. Le titre permettant de toucher 1USD en cas d’adoption du projet de loi a été lancé en mai 2009. Le 26 juin 2009, la Chambre des Représentants a adopté (à une courte majorité) le projet de loi (avec pour objectif de réduire de 17 % d’ici à 2020 par rapport à 2005 et de 80 % d’ici à 2050 les émissions de CO2). A ce moment, le titre était échangé pour 50¢, comme le montre le Graphique 1. Il restait à faire approuver la loi par le Sénat. Six mois plus tard, en signant l’accord de la conférence de Copenhague en décembre 2009, les Etats-Unis s’étaient davantage engagés, et un an plus tard (début 2011), le Sénat serait renouvelé lors des élections de mi-mandat. Le gouvernement de Barack Obama espérait faire approuver la session “lame duck” du congrès, les élections de mi-mandat devant probablement renforcer la majorité Républicaine. Mais le jeudi 22 juillet, le chef de file Démocrate au Sénat Harry Reid a annoncé le retrait du projet de loi, faute de majorité pour l’approuver. Ce jour-là, le titre s’échangeait pour 10¢.

Figure 1 : Evolution du prix du passage du contrat sur le passage du ‘American Clean Energy and Security Act’ devant le Sénat (2009-2010)

Source : auteur, données Meng (2017)

Voir un prix comme une probabilité

Avant toutes choses, rappelons qu’une propriété fondamentale des probabilités est la somme des probabilités d’un évènement et de son complémentaire doivent sommer à un. C’est assez logique si les probabilités sont vues comme des « fréquences », correspondant au nombre de cas favorables sur le nombre de cas possibles. Cette interprétation devient plus délicate quand on imagine le cas de variables aléatoires continues, et un passage par la théorie de la mesure s’impose (en particulier pour admettre qu’un « évènement de probabilité nulle » puisse malgré tout survenir). Mais historiquement, d’autres définitions des probabilités avaient été proposées, en particulier par Jacques Bernoulli, qui avait introduit une notion de « probabilité morale » (comme le rappelle Johnson (2017)) que l’on verra réapparaître dans les travaux de de Finetti sous le nom de « probabilité subjective », ou à la même époque, dans les travaux de Frank Ramsey. Pour ce dernier, une probabilité mesure un degré de croyance – « a degree of belief », qui pouvait se mesurer au travers de paris – « through betting odds » – et les prix que les joueurs sont disposés à payer donnent une information quant à ces croyances. Il devient alors possible de déduire des probabilités à partir de prix.

L’idée n’était pas nouvelle, puisqu’elle avait été émise par dans Van Rekeningh in Spelen van Gelucken publié par Christiaan Huyghens en 1655 (en latin sous le titre ‘De Ratiociniis in Aleæ Ludo‘), et utilisée quelques années plus tard, en 1671 par Wilhelmina de Witt dans le contexte d’annuités. Le prix d’un contrat versant une rente jusqu’au décès pouvant être vu à comme une moyenne pondérée d’annuités (à maturité fixe), en observant les prix des différents contrats d’assurance, on pouvait extraire des probabilités interprétées comme des probabilités de survie. Ces mesures de probabilités n’étaient plus construites à partir de fréquences observées, mais construites à partir de prix historiques observés. On retrouvera cette idée en mathématique financière presque trois siècles plus tard, avec le théorème fondamental de valorisation, et les prix d’Arrow-Debreu. Cette littérature est riche et abondante (on peut renvoyer à Charpentier (2014)), mais il existe également des résultats moins connus, sur la création d’un consensus à partir de probabilités, et les applications pour les paris.

Travaillant sur le marché de pari hippique, Eisenberg & Gale (1959) ont obtenu des résultats relativement généraux, qui peuvent se formaliser ainsi : considérons un individu i\in\{1,2,\cdots,I\} qui a la possibilité de parier sur un cheval j\in\{1,2,\cdots,J\} lors d’une course. Chaque joueur possède une somme totale b_i (que l’on normalise de telle sorte que b_1+b_2+\cdots+b_I=1, et qu’il place la somme \beta_{i,j} sur le cheval . Ces montants vérifient la contrainte \beta_{i,1}+\cdots+\beta_{i,J}=b_i , et on note \pi_j la somme total placée sur le cheval j, de telle sorte que \beta_{1,j}+\cdots+\beta_{I,j}=\pi_j. Compte tenu de la contrainte de budget, notons que \pi_1+\cdots+\pi_J=1, de telle sorte que les montants \pi_j peuvent être interprétés comme des probabilités. Supposons maintenant que chaque joueur ait des croyances (que chaque cheval gagne) matérialisées par un vecteur de probabilités \mathbf{p}_i=(p_{i,1},\cdots,p_{i,J}). Eisenberg Gale (1959) définissent alors une situation d’équilibre, dès lors que si \beta_{i,j}>0 alors on doit avoir

p_{i,j}=\pi_j \cdot \max_s\left\{\frac{p_{i,s}}{\pi_s} \right\}" Il montre qu’un tel équilibre existe, et qu’il est unique. On retrouve ici un résultat connu en mathématique financières sous le nom de loi du prix unique, même si les croyances ne sont pas explicitement mentionnées. Et c’est ce résultat qui permet d’interpréter les prix sur un marché prédictif comme des probabilités.

En s’inspirant du modèle de Eisenberg Gale (1959), Manski (2005) obtient une interprétation sur un marché prédictif. Dans ces marchés, on ne peut pas parier sur un cheval dans une course, mais sur la survenance d’un évènement. Considérons un contrat offrant 1$ si un évènement A survient (élection d’un nouveau pape, départ d’un pays de la zone euro, etc), proposé au prix \pi_A. S’il existe un contrat sur l’évènement complémentaire \bar A (on va supposer ici que c’est le cas), il devrait être proposé au prix \pi_{\bar A}=1-\pi_A, sinon un arbitrage serait possible[i]. Considérons maintenant un agent i croyant avec probabilité [lap_i>\pi_A, cette personne a intérêt à investir tout son argent pour acheter son titre, sinon, il a intérêt à investir dans l’autre titre, associé à \bar A. Au niveau agrégé la demande pour les titres de type A et \bar A sera respectivement

\frac{1}{\pi_A}\sum_i b_i P(p_i>\pi_A) et
<br /> \frac{1}{\pi_{\bar A}}\sum_i b_i P[p_i<\pi_A] et on aura un équilibre si \sum_i b_i=\frac{1}{\pi_A}\sum_i b_iP[p_i>\pi_A]=\frac{1}{\pi_{\bar A}}\sum_i b_iP[p_i<\pi_A]Si on suppose que les croyances (p_i) et les richesses (b_i) sont indépendantes, on peut interpréter les prix \pi_A et \pi_{\bar A} comme des probabilités, \pi_{A}=P[p_i>\pi_A] et \pi_{\bar A}=P[p_i<\pi_A] si les agents sont homogènes. Wolfers & Zitzewitz (2007) propose une interprétation différente, plus économique, dans le contexte d’agents ayant une utilité logarithmique.

Marché et pari en ligne

La situation que nous venons de décrire est celle d’un marché parfait, où un équilibre se fait entre les parieurs, et cet équilibre donne des prix interprétables comme des probabilités de gain (on parle éventuellement de « pari mutuel »). Mais la pratique des paris en ligne est souvent différente, car les prix sont proposés par un bookmaker. Comme le nom l’indique, ces personnes – à la base – enregistrent les paiements, mais en pratique, ce sont aussi eux qui fixent les prix. Par exemple lors des dernières élections présidentielles en France, plusieurs sites proposaient de parier sur l’issue du second tour.

Tableau 1 : côtes lors des élections présidentielles de 2017

côtes (sites de paris en ligne) Intention de votes (instituts de sondage)
bet365 comeon betsafe FT Bloomberg BBC
Macron 1/8 1/8 1/7 64% 60.5% 61%
Le Pen 5/1 9/2 9/2 36% 39.5% 39%

Source : http://www.onlinebettingsites.com/2017-french-election/ (données fin avril 2017)

Comme le montre le Tableau 1, les côtes offertes par les bookmakers sont différentes entre les sites de paris. Une cote de 1/8 signifie qu’en pariant 8€, mon gain net sera de 1€ si le candidat l’emporte le jour de l’élection (et je perds ma mise si le candidat perd). La « probabilité » associée est ici \pi_{A}=8/9\sim88.89\% pour E. Macron (mon gain brut étant de 9/8€ pour 1€ misé) et \pi_{\bar A}=1/6\sim16.67\%. Le soucis – on retrouve ici la discussion initié par Jacques Bernoulli – est que ces « probabilités » ne somment pas à 1, puisque 8/9+1/6\sim1.055\%. La différence vient du fait qu’il ne s’agit pas un « juste » prix, le bookmaker s’assurant un rendement certain (de l’ordre de 5,5% de la mise initiale, quel que soit le site choisi). Les vraies probabilités induites par ces cotes sont ici respectivement \pi_A=16/19\sim 84.2\% et \pi_{\bar A}=3/19\sim 15.8\%. Les sites de paris en ligne donneraient une probabilité de gain pour Emmanuel Macron de 84%, alors que les sondages indiquent une intention de vote de 62%, pour les instituts de sondage.

A-t-on réellement une probabilité, ou une fourchette ?

On le voit, les grandeurs données ne sont pas identiques, l’une indiquant une croyance sur une probabilité de gagner, et l’autre une proportion de votes que devrait avoir un candidat. Dans le modèle de Manski, les probabilités induites par la moyenne des croyances et les prix coïncident si on suppose que les agents ont une utilité logarithmique. Comme le montre Wolfers & Zitzewitz (2006), si on suppose que les agents ont une aversion au risque, ça ne sera plus le cas. En particulier, l’article suggère de considérer des bornes pour \pi_A de la forme [\pi_{A}^2,2\pi_A-\pi_A^2, représentées sur la Figure 2, où est appelé « prix de marché » pour reprendre une notation proche de la Figure 1.

Figure 2 : Bornes sur les croyances moyennes

Faut-il utiliser les prédictions des marchés prédictifs ?

On l’a vu lors de la dernière élection en France, beaucoup s’interrogent sur la place pris par les sondages, en particulier en notant qu’il ne s’agissait pas d’une photographie neutre, mais que les électeurs pouvaient choisir leur candidat précisément à la lecture des sondages (la fameuse théorie du « vote utile » en particulier). Le danger à utiliser les marchés prédictifs comme outils de prévision est que ces dangers sont manipulables, comme l’a montré Schröder (2008). Un autre danger – et cela a été noté lors de la dernière élection présidentielle aux Etats-Unis – les côtes sur les sites en ligne de pari reflètent les croyances d’une certaine « élite » (disons les personnes souvent relativement éduquées, et la classe supérieure, suivant la politique), contrairement aux sondages, qui sont supposés être conduits sur un échantillon représentatif des électeurs, et sont supposer refléter les choix du « peuple » (et non pas une opinion ou une croyance). Si les parieurs tiennent comptent de ces sondages, les prix reflètent leurs croyances. Une dernière critique que l’on peut faire est que la majorité des marchés prédictifs représentent des volumes de transactions très faibles.

Les sondages traditionnels font face aujourd’hui à une défiance croissante. Il devient difficile de constituer un échantillon représentatif, le taux de non-réponse atteint des proportions jamais atteint (et nul ne sait si cette non-réponse est corrélée avec des choix électoraux). Et la confusion règne, entre les sondages classiques, et les « sondages en lignes » hébergés sur divers sites d’information. Mais pour reprendre l’anecdote raconté par Scott Armstrong dans Knowledge@Wharton en 2004, si dans une réunion du conseil exécutif d’une grande entreprise on demande à tous les participants de faire des prévisions quant aux résultats pour l’année prochaine, ou si on demande aux participants de parier sur les résultats de l’année prochaine, quels prévisions seront les plus crédibles ?

Références

Arrow, K.J. & Debreu, G. Existence of an equilibrium for a competitive economy, Econometrica, vol. 22,‎ 1954

Baumgartner, F. Behind locked doors: a history of papal elections. 2003

Charpentier, A. L’efficience des marchés : hypothèse de modèle ou fait stylisé ? Risques, 96. 2014

Frank, E., Verbeek, E. & Nüesch, S. Inter-Market Arbitrage in Betting. Economica, 2012.

Johnson, T. Ethics in Quantitative Finance: a pragmatic theory of markets. Wiley, 2017.

Knowledge@Wharton, How Credible Are Polls? Is There a Better Way to Predict Outcomes in Politics and Business, 3 Novembre 2004,

Manski, C.F. Interpreting the Predictions of Prediction Markets. Working Paper 10359, 2004

Meng, K.C. Using a Free Permit Rule to Forecast the Marginal Abatement Cost of Proposed Climate Policy. American Economic Review vol 107, 2017

Polgreen. Use of Prediction Markets to Forecast Infectious Disease Activity. Healthcare Epidemiology, 2007

Ramsey, F.P. Truth and probability, dans The Foundations of Mathematics and other Logical Essays, 1931 .

Rhode, P.W. & Strumpf, K. (2008) Historical Election Betting Markets: an International Perspective.

Schröder, J. Manipulations in Prediction Markets. Universitätsverlag Karlsruhe, 2009.

The Economist, The Future of Futurology, 15 Novembre 2007

Wolfers, J. & Zitzewitz, E. Interpreting Prediction Market Prices as Probabilities. Working Paper 12200, 2006

Wood, T. Do betting markets outperform election polls? Hardly. The Washington Post, Août 2016

[i] Comme l’a montré Franck et al. (2012), comme il existe aujourd’hui plusieurs sites de paris (plusieurs marchés), il peut exister des arbitrages entre les marchés, rarement au sein d’un marché.

Modèle de régression et interaction(s) entre facteurs

Dans un modèle de régression, on veut écrire

Quand on se limite à un modèle linéaire, on écrit

ou encore

Mais on de doute que l’on rate quelque chose… en particulier, on va rater toutes les interactions possibles. On peut croiser les variables, et supposer que

qui peut s’étendre d’avantage, à l’ordre 3,

voire davantage.

Supposons que nos variables  soient ici qualitatives, et plus précisément binaires. Prenons un exemple simple, avec des données (classiques) en risque de crédit1. On peut trouver la base via

library(evtree)
db=GermanCredit

ou encore directement

myVariableNames = c("checking_status","duration","credit_history",
"purpose","credit_amount","savings","employment","installment_rate",
"personal_status","other_parties","residence_since","property_magnitude",
"age","other_payment_plans","housing","existing_credits","job",
"num_dependents","telephone","foreign_worker","class")

GermanCredit = read.table(
"http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data",
header=FALSE,col.names=myVariableNames)

Retenons pour commencer trois variables explicatives,

db=data.frame(Y=GermanCredit$class-1,
X1=GermanCredit$checking_status%in%c("A12","A13"),
X2=GermanCredit$credit_history%in%c("A30","A31"),
X3=GermanCredit$savings%in%c("A61","A62"))
reg=glm(Y~X1+X2+X3,data=db,family=binomial)
summary(reg)

La régression sans interaction donne ici

Call:
glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5431  -0.8421  -0.6295   1.3994   1.9999  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.8544     0.1699 -10.915  < 2e-16 ***
X1TRUE        0.3363     0.1496   2.249   0.0245 *  
X2TRUE        1.3462     0.2347   5.735 9.76e-09 ***
X3TRUE        1.0001     0.1787   5.596 2.19e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1221.7  on 999  degrees of freedom
Residual deviance: 1143.6  on 996  degrees of freedom
AIC: 1151.6

Number of Fisher Scoring iterations: 4

Il existe plusieurs interactions possibles ici (limitons nous aux paires). C’est ce que l’on observe quand on fait la régression

reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial)
summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3, family = binomial, 
    data = db)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5369  -0.8281  -0.6439   1.3954   1.9638  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -1.77109    0.20070  -8.825  < 2e-16 ***
X1TRUE         0.30296    0.33737   0.898 0.369186    
X2TRUE         0.88353    0.54255   1.628 0.103421    
X3TRUE         0.87709    0.22583   3.884 0.000103 ***
X1TRUE:X2TRUE -0.37917    0.49343  -0.768 0.442225    
X1TRUE:X3TRUE  0.09178    0.37278   0.246 0.805522    
X2TRUE:X3TRUE  0.80923    0.58185   1.391 0.164293    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1221.7  on 999  degrees of freedom
Residual deviance: 1141.0  on 993  degrees of freedom
AIC: 1155

Number of Fisher Scoring iterations: 4

On peut faire un dessin pour visualiser les interactions : on a trois sommets (nos trois variables), et on visualiser les interactions

indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",
xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

ce qui donne ici, pour nos trois variables

Ce modèle pourrait sembler incomplet, car on ne regarde que les interactions entre les modalités, par paires. En fait, c’est parce qu’il manque (visuellement) les variables non-croisées. On peut les rajouter si on veut (au risque de surcharger le dessin)

cercle=function(c,r,cl) lines(c[1]+r*cos(seq(0,2*pi,length=501)),
c[2]+r*sin(seq(0,2*pi,length=501)),col=cl)

reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=db,family=binomial)
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

soit ici

Si on change le ‘sens‘ de nos variables (en recodant a l’envers, en permutant les vrais et les faux), on obtient le graphique suivant

dbinv=db
dbinv[,2:k]=1-dbinv[,2:k]
reg=glm(Y~X1+X2+X3+X1:X2+X1:X3+X2:X3,data=dbinv,family=binomial)
indices=cbind(c(1,2,3),c(1,1,2),c(2,3,3))
k=3
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

qui peut alors être comparé au graphique précédant

Avec 5 variables, on augmente les interactions possibles… même si beaucoup risquent d’être non-significatifs. On peut déjà se focaliser sur les paires possibles d’interactions croisées. Pour simplifier le code, on va utiliser deux fonctions locales,

vrepeach=function(x,e){
v=NULL
for(i in 1:length(e)){v=c(v,rep(x[i],each=e[i]))}
return(v)}
vreplength=function(x,l){
v=NULL
for(i in 1:length(l)){v=c(v,x[l[i]:length(x)])}
return(v)}

et ensuite, on adapte le code précédant

indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

ce qui donne un schéma plus complexe,

On peut aussi prendre juste 2 variables, prenant 3 et 4 modalités respectivement. On va extraire deux variables indicatrices pour la première (la modalité restante sera la modalité de référence) et trois pour la seconde,

db=data.frame(Y=GermanCredit$class-1,
X1=GermanCredit$checking_status=="A12",
X2=GermanCredit$checking_status=="A13",
X3=GermanCredit$checking_status=="A14",
X4=GermanCredit$employment%in%c("A72","A73"),
X5=GermanCredit$employment%in%c("A74","A75"))
k=5
indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}
for(i in 1:k){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

On voit que plusieurs interactions ne sont alors plus possibles, sur la partie gauche (les trois modalités de la même variable) et sur la partie droite

On peut d’ailleurs simplifier les graphs, en ne visualisant que les interactions significatives.

indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
if(summary(reg)$coefficients[1+k+i,4]<.1){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}}
for(i in 1:k){
if(summary(reg)$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

soit ici

Ici, une seule interactions croisée est significative, et presque toutes les variables le sont. Et si on reprend le modèle avec 5 facteurs,

db=data.frame(Y=GermanCredit$class-1,X1=GermanCredit$checking_status%in%c("A12","A13"),
X2=GermanCredit$credit_history%in%c("A30","A31"),
X3=GermanCredit$savings%in%c("A61","A62"),
X4=GermanCredit$employment%in%c("A71","A72"),
X5=GermanCredit$other_payment_plans=="A143")

indices=cbind(1:(k*(k-1)/2),vrepeach(1:(k-1),(k-1):1),vreplength(2:k,1:(k-1)))
formule="Y~1"
for(i in 1:k) formule=paste(formule,"+X",i,sep="")
for(i in 1:nrow(indices)) formule=paste(formule,"+X",indices[i,2],":X",indices[i,3],sep="")
reg=glm(formule,data=db,family=binomial)
theta=pi/2+2*pi*(0:(k-1))/k
sommetX=cos(theta)
sommetY=sin(theta)
plot(sommetX,sommetY,cex=1,axes=FALSE,xlab="",ylab="",xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
for(i in 1:nrow(indices)){
if(!is.na(coefficients(reg)[1+k+i])){
if(summary(reg)$coefficients[1+k+i,4]<.1){
segments(sommetX[indices[i,2]],sommetY[indices[i,2]],
sommetX[indices[i,3]],sommetY[indices[i,3]],col="grey")
text(mean(sommetX[indices[i,2:3]]),mean(sommetY[indices[i,2:3]]),
trunc(10000*coefficients(reg)[1+k+i])/10000)
}}}
for(i in 1:k){
if(summary(reg)$coefficients[1+i]<.1){
cercle(c(cos(theta)[i]*1.18,sin(theta)[i]*1.18),.18,"grey")
text(cos(theta)[i]*1.35,sin(theta)[i]*1.35,
trunc(10000*coefficients(reg)[1+i])/10000)
}}
points(sommetX,sommetY,cex=6,pch=19,col="yellow")
points(sommetX,sommetY,cex=6,pch=1)
text(sommetX,sommetY,1:k)

on obtient

Je ne sais pas si mes graphiques sont pertinents, ou pas. Mais je trouve ça joli. En fait, je suis tombé un peu par hasard2 sur les Tables de Taguchi, développées par Gen’ichi Taguchi (田口 玄一). Le soucis est que je n’ai rien compris… Enfin, disons que je croyais comprendre, puis j’ai continué à faire des dessins… Si quelqu’un pourrait m’expliquer sur mon exemple les graphiques de Taguchi, je suis preneur ! car je doute que ce soit ce que je fais depuis tout à l’heure…

1. Cette base est largement utilisée dans le quatrième chapitre de Computational Actuarial Science with R, à paraître dans les mois à venir.

2.En l’occurence, le hasard est @Benavent qui a suscité ma curiosité ce matin en me parlant de ces tables, dont je n’avais alors jamais entendu parlé ! J’avais même lu rapidement Taniguchi (谷口 ジロー) et je ne voyais pas le rapport avec les statistiques….

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 !

Sondages et prévisions de séries temporelles

Ce soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant

sondage=read.table("http://freakonometrics.blog.free.fr/
public/data/sondage1ertour2012.csv",
header=TRUE,sep=";",dec=",")
sondage[36,1]="11/03/12"
sondeur=unique(sondage$SONDEUR)
sondage$date=as.Date(as.character(sondage$DATE),
"%d/%m/%y")

(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
couleur=c(CL[1],CL[6],"orange","grey")
datefin=as.Date("2012/04/12")
k=1
plot(sondage$date,sondage[,k+2],col=couleur[k],
pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40),
xlim=c(min(sondage$date),datefin))
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
}

Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste

rl=lowess(sondage$date,sondage[,k+2])

Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser),

D=data.frame(date=min(sondage$date)+0:148)
library(splines)
k=1
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[3],border=NA)
lines(D$date,prl,lwd=2,col=CL[1])
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,)
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
}

avec, là encore, la commande suivante pour lisser

rs=lm(sondage[,k+2]~bs(date),data=sondage)

Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets),

http://freakonometrics.hypotheses.org/files/2015/12/adic10.gif

On a ici la courbe suivante

avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

k=2
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[4],border=NA)
lines(D$date,prl,lwd=2,col=CL[6])
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7,)
for(k in c(1,3:4)){
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7)
}

Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête,

Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“.

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait…

Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via

D=data.frame(date=min(sondage$date)+0:148)
prl=predict(rs,newdata=D)

On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections)

http://freakonometrics.hypotheses.org/files/2015/12/ic-sondage01.gif

On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

library(mnormt)
simulation=function(S){
proba=c(S,100-sum(S))/100
variance= -proba%*%t(proba)
diag(variance)=proba*(1-proba)
variance=variance/1000
return(rmnorm(1,proba[1:4],variance[1:4,1:4]))}
simulsondages=function(M){
Z=rep(NA,ncol(M))
for(i in 1:nrow(M)){
Z=rbind(Z,simulation(M[i,]))
}
return(Z[-1,])}
prediction4=matrix(NA,nrow(D),4)
for(k in 1:4){
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prediction4[,k]=predict(rs,newdata=D)
}

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

set.seed(1)
S4=simulsondages(prediction4)
S100=100*S4
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}

(les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

library(forecast)
k=1
X=S100[,k]
ETS=ets(X)
F=forecast(ETS,h=60)
fdate=max(D$date)+1:60
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[3],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[2],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[1])

Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes

ETS=ets(X)
F=forecast(ETS,h=60)

Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

k=2
X=S100[,k]
ETS = ets(X)
F=forecast(ETS,h=60)
k=1
plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[4],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[5],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[6])

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

set.seed(1)
for(sim in 1:20){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
lines(fdate,as.numeric(F1$mean),col=CL[1])
lines(fdate,as.numeric(F2$mean),col=CL[6])}

Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour

set.seed(1)
VICTOIRE=rep(NA,1000)
for(sim in 1:1000){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
VICTOIRE[sim]=(F1$mean[30]>F2$mean[30])}
mean(VICTOIRE)

Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes !

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à….