Category Archives: Statistics

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
X = rmnorm(n,rep(0,3),Sigma)
df = data.frame(X1=X[,1],X2=X[,2],X3=X[,3],X4=rnorm(n),
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,
VLASSO[[s]] = as.matrix(lasso$beta)

To visualize confidence bands, one can compute quantiles

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


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,

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("",header=TRUE,sep=";")
   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,

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

Now, consider bootstrap samples.

for(s in 1:100){

We get here

The interpretation here is much more difficult

What about the order ?

for(s in 1:100000){
        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,

g=function(y) mean(apply(N,1,function(x) paste(x,collapse="")==y))

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}\rbracewhere\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
X = rmnorm(n,rep(0,3),Sigma)
df = data.frame(X1=X[,1],X2=X[,2],X3=X[,3],X4=rnorm(n),
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

colrs = c(brewer.pal(8,"Set1"))[c(1,4,5,2,6,3,7,8)]

The first regression we can run is a non-standardized one

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))<200)

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

idx = which(apply(lasso$beta,1,function(x) sum(x==0))<length(vlambda))

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\rbracethen \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\rbracewhere\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,


To get a convex regression, use

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

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))

Let us add vertical lines for the knots


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>=s)

then we can use it direcly in our regression model

reg = glm(PRONO~INSYS+pos(INSYS,15)+

The output of the regression is here

               Estimate Std. Error z value Pr(>|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")

Using bs() linear splines

Using the GAM function, things are slightly different. We will use here so called b-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)

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),
              Estimate Std. Error z value Pr(>|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…


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>=s)
reg = glm(PRONO~poly(INSYS,2)+pos2(INSYS,15)+pos2(INSYS,25),
                Estimate Std. Error z value Pr(>|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")

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)

If we run R code, we get

reg = glm(PRONO~bs(INSYS,knots=c(15,25),
               Estimate Std. Error z value Pr(>|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")

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


and the prediction will now be

reg = glm(PRONO~bs(INSYS,knots=c(15,25),
u = seq(5,55,length=201)
v = predict(reg,newdata=data.frame(INSYS=u),type="response")

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

   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")

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)

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")

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")

Additive models

Consider now the second dataset, with two variables. Consider here a model like
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)
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…


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)
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

 (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", 
v_beta = t(Vectorize(simu)(1:1000))

Ooops. There is a problem here. Clearly, we cannot rely on numerical optimization here. We can think about using another optimization routine

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]]
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

logit = function(x){1/(1+exp(-x))}
logLik = function(beta, X, y){
 -sum(y*log(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)

With our OLS starting point, we obtain

opt0 = optim_second(beta_init,500)
[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)
[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})

 for(s in 1:9){

Observe that here, I use only ten iterations of the algorithm !

                [,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

 for(s in 1:9){
                [,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\rbraceThe 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

  (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)))
              Estimate Std. Error t value Pr(>|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")))
              Estimate Std. Error z value Pr(>|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)
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))

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)
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)
contour(u,u,v,levels = .5,add=TRUE)
reg = glm(y~x1+x2,data=df,family=binomial)

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("",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, NaytaData (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 (NaytaData kindly sent them to me tonight).

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.


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 »


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.

Sigma=matrix(c(1,r,r,1), 2, 2)

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

idx=which(substr(names(reg1$coefficients), 1,2)=="X1")

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)


My main concern here is the “reference” that is considered. Should A be the reference? Why not B


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


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


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.


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)


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.


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


Since we focus on residuals, those are centered, and we have an easy interpretation of respective values


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

for(k in 1:5){

We can plot that sequence of values, including some confidence intervals (that would be related to significance with respect to all other regions)


Looking at values does not give intuitive results, but I have the feeling that it is easy to explain what we plot (we compare each region to “the rest of the world”), and the ordering of a seems to be consistent with \alpha (but I could not prove it).

Here are some ideas I got. I should be able to provide other graphs, but I would love to discuss with anyone on that topics, to find a proper and nice way to visualize effects of a categorical explanatory variable in a regression model (that can be a logistic one). Comments are open…

Holt-Winters with a Quantile Loss Function

Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord & Snyder (2008)),

when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.

Furthermore, I like it because I think it has nice pedagogical features. Consider simple exponential smoothing, L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1} where \alpha\in(0,1) is the smoothing weight. It is locally constant, in the sense that {}_{t}\hat Y_{t+h} = L_{t}

 SimpleSmooth = function(a){
  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}

When using the standard R function, we get

hw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1])
[1] 0.2465579

Of course, one can replicate that optimal value

     for(t in 2:T){
         erreur[t]=X[t]-L[t-1] }
[1] 0.2464844

Here, the optimal value for \alpha is the one that minimizes the one-step prediction, for the \ell_2 loss function, i.e. \sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2 where here {}_{t-1}\hat Y_t = L_{t-1}. But one can consider another loss function, for instance the quantile loss function, \ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0}). The optimal coefficient is then obtained using

loss=function(e) e*(tau-(e<=0)*1)
  for(t in 2:T){
  erreur[t]=X[t]-L[t-1] }

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

abline(h= hw$alpha,lty=2,col="red")

Note that the optimal \alpha is decreasing with \tau. I wonder how general this result can be…

Of course, one can consider more general exponential smoothing, for instance the double one, with L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}]andB_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}so that the prediction is now {}_{t}\hat Y_{t+h} = L_{t}+hB_t (it is now locally linear – and no longer constant).


The code to compute the smoothed series is the following

DoubleSmooth = function(a,b){
  L[1]=X[1]; B[1]=0
  for(t in 2:T){
  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }

Here also it is possible to replicate R using the \ell_2 loss function

     L[1]=X[1]; B[1]=X[2]-X[1]
     for(t in 2:T){
         erreur[t]=X[t]-(L[t-1]+B[t-1]) }
[1] 0.41904510 0.05988304

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

loss=function(e) e*(tau-(e<=0)*1)
  L[1]=X[1]; B[1]=X[2]-X[1]
  for(t in 2:T){
   erreur[t]=X[t]-(L[t-1]+B[t-1]) }

and we can plot those values on a graph

abline(v= hw$alpha,lwd=.4,lty=2,col="red")
abline(h= hw$beta,lwd=.4,lty=2,col="red")

(with \alpha on the x-axis, and \beta on the y-axis). So here, it is extremely simple to change the loss function, but so far, it should be done manually. Of course, one do it also for the seasonal exponential smoothing model.

Justice “actuarielle”, algorithmes… et données

Il y a un peu plus d’un an, Virginie Gautron m’envoyait plein de documents sur la “justice actuarielle”, concept que je découvrais alors. Pour comprendre un peu de quoi il s’agit, je peux renvoyer vers pénalité et gestion des risques : vers une justice « actuarielle » en Europe ? qui dresse un état de l’art, en francais. J’avoue avoir un peu mis de côté ensuite (faute de temps), et puis l’autre jour, j’ai dévoré surveillé et punir de Nicolas Bourgoin, qui contient un chapitre sur le sujet.

Sans vouloir tomber dans le cliché des deux cultures (à la CP Snow), il faut bien admettre que pour un économètre, lire des écrits de juristes parlait de modèles prédictifs (ou actuariels) est une expérience marquante. Et ce qui me trouble (j’ai d’ailleurs le même soucis avec la psychologie – c’est un sujet proche) c’est que l’accent est souvent mis sur le modèle, l’algorithme, ce truc qui semble fasciner les juristes, alors qu’il me semble que le coeur du problème, ce sont les données. On nous dit, par exemple, que les algorithmes sont racistes (c’est un peu la thèse de Cathy O’Neil dans Weapons of Math Destruction), mais le problème ce n’est pas le modèle. C’est le fait qu’il y a un énorme biais de sélection dans les données ! Et selon l’adage “garbage in, garbage out“… dit de manière plus policée, je pourrais tenter un “purgamentum init, exit purgamentum” pour parler comme un juriste. Ou dit autrement, si on entre des mauvaises données en entrée, faut pas espérer grand chose en sortie…

Pour comprendre le biais de sélection, on peut relire Mostly Harmless Econometrics de Joshua Angrist et Jörn-Steffen Pischke. En particulier, dans l’introduction, ils posent la question simple suivant : vaut-il mieux, suite à une opération, renvoyer les gens chez eux, ou bien les laisser en chambre de convalescence, à l’hôpital ? Le soucis est que les données sont biaisées : oui, si on regarde les données brutes, ceux qui sont rentré chez eux se sont remis plus vite… mais peut être a-t-on autorisé les gens à rentrer justement parce qu’ils allaient bien ? Et on retrouve ce biais de sélection un peu partout. Probablement même encore davantage dans les problèmes relatifs à la justice.

Prenons un exemple simple. On va s’intéresser à la probabilité d’être coupable d’un crime, et d’être condamné par la justice. C’est ma variable Y. Pour avoir un modèle économétrique, il me faut des variables explicatives. J’ai une variable X_3 et une autre X_4, la première influençant vraiment la probabilité de commettre un crime, et l’autre étant juste du bruit. Maintenant, entrons dans le coeur du modèle. Supposons que j’ai une variable corrélée à ma richesse \tilde{X}_2, comme le lieu d’habitation, un critère racial, ou pour faire simple le diplôme. Notons X_1 le niveau de diplôme de la personne (mais encore une fois, on peut imaginer n’importe quoi, l’idée est que cette variable est corrélée à la richesse). Et surtout, cette variable n’est pas du tout une variable causale dans mon modèle. La probabilité de commettre un crime n’est pas influencée par X_1. Une personne ayant un diplôme de second cycle universitaire a autant de chance de commettre un crime qu’une personne ayant raté le brevet. On va même supposer que la probabilité de commettre un crime soit indépendante de la richesse. En fait, la richesse me permet de payer un bon avocat. Et avoir un bon avocat permet de diminuer la probabilité d’être condamné. Pour faire simple, j’ai deux niveaux de richesse, qui me donne soit un avocat commis d’office, soit un avocat qui s’impliquera peut être davantage (je simplifie à des fins pédagogiques). Bref, mon modèle aura la forme suivante : la probabilité d’être condamné pour un crime s’écrit p=(0.025+0.05\cdot X_3)\cdot[k\cdot \mathbb{1}(X_2=R)+\mathbb{1}(X_2=P)] Autrement dit, ma probabilité est une fonction croissante (linéaire) de X_3, et si je suis riche (ce qu’on note X_2=R), j’ai moins de chance d’être condamné, avec un facteur k\in[0,1], toutes choses étant égales par ailleurs, ceteris paribus, histoire de parler latin comme les juristes.

Pour simuler un jeu de données, on considère une population d’un million de personnes, et je suppose que ma richesse est très corrélée avec mon niveau de diplôme. Pour info, la classe A ce sont les gens les plus diplômés, et je suppose qu’un bon avocat divise par trois la chance d’être condamné.

 X1=cut(vectX[,1],breaks=c(-10,1,-1,10), c("C","B","A"))
 Y1=rbinom(n,size = 1,prob=(.025+.05*X3)*(k*(X2=="R")+(X2=="P")))

On peut faire un peu de statistiques pour regarder les caractéristiques des gens condamnés (uniquement)


Par exemple si je regarde la proportion de gens les plus diplômés: ils sont 15% dans la population totale, et 5% parmi les personnes condamnées

[1] 0.158734
[1] 0.05016072
data:  B0$X1 == "A" and B1$X1 == "A"
t = 219.5, df = 1283500, p-value < 2.2e-16
95 percent confidence interval:
 0.1076038 0.1095428
sample estimates:
 mean of x  mean of y 
0.15873400 0.05016072

Autrement dit, on a une différence significative entre ces deux pourcentages. Pareil si je regarde la valeur moyenne de X_3, qui était de 50% dans la population et de 58% dans la population condamnée

data:  B0$X3 and B1$X3
t = -163.46, df = 844940, p-value < 2.2e-16
95 percent confidence interval:
 -0.08455430 -0.08255058
sample estimates:
mean of x mean of y 
0.5000094 0.5835618

alors que pour X_4, on a 50% dans la population totale, et 50% dans la population condamnée. Mais on s’y attendait….

data:  B0$X4 and B1$X4
t = -0.81161, df = 48645, p-value = 0.417
95 percent confidence interval:
 -0.003869879  0.001603452
sample estimates:
mean of x mean of y 
0.5001326 0.5012659

Passons à la seconde phase. Il semble y avoir un vrai intérêt dans la littérature sur les modèles actuariels pour la construction de modèles liés à la probabilité de récidive. Là encore, je vais faire une hypothèse forte : faire de la prison n’influence en rien la probabilité de récidiver. Autrement dit, pour cette population déjà condamnée, je vais refaire tourner mes simulations, avec une probabilité de commettre un crime qui dépend seulement de X_3, mais avec toujours notre effet richesse qui permet de payer un avocat qui fera baisser la probabilité d’être condamné.

 Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P")))

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

 Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P")))

Je fais ici un score logistique car c’est le plus logique.

            Estimate Std. Error z value  Pr(|z|)    
Inter.) -4.04753    0.11592 -34.916    2e-16 ***
X1C      0.61311    0.10914   5.618 1.93e-08 ***
X1B      0.68807    0.09977   6.896 5.34e-12 ***
X3       1.00610    0.08052  12.495    2e-16 ***
X4      -0.13562    0.07286  -1.861   0.0627 .

Premier élément à observer, la variable X_4 est à la limite d’être significative… Certains auteurs sur ces thèmes diraient que c’est un signal faible. Pourtant, en l’occurence, c’est du bruit ! Juste du bruit !! Et parmi les signaux forts, très clairement significatifs, on retrouve X_1, notre niveau de diplôme. Ici, si on est faiblement qualifié (niveaux B ou C) alors on a une chance de récidive beaucoup plus importante ! C’est très fortement significatif ! Alors que pour rappel, la probabilé de commettre un crime, et de récidiver n’est aucunement causé par cette variable ! C’est juste une variable corrélée avec la variable richesse, qui permet de faire baisser la probabilité d’être condamné.

Plus intéressant, si je regarde la valeur du coefficient, en fonction de k, autrement dit la valeur d’un bon avocat (que l’on peut s’offrir si on est riche) on obtient la courbe suivante.

Autrement, sur la partie droite, même si un bon avocat ne baisse que de 5% la chance d’être condamné, le niveau de diplôme semble être significativement liée à la probabilté de récidive. Avec un pic autour de 30%. Si les avocats ont un pouvoir incroyable, par exemple en divisant par 20 la probabilité d’être condamné, étrangement, la variable de diplôme redevient très peu significative, mais je pense que c’est parce qu’on n’a alors plus personne de très diplômé dans la base.

On pourrait me reprocher d’avoir une corrélation trop forte entre richesse et niveau de diplôme, ici. En mettant une corrélation à 70% (dans mon modèle Gaussien latent), j’obtiens

qui donne un effet moindre, mais là encore, si un bon avocat permet de réduire la probabilité d’être condamné par un facteur compris entre 5% et 90%, alors la variable de diplôme semblera significative.

Mon modèle n’est pas faux ici. C’est essentiellement que j’ai un biais dans mes données. Et ne pas en tenir compte fait que le modèle n’a aucun sens, et ne peut pas être interprété sans commettre d’erreurs (et ici je cumule les âneries, en voyant un signal faible alors qu’il y a du bruit, et en voyant un signal fort pour une variable qui n’a aucun lien causal). J’insiste lourdement, mais le diplôme n’intervient aucunement ici. Et on peut remplacer le diplôme par toute variable corrélée à la richesse, comme locataire et propriétaire, etc.

Bref, l’utilisation de modèles économétrique dans des problèmes judiciaires n’est pas simple ! Je parle ici de biais de sélection, mais la formulation qui me semble la plus juste serait une mauvaise spécification du modèle (on n’utilise pas les bonnes variables, mais ces dernières sont peut-être non-observables – ou difficiles à justifier éthiquement : peut-on demander le revenu de la personne pour faire un score de récidive, et avoir un score qui dépende effectivement du revenu ?). Il y a aussi un problème d’identification, car on a du mal à savoir ce qui influence le fait de commettre un crime, et la condamnation. Ce qu’on observe, ce sont les condamnations, pas les actes (je ne parle pas des condamnations à tort ici, que j’avais pu évoquer dans un autre billet, sur le principe de précaution). Mais j’essayerais de revenir sur ce point après les vacances….

Dernière anecdote pour finir. Après avoir découvert ces modèles via Virginie Gautron, j’avais discuté avec une amie, présidente de tribunal, qui me parlait de la nullité des experts qu’elle avait sous la main. Et elle voulait faire un peu de ménage. Bref, elle me demandait si je pouvais faire un modèle attribuant un score à un expert. Et là encore, ma réponse a été que oui, faire un modèle c’est facile, à condition d’avoir des données. Comment puis-je avoir un score Y d’un expert? Qui sait ce qu’est un bon expert? C’est quelqu’un contre qui on ne demande pas de contre-expertise? Malheureusement, la justice ne fonctionne pas ainsi, et la demande n’a rien à voir avec la qualité de l’expert. Et la qualité, vu du point de vue de la juge, est une notion qui peut être très différente de celle des deux parties qui s’opposent… Bref, je pourrais facilement faire un modèle. Mais en avoir un qui ait du sens, c’est une autre paire de manches !

The myth of interpretability of econometric models

There are important discussions nowadays about data modeling, to choose between the “two cultures” (as mentioned in Breiman (2001)), i.e. either econometrics models or machine/statistical learning models. We did discuss this issue recently in Econométrie et Machine Learning (so far only in French) with Emmanuel Flachaire and Antoine Ly. One argument often used by econometricians is the interpretability of econometric models. Or at least the attempt to get an interpretable model.

We also have this discussion in actuarial science, for instance in ratemaking (or insurance pricing). Machine learning based models usually perform better (for some a priori chosen metric), but actuaries claim that econometric models are more easily interpretable. In actuarial literature, we assume that claim frequency Y is driven by some non-observable risk factor \Theta, and therefore, we do have heterogeneous risks in our portfolio. And, it can be seen as legitimate to differentiate prices. Assume that this risk factor \Theta is strongly correlated with X_1, the age of the driver. Because in our portfolio, old drivers tend to have more accidents. Here, we could pretend to have a “causal story” (as defined in Freedman (2009)) because of a possible interpretation of the model. So it is natural here to consider a regression model of Y on X_1 to derive our actuarial pricing model. But assume that, possibly, risk factor \Theta is also strongly correlated with X_2, that can be related to spatial features (say latitude, which denoted a north/south position). Because in our portfolio, drivers living in the south tend to have more accidents (reads are known to be more dangerous there). Here, we could pretend to have a second “causal story”.

Of course, since \Theta is strongly correlated with X_1 and X_2, it means that X_1 and X_2 are strongly correlated. Here also, this correlation can be interpreted (not in a causal way as previously, but still), since we know that old people like to live in southern regions. So, what should we do here ? Let us run some simulations to  illustrate.


Our first idea was to consider a model where Y is “explained” by the first variable X_1,

         Estimate Std. Error z value Pr(>|z|)    
(Inter.) -2.97778    0.01544 -192.88   <2e-16 ***
X1        0.97926    0.01092   89.64   <2e-16 ***

As expected, our variable is “significant”, but also, probably more interesting, X_2, has no impact on the residuals

          Estimate Std. Error t value Pr(>|t|)
(Inter.) 0.0003618  0.0031696   0.114    0.909
X2       0.0028601  0.0031467   0.909    0.363

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

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

         Estimate Std. Error z value Pr(>|z|)    
(Inter.) -2.97724    0.01544 -192.81   <2e-16 ***
X2        0.97915    0.01093   89.56   <2e-16 ***

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

          Estimate Std. Error t value Pr(>|t|)
(Inter.) 0.0004863  0.0031733   0.153    0.878
X1       0.0027979  0.0031504   0.888    0.374

The story is similar here. If we correct from the spatial pattern, claims frequency does not depend on the age of the driver.

So, what should we do now? We do have two models, and each of them is as interpretable as the other one. Note that we can not use any statistical tool to distinguish the two: they are comparable

[1] 51013.39
[1] 51013.15

Why not incorporate the two explanatory variables X_1 and X_2, at the same time, in our regression model, and let “the model” decide what to do…?

         Estimate Std. Error  z value Pr(>|z|)    
(Inter.) -2.98132    0.01547 -192.723    2e-16 ***
X1        0.49310    0.06226    7.920 2.38e-15 ***
X2        0.49375    0.06225    7.931 2.17e-15 ***

It looks like we completely lost the interpretability of the model, since our two explanatory variables are (strongly) correlated. Actually, instead of saying “use one, and drop the other one (since it brings no further information)”, it says “use both, each one will explain half of the variable”. Strange interpretation, isn’t it?  So why not try some LASSO here?


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

A la recherche des homonymes…

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Les lois statistiques, du micro au macro

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

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

Norme, convention et aspects éthiques

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

I Got The Feelin’

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

Here is a function to extract information from the website,

> library(XML)
> extractbpm = function(VBP,P){
+ url=paste("",VBP,"&genre=&page=",P,sep="")
+ download.file(url,destfile = "page.html")
+ tables=readHTMLTable("page.html")
+ return(tables)}

For instance

> extractbpm(115,13)
Artist Title
1 Eros Ramazzotti y Claudio Guidetti Dimelo A Mi
2 Everclear Volvo Driving Soccer Mom
3 Evils Toy Dear God
4 Expose In Walked Love
5 Fabolous ft. 2 Chainz When I Feel Like It
6 Fabolous ft. 2 Chainz When I Feel Like It
7 Fabolous ft. 2 Chainz When I Feel Like It
8 Fanny Lu Fanfarron
9 Featurecast Ain't My Style
10 Fem 2 Fem Obsession
11 Fernando Villalona Mi Delito
12 Fever Ray Triangle Walks
13 Firstlove Freaky
14 Fito Blanko Pegadito Suavecito
15 Flechazo Del Norte Mariposa Traicionera
16 Fluke Switch/Twitch
17 Flyleaf Something Better
18 FM Static The Next Big Thing
19 Fonseca Eres Mi Sueno
20 Fonseca ft. Maffio & Nayer Eres Mi Sueno
21 Francesca Battistelli Have Yourself A Merry Little Christmas
22 Frankie Ballard Young & Crazy
23 Frankie J. More Than Words
24 Frank Sinatra The Hucklebuck
25 Franz Ferdinand The Dark Of The Matinée
Mix BPM Genre Label Year
1 — 115 — Sony 2009
2 — 115 — Capitol Records 2003
3 — 115 — — —
4 — 115 — Arista Records 1994
5 Explicit 115 Urban Def Jam/Island Def Jam 2013
6 — 115 Urban Def Jam/Island Def Jam 2013
7 Radio Edit 115 Urban Def Jam/Island Def Jam 2013
8 — 115 Latin Pop Universal Latino 2011
9 Psychemagik Dub 115 — Jalapeno 2012
10 — 115 — Critique Records 1993
11 — 115 — Mt&vi Records/caminante Records 2001
12 Rex The Dog Remix 115 — Little Idiot/Mute 2012
13 — 115 — Jwp Music 2000
14 — 115 Merengue Mambo Crown Loyalty 2012
15 — 115 — Hacienda 2010
16 Album Version 115 — One Little Indian Records 2004
17 — 115 Alternative A&M/Octone 2013
18 — 115 — Tooth & Nail Records 2007
19 — 115 Merengue Mambo 10 2012
20 Urban Version 115 — 10 2012
21 — 115 — Word/Fervent/Warner Bros 2009
22 — 115 Country Warner Bros 2015
23 Mynt Rocks Radio Edit 115 — Columbia 2005
24 — 115 Jazz Columbia 1950
25 — 115 New Wave — 2004

We have here one of the few old songs, a 1950 tune by Frank Sinatra. If we scrap the website, with a simple loop (where the bpm is from 40 to 200). Start with

> vbp=40
> p=1

and then, a loop based on

> while(vbp<=2017){
+ F=extractbmp(vbp,p)
+ if(length(F)==1){
+ BASE=rbind(BASE,F[[1]][,c("Artist","Title","BPM","Year")])
+ p=p+1}
+ if(length(F)==0){
+ vbp=vbp+1
+ p=1}}

Then we should clean the dataset


and we end up with almost 50,000 tunes.

col="light blue")

Over the past 20 years, it looks like speed of tunes has declined (let us forget tunes of 2017, clearly we have a problem here…)


which is confirmed with a (smoothed) regression


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