# Visualizing effects of a categorical explanatory variable in a regression

Recently, I’ve been working on two problems that might be related to semiotic issues in predictive modeling (i.e. instead of a standard regression table, how can we plot coefficient values in a regression model). To be more specific, I have a variable of interest $Y$ that is observed for several individuals $i$, with explanatory variables $\mathbf{x}_i$, year $t$, in a specific region $z_i\in\{A,B,C,D,E\}$. Suppose that we have a simple (standard) linear model (forget about time here) $$y_i=\beta_0+\beta_1x_{1,i}+\cdots+\beta_kx_{k,i}+\sum_j \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i$$

Let us forget the temporal effect to focus on the spatial effect today. And consider some simulated dataset. There will be only one (continuous) explanatory variable. And I will generate correlated covariates, just to be more realistic.

 n=1000 library(mnormt) r=.5 Sigma <- matrix(c(1,r,r,1), 2, 2) set.seed(1) X=rmnorm(n,c(0,0),Sigma) X1=cut(X[,1],c(-100,quantile(X[,1],c(.1,.4,.7,.85)), 100),labels=LETTERS[1:5]) X2=X[,2] Y=5+X[,1]-X[,2]+rnorm(n)/2 db=data.frame(Y,X1,X2)

Here we have $$y_i=\beta_0+\beta_1x_{1,i}+\sum_{j\in\{A,B,C,D,E\}} \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i$$ The goal here is to get to graph to visualize the vector $\hat\alpha=(\hat\alpha_A,\cdots,\hat\alpha_E)$. Let us run the linear regression

 reg1=lm(Y~X1+X2,data=db) idx=which(substr(names(reg1$coefficients), 1,2)=="X1") v1=reg1$coefficients[idx] names(v1)=LETTERS[2:5] barplot(v1,col=rgb(0,0,1,.4))

Note that it is possible to add some sort of “confidence interval” to discuss significance (or to avoid to spend hours discussing differences in bar heights that are not significantly different)

 library(Hmisc) sv1=summary(reg1)$coefficients[idx,2] (bp1=barplot(v1,ylim=range(c(0,v1+2*sv1)))) errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE) My main concern here is the “reference” that is considered. Should $A$ be the reference? Why not $B$  db$X1=relevel(db$X1,"B") reg1=lm(Y~X1+X2,data=db) idx=which(substr(names(reg1$coefficients),1,2)=="X1") v1=reg1$coefficients[idx] names(v1)=LETTERS[c(1,3:5)] library(Hmisc) sv1=summary(reg1)$coefficients[idx,2] (bp1=barplot(v1) errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

Why not the smallest one? Why not the largest one?… What if there is no simple way to choose. Furthermore, let us get back to the original point, which is that there might be some temporal aspects. More precisely, we can have $\hat\alpha^{(t)}=(\hat\alpha_A^{(t)},\cdots,\hat\alpha_E^{(t)})$. If we have also $\hat\alpha^{(t+1)}$ and we get another plot, how do we interpret it. If for $E$ the bar is taller, it means that relative to $A$, the difference has increased. I have the feeling that the interpretation is more complicated because we do not see, on that graph, changes in $\hat\alpha^{(t)}_A$.

Let us try something else. First, let us get back to the original setting

 db$X1=relevel(db$X1,"A")

Consider here the regression without the intercept, so that all values remain

 reg1=lm(Y~0+X1+X2,data=db) idx=which(substr(names(reg1$coefficients),1,2)=="X1") v1=reg2$coefficients[idx] names(v1)=LETTERS[1:5] barplot(v1)

It can be hard to read, especially if $Y$ takes (very) large values, and you think that barplots should start at 0. But still, having those 5 values is nice. Why not rescale that graph?

A natural idea my be to consider the case where no spatial component is considered, and to look at the difference with that reference.

 reg1=lm(Y~1+X2,data=db) reg2=lm(Y~0+X1+X2,data=db) idx=which(substr(names(reg2$coefficients),1,2)=="X1") v1=reg2$coefficients[idx] v2=v1-reg1$coefficients["(Intercept)"] barplot(v2,col=rgb(0,0,1,.4)) sv2=summary(reg2)$coefficients[idx,2] (bp2=barplot(v2,ylim=range(c(v2-2*sv2,v2+2*sv2)))) errbar(bp2[,1],v2,v2-2*sv2,v2+2*sv2,add=TRUE)

I like that graph, I should admit it. Now, I still have some remaining questions. For instance, can we insure that when only the intercept is considered, the value of $\hat\beta_0$ is somewhere between $\hat\beta_A,\cdots,\hat\beta_E$? Is it possible that $\hat\beta_A-\hat\beta_0,\cdots,\hat\beta_E-\hat\beta_0$ are all positive? In that case, I would find that hard to interpret.

Actually, if I really want values that can be seen as compared to some average, why not consider a (weighted) average of $\hat\beta_A,\cdots,\hat\beta_E$? (weights being here proportion in each class, in each region)

 w=table(db$X1) v3=v1-sum(w*v1)/sum(w) (bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3)))) errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE) I like that one. But what if, instead of normalizing at the end, we normalize the original dependent variable. By “normalize”, I mean “rescale”, to have a centered variable.  db$Y0=db$Y-mean(db$Y) reg3=lm(Y0~0+X1+X2,data=db) sv3=summary(reg3)$coefficients[idx,2] (bp3=barplot(v3,ylim=range(c(v3-2*sv3,v3+2*sv3)))) errbar(bp3[,1],v3,v3-2*sv3,v3+2*sv3,add=TRUE) This one is nice, because it is extremely simple to explain. But what if instead of a linear regression, we add a logistic one (with $Y\in\{0,1\}$)? or a Poisson regression… So maybe it cannot be the best solution here. Let us try something else… In insurance ratemaking, people like to use “zonier“. It is a two-stage regression. The idea is to run a regression without any spatial components, first. Then, consider the regression of residuals on spatial variables. Here, it would be something like  reg1=lm(Y~1+X2,data=db) reg2=lm(Y~0+X1+X2,data=db) Since we focus on residuals, those are centered, and we have an easy interpretation of respective values  sv4=summary(reg4)$coefficients[idx,2] v4=reg4$coefficients (bp4=barplot(v4,names.arg=LETTERS[1:5]))) errbar(bp4[,1],v4,v4-2*sv4,v4+2*sv4,add=TRUE) I guess that it can also be use in generalized linear models, with Pearson (or deviance) residuals. Another possible idea can be the following. Again, the goal is not to have the true values, but to visualize on a graph how regions can be different. Here, all of them are significantly different. And in region $A$, $Y$ is smaller, ceteris paribus (other things equal in the sense that we have taken into account $x_1$). And in region $E$ it is larger. Here, the graph helps to “see” those differences. Why not consider a completely different graph. What if we plot vector $a$ instead of $\alpha$, where $a_A$ can be interpreted as the value of the coefficient if we consider region $A$ against “not region $A$“. What if we consider 5 regressions where dichotomous versions of $Z$ are considered : $Z_j=\mathbf{1}_{Z=j}$.  v5=sv5=rep(NA,5) names(v5)=LETTERS[1:5] for(k in 1:5){ reg=lm(Y~I(X1==LETTERS[k])+X2,data=db) v5[k]=reg$coefficients[2] sv5[k]=summary(reg)$coefficients[2,2]} We can plot that sequence of values, including some confidence intervals (that would be related to significance with respect to all other regions)  (bp5=barplot(v5,ylim=range(c(v5-2*sv5,v5+2*sv5)))) errbar(bp5[,1],v5,v5-2*sv5,v5+2*sv5,add=TRUE) Looking at values does not give intuitive results, but I have the feeling that it is easy to explain what we plot (we compare each region to “the rest of the world”), and the ordering of $a$ seems to be consistent with $\alpha$ (but I could not prove it). Here are some ideas I got. I should be able to provide other graphs, but I would love to discuss with anyone on that topics, to find a proper and nice way to visualize effects of a categorical explanatory variable in a regression model (that can be a logistic one). Comments are open… # Holt-Winters with a Quantile Loss Function Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord & Snyder (2008)), when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class. Furthermore, I like it because I think it has nice pedagogical features. Consider simple exponential smoothing, $$L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1}$$ where $\alpha\in(0,1)$ is the smoothing weight. It is locally constant, in the sense that ${}_{t}\hat Y_{t+h} = L_{t}$  library(datasets) X=as.numeric(Nile) SimpleSmooth = function(a){ T=length(X) L=rep(NA,T) L[1]=X[1] for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]} return(L) } plot(X,type="b",cex=.6) lines(SimpleSmooth(.2),col="red") When using the standard R function, we get  hw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1]) hw$alpha [1] 0.2465579

Of course, one can replicate that optimal value

 V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(erreur^2)) } optim(.5,V)$par [1] 0.2464844 Here, the optimal value for $\alpha$ is the one that minimizes the one-step prediction, for the $\ell_2$ loss function, i.e. $$\sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2$$ where here ${}_{t-1}\hat Y_t = L_{t-1}$. But one can consider another loss function, for instance the quantile loss function, $$\ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0})$$. The optimal coefficient is then obtained using  HWtau=function(tau){ loss=function(e) e*(tau-(e<=0)*1) V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(loss(erreur))) } optim(.5,V)$par }

Here is the evolution of $\alpha^\star_\tau$ as a function of $\tau$ (the level of the quantile considered).

 T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(T,HW,type="l") abline(h= hw$alpha,lty=2,col="red") Note that the optimal $\alpha$ is decreasing with $\tau$. I wonder how general this result can be… Of course, one can consider more general exponential smoothing, for instance the double one, with $$L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}]$$and$$B_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}$$so that the prediction is now ${}_{t}\hat Y_{t+h} = L_{t}+hB_t$ (it is now locally linear – and no longer constant).  hw=HoltWinters(X,gamma=FALSE,l.start=X[1]) hw$alpha alpha 0.4200241 hw$beta beta 0.05973389 The code to compute the smoothed series is the following  DoubleSmooth = function(a,b){ T=length(X) L=B=rep(NA,T) L[1]=X[1]; B[1]=0 for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] } return(L+B) } Here also it is possible to replicate R using the $\ell_2$ loss function  V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(erreur^2)) } optim(c(.5,.05),V)$par [1] 0.41904510 0.05988304

(up to numerical optimization approximation, I guess). But here also, a quantile loss function can be considered

 HWtau=function(tau){ loss=function(e) e*(tau-(e<=0)*1) V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(loss(erreur))) } optim(c(.5,.05),V)$par } and we can plot those values on a graph  T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(HW[1,],HW[2,],type="l") abline(v= hw$alpha,lwd=.4,lty=2,col="red") abline(h= hw$beta,lwd=.4,lty=2,col="red") points(hw$alpha,hw$beta,pch=19,col="red") (with $\alpha$ on the $x$-axis, and $\beta$ on the $y$-axis). So here, it is extremely simple to change the loss function, but so far, it should be done manually. Of course, one do it also for the seasonal exponential smoothing model. # Justice “actuarielle”, algorithmes… et données Il y a un peu plus d’un an, Virginie Gautron m’envoyait plein de documents sur la “justice actuarielle”, concept que je découvrais alors. Pour comprendre un peu de quoi il s’agit, je peux renvoyer vers pénalité et gestion des risques : vers une justice « actuarielle » en Europe ? qui dresse un état de l’art, en francais. J’avoue avoir un peu mis de côté ensuite (faute de temps), et puis l’autre jour, j’ai dévoré surveillé et punir de Nicolas Bourgoin, qui contient un chapitre sur le sujet. Sans vouloir tomber dans le cliché des deux cultures (à la CP Snow), il faut bien admettre que pour un économètre, lire des écrits de juristes parlait de modèles prédictifs (ou actuariels) est une expérience marquante. Et ce qui me trouble (j’ai d’ailleurs le même soucis avec la psychologie – c’est un sujet proche) c’est que l’accent est souvent mis sur le modèle, l’algorithme, ce truc qui semble fasciner les juristes, alors qu’il me semble que le coeur du problème, ce sont les données. On nous dit, par exemple, que les algorithmes sont racistes (c’est un peu la thèse de Cathy O’Neil dans Weapons of Math Destruction), mais le problème ce n’est pas le modèle. C’est le fait qu’il y a un énorme biais de sélection dans les données ! Et selon l’adage “garbage in, garbage out“… dit de manière plus policée, je pourrais tenter un “purgamentum init, exit purgamentum” pour parler comme un juriste. Ou dit autrement, si on entre des mauvaises données en entrée, faut pas espérer grand chose en sortie… Pour comprendre le biais de sélection, on peut relire Mostly Harmless Econometrics de Joshua Angrist et Jörn-Steffen Pischke. En particulier, dans l’introduction, ils posent la question simple suivant : vaut-il mieux, suite à une opération, renvoyer les gens chez eux, ou bien les laisser en chambre de convalescence, à l’hôpital ? Le soucis est que les données sont biaisées : oui, si on regarde les données brutes, ceux qui sont rentré chez eux se sont remis plus vite… mais peut être a-t-on autorisé les gens à rentrer justement parce qu’ils allaient bien ? Et on retrouve ce biais de sélection un peu partout. Probablement même encore davantage dans les problèmes relatifs à la justice. Prenons un exemple simple. On va s’intéresser à la probabilité d’être coupable d’un crime, et d’être condamné par la justice. C’est ma variable $Y$. Pour avoir un modèle économétrique, il me faut des variables explicatives. J’ai une variable $X_3$ et une autre $X_4$, la première influençant vraiment la probabilité de commettre un crime, et l’autre étant juste du bruit. Maintenant, entrons dans le coeur du modèle. Supposons que j’ai une variable corrélée à ma richesse $\tilde{X}_2$, comme le lieu d’habitation, un critère racial, ou pour faire simple le diplôme. Notons $X_1$ le niveau de diplôme de la personne (mais encore une fois, on peut imaginer n’importe quoi, l’idée est que cette variable est corrélée à la richesse). Et surtout, cette variable n’est pas du tout une variable causale dans mon modèle. La probabilité de commettre un crime n’est pas influencée par $X_1$. Une personne ayant un diplôme de second cycle universitaire a autant de chance de commettre un crime qu’une personne ayant raté le brevet. On va même supposer que la probabilité de commettre un crime soit indépendante de la richesse. En fait, la richesse me permet de payer un bon avocat. Et avoir un bon avocat permet de diminuer la probabilité d’être condamné. Pour faire simple, j’ai deux niveaux de richesse, qui me donne soit un avocat commis d’office, soit un avocat qui s’impliquera peut être davantage (je simplifie à des fins pédagogiques). Bref, mon modèle aura la forme suivante : la probabilité d’être condamné pour un crime s’écrit $$p=(0.025+0.05\cdot X_3)\cdot[k\cdot \mathbb{1}(X_2=R)+\mathbb{1}(X_2=P)]$$ Autrement dit, ma probabilité est une fonction croissante (linéaire) de $X_3$, et si je suis riche (ce qu’on note $X_2=R$), j’ai moins de chance d’être condamné, avec un facteur $k\in[0,1]$, toutes choses étant égales par ailleurs, ceteris paribus, histoire de parler latin comme les juristes. Pour simuler un jeu de données, on considère une population d’un million de personnes, et je suppose que ma richesse est très corrélée avec mon niveau de diplôme. Pour info, la classe $A$ ce sont les gens les plus diplômés, et je suppose qu’un bon avocat divise par trois la chance d’être condamné.  n=1000000 r=.95 k=1/3 S=matrix(c(1,r,r,1),2,2) library(mnormt) set.seed(1) vectX=rmnorm(n,varcov=S) X1=cut(vectX[,1],breaks=c(-10,1,-1,10), c("C","B","A")) X1=relevel(X1,"A") X2=cut(vectX[,2],breaks=c(-10,1,10),c("P","R")) X3=runif(n) X4=runif(n) Y1=rbinom(n,size = 1,prob=(.025+.05*X3)*(k*(X2=="R")+(X2=="P"))) B0=data.frame(X1,X2,X3,X4,Y1) On peut faire un peu de statistiques pour regarder les caractéristiques des gens condamnés (uniquement)  B1=data.frame(X1=X1[Y1==1],X2=X2[Y1==1],X3=X3[Y1==1],X4=X4[Y1==1]) Par exemple si je regarde la proportion de gens les plus diplômés: ils sont 15% dans la population totale, et 5% parmi les personnes condamnées  mean(B0$X1=="A") [1] 0.158734 mean(B1$X1=="A") [1] 0.05016072 t.test(B0$X1=="A",B1$X1=="A") data: B0$X1 == "A" and B1$X1 == "A" t = 219.5, df = 1283500, p-value < 2.2e-16 95 percent confidence interval: 0.1076038 0.1095428 sample estimates: mean of x mean of y 0.15873400 0.05016072 Autrement dit, on a une différence significative entre ces deux pourcentages. Pareil si je regarde la valeur moyenne de $X_3$, qui était de 50% dans la population et de 58% dans la population condamnée  t.test(B0$X3,B1$X3) data: B0$X3 and B1$X3 t = -163.46, df = 844940, p-value < 2.2e-16 95 percent confidence interval: -0.08455430 -0.08255058 sample estimates: mean of x mean of y 0.5000094 0.5835618 alors que pour $X_4$, on a 50% dans la population totale, et 50% dans la population condamnée. Mais on s’y attendait….  t.test(B0$X4,B1$X4) data: B0$X4 and B1$X4 t = -0.81161, df = 48645, p-value = 0.417 95 percent confidence interval: -0.003869879 0.001603452 sample estimates: mean of x mean of y 0.5001326 0.5012659 Passons à la seconde phase. Il semble y avoir un vrai intérêt dans la littérature sur les modèles actuariels pour la construction de modèles liés à la probabilité de récidive. Là encore, je vais faire une hypothèse forte : faire de la prison n’influence en rien la probabilité de récidiver. Autrement dit, pour cette population déjà condamnée, je vais refaire tourner mes simulations, avec une probabilité de commettre un crime qui dépend seulement de $X_3$, mais avec toujours notre effet richesse qui permet de payer un avocat qui fera baisser la probabilité d’être condamné.  set.seed(2) n1=nrow(B1) Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P"))) B2=cbind(B1,Y2) reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

J’ai construit ici mon modèle permettant de calculer une probabilité de récidive. Le modèle est construit sur le niveau de diplôme $X_1$ (qui n’influence en rien, encore une fois, la probabilité de commettre un crime), la variable $X_3$ (qui effectivement joue un rôle) et la variable de bruit $X_4$.

 set.seed(2) n1=nrow(B1) Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P"))) B2=cbind(B1,Y2) reg=glm(Y2~X1+X3+X4,data=B2,family=binomial) Je fais ici un score logistique car c’est le plus logique.  summary(reg) Coefficients: Estimate Std. Error z value Pr(|z|) Inter.) -4.04753 0.11592 -34.916 2e-16 *** X1C 0.61311 0.10914 5.618 1.93e-08 *** X1B 0.68807 0.09977 6.896 5.34e-12 *** X3 1.00610 0.08052 12.495 2e-16 *** X4 -0.13562 0.07286 -1.861 0.0627 . Premier élément à observer, la variable $X_4$ est à la limite d’être significative… Certains auteurs sur ces thèmes diraient que c’est un signal faible. Pourtant, en l’occurence, c’est du bruit ! Juste du bruit !! Et parmi les signaux forts, très clairement significatifs, on retrouve $X_1$, notre niveau de diplôme. Ici, si on est faiblement qualifié (niveaux $B$ ou $C$) alors on a une chance de récidive beaucoup plus importante ! C’est très fortement significatif ! Alors que pour rappel, la probabilé de commettre un crime, et de récidiver n’est aucunement causé par cette variable ! C’est juste une variable corrélée avec la variable richesse, qui permet de faire baisser la probabilité d’être condamné. Plus intéressant, si je regarde la valeur du coefficient, en fonction de $k$, autrement dit la valeur d’un bon avocat (que l’on peut s’offrir si on est riche) on obtient la courbe suivante. Autrement, sur la partie droite, même si un bon avocat ne baisse que de 5% la chance d’être condamné, le niveau de diplôme semble être significativement liée à la probabilté de récidive. Avec un pic autour de 30%. Si les avocats ont un pouvoir incroyable, par exemple en divisant par 20 la probabilité d’être condamné, étrangement, la variable de diplôme redevient très peu significative, mais je pense que c’est parce qu’on n’a alors plus personne de très diplômé dans la base. On pourrait me reprocher d’avoir une corrélation trop forte entre richesse et niveau de diplôme, ici. En mettant une corrélation à 70% (dans mon modèle Gaussien latent), j’obtiens qui donne un effet moindre, mais là encore, si un bon avocat permet de réduire la probabilité d’être condamné par un facteur compris entre 5% et 90%, alors la variable de diplôme semblera significative. Mon modèle n’est pas faux ici. C’est essentiellement que j’ai un biais dans mes données. Et ne pas en tenir compte fait que le modèle n’a aucun sens, et ne peut pas être interprété sans commettre d’erreurs (et ici je cumule les âneries, en voyant un signal faible alors qu’il y a du bruit, et en voyant un signal fort pour une variable qui n’a aucun lien causal). J’insiste lourdement, mais le diplôme n’intervient aucunement ici. Et on peut remplacer le diplôme par toute variable corrélée à la richesse, comme locataire et propriétaire, etc. Bref, l’utilisation de modèles économétrique dans des problèmes judiciaires n’est pas simple ! Je parle ici de biais de sélection, mais la formulation qui me semble la plus juste serait une mauvaise spécification du modèle (on n’utilise pas les bonnes variables, mais ces dernières sont peut-être non-observables – ou difficiles à justifier éthiquement : peut-on demander le revenu de la personne pour faire un score de récidive, et avoir un score qui dépende effectivement du revenu ?). Il y a aussi un problème d’identification, car on a du mal à savoir ce qui influence le fait de commettre un crime, et la condamnation. Ce qu’on observe, ce sont les condamnations, pas les actes (je ne parle pas des condamnations à tort ici, que j’avais pu évoquer dans un autre billet, sur le principe de précaution). Mais j’essayerais de revenir sur ce point après les vacances…. Dernière anecdote pour finir. Après avoir découvert ces modèles via Virginie Gautron, j’avais discuté avec une amie, présidente de tribunal, qui me parlait de la nullité des experts qu’elle avait sous la main. Et elle voulait faire un peu de ménage. Bref, elle me demandait si je pouvais faire un modèle attribuant un score à un expert. Et là encore, ma réponse a été que oui, faire un modèle c’est facile, à condition d’avoir des données. Comment puis-je avoir un score $Y$ d’un expert? Qui sait ce qu’est un bon expert? C’est quelqu’un contre qui on ne demande pas de contre-expertise? Malheureusement, la justice ne fonctionne pas ainsi, et la demande n’a rien à voir avec la qualité de l’expert. Et la qualité, vu du point de vue de la juge, est une notion qui peut être très différente de celle des deux parties qui s’opposent… Bref, je pourrais facilement faire un modèle. Mais en avoir un qui ait du sens, c’est une autre paire de manches ! # The myth of interpretability of econometric models There are important discussions nowadays about data modeling, to choose between the “two cultures” (as mentioned in Breiman (2001)), i.e. either econometrics models or machine/statistical learning models. We did discuss this issue recently in Econométrie et Machine Learning (so far only in French) with Emmanuel Flachaire and Antoine Ly. One argument often used by econometricians is the interpretability of econometric models. Or at least the attempt to get an interpretable model. We also have this discussion in actuarial science, for instance in ratemaking (or insurance pricing). Machine learning based models usually perform better (for some a priori chosen metric), but actuaries claim that econometric models are more easily interpretable. In actuarial literature, we assume that claim frequency $Y$ is driven by some non-observable risk factor $\Theta$, and therefore, we do have heterogeneous risks in our portfolio. And, it can be seen as legitimate to differentiate prices. Assume that this risk factor $\Theta$ is strongly correlated with $X_1$, the age of the driver. Because in our portfolio, old drivers tend to have more accidents. Here, we could pretend to have a “causal story” (as defined in Freedman (2009)) because of a possible interpretation of the model. So it is natural here to consider a regression model of $Y$ on $X_1$ to derive our actuarial pricing model. But assume that, possibly, risk factor $\Theta$ is also strongly correlated with $X_2$, that can be related to spatial features (say latitude, which denoted a north/south position). Because in our portfolio, drivers living in the south tend to have more accidents (reads are known to be more dangerous there). Here, we could pretend to have a second “causal story”. Of course, since $\Theta$ is strongly correlated with $X_1$ and $X_2$, it means that $X_1$ and $X_2$ are strongly correlated. Here also, this correlation can be interpreted (not in a causal way as previously, but still), since we know that old people like to live in southern regions. So, what should we do here ? Let us run some simulations to illustrate.  set.seed(123) n=1e5 Theta=rnorm(n) X1=Theta+rnorm(n)/8 X2=Theta+rnorm(n)/8 L=exp(-3+Theta) Y=rpois(n,L) B=data.frame(Y,X1,X2) Our first idea was to consider a model where $Y$ is “explained” by the first variable $X_1$,  g1=glm(Y~X1,data=B,family=poisson) summary(g1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Inter.) -2.97778 0.01544 -192.88 <2e-16 *** X1 0.97926 0.01092 89.64 <2e-16 *** As expected, our variable is “significant”, but also, probably more interesting, $X_2$, has no impact on the residuals  B$e1=residuals(g1,type="pearson") g1e=lm(e1~X2,data=B) summary(g1e)   Coefficients: Estimate Std. Error t value Pr(>|t|) (Inter.) 0.0003618 0.0031696 0.114 0.909 X2 0.0028601 0.0031467 0.909 0.363

The interpretation is that once we corrected claim frequency for the age of the drivers, there is no spatial effect here. So, a good model should be based only on the age of the drivers.

But we can also consider the other story. We can consider a model where $Y$ is “explained” by the second variable $X_2$,

 g2=glm(Y~X2,data=B,family=poisson) summary(g2)   Coefficients: Estimate Std. Error z value Pr(>|z|) (Inter.) -2.97724 0.01544 -192.81 <2e-16 *** X2 0.97915 0.01093 89.56 <2e-16 ***

Here also we have a valid model, that can be interpreted, and here also $X_1$, has no impact on the residuals

 B$e2=residuals(g2,type="pearson") g2e=lm(e2~X1,data=B) summary(g2e) Coefficients: Estimate Std. Error t value Pr(>|t|) (Inter.) 0.0004863 0.0031733 0.153 0.878 X1 0.0027979 0.0031504 0.888 0.374 The story is similar here. If we correct from the spatial pattern, claims frequency does not depend on the age of the driver. So, what should we do now? We do have two models, and each of them is as interpretable as the other one. Note that we can not use any statistical tool to distinguish the two: they are comparable  AIC(g1) [1] 51013.39 AIC(g2) [1] 51013.15 Why not incorporate the two explanatory variables $X_1$ and $X_2$, at the same time, in our regression model, and let “the model” decide what to do…?  g=glm(Y~X1+X2,data=B,family=poisson) summary(g) Coefficients: Estimate Std. Error z value Pr(>|z|) (Inter.) -2.98132 0.01547 -192.723 2e-16 *** X1 0.49310 0.06226 7.920 2.38e-15 *** X2 0.49375 0.06225 7.931 2.17e-15 *** It looks like we completely lost the interpretability of the model, since our two explanatory variables are (strongly) correlated. Actually, instead of saying “use one, and drop the other one (since it brings no further information)”, it says “use both, each one will explain half of the variable”. Strange interpretation, isn’t it? So why not try some LASSO here?  library(glmnet) fit=glmnet(x=as.matrix(B[,c("X1","X2")]), y=B$Y,family="poisson") plot(fit,xvar="lambda")

Here also, it says that we either keep both, or none. So it cannot be used for variable selection (which is an important motivation to use LASSO technique). So, what should be do if we several interpretable models, but no way to choose? Because usually, we claim that we prefer to use a model with an interpretation. But what should be done here?

# A la recherche des homonymes…

Il y a quelques mois, Baptiste Coulmont m’avait contacté avec une question passionnante (comme chaque fois qu’il me contacte).

Baptiste voulait travailler sur la proportion de personnes qui ont un homonymes dans une population de taille donnée, ou sur la probabilité de ne pas avoir d’homonymes dans un bureau de vote, par exemple. Or ce dernier problème n’est pas sans rappeler le “paradoxe des anniversaires“.

Dans un groupe de 23 personnes, il y a plus d’une chance sur deux pour que deux personnes de ce groupe aient leur anniversaire le même jour. Ou dit autrement, il est plus surprenant de ne pas avoir deux personnes qui ont leur anniversaire le même jour que d’avoir deux personnes qui ont leur anniversaire le même jour.

On parle de “paradoxe” parce que le nombre de jours possibles est relativement grand ($k=365$), mais que la taille du groupe est elle relativement petite ($n=23$). Ici c’est un peu la même chose, sauf qu’on est sur des nombres beaucoup plus grands: on a des milliers de noms possibles (en fait, on s’intéresse au couple prénom+nom, ce qui augmente très largement le nombre de possibilités, que l’on notera $k$ ici), on considère des groupes relativement grands (la taille du groupe sera elle notée $n$), et surtout, l’hypothèse (implicite) d’uniformité dans la distribution des dates ne marche pas sur les noms de famille, et les prénoms. On peut malgré tout tenter de formaliser le problème.

Considérons une collection de $k$ éléments, notée $\mathcal{X}=\lbrace x_1,\cdots,x_k\rbrace$, et faisons $n$ tirages avec remise. Les tirages obtenus sont notés ici $X_1,\cdots,X_n$ et correspondent à des variables aléatoires iid, suivant des lois multinomiales sur $\mathcal{X}$ avec comme vecteur de probabilités $\mathbf{p}=(p_1,\cdots,p_k)$. Dans le cas des anniversaires, $\mathcal{X}$ correspond aux dates, avec $k=365$, et $\mathbf{p}$ est uniforme sur $\mathcal{X}$ – ou encore $\mathbf{p}=\mathbf{1}/k$ . Le problème classique est de calculer la probabilité que dans un groupe de taille donnée tout le monde ait des dates d’anniversaire différentes

$$\mathbb{P}[\forall j\neq i, X_i\neq X_j]$$

mais on peut aussi calculer (comme le voulait Baptiste) la proportion de personne ayant un alter-ego (une personne dans le groupe qui a son anniversaire le même jour)
$$\frac{1}{n}\sum_{i=1}^n \mathbf{1}(\exists j\neq i,X_i=X_j)$$

On a donc tenté d’approcher le second, quand on a un groupe de taille $n$. Notons $Z_i$ le nombre de personnes (parmi les $n-1$ restantes) qui partagent la même date d’anniversaire que $i$,
$$Z_i=\sum_{j\in\{1,\cdots,n\}\backslash\{i\}} \mathbf{1}(X_j \neq X_i)$$
Aussi, la personne $i$ a un alter-ego si $Z_i\geq 1$. Et donc, la proportion de personnes avec un alter-ego dans le groupe s’écrit
$$P_{n,k}=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(Z_i\geq 1)$$
Si $\mathbf{p}$ correspond à la loi uniforme sur $\mathcal{X}$, alors
$$1-P_{n,k}=\binom{k}{n}\frac{1}{k^n}\sim \exp\left[\binom{n+1}{2}\frac{1}{k}\right]\sim e^{-n^2/2k}$$
comme le montre Chatterjee et al. (2004). Je pourrais aussi renvoyer aux articles de Inoue & Aki (2008) et plus récemment de  Cortina Borja (2013) pour des revues de littérature sur les aspects computationnels.

Mais la plupart des résultats sont obtenu dans le cas où $\mathbf{p}$ correspond à une distribution uniforme. Les propriétés plus générales sont beaucoup plus rares(on pourra lire Munford (1977), Das Gupta (2004),  Inoue & Aki (2008) ou Nunnikhoven (1992) pour des tentatives. Utiliser des méthodes de Monte Carlo est probablement le plus naturel (si on n’est pas trop pressé) en se donnant un ensemble $\mathcal{X}$ (qui sont les dates dans le problème des anniversaires et les noms dans le cas des homonymes) de taille $k$, et une distribution $\mathbf{p}$ sur $\mathcal{X}$. On a ainsi pu approcher $P_{n,k}$, la proportion de personnes ayant un alter-ego dans une population de taille $n$: sur les figures ci-dessous, on a cette proportion en fonction de  $n$ (en abscisse), avec des courbes différentes pour différents $k$. Histoire d’être un peu réaliste, on va chercher le taux de personnes ayant un homonyme dans une popoluation de plusieurs centaine de milliers de personnes, et quand aussi des milliers, des dizaines de milliers, de possibilités.

Dans le cas d’une distribution uniforme, la proportion évolue comme sur la figure ci-dessous

On a tenté de quitter l’indépendance en utilisant une loi paramétrique, plus précisément une loi de type Pareto, ou Zipf, pour $\mathbf{p}$. En fait, c’est le type de distribution qu’avait observé Li (2012) sur la distribution des prénoms, avec $p_i \propto i^{-\alpha}$ (ou plus récemment, Höhle (2017)). On a alors l’évolution suivante

On est ici en échelle logarithmique pour la taille $n$ parce que c’est naturel. Mais pour l’instant, on était resté sur $P_{n,k}$. Baptiste insistait pour passer aussi à une échelle logarithmique sur l’axe des ordonnées, mais ca n’etait pas concluant. On a alors tenté deux autres représentations, plus classiques quand on travaille sur des probabilités (ou des proportions) : représenter $\text{logit}(P_{n,k})$ en fonction du logarithme de $n$ (où $\text{logit}(u)=\log[u/(1-u)]$) et une version probit, à savoir $\Phi^{-1}(P_{n,k})$ (où ici $\Phi$ est la fonctionc cumulative de la loi normale centrée et réduite).

En grande dimension, on a l’impression d’avoir une relation linéaire entre la transformée probit de la proportion de personnes ayant un homonyme et le logarithme de la taille de la population considérée

$$\Phi^{-1}(P_{n,k}) \sim a_{k}+b_{k} \log[n]$$

C’est joli, mais on se demandait comment utiliser ca sur nos données, pour travailler vraiment sur les homonymes. Le soucis, c’est qu’on n’a pas de fichier national, avec les noms et les prénoms. Alors Baptiste a sorti de son sac les listes électorales parisiennes. Le soucis est que la base est (relativement) petite. Le nombre de couple est très limité. Au lieu de faire des tirages avec remise dans les couples (“prénom”, “nom”) – encore une fois, la base est trop petite – on a voulu tenter de tirer indépendemment le “prénom” et le “nom”. On y gagne en ayant à disposition une base plus grande, mais on introduit un biais colossal, les “prénom” et “nom” étant assez corrélés. Dans notre cas, on a $k_1$ “prénom” et $k_2$ “prénom”. Il y a quelques années Chalabi & Flowers (2014)} (repris plus en détails sur ce blog) était revenu sur la non indépendance des noms et des prénoms, en s’interrogeant sur le nom le plus populaire aux Etats-Unis. Selon les auteurs “Michael Smith” n’était probablement pas le nom le plus populaire, même si “Michael” et “Smith” étaient respectivement les “prénom” et “nom” les plus populaires. En France “Jean” et “Marie” sont les prénoms les plus populaires pour les garcons et les filles, respectivement, et “Martin” est le nom le plus populaire (au moins sur les listes électorales parisiennes), mais en fait, le nom le plus populaire est “Thi Nguyen”, avec 350 personnes à Paris, alors que seulement 84 “Jean Martin” sont observés, le troisième plus haut score. Ah oui, je dois préciser qu’on se limite ici au premier prénom, en cas de prénoms composés. Tout comme Chalabi & Flowers (2014), on a fait un test du chi-deux pour tester l’indépendance, sur les 20 “noms” et “prénoms” les plus populaires. Sur la figure ci-dessous, on a les résidus normalisés du test. Par exemple “Michel Michel”est en rouge, car c’est beaucoup moins donné que ce qu’on aurait si  les deux composantes étaient indépendantes

Bref, on n’a pas trop le droit de tirer le “prénom” et le “nom” indépendament. Mais on va quand même le faire, pour voir. Ici $\mathcal{X}$ est l’ensemble $\mathcal{X}_1\times\mathcal{X}_2$ avec $k_1$ choix pour le “prénom” et $k_2$ pour le “nom”.

Soit $\mathbf{p}=(p_{i_1,i_2})$ la probabilité empirique sur $\mathcal{X}$. On imagine que $\mathbf{p}$ contient énormément de 0, car beaucoup de couples ne sont pas observés. Pour les simulations, on va noter $\mathbf{p}^{\perp}$ la pobabilité jointe sous hypothèse d’indépendance
$$p^{\perp}_{i_1,i_2}=\left(\frac{1}{k_2}\sum_{j=1}^{k_2} p_{i_1,j}\right)\cdot\left(\frac{1}{k_1}\sum_{j=1}^{k_1} p_{j,i_2}\right)=p_{i_1,\cdot}p_{\cdot,i_2}$$

Lors de nos simulations, on va considérer deux cas : on dire des $\mathbf{x}=(x_1,x_2)$ soit suivant la distribution $\mathbf{p}$, soit suivant $\mathbf{p}^{\perp}$. Bon, j’avais formalisé pour clarifier un peu dans ma tête et aussi pour impressionner Baptiste, mais techniquement, dans un cas on tire une paire observée – genre “Baptiste Coulmont” – et dans le second cas, on tire d’un côté un prénom – genre “Baptiste” – et de l’autre un nom de famille – genre “Charpentier” (même s’il n’existe aucun “Baptiste Charpentier”, on dit juste qu’il serait possible qu’il y en ait un).

Bref, on tirant des échantillons (plusieurs milliers de taille donnée, puis en moyennant), on obtient la figure ci-dessous (avec une échelle logarithme pour l’axe des ordonnées)

La courbe noir sur-estime sûrement à cause du faible nombre de choix possibles (on va tirer jusqu’à 50 millions de noms). En revanche, l’indépendance n’étant pas réaliste, la courbe rouge sous-estime probablement la proportion. Mais le plus intéressant est le graphique ci-dessous, obtenu en représentant la version probit de la probabilité

La courbe rouge est incroyablement linéaire ! Bien plus linéaire que ce qu’on espérait avoir suite à nos simulations dans le cas de la loi de Pareto.

En fait, on note que si on se limite au cas où $n$ n’est pas trop grand, la courbe noire semble aussi linéaire. Ensuite on perd la linéarité parce qu’on tire avec remise dans une base trop petite, mais au début (lorsque $n$ n’est pas trop grand), on pourrait accepter l’hypothèse d’une relation linéaire entre la transformée probit, et le logarithme de $n$. Bref, on estime une relation linéaire, et on extrapole pour les valeurs de $n$ plus grandes

Si on revient sur la représentation plus classique, de la proportion de personnes ayant un homonyme en fonction de $n$, on obtient

Autrement dit, en France, on peut estimer que 80% des francais ont un homonyme.

Deux points pour poursuivre, sur la robustesse de l’approche. Tout d’abord, si on change un peu la taille de l’intervalle des $n$ pour lesquels la courbe noire semble linéaire, la pente change peu

et on a un ordre de grandeur assez proche. Bref, ce 80% semble assez robuste,

Ensuite, comme cette relation linéaire me surprenait, on a aussi tenté sur des données américaines… et on a exactement le même genre de résultats…

Pour aller plus loin, on a tout mis sur github (avec des graphs en ggplot2 pour faire plus joli…) et une version (en anglais) est en ligne sur arxiv.

# C’est normal ! (partie 2) et si la normalité n’existait pas?

Avec un peu de retard, je vais tenter de poursuivre le billet sur la normalité, mis en ligne en avril dernier. Mon point de départ reprenait une discussion autour de la « moyenne » comme outil central pour définir la norme, ainsi que des idées que l’on peut retrouver dans le livre de Todd Rose, The End of Averages. Pour compléter mon précédant billet, je pourrais reprendre un exemple mentionné dans ce livre, revenant sur des problèmes rencontrés par l’Armée Américaine dans les années 1950. Lors de la conception des postes de pilotage des avions de chasse, les ingénieurs avaient utilisé les dimensions de plus de 4000 pilotes pour positionner de manière optimale le siège par rapport aux pédales, le manche à balai, la hauteur du pare-brise, mais aussi la forme du siège, du casque, etc. Ces mesures ont permis de calculer les mensurations du pilote « médian » dans une dizaine de dimensions. Par exemple la taille moyenne des pilotes était de 179cm, ce qui a permis de définir la taille d’un pilote moyen entre 175 et 185cm. Si une majorité des pilotes est de taille moyenne, parmi les 4000 pilotes, aucun (parmi les 4000 !) n’était « moyen » dans toutes les dimensions. Comme l’affirmait Daniels (1952), « concevoir un poste de pilotage pour le pilote moyen revenait en fait à n’en concevoir pour aucun ». Je faisais ensuite le lien avec le concept de malédiction de la dimension. Mais je voulais revenir aujourd’hui sur le sens (plus philosophique) de la « normalité ».

De la normalité statistique à la norme, de la régularité à la règle

D’un point de vue empirique, descriptif, être dans la norme, ne signifie rien d’autre qu’être dans la moyenne, ne pas se détacher trop de cette moyenne. On aura alors tendance à définir la norme comme la fréquence de ce qui se produit le plus souvent, comme l’attitude la plus fréquemment rencontrée ou de la préférence la plus régulièrement manifestée. Mais cette normalité n’est pas la normativité, et « être dans la norme », être exemplaire, relève alors d’une dimension différente, qui se rattache cette fois non plus à une description du réel mais à une identification de ce vers quoi il doit tendre. On passe donc du registre de l’être à celui du devoir-être, du « is » au « ought » pour reprendre la terminologie de Hume (1739). Il est en effet difficile d’envisager le modèle (ou la normalité) sans glisser vers ce second sens qu’on peut trouver au concept de norme, et qui déploie quant à lui une dimension à proprement parler normative. Cette vision mène à une confusion entre normes et lois, même si toute la normativité n’est pas épuisée par les lois. Hume constate ainsi que, dans tous les systèmes de morale, les auteurs passent de constat de faits, c’est-à-dire énonciatifs de type « il y a », à des propositions qui comportent une expression normative, comme « il faut », « on doit ». Ce que Hume conteste, c’est le passage d’un type d’affirmation à un autre : pour lui, ce sont là deux types d’énoncés qui n’ont rien à voir les uns avec les autres, et qu’on ne peut donc pas enchaîner logiquement les uns avec les autres, en particulier d’une norme empirique vers une règle normative. Pour Hume, une affirmation qui ne serait pas normatives ne peut pas donner lieu à une conclusion de type normatif. Cette affirmation de Hume a suscité nombre de commentaires et d’interprétations, en particulier parce qu’énoncé telle quel, elle semble être un obstacle à toute tentative de naturalisation de la morale (comme le détaillent Mc Intyre (1959) ou Rescher (1990)). En ce sens, on trouve la distinction forte entre la norme dans la régularité (normalité) et la règle (normativité).

Les lois statistiques, du micro au macro

La loi statistique relève de ce qui est « est » parce qu’on l’a observé ainsi (par exemple « les hommes sont plus grands que les chiens »). La loi humaine (divine, ou judiciaire) relève de ce qui est « est » parce qu’on l’a décrété, et donc « doit être » (« les hommes sont libres et égaux » ou « l’homme est bon »). Enfin, la loi physique relève de ce qui « est » parce qu’on peut le montrer (« les planètes sont attirées entre eux »), souvent dans le cadre d’hypothèses. On voit que les trois concepts peuvent être liés. Par exemple la « loi de Kepler » a été établie historiquement à l’aide d’observations (et rentrait historiquement dans la première catégorie), avant d’être démontrée dans le cadre du modèle Copernicien (et de passer alors dans la troisième). On peut d’ailleurs associer un concept d’équilibre à cette loi, cette « norme ». Toutefois, comme le souligne Hilpinen (1971), les lois probabilistes posent toutefois bon nombre de questions, il suffit de penser aux lancers de dés, ou aux attentes : que signifie « il est normal d’attendre 5 minutes le bus à l’arrêt », ou plus dérangeant d’un point de vue éthique, « il est normal qu’une personne placée en détention provisoire soit emprisonnée 18 mois » ?

On peut voir la « norme » comme une régularité de cas, observée à l’aide de fréquences (ou de moyennes), par exemple, sur la taille des individus, la durée du sommeil, autrement dit des données qui constituent la description d’individus. Les données anthropométriques ont ainsi permis de définir par exemple une taille moyenne des individus dans une population donnée, en fonction de leur âge ; par rapport à cette taille moyenne, un écart de 20% en plus ou en moins détermine le gigantisme ou le nanisme. Si l’on pense aux accidents de la route, il peut être considéré comme « anormal » d’avoir un accident de la route une année donné, à un niveau individuel (micro), car la majorité des conducteurs n’ont pas d’accident. Néanmoins, du point de vue de l’assureur (macro), la norme est que 10 % des conducteurs aient un accident. Il serait donc anormal que personne n’ait d’accident. C’est l’argument que l’on retrouve dans Durkheim (1897). De l’acte singulier qu’est le suicide, s’il est considéré du point de vue de l’individu qui le commet, Durkheim tente de le voir comme un acte social, relevant alors d’une réelle régularité, au sein d’une société donnée. Dès lors, le suicide devient, selon Durkheim, un phénomène « normal ». Les statistiques permettent alors de quantifier la tendance au suicide dans une société donnée, dès lors que l’on observe non plus l’irrégularité qui apparaît dans la singularité d’une histoire individuelle, mais une « normalité sociale » du suicide.

Norme, convention et aspects éthiques

Si on prend un point de vue évolutionniste, ce qui est « normal » est ce qui est le plus capable de s’adapter, de répondre aux besoins, de fournir un modèle pour la résolution des situations (la nature faisant disparaître l’anormalité), et la normalité tend vers la normativité, et il devient difficile de faire la distinction entre les deux aspects. En fait, David Hume aborde ce point dans l’exemple bien connu des rameurs, qui montent dans le même bateau pour traverser un fleuve et qui rament en cadence (cet exemple est longuement discuté dans Mackie (1980)). Les deux rameurs ajustent peu à peu leurs coups de rame, l’un par rapport à l’autre, et il n’est pas nécessaire d’obtenir un accord explicite (qui formulerait la « norme ») qu’ils respecteraient. La loi, qui consiste à imposer une norme peut être utile en cas de conflit (si un des rameurs refuse de ramer, ou deux rameurs de capacités physiques très différentes), mais bien souvent, il n’est pas nécessaire de formuler explicitement cette norme inhérente à leur conduite. L’observateur extérieur observera une régularité (lorsque le rythme de croisière sera atteinte) qu’il pourra modéliser, mais ce rythme « normal » observé n’est pas forcément imposé par une loi. Dans le cas de rameurs, on retrouve la notion d’équilibre évoquée précédemment. Construire un modèle, c’est extraire le signal du bruit (pour reprendre la distinction de Silver (2015)), c’est chercher une norme, au sens statistique. Mais cela va plus loin s’il l’on construit un modèle prédictif, la réalité devant alors se conformer au modèle, comme l’espèrent souvent les économètres (et autre modélisateurs).

# I Got The Feelin’

Last week, I’ve been going through my CD collection, trying to find records I haven’t been listing for a while. And I got the feeling that music I listen to nowadays is slower than the one I was listening to in my 20’s. I was wondering if that was an age issue, or it was simply the fact that music in the 90s was “faster” than the one released in 2015. So I tried to scrap the BPM database to get a more appropriate answer about this “feeling” I have. I did extract two information: the BPM (beat per minute) and the year (of release).

Here is a function to extract information from the website,

> library(XML)
> extractbpm = function(VBP,P){
+ url=paste("https://www.bpmdatabase.com/music/search/?artist=&title=&bpm=",VBP,"&genre=&page=",P,sep="")
+ return(tables)}

For instance

> extractbpm(115,13)
$track-table Artist Title 1 Eros Ramazzotti y Claudio Guidetti Dimelo A Mi 2 Everclear Volvo Driving Soccer Mom 3 Evils Toy Dear God 4 Expose In Walked Love 5 Fabolous ft. 2 Chainz When I Feel Like It 6 Fabolous ft. 2 Chainz When I Feel Like It 7 Fabolous ft. 2 Chainz When I Feel Like It 8 Fanny Lu Fanfarron 9 Featurecast Ain't My Style 10 Fem 2 Fem Obsession 11 Fernando Villalona Mi Delito 12 Fever Ray Triangle Walks 13 Firstlove Freaky 14 Fito Blanko Pegadito Suavecito 15 Flechazo Del Norte Mariposa Traicionera 16 Fluke Switch/Twitch 17 Flyleaf Something Better 18 FM Static The Next Big Thing 19 Fonseca Eres Mi Sueno 20 Fonseca ft. Maffio & Nayer Eres Mi Sueno 21 Francesca Battistelli Have Yourself A Merry Little Christmas 22 Frankie Ballard Young & Crazy 23 Frankie J. More Than Words 24 Frank Sinatra The Hucklebuck 25 Franz Ferdinand The Dark Of The Matinée Mix BPM Genre Label Year 1 — 115 — Sony 2009 2 — 115 — Capitol Records 2003 3 — 115 — — — 4 — 115 — Arista Records 1994 5 Explicit 115 Urban Def Jam/Island Def Jam 2013 6 — 115 Urban Def Jam/Island Def Jam 2013 7 Radio Edit 115 Urban Def Jam/Island Def Jam 2013 8 — 115 Latin Pop Universal Latino 2011 9 Psychemagik Dub 115 — Jalapeno 2012 10 — 115 — Critique Records 1993 11 — 115 — Mt&vi Records/caminante Records 2001 12 Rex The Dog Remix 115 — Little Idiot/Mute 2012 13 — 115 — Jwp Music 2000 14 — 115 Merengue Mambo Crown Loyalty 2012 15 — 115 — Hacienda 2010 16 Album Version 115 — One Little Indian Records 2004 17 — 115 Alternative A&M/Octone 2013 18 — 115 — Tooth & Nail Records 2007 19 — 115 Merengue Mambo 10 2012 20 Urban Version 115 — 10 2012 21 — 115 — Word/Fervent/Warner Bros 2009 22 — 115 Country Warner Bros 2015 23 Mynt Rocks Radio Edit 115 — Columbia 2005 24 — 115 Jazz Columbia 1950 25 — 115 New Wave — 2004 We have here one of the few old songs, a 1950 tune by Frank Sinatra. If we scrap the website, with a simple loop (where the bpm is from 40 to 200). Start with BASE=NULL > vbp=40 > p=1 and then, a loop based on > while(vbp<=2017){ + F=extractbmp(vbp,p) + if(length(F)==1){ + BASE=rbind(BASE,F[[1]][,c("Artist","Title","BPM","Year")]) + p=p+1} + if(length(F)==0){ + vbp=vbp+1 + p=1}} Then we should clean the dataset BASE=BASE[-duplicated(BASE),] BASE=BASE[-which(BASE$Year=="—"),]
BASE$y=as.numeric(as.character(BASE$Year))
BASE$bpm=as.numeric(as.character(BASE$BPM))
BASE=BASE[BASE$y>=1940,] and we end up with almost 50,000 tunes. boxplot(BASE$bpm~as.factor(BASE$y), col="light blue") Over the past 20 years, it looks like speed of tunes has declined (let us forget tunes of 2017, clearly we have a problem here…) library(mgcv) plot(BASE$y,BASE$bpm) reg=gam(bpm~s(y),data=BASE) B=data.frame(y=1950:2017) p=predict(reg,newdata=B) lines(B$y,p,lwd=3,col="red")

which is confirmed with a (smoothed) regression

p2=predict(reg,newdata=B,se.fit=TRUE)
plot(B$y,p2$fit,lwd=3,col="red",type="l",ylim=c(90,140))
lines(B$y,p2$fit+p2$se.fit) lines(B$y,p2$fit-p2$se.fit)

even when incorporating the confidence band. Bumps are probably related to smoothing parameters, but indeed, it looks like the average speed of music tune has decreased, from 110-115 in the 90’s to less than 100 nowadays. Now to be honest, I would love to have access to personal information from itunes, deezer or spotify, to get a better understanding (eg when in the week, in the day, do we like to listen to faster music for instance). But so far, I could not have access to such data. Too bad…

# Matching, Optimal Transport and Statistical Tests

To explain the “optimal transport” problem, we usually start with Gaspard Monge’s “Mémoire sur la théorie des déblais et des remblais“, where the the problem of transporting a given distribution of matter (a pile of sand for instance) into another (an excavation for instance). This problem is usually formulated using distributions, and we seek the “optimal” transport from one distribution to the other one. The formulation, in the context of distributions has been formulated in the 40’s by Leonid Kantorovich, e.g. from the distribution on the left to the distribution on the right.

Consider now the context of finite sets of points. We want to transport mass from points $\{A_1,\cdots,A_4\}$ to points $\{B_1,\cdots,B_4\}$. It is a complicated combinatorial problem. For 4 points, there are only 24 possible transfer to consider, but it exceeds 20 billions with 15 points (on each side). For instance, the following one is usually seen as inefficient

while the following is usually seen as much better

Of course, it depends on the cost of the transport, which depends on the distance between the origin and the destination. That cost is usually either linear or quadratic.

There are many application of optimal transport in economics, see eg Alfred’s book Optimal Transport Methods in Economics. And there are also applications in statistics, that what I’ve seen while I was discussing with Pierre while I was in Boston, in June. For instance if we want to test whether some sample were drawn from the same distribution,

set.seed(13)
npoints <- 25
mu1 <- c(1,1)
mu2 <- c(0,2)
Sigma1 <- diag(1, 2, 2)
Sigma2 <- diag(1, 2, 2)
Sigma2[2,1] <- Sigma2[1,2] <- -0.5
Sigma1 <- 0.4 * Sigma1
Sigma2 <- 0.4 *Sigma2
library(mnormt)
X1 <- rmnorm(npoints, mean = mu1, Sigma1)
X2 <- rmnorm(npoints, mean = mu2, Sigma2)
plot(X1[,1], X1[,2], ,col="blue")
points(X2[,1], X2[,2], col = "red")

Here we use a parametric model to generate our sample (as always), and we might think of a parametric test (testing whether mean and variance parameters of the two distributions are equal).

or we might prefer a nonparametric test. The idea Pierre mentioned was based on optimal transport. Consider some quadratic loss

ground_p <- 2
p <- 1
w1 <- rep(1/npoints, npoints)
w2 <- rep(1/npoints, npoints)
C <- cost_matrix_Lp(t(X1), t(X2), ground_p)
library(transport)
library(winference)
a <- transport(w1, w2, costm = C^p, method = "shortsimplex")

then it is possible to match points in the two samples

nonzero <- which(a$mass != 0) from_indices <- a$from[nonzero]
to_indices <- a$to[nonzero] for (i in from_indices){ segments(X1[from_indices[i],1], X1[from_indices[i],2], X2[to_indices[i], 1], X2[to_indices[i],2]) } Here we can observe two things. The total cost can be seen as rather large > cost=function(a,X1,X2){ nonzero <- which(a$mass != 0)
naa=a[nonzero,]
d=function(i) (X1[naa$from[i],1]-X2[naa$to[i],1])^2+(X1[naa$from[i],2]-X2[naa$to[i],2])^2
sum(Vectorize(d)(1:npoints))
}
> cost(a,X1,X2)
[1] 9.372472

and the angle of the transport direction is alway in the same direction (more or less)

> angle=function(a,X1,X2){
nonzero <- which(a$mass != 0) naa=a[nonzero,] d=function(i) (X1[naa$from[i],2]-X2[naa$to[i],2])/(X1[naa$from[i],1]-X2[naa$to[i],1]) atan(Vectorize(d)(1:npoints)) } > mean(angle(a,X1,X2)) [1] -0.3266797 > library(plotrix) > ag=(angle(a,X1,X2)/pi)*180 > ag[ag<0]=ag[ag<0]+360 > dag=hist(ag,breaks=seq(0,361,by=1)-.5) > polar.plot(dag$counts,seq(0,360,by=1),main=”Test Polar Plot”,lwd=3,line.col=4)

(actually, the following plot has been obtain by generating a thousand of sample of size 25)

In order to have a decent test, we need to see what happens under the null assumption (when drawing samples from the same distribution), see

Here is the optimal matching

Here is the distribution of the total cost, when drawing a thousand samples,

VC=rep(NA,1000)
VA=rep(NA,1000*npoints)
for(s in 1:1000){
X1a <- rmnorm(npoints, mean = mu1, Sigma1)
X1b <- rmnorm(npoints, mean = mu1, Sigma2)
ground_p <- 2
p <- 1
w1 <- rep(1/npoints, npoints)
w2 <- rep(1/npoints, npoints)
C <- cost_matrix_Lp(t(X1a), t(X1b), ground_p)
ab <- transport(w1, w2, costm = C^p, method = "shortsimplex")
VC[s]=cout(ab,X1a,X1b)
VA[s*npoints-(0:(npoints-1))]=angle(ab,X1a,X1b)
}
plot(density(VC)

So our cost of 9 obtained initially was not that high. Observe that when drawing from the same distribution, there is now no pattern in the optimal transport

ag=(VA/pi)*180
ag[ag<0]=ag[ag<0]+360
dag=hist(ag,breaks=seq(0,361,by=1)-.5)
polar.plot(dag$counts,seq(0,360,by=1),main="Test Polar Plot",lwd=3,line.col=4) Nice isn’t it? I guess I will spend some time next year working on those transport algorithm, since we have great R packages, and hundreds of applications in economics… # 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 (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 $\boldsymbol{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 monte 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é $p_i$ que l’évènement $A$ va survenir. Si $p_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 $\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]$

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

# Peut-on vraiment prévoir la probabilité de gagner une élection présidentielle ?

Cet article est co-signé avec Avner Bar-Hen, Conservatoire national des arts et métiers (CNAM)

L’élection du Président de la République française constitue un exemple important de mode de scrutin à deux tours. L’issue des deux tours est évidemment liée, et les résultats du premier permettent d’avoir une idée de l’issue du second tour.

Il est en effet plausible que les électeurs ayant voté au premier tour pour l’un des candidats qualifiés au second tour confortent leur choix. En revanche, le vote des électeurs ayant accordé au premier tour leur suffrage à des candidats éliminés du second tour est entaché de plus d’incertitude. Le choix de ces électeurs n’est pas totalement aléatoire car il est dépend fortement de la proximité idéologique et partisane entre le candidat pour lequel ils ont voté au premier tour et chacun des candidats qualifiés.

## Reports de vote

Les choix électoraux au second tour de scrutin des électeurs ayant voté au premier tour pour un candidat éliminé sont appelés les « reports de vote ». Ils sont généralement appréhendés de manière agrégée et sous forme de proportions : telle fraction de l’électorat de premier tour de tel candidat se reporte sur tel candidat au second tour.

La connaissance de ces reports de votes comporte principalement deux intérêts. D’une part, les reports de votes constituent un élément de science politique à part entière. Ils renseignent au sujet du positionnement respectif des électorats, de l’ordre de priorité qu’ils attribuent aux différentes thématiques du débat public, de leur fidélité partisane et de leur degré d’information politique. À ce titre, ils sont susceptibles d’intéresser les formations politiques et les candidats eux-mêmes dans l’élaboration ex ante de leurs stratégies électorales et dans leur évaluation ex post. D’autre part, connaître les reports de votes permet de construire des prédictions des résultats du second tour à l’issue du premier tour.

Les estimations de reports de votes publiées, notamment dans la presse, résultent le plus souvent de sondages effectués directement auprès des électeurs. L’estimation produite par le sondage va dépendre de l’échantillon tiré, avec pour corollaire l’existence d’un inévitable écart entre le paramètre et son estimation. Cet écart est quantifié en terme de probabilité associée à un sondage. Cependant les tailles d’échantillon utilisées et le très faible taux de réponse génèrent de telles incertitudes qu’il est difficile de conclure quoi que ce soit à partir de ces sondages.

Pour résumer rapidement, l’évolution (positive) de la présentation des sondages, et en particulier la présentation de termes de probabilité, doit s’accompagner d’une évolution des méthodes et en particulier des tailles d’échantillon. Il est important de rappeler que les intervalles de confiance n’ont aucune valeur prédictive (et donc aucune validité testable ex-post). Par exemple, comment lire le graphique ci-dessous sur les intentions de vote au premier tour :

Doit-on comprendre que fin janvier, 16 % des électeurs affirmaient vouloir voter pour Benoît Hamon, et qu’avec 95 chances sur 100 cette proportion était comprise entre 15 % et 17 % ? Un mois plus tard, il était à 14 %. A-t-il perdu des intentions de votes ? ou bien peut-on imaginer que des électeurs indécis fin janvier (et non comptabilisés) aient décidé de voter pour une autre alternative ? Les pourcentages ne sont pas, ici, relatifs à la même quantité.

Mais revenons à l’estimation des reports de voix, outils-clés pour prévoir la probabilité de gagner au second tour. Les estimations produites à l’aide d’outils mathématiques, à partir des résultats de scrutin par bureaux de vote, restent rares dans la presse mais c’est une source importante d’études en science politique. Essayons de voir le principe.

En démocratie, le vote de chaque électeur est confidentiel. À chaque tour d’une élection, seul le résultat agrégé des votes est connu, et le niveau le plus fin d’agrégation est donné par le bureau de vote. Les reports d’un vote à un autre sont inconnus, et leur estimation à partir des résultats agrégés nécessite de définir un modèle statistique. En effet, la décision individuelle de chaque électeur n’obéit à aucune fonction déterministe.

On résume souvent ce principe sous la formule « le vote est privé, les votes sont publics ». Néanmoins ce principe ouvre grand la porte au paradoxe de Simpson (on parle aussi d’erreur écologique) : à cause de la distribution hétérogène de l’échantillon, regrouper les données pointe une tendance qui peut être fausse, et qui disparaît si on analyse les données en séparant selon le facteur de confusion. L’utilisation de données socio-démographiques est donc importante.

## Estimer la probabilité qu’un candidat gagne une élection

Tout d’abord, notons que l’élection mobilise un corps électoral composé de N électeurs répartis dans K bureaux de vote, à raison de Nk électeurs dans le bureau k. Le nombre d’électeurs dans chaque bureau demeure identique lors des deux tours. Pour chaque tour et chaque électeur, le vote peut être vu comme la sélection d’une unique modalité dans un ensemble de choix commun à tous les électeurs. Pour l’élection présidentielle de 2017, il y a I=13 choix possibles au premier tour : s’abstenir, voter blanc ou pour un des onze candidats ; et J = 4 choix possibles au second : s’abstenir, voter blanc ou pour les candidats Macron ou Le Pen.

Étant donné l’hypothèse de stabilité du corps électoral, on peut construire pour chaque bureau k, un tableau croisé (on parle de tableau de contingence) constitué des nombres Nk(i,j) d’électeurs ayant opté pour le choix i au premier tour et j au second, à l’issu du dépouillement du second tour. De ces tableaux, seules les marges suivant les lignes Nk(i,) et suivant les colonnes Nk(,j) sont connues et correspondent respectivement aux résultats du premier et du second tour. Les conditions de vote (isoloir, bulletin sous enveloppe, dépouillement à la fin…) font que lors du scrutin, les électeurs sont globalement ignorants des choix des autres électeurs, à l’exception de l’abstention qui fait l’objet d’une estimation communiquée et révisée régulièrement en cours de journée. Par conséquent, il semble relativement raisonnable de considérer que les décisions de vote sont prises simultanément et qu’elles sont donc indépendantes entre électeurs.

Les probabilités de report du candidat i au candidat j diffèrent a priori entre électeurs d’un même bureau et entre bureaux. Même si le résultat final était connu, les différentes manières d’estimer ces probabilités soulèvent des problèmes complexes, résolus en mobilisant de nombreuses techniques algorithmiques. Estimer la probabilité qu’un électeur ayant fait un des 13 choix possibles (par exemple l’abstention) au premier tour en fasse un des 4 possibles au second (par exemple l’abstention) est un exercice compliqué mais indispensable si on souhaite vraiment estimer la probabilité de victoire d’un candidat. Et à part enfoncer des portes ouvertes, on ne peut pas dire grand-chose : si quelqu’un dit qu’une pièce est pipée, et à 80 % de chances de tomber sur « pile », le fait qu’elle tombe une fois sur « pile » ne signifiera pas qu’il a « raison », pas plus que tomber sur « face » ne signifiera qu’il a « tort ».

En statistique, c’est la répétition qui permet de savoir si l’estimation d’une probabilité est pertinente. Or c’est le propre d’une élection que de ne se produire qu’une fois.

Pour aller plus loin : Søren Risbjerg Thompsen, « Danish Elections 1920-79 : A logit Approach to Ecological Analysis and Inference » (1987) et Fendrich, L’Hour et Rateau, « Reports de votes entre les deux tours d’une élection présidentielle : estimation statistique et sociologie électorale » (2014).

Avner Bar-Hen, Professeur du Cnam, Conservatoire national des arts et métiers (CNAM) et Arthur Charpentier, Maître de Conférences, Université Rennes 1

La version originale de cet article a été publiée sur The Conversation.

# Proportion of people alive in 1945 that are still alive

In demography, we like to use life tables to estimate the probability that someone born in 1945 (say) is still alive nowadays.  But another interesting quantity might be the probability that someone alive in 1945 is still alive nowadays.

The main difference is that we do not know when that person, alive in 1945, was born. Someone who was old in 1945 is very unlikely still alive in 2017. To compute those probabilities, we can use datasets from http://www.mortality.org/hmd/. More precisely, we need both death and birth data. I assume that datasets (text files) were downloaded (it is necessary to register – for free – to get the data).

In the death dataset, there is a “110+” for people older than 110 years. For convenience, let us cap our observations at 110 years old,

D$Age=as.numeric(as.character(D$Age))
D$Age[is.na(D$Age)]=110

Consider now a first function that will return, for people born in 1930 (say) two informations

• the number of people (here, let us consider women only) born in 1930 (from the birth database)
• the number of death of people of age 0 in 1930, people of age 1 in 1931, people of age 2 in 1932, etc…

The code is simple

nb=function(y=1930){
debut=1816
MatDFemale=matrix(D$Female,nrow=111) colnames(MatDFemale)=debut+0:198 cly=y-debut+1:111 deces=diag(MatDFemale[,cly[cly%in%1:199]]) return(c(B$Female[B\$Year==y],deces))}

We have a single number for the number of births, and then a vector for the number of deaths. Consider now another function. Consider the people born in 1930. We want to get two numbers : the number of people still alive in 1945 (say), and the number of people still alive nowadays. The ratio will be the proportion of people born in 1930 that were alive in 1945, that are still alive in 2015.

pop=function(ne=1930,an=1945){
comptage=nb(ne)
s=0
if(an>ne) s=sum(comptage[seq(2,1+an-ne)])
p1=max(comptage[1]-s,0)
p2=max(p1-sum(comptage[seq(2+an-ne,length(comptage))]),0)
c(p1,p2)
}

Then, for a given year (say 1945), to get the proportion of people alive in 1945 that are still alive today, we need to count how many people born in 1944 were still alive in 1945, and in 2015, but also born in 1943, 1942, etc, And we simply consider the ratio of the total number of people alive in 2015 over the total number of people alive in 1945

ptn=function(y=1945){
V=Vectorize(function(x) pop(x,y))(1816:y)
sum(V[2,!is.na(V[2,])])/sum(V[1,!is.na(V[1,])])
}

Hence, 22% of those alive in 1945 are still alive in 2015,

> ptn(1945)
[1] 0.2209435

Actually, instead of looking only at 1945, it is possible to get a plot

P=Vectorize(ptn)(1900:2010)
plot(1900:2010,P,type="l",ylim=0:1)

For instance,

> ptn(1975)
[1] 0.6377413

i.e. 63.7% of those alive in 1975 are stil alive 40 years after. That is a rather interesting function, I was surprised that I couldn’t find it is standard demographical R package…

# Débat sur les prévisions

Il y a quelques semaines, la revue Risque organisait un débat sur la faillite des prévisions, auquel j’avais été convié, comme j’en parlais dans un précédant billet. La retranscription du débat va être publiée dans le prochain numéro de la revue.

There are now more than 10,000 R packages available from CRAN, much more if you include those available only on github. So, to be honest, it become difficult to know all of them. But sometimes, you discover a nice function in one of them, and that is really awesome. Consider for instance some (standard) censored lifetime data,

n=10000
idx=sample(1:4,size=n,replace=TRUE)
pd=LETTERS[idx]
lambda=1+(idx-1)/3
t=rexp(n,lambda)
x=rexp(n)
c=t>x
y=pmin(t,x)
df=data.frame(time=y,status=c,product=pd)

(yes, I will generate them here). Consider Kaplan-Meier estimator of the survival function,

library(survival)
km.base = survfit( Surv(time,status) ~ 1  , data = df )
plot(km.base)

This week end, Anat (currently finishing the Data Science for Actuaries program) made me discover a nice R function, to add information to that graph (well, not that graph, since it will be a ggplot version, but the same survival distribution plot)

library(ggplot2)
library(survminer)
ggsurvplot(km.base, main = "", color = "blue" , censor = FALSE, xlim = c(0,3), risk.table = TRUE ,
risk.table.col = "blue" , risk.table.height = 0.2, risk.table.title = "" , legend.labs = "All" , legend.title = "" , break.time.by = 1, xlab = "" , ylab = "")

This is more interesting when we have different lifetimes

km.prod = survfit( Surv(time,status) ~ product  , data = df )
ggsurvplot(km.prod, main = "", censor = FALSE, xlim = c(0,3), risk.table = TRUE , risk.table.col = "strata" , risk.table.height = 0.3, risk.table.title = "" , legend.labs = LETTERS[1:4] , legend.title = "" , break.time.by = 1, xlab = "" , ylab = "")

or, with a different time granularity

ggsurvplot(km.prod, main = "", censor = FALSE, xlim = c(0,3), risk.table = TRUE , risk.table.col = "strata" , risk.table.height = 0.3, risk.table.title = "" , legend.labs = LETTERS[1:4] , legend.title = "" , break.time.by = .5, xlab = "" , ylab = "")

Nice, isn’t it?

# Sondages, intentions de vote et “vote utile”

Il y a plusieurs semaines, Ashley Kirk et Patrick Scott publiaient French presidential election: How the polls are shaping up in the race to become president. Force est de constater qu’ils avaient raison: les sondages prennent une place essentielle dans cette période électorale.

Avant, une précaution d’usage s’impose, je pense. Les instituts de sondages le rappellent sans cesse, et ils ont raison : leur métier n’est pas véritablement de faire des sondages d’intention de vote. L’essentiel de leurs revenus vient des études de marché,
et d’opinion, généralement pour le compte d’entreprises privées. Néanmoins, les élections sont une fabuleuse vitrine publicitaire. On peut alors imaginer que l’institut qui aura “prédit” le scrutin avec la plus grande précision devrait s’assurer une rente pour les mois – voire les années – à venir. Cette vitrine publicitaire serait importante car elle empêcherait aux instituts de prendre parti : ils ont intérêt à prédire aussi précisément que possible, probablement plus que chercher à influencer les votes. Mais est-ce aussi simple ?

Depuis une quizaine d’années, on voit régulièrement ressurgir le concept de « vote utile ». Ce concept est particulièrement important en France, pour les élections présidentielles, compte tenu du système de vote : un vote en deux tours, avec seulement les deux meilleurs au second tour. Qui dit “meilleur” dit positionnement relatif, et donc importance du vote des autres. Quelque part, mon vote est important, mais celui des autres l’est tout autant (sinon plus). Pour voter de manière rationnelle, il faut savoir comment les autres vont voter. Et comme souvent dans un jeu statégique – et les élections en sont – l’information joue un rôle essentiel.

Lorsqu’un élécteur ne sait pas quel sera le vote des autres, il aura – a priori – intérêt à voter pour son candidat préféré, Mais si ce même électeur a des croyances quant au vote des autres, il peut avoir intérêt à modifier stratégiquement son choix, en votant pour son second, voire son troisième choix : on parlera alors de « vote utile ». Sans information sur les choix des autres (ou leurs stratégie), il ne peut pas y avoir de comportement stratégique. Et force est de constater que les “sondages” (terme générique pour désigner les estimations d’intension de vote) sont devenus des instruments essentiels pour déterminer justement quel sera le choix des autres électeurs.

Les sondages sont devenus l’élément central du vote stratégique, surtout s’ils sont présenté comme fournissant une information précise, voire “statistiquement rigoureuse” des intentions de vote. Les sondages deviennent alors auto-réalisateurs s’ils poussent à voter pour des candidats qui affichent le plus de chances de gagner, ce qu’on appelle communément le « vote utile ». Mais on peut aussi imaginer des stratégies inverses : si mon candidat préféré a 40% d’intention de votes, dans les derniers sondages, a-t-il vraiment besoin de ma voix ? Ne pourrais-je pas éviter de perdre du temps en allant voter ?

Bref, j’ai encore du mal à comprendre le rôle exact des sondages dans les élections, mais ils sont loin d’être neutres. Et j’ai l’impression que comprendre ce mécanisme est plus important pour quantifier les “marges d’erreurs” que toutes les méthodes statistiques que l’on peut imaginer (dont je parlais dans un précédent billet). Et les sondages en période électorale semblent relever davantage de la théorie des jeux que des problèmes d’inférence statistique.

# C’est normal ! (partie 1) et si la normalité n’existait pas?

Il y a quelques semaines, je discutais avec une collègue juriste (alors que j’essayais de glaner des statistiques juridiques), et alors que l’on parlait alors de la lenteur de l’instruction, ou de l’inversion de la charge de la preuve (je ne me souviens plus) j’ai été supris qu’elle me dise que “c’est normal“. Je sais qu’un juriste et un statisticien ne donnent pas forcément le même sens aux mots, et cette phrase m’a dérangé parce qu’une situation qui pourrait sembler “normale” (car observée régulièrement) n’est pas pour autant “juste” (c’est le paradoxe is/ought de Hume, mais j’en parlerais dans un autre billet).

Le point de départ est de comprendre ce qu’est la normalité “empirique”, telle qu’observée dans une population. Ce qui pourrait être “normal” pour un statisticien. Pour commencer, je pensais revenir sur un exemple raconté dans The End of Average de Todd Rose, qui essaye de montrer, exemples à l’appui, que l’homme moyen n’existe pas.

• L’homme moyen de Quételet

Au XIXème siècle, si plusieurs astronomes mesuraient la vitesse d’un même objet céleste, ils obtenaient (souvent) plusieurs mesures différentes. Pour savoir “laquelle” utiliser dans leurs calculs, l’idée d’utiliser “la méthode des moyennes” s’est rapidement imposée – comme le rappelle Stahl (2006), et surtout Sheynin (1973) – cette “moyenne” ayant une précision plus grande que n’importe quelle autre grandeur (ou dirait aujourd’hui “statistique”).

Adolphe Quetelet fut, semble-t-il, le premier à appliquer ce calcul de moyennes à des mesures humaines, introduisant son fameux concept d'”homme moyen”. Comme j’en avais parlé dans un précédant billet, la moyenne est une grande particulière, dont le sens n’est pas forcément clair. Si on définit la moyenne à l’aide d’une minimisation d’erreur quadratique, on a une interprétation en terme de prévision (on retrouve ici la notion de mesure elicitable dont je parlais dans mon dernier cours): la taille moyenne est la taille de “devrait” mesurer une personne tirée au hasard (à une variation aléatoire – et imprévisible – près). En 1846 dans une lettre (publiée dans l’ouvrage Lettres sur la théorie des probabilités, appliquée aux sciences morales) Adolphe Quételet utilise l’image de la statue du gladiateur pour expliquer ce que peut être l’homme moyen.

• L’interprétation de Francis Galton

Ce homme moyen a beaucoup plu à l’époque, en particulier au sein de l’école anglaise eugéniste, dirigée à l’époque par Francis Galton, même si ce dernier s’intéresse surtout aux déviations par rapport à cette norme (déviation vers le haut et déviation vers le bas). Comme le rappelle Bulmer (2004),the deviations from that average—upwards towards genius, and downwards towards stupidity—must follow the law that governs deviations from all true averages“. Les travaux de Galton ont visé à comprendre ces déviations. Si Florence Nightingale affirmait que “the Average Man is God’s Will“, Galton de son côté s’intéressait davantage au caractère héréditaire. Mais cet “homme moyen” a-t-il pour autant du sens ?

• L’être humain moyen n’existe pas

Une anecdote intéressante est celle de deux statues, à Cleveland, celles de Norma et Normann. L’artiste Abram Belskie et l’obstetricien Robert Latou Dickinson ont réalisé ensemble ces statues, en 1943. La particularité est qu’aucun modèle n’a été représenté. En fait, il s’agissait de représenter une femme et un homme qui avaient les mensurations moyennes de l’époque.

Une fois ces statues réalisées, un concours a été organisé pour trouver qui ces statues pouvaient bien représenter. Plusieurs milliers de personnes de l’Ohio ont envoyé leurs mensurations, mais aucun ne correspondaient à celles des statues. Certes, plusieurs centaines avaient la même taille. Plusieurs centaines avaient le même tour de poitrine. Mais aucune n’avait toutes les bonnes mesures. Car comme l’explique Todd Rose, l’homme n’est pas unidimensionnel: c’est sur plusieurs dimensions qu’on le mesure.

Et chercher à le résumer en une grandeur unidimensionnel est beaucoup trop réducteur. C’est ce qu’il montre dans son livre sur les tests d’intelligence, par exemple, où un même QI peut être associé à deux personnes très différentes.

Pareil pour décider de recruter quelqu’un, se focaliser sur un seul indicateur n’a pas de sens. Le soucis quand on travaille dans un contexte multivarié, c’est que la moyenne perd de son sens. Pour reprendre le titre d’un billet mis en ligne voilà 6 mois, être moyen peut être extraordinaire.

• La malédiction de la dimension (en grande dimension, l’espace est très vide…)

En fait, ce problème est bien connu par les statisticiens, sous le nom de “fléau de la dimension“. Prenons un exemple simple: supposons qu’une grandeur d’intérêt suive une loi normale $\mathcal{N}(\mu,\sigma^2)$, par exemple le poids, la taille, le tour de poitrine, etc. On pourrait dire que la norme, c’est se trouver dans un intervalle $[\mu\pm 1.5 \cdot\sigma]$. Si on a une loi normale, cette situation survient dans 85% des cas,

> (1-2*pnorm(-1.5))
[1] 0.8663856

Et 15% seront vues comme “anormales”. Elles peuvent être anormalement petites, ou anormalement grandes. C’est le dessin ci-dessous: on regarde ici seulement une dimension

On peut maintenant regarder deux dimensions, le poids et la taille, par exemple. La norme serait ici que dans les deux dimensions, on soit dans l’intervalle $[\mu\pm 1.5 \cdot\sigma]$. Si les grandeurs sont indépendantes, la probabilité que les deux grandeurs soient “normale” est de 75%

> (1-2*pnorm(-1.5))^2
[1] 0.750624

En dimension deux, 75% des observations sont “normales”, et 25% sont “anormales”

En dimension 3, on passe à 65%

> (1-2*pnorm(-1.5))^3
[1] 0.6503298

pour 35% d’observations “anormales” (plus du tiers)

Etc. En dimension cinq, on passe en dessous de 50%

> (1-2*pnorm(-1.5))^5
[1] 0.4881532

autrement dit, être dans la norme dans les 5 dimensions, ce n’est plus le cas de la majorité. Et en dimension vingt, ceux qui sont “normaux” sont plutôt atypiques, avec une proportion de l’ordre de 5%,

> (1-2*pnorm(-1.5))^20
[1] 0.0567838

Bref, la normalité est un concept particulièrement étrange sur le plan empirique, car intuitivement associé à l’idée d’une majorité. Alors que ce n’est pas le cas, la normalité étant justement atypique.