# 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 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… # 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é “surveiller 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 parlant 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 cœur 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")  0.158734 mean(B1$X1=="A")  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 &lt; 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 &lt; 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 ! # 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="") + download.file(url,destfile = "page.html") + tables=readHTMLTable("page.html") + 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[][,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… # Rupture Detection There are some graphs that you cannot forget. One graph that I found puzzling was mentioned on Andrew Gelman’s blog, a few years back, and was related to rupture detection What I remember from this graph is that if you want to get a rupture, you can easily find one… # Regression with Splines: Should we care about Non-Significant Components? Following the course of this morning, I got a very interesting question from a student of mine. The question was about having non-significant components in a splineregression. Should we consider a model with a small number of knots and all components significant, or one with a (much) larger number of knots, and a lot of knots non-significant? My initial intuition was to prefer the second alternative, like in autoregressive models in R. When we fit an AR(6) model, it’s not really a big deal if most coefficients are not significant (but the last one). It’s won’t affect much the forecast. So here, it might be the same. With a larger number of knots, we should be able to capture small bumps that we’ll never capture with a smaller number. Here is what a have with a small number of knots, and cubic splines and with a larger number of knots In order to understand what’s going on, consider a simple model, with the two splines above, in red > set.seed(1) > library(splines) > x=seq(0,1,by=.01) > v=bs(x,10) > x2=v[,2] > x10=v[,10] > set.seed(1) > y=1+3*x2+5*x10+rnorm(length(x))/4 > y_test=1+3*x2+5*x10+rnorm(length(x))/4 Note that here I have generated two sets of data, one to train a model, and one to test it. Here, the data looks like that > plot(x,y) It is based on two splines, > lines(df$x,1+3*x2+5*x10) If we use a spline model with 10 degrees of freedom, we get

> df=data.frame(x,y)
> reg=lm(y~bs(x,10),data=df)
> summary(reg)

Coefficients:
Estimate Std. Er t value Pr(>|t|)
(Intercept)  0.91671 0.17068   5.371 6.08e-07 ***
bs(x, 10)1   0.20485 0.32696   0.627    0.533
bs(x, 10)2   3.15593 0.22534  14.005  < 2e-16 ***
bs(x, 10)3   0.04847 0.25075   0.193    0.847
bs(x, 10)4   0.09373 0.21597   0.434    0.665
bs(x, 10)5   0.11624 0.22939   0.507    0.614
bs(x, 10)6   0.24829 0.22293   1.114    0.268
bs(x, 10)7  -0.06825 0.23498  -0.290    0.772
bs(x, 10)8   0.19633 0.26241   0.748    0.456
bs(x, 10)9   0.27557 0.26976   1.022    0.310
bs(x, 10)10  4.78134 0.24116  19.826  < 2e-16 ***

which makes sense, from what we have generated. Indeed, most of the components are not significant, but the second and the tenth. We can actually test that all those components are null (at the same time)

> A=matrix(0,8,11)
> colnames(A)=names(coefficients(reg))
> A[1,2]=A[2,4]=A[3,5]=A[4,6]=A[5,7]=
+ A[6,8]=A[7,9]=A[8,10]=1
> b=rep(0,8)
> linearHypothesis(reg, A,b)
Linear hypothesis test

Hypothesis:
bs(x, 10)1 = 0
bs(x, 10)3 = 0
bs(x, 10)4 = 0
bs(x, 10)5 = 0
bs(x, 10)6 = 0
bs(x, 10)7 = 0
bs(x, 10)8 = 0
bs(x, 10)9 = 0

Model 1: restricted model
Model 2: y ~ bs(x, 10)

Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     98 4.8766
2     90 4.6196  8   0.25701 0.6259  0.754

and yes, those coefficients are not significant.

> yp10=predict(reg)
> lines(df$x,yp10,col="red") # I Fought the (distribution) Law (and the Law did not win) A few days ago, I was asked if we should spend a lot of time to choose the distribution we use, in GLMs, for (actuarial) ratemaking. On that topic, I usually claim that the family is not the most important parameter in the regression model. Consider the following dataset > db <- data.frame(x=c(1,2,3,4,5),y=c(1,2,4,2,6)) > plot(db,xlim=c(0,6),ylim=c(-1,8),pch=19) To visualize a regression model, use the following code > nd=data.frame(x=seq(0,6,by=.1)) > add_predict = function(reg){ + prd1=predict(reg,newdata=nd,se.fit = TRUE,type="response") + y1=prd1$fit
+ y1_upp=prd1$fit+prd1$residual.scale*1.96*
prd1$se.fit + y1_low=prd1$fit-prd1$residual.scale*1.96* prd1$se.fit
+ polygon(c(nd$x,rev(nd$x)),c(y1_upp,
rev(y1_low)),col="light green",angle=90,
density=40,border=NA)
+ lines(nd$x,y1,col="red",lwd=2) + } For instance, with a Poisson regression (with a log link function) we get > plot(db) > reg1=glm(y~x,family=poisson(link="log"), + data=db) > add_predict(reg1) while, with a Gaussian regresion (but still with a log link function), we get > plot(db) > reg2=glm(y~x,family=gaussian(link="log"), + data=db) > add_predict(reg2) If we just care about the expected value of our prediction, the output is more or less the same > plot(db) > lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd, + type="response"),col="blue",lwd=1.5) So, indeed, forget about the (distribution) law when running a GLM. Not convinced? Consider – on the same dataset – a Poisson regression (with an identity link function this time) > plot(db) > reg1=glm(y~x,family=poisson(link="identity"), + data=db) > add_predict(reg1) while, with a Gaussian regresion (but still with an identity link function), we get > plot(db) > reg2=glm(y~x,family=gaussian(link="identity"), + data=db) > add_predict(reg2) Again, if we just plot the expected value of our prediction, the output is more or less the same > plot(db) > lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd, + type="response"),col="blue",lwd=1.5) So clearly, the simplistic message you should not care too much about the (distribution) law seems to be valid… # Classification with Categorical Variables (the fuzzy side) The Gaussian and the (log) Poisson regressions share a very interesting property, i.e. the average predicted value is the empirical mean of our sample. > mean(predict(lm(dist~speed,data=cars)))  42.98 > mean(cars$dist)
 42.98

One can prove that it is also the prediction for the average individual in our sample

> predict(lm(dist~speed,data=cars),
+ newdata=data.frame(speed=mean(cars$speed))) 42.98 The geometric interpretation is that the regression line passes through the centroid, > plot(cars) > abline(lm(dist~speed,data=cars),col="red") > abline(h=mean(cars$dist),col="blue")
> abline(v=mean(cars$speed),col="blue") > points(mean(cars$speed),mean(cars$dist)) But in all other cases, it is no longer the case. Consider for instance the case of a logistic regression. And to ask for something even more complicated, consider the case where we have only categorical explanatory variables. In that context, it is more difficult to get a prediction for the “average individual”. Unless we consider some fuzzy interpretation of the regression. # Regression Models, It’s Not Only About Interpretation Yesterday, I did upload a post where I tried to show that “standard” regression models where not performing bad. At least if you include splines (multivariate splines) to take into accound joint effects, and nonlinearities. So far, I do not discuss the possible high number of features (but with boostrap procedures, it is possible to assess something related to variable importance, that people from machine learning like). But my post was not complete: I was simply plotting the prediction obtained by some model. And it “looked like” the regression was nice, but so were the random forrest, the -nearest neighbour and boosting algorithm. What if we compare those models on new data? # On Some Alternatives to Regression Models When you start discussing with people in machine learning, you quickly hear something like “forget your econometric models, your GLMs, I can easily find a machine learning ‘model’ that can beat yours”. I am usually very sceptical, especially when I hear “easily” or “always“. I have no problem about the fact that I use old econometric models, but I had the feeling that things aren’t that easy. I can understand that we might have problems when we do have a lot of features (I am still working on that, I’ll get back to this point soon), but I have the feeling that I can still capture interactions, and non-linearities with standard econometric models as well as any machine learning algorithm. Just to illustrate, consider the following ‘model where is (just to illustrate) > n <- 5000 > rtf <- function(x1, x2) { sin(x1+x2)/(x1+x2) } > xgrid <- seq(1,6,length=31) > ygrid <- seq(1,6,length=31) > zgrid <- outer(xgrid,ygrid,rtf) > persp(xgrid,ygrid,zgrid,theta=30, phi=30, + col="green", ticktype="detailed",shade=TRUE) # Tukey and Mosteller’s Bulging Rule (and Ladder of Powers) When discussing transformations in regression models, I usually briefly introduce the Box-Cox transform (see e.g. an old post on that topic) and I also mention local regressions and nonparametric estimators (see e.g. another post). But while I was working on my ACT6420 course (on predictive modeling, which is a VEE for the SOA), I read something about a “Ladder of Powers Rule” also called “Tukey and Mosteller’s Bulging Rule“. To be honest, I never heard about this rule before. But that won’t be the first time I learn something while working on my notes for a course ! The point here is that, in a standard linear regression model, we have But sometimes, a linear relationship is not appropriate. One idea can be to transform the variable we would like to model, , and to consider This is what we usually do with the Box-Cox transform. Another idea can be to transform the explanatory variable, , and now, consider, For instance, this year in the course, we considered – at some point – a continuous piecewise linear functions, It is also possible to consider some polynomial regression. The “Tukey and Mosteller’s Bulging Rule” is based on the following figure. and the idea is that it might be interesting to transform and at the same time, using some power functions. To be more specific, we will consider some linear model for some (positive) parameters and . Depending on the shape of the regression function (the four curves mentioned on the graph above, in the four quadrant) different powers will be considered. To be more specific, let us generate different models, and let us look at the associate scatterplot, > fakedataMT=function(p=1,q=1,n=99,s=.1){ + set.seed(1) + X=seq(1/(n+1),1-1/(n+1),length=n) + Y=(5+2*X^p+rnorm(n,sd=s))^(1/q) + return(data.frame(x=X,y=Y))} > par(mfrow=c(2,2)) > plot(fakedataMT(p=.5,q=2),main="(p=1/2,q=2)") > plot(fakedataMT(p=3,q=-5),main="(p=3,q=-5)") > plot(fakedataMT(p=.5,q=-1),main="(p=1/2,q=-1)") > plot(fakedataMT(p=3,q=5),main="(p=3,q=5)") If we consider the South-West part of the graph, to get such a pattern, we can consider or more generally where and are both larger than 1. And the larger and/or , the more convex the regression curve. Let us visualize that double transformation on a dataset, say the cars dataset. > base=cars ­> names(base)=c("x","y") > MostellerTukey=function(p=1,q=1){ + regpq=lm(I(y^q)~I(x^p),data=base) + u=seq(min(min(base$x)-2,.1),max(base$x)+2,length=501) + par(mfrow=c(1,2)) + plot(base$x,base$y,xlab="X",ylab="Y",col="white") + vic=predict(regpq,newdata=data.frame(x=u),interval="prediction") + vic[vic<=0]=.1 + polygon(c(u,rev(u)),c(vic[,2],rev(vic[,3]))^(1/q),col="light blue",density=40,border=NA) + lines(u,vic[,2]^(1/q),col="blue") + lines(u,vic[,3]^(1/q),col="blue") + v=predict(regpq,newdata=data.frame(x=u))^(1/q) + lines(u,v,col="blue") + points(base$x,base$y) + + plot(base$x^p,base$y^q,xlab=paste("X^",p,sep=""),ylab=paste("Y^",q,sep=""),col="white") + polygon(c(u,rev(u))^p,c(vic[,2],rev(vic[,3])),col="light blue",density=40,border=NA) + lines(u^p,vic[,2],col="blue") + lines(u^p,vic[,3],col="blue") + abline(regpq,col="blue") + points(base$x^p,base$y^q) + } For instance, if we call > MostellerTukey(2,1) we get the following graph, On the left, we have the original dataset, and on the right, the transformed one, , with two possible transformations. Here, we did only consider the square of the speed of the car (and only one component was transformed, here). On that transformed dataset, we run a standard linear regression. We add, here, a confidence tube. And then, we consider the inverse transformation of the prediction. This line is plotted on the left. The problem is that it should not be considered as our optimal prediction, since it is clearly biased because . But quantiles associated with a monotone transformation are the transformed quantiles. So confidence tubes can still be considered as confidence tubes. Note that here, it could have be possible to consider another transformation, with the same shape, but quite different > MostellerTukey(1,.5) Of course, there is no reason to consider a simple power function, and the Box-Cox transform can also be used. The interesting point is that the logarithm can be obtained as a particular case. Furthermore, it is also possible to seek optimal transformations, seen here as a pair of parameters. Consider > p=.1 > bc=boxcox(y~I(x^p),data=base,lambda=seq(.1,3,by=.1))$y
> for(p in seq(.2,3,by=.1)) bc=cbind(bc,boxcox(y~I(x^p),data=base,lambda=seq(.1,3,by=.1))$y) > vp=boxcox(y~I(x^p),data=base,lambda=seq(.1,3,by=.1))$x
> vq=seq(.1,3,by=.1)
> library(RColorBrewer)
> blues=colorRampPalette(brewer.pal(9,"Blues"))(100)
> image(vp,vq,bc,col=blues)
> contour(vp,vq,bc,levels=seq(-60,-40,by=1),col="white",add=TRUE) The darker, the better (here the log-likelihood is considered). The optimal pair is here

> bc=function(a){p=a;q=a; as.numeric(-boxcox(y~I(x^p),data=base,lambda=q)$y)} > optim(c(1,1), bc,method="L-BFGS-B",lower=c(0,0),upper=c(3,3))$par
 0.5758362 0.3541601

$value  47.27395 and indeed, the model we get is not bad, Fun, ins’t it? # Modèles de prévision, fourre-tout Dans le dernier cours de modèle de prévision, la semaine passée, nous avions passé un peu de temps sur l’étude des points aberrants et des points influents. Tout est expliqué dans les slides (avec les codes) donc je ne reviendrais pas dessus. Je pourrais juste évoquer quelques lignes de codes utilisées pour voir l’impact d’une observation sur la régression, en enlevant l’observation de la base, et en regardant ce que ça donne sur la prévision > plot(cars) > text(cars[c(23,49),1],cars[c(23,49),2]+4,c(23,49)) > abline(lm(dist~speed,data=cars)) > abline(lm(dist~speed,data=cars[-23,]),col="red") > abline(lm(dist~speed,data=cars[-49,]),col="blue") On a ensuite vu les problèmes de choix de modèles, et les méthodes de pénalisation (pénalisation du R2, et de la log-vraisemblance). On a ensuite vu les méthodes automatiques de choix de modèle, > US=read.table("http://freakonometrics.free.fr/US.txt",sep=";") > reg=lm(Murder~.,data=US) > regs=step(reg) Start: AIC=63.01 Murder ~ Population + Income + Illiteracy + Life.Exp + HS.Grad + Frost + Area Df Sum of Sq RSS AIC - Income 1 0.236 128.27 61.105 - HS.Grad 1 0.973 129.01 61.392 <none> 128.03 63.013 - Area 1 7.514 135.55 63.865 - Illiteracy 1 8.299 136.33 64.154 - Frost 1 9.260 137.29 64.505 - Population 1 25.719 153.75 70.166 - Life.Exp 1 127.175 255.21 95.503 Step: AIC=61.11 Murder ~ Population + Illiteracy + Life.Exp + HS.Grad + Frost + Area Df Sum of Sq RSS AIC - HS.Grad 1 0.763 129.03 59.402 <none> 128.27 61.105 - Area 1 7.310 135.58 61.877 - Illiteracy 1 8.715 136.98 62.392 - Frost 1 9.345 137.61 62.621 - Population 1 27.142 155.41 68.702 - Life.Exp 1 127.500 255.77 93.613 Step: AIC=59.4 Murder ~ Population + Illiteracy + Life.Exp + Frost + Area Df Sum of Sq RSS AIC <none> 129.03 59.402 - Illiteracy 1 8.723 137.75 60.672 - Frost 1 11.030 140.06 61.503 - Area 1 15.937 144.97 63.225 - Population 1 26.415 155.45 66.714 - Life.Exp 1 140.391 269.42 94.213 Enfin, on a aussi vu le passage au logarithme sur la variable d’intérêt, avec un modèle log-normal, > plot(cars$speed,cars$dist) > regln=lm(log(dist)~speed,data=cars) > summary(regln) Call: lm(formula = log(dist) ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -1.46604 -0.20800 -0.01683 0.24080 1.01519 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.67612 0.19614 8.546 3.34e-11 *** speed 0.12077 0.01206 10.015 2.41e-13 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4463 on 48 degrees of freedom Multiple R-squared: 0.6763, Adjusted R-squared: 0.6696 F-statistic: 100.3 on 1 and 48 DF, p-value: 2.413e-13 > beta=coefficients(regln) > for(x in seq(2,30,by=.1)){ + zhat=beta+beta*x + yhat=exp(zhat+1/2*0.4463^2) + points(x,exp(zhat),col="red") + points(x,yhat,col="blue") + } > abline(lm(dist~speed,data=cars)) Pour comparer les modèles (linéaire versus log-linéaire), on avait suggéré de comparer les résidus des deux modèles, > U=cars$dist - exp(predict(regln)+1/2*0.4463^2)

> sd(U)
 17.12949
> sd(residuals(lm(dist~speed,data=cars)))
 15.22184

Ce soir, on fini cette première partie du cours, sur les données individuelles.

# Les Arbres de Classification

J’animerai une formation lundi 28 de 14:00 à 16:00 au local N-6320 de l’UQAM sur le thème introduction aux arbres de classification. Cette formation est organisée dans le cadre des séminaires en méthodes d’analyses quantitatives et qualitatives qui se tiennent régulièrement depuis un peu plus d’un mois. animé par le collectif pour le développement et les applications en mesure et évaluation (Cdame). Les slides sont disponibles en pdf (il y a quelques animations, qui ne passent qu’avec Acrobat)

La base utilisée tout au long des exposés est la suivante
> MYOCARDE=read.table("http://freakonometrics.free.fr/saporta.csv",head=TRUE,sep=";") # Elections Sans Bipartisme

Hier, sur Twitter, @JF_Godbout partageait un joli graphique relatif aux élections québécoises, avec les nombres de votes obtenus (ici en pourcentage des votes totaux) et le pourcentage de sièges que cela donne, Il faut dire qu’hier, c’était jour d’élection. Ce sont des élections à un tour, avec plusieurs partis (disons 4 grands partis si on se limite aux élections les plus récentes). Jusqu’à lundi, le PQ (Parti Québécois) avait une majorité (relative), et depuis lundi, le PLQ (Parti Libéral) est revenu au pouvoir. Et il y a à côté de ces deux grands partis d’autres partis, que l’on va regrouper. @JF_Godbout m’ayant envoyé ses données, j’ai pu reproduire les graphiques. La grande différence ici, entre nos approches est – je pense – la méthode de lissage: alors que  @JF_Godbout faisait des régressions locales, j’ai opté pour des splines.

reg=glm(Seat~bs(Election),data=d,family=binomial)
yp=predict(reg,newdata=data.frame(Election=1970:2015),type="response")

J’utilise des splines parce que ça se généralise facilement à d’autres lois, et qu’on peut aussi les utiliser quand on régresse sur plusieurs variables continues. Si on rajoute les résultats de l’élection de lundi, on obtient le graphique suivant Ce qui me gêne dans cette modélisation, c’est qu’on considère nos trois dernières régressions comme indépendantes…. alors que ce n’est pas le cas: on a trois modalités possibles, et la somme doit forcément faire 1. D’ailleurs, quand on somme les prédictions, on dépasse 1. Avec deux partis, notre régression logistique nous garantie d’avoir une somme à un, mais plus ici. A moins d’utiliser un modèle multinomial. Si on oublie le lissage, l’idée est d’écrire

qui est tout simplement l’extension du modèle logistique. La difficulté technique est que, sous R, quand on fait une régression logistique, la variable d’intérêt doit être une variable prenant trois facteurs, et non pas un vecteur de probabilités. Au lieu d’utiliser des taux, ou des pourcentages, je devrais utiliser des nombres d’observations de chaque modalités. On va donc crééer une fausse base

> vote=data.frame(
+ plq=plq$Vote, + pq=pq$Vote,
+ autre=autre$Vote, + Election=pq$Election)
> vote
plq     pq  autre Election
1  0.5465 0.3022 0.1513     1972
2  0.3378 0.4137 0.2485     1976
3  0.4608 0.4926 0.0466     1981
4  0.5599 0.3869 0.0532     1985
5  0.4995 0.4016 0.0989     1989
6  0.4440 0.4475 0.1085     1994
7  0.4355 0.4287 0.1358     1998
8  0.4599 0.3324 0.2077     2003
9  0.3308 0.2835 0.3857     2007
10 0.4208 0.3517 0.2275     2008
11 0.3120 0.3195 0.3685     2012
12 0.4140 0.2538 0.3312     2014
> M=t(round(vote[,1:3]*1000))
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
plq    546  338  461  560  500  444  436  460  331   421   312   414
pq     302  414  493  387  402  448  429  332  284   352   320   254
autre  151  249   47   53   99  108  136  208  386   227   368   331

On réplique alors autant de fois le nom du parti qu’on a de (pseudo) votes

> Y=X=NULL
> for(i in 1:3){
+ for(j in 1:ncol(M)){
+ Y=c(Y,rep(c("PLQ","PQ","Autre")[i],M[i,j]))
+ X=c(X,rep(pq$Election[j],M[i,j])) + }} > basev=data.frame(Y=as.factor(Y),X=X) C’est un peu artificiel, j’en conviens, mais ça marche… On peut alors faire la régression > library(nnet) > regv=multinom(Y~bs(X),data=basev) # weights: 15 (8 variable) initial value 13186.643301 iter 10 value 12241.931805 final value 12221.591172 converged et des prévisions > ypv=predict(regv,newdata= + data.frame(X=1973:2014),type="probs") > ypv Autre PLQ PQ 1 0.18426281 0.4705415 0.3451957 2 0.16051060 0.4767408 0.3627486 3 0.14149153 0.4805514 0.3779571 4 0.12631826 0.4826100 0.3910718 5 0.11426236 0.4834165 0.4023212 6 0.10474046 0.4833543 0.4119053 7 0.09729304 0.4827117 0.4199953 8 0.09156246 0.4817016 0.4267360 9 0.08727323 0.4804776 0.4322491 10 0.08421535 0.4791477 0.4366370 11 0.08223095 0.4777835 0.4399855 12 0.08120364 0.4764292 0.4423672 13 0.08105030 0.4751066 0.4438431 14 0.08171457 0.4738203 0.4444651 Et on notera que cette fois, nos prévisions somment à un. Cela dit, visuellement, on obtient des choses proches de ce qu’on avait auparavant, mais intellectuellement, je trouve plus logique d’utiliser autre chose que des régressions logistiques quand on est dans une élection où plusieurs partis sont en compétition. # Le Vote par Procuration en France La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“. Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition. • Modèle de régression, et recherche de variables explications Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que où est le nombre de procurations dans le bureau de vote , et où est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que , la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr) > bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";") > bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";") > bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";") > bv$BV=bv$BVCOM > baset1=merge(bt1,bv,by="BV") > baset2=merge(bt2,bv,by="BV") > baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT Si on suppose que est fonction le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote), > variable="PROPRIO" > reference="LOGEMENT" > baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference] il est légitime de tenter une régression logistique, voire un lissage par splines, si on pense que le lien peut ne pas être linéaire, Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques) > b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001) > pt1=predict(regt1,newdata=ND,se=TRUE,type="response") > pt2=predict(regt2,newdata=ND,se=TRUE,type="response") > library(RColorBrewer) > CL=brewer.pal(6, "RdBu") > plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration", + ylim=c(0,.1)) > for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0) + /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL)

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

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

(en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration). On peut faire la même chose sur le taux de logement HLM, dans le quartier, Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de  laisse à désirer (toute suggestion est la bienvenue !).

> summary(regt1l)

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

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

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

(Dispersion parameter for binomial family taken to be 1)

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

Number of Fisher Scoring iterations: 4

> summary(regt2l)

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

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

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

(Dispersion parameter for binomial family taken to be 1)

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

Number of Fisher Scoring iterations: 4

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

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

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

(j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre) Cette fois, on a une sortie de régression plus classique,

> summary(regt1l)

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

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

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

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

> summary(regt2l)

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

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

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

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

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

• Visualiser les bureaux de vote à Paris

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

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

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

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES)) > colnames(elec)=c("BV","PROCURATIONS","EXPRIMES") > elec=as.data.frame(elec) > elec$BV=bt1$BV Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur, > m=match(paris$BUREAU,elec$BV) > plotvar=100*elec$PROCURATIONS/elec$EXPRIMES > nclr=7 > plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] > class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1) > colcode=findColours(class, plotclr) Reste à conclure, en faisant une visualisation graphique de nos données > par(mar=c(1,1,1,1)) > plot(paris,col=colcode,border=colcode) > legend(656274.9, 6867308,legend=names(attr(colcode,"table")), + fill=attr(colcode, "palette"), cex=1, bty="n", + title="Frequence procurations (%)") Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives, > regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+ + I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1) Dans ce cas, la visualisation des résidus donne > m=match(paris$BUREAU,elec$BV) > plotvar=100*residuals(regt1) > nclr=7 > plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] > class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1) > colcode=findColours(class, plotclr) > par(mar=c(1,1,1,1)) > plot(paris,col=colcode,border=colcode) > legend(656274.9, 6867308,legend=names(attr(colcode,"table")), + fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus") Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux. # Precision, with Imprecise Words This morning, after my course on extreme values, some students did show me a question they got from practicals they were suppose to work on, with undergraduate students : To be more specific, they wanted some feedback about point B. Now, let’s make it clear : I have no idea what “precision” and “variation” could mean… But let’s try and see if we can get something usefull, that might help to understand the question. In order to illustrate, consider the following regression model, > plot(cars,pch=19,col="black",cex=.8) > abline(lm(dist~speed,data=cars),lty=2) Here is the summary table of the linear regression model > summary(lm(dist~speed,data=cars)) Estimate Std. Error t value Pr(>|t|) (Intercept) -17.5791 6.7584 -2.601 0.0123 * speed 3.9324 0.4155 9.464 1.49e-12 *** My first idea was that “variation of the X’s” should be related to the “variance” of the explanatory variable. But it is stupid. For instance, If we transform the explanatory variable, say with a multiplicative factor of 100, then the variance of X will be 10,000 times larger. And the regression will be the same > cars100=cars > cars100$speed=100*cars$speed > plot(cars100,pch=19,col="black",cex=.8) > abline(lm(dist~speed,data=cars100),lty=2) in the sense that > summary(lm(dist~speed,data=cars10)) Estimate Std. Error t value Pr(>|t|) (Intercept) -17.57909 6.75844 -2.601 0.0123 * speed 0.39324 0.04155 9.464 1.49e-12 *** And similarly, divide by 100. So, I guess using some affine transformation of the explanatory variable is clearly not the way we should get a variable with more “variability”. Let us try something else. And keep in mind the following quantities, > var(cars$speed)
 27.95918
> sd(cars$speed)/mean(cars$speed)
 0.3433535

with the variance, and the coefficient of variation. Consider the following modified dataset,

> carsg=cars
> carsg$speed=8 > carsg$speed=25
> carsg$speed=24 > carsg$speed=12

Four values were changed, here. Observe that, somehow, there is more variability

> var(carsg$speed)  31.84694 > sd(carsg$speed)/mean(carsg$speed)  0.3640845 But if we consider the output of the regression model, we get > summary(lm(dist~speed,data=carsg)) Estimate Std. Error t value Pr(>|t|) (Intercept) -18.5681 5.3621 -3.463 0.00113 ** speed 3.9708 0.3254 12.201 2.55e-16 *** It look like we got here a more precision on the slope, with a smaller variance, and a larger Student-t-value. But what if we consider the following transformation, > carsg=cars > carsg$speed=5
> carsg$speed=25 > carsg$speed=25
> carsg$speed=7 Again, we have more variability here, on the explanatory variable, > var(carsg$speed)
 32.9898
> sd(carsg$speed)/mean(carsg$speed)
 0.3754036

But this time,

> summary(lm(dist~speed,data=carsg))

Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.5078     8.0498  -0.187    0.852
speed         2.9077     0.4932   5.896 3.61e-07 ***

the estimator of the slope has more variance, and we have a smaller Student-t-value. So here, if we increase the “variability” of X, we get get… almost anything. The intuition about those two transformations is relatively simple. In the first case, I have moved observations that were far away from the regression line – but in the center of the distribution, and I put them closer to the regression line, but more on the border of the sample (to increase the variance) (I would not call them outliers since outliers are defined as observations far away from the model, but on Y, not on X). In the second case, I did exactly the opposite. I am not sure if I understood correctly this sentence. But it looks like it is incorrect. Since there is only one false statement here, I will go for this one. What do you think?