Category Archives: Statistics

Les sondages avaient tout faux en 2016? Ah bon?

Allez, un dernier billet sur les sondages, d’autant plus que le débat de demain soir devrait porter (plus spécifiquement) sur les faillites sur sondages anglais et américains, de 2016.

  • Alors comme ça les sondages seraient “faux” ?

Je n’aime pas cette terminologie. J’ai l’impression que la confusion vient en grande partie du fait qu’on analyse, au final, un vote aux élection comme un référendum. C’est pour ça qu’on a tendance à mettre dans le même parnier le référendum britannique et les élections américaines pour parler de la “faillite des sondages“. Un point amusant est que si classiquement, on essaye de quantifier le nombre de voix obtenues, aujourd’hui, on communique davantage sur la “probabilité de gagner”. Un peu comme en météo quand on parle de la probabilité qu’il pleuve (PDP, qui est un concept qui me dépasse complètement, comme j’en ai parlé voilà plusieurs années sur ce blog).

Mais revenons sur ce qui s’est passé, un instant. La veille du vote américain, on nous indiquait que la probabilité que Donald Trump gagne était de l’ordre de 15%. Ce qui est la probabilité de tomber sur 6 en lançant un dé. Ou de tomber sur un somme faisant 6 en lançant deux dés. Comment dire… Je me souviens de la réaction de mon fils, il y a plusieurs années, qui s’énervait lorsqu’on jouait à Serpents et Echelles (ou au jeu de l’oie, ou aux petits chevaux) et que, par hasard, je faisais un 6, et que c’est effectivement ce dont j’avais besoin pour gagner. “Tu as triché papa !”. Ai-je besoin de faire un billet de blog pour expliquer qu’avoir un 6 en lançant un dé, “ça arrive” ? Comme je le disais dans un billet l’autre jour, on retrouve ici des débats classiques en gestion des risques. Si je prédis le décès de ma belle-mère avec 1 chance sur 1000 ou même 10000, et qu’elle décède, avais-je raison, ou pas ? C’est la multitude des prévisions qui permet de juger de la qualité de la prévision. Si dans une population importante, la proportion de décès est de 1 sur 1000, alors mon estimation était bonne. Mais les élections, on ne les observe qu’une fois…

Continuons un instant notre réflexion… Ce matin, en partant, la météo m’annonçait une “probabilité de pluie de 15%”. Que dois-je faire de cette information ? Dois-je arrêter de la consulter, le matin, car elle s’est trompé puisqu’il a plu. S’est-elle d’ailleurs trompée ? Et ai-je d’ailleurs le droit de comparer la météo et un vote ? Si on y réfléchi un peu non. Si la météo m’indique une “probabilité de pluie de 15%”, je peux légitimement décider de remettre mon projet de faire en longue balade en vélo toute la journée à plus tard. Et ma décision n’aura pas d’impact sur le fait qu’il pleuve ou pas (sauf si, comme moi, vous avez observé que justement, il pleut précisément les jours où vous n’avez ni parapluie ni manteau imperméable…). Mais pour les votes ? Si le jour du vote on m’annonce que mon candidat a 15% de chances de gagner, et que je décide de ne pas sortir pour aller voter, ça n’est pas la même chose ! Le fait qu’il gagne, ou pas, dépend en grande partie de mon action ! Les gens peuvent changer d’avis, et prendre des décisions en fonction de ce qu’on leur annonce. Et ça aura un impact sur l’élection (si beaucoup de monde raisonne ainsi). Alors que pour la pluie….

  • Une “probabilité de gagner” ?

Prévoir l’avenir c’est dur. Ce que j’aime bien sur certains sites (je pense à 538) c’est par exemple la nuance qui est faite entre les forecast (prévision pour un scrutin à venir, dans plusieurs semaines) et les “now-cast“, correspondant à une projection de ce qui se passerait lors d’une hypothétique élection qui se tiendrait lors du sondage. Mine de rien, c’est une nuance assez fondamentale, je pense. Maintenant, un point très important est, me semble-t-il, l’importance donnée à cette “probabilité de gagner”.

Dans mon prédédant billet, je faisais quelques simulations pour montrer que la probabilité de gagner est un concept difficile à comprendre, car il reflète essentiellement une incertitude. Au lieu de faire des simulations, faisons plutôt quelques calculs. Dans le cas du référendum anglais (histoire de simplifier un peu l’interprétation), supposons qu’un sondage donne “Remain” à 52%. Si on veut un intervalle de confiance à 95%, on a (en gros) un intervalle de ±3 points, soit [49%,55%]. Dans la terminologie statistique, on aura une loi normale, centrée sur 52%, avec un écart-type de 1.5%. La probabilité d’être au dessus de 50% est ici de 91%

> 1-pnorm(50,52,1.5)
[1] 0.9087888

ce qui donne une cote à 10 contre 1

> (1-pnorm(50,52,1.5))/pnorm(50,52,1.5)
[1] 9.963564

Si on avait obtenu 51%, c’est à dire un intervalle de confiance [48%,54%], on aurait une côte à 3 contre 1.

> (1-pnorm(50,51,1.5))/pnorm(50,51,1.5)
[1] 2.960513

Autrement dit, si entre deux jours consécutifs, les sondages passent de 51% à 52% (un hausse de un point), la côté passe de 3 contre 1 à 10 contre 1 ! (je parle ici de côte car beaucoup ont utilisé les paris en ligne pour prévoir la victoire de tel ou tel camp). Autrement dit, les paris (ou les probabilités qu’un camp l’emporte) sont très très sensibles à de très faibles variation dans le modèle. Et ceci n’est pas sans poser de problème quand on sait comment fonctionne le lissage des sondages (très fort au début puis de plus en plus faible au fur et à mesure que le jour de l’élection approche): la probabilité de gagner devient incroyablement volatilé en fin de course, pouvant passer de 5 contre 1 pour un candidat à 5 contre 1 pour l’autre, si deux sondages varient de 3 points (51.5% pour le premier, 48.5% pour le second).

Historiquement, les sondages donnaient une proportion, genre 52%. Puis il y a quelques années, l’intervalle de confiance est apparu. Et là, j’ai l’impression que les journalistes ont du comprendre qu’ils y perdaient. Par ce que dire que la proportion est entre 49% et 55% (avec une probabilité de 95%), ça n’est pas super sexy. La parade semble être de parler de “la probabilité de gagner” qui permet d’avoir des chiffres très variables, très grands ou très petits. Bref, des chiffres télégéniques, presque clairs (“probabilité de gagner”, tout le monde comprend, hein ?), et souvent énorme, par exemple ici 10 contre 1 (ou 91% dans une terminologie plus probabiliste). Mais si on essayait de mettre un intervalle de confiance sur cette probabilité de gagner ? A la louche, quand on a une fréquence de l’ordre de 52%, la probabilité de gagner est de 91%, mais l’intervalle de confiance est entre 25% et 100%

> 1-pnorm(50,52+c(-2,+2)*1.5,1.5)
[1] 0.2524925 0.9995709

Avec une fréquence de 50.5% obtenue lors d’un sondage, j’aurais un intervalle de confiance à 95% entre 5% et 99%

> 1-pnorm(50,50.5+c(-2,+2)*1.5,1.5)
[1] 0.04779035 0.99018467

Bref, cette probabilité de gagner est assez incertaine… non ? Dit autrement, cette histoire de “probabilité de gagner” est une belle arnaque…

  • Que s’est-il passé en Angleterre ?

Oui, les sondages se sont trompé, dans le sens où il y a eu une différence (relativement) importante entre les résultats de l’élection et des sondages. Dans le premier cas, “leave” à gagner alors que les sondages semblaient avoir tendant à croire que “remain” aurait une majorité. Bref, globablement, on prévoit 0 et on obtient 1. C’est une grosse erreur. Maintenant, si on regarde un peu en détails, comme toujours il peut y avoir plusieurs explications à cette différences entre des sondages et les résultats finaux, comme le notait par exemple Andrew Gelman

  1. les personnes qui ont répondu aux sondages n’étaient pas un “échantillon représentatif” des électeurs (et peu importe la raison au fond, mais si les deux populations diffèrent, il ne faut pas être surpris d’avoir une différence)
  2. les réponses aux sondages permettent mal de mesurer les intensions de votes (par exemple si le taux d’indécision est trop fort, et non aléatoire)
  3. la décision a été prise le dernier jour
  4. la participation électorale a été différente de celle estimée (par exemple des personnes qui voulaient voter “remain” qui se sont abstenus, persuadé de la victoire de leur camp)
  5. la faute à pas de chance, ce qu’on appelle en statistique l’incertitude d’échantillonage.

J’avais fait un billet sur ce dernier point, rappelant qu’il est beaucoup plus important que ce qu’on apprend souvent en statistique, justement à cause des autres effets. Et c’est probablement ce qui s’est passé en Angleterre, un mélange de toutes ces causes. Ah si, peut être juste un point que je pourrais mentionner ici : le type de sondage utilisé a eu un impact considérable, les sondages par téléphone mettant le “remain” bien plus haut que les sondages en ligne.

  • Que s’est-il passé aux Etats-Unis ?

Comme je l’ai déjà dit, prévoir le vote aux Etats-Unis est compliqué, compte tenu du mode de scrutin. Mais au moins, les modèles utilisés sont clairement décrits (on pourra lire le guide publié par 538 en juin 2016). La modélisation permettant de faire une prévision se fait (en gros) en quatre temps,

  1. Collecter les sondages, pondérer et agréger. Par exemple, un sondage local (par état) est très différent d’un sondage national. Et un sondage sur toute la population n’est pas la même chose qu’un sondage sur des électeurs probables.
  2. Ajuster les sondages. En particulier, il peut être intéressant de lisser les résultats pour un même institut de sondage, en faisant des moyennes sur les sondages précédants. La subtilité est de faire varier le paramètre de lissage, avec des poids très conservateurs très en amont de l’élection (on lisse beaucoup, avec peu de variabilité entre deux dates) et au contraire, des poids plus faibles pour le passé au fur et à mesure que l’on se rapproche du jour de l’élection (on lisse peu, et les résultats deviennent très volatiles, très incertains)
  3. Combiner avec des données démographiques et économiques. Il est important de tenir compte de variables démographiques, aux Etats Unis, pour comprendre les changements possibles entre les élections de 2016 et celles de 2012. Par exemple en tenant compte de concentration éthniques, ou de changement d’âge moyen dans une zone géographique. Ensuite, il est possible de tenir compte de variables économiques, en particulier des niveaux d’emploi dans la région, le revenu personnel, des niveaux de consommation (personnels), ainsi que des niveaux d’indices boursiers (comme le note 538). On peut aussi tenir compte d’une information importante qui est le nombre de jours avant l’élection. Lauderdale & Linzer (2014) donne pas mal d’information sur les variables prédictives, justement.
  4. Projeter, par simulations. En simulant des scénarios, de manière corrélée entre les états, on peut arriver à des scénarios que l’on peut ensuite agréger, pour avoir des résultats nationaux.

Ce genre de modèle s’est relativement bien comporté, dans une marge d’erreur raisonnable. Mais le fait est que quelques pourcentages d’erreur ont donné des résultats très différents de ceux “attendus”. Ce qui est intéressant, c’est qu’avec 48% des voix en 2012, Mitt Romney a perdu l’élection présidentielle. Alors que Donald Trump l’a emporté, avec lui aussi près de 48% des voix. Mais Donald Trump a eu un résultat plus important que celui annoncé par les sondages dans des états clés (swing states). Un autre point important a été le fort changement de vote, jusqu’au bout. Plus de 15% des gens ont admis avoir décidé du candidat pour lequel ils voteraient au cours de la dernière semaine. D’autres sondeurs ont noté des fortes évolutions dans les non-réponses: beaucoup de personnes qui ont finalement voté pour Donald Trump ont longtemps refusé de l’avouer, et figuraient comme “non réponse”. Enfin, et c’est je pense un point important (j’avais tenté un billet sur le sujet voilà quelques années), les élections américaines ne sont pas des élections à deux candidats. Et mal prendre en compte Gary Johnson a probablement induit des erreurs dans les sondages. S’il était connu que Gary Johnson et Donald Trump chassaient le même électorat, beaucoup ont surestimé son score (estimé à 5% au niveau national, mais finalement à 3%), la différence s’étant probablement reporté sur Donald Trump

Après, il faut constater qu’il est difficile de construire un modèle, basé sur la rationalité des gens pour prévoir un comportement électoral. Quand 538 propose d’intégrer un indice économique, c’est parce qu’une hausse soudaine du chômage devrait avoir un impact sur le vote de certaines personnes (reste à savoir lequel, mais le choix des variables explicatives est un exercice en soi, comme le montre l’explosion des modèles sparse, permettant de trouver les bonnes variables prédictives, et non pas de rajouter davantage de bruit). Mais comment comprendre que Donald Trump, candidat ayant obtenu l’approbation du Ku Klux Klan, ait obtenu plus de voix que Mitt Romney en 2012 auprès des électeurs noirs, hispaniques et asiatiques ? Peut-on encore penser que les électeurs sont rationnels ?

Déclarer ou révéler, retour sur les sondages

Voilà une semaine, je mettais en ligne un court billet “les sondages et le bruit” sur les sondages, et leur biais, en préparation pour une intervention en milieu de semaine, sur Paris. Mais cette semaine, je voulais continuer à assembler mes notes de lecture…

  • Déclarer ou révéler

Bon, on va commencer par enfoncer des portes ouvertes, mais il est toujours bon de relire “l’opinion publique n’existe pas” de Bourdieu.

Toute enquête d’opinion suppose que tout le monde peut avoir une opinion ; ou, autrement dit, que la production d’une opinion est à la portée de tous. Quitte à heurter un sentiment naïvement démocratique, je contesterai ce premier postulat. Deuxième postulat : on suppose que toutes les opinions se valent. Je pense que l’on peut démontrer qu’il n’en est rien et que le fait de cumuler des opinions qui n’ont pas du tout la même force réelle conduit à produire des artefacts dépourvus de sens. Troisième postulat implicite : dans le simple fait de poser la même question à tout le monde se trouve impliquée l’hypothèse qu’il y a un consensus sur les problèmes, autrement dit qu’il y a un accord sur les questions qui méritent d’être posées. Ces trois postulats impliquent, me semble-t-il, toute une série de distorsions qui s’observent lors même que toutes les conditions de la rigueur méthodologique sont remplies dans la recollection et l’analyse des données.

Cela dit, j’irais probablement plus loin, en disant que le sondage crée une opinion, au sens où la personne interrogée se force à chercher la réponse à une question qu’elle ne s’était peut être jamais posée. Ou pas encore posée, dans le cas des sondages électoraux. Mais c’est le principe même des sondages, les personnes interrogée devant déclarer pour qui elles comptent voter le jour de l’élection (avec toutes sortes d’autres questions).

D’ailleurs, un des problèmes de cette “déclaration” est le biais de “désirabilité sociale“: les personnes interrogées aimer se présenter sous un jour favorable à l’enquêteur. Le vote pour Trump a été comme le vote pour le Front National (historiquement), quelquechose d’inavouable. Le soucis, c’est qu’il est difficile d’estimer ce biais. En France, une correction était possible, en tenant compte de la différence entre les sondages, et les résultats sur des élections passées. Mais aux Etats-Unis, il n’y avait pas de précédant.

Et si on va plus loin, on peut se demander s’il est pertinent d’utiliser des enquêtes sur les intentions de votes comme outils de prévision électorale. Ne pourrait-on pas y voir pour les personnes interrogées, de “mini-élection”, avec moins d’impact? Ne peut-on pas imaginer l’électeur socialiste, qui se dirait “et si je disais que je vais voter pour Jean-Luc Mélanchon, ça permettra peut être de tirer un peu plus vers la gauche, même si le jour J je voterais probablement pour Benoît Hamon” (le fameux vote “utile”). Compte tenu de l’importance des sondages pour les journalistes politiques, ne peut-on imaginer que les personnes interrogées jouent un jeu un peu plus subtile que juste réponse à une question sur leurs intentions de vote ? C’est un biais bien connu en économie expérimentale, lorsque l’on cherche à connaître les préférences des agents, entre les préférences déclarées et celles qui sont révelées.

  • Les sondeurs sont ils (vraiment) si mauvais ?

Les erreurs sont souvent bien connues. On peut penser aux élections en Angleterre de juin 1970, lorsque Harold Wilson perdit face à Edward Heathà la surprise générale (la légende prétend qu’il avait même vendu sa maison, tellement il était certain de rester au 10 Downing Street). En France, le choc d’avril 2002 reste dans la mémoire collective. Comment les instituts de sondages n’avaient pas vu venir le second tour entre Jacques Chirac et Jean-Marie Le Pen? C’est le repproche qui avait été fait, à l’époque, où les sondeurs ont été tenu pour responsable. Etrangement, retrospectivement, le risque était présent. Mais personne n’ausait l’admettre. Plus récement, en juin 2016, c’était le choc du Brexit, et en décembre 2016, c’était l’élection de Donald Trump. Si beaucoup on critiqué les sondages suite à son élection, il convient de relire certains articles publiés à cette époque. Par exemple, mi-octobre, “pourquoi des sondages donnent encore Trump gagnant ?“, parmi tant d’autres (oui, nous avons une mémoire très sélective).

Ce qui semble être reproché aux instituts de sondages, c’est leur opacité, et leur manipulation des chiffres. Mais il faut leur donner raison. Comme le dit la loi, “les listes électorales peuvent être consultées par tout électeur (même inscrit sur une autre commune) à condition de ne pas les utiliser à des fins commerciales“. Autrement dit, les instituts de sondages ne peuvent pas les utiliser pour constituer un plan de sondage. Ne pouvant faire de vrai sondage aléatoire, il faut donc ruser, et on va échantillonner par quotas. C’est ce qu’on appelle du redressement. Pour faire simple, on applique des poids pour se rapprocher de la “population de base”. Qu’on ne connait pas. Idéalement, on aimerait bien connaître ces gens inscrits sur les listes. Mieux encore, connaître ceux qui iront effectivement voter (car ce sont vraiment ces personnes qui nous intéressent). Si aucun homme ne va voter, le jour J, il conviendrait de faire un échantillon de sondage avec uniquement des femmes. Encore une fois, comprendre l’abstention est essentiel.

Pour faire les quotas, on utilise des données géographiques, mais aussi d’âge, de genre, et de catégorie soci-professionnelle. En fait, cette dernière est (malheuresement) seulement une catégorie profesionnelle. On pensera à cette merveilleuse catégorie des “inactifs” qui correspond à des réalités très différentes. On utilise très peu les données de niveau d’étude par exemple, comme le rappelle Avner dans les sondages sont-ils devenus fous, les personnes faiblement diplômées sont faiblement présentes dans les sondages,

Il devient nécessaire de faire des redressements. Mais ces derniers n’ont de sens que s’il bien faits. Certains préconisent d’ailleurs des catégories “socio-affinitaires“, et non plus socio-professionnelles.

  • On veut des chiffres !

Allez, quelques chiffres, que j’ai découvert dans l’article les sondages sont-ils devenus fous d’Avner Bar-Hen, qui permettent d’avoir un ordre de grandeur intéressant. Pour faire une enquête, on peut avoir envie d’avoir 5000 personnes qui répondent à une liste de questions. Mais combien de personnes faut-il solliciter, par téléphone, pour avoir 5000 réponses. La réponse est simple, plus de 15 fois plus : pour avoir 5240 réponses, il a été nécessaire de composer 84000 coups de fils,

  • Des intensions de vote à une probabilité de gagner

Les sondages reposent sur la théorie de l’échantillonnsages, et il y a une part (importante) d’incertitude, comme on le voyait la dernière fois. Et comme le savent les statisticiens, il est possible de mal faire son travail, mais d’avoir de la chance, et de faire une bonne prévision, mais aussi de faire son travail correctement, et de se tromper lourdement.

Ce qui m’a intéressé lors des dernières élections américaines, c’est qu’on ne regarde pas le nombre de voix (prédit), mais la probabilité de gagner, pour tel ou telle candidat. Ce qui a du sens quand on pense au mode de scrutin. Mais traduit dans un contexte plus français (où on se place entre les deux tours), la probabilité de gagner est une information sur l’incertitude du sondage.

> simu=function(i){
+ n=1+2000
+ vote=sample(c("A","B"),size=n,
+ replace=TRUE,prob=c(.515,.485))
+ france1=mean(vote=="A")
+ n=1+200
+ vote=sample(c("A","B"),size=n,
+ replace=TRUE,prob=c(.515,.485))
+ france2=mean(vote=="A")
+ prop=c(france1,france2)
+ names(prop)=c("FR1","FR2")
+ return(prop)
+ }
> Sbi=Vectorize(simu)(1:1e4)
> plot(density(Sbi[1,]))

Ici, le candidat rouge a 15% de chances de gagner, contre 85% pour le candidat bleu. Avec 2000 personnes interrogées. Mais si on a seulement 200 personnes interrogées, on passe à 40% de chances de gagner pour A.

> plot(density(Sbi[2,]))

Autrement dit, la probabilité de gagner est, à mon sens, essentiellement une information liée à l’incertitude qui entoure le scrutin.

  • Quid du… “big data” ?

Allez, mettons enfin les pieds dans le plat… pourquoi pas utiliser “le big data”? Oui, le “big data” peut tout, je pense qu’il peut même assurer le retour de l’être aimé, alors il doit bien pouvoir prédire qui va gagner une élection…

Il y a 15 ans, internet n’était pas considéré comme un média. Personne n’osait trop faire des sondages en ligne*. Cette époque est révolue*. Au point que bon nombre se sentent aujourd’hui observés, et adoptent des stratégies, des comportements, qu’il convient de décoder en prenant des précautions.

L’écoute de la toile permet d’éviter l’écueil du déclaratif, pour observer, et induire des tendances, des comportements. Il serait intéressant de voir ce qu’en pensent les data scientists de chez Facebook par exemple. Comme le notait Stéphane Rozès, “si la France est une grande consommatrice de sondages, l’idée même qu’un statisticien puisse annoncer en amont le résultat de l’élection est baroque et attentatoire à l’imaginaire politique français.C’était en 2013. Depuis, comme le montre the pulse of the people, l’analyse des tweets (par exemple) permet d’observer beaucoup de choses. La difficulté est que si une équipe électorale sait ce qui est observé, il est aisé de créer des robots pour manipuler l’analyse. Le jeu du chat et de la souris peut commencer.

Mais il n’en reste pas moins que l’analyse des traces laisses sur la toile est importante, et les instituts de sondages le savent. Un “like” est une information intéressante, de même que savoir les pages vues, et les liens partagés. Cela dit, si on pousse un peu le raisonnement, comme le faisait Andrew Gelman dans un article paru au début du mois dans Science, si on était capable de prédire parfaitement le résultat d’une élection, cette dernière devient superflue, non ?

* ce qui est intéressant, c’est que dans le cas du Brexit, ce sont justement les sondages en ligne qui avaient “prévu” la victoire du leave, contrairement aux sondages téléphoniques qui donnaient le remain en tête….

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.

Prédire ou Prévoir, réflexions (probablement) personnelles

Il y a quelques années, j’avais fait un billet où je m’énervais contre une phrase que j’avais entendu (“l’homme qui avait prédit la crise). Le point de départ était surtout un énervement lié au fait que beaucoup de journalistes vivent avec des oeillères (par exemple en associant “économiste” à “type à la télé qui parle d’économie“), et ne connaissent pas la diversité de la communauté des économistes (surtout que la pluralité n’est pas vraiment une vertue présente dans les grands médias, la manière dont sont présentés les différents candidats à la future élection présidentielle peut en être une illustration). Mais je m’écarte…. en fait, un de mes points portait sur l’utilisation du terme “prédit“.

Avec un peu de recul, je ne sais pas si la langue française à changé depuis, mais j’aurais envie de revenir sur ma distinction entre “prévoir” et “prédire“. Prévoir et prédire ont deux racines latines distinctes, prévoir venant du latin praevidere, et prédire de praedicere. Beaucoup notent que “prédire” a longtemps eu le sens de “prophétiser“, où l’annonce (de ce qui doit arriver) est souvent d’inspiration divine. En anglais, c’est presque plus compliqué (car il y a plus de synonymes) mais la même confusion existe. On peut retrouver “foretell” construit sur “fore” et “tell” (qui est à relier à une notion de prémonition), “forespell” construit sur “fore” et “spell” (on retroue la notion de sort, proche du néerlandais “voorspellen”), “foresee” construit sur “fore” et “see”, ou encore  “forecast“, dont on reparlera très bientôt. A côté de ces mots, construits sur le préfixe “fore”, il existe le mot “predict“, tiré du latin praedicere, construit sur “prae” ‎(avant, que l’on retrouvé dans “pré” en français) et “dicere” ‎(dire). Comme le note le dictionnaire, “predicted” a un sens proche de “foretold” mais “without the mystical flavour”. Bref, “prédit” a longtemps eu une connotation mystique, divine.

Pourtant, les deux semblent parfois être interchangeables. Au point que fin 2006, Le Devoir écrivait « parce que cette administration n’a pas su prévoir ce qui était prévisible, nous en sommes rendus à siéger à 6 h un lundi matin! ». “prévoir ce qui était prévisible” vaste programme, non ?

Si on reste sur l’anglais, on peut tenter une distinction. Par exemple “forecast“, c’est à dire “prévoir“, est associé aux ouvrages suivants, sur amazon.com

Autrement dit, “forecast” (ou “prévision“) correspond à un aspect temporel. On cherche à prévoir, à un instant donné https://latex.codecogs.com/gif.latex?t la valeur d’une variable https://latex.codecogs.com/gif.latex?x à une date ultérieur, disons https://latex.codecogs.com/gif.latex?t+h. Classiquement, on cherche (et ça sera notre “prévision”)

https://latex.codecogs.com/gif.latex?\mathbb{E}_t[X_{t+h}]=\mathbb{E}[X_{t+h}|x_t,x_{t-1},\cdots]

C’est pour cela que bon nombre de cours de séries temporelles, ou séries chronologiques sont des cours de “modèles de prévision“.

Mais on commence à voir arriver des ouvrages de “méthodes prédictives“, “predictive modeling“, proposant des “prédictions

Si on regarde de plus près, ces techniques sont utilisées non pas dans un contexte de séries temporelles, mais de données individuelles. Ici, on cherche à prévoir une variable d’intérêt https://latex.codecogs.com/gif.latex?y à l’aide de covariables (variables dites explicatives https://latex.codecogs.com/gif.latex?x_1,\cdots,x_k. La prédiction sera classiquement

https://latex.codecogs.com/gif.latex?\mathbb{E}_{\boldsymbol{x}}[Y]=\mathbb{E}[Y|x_1,\cdots,x_k]

Sur le graphique du taux d’imposition, en France

Voilà plusieurs années (en 2011) Landais, Piketty et Saez publiaient le graphique suivant, présentant le sytème fiscal, en France, avec l’évolution du taux global d’imposition, en fonction du revenu.

https://freakonometrics.hypotheses.org/files/2017/02/revol-fiscale.png

Il y a deux ans, Yves Guiard critiquait ce graphique (rappelant que “la conversion numérique/graphique obéit à certaines règles“) et proposait une version “correcte“, sous la forme

J’ai déjà dit – voilà plusieurs semaines – que cette vision dogmatique des règles de représentation graphique ne me plaisait pas. Mais le plus simple est peut être de reprendre la réponse de Thomas Piketty, sur cette critique (je me retrouve totalement dans le propos),

(je suis désolé de découvrir cette discussion aussi tard, Béatrice m’a posté un commentaire ce matin, qui pointait vers ces deux articles)

Mais essayons de regarder un peu plus en détails. Oui, si on utilise la représentation standard de la courbe, on retrouve le graphique suivant

avec les plus pauvres à gauche, et les plus riches à droites. On a ici un graphique où le même pas sur l’axe des abscisses correspond au même nombre d’individus. Les 10% les plus pauvres (bande verticale la plus à gauche) ont un taux moyen d’imposition de 40.46% alors qu’il est de 48.09% pour les 10% qui sont juste au dessus du salaire médian (bande juste à droite de 50%). Je ne reviendrais pas sur le débat de commencer le graphique à 0% (ce débat aurait du sens, j’en conviens, mais commencer à 0% et mettre du texte entre l’axe des abscisses et la courbe, c’est un peu tricher), mais si on respecte cette règle, on obtient

Mais ce n’est pas cette visualisation qui est retenue dans le livre de Landais, Piketty et Saez, mais plutôt quelque chose qui ressemble au graphique ci-dessous (je commence à 30% pour éviter d’avoir un billet trop long)

Ce graphique est (je pense) plus intéressant que le premier. Bon, bien sûr, il se rapproche de la courbe parfaite que tout statisticien qui se respecte essaye reproduire

mais surtout, il permet de tenir compte de l’hétérogénéité qui existe dans ce dernier décile. Histoire d’illustrer mon propos, je voudrais proposer une représentation alternative de celle de Landais, Piketty et Saez. Pour se faire, on va utiliser une information complémentaire (en ligne sur le site) à savoir le revenu moyen, par centile, dans la population (voire plus fin dans le centile le plus élevé)

L’idée sera d’utiliser en abscisse non pas un pourcentage de personnes, mais un pourcentage de richesse (un peu comme avec la courbe de Lorenz),

Par exemple sur le graphique ci dessus, la bande la plus à gauche ne représente pas 10% des français, mais la fraction de français les plus pauvres, qui totalisent 10% de la richesse totale. La bande la plus à droite, ce sont les français les plus riches, qui accumulent 10% de la richesse totale. On peut d’ailleurs être un peu plus éclairant (je pense) en mettant des bandes non pas tous les 10% du revenu cumulé, mais faire l’équivalent du premier graphique, avec 10% des français. Bien entendu, les 10% les plus pauvres possède un très faible pourcentage de la richesse totale. En fait, les 30% les plus pauvres possèdent, tous ensemble, 5% de la richesse totale. Alors que les 10% les plus riches possèdent un tiers de la richesse. C’est ce qu’on peut voir sur le graphique ci-dessous

Ce graphique est – je pense – juste, au sens des normes sémiologiques, et permet de faire passer le message du graphique initial.

Après, on peut tenter une vision qui irait encore plus loin, avec non plus un taux en abscisse, mais un niveau de richesse. Aussi, le taux d’imposition, en fonction du revenu brut mensuel (sur une échelle logarithmique) ressemble au graphique ci-dessous

Forecasting Natural Catastrophes (is rather difficult)

Following my previous post, I wanted to spend more time, on the time series with “global weather-related disaster losses as a proportion of global GDP” over the time period 1990-2016 that Roger Pilke sent me last night.

db=data.frame(year=1990:2016,
ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))

In my previous post, I spend some time explaining that we should provide some sort of ‘confidence interval’ when we try to predict a pattern. That was what we call ‘model uncertainty’. But there are two (important) issues that I did not mention. (1) it is a time series, so why not use techniques dedicated to time series objects ? (2) we do not really care actually about ‘model uncertainty’ (unless we want to assess if a decreasing trend is significant, or not), and we care more about real prediction uncertainty: in the next ten years, what could be the range for the this ratio, with some given probability (say 95%)? Could we say that with 95% chance the global weather-related disaster losses as a proportion of global GDP should be (each year) within 0 and 0.35 or 0 and 0.7?

A first idea might be to use exponential smoothing techniques (without a seasonal component here).

ratio=ts(db$ratio,start=1990,frequency=1)
plot(ratio,xlim=c(1990,2030))
hw=HoltWinters(ratio,gamma=FALSE)
phw=predict(hw,n.ahead=15,prediction.interval = TRUE)
plot(hw,phw,xlim=c(1990,2030))
polygon(c(2017:2031,rev(2017:2031)), c(phw[,2],rev(phw[,3])),border=NA,col=rgb(0,0,1,.2))

The decreasing trend is coming from the fact that exponential smoothing is here a linear regression, with weight exponentially decaying with time (the older, the smaller the weight). But we cannot use that prediction, since the ratio cannot (obviously) be negative. So why not consider, here, the logarithm of the ratio

plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
hw=HoltWinters(log(ratio),gamma=FALSE)
phw=predict(hw,n.ahead=15,prediction.interval = TRUE)
abline(v=2016,lty=2,col="grey")
lines(2017:2031,exp(phw[,2]),col="blue")
lines(2017:2031,exp(phw[,3]),col="blue")
lines(c(1992:2016,2017:2031),c(exp(hw$fitted[,1]),exp(phw[,1])),col="red")
polygon(c(2017:2031,rev(2017:2031)),exp(c(phw[,2],rev(phw[,3]))),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

The confidence band is huge, here. What if we consider some ARIMA model here?

fit=auto.arima(ratio)
farma=forecast(fit,15)
farma=cbind(as.numeric(farma$fitted)[1:15],as.numeric(farma$lower[,1]),as.numeric(farma$upper[,1]),as.numeric(farma$lower[,2]),as.numeric(farma$upper[,2]))
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
lines(2017:2031,farma[,4],col="blue")
lines(2017:2031,farma[,5],col="blue")
lines(2017:2031,farma[,1],col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,4],rev(farma[,5])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Here, there is an intercept, but no dynamics for the time series (which is considered, here, as a pure white noise). We get exactly the same if we consider the average value of the series

fit=lm(ratio~1,data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Here, we get back to my previous post, if we want to consider a possible trend (and not only an intercept)

fit=lm(ratio~year,data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Again, the confidence region is not based on inference related error, but on model uncertainty: we try to visualize where future observations might be with (say) 95% chance. Note we can also consider (why not?) a quadratic regression

fit=lm(ratio~poly(year,2),data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

I am usually not a huge fan of those polynomial regression, but recently, I’ve seen that a lot in economic papers (like “if it’s not linear, add a squared version of the explanatory variable”, which is a rather odd strategy, I’ll publish some posts on that issue this year).

Here again, it might be more clever to consider a logarithmic transformation of the ratio, to insure that the ratio remains positive

fit=lm(log(ratio)~year,data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(exp(pf+s^2/2),exp(pf-1.96*s),exp(pf+1.96*s))
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(exp(predict(fit)+s^2/2),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Observe that future trend is mainly driven by the three latest observations, that were rather low (compared with older observations). What if we remove them?

dbna=db
db$ratio[25:27]=NA
fit=lm(ratio~1,data=dbna)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016-3,lty=2,col="grey")
ndb=data.frame(year=2014:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2014:2031,farma[,2],col="blue")
lines(2014:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red")
polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

More funny, if we consider a quadratic regression, we obtain an increasing trend for the future

fit=lm(ratio~poly(year,2),data=dbna)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016-3,lty=2,col="grey")
ndb=data.frame(year=2014:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2014:2031,farma[,2],col="blue")
lines(2014:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red")
polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

As we can see, it is rather difficult to get relevant prediction for the future, based on 25 observations…. If anyone has a suggestion, comments are open…

 

What is a Linear Trend, by the way?

I had a very strange discussion on twitter (yes, another one), about regression curves. I think it started with a tweet based on some xkcd picture (just for fun, because it was New Year’s Day)

There were comments on that picture, by econometricians, mainly about ‘significant’ trends when datasets are very noisy. And I mentioned a graph that I saw earlier, a couple of days ago

Let us reproduce that graph (Roger kindly sent me the dataset)

db=data.frame(year=1990:2016,
ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))
library(ggplot2)

The graph is here (with the same aesthetic conventions as Roger’s initial graph, i.e. using some sort of barplot)

ggplot(db, aes(year, ratio)) +
geom_bar(stat="identity") +
stat_smooth(method = "lm", se = FALSE)

My point was that we miss the ‘confidence band’ of the regression

In R, at least, it is quite natural to get (and actually, it is the default version of the graph function)

ggplot(db, aes(year, ratio)) +
geom_bar(stat="identity") +
stat_smooth(method = "lm", se = TRUE)

It is hard to claim that the ‘regression line’ is significant (in the sense “significantly non horizontal”). To be more specific, if we look at the output of the regression model, we get

summary(lm(ratio~year,data=db))

Coefficients:
Estimate    Std. Error t value Pr(>|t|)
(Intercept) 9.158531 4.549672  2.013 0.055 .
year       -0.004457 0.002271 -1.962 0.061 .
---
Signif. codes: 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(which is exactly what Roger used in his graph to plot his red straight line). The p-value of the estimator of the slope, in a linear regression model is here 6%. But I found Roger’s point puzzling

See also

First of all, let us get back to a more standard graph, with a scatterplot, and not bars,

ggplot(db, aes(year, ratio)) +
stat_smooth(method = "lm") +
geom_point()

Here, we observe points https://latex.codecogs.com/gif.latex?\{y_{1990},y_{1991},\cdots,y_{2016}\}. In order to draw that blue line, we assume (Econometrics 101, actually) that those observations are realizations of random variables https://latex.codecogs.com/gif.latex?\{Y_{1990},Y_{1991},\cdots,Y_{2016}\}. Randomness here does not come from a survey, or from ‘balls in an urn’. Randomness is because hurricanes and floods are themselves seen are realizations of random events. Yes, there might be measurement errors, but that’s not where randomness comes from (here). When we talk about ‘randomness’, it should be related to ‘model error’ i.e. the error we make if we consider a linear model (here), that is

https://latex.codecogs.com/gif.latex?Y_t=\beta_0+\beta_1t+\varepsilon_t

Even if observations are not obtained from balls in an urn, there is some kind of randomness here. Randomness means that we might have errors (random errors) around the estimated value (that is on the blue curve), https://latex.codecogs.com/gif.latex?y_t=\widehat{y}_t+\widehat{\varepsilon}_t. One might consider a nonlinear model to reduce the error,

ggplot(db, aes(year, ratio)) +
geom_point() +
geom_smooth()

but in the case, the danger is to overfit

So yes, when we fit a linear model, there is always some kind of randomness, and it is possible to get a ‘confidence band’, that will be very useful for predictions (e.g. for reinsurance purpose here).

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 ?

Les gens bien élevés font commencer les barres en zéro

Il y a quelques jours, je mentionnais sur twitter un graphique mis en ligne sur imgur.com, au sujet des trois dernières élections présidentielles américaines,

Le premier point est que le graphique a été fait alors que les dépouillements n’étaient pas encore finis. 72 heures après la clôture des bureaux de votes, plusieurs millions de bulletins n’avaient pas été dépouillés, comme me l’ont fait noter François Gourio et Tom Roud. Comme on visualise ici le nombre de votants, et pas le pourcentage de votants, c’est génant… Le second point avait été souligné par Marin Dacos, et par plusieurs personnes, toujours sur twitter, qui disaient (je résume un peu)

On ne peut pas faire des diagrammes en barres qui ne commencent pas en zéro ! Sacrilège !

C’est en effet une règle qui semble bien établie. Enfin, quand je dis ‘règle‘, on pourrait penser qu’il s’agit d’une ‘suggestion‘. Des sortes de ‘règles de savoir vivre’, à la Nadine de Rothschild,

Quand on est bien élevé, il est de bon ton de faire commencer ses diagrammes en barres en zéro.

Et cette règle semble importante, car c’est la première que donne Nathan Yau dans ses règles pour faire des graphiques

A bar chart’s baseline must start at zero.

L’intuition est que dans un diagramme en barres, la longueur des barres a un sens.

The bar chart relies on length to show data. Shorter bars represent lower values, and longer bars represent greater values. Compare bar lengths to compare values. That’s how it works. When you shift the baseline, you distort the visual.

Quittons nos élections un instant. Regardons par exemple les temps mis au 100 mètres lors des Championnats du Monde d’Athlétisme.

> year=c(1983,1987,seq(1991,2015,by=2))
> time=c(10.07,9.93,9.86,9.87,9.97,9.86,9.80,9.82, 10.07,9.88,9.85,9.58,9.92,9.77,9.79)
> names(time)=year
> barplot(time,col="darkgreen",beside=TRUE)

Oui, on a compris, les barres ont une interpréation physique, et il ne faut pas les raccourcir ! C’est le temps mis par un sportif pour courrir 100 mètres. Mais force est de constater qu’on ne voit pas grand chose… Allez, souvenez vous en 2009,

Sur notre diagramme en barre, on ne voit rien d’exceptionnel… Allez, osons faire commencer l’axe des ordonnées à 9 secondes,

> barplot(time,ylim=c(9,10.2),col="darkgreen")

Ok, je suis un hérétique… Certes, je mesure ici le temps au delà de 9 secondes, ce qui ne veut pas dire grand chose… J’opère ici la “distortion” dont parle Nathan Yau…. mais j’ai quand même l’impression qu’on voit mieux l’exploit de Bolt (sans laquelle on ne voit absolument rien). Cela dit, on pourrait même tenter quelque chose d’encore plus violent

> barplot(time,ylim=c(min(time),10.2), col="darkgreen")

C’est ridicule comme graphique n’est-ce pas ? Cela dit, à peu de chose près, ce graphique peut se lire en translatant, et en disant qu’il s’agit le temps mis en plus, par rapport au record de Bolt en 2009,

et ce graphique, lui, a du sens et suit la règle de Nathan Yau : il commence en 0.

Sours R, naturellement, un diagramme en baton commencera en 0  mais ça ne sera pas le cas pour une courbe,

Ici, le graphique commence à peu près au minimum de nos observations. Et on comprend avec ce graphique le caractère exceptionnel du temps mis par Bolt, en 2009. Faire commencer le graphique en 0 est ridicule,

C’est presque pire lorsqu’on représente l’évolution du record du monde, au 100 mètres,

Ca n’a pas d’allure ! On ne voit rien sur ce graphique. Alors qu’en changeant l’axe des ordonnées, ça devient plus clair !

Bref, je suis surpris qu’on passe pour un hérétique quand on dit que ce n’est pas génant de ne pas commencer en 0. Ca dépend du message que l’on souhaite faire passer avec le graphique.

En fait, j’ose affirmer que commencer en 0 n’a pas sens ici, car il est impossible qu’un être humain court le 100 mètres en 0 secondes. Cette référence n’a pas lieu d’être ! Mais supposons qu’on nous dise que pour des raisons morphologiques, on pourrait imaginer qu’un homme court, un jour, le 100 mètres en 9 secondes, mais qu’il serait physiologiquement impossible de passer en dessous. Dans ce cas, j’ai envie de croire que la réprésentation en commençant à 9 secondes aurait du sens…

On voit ici non seulement la différence entre Bolt en 2009 et les autres champions du monde, mais on voit aussi la distance qui sépare ces athlètes d’un record absolu. Je suis désolé, mais je ne trouve pas ce graphique dérangeant. Certes, il convient de lire (rapidement) l’axe des ordonnées, pour éviter d’affirmer que Bolt aurait mis 2 fois moins de temps en 2009 que le champion du monde de 2003 pour courrir le 100 mètres (car la barre semble deux fois plus haute).

Revenons maintenant à nos élections américaines, en reprenant les nombres de votes obtenus par candidats en 2008, 2012 et 2016 (avec des chiffres mis à jour)

> votes=data.frame(red=c(59948323,60933504, 62213384),blue=c(69498516,65915795,64154610))
> MV=as.matrix(t(votes))
> colnames(MV)=c(2008,2012,2016)
> barplot(MV/1e6, main="",
+         xlab="", col=c("red","darkblue"),
+         beside=TRUE)

Là encore, j’aurais presque envie d’affirmer que commencer en 0 n’a pas grand sens. Tout d’abord parce qu’il est dur d’imaginer que le candidat des Républicains comme le candidat des Démocrates ait un score proche de zéro. Ensuite, il est aussi difficile d’imaginer qu’un candidat ait deux fois plus de vote que l’autre. Car on peut comparer des barres quand une des deux fois plus grande que l’autre, pas quand elle est 14% plus grande… ça ne se voit pas à l’oeil nu… disons que la différence entre 14% et 18% n’est pas visible. Aussi, quand on connait le contexte, on sait qu’un candidat ne peut pas avoir deux fois plus que l’autre… Donc quand on voit un diagramme à baton de ce genre, on va regarder les valeurs des ordonnées.

> barplot(MV/1e6, main="",
+         xlab="", col=c("red","darkblue"),
+         beside=TRUE,ylim=c(52,72))

J’ai envie de dire qu’on pourrait accepter de fausser la représentation, parce que ce n’est pas ce qu’on veut illustrer avec ce graphique. Ce qu’on veut montrer c’est que le nombre de votes rouges varie augmente un peu, et que le nombre de votes bleus a baissé. Avec des courbes (là on a le droit), on aurait la représentation suivantes

mais je trouve le diagramme en barres plus joli. Et je trouve dommage de dire que le diagramme avec les barres juste au dessus est faussé, alors que le diagramme avec les courbes ne l’est pas… Si on fait commencer en zéro, comme le dit la règle, on ne voit plus rien. Sauf qui ça se joue dans un mouchoir de poche,

Là encore, si on veut respecter la règle du zero, et faire malgré tout des barres, peut être faut-il représenter les variations du nombre de votants (en millions), d’une élection sur l’autre

mais je n’aime pas trop…. Là encore, le graphique le plus parlant reste celui mentionné dans mon tweet

(la nuance est peut être que je ne mets pas ici de barre horizontale en bas). On peut voir le message qu’on veut faire passer avec le graphique de manière beaucoup plus claire. Certains diraient que je fausse l’interprétation, mais je trouve l’affirmation excessive. Il s’agit plutôt de déformer pour s’assurer justement qu’un message (juste et vrai) passe mieux. C’est ce qu’on voyait avec l’exemple du 100 mètres : sans déformation, le seul message que l’on retient est que le record de Bolt est anecdotique.

J’aurais presque envie de croire que la règle ne s’applique pas si la valeur zéro n’est pas atteignable (et n’aurait pas de sens). Mais je ne demande pas mieux que d’être convaincu que même ici, la règle de commencer les barres en zéro doit s’appliquer !

Dominance stochastique, ordre 1 ou ordre 2 ?

Mercredi dernier, on commençait le cours avec une présentation rapide du test de Wilcoxon-Mann-Whitney (parfois appelé test de la somme des rangs, de Wilcoxon). Ce test a été proposé par Frank Wilcoxon en 1945, et (presque en même temps) par Henry Mann et Donald Whitney (en 1947), dans un article intitulé On a Test of Whether one of Two Random Variables is Stochastically Larger than the Other. Le titre est clair. Ce test (non paramétrique) a vocation à tester l’hypothèse selon laquelle la distribution dans deux échantillons est la même. Ou, pour reprendre les notations du derniers billets, la distribution est la même dans deux groupes.

Autrement dit, on va chercher à tester, ici, https://latex.codecogs.com/gif.latex?H_0:F=G, au sens où https://latex.codecogs.com/gif.latex?H_0:\forall%20x,F(x)=G(x). Dans le test de Kolmogorov Smirnov, l’hypothèse alternative est https://latex.codecogs.com/gif.latex?H_1:\exists%20x,F(x)\neq%20G(x).

Ici, comme le notent d’ailleur Mann et Whitney, “the test is shown to be consistent with respect to the class of alternatives https://latex.codecogs.com/gif.latex?F(x)%3EG(x) for every https://latex.codecogs.com/gif.latex?x“. Cette propriété est un peu particulière, et correspond à la notion de dominance stochastique à l’ordre 1. https://latex.codecogs.com/gif.latex?F\preceq_1%20Gcorrespond à https://latex.codecogs.com/gif.latex?F(x)\geq%20G(x) https://latex.codecogs.com/gif.latex?\forall%20x, ou encore, pour tout fonction https://latex.codecogs.com/gif.latex?u(\cdot) croissante, si https://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)] avec https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G.

Mais regardons un peu le test. On considère deux populations https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_{n_x}\} et https://latex.codecogs.com/gif.latex?\boldsymbol{y}=\{y_1,\cdots,y_{n_y}\}. L’idée est de réordonner les https://latex.codecogs.com/gif.latex?(n_x%20+%20n_y) éléments de https://latex.codecogs.com/gif.latex?\boldsymbol{z}=\boldsymbol{x}\cup\boldsymbol{y} par ordre croissant. On note https://latex.codecogs.com/gif.latex?\boldsymbol{r} le vecteur des rangs associés. Soit https://latex.codecogs.com/gif.latex?S_x la somme des https://latex.codecogs.com/gif.latex?n_x (premiers) rangs des éléments de https://latex.codecogs.com/gif.latex?\boldsymbol{x}, i.e. https://latex.codecogs.com/gif.latex?S_x=r_1+\cdots+r_{n_x}.

Comme l’ont montré Wilcoxon, Mann et Whitney, sous https://latex.codecogs.com/gif.latex?H_0, https://latex.codecogs.com/gif.latex?S_x suit une “loi connue”, où classiquement, “connue” veut dire tabulée pour de petits échantillons, et pour des échantillons plus grands.

Histoire d’illustrer, considerons deux échantillons, tirés suivant deux lois https://latex.codecogs.com/gif.latex?F et https://latex.codecogs.com/gif.latex?G avec https://latex.codecogs.com/gif.latex?F\preceq_1%20G, par exemple deux lois normales de moyennes différentes (et de même variance)

> u=seq(-3,4,length=251)
> plot(u,pnorm(u,0,1),col="red",type="l")
> lines(u,pnorm(u,.5,1),col="darkblue")

Effectivement, nos deux fonctions de répartition sont ordonnées, et on est bien dans une situation de dominance stochastique à l’ordre 1. Les densité sont ici

>plot(u,dnorm(u,0,1),col="red",type="l")
>lines(u,dnorm(u,.5,1),col="darkblue")

On peut générer 10,000 paires d’échantillons, de taille 100 (respectivement), et regarder la distribution de https://latex.codecogs.com/gif.latex?S_x (et celle de https://latex.codecogs.com/gif.latex?S_y).

> P=matrix(NA,1e4,3)
> for(i in 1:1e4){
+  x=c(rnorm(100,0,1),rnorm(100,.5,1))
+  P[i,1:2]=c(sum(rank(x)[1:100]),sum(rank(x)[101:200]))
+  x0=c(rnorm(100,0,1),rnorm(100,0,1))
+  P[i,3]=  sum(rank(x0)[1:100])
+}
> plot(density(P[,1]),col="red"))
> lines(density(P[,2]),col="darkblue")
> lines(density(P[,3]),lty=2)

Sous https://latex.codecogs.com/gif.latex?H_0, la somme aurait eu la distribution avec les pointillés. On en est manifestement assez loin, autrement dit, on va rejeter ici l’hypothèse de lois égales, https://latex.codecogs.com/gif.latex?H_0:F=G.

Intuitivement, c’est normal, car effectivement, les rangs dans un cas sont – en moyenne – plus faibles que pour l’autre,

Au lieu de simuler des lois dominées, à l’ordre 1, on peut tenter des lois dominées à l’ordre 2. Rappelons que https://latex.codecogs.com/gif.latex?F\preceq_2%20G signifie

https://latex.codecogs.com/gif.latex?\int_{-\infty}^x%20F(t)dt%20\geq\int_{-\infty}^x%20G(t)dt,\forall%20x

ou encorehttps://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)]

avec https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G. où https://latex.codecogs.com/gif.latex?u(\cdot) est unne fonction croissante et concave.

En gestion des risques (ou sur les inégalités), on va supposer en plus que https://latex.codecogs.com/gif.latex?\mathbb{E}[X]=\mathbb{E}[Y], ce qui correspond à une notion de “mean preserving spread” (on peut relire le document sur la notion de
single crossing, par Alain Chateauneuf). Une caractérisation est

https://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)]https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G. où https://latex.codecogs.com/gif.latex?u(\cdot) est concave. Une caractérisation alternative est https://latex.codecogs.com/gif.latex?Y=X+\varepsilon avec https://latex.codecogs.com/gif.latex?\mathbb{E}[\varepsilon\vert%20X]=0 (on pourra relire Rothschild et Stiglitz, Increasing Risk: A Definition)

Dans le cas des lois normales, on a cette notion de ‘mean preserving spread’ lorsque les moyennes sont identiques, mais pas les variances. On a alors plus de dispersion,

> u=seq(-4,4,length=251)
> plot(u,pnorm(u,0,1),col="red",type="l")
> lines(u,pnorm(u,0,1.5),col="darkblue")

Effectivement, nos deux fonctions de répartition sont ordonnées, et on est bien dans une situation de dominance stochastique à l’ordre 2, avec du single crossing. Les densités sont ici

> plot(u,dnorm(u,0,1),col="red",type="l")
> lines(u,dnorm(u,0,1.5),col="darkblue")

On peut générer 10,000 paires d’échantillons, de taille 100 (respectivement), et regarder la distribution de https://latex.codecogs.com/gif.latex?S_x (et celle de https://latex.codecogs.com/gif.latex?S_y).

> P=matrix(NA,1e4,3)
> for(i in 1:1e4){
+  x=c(rnorm(100,0,1),rnorm(100,0,1.5))
+  P[i,1:2]=c(sum(rank(x)[1:100]), sum(rank(x)[101:200]))
+  x0=c(rnorm(100,0,1),rnorm(100,0,1))
+  P[i,3]=  sum(rank(x0)[1:100])
+}
> plot(density(P[,1]),col="red"))
> lines(density(P[,2]),col="darkblue")
> lines(density(P[,3]),lty=2)

Autrement dit, ici les sommes des rangs sont identiques, avec les deux échantillons,  et donc on aura tendance à accepter ici l’hypothèse de lois égales, https://latex.codecogs.com/gif.latex?H_0:F=G. L’idée naturelle est que avec les dans le premiers cas, effectivement, on avait plus de petits rangs avec le petit échantillon. Mais plus ici: le fait d’avoir plus de dispersion fait qu’on a en même temps plus de petits, et plus de grands. Et en moyenne, ça ne change rien

Modéliser et prévoir, entre ‘avoir faux’ et ‘avoir pas de chance’

Depuis plusieurs heures, il y a une phrase qui revient sans cesse, et qui m’agace au plus haut point, sur “le modèle est faux…“. C’est une phrase que j’ai entendu pas plus tard qu’au séminaire ce midi, lorsqu’un collègue a dit que très clairement les modèles étaient faux puisqu’ils n’avaient pas prédit le vainqueur des élections américaines. Par exemple, le jour des élections, Huffington Post annonçait que Donald Trump avait 1,6% de chances d’être élu…

et pourtant, il a gagné. C’est donc pourtant que le modèle était “faux”, c’est évident, non ?

Pourtant, tout ceux qui ont fait un minimum de statistique dans leurs vies le savent… les évènements rares, ça peut se produire,

> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "clinton"
> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "clinton"
> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "trump"

Ce type de raisonnement me dérange car il dénote une mauvaise compréhense de ce qu’est une probabilité. Autre exemple: supposons que j’annonce que ma belle mère a une chance sur un million de mourrir lorsqu’elle prendra la route pendant les prochaines vacances, et qu’elle meurt effectivement cet hiver d’un accident de voiture. Comme pour les élections américaines, on sera surpris, sous le choc. Mais est-ce que ça veut dire que mon modèle était faux ? Non…Si j’avais annoncé qu’elle avait 85% de ne pas passer les vacances en prenant la voiture, j’aurais eu un modèle des plus douteux… Le problème ici c’est qu’on avance des probabilités sur des évènements qu’on ne pourra observer qu’une fois. On ne peut pas le tester suffisement pour savoir s’il est juste, en moyenne.

Cela dit, c’est un problème assez profond, car il est très dur de savoir ce qui relève effectivement de la malchance… C’est en tous cas comme ça que je comprends le strip ci-dessous

Et je trouve dommage que sous prétexte qu’un évènement rare se réalise, on blâme la personne en charge de la prévision… Il y a manifestement une incompréhension profonde de ce que peut être le hasard. C’est triste…

Ce que la courbe ROC (et l’AUC) ne raconte pas

En préparant une intervention pour mardi prochain, j’épluchais les résultats renvoyés pour un exercice, et j’ai eu un résultat assez étrange avec un modèle de classification. J’avais donné la même base cet automne à l’ensae, et j’avais donc près d’une trentaine d’autres modèles, pour comparer (disons plutôt que sur la même base de test, j’ai une trentaine de prévisions). Les observations noires sont celles obtenues cet automne (le trait correspond aux meilleurs AUC sur la base de test), et les observations rouges sont celles que j’ai obtenu pour l’intervention de mardi (là encore, le trait vertical correspond aux meilleurs modèles, au sens du AUC), sur une observation de la base de test,

Ce sont les probabilités prédites (mais j’ai enlevé l’échelle).

Pour presque toutes mes observations, les poids rouges sont bien au dessus des autres… Mais cela n’enlève en rien le fait que l’AUC obtenu (pour les deux modèles rouges) est très bon. C’est effectivement un résultat important (et connu) : le critère AUC (mais plus généralement la courbe ROC en entier) n’indique en rien si la valeur prédite est bonne, ou pas. Il nous dit juste si l’ordre obtenu est correct. Si les valeurs les plus importantes sont effectivement les valeurs pour lesquelles on a un 1, l’AUC sera très bon.

C’est ce qu’on peut observer sur le petit exemple ci-dessous. Considérons un modèle logistique simulé assez simple,

> n=1e3
> set.seed(1)
> x1=rnorm(n)
> x2=runif(n)
> u=-3+x2+x1
> p=exp(u)/(1+exp(u))
> y=rbinom(n,prob=p,size=1)
> library(ROCR)
> df=data.frame(y,x1,x2)
> mean(df$y)
[1] 0.116
> reg=glm(y~.,data=df,family=binomial)
> p=predict(reg,type="response")
> mean(p)
[1] 0.116
> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")

L’AUC est ici

> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

et la courbe ROC est la suivante,

> plot(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="blue")

Supposons que l’on change la constante du modèle logistique,

> reg$coefficients[1]=0

Dans ce cas, notre prévision est assez mauvaise, car la probabilité moyenne prédite est ici

> u=reg$coefficients[1]+reg$coefficients[2]*
+ df$x1+reg$coefficients[3]*df$x2
> p=exp(u)/(1+exp(u))
> mean(p)
[1] 0.6060676

(on est loin des 11,6% de 1 dans la base). Pourtant le AUC est bon

> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")
> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

(c’est en fait la même valeur qu’auparavant, ce qui a du sens puisque la courbe ROC est identique)

> lines(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="red")

Autrement dit, ces outils, classiquement utilisés pour juger la qualité d’un classifieur, ne permettent en aucun cas de dire que la probabilité prédite a du sens. Ces critères permettent juste de dire qu’on identifie assez bien les personnes qui ont le plus de chance d’avoir la réponse 1. Ce qui n’est pas si mal… mais c’est un autre problème que celui d’avoir une probabilité qui soit pertinente.

Estimer une Proportion sur une Question Sensible, par Enquête

Comme le rapportait Slate hier (dans un article qui, malheureusement ne mérite pas d’être lu), un joli sujet a été proposé au baccalauréat au lycée français de Pondichéry, en mathématiques.

Ce sujet est d’autant plus intéressant qu’il revient sur une méthode assez classique pour questionner les gens sur des questions sensibles (ici sur le téléchargement pirate, mais on peut poser la même question sur la fraude, ou des pratiques sexuelles).

Comme le notait , cette technique avait été introduite initialement par Stanley Warner, voilà un demi-siècle,

Cette technique est particulièrement intéressante, tant d’un point de vu intellectuel que pratique. Plusieurs sur Twitter me faisaient remarquer que, contrairement à ce qui était raconté dans Slate, cette technique est effectivement utilisée.

Quand j’enseignais cette technique, je prenais la version simple, avec une pièce de monnaie. Supposons qu’une proportion  de personnes fraudent le fisc. Et qu’une proportion  aime Prince, va à la piscine une fois par mois, a un diplôme de l’enseignement supérieur, etc. Peut importe, la seconde proportion doit pouvoir être estimée sans problème avec une autre enquête (voire par le recensement). Bref, on a un estimateur de . Ici, les personnes interrogées tirent une pièce et répondent soit à la question “avez-vous déjà fraudé au Fisc?” s’ils tombent sur face, et à la question “aimez-vous Prince?” s’ils tombent sur pile. On pourrait dire qu’ils tirent une pièce dans leur tête, mais il faut s’assurer que la probabilité de tomber sur face est effectivement 1/2.

La magie est de noter que la proportion de personnes qui fraudent le fisc peut s’écrire

On peut avoir, par une autre enquête un estimateur du second terme, . Et en regardant des personnes qui ont répondu “oui” à une question inconnue, on a un estimateur pour

Et par différence, on en déduit un estimateur de . Les maths sont dans l’article de Stanley Warner, mais il faut reconnaître que la méthode est vraiment jolie…. Dans le cas du sujet du bac, la mise en oeuvre est un peu différente, mais l’idée est un peu la même….

(si certains ont des exemples pratiques d’application, les commentaires sont ouverts, ou sinon par mail pour joindre des documents en pièce jointe)