Category Archives: Statistics

Régression sur des variables corrélées, un peu de géométrie

Depuis quelques jours, je mets en ligne des billets sur la régression sur des variables corrélés, en commençant par du deuxième effet kiss-cool suivi d’un court billet sur des régressions en cascade. Oui, depuis le premier, mon collègue Olivier ne cesse de demander des compléments. Alors ce soir, on va faire des dessins…

On a ici 3 vecteurs, \color{red}{\vec{x_1}}, \color{green}{\vec{x_2}} et \color{blue}{\vec{y}}. Pour rappel, la longueur de mes vecteurs, c’est la variance de la variable associée, \|{\color{red}{\vec{x_1}}}\|=\text{Var}({\color{red}{x_1}}). Et commençons par supposer que les variables \color{red}{x_1} et \color{green}{x_2} sont indépendantes (ou plus simplement non-corrélées) soit {\color{red}{\vec{x_1}}}\perp\color{green}{\vec{x_2}}. C’est le dessin ci-dessous,

La régression linéaire multiple est ici équivalente aux régressions simples, comme on en parlait dans les billets précédant. Géométriquement, on voit que les coordonnées de la projection de \color{blue}{\vec{y}} sur l’espace engendré par \color{red}{\vec{x_1}} et \color{green}{\vec{x_2}} (c’est à dire le plan (x,y)) correpond au point noir. C’est la régression multiple. Les deux régressions simples sont obtenus en projetant respectivement sur sur l’espace engendré par \color{red}{\vec{x_1}} (c’est à dire la droite (x)) et sur \color{green}{\vec{x_2}} (c’est à dire la droite (y)).

Supposons maintenant que \color{red}{\vec{x_1}} et \color{green}{\vec{x_2}} ne sont plus orthogonaux. Plus précisément, on va garder \color{red}{\vec{x_1}} et \color{blue}{\vec{y}} inchangés, et on va changer , \color{green}{\vec{x_2}}, tout en gardant la norme \|\color{green}{\vec{x_2}}\| inchangée. C’est ce qu’on a sur le schéma ci-dessous

On avait vu qu’on ne pouvait plus faire deux régressions indépendantes, mais qu’on pouvait garder la première, et en faire ensuite une seconde, mais en faisant un peu attention. Pour la première régression, on va projeter \color{blue}{\vec{y}} sur \color{red}{\vec{x_1}} et ca c’est facile, rien n’a changé,

Pour la seconde étape, rappelons qu’on doit projeter à la fois \color{blue}{\vec{y}} et \color{green}{\vec{x_2}} sur l’orthogonal de \color{red}{\vec{x_1}}, c’est à dire ici le plan (y,z), en rouge ci-dessous

(je renvois au précédant billet pour les formules algébriques). On obtient alors les projections suivantes, qu’on pourrait noter \Pi_{\color{red}{x_1^{\perp}}}\color{blue}{\vec{y}} et \Pi_{\color{red}{x_1^{\perp}}}\color{green}{\vec{x_2}}

On le voit, si on regarde sur le mur du fond, en rouge, la projection de \Pi_{\color{red}{x_1^{\perp}}}\color{blue}{\vec{y}} sur \Pi_{\color{red}{x_1^{\perp}}}\color{green}{\vec{x_2}} est un peu plus grande, autrement dit, le coefficient de la régression simple sera un peu plus grand. En fait, si on veut faire des calculs, on voit que le sinus de l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle va intervenir ici… Plus l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle sera petit, plus les variables seront corrélées, et plus la correction sera importante. Graphiquement, la correction dont je parle, c’est le petit segment, à gauche sur le dessin

Mais simulons quelques données pour visualiser tout ça…

library(mnormt)
n=10000
r=0
set.seed(1)
Z=cbind(rnorm(n)*sqrt(2),rnorm(n)*sqrt(3))
df=data.frame(y=Z[,1]+Z[,2]+rnorm(n)*sqrt(3),x1=Z[,1],x2=Z[,2])

On simule ici des variables explicatives indépendamment. Et on peut vérifier que les distances sont bien celles de notre dessin

var(df$x1)
[1] 2.049731
var(df$x2)
[1] 2.810419
var(df$y)
[1] 8.066665

et en plus, les variables sont quasiment orthogonales

cor(df$x1,df$x2)
[1] 0.004845078

Si on fait la régression multiple

reg=lm(y~0+x1+x2,data=df)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  0.99904    0.01219   81.94   2e-16 ***
x2  1.00090    0.01017   98.40   2e-16 ***

On a des coefficients qui sont identiques (à très peu de choses près, mais le soucis vient du fait que j’ai enlevé la constante du modèle pour coller au mieux avec mon dessin)

reg1=lm(y~0+x1,data=df)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  1.00489    0.01711   58.75   2e-16 ***

et pour la seconde régression aussi

reg2=lm(y~0+x2,data=df)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x2  1.00496    0.01315   76.42   2e-16 ***

On a déjà parlé de ce résultat longuement (je renvoie au précédant billet, des régressions en cascade). Regardons maintenant ce qui se passe si on garde nos observations \color{blue}{y_i}, \color{red}{x_{1,i}} mais que maintenant, la variable \color{red}{x_{2,i}} se retrouve corrélée avec \color{red}{x_{1,i}} , disons avec une corrélation de l’ordre de 0.4

r=.4
df$x2=r*df$x1+sqrt(1-r^2)*df$x2
cor(df$x1,df$x2)
[1] 0.3461496

Faire la régression multiple  donne ici

reg=lm(y~0+x1+x2,data=df)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1   0.5622     0.0130   43.26   2e-16 ***
x2   1.0921     0.0111   98.40   2e-16 ***

ce qui est différent des régressions simples, comme la première

reg1=lm(y~0+x1,data=df)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  1.00489    0.01711   58.75   2e-16 ***

Ce que nous dit le théorème de Frisch-Waugh, c’est qu’on peut faire des régressions en cascades: on peut garder cette première régression, mais pour la seconde, on va régresser les projections sur l’orthogonal à \color{red}{\vec{x_1}} (le mur en rouge au fond dans notre dessin). Ici, on obtient

reg2bis=lm(residuals(lm(y~x1,data=df))~residuals(lm(x2~x1,data=df)))
summary(reg2bis)
 
Coefficients:
           Estimate Std. Error t value Pr(|t|)      
residuals(2) 1.0921     0.0111    98.4   2e-16 ***

qui donne effectivement un coefficient de régression un peu supérieur à celui que nous avions auparavant, comme nous l’observions sur le dessin. En fait on peut remarquer que la différence relative s’obtient très facilement avec la corrélation

1/sqrt(1-r^2)
[1] 1.091089

et pour rappel, d’un point de vue géométrique, la corrélation est un cosinus d’un angle, r=\cos(\langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle). Autrement dit \sqrt{1-r^2}=\sin(\langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle), c’est à dire que la correction est directement liée au sinus de l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle, comme mentionné auparavant… amusant, non?

Des régressions en cascade

Cette fin de semaine, je mettais en ligne un court billet du deuxième effet kiss-cool où je rappelais que quand on fait une régression sur plusieurs variables explicatives corrélées, ce n’est pas équivalent à faire plusieurs régressions simple. C’est ce que disait le théorème de Frisch-Waugh, publié il y a près de 90 ans. Mais je glissais qu’il était possible de faire des régressions en cascade, sans aller beaucoup plus loin. Il faut relire l’article de Michael Lovell qui donner une interprétation géométrique de ce dessin. Mais programmons le, pour le voir… Considérons une base avec 3 variables explicatives, correspondant à des nombres d’incendies à Chicago

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

La régression multiple s’écrit iciy_i=\beta_0+\beta_1 x_{1,i}+\beta_2 x_{2,i}+\beta_3 x_{3,i}+\varepsilon_iOn peut estimer ces paramètres par moindre carrés,

reg=lm(Fire~.,data=chicago)
summary(reg)
 
Coefficients:
            Estimate Std. Error t value  Pr(|t|)    
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X_1         -0.62764    5.28130  -0.119 0.905953    
X_2          0.22378    0.06161   3.632 0.000744 ***
X_3         -1.55059    0.38195  -4.060 0.000204 ***

et on a alors la prévision suivante
\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 x_{1,i}+\widehat{\beta}_2 x_{2,i}+\widehat{\beta}_3 x_{3,i}En fait, je peux même visualiser cette prévision.

BC = cbind(x0=coefficients(reg)[1],
           x1=coefficients(reg)[2]*chicago[,2],
           x2=coefficients(reg)[3]*chicago[,3],
           x3=coefficients(reg)[4]*chicago[,4],
           yp=predict(reg),y=chicago$Fire)
colr = 1:4
dessin_fire = function(i,BC=BC){
plot(0:3,0:3,xlim=c(0,35),col="white",ylim=c(-1,3),xlab="",ylab="",axes=FALSE)
abline(v=BC[i,1],lty=2)
rect(BC[i,1],3,BC[i,1]+BC[i,2],2,col=colr[1],border=NA)
arrows(BC[i,1],2.5,BC[i,1]+BC[i,2],2.5,lwd=2,length=.1,col="white")
rect(BC[i,1]+BC[i,2],2,BC[i,1]+BC[i,2]+BC[i,3],1,col=colr[2],border=NA)
arrows(BC[i,1]+BC[i,2],1.5,BC[i,1]+BC[i,2]+BC[i,3],1.5,lwd=2,length=.1,col="white") 
rect(BC[i,1]+BC[i,2]+BC[i,3],1,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],0,col=colr[3],border=NA)
segments(BC[i,5],0,BC[i,5],1,lwd=3)
arrows(BC[i,1]+BC[i,2]+BC[i,3],.5,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],.5,lwd=2,length=.1,col="white") 
abline(v=BC[i,1],lty=2)
rect(BC[i,5],-1,BC[i,6],0,col=colr[4],density=20,border=NA)
segments(BC[i,6],0,BC[i,6],-1,lwd=3,col=colr[4])
text(26,2.5,expression(X[1]),pos=4,col=colr[1])
text(26,1.5,expression(X[2]),pos=4,col=colr[2])
text(26,.5,expression(X[3]),pos=4,col=colr[3])
text(26,-.5,expression(epsilon),pos=4,col=colr[4])
axis(1)}

Par exemple, pour la douzième observation de ma base

dessin_fire(12)

on obtient

On lit de haut en bas : on commence par constante, vers 22. Puis on rajoute trois fois rien, parce que x_1 était non-significative dans notre régression. Puis on rajoute un petit quelque chose à cause de x_2, pour passer à 25, et x_3 nous fait plonger de plus de 20 points, pour finir vers 3.67.Le trait noir est la prévision que l’on obtient, à partir des trois variables. En bleu, on peut même voir l’erreur.

On peut noter que x_1 est non-significative, tout en étant relativement corrélée avec y

cor(chicago[,1:2])
          Fire       X_1
Fire 1.0000000 0.3773486
X_1  0.3773486 1.0000000

Car si x_1 est non-significative, c’est parce que cette variable est très corrélée avec une autre variable explicative. En fait, on peut faire une série de régressions en cascade, en commençant par x_1, et plus spécifiquement, on peut écrire\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_1 x_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,2}+\widehat{b}_2 \tilde{x}_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}où le premier terme (1) est obtenu en régression simplement y sur x_1 (oui, on fait juste une régression simple),  y_i=b_0+b_1 x_{1,i}+\eta_ipuis on apporte une petite correction, pour tenir compte de ce que nous apprend x_2 une fois la première régression effectuée. C’est exactement ce que raconte le théorème de Frish-Waugh : on va projeter y sur x_1^{\perp} (car c’est justement ce que x_1 n’explique pas), et x_2 sur x_1^{\perp}  et faire la régression de ces deux projections, soit, pour (2) \Pi_{x_1^{\perp}}y_i=b_{0,2}+b_2 \Pi_{x_1^{\perp}}x_{2,i}+\eta_{2,i}puis pour (3) \Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i}En terme informatique, cela donne

reg1=lm(Fire~X_1,data=chicago)
reg2=lm(residuals(lm(Fire~X_1,data=chicago))~residuals(lm(X_2~X_1,data=chicago)))
reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago)))

soit ici

BC123=cbind(x0=coefficients(reg1)[1],
            x1=coefficients(reg1)[2]*chicago[,2],
            x2=coefficients(reg2)[1]+coefficients(reg2)[2]*residuals(lm(X_2~X_1,data=chicago)),
            x3=coefficients(reg3)[1]+coefficients(reg3)[2]*residuals(lm(X_3~X_1+X_2,data=chicago)),
            yp=predict(reg),y=chicago$Fire)

En adaptant la fonction de dessin précédante on obtient

dessin_fire_123(12)

Autrement dit, on commence avec une régression simple, en rouge: on part d’une constante vers 3.5, puis on augemente notre prévision, en tenant compte de x_1. Puis on corrige avec x_2 (ou plutôt la partie de x_2 non expliquée par x_1). sur ce qui n’avait pas été expliqué auparavant, ce qui pousse à revoir un peu à la baisse notre prévision, puis on tient compte de x_3 (ou là encore, sa projection sur l’orthogonal de x_1 et x_2). On note que la prévision est très exactement la même qu’auparavant.

Mais on peut aller plus loin… pourquoi commencer par x_1? On peut aussi commencer par x_2. On considère alors\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_2 x_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,1}+\widehat{b}_1 \tilde{x}_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}où le premier terme (2) est obtenu en régression simplement y sur x_2y_i=b_0+b_2 x_{2,i}+\eta_ipuis on apporte une petite correction, pour tenir compte de ce que nous apprend x_1 une fois la première régression effectuée, (1) \Pi_{x_2^{\perp}}y_i=b_{0,1}+b_1\Pi_{x_2^{\perp}}x_{1,i}+\eta_{1,i}puis pour (3) \Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i} Soit

reg1=lm(Fire~X_2,data=chicago)
reg2=lm(residuals(lm(Fire~X_2,data=chicago))~residuals(lm(X_1~X_2,data=chicago)))
reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago)))
BC213=cbind(x0=coefficients(reg1)[1],
            x1=coefficients(reg1)[2]*chicago[,3],
            x2=coefficients(reg2)[1]+coefficients(reg2)[2]*residuals(lm(X_1~X_2,data=chicago)),
            x3=coefficients(reg3)[1]+coefficients(reg3)[2]*residuals(lm(X_3~X_1+X_2,data=chicago)),
            yp=predict(reg),y=chicago$Fire)

et visuellement, on obtient

dessin_fire_213(12)

Cette fois, la constante est un peu plus élevée, et on commence par faire une régression seulement avec x_2, pour arriver un peu en dessous de 10. Puis on corrige. On note que cette seconde correction nous ramène… très exactement comme dans le cas précédant. Mais si on pense en projection successives, on ne doit pas être surpris.

On peut récapituler avec une autre prévision, soit avec le modèle de régresion multiple (mais là encore, x_1 n’explique pas grand chose)

dessin_fire(15)

soit en commencant par une régression simple sur x_1 puis en faisant des régression en cascade,

dessin_fire_123(15)

soit, alternativement, en commencant par une régression simple sur x_2

dessin_fire_213(15)

Comme annoncé, cette trois approches sont équivalentes, et donnent très exactement la même prévision.

Regression discontinuity model for TV series

In September, we are usually happy to see our favorite TV series back on air… Or not? Because, admit it, if we are happy to see those characters back, most of the time, we are disappointed, too. So why not look at the data, to confirm this feeling? Nazareno Andrade shared some nice codes to get IMDB ratings in a nice csv file (you can either use the large csv file, or run your own codes)

download.file("https://github.com/nazareno/imdb-series/raw/master/data/series_from_imdb.csv",
destfile="series_from_imdb.csv")
base = read.csv("series_from_imdb.csv")

It is a large dataset, with more than 64,000 episodes of almost 890 TV series,

str(base)
'data.frame':	64018 obs. of  18 variables:
 $ series_name: Factor w/ 889 levels "'Allo 'Allo!",..: 137 137 137 137 137 137 137 137 137 137 ...
 $ episode    : Factor w/ 54090 levels "-30-","¡Viva los muertos!",..: 32314 7446 16 7176 17748 9562 1379 36218 17845 5553 ...
 $ series_ep  : int  1 2 3 4 5 6 7 8 9 10 ...
 $ season     : int  1 1 1 1 1 1 1 2 2 2 ...
 $ season_ep  : int  1 2 3 4 5 6 7 1 2 3 ...
 $ user_rating: num  8.9 8.7 8.7 8.2 8.3 9.2 8.8 8.7 9.2 8.3 ...

Just pick a TV series, for instance Dan Harmon’s Community,

sbase = base[base$series_name=="Community",]

We can plot the evolution of the rating over the 110 episodes.

sbase=sbase[!duplicated(sbase[,c(1,2,4,5)]),]
sbase$series_ep=1:nrow(sbase)

()since there could be some problem with the data (such as duplicates, let us clean it quickly)

plot(sbase$series_ep,sbase$UserRating,xlab=sbase$series_name[1])
idx=c(0,which(diff(sbase$season)!=0),nrow(sbase))
abline(v=idx+.5,lty=2,col=colr[2])
a = unique(sbase$season)
for(u in a){
  ssbase = sbase[sbase$season==u,]
  reg = lm(UserRating~series_ep,data=ssbase)
  lines(ssbase$series_ep,predict(reg),col=colr[3],lwd=2)
}

The vertical lines are here to visualize the seasons. On issue is that the lenght can vary with time. Consider Linwood Boomer’s Malcom in The Middle,

sbase = base[base$series_name=="Malcolm in the Middle",]

or Craig Thomas and Carter Bays’s How I Met Your Mother,

sbase = base[base$series_name=="How I Met Your Mother",]

On those two, the evolution is rather stable. Look at AMC’s The Walking Dead,

sbase = base[base$series_name=="The Walking Dead",]

Now, look at Howard Gordon and Alex Gansa’s Homeland,

sbase = base[base$series_name=="Homeland",]

There is an issue here with the last episode of season4, “Long Time Coming“, that has a very poor rating. If we remove that point, we get the thin line. Note that the regression line is always increasing. For Michael Hirst’s Vickings, we have

sbase = base[base$series_name=="Vicking",]

If we look more carefully on the previous graph, for five seasons (out of six), we have a positive slope. Well, to be honest, it is not significantly positive most of the time, but still. Out of 80 shows, and a total of 583 seasons, the slope is postive 75% of the time (433) and negative 25% of the time (150).

BASE = NULL
L80 = unique(base$series_name)
for(j in 1:length(L)){
sbase=base[base$series_name==L[j],]
sbase=sbase[!duplicated(sbase[,c(1,2,4,5)]),]
sbase=sbase[sbase$season>0,]
sbase$series_ep=1:nrow(sbase)
a=unique(sbase$season)
a=a[!is.na(a)]
for(u in a){
  ssbase=sbase[sbase$season==u,]
  reg=lm(UserRating~series_ep,data=ssbase)
  pente = NA
  if((!is.na(coefficients(reg)[2]))&(!is.na((summary(reg)$coefficients[2,4])))){
  if((summary(reg)$coefficients[2,4]<.05)&(coefficients(reg)[2]>0)) pente="positive"
  if((summary(reg)$coefficients[2,4]<.05)&(coefficients(reg)[2]<0)) pente="negative" sdf=data.frame(nom=sbase$series_name[1],season=u,slope=coefficients(reg)[2],inf=confint(reg)[2,1],sup=confint(reg)[2,2],signe=pente) BASE=rbind(BASE,sdf)} }} str(BASE) 'data.frame': 583 obs. of 6 variables: $ nom : Factor w/ 80 levels "Friends","Game of Thrones",..: 1 1 1 1 1 1 1 1 1 1 ... 
 
mean(BASE$slope>0)
[1] 0.7427101
table(BASE$signe)
negative positive 
      15      144

Most of the time, the slope is not significant. To be more specific, 72% of the time, the slope is not significant. But when it is, 90% of the time, it is positive (144 seasons). Let us look at other TV series, for instance Joel Surnow and Robert Cochran’s 24,

sbase = base[base$series_name=="24",]

Álex Pina’s La Casa de Papel,

sbase = base[base$series_name=="La Casa de Papel",]

Steven Knight’s Peaky Blinders,

sbase = base[base$series_name=="Peaky Blinders",]

or David Simon’s The Wire,

sbase = base[base$series_name=="The Wire",]

The slope is increasing over almost all seasons. But a major drawback is that when we get back to our show, for a new season, we usually get disapointed. More specifically, we can quantify the difference in red below

that can be estimated using

sbase12 = sbase[sbase$season%in%c(a[ij],a[ij+1]),]
seuil = sbase12$series_ep[which(diff(sbase12$season)!=0)]+.5
s = function(x) (x-seuil)*(x>seuil)
reg = lm(UserRating~series_ep+s(series_ep)+I(series_ep>seuil),data=sbase12)

Here we have

summary(reg)
Coefficients:
                         Estimate Std. Error t value  Pr(|t|)    
(Intercept)               8.45000    0.16338  51.719    2e-16 ***
series_ep                 0.10000    0.03235   3.091 0.008598 ** 
s(series_ep)              0.02000    0.04218   0.474 0.643291    
I(series_ep)TRUE.        -1.01778    0.20486  -4.968 0.000257 ***

so the drop of 1 point (out of 10) cannot be claimed as being significant. That is the idea of regression discontinuity.

If we loop again over all our series, we have 485 pairs of consecutive seasons. As expected, in 75% of the casse, from season t-1 to season t, we observe a negative rupture. As previously, in 70% of the cases, it is not significat (with linear models before and after), and when it is significant, it is negative in 96% of the cases ! But an alternative can be to use nonparametric models, on both sides.

To illustrate, consider David Benioff and D. B. Weiss’s Game of Thrones,

sbase = base[base$series_name=="Game of Thrones",]

but let us remove the last season (no spoiler here, but clearly not worst watching)

Consider for instance the drop between season 1 and season 2,

library(rdd)
sbase12=sbase[sbase$season%in%c(1,2),]
lmr=RDestimate(UserRating~series_ep,data=sbase12,cutpoint=mean(range(sbase12$series_ep)))
plot(lmr)

This is very consistent with what we observed with our linear regressions actually,

seuil=10.5
s = function(x) (x-seuil)*(x>seuil)
reg = lm(UserRating~series_ep+s(series_ep)+I(series_ep>seuil),data=sbase12)
summary(reg)
 
Coefficients:
                         Estimate Std. Error t value  Pr(|t|)    
(Intercept)               8.70000    0.15458  56.281    2e-16 ***
series_ep                 0.07273    0.02491   2.919  0.01003 *  
s(series_ep)              0.01455    0.03523   0.413  0.68520    
I(series_ep)TRUE         -0.94000    0.20316  -4.627  0.00028 ***

Here, the drop of one point is significant…

So, your favorite show had an outstanding finale ? and you can’t wait to watch the new season… Well, statistically, it’s very likely that you will be disapointed by the first episode of the forthcoming season…

Du deuxième effet kiss-cool (régression multiple, scoring et évaluation)

Lorsque j’étais petit (il y a fort longtemps, à une époque où je regardais pas mal la télévision) il y avait une publicité pour les pastilles kiss cool,

Et quand je présente la régression multiple à mes étudiants, je ne peux m’empêcher d’y penser… Mais avant d’aller plus loin sur le parallèle, faisons un peu de mathématiques.

Les techniques de régression permettent d’avoir des jolis théorèmes, souvent d’une portée incroyablement générale (moyennant quelques petites hypothèses techniques). Par exemple le théorème de Frisch-Waugh, en régression multiple, dont j’ai déjà parlé dans des vieux billets. Un des corolaires est que lorsque les variables explicatives dans un modèle de régression sont orthogonales, la régression multiple correspond à une collection de régressions simples (autrement dit, les estimateurs par moindres carrés coïncident). Formellement, si on considère le modèley_i=\beta_0+\beta_1x_{1,i}+\beta_2x_{2,i}+\varepsilon_i(avec les hypothèses usuelles des modèles de régression) alors, si les variables x_1 et x_2 sont non-corrélées, \widehat{\beta}_1 coïncide avec \widehat{b}_1 dans le modèley_i=b_0+b_1x_{1,i}+\eta_iOn peut faire une petite simulation pour confirmer (sinon, bien entendu, on peut regarder la démonstration qui se trouve dans tous les livres d’économétrie, qui est d’ailleurs un simple résultat d’algèbre linéaire, ou de géométrie, avec des projections successives sur des sous-espaces orthogonaux – même si c’est à Michael Lovell que l’on doit l’approche géométrique).

library(mnormt)
r = 0
S = matrix(c(1,r,r,1),2,2)
n = 1000
set.seed(1)
X = rmnorm(n,c(0,0),S)
E = rnorm(n,0,.3)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))

Petite note technique: je vais centrer les variables, histoire de ne pas avoir à garder la constante dans mon modèle (qui va compliquer les notations, et potentiellement embrouiller un peu le billet). La constante est nulle ici, on le voit,

reg = lm(Y~X1+X2,data=base)
summary(reg)
 
Coefficients:
              Estimate Std. Error t value  Pr(|t|)    
(Intercept)  7.449e-18  9.777e-03     0.0        1    
X1           1.012e+00  9.171e-03   110.3    2e-16 ***
X2          -1.988e+00  9.719e-03  -204.6    2e-16 ***

Bref, je régresse sans constante

reg = lm(Y~0+X1+X2,data=base)
summary(reg)
 
Coefficients:
    Estimate Std. Error t value  Pr(|t|)    
X1  1.011520   0.009166   110.3   2e-16 ***
X2 -1.988321   0.009714  -204.7   2e-16 ***

Maintenant, on va regarder les deux régressions simple

reg1 = lm(Y~0+X1,data=base)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X1  1.01300    0.06006   16.86   2e-16 ***

quand on régresse juste sur la première variable, et pour la seconde, on obtient

reg2 = lm(Y~0+X2,data=base)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X2 -1.98916    0.03528  -56.39   2e-16 ***

Autrement dit, nos estimateurs sont très proches (si on avait laissé la constante, ils coïncideraient).

Maintenant, le gros soucis est que ce résultat n’est plus valide lorsque les variables explicatives sont corrélées. Le théorème de Frisch-Waugh nous explique comment ces estimateurs divergent, mais le point ici est que si les variables sont corrélées, utiliser les régressions simples donne deux estimateurs biaisés des vrais paramètres (du modèle multiple). Recommençons l’exercice précédant

r=.9
S=matrix(c(1,r,r,1),2,2)
set.seed(1)
X=rmnorm(n,c(0,0),S)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))
reg = lm(Y~X1+X2,data=base)
summary(reg)
reg = lm(Y~0+X1+X2,data=base)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X1  0.98740    0.02205   44.79   2e-16 ***
X2 -1.97321    0.02229  -88.54   2e-16 ***

(on retrouve des valeurs proches de celles utilisées pour simuler nos données, donc tout va bien). En revanche, pour les régressions simples, on obtient

reg1 = lm(Y~0+X1,data=base)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X1 -0.78784    0.02726   -28.9   2e-16 ***

et

reg2 = lm(Y~0+X2,data=base)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X2 -1.06543    0.01607  -66.31   2e-16 ***

Autrement dit, \widehat{b}_1\neq \widehat{\beta}_1 (et pareil pour le second). Ce qui signifie que si on construit une prévision à partir de ce modèle, \widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}, on sera potentiellement très loin de la “bonne” prévision \widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i} (qui sera sans biais, etc, je renvoie ici vers n’importe quel cours d’économétrie). On peut le voir sur un dessin,

Yp=reg1$coefficients[1]*base$X1+reg2$coefficients[1]*base$X2
plot(base$Y,predict(reg),ylim=range(Yp),col=rgb(0,0,1,.5),cex=.7,xlab="observé",ylab="prédit")
abline(a=0,b=1,lty=2)
points(base$Y,Yp,col=rgb(1,0,0,.5),cex=.7)
abline(lm(predict(reg)~base$Y),col="blue")
abline(lm(Yp~base$Y),col="red")

En bleu, on a les prévisions avec le modèle linéaire multiple, et en rouge, en faisant deux régressions indépendantes… Si on regarde sur la droite, le modèle rouge sur-valorise, ou disons sur-estime largement, alors qu’il sous-valorise à gauche. La différence entre les deux droites s’interprète ici comme un biais. Sur le graphique ci-dessous, on peut visualiser la distribution des \widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}, en rouge, et les \widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i} , en bleu. Cet excès de dispersion, de variance, qu’on observe sur les points rouges, j’interprète ça comme de la polarisation

plot(density(Yp),col="red")lwd=2)
lines(density(predict(reg)),col="blue",lwd=2)

En fait, ce que raconte le théorème de Frisch-Waugh (et je renvoie à mon précédant billet pour plus de détails), c’est qu’on a le droit de faire plusieurs régressions, mais en cascade ! et surtout pas indépendamment : je peux expliquer y avec la première variable x_1, et ensuite régresser le résidu (ce qu’on n’a pas pu expliquer) sur la seconde x_2. Cette méthode donnera la même prévision que le modèle multiple.

On peut aller un peu plus loin: on peut jouer sur la valeur de la corrélation, pour mesure l’écart entre les deux prévisions (ici je prévois pour une observation au hasard, la 78ème)

comp=function(r){
S=matrix(c(1,r,r,1),2,2)
set.seed(1)
X=rmnorm(n,c(0,0),S)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))
reg = lm(Y~0+X1+X2,data=base)
reg1 = lm(Y~0+X1,data=base)
reg2 = lm(Y~0+X2,data=base)
y1=predict(reg)
y2=reg1$coefficients[1]*base$X1+reg2$coefficients[1]*base$X2
c(y1[78],y2[78],(y2[78]-y1[78])/y1[78])}
vR=seq(0,.98,by=.02)
vc=Vectorize(comp)(vR)
plot(vR,vc[3,]*100,ylab="Différence relative (%)",xlab="Corrélation",type="l")

On prédit ici pour une observation avec un large y_i, ce qui correspondait à la partie de droite du graphique précédant (avec les points rouges et bleus).

C’est ce que j’appelle le deuxième effet kiss-cool. Quand on est dans le premier cas, avec les variables indépendantes (corrélation nulle, c’est à dire à gauche sur ma figure), on explique ce qu’on peut avec la première variable, et on rajoute l’impact de la seconde. Et \widetilde{y}_i\sim\widehat{y}_i. Le soucis ici est que je n’ai pas le droit de considérer deux modèles indépendants lors que les variables sont très corrélées. Une partie de l’explication fournie par la seconde variable était déjà inclue dans la première. Par exemple avec deux variables très (positivement) corrélées, la prévision qu’on obtient en ajoutant les deux effets estimés indépendamment avec deux régressions simple \widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}, on sur-estime de 50% à 70% la “vraie” prévision \widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i}.

Tous les chercheurs savent savent ça… on parle ici de résultats du tout premier cours de modèles linéaires. Et malgré tout, en pratique, on continue à utiliser cette seconde méthode. Un exemple bien connu est celui de l’évaluation (des étudiants, des chercheurs, peu importe). Par exemple, quand on évalue un dossier de financement pour un chercheur, on nous demande de mettre un score

  • pour les publications scientifiques (nombre, qualité, etc)
  • pour l’encadrement d’étudiants (nombre, niveau, etc)
  • pour la qualité de l’environnement (prestige du labo, etc)
  • etc

Et à la fin, on somme tout. Mais on le voit, ces variables sont très très corrélées: si vous êtes dans un labo prestigieux, vous attirez beaucoup de candidatures d’étudiants (et des bons), et avoir beaucoup d’étudiants va permettre d’avoir plus de publications (si on ajoute son nom comme co-auteur). Bref, on est typiquement dans un modèle à double (voire triple) effet kiss-cool. Quelqu’un dans un bon labo aura un bon score sur le troisième item, mais aussi un bon sur le nombre d’étudiants, et aussi sur les publications. Ajouter ces scores est stupide, car on a une spirale infernale (les bons sont sur-évalués, et les moins bon, sous-évalués), c’est ce que racontait mon premier dessin, avec les points rouges et bleus. C’est un effet clivant de polarisation forte.

Si on voulait faire les choses proprement, ce que dit le le théorème de Frisch-Waugh, c’est que les scores devraient être attribués en corrigeant de la corrélation entre les variables

  • on peut commencer par calculer un score pour la qualité de l’environnement (prestige du labo, etc)
  • à environnement donné, calculer un score pour l’encadrement d’étudiants (nombre, niveau, etc)
  • à environnement donné, et à encadrement d’étudiants donné, calculer un score pour les publications
  • etc

C’est comme la situation que je voyais en France, où on pouvait avoir une variable qui tenait compte d’un prestige du chercheur (par exemple, être chercheur CNRS donnait un bonus) et une autre sur le dossier de publications. Sauf que les deux sont corrélés. Et la plupart des classements d’universités sont construits à partir de scores qui sont loin d’être indépendants.

Bref, tant que l’évaluation se fera en sommant des scores qui sont construits sur des critères souvent très corrélés, on polarise fortement la population.

Est-ce gênant? A priori oui. Car le message que cela envoie est qu’il existe deux classes, les bons, et les mauvais, alors qu’en réalité, le niveau est beaucoup plus homogène qu’il n’y paraît.  Un petit effet positif se retrouve démultiplié par le fait qu’il va se répercuter (positivement) sur plein d’autres variables. C’est mon effet kiss-cool. Mais on pourrait se dire que c’est un problème de distribution des notes finales. Si l’ordre est préservé, on pourrait se dire que ce n’est pas très grave. Malheureusement, ce n’est pas le cas.

Si on quitte un instant le cas de la corrélation très forte, les rangs des prédictions (c’est à dire les rangs des chercheurs une fois donnés les notes \widehat{y}_i ou \widetilde{y}_i ) sont moins corrélés si la corrélation sous-jacente est importante (mais pas trop)

comp=function(r){
S=matrix(c(1,r,r,1),2,2)
set.seed(1)
X=rmnorm(n,c(0,0),S)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))
reg = lm(Y~0+X1+X2,data=base)
reg1 = lm(Y~0+X1,data=base)
reg2 = lm(Y~0+X2,data=base)
y1=predict(reg)
y2=reg1$coefficients[1]*base$X1+reg2$coefficients[1]*base$X2
cor(y1,y2,method="spearman")}
vR=seq(0,.98,by=.02)
vc=Vectorize(comp)(vR)
plot(vR,vc,ylab="Corrélation de rangs",xlab="Corrélation",type="l")

Autrement dit, si les variables x_1 et x_2 sont très peu corrélées, on a les mêmes rangs  (globalement). En revanche, si la corrélation entre les variables x_1 et x_2 augmente, le rang des \widetilde{y}_i est de moins en moins cohérent avec celui entre les \widehat{y}_i  (qui devrait être celui que l’on recherche).

Bref, il serait temps de comprendre enfin sérieusement les conséquences de ce joli papier, publié il y a presque 90 ans

Prédire le nombre de morts, suite

Il y a quelques jours, un tweet de Baptiste Coulmont m’inspirait un billet sur l’estimation de la surmortalité, en France. Ce soir, c’est un tweet d’Alexandre Blanchet, qui parlait de surmortalité, au Québec.

J’étais surpris, parce que la prévision est bien au dessus des années passés contrairement à ce qu’on pouvait observer en France. En repartant des codes d’Alexandre (data_prep.R puis predictions.R) j’ai remis un peu en forme le graphique, pour avoir la vision de long-terme, pour comprendre la tendance linéaire utilisée par Alexandre (qui utilise un modèle à la Buys-Ballot)

Je me disais qu’on pouvait tenter des modèles de lissages exponentiel, à la Holt-Winters,

library(forecast)
x=ts(df.l_training$deces,start = 2010, frequency = 12)
x.hw <- ets(x, model = "AAA")
autoplot(forecast(x.hw,24))

en version additive, ou multiplicative

x.hw2 <- ets(x, model = "MAM")
autoplot(forecast(x.hw2,24))

Je laisse à d’autre le soin de comparer aux nombres de morts observés depuis 2 mois, je reste fasciné pour ma part par la tendance à la hausse de la mortalité, depuis 10 ans… Et j’essayerais de comparer ces trois modèles dès que j’aurais finis mes corrections de copies de la session d’hiver.

Testing for Covid-19 in the U.S.

For almost a month, on a daily basis, we are working with colleagues (Romuald, Chi and Mathieu) on modeling the dynamics of the recent pandemic. I learn of lot of things discussing with them, but we keep struggling with the tests. Paul, in Montréal, helped me a little bit, but I think we will still have to more to get a better understand. To but honest, we stuggle with two very simple questions

  • how many people are tested on a daily basis ?

Recently, I discovered Modelling COVID-19 exit strategies for policy makers in the United Kingdom, which is very close to what we try to do… and in the document two interesting scenarios are discussed, with, for the first one, “1 million ‘reliable’ daily tests are deployed” (in the U.K.) and “5 million ‘useless’ daily tests are deployed”. There are about 65 millions unhabitants in the U.K. so we talk here about 1.5% people tested, on a daily basis, or 7.69% people ! It could make sense, but our question was, at some point, is that realistic ? where are we today with testing ? In the U.S. https://covidtracking.com/ collects interesting data, on a daily basis, per state.

url = "https://raw.githubusercontent.com/COVID19Tracking/covid-tracking-data/master/data/states_daily_4pm_et.csv"
download.file(url,destfile="covid.csv")
base = read.csv("covid.csv")

Unfortunately, there is no information about the population. That we can find on wikipedia. But in that table, the state is given by its full name (and the symbol in the previous dataset). So we new also to match the two datasets properly,

url="https://en.wikipedia.org/wiki/List_of_states_and_territories_of_the_United_States_by_population"
download.file(url,destfile = "popUS.html")
#pas contaminé 2/3 R=3
library(XML)
tables=readHTMLTable("popUS.html")
T=tables[[1]][3:54,c("V3","V4")]
names(T)=c("state","pop")
url="https://en.wikipedia.org/wiki/List_of_U.S._state_abbreviations"
download.file(url,destfile = "nameUS.html")
tables=readHTMLTable("nameUS.html")
T2=tables[[1]][13:63,c(1,4)]
names(T2)=c("state","symbol")
T=merge(T,T2)
T$population = as.numeric(gsub(",", "", T$pop, fixed = TRUE))
names(base)[2]="symbol"
base = merge(base,T[,c("symbol","population")])

Now our dataset is fine… and we can get a function to plot the number of people tested in the U.S. (cumulated). Here, we distinguish between the positive and the negative,

drawing = function(st ="NY"){
sbase=base[base$symbol==st,c("date","positive","negative","population")]
sbase$DATE = as.Date(as.character(sbase$date),"%Y%m%d")
sbase=sbase[order(sbase$DATE),]
par(mfrow=c(1,2))
plot(sbase$DATE,(sbase$positive+sbase$negative)/sbase$population,ylab="Proportion Test (/population of state)",type="l",xlab="",col="blue",lwd=3)
lines(sbase$DATE,sbase$positive/sbase$population,col="red",lwd=2)
legend("topleft",c("negative","positive"),lwd=2,col=c("blue","red"),bty="n")
title(st)
plot(sbase$DATE,sbase$positive/(sbase$positive+sbase$negative),ylab="Ratio of positive tests",ylim=c(0,1),type="l",xlab="",col="black",lwd=3)
title(st)}

Let us start with New York

drawing("NY")

As at now, 4% of the entiere population got tested… over 6 weeks…. The graph on the right is the proportion of people who tested positive… I won’t get back on that one here today, I keep it for our work. In New Jersey, we got about 2.5% of the entiere population tested, overall,

drawing("NJ")

Let us try a last one, Florida

drawing("FL")

As at today, it is 1.5% of the population, over 6 weeks. Overall, in the U.S. less than 0.1% people are tested, on a daily basis. Which is far from the 1.5% in the U.K. scenarios. Now, here come the second question,

  • what are we actually testing for ?

On that one, my experience in biology is… very limited, and Paul helped me. He mentioned this morning a nice report, from a lab in UC Berkeley

One of my question was for instance, if you get tested positive, and you do it again, can you test negative ? Or, in the context of our data, do we test different people ? are some people tested on a regular basis (perhaps every week) ? For instance, with antigen tests (Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) – also called molecular or PCR – Polymerase Chain Reaction – test) we test if someone is infectious, while with antibody test (using serological immunoassays that detect viral-specific antibodies — Immunoglobin M (IgM) and G (IgG) — also called serology test), we test for immunity. Which is rather different…

I have no idea what we have in our database, to be honest… and for the past six weeks, I have seen a lot of databases, and most of the time, I don’t know how to interpret, I don’t know what is measured… and it is scary. So, so far, we try to do some maths, to test dynamics by tuning parameters “the best we can” (and not estimate them). But if anyone has good references on testing, in the context of Covid-19 (for instance on specificity, sensitivity of all those tests) I would love to hear about it !

Estimer la surmortalité

Ce matin, Baptiste Coulmont m’envoyait un tweet avec un joli graphique, présentant le nombre quotidien de décès en France.

Comme les données sont en ligne, je me suis dit que je pourrais jouer un peu avec. Pour les plus curieux, on a la liste de tous les décès depuis… longtemps ! (on a plus de 25 millions de décès). La seule chose qui nous intéresse, c’est la date, alors on la récupère. Puis on compte, par jour, combien de décès on a eu,

D = read.csv("insee_deces.csv",header=TRUE)
Vecteur_Dates = as.character(D[,8])
TV = table(Vecteur_Dates)
B = data.frame(dateC=names(TV),nb=as.numeric(TV))
B$date=as.Date(B$dateC,"%Y-%m-%d")
B$year = format(B$date,"%Y")
B$month = format(B$date,"%m")
B$day = format(B$date,"%d")
B$deb=as.Date(paste(B$year,"01-01",sep="-"),"%Y-%m-%d")
B$fin=as.Date(paste(B$year,"12-31",sep="-"),"%Y-%m-%d")
B$dif=as.numeric(B$date-B$deb)/as.numeric(B$fin-B$deb)

Je crée ici une une variable qui me dit où je suis dans l’année en cours (0 pour le 1er janvier, et 1 pour pour 31 décembre). Ça me permet de contourner un petit problème technique : les années bissextiles. On peut ensuite représenter toutes les années,

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",xlim=0:1,ylim=c(1200,4000),col="white")
for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")

On veut ensuite avoir la tendance moyenne. Baptiste prenait la moyenne par jour (la courbe noire). Ici, on va lisser avec des splines

library(gam)
reg = gam(NormCpte~bs(dif,40),data=B2[B2$year != "2020",])
vx = seq(0,1,length=501)
vy = predict(reg,newdata=data.frame(dif=vx))
lines(vx,vy,col="blue",lwd=3)

ce qui donne la courbe bleue. C’est l’écart à cette courbe qui donne la surmortalité. Par exemple (prenons un exemple bien connu, et largement visible sur le graphique précédant) l’année 2003, et les quinze premiers jours d’août, correspondant (plus ou moins) à la période de canicule.

i = 2003
x = B2$dif[(B2$year == i)&((B2$month == "08"))&(B2$day %in% c(paste("0",1:9,sep=""),10:15))]
y = B2$NormCpte[(B2$year == i)&((B2$month == "08"))&(B2$day %in% c(paste("0",1:9,sep=""),10:15))]
yp = predict(reg,newdata=data.frame(dif=x))
e = y-yp
sum(e)
[1] 14294.39

Plusieurs chiffres sont mentionnés sur wikipédia, dont le rapport de l’INSERM qui annonçait 14,802 morts. On n’est pas trop loin… On peut visualiser cette surmortalité sur le graphique ci-dessous

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",col="white",xlim=0:1,ylim=c(1200,4000))
for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")
lines(vx,vy,col="blue")
lines(B2$dif[B2$year == 2003],B2$NormCpte[B2$year == 2003],col="red")
for(u in 1:length(x)) segments(x[u],y[u],x[u],yp[u],lwd=3,col="red")

Prenons un autre exemple qu’on voit bien sur le graphique de Baptiste, la grippe de 2017 (dite de 2016-2017). Si on prend les 45 premiers jours de 2017, on obtient les chiffres suivants

i = 2003
x=B2$dif[(B2$year == i)&((B2$month == "01")|
  ((B2$month == "02")&(B2$day %in% c(paste("0",1:9,sep=""),10:15))))]
y=B2$NormCpte[(B2$year == i)&((B2$month == "01")|
  ((B2$month == "02")&(B2$day %in% c(paste("0",1:9,sep=""),10:15))))]
yp=predict(reg,newdata=data.frame(dif=x))
yp = predict(reg,newdata=data.frame(dif=x))
e = y-yp
sum(e)
[1] 21177.33

Le rapport de Santé Publique France parle de 21,200 décès (pour la “surmortalité toute cause”) ce qui est là aussi comparable…

Maintenant, avant de conclure, et avant de me faire troller pendant quelques jours, je tiens à préciser que c’est juste un exercice de modélisation, et oui, j’ai exclu 2020. Je ne tiens pas à faire de comparaison avec la Covid-19, à qui 21,000 morts sont attribués, en France, ce soir. Si l’ordre de grandeur est le même, qu’on en me fasse pas dire qu’on est dans une situation comparable, et que la Covid-19 est finalement une simple grippe. Je ne me prononcerais pas sur ce point car je ne suis pas virologue. Mais surtout, on a 21,000 morts malgré des mesures incroyables de confinement. Je reviendrais probablement dans les jours qui viennent sur l’impact du contrôle par confinement, mais l’objet du court billet de ce soir était de voir comment on pouvait rapidement quantifier une surmortalité.

  • petite note technique : pour les 3 graphiques (mais pas pour les codes et les calculs) j’avais initialement normalisé les données. En effet, la population française a relativement pas mal changé en 20 ans, augmentant de 10%… Aussi, 1000 morts en 2003, ça n’est pas autant (relativement à la taille de la population) qu’en 2020.

Modeling Pandemics (3)

In Statistical Inference in a Stochastic Epidemic SEIR Model with Control Intervention, a more complex model than the one we’ve seen yesterday was considered (and is called the SEIR model). Consider a population of size N, and assume that S is the number of susceptible, E the number of exposed, I the number of infectious, and R for the number recovered (or immune) individuals, \displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-\beta {\frac {I}{N}}S\\[8pt]{\frac {dE}{dt}}&=\beta {\frac {I}{N}}S-aE\\[8pt]{\frac {dI}{dt}}&=aE-b I\\[8pt]{\frac {dR}{dt}}&=b I\end{aligned}}Between S and I, the transition rate is \beta I, where \beta is the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious subject. Between I and R, the transition rate is b (simply the rate of recovered or dead, that is, number of recovered or dead during a period of time divided by the total number of infected on that same period of time). And finally, the incubation period is a random variable with exponential distribution with parameter a, so that the average incubation period is a^{-1}.

Probably more interesting, Understanding the dynamics of ebola epidemics suggested a more complex model, with susceptible people S, exposed E, Infectious, but either in community I, or in hospitals H, some people who died F and finally those who either recover or are buried and therefore are no longer susceptible R.

Thus, the following dynamic model is considered\displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}\\[8pt]\frac {dE}{dt}&=(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}-\alpha E\\[8pt]\frac {dI}{dt}&=\alpha E+\theta\gamma_H I-(1-\theta)(1-\delta)\gamma_RI-(1-\theta)\delta\gamma_FI\\[8pt]\frac {dH}{dt}&=\theta\gamma_HI-\delta\lambda_FH-(1-\delta)\lambda_RH\\[8pt]\frac {dF}{dt}&=(1-\theta)(1-\delta)\gamma_RI+\delta\lambda_FH-\nu F\\[8pt]\frac {dR}{dt}&=(1-\theta)(1-\delta)\gamma_RI+(1-\delta)\lambda_FH+\nu F\end{aligned}}In that model, parameters are \alpha^{-1} is the (average) incubation period (7 days), \gamma_H^{-1} the onset to hospitalization (5 days), \gamma_F^{-1} the onset to death (9 days), \gamma_R^{-1} the onset to “recovery” (10 days), \lambda_F^{-1} the hospitalisation to death (4 days) while \lambda_R^{-1} is the hospitalisation to recovery (5 days), \eta^{-1} is the death to burial (2 days). Here, numbers are from Understanding the dynamics of ebola epidemics (in the context of ebola). The other parameters are \beta_I the transmission rate in community (0.588), \beta_H the transmission rate in hospital (0.794) and \beta_F the transmission rate at funeral (7.653). Thus

epsilon = 0.001 
Z = c(S = 1-epsilon, E = epsilon, I=0,H=0,F=0,R=0)
p=c(alpha=1/7*7, theta=0.81, delta=0.81, betai=0.588,
    betah=0.794, blambdaf=7.653,N=1, gammah=1/5*7,
    gammaf=1/9.6*7, gammar=1/10*7, lambdaf=1/4.6*7,
    lambdar=1/5*7, nu=1/2*7)

If \boldsymbol{Z}=(S,E,I,H,F,R), if we write \frac{\partial \boldsymbol{Z}}{\partial t} = SEIHFR(\boldsymbol{Z})where SEIHFR is

SEIHFR = function(t,Z,p){
  S=Z[1]; E=Z[2]; I=Z[3]; H=Z[4]; F=Z[5]; R=Z[6]
  alpha=p["alpha"]; theta=p["theta"]; delta=p["delta"]
  betai=p["betai"]; betah=p["betah"]; gammah=p["gammah"]
  gammaf=p["gammaf"]; gammar=p["gammar"]; lambdaf=p["lambdaf"]
  lambdar=p["lambdar"]; nu=p["nu"]; blambdaf=p["blambdaf"]
  N=S+E+I+H+F+R
  dS=-(betai*I+betah*H+blambdaf*F)*S/N
  dE=(betai*I+betah*H+blambdaf*F)*S/N-alpha*E
  dI=alpha*E-theta*gammah*I-(1-theta)*(1-delta)*gammar*I-(1-theta)*delta*gammaf*I
  dH=theta*gammah*I-delta*lambdaf*H-(1-delta)*lambdaf*H
  dF=(1-theta)*(1-delta)*gammar*I+delta*lambdaf*H-nu*F
  dR=(1-theta)*(1-delta)*gammar*I+(1-delta)*lambdar*H+nu*F
  dZ=c(dS,dE,dI,dH,dF,dR)
  list(dZ)}

We can solve it, or at least study the dynamics from some starting values

library(deSolve)
times = seq(0, 50, by = .1)
resol = ode(y=Z, times=times, func=SEIHFR, parms=p)

For instance, the proportion of people infected is the following

plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="",col="red")
lines(resol[,"time"],resol[,"H"],col="blue")

Modeling pandemics (2)

When introducing the SIR model, in our initial post, we got an ordinary differential equation, but we did not really discuss stability, and periodicity. It has to do with the Jacobian matrix of the system. But first of all, we had three equations for three function, but actually\displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}so it means that our problem is here simply in dimension 2. Hence\displaystyle {\begin{aligned}&X={\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&Y={\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\mu+\gamma)I\end{aligned}}and therefore, the Jacobian of the system is\begin{pmatrix}\displaystyle{\frac{\partial X}{\partial S}}&\displaystyle{\frac{\partial X}{\partial I}}\\[9pt]\displaystyle{\frac{\partial Y}{\partial S}}&\displaystyle{\frac{\partial Y}{\partial I}}\end{pmatrix}=\begin{pmatrix}\displaystyle{-\mu-\beta\frac{I}{N}}&\displaystyle{-\beta\frac{S}{N}}\\[9pt]\displaystyle{\beta\frac{I}{N}}&\displaystyle{\beta\frac{S}{N}-(\mu+\gamma)}\end{pmatrix}We should evaluate the Jacobian at the equilibrium, i.e. S^\star=\frac{\gamma+\mu}{\beta}=\frac{1}{R_0}andI^\star=\frac{\mu(R_0-1)}{\beta}We should then look at eigenvalues of the matrix.

Our very last example was

times = seq(0, 100, by=.1)
p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")

We can compute values at the equilibrium

mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"]
N=1
S = (gamma + mu)/beta
I = mu * (beta/(gamma + mu) - 1)/beta

and the Jacobian matrix

J=matrix(c(-(mu + beta * I/N),-(beta * S/N),
         beta * I/N,beta * S/N - (mu + gamma)),2,2,byrow = TRUE)

Now, if we look at the eigenvalues,

eigen(J)$values
[1] -0.024975+0.6318831i -0.024975-0.6318831i

or more precisely 2\pi/b where a\pm ib are the conjuguate eigenvalues

2 * pi/(Im(eigen(J)$values[1]))
[1] 9.943588

we have a damping period of 10 time lengths (10 days, or 10 weeks), which is more or less what we’ve seen above,

The graph above was obtained using

p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[1:1e5,"time"],resol[1:1e5,"I"],type="l",xlab="time",ylab="",lwd=3,col="red")
yi=resol[,"I"]
dyi=diff(yi)
i=which((dyi[2:length(dyi)]*dyi[1:(length(dyi)-1)])<0)
t=resol[i,"time"]
arrows(t[2],.008,t[4],.008,length=.1,code=3)

If we look carefully. at the begining, the duration is (much) longer than 10 (about 13)… but it does converge towards 9.94

plot(diff(t[seq(2,40,by=2)]),type="b")
abline(h=2 * pi/(Im(eigen(J)$values[1]))

So here, theoretically, every 10 weeks (assuming that our time length is a week), we should observe an outbreak, smaller than the previous one. In practice, initially it is every 13 or 12 weeks, but the time to wait between outbreaks decreases (until it reaches 10 weeks).

Modeling pandemics (1)

The most popular model to model epidemics is the so-called SIR model – or Kermack-McKendrick. Consider a population of size N, and assume that S is the number of susceptible, I the number of infectious, and R for the number recovered (or immune) individuals, \displaystyle {\begin{aligned}&{\frac {dS}{dt}}=-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-\gamma I,\\[6pt]&{\frac {dR}{dt}}=\gamma I,\end{aligned}}so that \displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}which implies that S+I+R=N. In order to be more realistic, consider some (constant) birth rate \mu, so that the model becomes\displaystyle {\begin{aligned}&{\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\gamma+\mu) I,\\[6pt]&{\frac {dR}{dt}}=\gamma I-\mu R,\end{aligned}}Note, in this model, that people get sick (infected) but they do not die, they recover. So here, we can model chickenpox, for instance, not SARS.

The dynamics of the infectious class depends on the following ratio:\displaystyle{R_{0}={\frac {\beta }{\gamma +\mu}}} which is the so-called basic reproduction number (or reproductive ratio). The effective reproductive ratio is R_0S/N, and the turnover of the epidemic happens exactly when R_0S/N=1, or when the fraction of remaining susceptibles is R_0^{-1}. As shown in Directly transmitted infectious diseases:Control by vaccination, if S/N<R_0^{-1} the disease (the number of people infected) will start to decrease.

Want to see it  ? Start with

mu = 0
beta = 2
gamma = 1/2

for the parameters. Here,  R_0=4. We also need starting values

epsilon = .001
N = 1
S = 1-epsilon
I = epsilon
R = 0

Then use the ordinary differential equation solver, in R. The idea is to say that \boldsymbol{Z}=(S,I,R) and we have the gradient \frac{\partial \boldsymbol{Z}}{\partial t} = SIR(\boldsymbol{Z})where SIR is function of the various parameters. Hence, set

p = c(mu = 0, N = 1, beta = 2, gamma = 1/2)
start_SIR = c(S = 1-epsilon, I = epsilon, R = 0)

The we must define the time, and the function that returns the gradient,

times = seq(0, 10, by = .1)
SIR = function(t,Z,p){
S=Z[1]; I=Z[2]; R=Z[3]; N=S+I+R
mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"]
dS=mu*(N-S)-beta*S*I/N
dI=beta*S*I/N-(mu+gamma)*I
dR=gamma*I-mu*R
dZ=c(dS,dI,dR)
return(list(dZ))}

To solve this problem use

library(deSolve)
resol = ode(y=start_SIR, times=times, func=SIR, parms=p)

We can visualize the dynamics below

par(mfrow=c(1,2))
t=resol[,"time"]
plot(t,resol[,"S"],type="l",xlab="time",ylab="")
lines(t,resol[,"I"],col="red")
lines(t,resol[,"R"],col="blue")
plot(t,t*0+1,type="l",xlab="time",ylab="",ylim=0:1)
polygon(c(t,rev(t)),c(resol[,"R"],rep(0,nrow(resol))),col="blue")
polygon(c(t,rev(t)),c(resol[,"R"]+resol[,"I"],rev(resol[,"R"])),col="red")

We can actually also visualize the effective reproductive number is R_0S/N, where

R0=p["beta"]/(p["gamma"]+p["mu"])

The effective reproductive number is on the left, and as we mentioned above, when we reach 1, we actually reach the maximum of the infected,

plot(t,resol[,"S"]*R0,type="l",xlab="time",ylab="")
abline(h=1,lty=2,col="red")
abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen")
points(max(t[resol[,"S"]*R0&gt;=1]),1,pch=19)
plot(t,resol[,"S"],type="l",xlab="time",ylab="",col="grey")
lines(t,resol[,"I"],col="red",lwd=3)
lines(t,resol[,"R"],col="light blue")
abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen")
points(max(t[resol[,"S"]*R0&gt;=1]),max(resol[,"I"]),pch=19)

And when adding a \mu parameter, we can obtain some interesting dynamics on the number of infected,

times = seq(0, 100, by=.1)
p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")

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&lt;=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…

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(&gt;|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(&gt;|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(&gt;|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(&gt;|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 !

Le R² pour justifier la causalite…

Ce soir, Louis (@LouisdeCharson) me signalait un article assez délirant où, dans une entrevue, un monsieur souvent présenté comme un chercheur (je n’ai pas réussi à trouver où) disait une phrase assez étonnante,

“… une variable ancienne explique 85 % de la variation d’une variable nouvelle. On peut donc parler de causalité”

Histoire de plaisanter un peu, je l’ai pris au mot, et j’ai ressorti un vieux jeu de données que j’avais utilisé dans un précédant billet, avec le traffic cycliste à Helsinki, en fonction de la température. Pour rester dans la classe des modèles linéaires, j’ai enlevé quelques jours d’hiver, et j’ai tenté de régresser la température le jour j+1 (la nouvelle variable comme dit dans la discussion) sur le nombre de cyclistes la veille, autrement dit le jour j

Si on regarde la régression, on obtient

lm(formula = temp ~ cyclists, data = df0)
Coefficients:
Estimate Std. Error t value Pr(&gt;|t|)
(Intercept) -4.170e+00  4.052e-01  -10.29   &lt;2e-16 ***
cyclists     1.066e-03  3.558e-05   29.96   &lt;2e-16 ***
 
Residual standard error: 2.833 on 212 degrees of freedom
Multiple R-squared: 0.809, Adjusted R-squared: 0.8081

autrement dit, 81\% de la variation de la température (moyenne) le  jour j+1 est expliqué par le nombre de cyclistes le jour j, c’est la définition du R^2 (“the coefficient of determination, denoted R2, is the proportion of the variance in the dependent variable that is predictable from the independent variable“). Et le monsieur en déduit qu’il y a une relation causale, autrement dit, dans notre exemple, le nombre de cyclistes sur la route le jour j cause la température le jour j+1. Loin de moi l’idée de me poser comme un expert sur les modèles causaux et les modèles climatiques, je doute que ce soit le cas… sinon la solution pour stopper le réchauffement climatique est toute trouvée : il faut diminuer le nombre de cyclistes !

 

Des modeles prédictifs en assurance

Cet post est aussi en ligne sur https://hal.archives-ouvertes.fr/hal-02350006, il a été coécrit avec Laurence Barry et Ewen Gallic.

Les compagnies d’assurance émettent des contrats qui prévoient des paiements d’indemnités en cas de survenance d’évènements aléatoires (accident, maladie, décès, etc.). En contrepartie, l’assuré doit s’acquitter d’une prime, dont le montant est déterminé ex-ante, avant le début de la période de couverture. Cette prime se décompose en deux termes : une prime pure (destinée à couvrir les pertes anticipées) et un chargement (incluant des commissions à des agents, divers frais, mais aussi couvrant contre le risque de variabilité des pertes). La prime pure est souvent calculée par classe de risque, et une classification est alors nécessaire.

Assurer une population hétérogène, ou l’importance de la classification

Le regroupement des risques selon diverses informations telles l’âge de l’assuré, son état de santé ou encore sa profession constitue ce que l’on appelle la classification des risques. Cette pratique de segmentation se justifie (à des fins d’admissibilité mais aussi de tarification) par la supposition que les risques sont placés dans des groupes relativement homogènes, au sein desquels les probabilités de survenance sont similaires. Pour Schauer (2006), cette « généralisation », qui vise à voir l’individu sous le prisme de sa classe de risque, de généraliser son comportement à partir de quelques variables explicatives, est probablement la raison d’être de l’actuaire : « to be an actuary is to be a specialist  in  generalization,  and  actuaries  engage  in  a form of decision-making that is sometimes called actuarial ». Statistiquement on cherche une méthode de classification aussi « discriminatoire » que possible[1], en gardant en mémoire que la discrimination est interdite, ce qui rend l’exercice périlleux, et souvent critiqué (nous y reviendrons plus loin).

Les assureurs évoquent souvent deux arguments pour justifier une segmentation. Le premier est qu’elle serait rendue économiquement nécessaire par la concurrence ; ne pas classifier conduit à une anti-sélection, les risques importants restant seuls chez les assureurs qui ne segmentent pas. Dans une telle situation, l’équilibre de marché ne serait pas possible puisque les risques faibles seraient chez un concurrent ayant segmenté. Si le facteur de risque était observable, tant par les assurés que les assureurs, il y aurait un phénomène d’auto-sélection, les assurés à risque faible ayant les polices les moins chères. Cette situation constitue un équilibre séparant de Nash. Mais si le facteur de risque n’est pas observable, un équilibre sous-optimal peut être atteint, résultant d’une externalité négative de cette information non-accessible, à la manière de Wilson (1977), tel que décrit dans Cummins et al. (1982) dans le cas des contrats d’assurance-vie. Cela dit, Kleindorfer & Kunreuther (1980) montrent qu’accéder à davantage d’information ne conduit pas nécessairement à une amélioration du bien-être des consommateurs. De plus si la classification n’est pas autorisée, l’équilibre est maintenu, les risques faibles subventionnant les risques élevés.

Le second argument avancé pour justifier une segmentation est que cette dernière (et le fait, par conséquent, d’ajuster les primes au risque) serait juste et équitable. Mais cette vision de l’équité n’a pas toujours été de mise et semble portée par les développements techniques. Ainsi la classification est devenue de plus en plus fine, multipliant les classes de risque et conduisant à des tarifs « personnalisés ». En plus des avancées statistiques, des facteurs économiques pourraient justifier cette sophistication : la concurrence de plus en plus forte sur certaines branches.

Incertitude en assurance

Il y a plusieurs manières de caractériser l’incertitude en assurance. Comme bien souvent quand on fait des prévisions, il convient de distinguer l’incertitude associée à l’estimation des probabilités et l’incertitude réelle sur le résultat (aléa de l’évènement). Pour la seconde notion, Hacking (1975) parle de probabilité structurelle, et c’est celle qui est souvent utilisée pour introduire les concepts de probabilité, par exemple avec des dés ou des jeux de cartes : les probabilités sont connues, seule l’issue du jeu est incertaine. Par exemple je sais que la probabilité d’avoir 6 en lançant un dé est 1/6 (compte tenu de la géométrie du cube).

D’un point de vue statistique, la probabilité se mesure quand on peut observer une fréquence, c’est-à-dire une répétition de risques semblables. Les statisticiens ont ainsi défini une notion de probabilité empirique, basée sur la répétition[2]. Si, en lançant mille dés j’obtiens 173 fois la face 6, la probabilité empirique d’avoir 6 est de 17,3%. La loi des grands nombres nous assure que cette fréquence va tendre vers la vraie valeur en répétant l’expérience, et le théorème central limite permet d’en contrôler les fluctuations. C’est la première incertitude dont nous parlions au début de cette section, que nous appellerions l’erreur d’estimation.

On peut enfin mentionner deux notions supplémentaires ; tout d’abord, les probabilités conditionnelles. Cette idée est introduite en assurance par de Moivre, ou de Witt, lorsqu’ils notaient que pour estimer une probabilité de décès, il fallait considérer des personnes de même âge. C’est cette idée que l’on retrouve quand on considère une classification : on veut des risques homogènes, similaires, sans être pour autant identiques. La probabilité que l’on obtient est alors conditionnelle à ce facteur commun qui caractérise la classe observée. Dans notre exemple des dés, cela revient à dire qu’il ne faut pas lancer mille dés, mais mille fois le même dé – ou à défaut des dés semblables.

Enfin, les probabilités subjectives ont été formalisées par Bruno de Finetti et Leonard Savage (ainsi que plus philosophiquement par Frank P. Ramsey) pour comprendre et modéliser la prise de décision. Elles sont relativement populaires en économie de l’incertain, mais difficile à mettre en œuvre dans un contexte de valorisation de contrats d’assurance automobile ou habitation. Il s’agit d’un jugement, qui ne peut être confronté à la réalité, mais envisageable pour l’assurance de risques encore mal connus (McGrayne (2012) évoque ainsi les premiers contrats d’assurance aviation). Une approche bayésienne consiste alors à combiner cette probabilité subjective avec la probabilité comme fréquence observée d’un phénomène : partant d’une croyance a priori, on affine l’estimation par une mise à jour progressive en répétant les expériences. Classiquement, la probabilité d’avoir la face 6 sera une moyenne entre notre croyance (1 chance sur 6) et une probabilité dite historique, obtenue en faisant quelques lancés (3 sur 20 lancers, par exemple). Les poids attribués aux deux dépendant du nombre d’expériences effectuées : on donnera plus de crédit à l’expérience si on fait mille lancés que si on en fait soixante.

Incertitude sur le résultat, ou aléa fondamental

Les probabilités prédictives, utilisées pour calculer la prime d’un contrat d’assurance, sont la première étape d’un problème de classification. Un outil classique pour juger de la pertinence d’un classifieur est la courbe ROC, décrite dans Kuhn (2018)) : on compare la probabilité individuelle (a priori, telle que résultant du modèle de classification) à un seuil, compris entre 0 et 1; si la probabilité est inférieure au seuil, l’estimation est que la personne survit, sinon qu’elle décède.

On compare ensuite cette estimation aux réalisations (ex-post) de survie et de décès. Pour chaque seuil, on peut considérer la matrice classique dite matrice de confusion de théorie de la décision : elle consiste à répartir les observations suivant le résultat observé (en colonne) et l’estimation résultant du modèle en ligne (en fonction de la probabilité estimée pour l’individu et le seuil que l’on s’est fixé). On peut ainsi partager la population entre les classements corrects, et les erreurs (dont les « faux positifs » si la personne a survécu malgré une probabilité estimée de décès supérieure au seuil, et les « faux négatifs » si la personne décède malgré une probabilité estimée inférieure au seuil).

Figure 1: Courbe ROC et classification pour un seuil de probabilité valant 1.5%.

La courbe ROC est obtenue en faisant varier le seuil. Chaque seuil correspond à un point de la courbe, rapportant graphiquement les taux de faux positifs (en abscisse) et de vrais positifs (en ordonnée), comme sur la Figure 1.

Considérons un groupe de 1000 assurés, où 20 personnes sont décédées l’an passé. Supposons un modèle dans lequel on admet que la population est parfaitement homogène, la probabilité estimée de décès est de 2% pour tout le monde. Dans ce cas pour tout seuil supérieur à 2%, on estimera que la totalité de la population survit : on aura un taux de faux positifs de 0% et un taux de vrais positifs de 0%, d’où un point (0,0) sur le graphe. A l’inverse pour tout seuil inférieur à 2%, on estimera que la totalité de la population décède : on aura un taux de faux positifs de 100% et un taux de vrais positifs de 100%, d’où un point (1,1) sur le graphe. La courbe de ROC de ce modèle uniforme à 2% est donc la diagonale du carré sur la figure 1.

Mais on peut aussi imaginer qu’il existe un peu d’hétérogénéité avec, par exemple, une probabilité de décès de 1% pour une moitié de la population et de 3% pour l’autre moitié, ou encore que le modèle produit des probabilités comprises entre 1% et 3% de façon non dichotomique. Les données simulées pour construire la courbe noire sur la Figure 1 suppose que la population a des probabilités de décès variables, comprises entre 1% et 3%, obtenues par une régression logistique. Comme le montre le tableau de droite, on commet des erreurs, et comme le montre celle de gauche, la nature de celle-ci varie en fonction du seuil choisi, qui modifie les taux de faux positifs et de faux négatifs.

Le cas extrême serait celui où le modèle aurait correctement attribué une probabilité de 100% aux 20 personnes qui sont effectivement décédées. C’est la courbe rouge sur la Figure 1. Ce partage est possible ex-post, une fois réalisation de l’aléa : a posteriori, il y a une certitude de décès pour ceux qui sont effectivement morts. Mais cela n’a cependant pas grande réalité dans l’assurance, à moins d’imaginer que l’actuaire serait un oracle, qui saurait avec certitude qui va mourir, et qui va survivre. La réalité est plutôt celle de la situation intermédiaire entre la courbe rouge et la diagonale, avant d’arriver dans la région hachurée, où le taux d’erreur est faible, mais pas nul : on ne peut pas prédire, avec certitude, qui va décéder. L’assurance n’est possible que si cette borne supérieure n’est pas trop élevée.

Incertitude statistique, données et modèles

Une question fondamentale pour la survie de l’assurance est de savoir où se situe cette borne supérieure : jusqu’où peut-on aller, entre les deux cas extrêmes (population homogène avec une probabilité de 2% pour tous, et une population très discriminée, avec 2% de la population ayant 100% de chances de mourir, et l’autre 0%) ? Et de quoi cette borne dépend-elle ? En particulier, des modèles plus complexes, tels que les réseaux de neurones très profonds permettent-ils vraiment d’améliorer la prévision ? Et l’enrichissement de données, tel qu’on l’observe grâce aux objets connectés et la fusion avec toutes sortes d’informations externes, va-t-il déplacer la borne supérieure vers le haut ?

Si l’apprentissage profond – voir Goodfellow et al. (2018) – permet d’avoir des classifieurs d’images avec un taux d’erreur proche de 0%, il est difficile d’imaginer qu’il sera possible de prévoir, presque un an à l’avance (à la signature du contrat), qui décèdera dans l’année, qui aura la grippe, qui aura un dégât des eaux, etc. Les modèles plus complexes permettent d’améliorer les prévisions, en tenant compte de non-linéarités, d’effets croisés entre les variables tarifaires, mais pas au point de faire disparaître l’aléa. Et tant que l’assurance est envisagée ex-ante (la prime est fixée au début de la période de couverture), il est difficile d’imaginer que rajouter de l’information fera aussi disparaître l’aléa. C’est d’ailleurs le cas pour les tests génétiques qui n’expliquent qu’une (petite) partie du risque de cancer, par exemple. Et rajouter des données revient souvent à rajouter du bruit, ce qui rend le travail d’analyse plus complexe. Cependant, force est de constater que des modèles plus complexes et des données plus riches tendent effectivement à « améliorer » la prévision, en remontant la courbe ROC vers le haut. Mais se pose-t-on les bonnes questions ? Que signifie vraiment une borne très éloignée du cas homogène, sur la diagonale ?

Homogénéité, équité et causalité

Comme nous l’avons vu, la tarification en assurance repose sur une répartition des risques (des contrats) en catégories, au sein desquelles la distribution des pertes peut être estimée, afin de fixer un niveau de primes. La répartition se fait à partir des caractéristiques de l’assuré, et du bien assuré. En retraçant l’histoire de l’assurance, Ewald (1986) montre que les mécanismes de prévoyance se sont mis en place en déplaçant la charge des accidents du travail sur la société : on abandonne l’idée d’une responsabilité individuelle de l’accident en faveur de la solidarité. L’assurance distingue « entre le dommage que subit tel ou tel individu — c’est affaire de chance ou de malchance — et la perte liée au dommage dont l’attribution est, quant à elle, toujours collective et sociale ». Ce principe de solidarité sociale, de mutualisation des risques, fait que le risque (en assurance) est toujours pensé collectivement.

Aujourd’hui, les tarifs sont considérés comme « justes », ou « actuariellement équitables » si chaque prime correspond à la perte attendue (pour ne pas dire « espérée », au sens mathématique) pour chaque assuré. Dans cette perception de l’équité, une hypothèse essentielle est que les classes soient « homogènes ». En effet, dans l’hypothèse inverse, les personnes les moins risquées subventionnent les personnes les plus risquées, ce qui est perçu comme socialement injuste.

On peut décrire cette version de l’équité actuarielle à l’aide de la formule de décomposition de la variance. La variance globale se décompose en effet en deux termes, la variance inter-classes et la variance intra-classes : l’ « équité actuarielle » vise à ce que  les classes de risque soient relativement distinctes les unes des autres, donc une variance inter-classes forte, accompagnée d’une homogénéité des classes, donc une variance intra-classes faible. D’un point de vue statistique, chercher à augmenter l’une est équivalente à faire diminuer l’autre. Cette mécanique n’est pas toujours claire pour des observateurs non avertis ; ainsi dans l’affaire Manhart, un des cas les plus documentés sur la discrimination par le genre en assurance, le juge Stevens affirme : « we focus on fairness to individuals rather than on fairness to classes […] even a true generalization about a class is an insufficient reason for disqualifying an individual to whom the generalization does not apply» (cité dans Anzalone (2016)). Autrement dit, pour la justice, un critère statistique de type « true generalization » ne peut être opposable à un individu.

Une autre critique importante, que l’on retrouve dans la « gender directive », est le lien entre discrimination et causalité. En effet, statistiquement, les actuaires vont chercher des facteurs de classification fortement corrélés avec la sinistralité. Mais il est possible que ces facteurs ne soient qu’un proxy de la vraie variable causale, restée elle inobservée, conduisant à une mauvaise estimation du risque pour certains. Comme le notent Antonio et Charpentier (2017), le genre a ainsi été longtemps utilisé en assurance automobile car très corrélé avec des variables associées au style de conduite et à d’autres variables historiquement non observables (mais qui le sont aujourd’hui grâce aux objets connectés, comme le kilométrage, les heures de conduites, les types de routes utilisés, etc).

Ce lien avec les mécanismes causaux est d’ailleurs relativement profond, et Hacking (1975) y voit une connexion avec la « révolution probabiliste » : on peut assez facilement mettre en évidence des corrélations, mais les causes, si elles existent, nous restent plus opaques. Laplace au début du 19e siècle déclare ainsi que « la probabilité est relative en partie à nos connaissances, en partie à notre ignorance », liant les probabilités à la fois à une vision newtonienne déterministe du monde et à notre incapacité à le connaitre parfaitement. Cette dernière composante fait que l’on ne peut pas annoncer la date exacte du décès d’un individu, mais statistiquement, dans un groupe homogène, on peut prédire le nombre de décès au cours d’une année. Et pour revenir à la relation causale, le tabagisme par exemple ne cause pas forcément une mort prématurée mais fumer sera vu comme dangereux car il augmente la probabilité de décès pendant une période donnée. Ainsi nous montre Hacking (1975), la causalité est pensée aujourd’hui dans un contexte probabiliste, et non plus déterministe.

Abraham, K.S. (1985) Efficiency and fairness in insurance risk classification, Virginia Law Review 71: 403-451.

Anderson, A.W. (1978). A Critique of the Manhart Brief. The Actuary, 12:5.

Antonio, K. & Charpentier, A. (2017). La tarification par genre en assurance, corrélation ou causalité ? Risques, 109.

Anzalone, C.A. (2016). U.S. Supreme Court Cases on Gender and Sexual Equality. Routledge.

Bailey, H., Hutchison, T. & Narber, G. (1975) The regulatory challenge to life insurance classification, Drake Law Review Insurance Law Annual 4: 779-827

Barry L. (2019). Justice ou justesse? L’équité de l’assurance. Working paper #15, chaire PARI.

Charpentier, A. & Denuit, M. (2004). Mathématiques de l’Assurance Non-Vie : Principes Généraux de Théorie du Risque. Economica.

Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press.

Cummins, J.D., Smith, B.D., Vance, R.N. & VanDerhai, J.L. (1982). Risk Classification in Life Insurance. Kluwer-Nijhoff Publishing.

Ewald F. (1986). L’État providence. Grasset.

Fisher, R. A. (1936).  The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics. 7 (2): 179–188.

Frézal, S. & Barry, L. (2019). Fairness in Uncertainty: Some Limits and Misinterpretations of Actuarial Fairness, Journal of Business Ethics.

Goodfellow, U., Bengio, Y. & Courville, A. (2018) L’apprentissage profond. Massot Edition.

Hacking, I. (1975) The Emergence of Probability. Cambridge University Press

Hoy, M. (1982) Categorizing risks in the insurance industry, The Quarterly Journal of Economics 97: 321-336.

Kleindorfer, P. & Kunreuther, H. (1980) Misinformation and Equilibrium in Insurance Markets, in Economic Analysis of Regulated Markets, Jörg Finsinger Editor, Springer Verlag, 67-90

Kuhn, M. & Johnson, K. (2018). Applied Predictive Modeling. Springer Verlag.

Martin,G.D. (1976) Gender Discrimination in Pension Plans, Journal of Risk and Insurance 43.

McGrayne, S.B. (2012) The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, and Emerged Triumphant from Two Centuries of Controversy. Yale University Press.

Ramsey P.F. (1926). Truth and Probability.

Schauer, F. (2006) Profiles, Probabilities, and Stereotypes. Harvard University Press.

Von Mises, R. (1957). Probability, Statistics and Truth. Dover publications.

Wilson, C. (1977). A model of insurance markets with incomplete information. Journal of Economic Theory, 16:2, 167-207.

[1] Au sens statistique du mot, dans le sens introduit par Fisher (1936).

[2] Dans cette approche fréquentiste, et notamment pour Ronald Fisher et Richard von Mises, la probabilité d’un évènement unique (dit « one shot ») n’a pas de sens.

Principal Component Analysis: A Generalized Gini Approach

With Stéphane Mussard and Téa Ouraga, we recently uploaded on arxiv a paper Principal Component Analysis: A Generalized Gini Approach,

A principal component analysis based on the generalized Gini correlation index is provided. It is proven that the reduction dimensionality based on the generalized Gini correlation index, that relies on city-block distances, is robust to outliers.

Some codes are also available on a dedicated github repo.