# Histogramme et densité en échelle logarithmique

Il y a presque 20 ans, Paul-André Rosental, Gilles Postel-Vinay, Akiko Suwa-Eisenmann et Jérôme Bourdieu publiaient Migrations et transmissions inter-générationnelles dans la France
du XIXe et du début du XXe siècle qui présentait le graphique ci-dessous,

On est tombé sur ce graphique avec Ewen Gallic en écrivant notre article Using Collaborative Genealogy Data to Study Migration (qui paraîtra très bientôt dans le journal The History of the Family). Classiquement dans un article académique, on a besoin de se positionner par rapport à la littérature existante. Avec nos propres données, et nous avions construit l’histogramme suivant

histoire de montrer que nous retrouvions cette forme bimodale. Depuis le début, je suis pas très confortable avec ce graphique, mais on l’a laissé non pas pour ses vertus sémiologiques, mais parce qu’il valide notre approche, en obtenant des résultats comparables à d’autres, déjà établis dans la littérature.

Je repensais à tout ça avant hier, en mentionnant sur twitter le graphique suivant

qui contient quelque chose de similaire, à savoir une espèce d’histogramme en bleu, sauf que la largeur des bandes décroit ici de manière logarithmique…

Je suis très mal à l’aise face à ces graphiques, parce que je ne sais pas comment les interpréter.  Mais peut-être faut-il revenir à la base, pour comprendre ce qu’on fait.

On a un échantillon $\{x_1,\cdots,x_n\}$ disons de montant d’impôt payé. On a un échantillon comme on dit en statistique descriptive. En statistique mathématique, on suppose que les $x_i$ sont des réalisations de variables aléatoires. On pourrait ainsi dire que $x_i=X(\omega_i)$$X$ est une variable aléatoire définie sur un espace probabilisé $(\Omega,\mathcal{A},\mathbb{P})$. Pour rappel, $\Omega$, c’est l’univers, une espèce d’espace fondamental abstrait. Sur le site les-mathématiques, on nous dit que $\Omega$ est l’espace des observables, ce qui me déplait assez car pour moi, justement $\Omega$ est un espace abstrait. En théorie de la décision, si on reprend par exemple Probability and Uncertainty in Economic Modeling, on parle d’états de la nature (c’est aussi la terminologie que l’on retrouve dans States of Nature and the Nature of States). Et assez souvent, on oublie cet espace, grâce au théorème de transfert : en effet, si on a une variable aléatoire réelle, $X:\Omega \rightarrow \mathbb {R}$, alors $${\displaystyle \mathbb {E} \left[\varphi (X)\right]=\int _{\Omega }\varphi {\big (}X(\omega ){\big )}\mathbb {P} (\mathrm {d} \omega )=\int _{\mathbb {R} }\varphi (x)\mathbb {P} _{X}(\mathrm {d} x)}$$ce qui veut dire, de manière polie, qu’on ne se place jamais sur ces espace fondamental, mais on regarde juste leur “transfert” sur la droite réelle (ici, je vais me contenter du cas unidimensionnel), c’est à dire la “valeur” prise par la variable aléatoire (je reviendrai tout à l’heure sur ces écritures sous forme intégrales). C’est un abus de langage que l’on fait quand on dit que pour un lancer de dé, l’univers est $\Omega = \{1, 2, 3, 4, 5, 6\}$, correspondant aux valeurs des faces du dé. En fait, l’espace fondamental, des “états de la nature” peut être plus compliqué que ça, mais pour mieux comprendre, on va le “transférer” sur l’espace comptable $\{1, 2, 3, 4, 5, 6\}$. Mais je commence à m’écarter du sujet, surtout que je m’étais promis que je ne ferais pas (trop) de théorie de la mesure. L’autre raison est que je pense que ce n’est pas exactement la formalisation qu’on utilise en statistique. Je pense que $x_i=X_i(\omega)$ où les variables $X_1,\cdots,X_n$ sont des variables aléatoires indépendantes, et de même loi. Pourquoi ? Avec les notations précédentes, on dispose de l’échantillon suivant

que l’on peut voir comme une réalisation des variables aléatoires suivantes

Une statistique sera alors une fonction construit sur notre échantillon, $\hat{\theta}=h(x_1,\cdots,x_n)$ , par exemple la moyenne (empirique)

qui est alors un nombre. Mais ça peut être aussi une variable aléatoire, $\hat{\theta}=h(X_1,\cdots,X_n)$ , soit ici, pour la moyenne

La principale difficulté est qu’en statistique mathématique, on utilise la même notation, $\widehat{\theta}$ pour deux objets de nature très différentes. Et donc avoir des abus de langage, en pouvait dire en cours “ici la moyenne vaut 37.5%” et ensuite demander de “calculer la variance de la moyenne”. Mais là encore, je m’écarte. Revenons à notre variable aléatoire $X$. Souvent, on va essayer de décrire sa loi de probabilité.

Le plus simple est de passer par la fonction de répartition, définie sans ambiguïté par $F(x) = \mathbb{P}[\{\omega\in\Omega:X(\omega)\leq x\}]$ (ce qui a du sens cas la mesure de probabilité $\mathbb{P}$ est effectivement définie sur l’espace des états de la nature $\Omega$. Mais par abus de langage, on se contentera de noter $F(x) = \mathbb{P}[X\leq x]$ . Supposons maintenant que notre variable $X$ est absolument continue. Alors dans ce cas – ça correspond je pense au premier théorème de l’analyse$F(x)=\int _{-\infty}^{x}f(t) dt$, où $f$ est alors la densité de la variable aléatoire. Le point un peu subtile ici est que l’intégrale est ici définie par rapport à la mesure de Lebesgue (ce qui permet de définir proprement ce qu’est ce $dt$ à la fin de l’intégrale). Mais oublions ce point un instant.

On arrive enfin au premier point auquel je voulais arriver : pour les variables (absolument) continue, on peut “représenter” leur loi par leur densité. C’est le dessin ci-dessus, si $X$ est une variable qui suit une loi lognormale

u=seq(0,1000,length=251) v=dlnorm(u,5,1) plot(u,v,type="l",xlab="",ylab="")

L’interprétation est simple : comme la probabilité $\mathbb {P} (a se calcule alors par la relation suivante : $$\mathbb {P} \left(ala probabilité $\mathbb {P} (a se lit comme l’aire sous la courbe sur l’intervalle $[a, b]$.

Parmi les estimateurs usuels de la densité, on peut utiliser un histogramme. Même si c’est un objet que tout le monde manipule depuis ses cours au secondaire, la construction formelle de cet objet est un peu technique. En effet, le premier point est qu’il faut découper l’ensemble des valeurs prises par les observations $x_i$ – disons $(a,b]$ – en une partition, c’est à dire un ensemble d’intervalles disjoints qui vont recouvrir  $(a,b]$, $I_1,\cdots,I_k$. Le plus simple est d’avoir une partition régulière, $I_1=(a,a+(b-a)/k]$, $I_2=(a+(b-a)/k,a+2(b-a)/k]$, etc. Généralement, $I_j=(a+(b-a)(j-1)/k,a+(b-a)j/k]$. Et classiquement, on construit un histogramme en comptant le nombre d’observations dans chacun des intervalles,

x=rlnorm(500,5,1) hist(x,xlim=c(0,1000),breaks=seq(0,10000,by=100))

C’est l’histogramme classique, avec un pas de 100.

hist(x,xlim=c(0,1000),breaks=seq(0,10000,by=100),probability=TRUE,ylim=c(0,.005)) lines(u,v,col="red")

Pour avoir une densité (i.e. une “fonction en escalier qui va s’intégrer à 1”) il faut normaliser la hauteur (en divisant par le pas). On peut d’ailleurs superposer la densité de la loi lognormale.

Dans les graphiques que je présentais, on utilise en échelle logarithmique en abscisse. On peut le faire avec R en demandant

plot(u,v,log="x",type="l")

ce qui correspond aux graphiques précédant. On peut aussi, à la main, dire qu’on veut non plus voir $\{x,f(x)\}$ mais $\{\log(x),f(x)\}$.

plot(log(u),v,type="l")

L’aire précédente devient alors

Visuellement, on raisonne par comparaison. Quand on dit que la densité est une aire sous la courbe, en réalité, c’est la proportion de l’aire sous la courbe (par rapport à l’aire totale) qui nous parle. Mais autant le faire avec deux aires, pour mieux comprendre. A droite, en rouge, on a la probabilité d’avoir une valeur entre 200 et 600, soit ici environ 30%.

plnorm(600,5,1)-plnorm(200,5,1) [1] 0.3015131

et en bleu, on a la probabilité d’être entre 92 et 200, qui vaut là aussi environ 30%

plnorm(200,5,1)-plnorm(92,5,1) [1] 0.3010197

Ici les deux aires sont égales. On notera que ce n’est pas forcément évident, au premier coup d’œil, et l’asymétrie (vers la droite) laisse à croire que l’aire rouge est peut-être plus grande. Mais passons.

Considérons maintenant notre échelle logarithmique, en abscisse,

ou alors la version sur le logarithme des abscisses,

Comme on le notait auparavant, on aime ce graphique car on reconnait une forme de densité de “loi normale” (ce qui pourrait faire du sens, car on obtient une loi lognormale justement en prenant l’exponentielle d’une loi normale…). Cela dit, on a clairement un problème d’échelle ici, car l’aire totale sous la courbe vaut moins de 1%… autrement dit, on est loin de 1. Mais on peut quand même retrouver une loi normale,

lines(u2,dnorm(u2,4,1)/90,col="blue")

(on pourra noter d’ailleurs que non seulement, il a fallu diviser par 90, mais la moyenne est ici 4, et pas 5). Comme on a la densité d’une loi normale (à une transformation près), on peut facilement faire des calculs, en particulier, calculer les deux aires, rouges et bleues

range(log(c(u[I],rev(u[I])))) [1] 5.300315 6.396263 pnorm(range(log(c(u[I],rev(u[I])))),4,1) [1] 0.9032535 0.9917184 diff(pnorm(range(log(c(u[I],rev(u[I])))),4,1)) [1] 0.08846485 diff(pnorm(range(log(c(u[I2],rev(u[I2])))),4,1)) [1] 0.2019666

Aussi, l’aire bleue vaut ici plus du double de l’aire rouge. Autrement dit, si on interprète la courbe comme une densité sur le graphique ci-dessous

on a intuitivement envie de dire qu’il y a 2 fois plus de chance d’avoir une valeur entre 92 et 200 qu’entre 200 et 600. Ce qui est faux.

Si on essaye de comprendre un peu mieux, je vais revenir sur deux points… La première c’est que, aussi étrange que ça puisse paraître, c’est presque un coup de chance qu’en prenant une échelle logarithmique, on obtient une densité lognormale. On avait écrit un article avec Emmanuel Flachaire (Log-Transform Kernel Density Estimation of Income Distribution) qui revenait sur l’importance de cette transformation logarithmique. Mais ici, c’est un peu différent. On a en effet deux graphiques : le premier c’est $\{x,f(x)\}$ alors le second c’est $\{\log(x),f(x)\}$, soit, en faisant un changement de variable $y=\log(x)$ (ou $x=e^y$), c’est $\{y,f(e^y)\}$. Or dans le premier cas, on avait une loi lognormale, autrement dit $$f(x )=\frac {1}{x\sigma {\sqrt {2\pi }}}\exp \left(-\frac {(\log x-\mu )^{2}}{2\sigma ^{2}}\right)={\frac {1}{x}}\phi(\log(x);\mu ,\sigma^2 )$$$\phi$ est la densité de la loi normale centrée réduite. Si on regarde le second graphique, $\{y,f(e^y)\}$, on représente alors $\{y,e^{-y}\phi(y;\mu ,\sigma^2 )\}$. La magie va opérer parce qu’on peut rentrer le $e^{-y}$ dans la densité de la loi normale, ce qui donnera deux choses (1) une translation de la moyenne, (2) un facteur multiplicatif sur la densité. Ce sont les deux phénomème que nous avions observé ici. Si on détaille un peu $$e^{-y}\exp \left(-\frac {(y-\mu )^{2}}{2\sigma ^{2}}\right)=\exp \left(-\frac {(y-[\mu-\sigma^2] )^{2}}{2\sigma ^{2}}+\star\right)$$où on retrouve la translation de 1 de moyenne (centrée sur 4 et non plus sur 5, mais c’est logique puisqu’on avait pris $\sigma^2$ ayant pour valeur 1. Je laisse les plus courageux calculer le $\star$ qui va donner le facteur multiplicatif. Le second point est un peu plus technique, il est lié à un problème de mesure. Dans le premier cas, quand on faisait un calcul d’intégrale, on avait un $dx$ et dans le second cas, on calcule des aires avec un $dy$, mais ce n’est pas la bonne transformation. En effet, si $y=\log(x)$ (ou $x=e^x$) alors $dx=e^ydy$. Plus formellement, dans le premier cas, quand on calculait $\mathbb{P}[X\in[a,b]]$ on calculait l’intégrale $\int _{a}^{b}f(x) dx$. Dans le second cas, on calcule $\int _{\alpha}^{\beta}f(e^y) dy$ puisqu’on visualise la courbe $\{y,f(e^y)\}$. Mais si on fait un changement de variable propre $$\int _{\alpha}^{\beta}f(e^y) dy=\int _{a}^{b}f(x) \frac{dx}{x}$$qui n’est plus du tout l’intégrale qu’on cherche à calculer. Je pense que cette histoire de $dy$ correspondant à $x^{-1}dx$ pourrait avoir une interprétation en terme de changement de mesure (la mesure de référence n’est plus la mesure de Lebesgue, qui nous donne une relative uniformité sur l’axe des abscisses, et permet de faire un lien avec l’intégrale “classique” – au sens de Rieman).

Bref, cette histoire de mettre une échelle logarithmique en abscisse pour visualiser une densité est très perturbant. L’intuition que l’on peut en avoir est biaisée, et l’objet mathématique que l’on créé est complexe…

# The “probability to win” is hard to estimate…

Real-time computation (or estimation) of the “probability to win” is difficult. We’ve seem that in soccer games, in elections… but actually, as a professor, I see that frequently when I grade my students.

Consider a classical multiple choice exam. After each question, imagine that you try to compute the probability that the student will pass. Consider here the case where we have 50 questions. Students pass when they have 25 correct answers, or more. Just for simulations, I will assume that students just flip a coin at each question… I have $n$ students, and 50 questions

set.seed(1) n=10 M=matrix(sample(0:1,size=n*50,replace=TRUE),50,n)

Let $X_{i,j}$ denote the score of student $i$ at question $j$. Let $S_{i,j}$ denote the cumulated score, i.e. $S_{i,j}=X_{i,1}+\cdots+X_{i,j}$. At step $j$, I can get some sort of prediction of the final score, using $\hat{T}_{i,j}=50\times S_{i,j}/j$. Here is the code

SM=apply(M,2,cumsum) NB=SM*50/(1:50)

We can actually plot it

plot(NB[,1],type="s",ylim=c(0,50)) abline(h=25,col="blue") for(i in 2:n) lines(NB[,i],type="s",col="light blue") lines(NB[,3],type="s",col="red")

But that’s simply the prediction of the final score, at each step. That’s not the computation of the probability to pass !

Let’s try to see how we can do it… If after $j$ questions, the students has 25 correct answer, the probability should be 1 – i.e. if $S_{i,j}\geq 25$ – since he cannot fail. Another simple case is the following : if after $j$ questions, the number of points he can get with all correct answers until the end is not sufficient, he will fail. That means if $S_{i,j}+(50-i+1)< 25$ the probability should be 0. Otherwise, to compute the probability to sucess, it is quite straightforward. It is the probability to obtain at least $25-S_{i,j}$ correct answers, out of $50-j$ questions, when the probability of success is actually $S_{i,j}/j$. We recognize the survival probability of a binomial distribution. The code is then simply

PB=NB*NA for(i in 1:50){ for(j in 1:n){ if(SM[i,j]&gt;=25) PB[i,j]=1 if(SM[i,j]+(50-i+1)&lt;25) PB[i,j]=0 if((SM[i,j]&lt;25)&amp;(SM[i,j]+(50-i+1)&gt;=25)) PB[i,j]=1-pbinom(25-SM[i,j],size=(50-i),prob=SM[i,j]/i) }}

So if we plot it, we get

plot(PB[,1],type="s",ylim=c(0,1)) abline(h=25,col="red") for(i in 2:n) lines(PB[,i],type="s",col="light blue") lines(PB[,3],type="s",col="red")

which is much more volatile than the previous curves we obtained ! So yes, computing the “probability to win” is a complicated exercice ! Don’t blame those who try to find it hard to do !

Of course, things are slightly different if my students don’t flip a coin… this is what we obtain if half of the students are good (2/3 probability to get a question correct) and half is not good (1/3 chance),

If we look at the probability to pass, we usually do not have to wait until the end (the 50 questions) to know who passed and who failed

PS : I guess a less volatile solution can be obtained with a Bayesian approach… if I find some spare time this week, I will try to code it…

# Solving the chinese postman problem

Some pre-Halloween post today. It started actually while I was in Barcelona : kids wanted to go back to some store we’ve seen the first day, in the gothic part, and I could not remember where it was. And I said to myself that would be quite long to do all the street of the neighborhood. And I discovered that it was actually an old problem. In 1962, Meigu Guan was interested in a postman delivering mail to a number of streets such that the total distance walked by the postman was as short as possible. How could the postman ensure that the distance walked was a minimum?

A very close notion is the concept of traversable graph, which is one that can be drawn without taking a pen from the paper and without retracing the same edge. In such a case the graph is said to have an Eulerian trail (yes, from Euler’s bridges problem). An Eulerian trail uses all the edges of a graph. For a graph to be Eulerian all the vertices must be of even order.

An algorithm for finding an optimal Chinese postman route is:

1. List all odd vertices.
2. List all possible pairings of odd vertices.
3. For each pairing find the edges that connect the vertices with the minimum weight.
4. Find the pairings such that the sum of the weights is minimised.
5. On the original graph add the edges that have been found in Step 4.
6. The length of an optimal Chinese postman route is the sum of all the edges added to the total found in Step 4.
7. A route corresponding to this minimum weight can then be easily found.

For the first steps, we can use the codes from Hurley & Oldford’s Eulerian tour algorithms for data visualization and the PairViz package. First, we have to load some R packages

require(igraph) require(graph) require(eulerian) require(GA)

Then use the following function from stackoverflow,

make_eulerian = function(graph){ info = c("broken" = FALSE, "Added" = 0, "Successfull" = TRUE) is.even = function(x){ x %% 2 == 0 } search.for.even.neighbor = !is.even(sum(!is.even(degree(graph)))) for(i in V(graph)){ set.j = NULL uneven.neighbors = !is.even(degree(graph, neighbors(graph,i))) if(!is.even(degree(graph,i))){ if(sum(uneven.neighbors) == 0){ if(sum(!is.even(degree(graph))) &gt; 0){ info["Broken"] = TRUE uneven.candidates &lt;- !is.even(degree(graph, V(graph))) if(sum(uneven.candidates) != 0){ set.j &lt;- V(graph)[uneven.candidates][[1]] }else{ info["Successfull"] &lt;- FALSE } } }else{ set.j &lt;- neighbors(graph, i)[uneven.neighbors][[1]] } }else if(search.for.even.neighbor == TRUE &amp; is.null(set.j)){ info["Added"] &lt;- info["Added"] + 1 set.j &lt;- neighbors(graph, i)[ !uneven.neighbors ][[1]] if(!is.null(set.j)){search.for.even.neighbor &lt;- FALSE} } if(!is.null(set.j)){ if(i != set.j){ graph &lt;- add_edges(graph, edges=c(i, set.j)) info["Added"] &lt;- info["Added"] + 1 } } } (list("graph" = graph, "info" = info))}

Then, consider some network, with 12 nodes

g1 = graph(c(1,2, 1,3, 2,4, 2,5, 1,5, 3,5, 4,7, 5,7, 5,8, 3,6, 6,8, 6,9, 9,11, 8,11, 8,10, 8,12, 7,10, 10,12, 11,12), directed = FALSE)

To plot that network, use

V(g1)$name=LETTERS[1:12] V(g1)$color=rgb(0,0,1,.4) ly=layout.kamada.kawai(g1) plot(g1,vertex.color=V(newg)$color,layout=ly) Then we convert it to some traversable graph by adding 5 vertices eulerian = make_eulerian(g1) eulerian$info broken Added Successfull 0 5 1 g = eulerian$graph as shown below ly=layout.kamada.kawai(g) plot(g,vertex.color=V(newg)$color,layout=ly)

We cut those 5 vertices in two part, and therefore, we add 5 artificial nodes

A=as.matrix(as_adj(g)) A1=as.matrix(as_adj(g1)) newA=lower.tri(A, diag = FALSE)*A1+upper.tri(A, diag = FALSE)*A for(i in 1:sum(newA==2)) newA = cbind(newA,0) for(i in 1:sum(newA==2)) newA = rbind(newA,0) s=nrow(A) for(i in 1:nrow(A)){ Aj=which(newA[i,]==2) if(!is.null(Aj)){ for(j in Aj){ newA[i,s+1]=newA[s+1,i]=1 newA[j,s+1]=newA[s+1,j]=1 newA[i,j]=1 s=s+1 }}}

We get the following graph, where all nodes have an even number of vertices !

newg=graph_from_adjacency_matrix(newA) newg=as.undirected(newg) V(newg)$name=LETTERS[1:17] V(newg)$color=c(rep(rgb(0,0,1,.4),12),rep(rgb(1,0,0,.4),5)) ly2=ly transl=cbind(c(0,0,0,.2,0),c(.2,-.2,-.2,0,-.2)) for(i in 13:17){ j=which(newA[i,]&gt;0) lc=ly[j,] ly2=rbind(ly2,apply(lc,2,mean)+transl[i-12,]) } plot(newg,layout=ly2)

Our network is now the following (new nodes are small because actually, they don’t really matter, it’s just for computational reasons)

plot(newg,vertex.color=V(newg)$color,layout=ly2, vertex.size=c(rep(20,12),rep(0,5)), vertex.label.cex=c(rep(1,12),rep(.1,5))) Now we can get the optimal path n &lt;- LETTERS[1:nrow(newA)] g_2 &lt;- new("graphNEL",nodes=n) for(i in 1:nrow(newA)){ for(j in which(newA[i,]&gt;0)){ g_2 &lt;- addEdge(n[i],n[j],g_2,1) }} etour(g_2,weighted=FALSE) [1] "A" "B" "D" "G" "E" "A" "C" "E" "H" "F" "I" "K" "H" "J" "G" "P" "J" "L" "K" "Q" "L" "H" "O" "F" "C" [26] "N" "E" "B" "M" "A" or edg=attr(E(newg), "vnames") ET=etour(g_2,weighted=FALSE) parcours=trajet=rep(NA,length(ET)-1) for(i in 1:length(parcours)){ u=c(ET[i],ET[i+1]) ou=order(u) parcours[i]=paste(u[ou[1]],u[ou[2]],sep="|") trajet[i]=which(edg==parcours[i]) } parcours [1] "A|B" "B|D" "D|G" "E|G" "A|E" "A|C" "C|E" "E|H" "F|H" "F|I" "I|K" "H|K" "H|J" "G|J" "G|P" "J|P" [17] "J|L" "K|L" "K|Q" "L|Q" "H|L" "H|O" "F|O" "C|F" "C|N" "E|N" "B|E" "B|M" "A|M" trajet [1] 1 3 8 9 4 2 6 10 11 12 16 15 14 13 26 27 18 19 28 29 17 25 24 7 22 23 5 21 20 Let us try now on a real network of streets. Like Missoula, Montana. I will not try to get the shapefile of the city, I will just try to replicate the photography above. If you look carefully, you will see some problem : 10 and 93 have an odd number of vertices (3 here), so one strategy is to connect them (which explains the grey line). But actually, to be more realistic, we start in 93, and we end in 10. Here is the optimal (shortest) path which goes through all vertices. Now, we are ready for Halloween, to go through all streets in the neighborhood ! # « Dans toute statistique, l’inexactitude du nombre est compensée par la précision des décimales » Le statisticien et économiste Alfred Sauvy est resté dans les mémoires pour avoir inventé en 1952 le terme “tiers-monde”. Mais on lui a aussi attribué la paternité de la phrase suivante « dans toute statistique, l’inexactitude du nombre est compensée par la précision des décimales ». J’ai du l’entendre alors que j’étais étudiant, et depuis, elle me suit partout. J’y repensais l’autre jour, quand Mathieu Gallard mentionnait sur Twitter le graphique suivant (correspondant a la popularité du président de la république, en France, dans les 18 mois qui suivent l’élection). Le tweet disait (entre autres) “ est à ce stade de son quinquennat légèrement moins populaire que “. Pour rappel, ces courbes “de popularité” sont construites par un sondage, avec a chaque fois environ 1000 personnes interrogées (“sur la taille de l’échantillon on est toujours entre 950 et 970 interviews” me disait Mathieu). Bon, tous ceux qui ont des souvenirs de cours de stats se souviennent qu’avec 1000 personnes interrogées, 2 points de différence, c’est rarement significatif. Mais plus globalement, compte tenu de la marge d’erreur, je me suis demande pourquoi les courbes n’étaient pas lissées ? Ça éviterait les discussions stériles pour une variation de 2 points par exemple.. Si on reprend les données brutes (merci Mathieu), on a ici rate=read.csv2("http://freakonometrics.free.fr/satisfaction.csv") plot(rate[,2],type="b",ylim=c(.2,.75),xlim=c(0,20),pch=19) lines(rate[,3],type="b",col="red") lines(rate[,4],type="b",col="blue") lines(rate[,5],type="b",col="dark green") text(18.15,rate[17,2],"EM") text(18.15,rate[17,3],"FH",col="red") text(18.15,rate[17,4],"NS",col="blue") text(18.15,rate[17,5],"JC",col="dark green") ce qui donne le même que dans le tweet (même si je n’ose pas interpoler linéairement les valeurs manquantes – il y en a deux dans mon fichier) Prenons la courbe la plus récente, celle d’Emmanuel Macron (avec en plus les valeurs manquantes pour pimenter un peu) et rajoutons les intervalles de confiance ponctuels. plot(rate$EM,type="b",ylim=c(.2,.5)) p=rate$EM n=length(p) arrows(1:n,p-2/sqrt(1000)*sqrt(p*(1-p)),1:n,p+2/sqrt(1000)*sqrt(p*(1-p)),code=3,angle=90,length=.1,col="blue") Personnellement, j’aurais bien voulu (1) lisser tout ça, (2) rajouter quelque chose qui s’apparente a des bandes de confiance. Mais avec des erreurs de mesure (c’est comme ça qu’on peut interpréter le fait que les points viennent d’un sondage), je ne sais pas trop quoi faire. J’ai tenté la méthode suivante : le tire au hasard des points dans l’intervalle de confiance, puis je lisse sur ce nouveau jeu de points. Et je répète mille fois library(mgcv) Y=matrix(NA,73,1000) for(s in 1:1000){ x=(1:n)[!is.na(p)] pna=p[!is.na(p)] ps=rnorm(length(x),pna,1/sqrt(1000)*sqrt(pna*(1-pna))) b=data.frame(x=x,y=ps) reg=gam(y~s(x),data=b) yp=predict(reg,newdata=data.frame(x=seq(0,18,by=.25))) Y[,s]=yp if(s&lt;100) lines(seq(0,18,by=.25),yp,col="light blue") } lines(seq(0,18,by=.25),apply(Y,1,mean),col="red",lwd=2) lines(seq(0,18,by=.25),apply(Y,1,function(x) quantile(x,.95)),col="red",lty=2) lines(seq(0,18,by=.25),apply(Y,1,function(x) quantile(x,.05)),col="red",lty=2) On voit que notre courbe lissée est réaliste, voire même les pseudo-bandes de confiance autour. Pour obtenir ces trois courbes, on peut utiliser la fonction suivante courbe=function(j=1){ p=rate[,1+j] n=length(p) Y=matrix(NA,73,1000) for(s in 1:1000){ x=(1:n)[!is.na(p)] pna=p[!is.na(p)] ps=rnorm(length(x),pna,1/sqrt(1000)*sqrt(pna*(1-pna))) b=data.frame(x=x,y=ps) reg=gam(y~s(x),data=b) yp=predict(reg,newdata=data.frame(x=seq(0,18,by=.25))) Y[,s]=yp } data.frame( x=seq(0,18,by=.25), pred=apply(Y,1,mean), upr=apply(Y,1,function(x) quantile(x,.975)), lwr=apply(Y,1,function(x) quantile(x,.025))) } Sur les quatre colonnes de notre tableau, ça donne plot(rate[,2],type="b",ylim=c(0,.8),xlim=c(0,20),col="white") Y=courbe(4) polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col=rgb(0,1,0,.4),border=NA) lines(Y$x,Y$pred,col="dark green",lwd=2) text(18.65,Y$pred[73],"JC",col="dark green") Y=courbe(3) polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col=rgb(0,0,1,.4),border=NA) lines(Y$x,Y$pred,col="blue",lwd=2) text(18.65,Y$pred[73],"NS",col="blue") Y=courbe(2) polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col=rgb(1,0,0,.4),border=NA) lines(Y$x,Y$pred,col="red",lwd=2) text(18.65,Y$pred[73],"FH",col="red") Y=courbe(1) polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col="grey",border=NA) lines(Y$x,Y$pred,col="black",lwd=2) text(18.65,Y$pred[73],"EM",col="black") Pourquoi les instituts de sondages, qui produisent les courbes de popularité, ne montrent pas ce genre de courbes ? Elles sont – a mon avis – aussi justes que celles qu’ils fournissent, au centième près, jouant sur une précision que l’incertitude ne devrait pas autoriser… # (A brief) history of randomness, and simulation techniques Hearing “there is a 10% chance of rain today” or “the medical test has a positive predictive value of 75%” shows that the probabilities are now everywhere. A probability is a quantity that is difficult to grasp, but essential when trying to theorize and measure chance, or randomness. And if mathematical theory finally came very late, as Hacking (2006) points out, this did not prevent insurance from developing early enough, and from having the first (actuarial) mortality tables even before the “probability of death” or “life expectancy” had a mathematical basis. And in the same way, many techniques were invented to “generate randomness“, before the explosion of the so-called Monte Carlo methods, in parallel with the development of computing (and the fact that a machine could generate chance). Continue reading (A brief) history of randomness, and simulation techniques # On the robustness of LASSO Probably the last post on lasso, before the summer break… More specifically, I was wondering about the interpretation of graphs $\lambda\mapsto\widehat{\beta}_\lambda$. We use them for variable selection, but my major concern was about confidence intervals : how can we trust those lines ? As usual, a natural way is to use simulations on generated datasets. Consider for instance Sigma = matrix(c(1,.8,.2,.8,1,.4,.2,.4,1),3,3) n = 1000 library(mnormt) X = rmnorm(n,rep(0,3),Sigma) set.seed(123) df = data.frame(X1=X[,1],X2=X[,2],X3=X[,3],X4=rnorm(n), X5=runif(n), X6=exp(X[,3]), X7=sample(c("A","B"),size=n,replace=TRUE,prob=c(.5,.5)), X8=sample(c("C","D"),size=n,replace=TRUE,prob=c(.5,.5))) df$Y = 1+df$X1-df$X4+5*(df$X7=="A")+rnorm(n) One can use other simulations of datasets, and store the output vlambda = exp(seq(-8,1,length=201)) lasso = glmnet(x=X,y=df[,"Y"],family="gaussian",alpha=1, lambda=vlambda,standardize=TRUE) VLASSO[[s]] = as.matrix(lasso$beta)

To visualize confidence bands, one can compute quantiles

Q05=Q95=Qm=matrix(NA,9,201) for(i in 1:nrow(Q05)){ for(j in 1:ncol(Q05)){ v = unlist(lapply(VLASSO,function(x) x[i,j])) Q05[i,j] = quantile(v,.05) Q95[i,j] = quantile(v,.95) Qm[i,j] = mean(v) }}

and get get the graph

plot(lasso,col=colrs,"lambda"ylim=c(min(Q05),max(Q95))) colrs=c(brewer.pal(8,"Set1")) polygon(c(log(lasso$lambda),rev(log(lasso$lambda))), c(Q05[2,],rev(Q95[2,])),col=colrs[1],border=NA) polygon(c(log(lasso$lambda),rev(log(lasso$lambda))), c(Q05[5,],rev(Q95[5,])),col=colrs[2],border=NA) polygon(c(log(lasso$lambda),rev(log(lasso$lambda))), c(Q05[8,],rev(Q95[8,])),col=colrs[3],border=NA)

An alternative (more realistic on real data) is to use bootstrapped version of the dataset

id = sample(1:nrow(X),size=nrow(X),replace=TRUE) lasso = glmnet(x=X[id,],y=df[id,"Y"],family="gaussian",alpha=1, lambda=vlambda,standardize=TRUE)

So far, it looks it’s working very well. Now, what if we have a smaller dataset

n = 100

On simulated new samples, we get

while the bootstrap version is

There is more uncertainty, clearly, but the conclusion is not ambiguous here.

Now, what about real data. Consider the following

chicago = read.table("http://freakonometrics.free.fr/chicago.txt",header=TRUE,sep=";") tail(chicago) Fire X_1 X_2 X_3 42 4.8 0.152 19 13.323 43 10.4 0.408 25 12.960 44 15.6 0.578 28 11.260 45 7.0 0.114 3 10.080 46 7.1 0.492 23 11.428 47 4.9 0.466 27 13.731

with one variable of interest (the number of fires, per unhabitants) and 3 features. We can here use bootstrap to generate samples, and then fit a lasso regression. On the original dataset, the regression is

X = model.matrix(lm(Fire~.,data=chicago)) id = sample(1:nrow(X),size=nrow(X),replace=TRUE) vlambda = exp(seq(-4,2,length=201)) lasso = glmnet(x=X[id,],y=chicago[id,"Fire"],family="gaussian",alpha=1, lambda=vlambda,standardize=TRUE)

And if we just plot lines $\lambda\mapsto\widehat{\beta}_\lambda$ we get

Now, consider bootstrap samples.

for(s in 1:100){ id=sample(1:nrow(X),size=nrow(X),replace=TRUE) library(glmnet) vlambda=exp(seq(-4,2,length=201)) lasso=glmnet(x=X[id,],y=chicago[id,"Fire"],family="gaussian",alpha=1, lambda=vlambda,standardize=TRUE) plot(lasso,col=colrs,"lambda",lwd=.2,add=TRUE)}

We get here

The interpretation here is much more difficult

What about the order ?

N=matrix(NA,100000,4) for(s in 1:100000){ id=sample(1:nrow(X),size=nrow(X),replace=TRUE) library(glmnet) vlambda=exp(seq(-4,2,length=201)) lasso=glmnet(x=X[id,],y=chicago[id,"Fire"], family="gaussian",alpha=1, lambda=vlambda,standardize=TRUE) N[s,]=names(sort(apply(as.matrix(lasso$beta), 1,function(x) sum(x!=0))))} The ordering that was obtained on the original dataset was the same in 56% of the scenarios, mean(apply(N,1,function(x) paste(x,collapse="")=="(Intercept)X_1X_2X_3")) [1] 0.5693 We can look at all the cases, L=as.character(c(123,132,213,231,312,321)) Li=paste("(Intercept)X_",substr(L,1,1),"X_", substr(L,2,2),"X_",substr(L,3,3),sep="") g=function(y) mean(apply(N,1,function(x) paste(x,collapse="")==y)) vL=unlist(lapply(Li,g)) names(vL)=L barplot(vL,las=2,horiz=TRUE) # Standardization in LASSO The lasso regression is based on the idea of solving$$\widehat{\mathbf{\beta}}_{\lambda}=\text{argmin}\lbrace -\log\mathcal{L}(\mathbf{\beta}|\mathbf{x},\mathbf{y})+\lambda\|\mathbf{\beta}\|_{\ell_1}\rbrace$$where$$\Vert\mathbf{a} \Vert_{\ell_1}=\sum_{i=1}^d |a_i|$$for any $\mathbf{a}\in\mathbb{R}^d$. In a recent post, we’ve seen computational aspects of the optimization problem. But I went quickly throught the story of the $\ell_1$-norm. Because it means, somehow, that the value of $\beta_1$ and $\beta_2$ should be comparable. Somehow, with two significant variables, with very different scales, we should expect orders (or relative magnitudes) of $\widehat{\beta}_1$ and $\widehat{\beta}_2$ to be very very different. So people say that it is therefore necessary to center and reduce (or standardize) the variables. Consider the following (simulated) dataset Sigma = matrix(c(1,.8,.2,.8,1,.4,.2,.4,1),3,3) n = 1000 library(mnormt) X = rmnorm(n,rep(0,3),Sigma) set.seed(123) df = data.frame(X1=X[,1],X2=X[,2],X3=X[,3],X4=rnorm(n), X5=runif(n),X6=exp(X[,3]), X7=sample(c("A","B"),size=n,replace=TRUE,prob=c(.5,.5)), X8=sample(c("C","D"),size=n,replace=TRUE,prob=c(.5,.5))) df$Y = 1+df$X1-df$X4+5*(df$X7=="A")+rnorm(n) X = model.matrix(lm(Y~.,data=df)) Use the following colors for the graphs and the value of $\lambda$ library("RColorBrewer") colrs = c(brewer.pal(8,"Set1"))[c(1,4,5,2,6,3,7,8)] vlambda=exp(seq(-8,1,length=201)) The first regression we can run is a non-standardized one library(glmnet) lasso = glmnet(x=X,y=df[,"Y"],family="gaussian",alpha=1,lambda=vlambda,standardize=FALSE) We can visualize the graphs of $\lambda\mapsto\widehat{\beta}_\lambda$ idx = which(apply(lasso$beta,1,function(x) sum(x==0))&lt;200) plot(lasso,col=colrs,'lambda',xlim=c(-5.5,2.3),lwd=2) legend(1.2,.9,legend=paste('X',0:8,sep='')[idx],col=colrs,lty=1,lwd=2)

At least, observe that the most significant variables are the one that were used to generate the data.

Now, consider the case that we standardize the data

lasso = glmnet(x=X,y=df[,"Y"],family="gaussian",alpha=1,lambda=vlambda,standardize=TRUE)

The graphs of $\lambda\mapsto\widehat{\beta}_\lambda$

The graph is (strangely) very similar to the previous one. Except perhaps for the green curve. Maybe that categorical are not simular to continuous variables… Because somehow, standardisation of categorical variables might be not natural…

Why not consider some home-made function ? Let us transform (linearly) all variable in the $X$ matrix (except the first one, which is the intercept)

Xc = X for(j in 2:ncol(X)) Xc[,j]=(Xc[,j]-mean(Xc[,j]))/sd(Xc[,j])

Now, we can run our lasso regression on that one (with the intercept since all the variables are centered, but $y$)

lasso = glmnet(x=Xc,y=df$Y,family="gaussian",alpha=1,intercept=TRUE,lambda=vlambda) The plot is now plot(lasso,col=colrs,"lambda",xlim=c(-6.7,1.3),lwd=2) idx = which(apply(lasso$beta,1,function(x) sum(x==0))&lt;length(vlambda)) legend(.15,.45,legend=paste('X',0:8,sep='')[idx],col=colrs,lty=1,bty=&quot;n&quot;,lwd=2)

Actually, why not also center the $y$ variable, and remove also the intercept

Yc = (df[,"Y"]-mean(df[,"Y"]))/sd(df[,"Y"]) lasso = glmnet(x=Xc,y=Yc,family="gaussian",alpha=1,intercept=FALSE,lambda=vlambda)

Hopefully, those graphs are very consistent (and if we use those for variable selection, they suggest to use variables that were actually used to generate the dataset). And having qualitative and quantitative variable is not a big deal. But still, I do not feel confortable with the differences…

# Convex Regression Model

This morning during the lecture on nonlinear regression, I mentioned (very) briefly the case of convex regression. Since I forgot to mention the codes in R, I will publish them here. Assume that $y_i=m(\mathbf{x}_i)+\varepsilon_i$ where $m:\mathbb{R}^d\rightarrow \mathbb{R}$ is some convex function.

Then $m$ is convex if and only if $\forall\mathbf{x}_1,\mathbf{x}_2\in\mathbb{R}^d$, $\forall t\in[0,1]$, $$m(t\mathbf{x}_1+[1-t]\mathbf{x}_2) \leq tm(\mathbf{x}_1)+[1-t]m(\mathbf{x}_2)$$Hidreth (1954) proved that if$$m^\star=\underset{m \text{ convex}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-m(\mathbf{x_i})\big)^2\right\rbrace$$then $\mathbf{\theta}^\star=(m^\star(\mathbf{x_1}),\cdots,m^\star(\mathbf{x_n}))$ is unique.

Let $\mathbf{y}=\mathbf{\theta}+\mathbf{\varepsilon}$, then $$\mathbf{\theta}^\star=\underset{\mathbf{\theta}\in \mathcal{K}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-\theta_i)\big)^2\right\rbrace$$where$$\mathcal{K}=\{\mathbf{\theta}\in\mathbb{R}^n:\exists m\text{ convex },m(\mathbf{x}_i)=\theta_i\}$$. I.e. $\mathbf{\theta}^\star$ is the projection of $\mathbf{y}$ onto the (closed) convex cone $\mathcal{K}$. The projection theorem gives existence and unicity.

For convenience, in the application, we will consider the real-valued case, $m:\mathbb{R}\rightarrow \mathbb{R}$, i.e. $y_i=m(x_i)+\varepsilon_i$. Assume that observations are ordered $x_1\leq x_2\leq\cdots \leq x_n$. Here $$\mathcal{K}=\left\lbrace\mathbf{\theta}\in\mathbb{R}^n:\frac{\theta_2-\theta_1}{x_2-x_1}\leq \frac{\theta_3-\theta_2}{x_3-x_2}\leq \cdots \leq \frac{\theta_n-\theta_{n-1}}{x_n-x_{n-1}}\right\rbrace$$

Hence, quadratic program with $n-2$ linear constraints.

$m^\star$ is a piecewise linear function (interpolation of consecutive pairs $(x_i,\theta_i^\star)$).

If $m$ is differentiable, $m$ is convex if $$m(\mathbf{x})+ \nabla m(\mathbf{x})^{\text{T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})$$

More generally, if $m$ is convex, then there exists $\xi_{\mathbf{x}}\in\mathbb{R}^n$ such that $$m(\mathbf{x})+ \xi_{\mathbf{x}}^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})$$
$\xi_{\mathbf{x}}$ is a subgradient of $m$ at ${\mathbf{x}}$. And then $$\partial m(\mathbf{x})=\big\lbrace m(\mathbf{x})+ \xi^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y}),\forall \mathbf{y}\in\mathbb{R}^n\big\rbrace$$

Hence, $\mathbf{\theta}^\star$ is solution of $$\text{argmin}\big\lbrace\|\mathbf{y}-\mathbf{\theta}\|^2\big\rbrace$$$$\text{subject to }\theta_i+\xi_i^{\text{ T}}[\mathbf{x}_j-\mathbf{x}_i]\leq\mathbf{\theta}_j,~\forall i,j$$ and $\xi_1,\cdots,\xi_n\in\mathbb{R}^n$. Now, to do it for real, use cobs package for constrained (b)splines regression,

library(cobs)

To get a convex regression, use

plot(cars) x = cars$speed y = cars$dist rc = conreg(x,y,convex=TRUE) lines(rc, col = 2)

Here we can get the values of the knots

rc   Call: conreg(x = x, y = y, convex = TRUE) Convex regression: From 19 separated x-values, using 5 inner knots, 7, 8, 9, 20, 23. RSS = 1356; R^2 = 0.8766; needed (5,0) iterations

and actually, if we use them in a linear-spline regression, we get the same output here

reg = lm(dist~bs(speed,degree=1,knots=c(4,7,8,9,,20,23,25)),data=cars) u = seq(4,25,by=.1) v = predict(reg,newdata=data.frame(speed=u)) lines(u,v,col="green")

Let us add vertical lines for the knots

abline(v=c(4,7,8,9,20,23,25),col="grey",lty=2)

# Classification from scratch, logistic with splines 2/8

Today, second post of our series on classification from scratch, following the brief introduction on the logistic regression.

## Piecewise linear splines

To illustrate what’s going on, let us start with a “simple” regression (with only one explanatory variable). The underlying idea is natura non facit saltus, for “nature does not make jumps”, i.e. process governing equations for natural things are continuous. That seems to be a rather strong assumption, because we can assume that there is a fixed threshold to explain death. For instance, if patients die (for sure) if the “stroke index” exceeds a threshold, we might expect some discontinuity. Exceept that if that threshold is an heterogeneous (non-observable continuous) variable, then we get back to the continuity assumption.

The most simple model we can think of to extend the linear model we’ve seen in the previous post is to consider a piecewise linear function, with two parts : small values of $x$, and larger values of $x$. The most convenient way to do so is to use the positive part function $(x-s)_+$ which is the difference between $x$ and $s$ if that difference is positive, and $0$ otherwise. For instance $$\beta_1 x+\beta_2(x-s)_+$$ is the following piecewise linear function, continuous, with a “rupture” at knot $s$.

Observe also the following interpretation: for small values of $x$, there is a linear increase, with slope $\beta_1$, and for lager values of $x$, there is a linear decrease, with slope $\beta_1+\beta_2$. Hence, $\beta_2$ is interpreted as a change of the slope.

And of course, it is possible to consider more than one knot. The function to get the positive value is the following

pos = function(x,s) (x-s)*(x&gt;=s)

then we can use it direcly in our regression model

reg = glm(PRONO~INSYS+pos(INSYS,15)+ pos(INSYS,25),data=myocarde,family=binomial)

The output of the regression is here

summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.1109 3.2783 -0.034 0.9730 INSYS -0.1751 0.2526 -0.693 0.4883 pos(INSYS, 15) 0.7900 0.3745 2.109 0.0349 * pos(INSYS, 25) -0.5797 0.2903 -1.997 0.0458 *

Hence, the original slope, for very small values is not significant, but then, above 15, it become significantly positive. And above 25, there is a significant change again. We can plot it to see what’s going on

u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,type="l") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Using bs() linear splines

Using the GAM function, things are slightly different. We will use here so called b-splines,

library(splines)

We can define spline functions with support $(5,55)$ and with knots $\{15,25\}$

clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02") x = seq(0,60,by=.25) B = bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=1) matplot(x,B,type="l",lty=1,lwd=2,col=clr6)

as we can see, the functions defined here are different from the one before, but we still have (piecewise) linear functions on each segment $(5,15)$, $(15,25)$ and $(25,55)$. But linear combinations of those functions (the two sets of functions) will generate the same space. Said differently, if the interpretation of the output will be different, predictions should be the same

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=1), data=myocarde,family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.9863 2.0555 -0.480 0.6314 bs(INSYS,..)1 -1.7507 2.5262 -0.693 0.4883 bs(INSYS,..)2 4.3989 2.0619 2.133 0.0329 * bs(INSYS,..)3 5.4572 5.4146 1.008 0.3135

Observe that there are three coefficients, as before, but again, the interpretation is here more complicated…

v=predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Nevertheless, the prediction is the same… and that’s nice.

## Piecewise quadratic splines

Let us go one step further… Can we have also the continuity of the derivative ? Yes, and that’s easy actually, considering parabolic functions. Instead of using a decomposition on $x,(x-s_1)_+$ and $(x-s_2)_+$ consider now a decomposition on $x,x^{\color{red}{2}},(x-s_1)^{\color{red}{2}}_+$ and $(x-s_2)^{\color{red}{2}}_+$.

 pos2 = function(x,s) (x-s)^2*(x&gt;=s) reg = glm(PRONO~poly(INSYS,2)+pos2(INSYS,15)+pos2(INSYS,25), data=myocarde,family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 29.9842 15.2368 1.968 0.0491 * poly(INSYS, 2)1 408.7851 202.4194 2.019 0.0434 * poly(INSYS, 2)2 199.1628 101.5892 1.960 0.0499 * pos2(INSYS, 15) -0.2281 0.1264 -1.805 0.0712 . pos2(INSYS, 25) 0.0439 0.0805 0.545 0.5855

As expected, there are here five coefficients: the intercept and two for the part on the left (three parameters for the parabolic function), and then two additional terms for the part in the center – here $(15,25)$ – and for the part on the right. Of course, for each portion, there is only one degree of freedom since we have a parabolic function (three coefficients) but two constraints (continuity, and continuity of the first order derivative).

On a graph, we get the following

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2,xlab="INSYS",ylab="") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Using bs() quadratic splines

Of course, we can do the same with our R function. But as before, the basis of function is expressed here differently

 x = seq(0,60,by=.25) B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=2) matplot(x,B,type="l",xlab="INSYS",col=clr6)

If we run R code, we get

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=2),data=myocarde, family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 7.186 5.261 1.366 0.1720 bs(INSYS, ..)1 -14.656 7.923 -1.850 0.0643 . bs(INSYS, ..)2 -5.692 4.638 -1.227 0.2198 bs(INSYS, ..)3 -2.454 8.780 -0.279 0.7799 bs(INSYS, ..)4 6.429 41.675 0.154 0.8774

But that’s not really a big deal since the prediction is exactly the same

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Cubic splines

Last, but not least, we can reach the cubic splines. With our previous notions, we would consider a decomposition on (guess what) $x,x^2,x^{\color{red}{3}},(x-s_1)^{\color{red}{3}}_+,(x-s_2)^{\color{red}{3}}_+$, to get this time continuity, as well as continuity of the first two derivatives (and to get a very smooth function, since even variations will be smooth). If we use the bs function, the basis is the followin

B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=3) matplot(x,B,type="l",lwd=2,col=clr6,lty=1,ylim=c(-.2,1.2)) abline(v=c(5,15,25,55),lty=2)

and the prediction will now be

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=3), data=myocarde,family=binomial) u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Two last things before concluding (for today), the location of the knots, and the extension to additive models.

## Location of knots

In many applications, we do not want to specify the location of the knots. We just want – say – three (intermediary) knots. This can be done using

reg = glm(PRONO~1+bs(INSYS,degree=1,df=4),data=myocarde,family=binomial)

We can actually get the locations of the knots by looking at

attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 1L, knots = c(15.8, 21.4, 27.15), Boundary.knots = c(8.7, 54), intercept = FALSE) which provides us with the location of the boundary knots (the minumun and the maximum from from our sample) but also the three intermediary knots. Observe that actually, those five values are just (empirical) quantiles quantile(myocarde$INSYS,(0:4)/4) 0% 25% 50% 75% 100% 8.70 15.80 21.40 27.15 54.00

If we plot the prediction, we get

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2) If we get back on what was computed before the logit transformation, we clealy see ruptures are the different quantiles B = bs(x,degree=1,df=4) B = cbind(1,B) y = B%*%coefficients(reg) plot(x,y,type="l",col="red",lwd=2) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2)

Note that if we do specify anything about knots (number or location), we get no knots…

reg = glm(PRONO~1+bs(INSYS,degree=2),data=myocarde,family=binomial) attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 2L, knots = numeric(0), Boundary.knots = c(8.7,54), intercept = FALSE) and if we look at the prediction u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) actually, it is the same as a quadratic regression (as expected actually) reg = glm(PRONO~1+poly(INSYS,degree=2),data=myocarde,family=binomial) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) ## Additive models Consider now the second dataset, with two variables. Consider here a model like $$\mathbb{P}[Y|X_1=x_1,X_2=x_2]=\frac{\exp[\eta(x_1,x_2)]}{1+\exp[\eta(x_1,x_2)]}$$ where $$\exp[\eta(x_1,x_2)]=\beta_0+\color{red}{s_1(x_1)}+\color{blue}{s_2(x_2)}$$ $$\color{red}{s_1(x_1)}=\beta_{1,0}x_1+\beta_{1,1}(x_1-s_{11})_++\beta_{1,2}(x_1-s_{12})_+$$ and $$\color{blue}{s_2(x_2)}=\beta_{2,0}x_2+\beta_{2,1}(x_2-s_{21})_++\beta_{2,2}(x_2-s_{22})_+$$ It might seem a little bit restrictive, but that’s actually the idea of additive models. reg = glm(y~bs(x1,degree=1,df=3)+bs(x2,degree=1,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Now, if think about is, we’ve been able to get a “perfect” model, so, somehow, it seems no longer continuous…

persp(u,u,v,theta=20,phi=40,col="green"

Of course, it is… it is piecewise linear, with hyperplane, some being almost vertical.

And one can also consider piecewise quadratic functions

reg = glm(y~bs(x1,degree=2,df=3)+bs(x2,degree=2,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) Funny thing, we now have two “perfect” models, with different areas for the white and the black dots… Don’t ask me how to choose on that one. In R, it is possible to use the mgcv package to run a gam regression. It is used for generalized additive models, but here, we have only one variable, so it is difficult to see the “additive” part, actually. And to be more specific, mgcv is using penalized quasi-likelihood from the nlme package (but we’ll get back on penalized routines later on). But maybe I should also mention another smoothing tool before, kernels (and maybe also $k$-nearest neighbors). To be continued # Classification from scratch, logistic regression 1/8 Let us start today our series on classification from scratch The logistic regression is based on the assumption that given covariates $\mathbf{x}$, $Y$ has a Bernoulli distribution,$$Y|\mathbf{X}=\mathbf{x}\sim\mathcal{B}(p_{\mathbf{x}}),~~~~p_\mathbf{x}=\frac{\exp[\mathbf{x}^T\mathbf{\beta}]}{1+\exp[\mathbf{x}^T\mathbf{\beta}]}$$The goal is to estimate parameter $\mathbf{\beta}$. Recall that the heuristics for the use of that function for the probability is that$$\log[\text{odds}(Y=1)]=\log\frac{\mathbb{P}[Y=1]}{\mathbb{P}[Y=0]}=\mathbf{x}^T\mathbf{\beta}$$ ## Maximimum of the (log)-likelihood function The log-likelihood is here$$\log\mathcal{L} = \sum_{i=1}^n y_i\log p_i+(1-y_i)\log (1-p_i)$$ where $p_{i}=(1+\exp[-\mathbf{x}_i^T\mathbf{\beta}])^{-1}$. Numerical techniques are based on (numerical) gradient descent to compute the maximum of the likelihood function. The (negative) log-likelihood is the following function y = myocarde$PRONO X = cbind(1,as.matrix(myocarde[,1:7])) negLogLik = function(beta){ -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta))) }

We use the minus sign since standard optimization routines compute minima, not maxima. Now, to find the minimum of that function, we need a starting point to initiate the algorithm

beta_init = lm(PRONO~.,data=myocarde)$coefficients Why not start with the parameter of the OLS. Somehow, we might think that at least, sign should be ok for instance. Anyway, we need a starting point, and let us use that one. logistic_opt = optim(par = beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) Here, we obtain  logistic_opt$par (Intercept) FRCAR INCAR INSYS 1.656926397 0.045234029 -2.119441743 0.204023835 PRDIA PAPUL PVENT REPUL -0.102420095 0.165823647 -0.081047525 -0.005992238

Let us verify here that this output is valid. For instance, what if we change the value of the starting point (randomly)

simu = function(i){ logistic_opt_i = optim(par = rnorm(8,0,3)*beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) logistic_opt_i$par[2:3] } v_beta = t(Vectorize(simu)(1:1000)) plot(v_beta) par(mfrow=c(1,2)) hist(v_beta[,1],xlab=names(myocarde)[1]) hist(v_beta[,2],xlab=names(myocarde)[2]) Ooops. There is a problem here. Clearly, we cannot rely on numerical optimization here. We can think about using another optimization routine library(optimx) logit = function(mX, vBeta) { exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) } logLikelihoodLogitStable = function(vBeta, mX, vY) { -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta))) + (1-vY)*(-log(1 + exp(mX %*% vBeta)))) } likelihoodScore = function(vBeta, mX, vY) { return(t(mX) %*% (logit(mX, vBeta) - vY) ) } optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable, method = 'L-BFGS-B', gr = likelihoodScore, mX = X, vY = y, hessian=TRUE) The optimum is here attr(optimLogitLBFGS, "details")[[2]] [,1] 0.066680272 FRCAR 0.003080542 INCAR 0.079031364 INSYS -0.001586194 PRDIA 0.040500697 PAPUL -0.041870705 PVENT -0.014162756 REPUL 0.195632244 Let’s be honest here, I do not feel confortable with those techniques. So, what happened here ? Here, the technique we use is based on the following idea,$$\mathbf{\beta}_{new}=\mathbf{\beta}_{old} -\left(\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}\right)^{-1}\cdot \frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}$$The problem is that my computer does not know this first and second derivatives. So it will compute them using approximation techniques. Actually, it is possible to use functions dedicated to such computation library(numDeriv) library(MASS) logit = function(x){1/(1+exp(-x))} logLik = function(beta, X, y){ -sum(y*log(logit(X%*%beta)) + (1-y)*log(1-logit(X%*%beta))) } optim_second = function(beta, num_iter){ LL = vector() for(i in 1:num_iter){ grad = (t(X)%*%(logit(X%*%beta) - y)) H = hessian(logLik, beta, method = "complex", X = X, y = y) beta = beta - ginv(H)%*%grad LL[i] = logLik(beta, X, y) } result = list(beta, H) return(result) } With our OLS starting point, we obtain opt0 = optim_second(beta_init,500) opt0[[1]] [,1] [1,] 0.951074420 [2,] 0.018860280 [3,] 0.275428978 [4,] 0.144803636 [5,] -0.058535606 [6,] 0.001182178 [7,] -0.108651776 [8,] -0.002940315 But if we try with another starting point opt1 = optim_second(beta_init*runif(8),500) opt1[[1]] [,1] [1,] 0.052894794 [2,] 0.024718435 [3,] 0.167953661 [4,] 0.171662947 [5,] -0.057458066 [6,] -0.011361034 [7,] -0.107532114 [8,] -0.002679064 Clearly, some coefficients are rather close. But other aren’t. From my point of viezw, that is a major problem (keep in mind that we do not deal here with massive data ! There are only 7 explanatory variables, and only 71 observations). Why not try to be clever, and use the analytical values of those derivatives ? Even if some people claim the oppositive, sometimes, it can actually be usefull to do the maths, instead of considering only numerical values. ## Newton (or Fisher) Algorithm If you open any Econometrics textbooks (one can also try to derive it), you will get $$\frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}=\mathbf{X}^T(\mathbf{y}-\mathbf{p}_{old})$$ while$$\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=-\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}$$ Y=myocarde$PRONO X=cbind(1,as.matrix(myocarde[,1:7])) colnames(X)=c("Inter",names(myocarde[,1:7])) beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} Observe that here, I use only ten iterations of the algorithm !  beta[,8:10] [,1] [,2] [,3] XInter -10.187641685 -10.187641696 -10.187641696 XFRCAR 0.138178119 0.138178119 0.138178119 XINCAR -5.862429035 -5.862429037 -5.862429037 XINSYS 0.717084018 0.717084018 0.717084018 XPRDIA -0.073668171 -0.073668171 -0.073668171 XPAPUL 0.016756506 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187 The thing is that is seems to converge extremely fast. And it is rather robust ! Look at what we get if we change our starting point beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} beta[,8:10] [,1] [,2] [,3] XInter -10.187641586 -10.187641696 -10.187641696 XFRCAR 0.138178118 0.138178119 0.138178119 XINCAR -5.862429017 -5.862429037 -5.862429037 XINSYS 0.717084013 0.717084018 0.717084018 XPRDIA -0.073668172 -0.073668171 -0.073668171 XPAPUL 0.016756508 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187

Nice, isn’t it? Looks like we got our winner, don’t we? And one can use the inverse of the Hessian matrix to get standard deviations.

## Weighted Least-Squares

Let us go one step further. We’ve seen that we want to compute something like$$\mathbf{\beta}_{new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}$$(if we do substitute matrices in the analytical expressions) where $\mathbf{z}=\mathbf{X}\mathbf{\beta}_{old}+\mathbf{\Delta}_{old}^{-1}[\mathbf{y}-\mathbf{p}_{old}]$. But actually, that’s simply a standard least-square problem$$\mathbf{\beta}_{new} = \text{argmin}\left\lbrace(\mathbf{z}-\mathbf{X}\mathbf{\beta})^T\mathbf{\Delta}_{old}^{-1}(\mathbf{z}-\mathbf{X}\mathbf{\beta})\right\rbrace$$The only problem here is that weights $\mathbf{\Delta}_{old}$ are functions of unknown $\mathbf{\beta}_{old}$. But actually, if we keep iterating, we should be able to solve it : given the $\mathbf{\beta}$ we got the weights, and with the weights, we can use weighted OLS to get an updated $\mathbf{\beta}$. That’s the idea of iteratively reweighted least squares.

The algorithm will be

df = myocarde beta_init = lm(PRONO~.,data=df)$coefficients X = cbind(1,as.matrix(myocarde[,1:7])) beta = beta_init for(s in 1:1000){ p = exp(X %*% beta) / (1+exp(X %*% beta)) omega = diag(nrow(df)) diag(omega) = (p*(1-p)) df$Z = X %*% beta + solve(omega) %*% (df$PRONO - p) beta = lm(Z~.,data=df[,-8], weights=diag(omega))$coefficients }

and the output is here

 beta (Intercept) FRCAR INCAR INSYS PRDIA -10.187641696 0.138178119 -5.862429037 0.717084018 -0.073668171 PAPUL PVENT REPUL 0.016756506 -0.106776012 -0.003154187

which is almost what we’ve obtained before. Nice isn’t it ? Actually, here we also have standard deviations of estimators

summary( lm(Z~.,data=df[,-8], weights=diag(omega)))   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) -10.187642 10.668138 -0.955 0.343 FRCAR 0.138178 0.102340 1.350 0.182 INCAR -5.862429 6.052560 -0.969 0.336 INSYS 0.717084 0.503527 1.424 0.159 PRDIA -0.073668 0.261549 -0.282 0.779 PAPUL 0.016757 0.306666 0.055 0.957 PVENT -0.106776 0.099145 -1.077 0.286 REPUL -0.003154 0.004386 -0.719 0.475

## The standard glm function

Of course, it is possible to use an R built-in function to get our estimate

summary(glm(PRONO~.,data=myocarde,family=binomial(link = "logit")))   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -10.187642 11.895227 -0.856 0.392 FRCAR 0.138178 0.114112 1.211 0.226 INCAR -5.862429 6.748785 -0.869 0.385 INSYS 0.717084 0.561445 1.277 0.202 PRDIA -0.073668 0.291636 -0.253 0.801 PAPUL 0.016757 0.341942 0.049 0.961 PVENT -0.106776 0.110550 -0.966 0.334 REPUL -0.003154 0.004891 -0.645 0.519

## Application and visualisation

Let us visualize the prediction obtained from the logistic regression, on our second dataset

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x,x2=y,y=as.factor(z)) reg = glm(y~x1+x2,data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Here level curves – or iso-probabilities – are linear, so the space is divided in two (0 and 1, survival and death, white and black) by a straight line (or an hyperplane in higher dimension). Furthermore, since we have a linear model, if we change the cutoff (the threshold used to create the two classes), we obtain another straight line (or hyperplane) parallel to the first one.

Next time, we will introduce splines to smooth those continuous covariates… to be continued.

# Classification from scratch, overview 0/8

Before my course on « big data and economics » at the university of Barcelona in July, I wanted to upload a series of posts on classification techniques, to get an insight on machine learning tools.

According to some common idea, machine learning algorithms are black boxes. I wanted to get back on that saying. First of all, isn’t it the case also for regression models, like generalized additive models (with splines) ? Do you really know what the algorithm is doing ? Even the logistic regression. In textbooks, we can easily find math formulas. But what is really done when I run it, in R ?

When I started working on academia, someone told me something like « if you really want to understand a theory, teach it ». And that has been my moto for more than 15 years. I wanted to add a second part to that statement: « if you really want to understand an algorithm, recode it ». So let’s try this… My ambition is to recode (more or less) most of the standard algorithms used in predictive modeling, from scratch, in R. What I plan to mention, within the next two weeks, will be

I will use two datasets to illustrate. The first one is inspired by the cover of « Foundations of Machine Learning » by Mehryar Mohri, Afshin Rostamizadeh and Ameet Talwalkar. At least, with this dataset, it will be possible to plot predictions (since there are only two – continuous – features)

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x,x2=y,y=as.factor(z)) plot(x,y,pch=c(1,19)[1+z])

Here is some code to get a visualization of the prediction (here the probability to be a black point)

rmatrix_model = function(model){ u = seq(0,1,length=101) p = function(x,y) predict(model,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) return(v)} nice_graph=function(v){ u = seq(0,1,length=101) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10[c(1,10)],breaks=c(0,5,10)/10) points(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) } reg = glm(y~x1+x2,data=df,family=binomial) nice_graph(rmatrix_model(reg))

Note that colors are defined here as

clr10= c("#ffffff","#f7fcfd","#e5f5f9","#ccece6","#99d8c9","#66c2a4","#41ae76","#238b45","#006d2c","#00441b")

or with some nonlinear model

The second one is a dataset I got from Gilbert Saporta, about heart attacks and decease (our binary variable).

myocarde = read.table("http://freakonometrics.free.fr/myocarde.csv",head=TRUE, sep=";") myocarde$PRONO = (myocarde$PRONO=="SURVIE")*1 y = myocarde$PRONO X = as.matrix(cbind(1,myocarde[,1:7])) So far, I do not plan to talk (too much) on the choice of tunning parameters (and cross-validation), on comparing models, etc. The goal here is simply to understand what’s going on when we call either glm, glmnet, gam, random forest, svm, xgboost, or any function to get a predict model. # On the interpretation of a regression model Yesterday, (aka @NaytaData ) posted a nice graph on reddit, with bicycle traffic and mean air temperature, in Helsinki, Finland, per day, I found that graph interesting, so I did ask for the data ( kindly sent them to me tonight). df=read.csv("cyclistsTempHKI.csv") library(ggplot2) ggplot(df, aes(meanTemp, cyclists)) + geom_point() + geom_smooth(span = 0.3) But as mentioned by someone on twitter, the interpretation is somehow trivial : people get out on their bike when the weather is nice. The hotter, the more cyclists on the road. Which is interpreted here in a causal way… But actually, we can also visualize the data as follows, as suggested by Antoine Chambert-Loir  ggplot(df, aes(cyclists, meanTemp)) + geom_point() + geom_smooth(span = 0.3) The interpretation would be, somehow, that the more cyclists on the road, the hotter it is. Why not consider this causal interpretation here ? Like cyclists go so fast, or sweat so much, that they increase temperature… Of course, it is the standard (recurrent) discussion “correlation is not causality”, but in regression models, we like to tell a story, to pretend that we have some sort of a causal story. But we do not prove it. Here, we know that the first one is more credible than the second one, but how do we know that ? To go further, how can we use machine learning techniques to prove causal relationships ? How could a machine choose between the first and the second story ? # Les modèles prédictifs peuvent-il être justes ? Dans Nosedive (traduit par le titre Chute Libre en France), le premier épisode de la saison 3 de la série télévisée Black Mirror, on découvre la dystopie d’une une société régie par une « cote personnelle », une note, un score allant de 0 à 5. Dans ce monde, chaque personne note les autres, les mieux notés ayant accès à de meilleurs services (priorité dans les services, meilleurs taux, meilleurs prix, etc). Cette tendance à construire des scores dans toutes sortes de domaines (historiquement sur les crédits mais aujourd’hui sur des aspects criminels, voire civiques dans certains pays) ne va-t-elle pas déboucher sur un monde qui serait un concours de popularité sans fin ? Et comment serait-elle conciliable avec une justice sociale, a priori souhaitable ? # Les scores de crédit et les réseaux sociaux Un score de crédit est, d’un point de vue actuariel, une grandeur proportionnelle à la probabilité de ne pas honorer ses engagements en tant que créditeur. Ça peut être de plus pouvoir payer les échéances trois mois consécutifs, ou juste d’avoir un retard. Dans la vraie vie, comme toujours, c’est un peu plus compliqué. Aux Etats-Unis ou en Grande Bretagne, il n’est pas rare que les étudiants s’endettent sur des dizaines d’années pour avoir l’opportunité de suivre les cours qui les intéressent (même si la motivation est surtout d’obtenir un diplôme en fin de parcours). Mais surtout, dès qu’ils atteignent d’âge de 18 ans, des sociétés de notation de crédit vont surveiller tous leurs déplacements. Souvent à leur insu. Et si un jour, un crédit consommation ou hypothécaire est refusé, les raisons ne sont jamais motivées. Est-ce un retard dans un paiement de loyer ? Des amendes de bibliothèques oubliées ? Une facture d’eau impayée, vieille de plusieurs années ? Les sociétés de notation de crédit aux États-Unis, mais aussi en Chine, commencent à étudier la possibilité d’utiliser des données provenant de médias sociaux pour améliorer le score de crédit. Compter le nombre de fois où un utilisateur utilise le mot « gaspillé » (« wasted » en anglais) dans ce qu’il poste en ligne ne peut-il pas révéler une information quant au remboursement de dettes ? C’est en tout cas ce que prétend l’analyste de crédit américain FICO: « If you look at how many times a person says ‘wasted’ in their profile, it has some value in predicting whether they’re going to repay their debt (…) It’s not much, but it’s more than zero » (cité dans McLannahan (2015)). En Chine, le prêteur peer-to-peer Jubao a révélé qu’il était plus susceptible de donner des « bonus » aux emprunteurs s’ils étaient des amis Facebook avec des célébrités, tel que le raconte Botsman (2017). Pour l’instant, les sociétés de notation de crédit utilisent encore les données qu’elles connaissent bien (factures de services publics et cartes de crédit), mais elles imaginent que bien des informations intéressantes doivent être accessibles (d’une manière ou une autre) sur les réseaux sociaux. Mais les données sont encore rares, et difficiles à analyser. Quid de la composante sarcastique ou humoristique dans un tweet utilisant le mot « wasted » ? Comme souvent, la difficulté est que les données réellement pertinentes sont difficiles à obtenir. S’il est possible d’avoir des informations sur le payement du loyer quand un locataire passe par une agence, que faire pour les transactions entre deux particuliers ? Et si c’était possible, comment traiter le cas de colocataires ? Ne pas obtenir de crédit parce qu’un ancien colocataire n’a pas payé dans les temps devient dérangeant. D’autant plus s’il s’agit peut-être d’une facture de téléphone cellulaire réclamée abusivement par la compagnie de téléphonie, alors que l’abonnement avait été résilié. Mais le gros « malus » dans le score de crédit est bien souvent le fait de ne jamais avoir eu de carte de crédit. On pourrait penser qu’une personne qui n’a pas eu besoin d’une carte de crédit (et se contentait d’une carte de débit, permettant d’acheter chez un commerçant, comme la majorité des cartes bancaires en France) est le propre d’une personne prudente, qui n’a pas besoin de crédit pour des dépenses quotidiennes. Mais pour les établissements de crédit, cette personne n’est pas fiable car on ne la connait pas. Et c’est à elle de prouver qu’elle l’est (on revient à la pratique récurrente d’inversion de la charge de la preuve évoqué dans Charpentier(2016)). C’est étrangement ce qui se passe aujourd’hui quand on veut entrer sur le sol américain sans avoir de page Facebook. # Dans un monde de surveillance généralisée Et si les établissements de crédit n’étaient pas les seuls intéressés par notre vie ? Que serait un monde si, en plus de savoir si je paye mes factures à temps, certains cherchaient à connaître mes réseaux d’amis, à savoir quels journaux je lis, si je préfère acheter du lait entier ou du lait demi-écrémé ? Quand on visite le musée de la Stasi à Berlin, on découvre ce monde a existé, que 1 personne sur 63 était agent (ou indicateur) de la Stasi (en comptant les indicateurs occasionnels, la proportion peut atteindre une personne sur 6). Le musée décrit un panoptisme total, chacun étant observé en permanence, comme le décrit Foucault (1975). Mais ce cauchemar ne correspond-il pas à monde actuel, de surveillance permanente, plus ou moins consentie. Surveillance via les téléphones cellulaires (géolocalisation pour la fonction la plus courante, mais parfois aussi des enregistrements audio à l’insu de l’utilisateur par certaines applications), via les objets connectés, mais aussi les caméras de surveillance couplées à des algorithmes de reconnaissance faciale de plus en plus performant. Fin 2017, en Chine 170 millions de caméra étaient installées, et le cap des 300 millions devrait être atteint d’ici 2020. Lors d’une expérience tentée par la BBC[1], il a fallu 7 minutes pour retrouver le journaliste John Sudworth qui se promenait dans les rues. Le danger est que l’on ne sait jamais trop qui contrôle. De plus en plus de compagnies privées de sécurité se sont associées aux gouvernements. Les fournisseurs de messageries électroniques lisent nos messages pour détecter du spam, mais aussi pour revendre certaines informations. Par exemple, dans les règles de confidentialité annexées aux conditions générales d’utilisation de Gmail (Google) on lit « Nos systèmes automatisés analysent vos contenus (y compris les e-mails) afin de vous proposer des fonctionnalités personnalisées sur les produits, telles que (…) des publicités sur mesure ». Les assureurs envisagent de plus de plus l’installation de boitiers GPS dans les voitures, mais en passant par des prestataires externes. Au-delà de la propriété des données (évoqué dans Charpentier & Suire (2016)) on peut s’interroger sur leur revente, et leur utilisation. Savoir que quelqu’un se rend régulièrement dans un centre de transfusion sanguine est une information potentiellement intéressante, surtout couplée à d’autres. Depuis 2014, le gouvernement chinois travaille sur un système d’évaluation de ses propres citoyens, programmé pour être mis en place en 2020, comme le raconte Trujillo (2017). Ce « système de crédit social » vise à créer un « score citoyen » (pour reprendre l’expression de Galeon & Bergan (2017)), afin de prédire, et prévenir, les dangers potentiels, normalisant les conduites individuelles par des dispositifs panoptiques (par exemple la vidésurveillance), induisant des réflexes d’autodéfense et d’autocontrôle. Comme le disait Foucault (1975), il s’agit de « faire que la surveillance soit permanente dans ses effets, même si elle est discontinue dans son action; que la perfection du pouvoir tende à rendre inutile l’actualité de son exercice » (même souvent, aujourd’hui, elle est en plus continue dans son action). Certains de ces scores sont utilisés par la police pour savoir où patrouiller pour faire baisser la criminalité, comme PredPol. Mais quand on y regarde de plus près, les prédictions disent, en substance, que les crimes auront lieu (en majorité) dans les zones (historiquement) les plus criminogènes de la ville. La frontière entre la banalité et la tautologie est étroite. Et le réel danger est que bien souvent, les scores transforment les probabilités en quasi-certitudes, et le soupçon devient une preuve, comme le notait Supiot (2015). # Justice prédictive et méthodes actuarielles En juin 2010, un rapport de l’Académie de médecine préconisait d”« améliorer la pratique des expertises de dangerosité des criminels sexuels en enseignant et en diffusant les méthodes actuarielles ». Ces « méthodes actuarielles » sont tout simplement les techniques de scoring, de « profilage » tel que le défini le règlement européen relatif aux données personnelles du 27 avril 2016 (RGDP). Angèle Christin s’est intéressée aux algorithmes qui estiment la probabilité de récidive dans la justice pénale américaine. Comme elle l’a montré, ces techniques posent nombre de questions, en particulier des biais discriminatoires, l’opacité qui rend difficile les recours, et surtout la difficulté de comprendre ce qui est réellement calculé. Dans l’État de Virginie, un score entre 1 et 10 est renvoyé, convention reprise par Compas (Correctional Offender Management Profiling Alternative Sanctions ) qui offre en plus un code couleur qui prédit le risque de récidive violent. Il s’agit alors d’un outil d’aide à la décision, une machine ne pouvant placer une personne en détention seule (Christin et al. (2015)). Les conclusions d’un score prédictif dépendent de deux éléments clés : le modèle utilisé, et les données. Dans la majorité des cas aux Etats Unis, les codes des modèles restent opaques (et donc impossible à attaquer), et rares sont ceux qui ont vu les données utilisées pour calibrer ces modèles. Mais on peut se demander si les décisions de justice ne sont pas eux aussi relativement opaques ? Certes, les juges doivent motiver leurs décisions, ce qui les rend critiquables et attaquables, mais si le processus était si transparent, les issues d’un procès (humains) ne devraient-elles pas être alors davantage prévisibles ? Enfin, les différents biais sont assez simples à comprendre. Supposons qu’être riche permet d’avoir un bon avocat, et avoir un bon avocat permet de ne pas avoir certaines condamnations. Dans ce cas, une variable liée à la richesse (le type de véhicule possédé par exemple) sera liée positivement avec le fait de ne pas être coupable (reconnu coupagle), et fera baisser le score de dangerosité. L’autre danger dans les biais de sélection est qu’ils sont parfois complexes à comprendre, voire paradoxaux. Un exemple classique est celui illustré sur la Figure 1. Pendant la seconde guerre mondiale, il a été demandé à des ingénieurs et des statisticiens comment renforcer les bombardiers qui faisaient face au feu ennemi. Figure1 : Endroits endommagés des avions revenus (source : McGeddon 2016) Le statisticien Abraham Wald a commencé à collecter des données sur les impacts dans la carlingue, comme le raconte Mangel & Samaniego (1984). A la surprise générale, il a recommandé de blinder les endroits des appareils qui présentaient le moins de dommages. En effet, les avions utilisés dans l’échantillon présentait un biais important : seuls les avions revenus ont été pris en compte. S’ils ont pu revenir avec des trous au bout des ailes, c’est que ces parties sont suffisamment solides. Et comme aucun avion n’est revenu avec des trous au niveau des moteurs des hélices, ce sont ces parties qu’il convenait de renforcer. Un autre danger est celui où les relations causales sont inversées. Que penser de ce médecin qui prescrit un puissant neuroleptique à un patient mis en examen, de peur que la justice lui reproche de ne pas avoir vu la dangerosité de son patient, et qu’inversement, la justice s’appuie sur cette prescription pour prouver que le patient est dangereux ? Un algorithme mal conçu pourrait comprendre de travers le sens des relations causales. Mais les modèles prédictifs en matière judiciaire ne sont pas que du côté des juges. Lors d’un accident corporel sur la route, la loi Badinter (du 5 juillet 1985) prévoit un « droit à indemnisation » pour toute victime d’un accident de la circulation dans lequel est impliqué un véhicule terrestre à moteur. Lorsque la compagnie d’assurance du conducteur propose une indemnité, la victime fait une rapide analyse coût/bénéfice pour savoir s’il va au tribunal. Si elle ne construit pas formellement un modèle prédictif, elle tente de voir, à partir de quelques éléments à sa connaissance, les coûts de demander à un juge de trancher sur le montant de l’indemnité, et ses bénéfices (potentiels). Autre point important, les juristes appellent ces modèles « prédictifs » des modèles « actuariels ». Or la première fonction des actuaires était d’actualiser, de calculer la valeur du temps. Et le temps judiciaire a des conséquences souvent désastreuses. En quoi une décision humaine, imparfaite, prise au bout de 5 ans de procédure serait « meilleure » qu’une décision automatique prise en 15 jours par une machine ? Nombre de personnes qui ont connu des procédures de plusieurs années, aboutissant à un non-lieu, rêvent de procédures accélérées. Car le « temps perdu » a une valeur, les actuaires le savent bien. Que penser alors de cette efficacité des modèles algorithmiques ? La justice se doit d’être efficace, mais cette contrainte ne doit pas faire oublier l’objectif central qui est celui de rendre la justice. Que se passe-t-il si l’efficacité devient un objectif, pour ne pas dire le principal objectif ? Car c’est bien la question que posent les modèles prédictifs : quel est l’objectif que l’on cherche à maximiser ? Et comment le formule-t-on de manière simple ? # Aide à la décision, ou justification d’une prise de décision Aux Etats-Unis, nombre de juges se sont vu reproché de motiver un jugement à l’aide d’outils d’aide à la décision, ce qui laisse planer un doute sur la fonction réelle de ces outils. L’idée était initialement d’apporter une aide. Récemment, plusieurs systèmes mis en place dans les années passées ont été questionnés. Par exemple en Australie, le STMP (Suspect Targeting Management Plan) proposait d’identifier si des préadolescents doivent être surveillés, ou pas. Ce modèle ressemble à s’y méprendre à n’importe quel modèle actuariel, c’est à dire un outil d’évaluation et prédiction des risques, en se concentrant soit sur les récidivistes, soit sur les personnes suspectées de commettre un futur crime. Or un rapport récent montrait que son utilisation n’a eu « aucun impact observable sur la prévention du crime »[2]. Parallèlement aux Etats-Unis, l’outil Compas (Correctional Offender Management Profiling Alternative Sanctions) a été critiqué dans Dressel & Farid (2017) : « Les défenseurs de ces systèmes soutiennent que les données et l’apprentissage automatique avancé rendent ces analyses plus précises et moins biaisées que celles des humains. Cependant, nous montrons que le logiciel d’évaluation des risques Compas, largement utilisé, n’est pas plus précis ou juste que des prédictions faites par des personnes qui ont peu ou pas du tout d’expertise en matière de justice pénale ». En questionnant des personnes recrutées sur Internet, sans compétences en droit, il s’agissait de prévoir si des personnes allaient, ou pas commettre un autre crime dans les deux ans à venir. Compas s’est trompé dans 34.8% des cas, et les internautes dans 33% des cas. Cela dit, on peut se demander ici ce que “se tromper” signifie. En l’occurrence, on ne mesure pas ici la récidive des personnes, mais la condamnation pour récidive des personnes. Et si les modèles (ou les gens) ne s’étaient pas trompés, mais que les juges, en revanche, oui ? # Prévoir, et se tromper Et si un des soucis ne venait pas dans ce que l’on demande à un outil prédictif ? Prévoir, c’est (fondamentalement) établir une probabilité pour un fait futur. Comme cela avait été rappelé dans un débat sur les sondages et les élections, est-ce qu’on peut dire qu’on se trompe si on annonce qu’un événement peut se produire avec 5% de chances, et qu’il se produit effectivement ? Pour savoir si une technique de prévision est bonne, il faut collecter un ensemble de prévision, et les comparer aux observations. C’est ce que font les météorologues depuis une quinzaine d’années, et qui a été formalisé par Gneiting et al. (2007). Leur idée est qu’un ensemble de distributions prédictives est obtenu par un modèle $\{\hat F_t,\hat F_{t+1},..,\hat F_{t+h}\}$ et il convient de comparer ces distributions aux observations $\{y_t,y_{t+1},..,y_{t+h}\}$ – et non pas $\{\hat y_t,\hat y_{t+1},..,\hat y_{t+h}\}$. Il faut alors introduire une distance entre les distributions prédictives, et les observations. Dans un système physique, il est possible d’imaginer comprendre les différentes relations causales, et donc de prévoir. Mais dans les relations humaines (et la justice en est un exemple parfait), rien n’est aussi simple, aussi automatique que les lois de mécanique des fluides qui permettent de modéliser des phénomènes météorologiques. # Références Binet, Jacques-Louis, 2010, La prévention médicale de la récidive chez les délinquants sexuels. Académie de Médecine. Botsman, Rachel. 2017. Who Can You Trust?: How Technology Brought Us Together – and Why It Could Drive Us Apart. Portfolio Penguin Charpentier, Arthur & Suire, Raphaël 2016. Données et santé: valeurs, acteurs et santé. Risques, 107 Charpentier, Arthur. 2016. Les dérives du principe de précaution. Risques. 108 Christin, Agnèle, Rosenblat, Alex & Boyd, Danah 2015. Courts and Predictive Algorithms. Datacivilrights Dressel, Julie & Farid, Hany 2018. The accuracy, fairness, and limits of predicting recidivism. Science Adavances Foucault, Michel 1975 Surveiller et punir, naissance de la prison. Gallimard Galeon, Dom & Bergan, Brad 2017. China’s “Social Credit System” Will Rate How Valuable You Are as a Human. Futurism Gneiting, Balabdaoui & Raftery 2007. Probabilistic forecasts, calibration and sharpness. JRRS-B, 69, 243–268. Mangel, Marc & Samaniego, Francisco 1984. Abraham Wald’s work on aircraft survivability , Journal of the American Statistical Association, vol. 79, no 386,‎ 259–267 McLannahan, 2015 Being ‘wasted’ on Facebook may damage your credit score (Octobre 2015, Financial Times) Supiot, Alain 2015. La gouvernance par les nombres : cours au Collège de France, 2012-2014. Fayard. Trujillo, Elsa 2017. La Chine met en place un système de notation de ses citoyens pour 2020. Le Figaro, décembre 2017 [1] Décrit dans « In Your Face: China’s all-seeing state » http://www.bbc.com/news/av/world-asia-china-42248056/in-your-face-china-s-all-seeing-state [2] https://www.numerama.com/politique/300907-un-algorithme-teste-par-la-police-pour-anticiper-les-crimes-des-jeunes-inquiete-laustralie.html # 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…