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 ! Acheter un billet de train (pas trop cher) Hier, je suis tombé sur article qui discutait des prix des billets de train, en France (et du prix très élevé, a certaines dates, genre pendant les vacances d’hiver). Tous ceux qui ont l’habitude de prendre le train savent que le prix que l’on paye dépend du moment ou on achète le billet (et de la souplesse que l’on peut avoir sur l’heure (voire la date) du trajet). Cet été, dans le cadre du projet de la formation Data Science pour l’Actuariat, pour mon cours, Pierre proposait de moissonner le site https://www.oui.sncf/ pour suivre un peu l’évolution du prix des billets. Le principal soucis est que le site https://www.oui.sncf/ s’appuie sur du javascript pour les formulaires de saisie et pour l’affichage des résultats, ce qui empêche l’utilisation classique du package rvest, par exemple. J’avais évoqué dans un autre billet l’utilisation de wdman, pour scraper le site des incendies de forets. Ici, Pierre proposait de passer par casperjs, et je vais reprendre un peu ici sa stratégie: • on va utiliser casperjs, un émulateur de navigateur écrit en javascript. Il permet d’émuler un véritable navigateur (même moteur que google chrome) et de résoudre le javascript intégrés à la page • on va utiliser un petit code en bash pour lancer le code Pour le code en bash, c’est juste que je suis sous mac et linux. Sous mac, on peut faire pas mal de choses en bash… Tout passe par des variables, que l’on peut définir, et afficher, par exemple qui nous donne l’heure. Je peux aussi définir une variable, et l’incrémenter (pratique pour faire des boucles) Plus intéressant, on peut planifier des taches. Pour ça, on tape qui va ouvrir un éditeur, Je demande ici de lancer un code R, tous les heures, a 13:50, 14:50, 15:50, etc. Pour demander tous les jours a 13:50, je tape on va ensuite sauver l’instruction On voit que la commande sera lancée, tous les jours : elle est dans la liste des taches a faire notons qu’on peut lancer un code tout en passant des arguments : ici je lui dis quel objet manipuler (le premier argument) et le second sert pour créer un fichier (et pour le nommer). Bref, c’est assez facile de lancer automatiquement des codes, pour scraper. Sous windows, on passe par le planificateur de taches. Un premier script permet d’extraire les trains ouverts dans la journée, et le second, les trains ouverts depuis plus de 24 heures, Ensuite, on va utiliser http://docs.casperjs.org/en/latest/ pour coder notre émulateur de navigateur internet (le code est ici en ligne). On va ainsi créer plein de fichiers, contenant les informations que l’on veut ! Je vais passer un peu le retraitement, et juste présenter les informations qu’on peut en tirer. En particulier, Pierre avait stocke les données entre mars et juin dernier. Ici, on est juste sur quelques trajets de grande ligne, library(readr) library(rgdal) nomFichier = tempfile(fileext = ".zip") download.file("https://freakonometrics.free.fr/CarteFrance.zip", destfile = nomFichier, mode = "wb") unzip(zipfile = nomFichier, exdir = getwd()) download.file("https://freakonometrics.free.fr/LgTroncons.csv", destfile = "LgTroncons.csv", mode = "wb") download.file("https://freakonometrics.free.fr/CoordVilles.csv", destfile = "CoordVilles.csv", mode = "wb") fra0 = readOGR(dsn = paste(getwd(), "/CarteFrance", sep = ""), layer = "gadm36_FRA_0", verbose = F) LgTroncons = read_delim("LgTroncons.csv",";", escape_double = FALSE,locale = locale(decimal_mark = ","), trim_ws = TRUE) CoordVilles = read_delim("CoordVilles.csv",";", escape_double = FALSE,locale = locale(decimal_mark = ","), trim_ws = TRUE) NomsVilles = CoordVilles[CoordVilles$NOM_A_AFFICHER==1,] library(ggplot2) fr_df = fortify(fra0) ggp = ggplot() + geom_polygon(data=fr_df, aes(long, lat,group = group), fill = "#3A8EBA") ggp = ggp + geom_path(data = LgTroncons, aes(x = LONG, y = LAT, group = ID_TRONCON), colour= "#CC5500", lineend = "round", size=3) + geom_path(data = LgTroncons, aes(x=LONG, y=LAT, group = ID_TRONCON), colour="white", lineend = "round", size=1.75) ggp = ggp + geom_point(data = NomsVilles, aes(x=LONG, y=LAT), colour = "blue", fill = "white", shape=21, size = NomsVilles$PT_SIZE, stroke = NomsVilles$PT_STROKE) + theme_void() ggp = ggp + geom_text(data = NomsVilles, aes(x=LONG, y=LAT, label=NOM),hjust = NomsVilles$H_AJUST, vjust=NomsVilles$V_AJUST, colour = "white", fontface = "bold", size =3.25)+coord_fixed(1.47) ggp = ggp + ggtitle("Représentation des trajets étudiées") + theme(plot.title = element_text(hjust = 0.5, face="bold")) print(ggp)

On va travailler sur les trajets suivants,

Comparons ici les billets pour un trajet Pars-Rennes, a partir des informations moissonnées pendant 3 mois, le vendredi soir, en particulier 2 vendredi de juin 2018 (les 15 et 22 juin). Pour ces deux jours, il y a eu 6 trains, entre 17 et 20 heures. Pour le 15 juin, les trois premiers ont commencé avec un prix de 45 €. Pour le 22 juin, le premier a commence a 45 €, mais les deux suivant on été lancés a 33 €. Assez rapidement, les prix sont monte a 45 €.

On peut regarder l’evolution du prix

Si on regarde plusieurs destination, on observe des comportement très différentes,

• pour Le Mans, les prix montent très vite, commençant a 15 €, montant a 18 € 10 heures après le lancement, 21 € le lendemdemain, 27 € au bout d’une semaine. En un mois, les prix ont presque doublé.
• pour Rennes, on observe une evolution similaire, passant de 20 € a 25 € en quelques heures, et 40 € deux semaines après !
• pour Toulouse en revanche, le prix initial est plus haut, 43 €, monte de 3 € en 10 heures, 6 € en 16 heures, pour rester a 49 €

Mais pour toutes les destinations, les prix sont croissants.

soit graphiquement

On peut aussi faire une carte. Si on regarde les prix a l’ouverture, Lille, Le Mans et Rennes sont peu chères.

Et les plus fortes variations sur 10 heures sont observées sur Nantes et Bordeaux.

Amusant non ?

Monte Carlo techniques to create counterfactuals

In the previous STT5100 course, last week, we’ve seen how to use monte carlo simulations. The idea is that we do observe in statistics a sample $\{y_1,\cdots,y_n\}$, and more generally, in econometrics $\{(y_1,\mathbf{x}_1),\cdots,(y_n,\mathbf{x}_n)\}$. But let’s get back to statistics (without covariates) to illustrate. We assume that observations $y_i$ are realizations of an underlying random variable $Y_i$. We assume that $Y_i$ are i.id. random variables, with (unkown) distribution $F_{\theta}$. Consider here some estimator $\widehat{\theta}$ – which is just a function of our sample $\widehat{\theta}=h(y_1,\cdots,y_n)$. So $\widehat{\theta}$ is a real-valued number like . Then, in mathematical statistics, in order to derive properties of the estimator $\widehat{\theta}$, like a confidence interval, we must define $\widehat{\theta}=h(Y_1,\cdots,Y_n)$, so that now, $\widehat{\theta}$ is a real-valued random variable. What is puzzling for students, is that we use the same notation, and I have to agree, that’s not very clever. So now, $\widehat{\theta}$ is .

There are two strategies here. In classical statistics, we use probability theorem, to derive properties of $\widehat{\theta}$ (the random variable) : at least the first two moments, but if possible the distribution. An alternative is to go for computational statistics. We have only one sample, $\{y_1,\cdots,y_n\}$, and that’s a pity. But maybe we can create another one $\{y_1^{(1)},\cdots,y_n^{(1)}\}$, as realizations of $F_{\theta}$, and another one $\{y_1^{(2)},\cdots,y_n^{(2)}\}$, anoter one $\{y_1^{(3)},\cdots,y_n^{(3)}\}$, etc. From those counterfactuals, we can now get a collection of estimators, $\widehat{\theta}^{(1)}$,$\widehat{\theta}^{(2)}$, $\widehat{\theta}^{(3)}$, etc. Instead of using mathematical tricks to calculate $\mathbb{E}(\widehat{\theta})$, compute $$\frac{1}{k}\sum_{s=1}^k\widehat{\theta}^{(s)}$$That’s what we’ve seen last friday.

I did also mention briefly that looking at densities is lovely, but not very useful to assess goodness of fit, to test for normality, for instance. In this post, I just wanted to illustrate this point. And actually, creating counterfactuals can we a good way to see it. Consider here the height of male students,

Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] X=Davis$height[Davis$sex=="M"]

We can visualize its distribution (density and cumulative distribution)

u=seq(155,205,by=.5) par(mfrow=c(1,2)) hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) lines(u,dnorm(u,178,6.5),col="black") Xs=sort(X) n=length(X) p=(1:n)/(n+1) plot(Xs,p,type="s",col="blue") lines(u,pnorm(u,178,6.5),col="black")

Since it looks like a normal distribution, we can add the density a Gaussian distribution on the left, and the cdf on the right. Why not test it properly. To be a little bit more specific, I do not want to test if it’s a Gaussian distribution, but if it’s a $\mathcal{N}(178,6.5^2)$. In order to see if this distribution is relevant, one can use monte carlo simulations to create conterfactuals

hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) Y=rnorm(n,178,6.5) hist(Y,col=rgb(1,0,0,.3)) lines(density(Y),col="red",lwd=2) Ys=sort(Y) plot(Xs,p,type="s",col="white",lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(155,205)) polygon(c(Xs,rev(Ys)),c(p,rev(p)),col="yellow",border=NA) lines(Xs,p,type="s",col="blue",lwd=2) lines(Ys,p,type="s",col="red",lwd=2)

We can see on the left that it is hard to assess normality from the density (histogram and also kernel based density estimator). One can hardly think of a valid distance, between two densities. But if we look at graph on the right, we can compare the empirical distribution cumulative distribution $\widehat{F}$ obtained from $\{y_1,\cdots,y_n\}$ (the blue curve), and some conterfactual, $\widehat{F}^{(s)}$ obtained from $\{y_1^{(s)},\cdots,y_n^{(s)}\}$ generated from $F_{\theta_0}$ – where $\theta_0$ is the value we want to test. As suggested above, we can compute the yellow area, as suggest in Cramer-von Mises test, or the Kolmogorov-Smirnov distance.

d=rep(NA,1e5) for(s in 1:1e5){ d[s]=ks.test(rnorm(n,178,6.5),"pnorm",178,6.5)$statistic } ds=density(d) plot(ds,xlab="",ylab="") dks=ks.test(X,"pnorm",178,6.5)$statistic id=which(ds$x&gt;dks) polygon(c(ds$x[id],rev(ds$x[id])),c(ds$y[id],rep(0,length(id))),col=rgb(1,0,0,.4),border=NA) abline(v=dks,col="red")

If we draw 10,000 counterfactual samples, we can visualize the distribution (here the density) of the distance used a test statistic $\widehat{d}^{(1)}$, $\widehat{d}^{(2)}$, etc, and compare it with the one observe on our sample $\widehat{d}$. The proportion of samples where the test-statistics exceeded the one observed

mean(d&gt;dks) [1] 0.78248

is the computational version of the $p$-value

ks.test(X,"pnorm",178,6.5)   One-sample Kolmogorov-Smirnov test   data: X D = 0.068182, p-value = 0.8079 alternative hypothesis: two-sided

I thought about all that a couple of days ago, since I got invited for a panel discussion on “coding”, and why “coding” helped me as professor. And this is precisely why I like coding : in statistics, either manipulate abstract objects, like random variables, or you actually use some lines of code to create counterfactuals, and generate fake samples, to quantify uncertainty. The later is interesting, because it helps to visualize complex quantifies. I do not claim that maths is useless, but coding is really nice, as a starting point, to understand what we talk about (which can be very usefull when there is a lot of confusion on notations).

October, grant proposal season

In 2012, Danielle Herbert, Adrian Barnett, Philip Clarke and Nicholas Graves published an article entitled “on the time spent preparing grant proposals: an observational study of Australian researchers“, whose conclusions had been included in Nature under a more explicit title, “Australia’s grant system wastes time” ! In this study, they included 3700 grant applications sent to the National Health and Medical Research Council, and showed that each application represented 37 working days: “Extrapolating this to all 3,727 submitted proposals gives an estimated 550 working years of researchers’ time (95% confidence interval, 513-589)“. But in these times when I have to write my funding application, I find that losing 37 days of work is huge. Because it’s become the norm! And somehow, it’s sad.

Forget about the crazy idea that I would rather, in fact, spend more time doing my research. In fact, the thought I had this morning was that it is rather sad that in the Faculty of Science, mathematicians are asked to spend a considerable amount of time, comparable to that required of physicists or chemists, for often smaller amounts of funding… And I thought it could be easily verified. We start by retrieving the discipline codes

url="http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSC-ResultatsCSS_eng.asp" download.file(url,destfile = "GSC.html") library(XML) tables=readHTMLTable("GSC.html") GSC=tables[[1]]$V1 GSC=as.character(GSC[-(1:2)]) namesGSC=tables[[1]]$V2 namesGSC=as.character(namesGSC[-(1:2)])

We’re going to need a small function, to remove the $and other symbols that pollute the data (and prevent them from being treated as numbers) library(stringr) Correction = function(x) as.numeric(gsub('[$,]', '', x))

We will now read the 12 pages, and harvest (we will just take the 2017 data, but we could go back a few years before)

grants= function(gsc){ url=paste("http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSCDetail-ResultatsCSSDetails_eng.asp?Year=2017&amp;GSC=",gsc,sep="") download.file(url,destfile = "GSC.html") library(XML) tables=readHTMLTable("GSC.html") X=as.character(tables[[1]]$"Awarded Amount") A=as.numeric(Vectorize(Correction)(X)) return(c(median(A),mean(A),as.numeric(quantile(A,(1:99)/100)))) } M=Vectorize(grants)(GSC[1:12]) The average amounts of individual grants can be compared, barplot(M[2,]) In mathematics, the average grant amount is$24400. If we normalize by this quantity, we obtain

barplot(M[2,]/M[2,8])

In other words, the average amount of a (individual) grant in chemistry (to pay for students, conferences, etc.) is twice that in mathematics, 60% higher in physics than in maths…

We can also look at the median values (rather than the averages)

barplot(M[1,])

Here again, it is in mathematics that it is the weakest….

barplot(M[1,]/M[1,8])

in comparable proportions. If we think that the time spent writing should be proportional to the amount allocated, we should spend half as much time in math as in chemistry.

Cumulative functions can also be ploted,

plot(M[3:101,8],(1:99)/100,type="s",xlim=range(M)) lines(M[3:101,5],(1:99)/100,type="s",col="red") lines(M[3:101,4],(1:99)/100,type="s",col="blue")

with math in black, physics in red, and chemistry in blue. What is surprising is the bottom part: a “bad” researcher in chemistry or physics will earn more than the median researcher in mathematics…

Now that my intuition is confirmed, I have to go back, writing my proposal… and explain to my coauthors that I have to postpone some research projects because, well, you know…

Combining automatically factor levels in R

Each time we face real applications in an applied econometrics course, we have to deal with categorial variables. And the same question arise, from students : how can we combine automatically factor levels ? Is there a simple R function ?

I did upload a few blog posts, over the pas years. But so far, nothing satistfying. Let me write down a few lines about what could be done. And if some wants to write a nice R function, that would be awesome. To illustrate the idea, consider the following (simulated dataset)

n=200 set.seed(1) x1=runif(n) x2=runif(n) y=1+2*x1-x2+rnorm(n,0,.2) LB=sample(LETTERS[1:10]) b=data.frame(y=y,x1=x1, x2=cut(x2,breaks= c(-1,.05,.1,.2,.35,.4,.55,.65,.8,.9,2), labels=LB)) str(b) 'data.frame': 200 obs. of 3 variables: $y : num 1.345 1.863 1.946 2.481 0.765 ...$ x1: num 0.266 0.372 0.573 0.908 0.202 ... $x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ... table(b$x2)[LETTERS[1:10]]   A B C D E F G H I J 11 12 23 34 23 36 12 32 3 14

There is one (continuous) dependent variable $y$, one continuous covariable $x_1$ and one categorical variable $x_2$, with here ten levels. We can plot the data using

plot(b$x1,y,col="white",xlim=c(0,1.1)) text(b$x1,y,as.character(b$x2),cex=.5) The output of a linear regression yield the following predictions for(i in 1:10){ p=function(x) predict(lm(y~x1+x2,data=b),newdata=data.frame(x1=x,x2=LETTERS[i])) u=seq(-1,1.065,by=.01) v=Vectorize(p)(u) lines(u,v)} the slope for $x_1$ is the same, we simply add a different constant for each level. As we can see, some levels are very very close, so it seems legitimate to combine them into one single category. Here is the output of the linear regression, summary(lm(y~x1+x2,data=b)) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 0.843802 0.119655 7.052 3.23e-11 *** x1 1.992878 0.053838 37.016 &lt; 2e-16 *** x2A 0.055500 0.131173 0.423 0.6727 x2H 0.009293 0.121626 0.076 0.9392 x2F -0.177002 0.121020 -1.463 0.1452 x2B -0.218152 0.130192 -1.676 0.0955 . x2D -0.206970 0.121294 -1.706 0.0896 . x2G -0.407417 0.129999 -3.134 0.0020 ** x2C -0.526708 0.123690 -4.258 3.24e-05 *** x2J -0.664281 0.128126 -5.185 5.54e-07 *** x2E -0.816454 0.123625 -6.604 3.94e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2014 on 189 degrees of freedom Multiple R-squared: 0.8995, Adjusted R-squared: 0.8942 F-statistic: 169.1 on 10 and 189 DF, p-value: &lt; 2.2e-16 AIC(lm(y~x1+x2,data=b)) [1] -60.74443 BIC(lm(y~x1+x2,data=b)) [1] -21.16463 Here the reference category is “I”. And it looks like we could actually combine that category with several others. One strategy here would be to select all categories that seem to be not significantly different, and to run a (multiple) test library(car) linearHypothesis(lm(y~x1+x2,data=b), c("x2A = 0", "x2H = 0", "x2F = 0")) Hypothesis: x2A = 0 x2H = 0 x2F = 0 Model 1: restricted model Model 2: y ~ x1 + x2 Res.Df RSS Df Sum of Sq F Pr(&gt;F) 1 192 8.4651 2 189 7.6654 3 0.79971 6.5726 3e-04 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 It seems that we can combine those four categories together. Here, we can see what’s going on when we change the reference category (actually, loop on all categories) P=matrix(NA,nlevels(b$x2),nlevels(b$x2)) colnames(P)=rownames(P)=LETTERS[1:10] plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5), ylim=c(0,10.5)) text(1:10,0,LETTERS[1:10]) text(0,1:10,LETTERS[1:10]) for(i in 1:nlevels(b$x2)){ #levels(b$x2)=LETTERS[1:10] b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,3) P[LETTERS[i],names(p)]=p p=P[LETTERS[i],] idx=which(p&gt;.05) points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2) idx=which(p&gt;.1) points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)}

We are glad to see that it is symmetric : if “H” should be combined with “I”, “I” should also be combined with “H”.

Here black points are related with the 10% $p$-value, and white points the 5% $p$-value. This graph is actually hard to read… And actually, this reminds us of  Bertin (1967).

Here, we can predefine manually some ordering (we will see below how it might be automatised)

LETTERSord=c("I","A","H","F","B","D","G","C","J","E") P=matrix(NA,nlevels(b$x2),nlevels(b$x2)) colnames(P)=rownames(P)=LETTERSord plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5), ylim=c(0,10.5)) ct=c(3,3,2,1,1) abline(v=.5+c(0,cumsum(ct)),lty=2) abline(h=.5+c(0,cumsum(ct)),lty=2) text(1:10,0,LETTERSord) text(0,1:10,LETTERSord) for(i in 1:nlevels(b$x2)){ #levels(b$x2)=LETTERS[1:10] b$x2=relevel(b$x2,LETTERSord[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,3) P[LETTERSord[i],names(p)]=p p=P[LETTERSord[i],] idx=which(p&gt;.05) points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2) idx=which(p&gt;.1) points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2) } Here we get the following It looks like we have our combined categories… Actually, it is possible to use another strategy. We start from some level, say “A”. Then, we merge it with all non-significantly different levels. If “B” is not one of them, we use it as the new reference. Etc. for(i in 1:nlevels(b$x2)){ if(LETTERS[i]%in%levels(b$x2)){ b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,nchar(p)) idx=which(p&gt;.05) mix=c(LETTERS[i],names(p)[idx]) b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep="")) }}

The final categories are

table(b$x2) A+I+H B+D+F C+G E J 46 82 35 23 14 with the following regression output summary(lm(y~x1+x2,data=b)) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 0.86407 0.03950 21.877 &lt; 2e-16 *** x1 1.99180 0.05323 37.417 &lt; 2e-16 *** x2B+D+F -0.21517 0.03699 -5.817 2.44e-08 *** x2C+G -0.50545 0.04528 -11.164 &lt; 2e-16 *** x2E -0.83617 0.05128 -16.305 &lt; 2e-16 *** x2J -0.68398 0.06131 -11.156 &lt; 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2008 on 194 degrees of freedom Multiple R-squared: 0.8975, Adjusted R-squared: 0.8948 F-statistic: 339.6 on 5 and 194 DF, p-value: &lt; 2.2e-16 AIC(lm(y~x1+x2,data=b)) [1] -66.76939 BIC(lm(y~x1+x2,data=b)) [1] -43.68117 Which is consistent with the group we got before. But actually, if we change the order, we can get different combinations. For instance, if we go from “J” to “A”, instead of “A” to “J”, we obtain for(i in nlevels(b$x2):1){ #levels(b$x2)=LETTERS[1:10] if(LETTERS[i]%in%levels(b$x2)){ b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,nchar(p)) idx=which(p&gt;.05) mix=c(LETTERS[i],names(p)[idx]) b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep="")) }} table(b$x2)   E G+C I+A+B+D+F+H J 23 35 128 14

with different information criteria here

AIC(lm(y~x1+x2,data=b)) [1] -36.61665 BIC(lm(y~x1+x2,data=b)) [1] -16.82675

I guess it would be necessary to run randomly the order we go through the levels. Last, but not least, one can use regression trees (even if it not per se in the syllabus of the course). The problem is that there is another explanatory variable that might interphere. So I would suggest (1) to fit a linear model $y=\beta_0+\beta_1x_1+u_i$, to calculate the residuals, $\widehat{u}_i$ (2) to run a regression tree, to explain $\widehat{u}_i$ with categorical variable $x_2$ (I did explain how trees are build when the explanatory variable is a categorical one in a previous post)

library(rpart) library(rpart.plot) b$e=residuals(lm(y~x1,data=b)) arbre=rpart(e~x2,data=b) prp(arbre,type=2,extra=1) Observe that the leaves have the same groups as the one we got. arbre n= 200 node), split, n, deviance, yval * denotes terminal node 1) root 200 22.563500 7.771561e-18 2) x2=G,C,J,E 72 4.441495 -3.232525e-01 4) x2=J,E 37 1.553520 -4.578492e-01 * 5) x2=G,C 35 1.509068 -1.809646e-01 * 3) x2=I,A,H,F,B,D 128 6.366628 1.818295e-01 6) x2=F,B,D 82 2.983381 1.048246e-01 * 7) x2=I,A,H 46 2.030229 3.190993e-01 * I guess that it should be possible to put all that in an R function, to suggest combinations of level that might improve the regression. Traitement des valeurs manquantes : remplacer les NA par une constante ? Un rapide billet pour répondre à une question posée à la fin du cours de ce matin (ST5100), par Jean-Pierre Liégeois, jeune lecteur du Var (pour préserver un peu d’anonymat) Dans mon stage, quand on avait des valeurs manquantes, on me disait de remplacer par -1, puis de rajouter une indicatrice comme quoi la variable vaut -1. Ça permet de ne supprimer ni variable, ni observations. On peut faire ça ? Si je formalise un peu, on va simuler ici des données, disons $x_1$ et $x_2$, on génère ensuite des données suivant un modèle, de la forme $y=\beta_0+\beta_1x_1+\beta_2x_2+\varepsilon$. Une proportion $\alpha$ de $x_1$ sera transformée en NA. Ce que suggérais Jean-Pierre, c’est de remplacer les valeurs manquantes par -1, puis d’ajuster un modèle $$y=\beta_0+\beta_1x_1+\beta_{-1}\mathbf{1}(x_1=-1)+\beta_2x_2+\varepsilon$$. Côté code, c’est assez simple. Par défaut, la stratégie de R est de supprimer les valeurs manquantes. Si 50% des données de $x_1$ sont manquantes, la moitié des lignes sont supprimées n=1000 x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.05 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1) base$x1[indice]=NA reg=lm(y~x1+x2,data=base)

Au lieu de générer un unique échantillon, on va en simuler 10,000, et regarder la distribution de $\widehat{\beta}_1$,

m=10000 B=rep(NA,m) for(s in 1:m){ x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.5 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA reg=lm(y~x1+x2,data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white",xlab="missing values = 50%") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") Bien sur, avec un taux de valeurs manquantes plus faibles – disons $\alpha=5\%$ – on perd moins d’observations, et donc l’estimateur a une variance plus faible. Tentons maintenant la stratégie consistant à remplacer les valeurs manquantes par des valeurs numériques fixes, et de rajouter une indicatrice, B=rep(NA,m) for(s in 1:m){ x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.5 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg=lm(y~x1+x2+I(x1==(-1)),data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red")

Ce qui ne change pas grand chose, on en conviendra…  y compris si le taux de valeurs manquantes passe à 5%,

On peut se demander ce qui se passe si le shift n’est plus de 1 mais de 10 (a priori, c’est arbitraire, la variable $x_1$  pouvant être plus ou moins dispersée… -1 pour une variable entre 0 et 1, ou entre 0 et 1000, ça n’est pas tout à fait pareil). Mais non, par exemple avec toujours 5% de valeurs manquantes, on a

Si on regarde notre échantillon, en particulier le nuage de points  $(x_1,y)$, on observe

ici, les valeurs manquantes sont choisies au hasard, de manière totalement indépendante,

x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.3333333 indice=sample(1:n,size=round(n*alpha)) clr=rep("black",n) clr[indice]="red" plot(x1,y,col=clr)

(ici avec 1/3 de valeurs manquantes, en rouge). Mais on pourrait supposer que les valeurs manquantes sont les plus grandes valeurs de $x_1$, par exemple,

x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.3333333 indice=sample(1:n,size=round(n*alpha),prob = x1^3) clr=rep("black",n) clr[indice]="red" plot(x1,y,col=clr)

On peut se demander ce que ça donnerait sur l’estimateur $\widehat{\beta}_1$

Ça ne change pas grand chose, mais on a plus de variance, si on regarde bien. Dernier essai : que se passe-t-il si les variables $x_1$ et $x_2$ sont maintenant corrélées,

B=rep(NA,m) library(mnormt) r=.8 S = matrix(c(1,r,r,1),2,2) for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg=lm(y~x1+x2+I(x1==(-1)),data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") Cette fois, on a un estimateur biaisé (de l’ordre de 10% sur cet exemple numérique). Manifestement, cette technique n’est pas très concluante… Je pourrais ajouter que cette méthode ne revient pas à la première, même si la distribution des estimateurs est proche set.seed(1) x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg1=lm(y~x1+x2+I(x1==(-1)),data=base) coefficients(reg1) (Intercept) x1 x2 I(x1 == (-1))TRUE 1.0988005 1.7454385 -0.5149477 3.1000668 base$x1[indice]=NA reg2=lm(y~x1+x2,data=base) coefficients(reg2) (Intercept) x1 x2 1.1123953 1.8612882 -0.6548206 Comme je le disais (lors de la discussion qui a suivi le cours) une méthode plus prometteuse est l’imputation. L’idée est de prédire une valeur pour les $x_1$ manquants. On pourrait être tenté de prendre $\overline{x}_1$, la moyenne des $x_1$ observés. Mais on sait que les valeurs manquantes sont justement pour les grandes valeurs de $x_1$, ici, donc on doit pouvoir faire mieux ! On sait aussi que $x_1$ et $x_2$ sont corrélés ici. Corrélés positivement, en plus. Autrement dit, si $x_2$ est grand, on sait que le $x_1$ (non observé) devait être grand. Le plus simple est de faire un modèle linéaire, $x_1=\alpha_0+\alpha_2x_2+\eta_i$, calibré sur les valeurs non-manquantes. Puis on utilise $\widehat{x}_1=\widehat{\alpha}_0+\widehat{\alpha}_2x_2$ pour les valeurs manquantes. C’est simpliste, mais pourquoi pas ? On estime alors le modèle sur cette nouvelle base. for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA reg0=lm(x1~x2,data=base[-indice,]) base$x1[indice]=predict(reg0,newdata=base[indice,]) reg=lm(y~x1+x2,data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") sur l’exemple numérique, on obtient base$x1[indice]=NA reg0=lm(x1~x2,data=base[-indice,]) base$x1[indice]=predict(reg0,newdata=base[indice,]) reg3=lm(y~x1+x2,data=base) coefficients(reg3) (Intercept) x1 x2 1.1593298 1.8612882 -0.6320339 Cette méthode a au moins pu corriger du biais… Après, si on regarde attentivement, on a exactement la même valeur qu’avec la première méthode qui consiste à supprimer les lignes avec des valeurs manquantes ! summary(reg3) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 1.15933 0.06649 17.435 &lt; 2e-16 *** x1 1.86129 0.21967 8.473 &lt; 2e-16 *** x2 -0.63203 0.20148 -3.137 0.00176 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.051 on 997 degrees of freedom Multiple R-squared: 0.1094, Adjusted R-squared: 0.1076 F-statistic: 61.23 on 2 and 997 DF, p-value: &lt; 2.2e-16 summary(reg2) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 1.11240 0.06878 16.173 &lt; 2e-16 *** x1 1.86129 0.21666 8.591 &lt; 2e-16 *** x2 -0.65482 0.20820 -3.145 0.00172 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.037 on 797 degrees of freedom (200 observations deleted due to missingness) Multiple R-squared: 0.1223, Adjusted R-squared: 0.12 F-statistic: 55.5 on 2 and 797 DF, p-value: &lt; 2.2e-16 Au lieu de faire une régression linéaire, on peut utiliser une autre méthode d’imputation. Par exemple prendre la moyenne des $k$ valeurs de $x_1$ (observées) pour les $k$ individus ayant des $x_2$ les plus proches du $x_2$ de l’individu ayant $x_1$ manquant : kpp=function(i,basena,k=5){ x2=basena$x2[i] sb=basena[!is.na(basena$x1),] idx=rank(abs(sb$x2-x2)) mean(sb[which(idx&lt;=k),"x1"]) }

Sur notre base simulée on obtient

base$x1[indice]=NA base0=base for(j in indice) base0$x1[j]=kpp(j,base0,k=5) reg4=lm(y~x1+x2,data=base) coefficients(reg4) (Intercept) x1 x2 1.197944 1.804220 -0.806766

Si on regarde ce que ça donne sur nos 10,000 simulations, on a (c’est un peu long, car j’ai codé ça très rapidement, et pas du tout de manière optimale)

for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA base0=base for(j in indice) base0$x1[j]=kpp(j,base0,k=5) reg=lm(y~x1+x2,data=base0) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red")

Le biais semble ici plus faible que sans imputation… autrement dit, les méthodes d’imputation me semblent plus robustes que cette stratégie visant à remplacer des NA par une valeur arbitraire, et de rajouter une indicatrice dans la régression.