Richesse et espérance de vie

Ce matin, je découvrais un graphique de l’INSEE qui présentait les taux de mortalité par sexe, âge et niveau de vie, avec entre autres le graphique suivant

Comme souvent avec l’INSEE, on a accès aux données… pas celles au niveau individuel (malheureusement) mais au moins on peut retravailler la visualisation. En fait, les données sont même encore plus fines, puisque les niveaux de richesses sont définis avec des tranches de 5%, et en plus, on a le détail entre hommes et femmes

b = read.csv2("MORT-RICHESSE.csv")
plot(b[,1],b[,2]/1000,col="red",type="l",ylab="% de survivants",xlab="Age")
lines(b[,1],b[,3]/1000,col="red",type="l",lty=2)          
lines(b[,1],b[,4]/1000,col="blue",type="l",lty=1)           
lines(b[,1],b[,5]/1000,col="blue",type="l",lty=2)  
legend("bottomleft",c("Femmes 95-100%","Hommes 95-100%","Femmes 0-5%","Hommes, 0-5%"),
bty="n",col=c("red","blue","red","blue"),lty=c(2,2,1,1))

Je me demandais si on ne pouvait pas tenter une lecture inverse de ce graphique : sur ce graphique, assez naturellement, on regarde à pourcentage donné l’écart entre la courbe en trait plein (les pauvres) et celle en trait pointillé (les riches). Si c’est cette information qu’on veut avoir, on peut alors tenter de la visualiser. Pour ça, il faut inverser notre fonction de survie (je l’ai fait rapidement, avec une simple interpolation linéaire… je pense qu’on peut faire mieux)

inversef = function(p,k=2){
  y=1-b[,k]/100000
  idx=sum(y<=p) 
  y1=y[idx-1]
  y2=y[idx]
  w1=(y1-p)/(y1-y2)
  w2=(p-y2)/(y1-y2)
  w2*b[idx-1,1]+w1*b[idx,1]
}

Ensuite, on peut construire les inverses, et mieux, les différences entre les courbes des riches et des pauvres

diffF = function(p) inversef(p,3)-inversef(p,2)
diffH = function(p) inversef(p,5)-inversef(p,4)
u  = seq(.01,.99,by=.01)
vF = Vectorize(diffF)(u)
vH = Vectorize(diffH)(u)
plot(u*100,vF,col="red",type="l",xlab="Probabilité (%)",ylab="Nombre d'années",ylim=c(0,max(vF,vH)))
lines(u*100,vH,col="blue",type="l",lty=1)           
legend("topright",c("Femmes","Hommes"),
bty="n",col=c("red","blue"),lty=c(1,1))

Je ne suis pas très à l’aise avec le graphique. Tout d’abord parce que je ne sais pas ce que la richesse indique (un pauvre de 20 ans peut devenir un riche de 50 ans, non ?), la richesse était souvent liée à l’âge. Après l’axe des abscisses me semble aussi avoir une interprétation compliquée : quand on regarde à 10%, on regarde les pauvres et les riches qui sont morts relativement jeunes (relativement car je regarde le quantile de la fonction de survie à tranche de richesse donnée). Autrement dit, si je me place à 10%, je compare les jeunes riches morts très jeunes (10% des riches seulement sont morts plus jeunes – c’est l’interprétation d’un quantile) et les jeunes pauvres (mort à un âge que seulement 90% des pauvres ont dépassé), on observe une différence de 20 ans environ. J’ai aussi l’impression qu’on pourrait dire que pour la majorité des hommes, les riches vivent 12 ans de plus que les pauvres, soit le double des femmes (de l’ordre de 6 ans).

On pourrait bien sûr se contenter de calculer les différences entre les aires, ce qui donne une différence entre les espérances de vie à la naissance des pauvres et des riches (comme le fait l’INSEE)

et qu’on peut visualiser sur le graphique suivant

plot(b[,1],b[,2]/1000,col="white",type="l",ylab="% de survivants",xlab="Age")
polygon(c(b[,1],rev(b[,1])),c(b[,3]/1000,rev(b[,2]/1000)),col="red",border=NA)

et un calcul donne une différence de l’ordre de 8 ans

sum(b[,3]-b[,2])/100000
[1] 8.239346

mais la visualisation raconte bien plus qu’un simple calcul d’aire. Par exemple, le graphique ci-dessous donne exactement le même écart entre les espérances de vie des pauvres et des riches

diff = sum(b[,3]-b[,2])/1000
y1 = b[,2]/1000
for(i in 1:100){
  y1[i] = b[i,2]/1000+min(100-b[i,2]/1000,diff)
  diff  = diff-(y1[i]-b[i,2]/1000)
}
plot(b[,1],b[,2]/1000,col="white",type="l",ylab="% de survivants",xlab="Age")
polygon(c(b[,1],rev(b[,1])),c(y1,rev(b[,2]/1000)),col="red",border=NA) 
sum(b[,3]-b[,2])/100000
lines(b[,1],b[,2]/1000,col="red")

Ici, on dit que la moitié des femmes riches meurent vers 81 ans, et les autres meurent au même âge qu’une femmes pauvre (mais une femme pauvre qui vivrait longtemps. Les distributions sont vraiment différentes, et c’est ça que je cherche à visualiser… Parce que la densité de l’âge au décès ne me semble pas forcément très simple à analyser…

Comme toujours, les commentaires sont ouverts si certains ont des idées quant à la visualisation de ces données…

INF7100 – Initiation à la science des données et à l’intelligence artificielle

Cette été, le cours INF7100 Initiation à la science des données et à l’intelligence artificielle – devrait être offert pour la première fois. Le premier cours aura lien le mardi 28 avril, et on le donnera conjointement avec Marie-Jean Meurs (du département d’informatique) et Jean-Hugues Roy (de l’école des médias). J’essayerai de mettre en ligne des informations au fur et à mesure…

 

Testing for a causal effect (with 2 time series)

A few days ago, I came back on a sentence I found (in a French newspaper), where someone was claiming that

“… an old variable explains 85% of the change in a new variable. So we can talk about causality”

and I tried to explain that it was just stupid : if we consider the regression of the temperature on day t+1 against the number of cyclist on day t, the R^2 exceeds 80%… but it is hard to claim that the number of cyclists on specific day will actually cause the temperature on the next day…

Nevertheless, that was frustrating, and I was wondering if there was a clever way to test for causality in that case. A popular one is Granger causality (I can mention a paper we published a few years ago where we use such a test, Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media). To explain that test, consider a bivariate time series (just like the one we have here), \boldsymbol{z}_t=(x_t,y_t), and consider some bivariate autoregressive model
{\displaystyle {\begin{bmatrix}x_{t}\\y_{t}\end{bmatrix}}={\begin{bmatrix}c_{1}\\c_{2}\end{bmatrix}}+{\begin{bmatrix}a_{1,1}&\textcolor{red}{a_{1,2}}\\\textcolor{blue}{a_{2,1}}&a_{2,2}\end{bmatrix}}{\begin{bmatrix}x_{t-1}\\y_{t-1}\end{bmatrix}}+{\begin{bmatrix}u_{t}\\v_{t}\end{bmatrix}}}where \boldsymbol{\varepsilon}_t=(u_t,v_t) is some bivariate white noise, in the sense that (i) {\displaystyle \mathbb{E} (\boldsymbol{\varepsilon}_{t})=\boldsymbol{0}} (the noise is centered) (ii) {\displaystyle \mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t}^\top)=\Omega } , so the variance matrix is constant, but possibly non-diagonal (iii) {\displaystyle \mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t-h}^\top)=\boldsymbol{0} } for all h\neq 0. Note that we can use the simplified expression{\displaystyle {\boldsymbol{z}_t=\boldsymbol{c}+\boldsymbol{A}\boldsymbol{z}_{t-1}+\boldsymbol{\varepsilon}_t}}Now, Granger test is based on several quantities. With off-diagonal terms of matrix \Omega, we have a so-called instantaneous causality, and since \Omega is symmetry, we will write x\leftrightarrow y. With off-diagonal terms of matrix \boldsymbol{A}, we have a so-called lagged causality, with either \textcolor{blue}{x\rightarrow y} or \textcolor{red}{x\leftarrow y} (and possibly both, if both terms are significant).

So I wanted to try on my two-variable problem.

df = read.csv("cyclistsTempHKI.csv")
dfts = cbind(C=ts(df$cyclists,start = c(2014, 1,2), frequency = 365),
             T=ts(df$meanTemp,start = c(2014, 1,2), frequency = 365))
library(vars)

I now have “time series” objects, and we can fit a VAR model,

var2 = VAR(dfts, p = 1, type = "const")
coefficients(var2)
$C
         Estimate   Std. Error   t value      Pr(>|t|)
C.l1    0.8684009   0.02889424 30.054460 8.080226e-107
T.l1   70.3042012  20.07247411  3.502518  5.102094e-04
const 807.6394001 187.75472482  4.301566  2.110412e-05
 
$T
           Estimate   Std. Error   t value     Pr(>|t|)
C.l1   0.0003865391 6.257596e-05  6.177118 1.540467e-09
T.l1   0.6611135594 4.347074e-02 15.208241 6.086394e-42
const -1.6413074565 4.066184e-01 -4.036481 6.446018e-05

For instant, we can run a causality, to test if the number of cyclists can cause the temperature (on the next day)

causality(var2, cause = "C")
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 38.157, df1 = 1, df2 = 842, p-value = 1.015e-09

Here, we should clearly reject H_0, which is that there is no causal effect. Which is the way statisticians say that there should be some causal effect between the number of cyclist and the temperature…

So clearly, something is wrong here. Either it is some sort of superpower that cyclists are not aware of. Or this test that was used for forty years (Clive Granger even got a Nobel price for it) is not working. Or we missed something. Actually… I think we missed something here. Possibly because the series are not stationary. We can almost see it with

Phi = matrix(c(coefficients(var2)$C[1:2,1],coefficients(var2)$T[1:2,1]),2,2)
eigen(Phi)
eigen() decomposition
$values
[1] 0.9594810 0.5700335

where the highest eigenvalue is very close to one. But actually, we look here at the temperature…

plot(dfts)

so, at least, we should expect some seasonal unit root here. So let us use two techniques. The first one is a classical one-year difference, \Delta_{365}\boldsymbol{z}_t=\boldsymbol{z}_t-\boldsymbol{z}_{t-365}

var2 = VAR(diff(dfts,365), p = 1, type = "const")
coefficients(var2)
$C
          Estimate   Std. Error   t value     Pr(>|t|)
C.l1     0.8376424   0.07259969 11.537823 1.993355e-16
T.l1    42.2638410  28.58783276  1.478386 1.449076e-01
const -507.5514795 219.40240747 -2.313336 2.440042e-02
 
$T
         Estimate   Std. Error   t value     Pr(>|t|)
C.l1  0.000518209 0.0003277295 1.5812096 1.194623e-01
T.l1  0.598425288 0.1290511945 4.6371154 2.162476e-05
const 0.547828079 0.9904263469 0.5531235 5.823804e-01

The test on the fited VAR model yields

causality(var2, cause = "C") 
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 2.5002, df1 = 1, df2 = 112, p-value = 0.1167

i.e., with a 11% p-value, we should reject the assumption that the number of cyclists cause the temperature (on the next day), and actually, we should also reject the other way

causality(var2, cause = "T") 
$Granger
 
	Granger causality H0: T do not Granger-cause C
 
data:  VAR object var2
F-Test = 2.1856, df1 = 1, df2 = 112, p-value = 0.1421

Nevertheless, if we look at the instantaneous causality, this one makes more sense

$Instant
 
	H0: No instantaneous causality between: T and C
 
data:  VAR object var2
Chi-squared = 13.081, df = 1, p-value = 0.0002982

The second idea would be to use a one day difference, \Delta_{1}\boldsymbol{z}_t=\boldsymbol{z}_t-\boldsymbol{z}_{t-1} and to fit a VAR model on that one

VARselect(diff(dfts,1), lag.max = 4, type="const")
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      3      2      3

but on that one, a VAR(1) model – with only one lag – might not be sufficient. It might be better to consider a VAR(3)

var2 = VAR(diff(dfts,1), p = 3, type = "const")

and on that one, one more time, we should reject the causal effect of the number of cyclists on the temperature (on the next day)

causality(var2, cause = "C")  
$Granger
 
	Granger causality H0: C do not Granger-cause T
 
data:  VAR object var2
F-Test = 0.67644, df1 = 3, df2 = 828, p-value = 0.5666

and this time, there could be a (lagged) causal effect of the temperature on the number of cyclists

causality(var2, cause = "T")  
$Granger
 
	Granger causality H0: T do not Granger-cause C
 
data:  VAR object var2
F-Test = 7.7981, df1 = 3, df2 = 828, p-value = 3.879e-05
 
$Instant
 
	H0: No instantaneous causality between: T and C
 
data:  VAR object var2
Chi-squared = 55.83, df = 1, p-value = 7.905e-14

but nothing instantaneously… So it looks like Granger causality performs well on that one !

Lasso Regression (home made)

Again, this post is related to my MAT7381 course, where we will see that it is actually possible to write our own code to compute Lasso regression, \min\left\lbrace\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|_{\ell_2}^2+\lambda\|\mathbf{\beta}\|_{\ell_1}\right\rbraceWe have to define the soft-thresholding functionS(z,\gamma)=\text{sign}(z)\cdot(|z|-\gamma)_+=\begin{cases}z-\gamma&\text{ if }\gamma>|z|\text{ and }z<0\\z+\gamma&\text{ if }\gamma<|z|\text{ and }z<0 \\0&\text{ if }\gamma\geq|z|\end{cases}The R function would be

soft_thresholding = function(x,a){
sign(x) * pmax(abs(x)-a,0)
}

To solve our optimization problem, set\mathbf{r}_j=\mathbf{y} - \left(\beta_0\mathbf{1}+\sum_{k\neq j}\beta_k\mathbf{x}_k\right)=\mathbf{y}-\widehat{\mathbf{y}}^{(j)}
so that the optimization problem can be written, equivalently
\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p [\mathbf{r}_j-\beta_j\mathbf{x}_j]^2+\lambda |\beta_j|\right\rbrace
hence\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p \beta_j^2\|\mathbf{x}_j\|-2\beta_j\mathbf{r}_j^T\mathbf{x}_j+\lambda |\beta_j|\right\rbrace
and one gets
\beta_{j,\lambda} = \frac{1}{\|\mathbf{x}_j\|^2}S(\mathbf{r}_j^T\mathbf{x}_j,n\lambda)
or, if we develop
\beta_{j,\lambda} = \frac{1}{\sum_i x_{ij}^2}S\left(\sum_ix_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)
Again, if there are weights \mathbf{\omega}=(\omega_i), the coordinate-wise update becomes
\beta_{j,\lambda,{\color{red}{\omega}}} = \frac{1}{\sum_i {\color{red}{\omega_i}}x_{ij}^2}S\left(\sum_i{\color{red}{\omega_i}}x_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)
The code to compute this componentwise descent is

lasso_coord_desc = function(X,y,beta,lambda,tol=1e-6,maxiter=1000){
  beta = as.matrix(beta)
  X = as.matrix(X)
  omega = rep(1/length(y),length(y))
  obj = numeric(length=(maxiter+1))
  betalist = list(length(maxiter+1))
  betalist[[1]] = beta
  beta0list = numeric(length(maxiter+1))
  beta0 = sum(y-X%*%beta)/(length(y))
  beta0list[1] = beta0
  for (j in 1:maxiter){
    for (k in 1:length(beta)){
      r = y - X[,-k]%*%beta[-k] - beta0*rep(1,length(y))
      beta[k] = (1/sum(omega*X[,k]^2))*
        soft_thresholding(t(omega*r)%*%X[,k],length(y)*lambda)
    }
    beta0 = sum(y-X%*%beta)/(length(y))
    beta0list[j+1] = beta0
    betalist[[j+1]] = beta
    obj[j] = (1/2)*(1/length(y))*norm(omega*(y - X%*%beta - 
           beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta))
    if (norm(rbind(beta0list[j],betalist[[j]]) - 
             rbind(beta0,beta),'F') &lt; tol) { break } 
  } 
  return(list(obj=obj[1:j],beta=beta,intercept=beta0)) }

For instance, consider the following (simple) dataset, with three covariates

chicago = read.table("http://freakonometrics.free.fr/chicago.txt",header=TRUE,sep=";")

that we can “normalize” (or “standardize“)

X = model.matrix(lm(Fire~.,data=chicago))[,2:4]
for(j in 1:3) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
y = chicago$Fire
y = (y-mean(y))/sd(y)

To initialize the algorithm, use the OLS estimate

beta_init = lm(Fire~0+.,data=chicago)$coef

For instance

lasso_coord_desc(X,y,beta_init,lambda=.001)
$obj
[1] 0.001014426 0.001008009 0.001009558 0.001011094 0.001011119 0.001011119
 
$beta
          [,1]
X_1  0.0000000
X_2  0.3836087
X_3 -0.5026137
 
$intercept
[1] 2.060999e-16

and we can get the standard Lasso plot by looping,

Quantile Regression (home made, part 2)

A few months ago, I posted a note with some home made codes for quantile regression… there was something odd on the output, but it was because there was a (small) mathematical problem in my equation. So since I should teach those tomorrow, let me fix them.

Median

Consider a sample \{y_1,\cdots,y_n\}. To compute the median, solve\min_\mu \left\lbrace\sum_{i=1}^n|y_i-\mu|\right\rbracewhich can be solved using linear programming techniques. More precisely, this problem is equivalent to\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^na_i+b_i\right\rbracewith a_i,b_i\geq 0 and y_i-\mu=a_i-b_i, \forall i=1,\cdots,n. Heuristically, the idea is to write y_i=\mu+\varepsilon_i, and then define a_i‘s and b_i‘s so that \varepsilon_i=a_i-b_i and |\varepsilon_i|=a_i+b_i, i.e. a_i=(\varepsilon_i)_+=\max\lbrace0,\varepsilon_i\rbrace=|\varepsilon|\cdot\boldsymbol{1}_{\varepsilon_i>0}andb_i=(-\varepsilon_i)_+=\max\lbrace0,-\varepsilon_i\rbrace=|\varepsilon|\cdot\boldsymbol{1}_{\varepsilon_i<0}denote respectively the positive and the negative parts.

Unfortunately (that was the error in my previous post), the expression of linear programs is\min_{\mathbf{z}}\left\lbrace\boldsymbol{c}^\top\mathbf{z}\right\rbrace\text{ s.t. }\boldsymbol{A}\mathbf{z}=\boldsymbol{b},\mathbf{z}\geq\boldsymbol{0}In the equation above, with the a_i‘s and b_i‘s, we’re not far away. Except that we have \mu\in\mathbb{R}, while it should be positive. So similarly, set \mu=\mu^+-\mu^- where \mu^+=(\mu)_+ and \mu^-=(-\mu)_+.

Thus, let\mathbf{z}=\big(\mu^+;\mu^-;\boldsymbol{a},\boldsymbol{b}\big)^\top\in\mathbb{R}_+^{2n+2}and then write the constraint as \boldsymbol{A}\mathbf{z}=\boldsymbol{b} with \boldsymbol{b}=\boldsymbol{y} and \boldsymbol{A}=\big[\boldsymbol{1}_n;-\boldsymbol{1}_n;\mathbb{I}_n;-\mathbb{I}_n\big]And for the objective function\boldsymbol{c}=\big(\boldsymbol{0},\boldsymbol{1}_n,-\boldsymbol{1}_n\big)^\top\in\mathbb{R}_+^{2n+2}

To illustrate, consider a sample from a lognormal distribution,

n = 101 
set.seed(1)
y = rlnorm(n)
median(y)
[1] 1.077415

For the optimization problem, use the matrix form, with 3n constraints, and 2n+1 parameters,

library(lpSolve) 
X = rep(1,n) 
A = cbind(X, -X, diag(n), -diag(n))
b = y
c = c(rep(0,2), rep(1,n),rep(1,n))
equal_type = rep("=", n) 
r = lp("min", c,A,equal_type,b)
head(r$solution,1)
[1] 1.077415

It looks like it’s working well…

Quantile

Of course, we can adapt our previous code for quantiles

tau = .3
quantile(y,tau)
      30% 
0.6741586

The linear program is now\min_{q^+,q^-,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbracewith a_i,b_i,q^+,q^-\geq 0 and y_i=q^+-q^-+a_i-b_i, \forall i=1,\cdots,n. The R code is now

c = c(rep(0,2), tau*rep(1,n),(1-tau)*rep(1,n))
r = lp("min", c,A,equal_type,b)
head(r$solution,1)
[1] 0.6741586

So far so good…

Quantile Regression

Consider the following dataset, with rents of flat, in a major German city, as function of the surface, the year of construction, etc.

base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)

The linear program for the quantile regression is now\min_{\boldsymbol{\beta}^+,\boldsymbol{\beta}^-,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbracewith a_i,b_i\geq 0 and y_i=\boldsymbol{x}^\top[\boldsymbol{\beta}^+-\boldsymbol{\beta}^-]+a_i-b_i\forall i=1,\cdots,n and \beta_j^+,\beta_j^-\geq 0 \forall j=0,\cdots,k. So use here

require(lpSolve) 
tau = .3
n=nrow(base)
X = cbind( 1, base$area)
y = base$rent_euro
K = ncol(X)
N = nrow(X)
A = cbind(X,-X,diag(N),-diag(N))
c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N))
b = base$rent_euro
const_type = rep("=",N)
r = lp("min",c,A,const_type,b)
beta = r$sol[1:K] -  r$sol[(1:K+K)]
beta
[1] 148.946864   3.289674

Of course, we can use R function to fit that model

library(quantreg)
rq(rent_euro~area, tau=tau, data=base)
Coefficients:
(Intercept)        area 
 148.946864    3.289674

Here again, it seems to work quite well. We can use a different probability level, of course, and get a plot

plot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")),
     ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5)
sf=0:250
yr=r$solution[2*n+1]+r$solution[2*n+2]*sf
lines(sf,yr,lwd=2,col="blue")
tau = .9
r = lp("min",c,A,const_type,b)
tail(r$solution,2)
[1] 121.815505   7.865536
yr=r$solution[2*n+1]+r$solution[2*n+2]*sf
lines(sf,yr,lwd=2,col="blue")

And we can adapt the later to multiple regressions, of course,

X = cbind(1,base$area,base$yearc)
K = ncol(X)
N = nrow(X)
A = cbind(X,-X,diag(N),-diag(N))
c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N))
b = base$rent_euro
const_type = rep("=",N)
r = lp("min",c,A,const_type,b)
beta = r$sol[1:K] -  r$sol[(1:K+K)]
beta
[1] -5542.503252     3.978135     2.887234

to be compared with

library(quantreg)
rq(rent_euro~ area + yearc, tau=tau, data=base)
 
Coefficients:
 (Intercept)         area        yearc 
-5542.503252     3.978135     2.887234 
 
Degrees of freedom: 4571 total; 4568 residual

De la pratique de la régression

Depuis le début de la session, j’ai imposé une petite innovation, en donnant, environ une semaine sur deux, un petit exercice (obligatoire mais non noté) avant le cours, en vue de forcer à réfléchir (et de donner des éléments de réponse). Par exemple pour le premier cours, il fallait “prévoir” une valeur manquante, et le but était de montrer que, naturellement, on choisit la valeur moyenne.

Pour demain, j’avais posé un exercice un peu plus compliqué, sachant qu’on avait vu, lors du dernier cours, comme faire une régression linéaire, et qu’on avait fini en discutant les tests simples (en lien avec la significativité) et les tests multiples. Pour l’exercice, j’avais mis en ligne une petite base de données,

download.file("http://freakonometrics.free.fr/data3.RData","data3.RData")
load("data3.RData")
str(df)
'data.frame':	147 obs. of  10 variables:
 $ Y : num  11.72 15.91 14.19 11.15 8.31 ...
 $ X1: num  1.33 3.18 0.28 2.08 0.11 1.67 1.97 1.27 4.38 0.52 ...
 $ X2: num  3.66 3.75 3.32 2.68 4.97 2.98 4.56 1.78 2.83 6.36 ...
 $ X3: num  1.41 3.01 0.34 2.19 0.25 1.69 2.01 1.25 4.41 0.43 ...
 $ X4: num  -3.53 -4.46 -3.35 -7.54 -7.02 -2.53 -6.1 -5.99 -3.92 -5.84 ...
 $ X5: num  0.57 0.01 -0.7 1.62 -0.95 -1.37 1.18 -0.72 2.63 -1.63 ...
 $ X6: num  -0.82 1.2 3.03 -0.91 -1.6 1.77 1 -1.33 1.31 -0.7 ...
 $ X7: num  1.01 0.06 2.02 3.63 2.66 2.53 1.29 3.5 1.17 1.8 ...
 $ X8: num  8.31 8.52 9.78 7.34 7.26 ...
 $ X9: num  6.04 6.53 7.52 5.61 4.52 6.06 6.2 5.99 6.93 4.38 ...

Je vais mettre ici les questions que je posais, et donner des pistes de réflexions, non pas sur les réponses attendues (je n’attends rien de cet exercice à part une réflexion), mais sur la discussion que peut amener chacune des questions,

  • Faîtes un modèle linéaire pour expliquer Y en fonction des neuf variables explicatives. Combien de variables explicatives garderiez-vous ?

Commençons par faire une régression sur toutes les variables

summary(lm(Y~., data=df))
 
Coefficients:
             Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  8.046133   2.193448   3.668 0.000349 ***
X1           0.342293   0.915894   0.374 0.709186    
X2          -0.040479   0.103409  -0.391 0.696073    
X3           1.683875   0.897278   1.877 0.062693 .  
X4          -0.009254   0.062382  -0.148 0.882295    
X5          -1.085367   0.113840  -9.534  &lt; 2e-16 ***
X6           0.983207   0.111830   8.792 5.49e-15 ***
X7          -0.015646   0.087483  -0.179 0.858327    
X8           0.012165   0.094756   0.128 0.898033    
X9           0.172210   0.180605   0.954 0.342005    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.019 on 137 degrees of freedom
Multiple R-squared:  0.9228,	Adjusted R-squared:  0.9178 
F-statistic:   182 on 9 and 137 DF,  p-value: &lt; 2.2e-16

On a, dans la sortie, une dizaine de tests de Student qui sont évoqués, correspondant au test de l’hypothèse H_0:\beta_j=0 (contre l’hypothèse alternative (bilatérale) H_1:\beta_j\neq 0), dans un modèle de la forme y_i=\beta_0+\beta_1x_{1,i}+\dots+\beta_9x_{9,i}+\varepsilon_iC’est ce qu’on appelle le test de significativité de la variable x_j (oui, comme on l’a vu en cours, on peut dire le test parce que les autres tests classiques – Fisher, ou Wald – sont équivalent – sauf qu’au lieu de regarder t on regarde la statistique t^2 – ce qui présente l’avantage de voir la statistique de test comme une forme de distance à l’hypothèse H_0 : si c’est trop grand, on rejette…). Avec un seuil d’acceptation de l’ordre de 5%, on nous dit que 6 variables ne sont pas significatives. Mais gardons bien en mémoire que le test de significativité de x_j est fait ici en supposant que toutes les autres variables restent dans le modèle. Autre chose: avec un seuil d’acceptation de l’ordre de 10%, une des variables (la troisième) peut être vue comme significative.

Faisons un test multiple, pour savoir si on peut supprimer 6 des 9 variables explicatives possibles (faisons le à la main, inutile d’aller chercher un package pour le faire)

reg1 = lm(formula = Y ~ ., data = df)
reg0 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg0)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.11031    0.16788   54.27   &lt;2e-16 ***
X3           1.96784    0.07604   25.88   &lt;2e-16 ***
X5          -0.99296    0.08751  -11.35   &lt;2e-16 ***
X6           1.08761    0.05934   18.33   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: &lt; 2.2e-16 anova(reg0,reg1) Analysis of Variance Table Model 1: Y ~ X3 + X5 + X6 Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 Res.Df RSS Df Sum of Sq F Pr(&gt;F)
1    143 146.16                           
2    137 142.12  6    4.0338 0.6481 0.6916

Le test de Fisher nous dit qu’on peut accepter l’hypothèse que les 6 coefficients sont nuls – ici H_0:\beta_1=\beta_2=\beta_4=\beta_7=\beta_8=\beta_9=0(qui est un test multiple, contrairement au test de Student précédant). Mais il nous dit aussi, qu’individuellement, les trois variables restants semblent significatives. Donc j’aurais tendance à garder 3 variables explicatives (je ne parle pas de la constante : on garde toujours la constante – qui n’explique pas grand chose, sauf la valeur moyenne de y).

  • faites une prévision pour un individu dont on sait que X3=1, X5=1 et X6=8

(en réalité, la vraie question que j’ai posée contenait une typo ce qui la rendait vicieuse parce que ce n’est pas le modèle qu’on vient de construire… mais pour commencer, regardons cette question)

On vient de calibrer ce modèle, donc il suffit de faire une prévision,

predict(reg0, newdata=data.frame(X3=1, X5=1, X6=8))
       1 
18.78603

Mais ce n’était pas la vraie question…

  • faites une prévision pour un individu dont on sait que X1=3, X5=1 et X6=8

Je pense que pour répondre à cette question, il convient d’oublier tout ce qu’on vient de voir. On nous donne trois informations, et on va voir si on peut les exploiter. Autrement dit, on va commencer par regarder la régression sur ces trois variables,

reg156 = lm(formula = Y ~ X1+X5+X6, data = df)
summary(reg156)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.06688    0.17222   52.65   &lt;2e-16 ***
X1           1.98342    0.07800   25.43   &lt;2e-16 ***
X5          -1.01438    0.08947  -11.34   &lt;2e-16 ***
X6           1.07824    0.06039   17.86   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.026 on 143 degrees of freedom
Multiple R-squared:  0.9183,	Adjusted R-squared:  0.9166 
F-statistic: 535.8 on 3 and 143 DF,  p-value: &lt; 2.2e-16

Le modèle est ici bon, les trois variables étant significatives. Cela dit, je dis qu’il est “bon” mais on verra vendredi comment discuter davantage ce point… En tous cas, on peut tenter de faire une prévision, et on obtient

predict(reg156, newdata=data.frame(X1=3, X5=1, X6=8))
       1 
22.62867
  • faites une prévision pour un individu dont on sait que X1=3, X3=2, X5=1 et X6=8

Comme auparavant, on nous donne 4 informations… on va regarder le modèle avec les 4 variables

reg1356 = lm(formula = Y ~ X1+X3+X5+X6, data = df)
summary(reg1356)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.10458    0.17131  53.148   &lt;2e-16 ***
X1           0.16362    0.89028   0.184    0.854    
X3           1.80663    0.88052   2.052    0.042 *  
X5          -0.99558    0.08895 -11.192   &lt;2e-16 ***
X6           1.08649    0.05986  18.152   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.014 on 142 degrees of freedom
Multiple R-squared:  0.9207,	Adjusted R-squared:  0.9184 
F-statistic: 411.9 on 4 and 142 DF,  p-value: &lt; 2.2e-16

Cette fois, une des variables n’est pas significative (la première) donc on devrait l’enlever: on fait alors la régression juste sur les trois autres variables

reg356 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg356)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.11031    0.16788   54.27   &lt;2e-16 ***
X3           1.96784    0.07604   25.88   &lt;2e-16 ***
X5          -0.99296    0.08751  -11.35   &lt;2e-16 ***
X6           1.08761    0.05934   18.33   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: &lt; 2.2e-16

Cette fois, le modèle est “bon” et on peut alors faire la prévision

predict(reg356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.75388

A titre d’information, si on avait fait la prévision sans enlever la variable non significative, on aurait obtenu la valeur suivante

predict(reg1356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.90501

qui est globalement assez proche. On aura l’occasion d’en reparler en cours, quand on se demandera s’il est plus grave d’avoir une valeur non-significative dans la régression, ou d’oublier une variable importante… Mais encore une fois, le but de ces petits exercices est d’appliquer ce qu’on a vu en cours, et d’introduire des questions auxquelles j’apporterai des réponses au prochain cours !

Quelle responsabilité pour les algorithmes ?

Il y a quelques semaines, avec Rodolphe Bigot, maître de conférences à l’Université de Picardie Jules Verne, on avait commencé une réflexion sur le thème “Repenser la responsabilité, et la causalité“, et voici la suite…

Historiquement, les algorithmes se contentaient de fournir une aide à la décision, laissant à un être humain le rôle de prendre la décision, mais des expériences sont en cours, avec des systèmes autonomes, prenant des décisions, que ce soit les systèmes de conduite de voiture, ou les algorithmes de justice prédictive, comme le montre Huss et al. (2018). Cette autonomie, qui signifie fondamentalement la « faculté d’agir librement » désigne aussi l’idée « de se gouverner par ses propres lois ». Mais quelle est la responsabilité du décisionnaire dans le cas d’une prédiction qui entrainerait un préjudice ?

Continue reading Quelle responsabilité pour les algorithmes ?