Category Archives: Statistics

From multinomial regression to binary classification on some Siamese data

There are two kinds of people in the world: people who think there are two kinds of people in the world and people who don’t

(borrowed from Menand (2018)). Because things are always simpler when we face only binary choice, aren’t they? But consider here the case were multiple options are possible, and let us see if we cannot get back to simpler binary choices. Consider a collection of observations (y_i,\boldsymbol{x}_i) where y_i is some categorical variable, y_i\in\mathcal{A} where \mathcal{A}=\lbrace A_1,\cdots,A_\kappa \rbrace, with \kappa possible categories. Let \mathcal{I}_k=\lbrace i:y_i\in A_k \rbrace.

In a classical multinomial logistic regression, suppose that A_1 is the reference, then \mathbb{P}[Y=A_j|\boldsymbol{X}=\boldsymbol{x}]=\frac{\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_j]}{1+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_2]+\cdots+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]}With a lot a categories, and a small number of observations, inference can be complicated, and non-robust.

  • the Siamese dataset

The name Siamese I use, here, comes from Siamese Networks. Or sort of… As we say in French, it is an « histoire de l’homme qui a vu l’homme qui a vu l’ours » (story of the man who saw the man who saw the bear). A few years ago, a student tried to explain to me the idea of Siamese Networks and this is what I understood. I might be completely wrong, but the idea I got from it did make sense, in my mind at least. That is the story of that blog post…

The idea of the siamese algorithm will be to consider all pairs of observations, (y_i,\boldsymbol{x}_i) and (y_j,\boldsymbol{x}_j) :

  1. \tilde y_{i,j}=\boldsymbol{1}(y_i=y_j) indicating if individuals i and j are in the same category
  2. \tilde{\boldsymbol{x}}_{i,j} is a collection of p-1 variables,
  • \tilde {x}_{k:i,j}={x}_{k:i}-{x}_{k:j} if x_k is continuous, or \tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}| (we can use another metric, e.g. \tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|^2, and this is why I decided to use some GAM model in the logistic regression on the Siamese dataset)
  • \tilde {x}_{k:i,j}=({x}_{k:i},{x}_{k:j})\in\mathcal{X}_k\times\mathcal{X}_k if x_k is a categorical variables (taking values in the set \mathcal{X}_k), or \tilde {x}_{k:i,j}=\boldsymbol{1}({x}_{k:i}\neq{x}_{k:j})\in\{0,1\}

The original dataset was a n\times p matrix, and (if there are no categorical variable), it becomes a n(n-1)/2\times p matrix. The key point is that if the original variable y_i was multinomial, y_{i,j} is now binomial. For instance, if our initial dataset was the following, with two covariates, one continuous and one categorical

its siamese counterpart is the following

  • Classification step

On the dataset (\tilde y_{i,j},\tilde {x}_{k:i,j})_{i,j}, fit a logistic regression, \mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}]=\frac{\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}(or any classification model – CART, random forest, etc). But that is the easy part (unless n is large, because the siamese dataset has (roughly) n^2/2 rows). The difficult task is the prediction

  • Prediction step

Consider a new input variable \boldsymbol{x}_{\cdot}, and define its siamese version, \tilde{\boldsymbol{x}}_{\cdot}=(\tilde{\boldsymbol{x}}_{\cdot,j})_j, i.e. a database with n rows. Then compute
p_{\cdot,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{\cdot,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]} where p_{\cdot,j} is the probability that (y_j,\boldsymbol{x}_{j}) and (y_{\cdot},\boldsymbol{x}_{\cdot}) are in the same category, as well as p_{i,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{i,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}
Let \boldsymbol{p}_{\cdot}=(p_{\cdot,j}), and similarly \boldsymbol{p}_{i}=(p_{i,j}). Then several techniques can be used to predict y_{\cdot}.

  1. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{j=1,\cdots,n}{\text{argmax}}\{p_{\cdot,j}\}: the predicted class is the one of the observation the most likely to be other same class
  2. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{\ell=1,\cdots,k}{\text{argmax}}\{\overline{p}_{\ell}\}, where\overline{p}_{\ell} = \frac{1}{n_{\ell}}\sum_{j\in\mathcal{I}_s} \boldsymbol{1}(y_j =y_{\ell}),\text{ where }\mathcal{I}_s=\lbrace i:p_{\cdot,i}>s\rbraceconsider only probabilities sufficiently high to be considered, and predict the most important class (majority rule)
  3. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{\theta_i\} where \theta_i=\cos(\boldsymbol{p}_{\cdot},\boldsymbol{p}_{i})=\displaystyle{\frac{\boldsymbol{p}_{\cdot}\cdot\boldsymbol{p}_{i}}{\|\boldsymbol{p}_{\cdot}\|\|\boldsymbol{p}_{i}\|}}
  4. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{KL_{\cdot|i}\} and KL_{\cdot|i}=\displaystyle{\sum_{j=1}^n p_{\cdot,j}\log\frac{p_{\cdot,j}}{p_{i,j}}} (but one can select another metric)
  5. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{j\in\mathcal{J}}{\text{argmax}}\{p_{\cdot,j}\} and \mathcal{J} is a sample of k observations, chosen randomly, one in each group (one-shot procedure): the predicted class is the one of the observation the most likely to be other same class

Heuristically, it can be related to some k nearest neighbors strategy: we give the attribute that most neighbors have. The total distance is a weighted sum of the componentwise distances (for the logistic regression).

  • Simulation study

In order to test that technique, let us generate some multinomial model where y has 10 possible labels, with 6 (independent) covariates x_1,\cdots,x_6, and \mathbb{P}[Y=A_k|\boldsymbol{X}=\boldsymbol{x}]\propto \exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k] (where coefficients \boldsymbol{\beta}_k where generated randomly) for k\in\{1,2,\cdots,10\} (there were 10 categories) and with n=700 observations.

n=700
X1=rnorm(n)
X2=rnorm(n)
X3=rnorm(n)
X4=rnorm(n)
X5=rnorm(n)
X6=rnorm(n)
X=cbind(1,X1,X2,X3,sqrt(abs(X4)),X5*X1,X6)
k = 10
 PARAM = matrix(rnorm(k*6),k,6)
 PARAM[,1]=PARAM[,1]-1
 PARAM=cbind(PARAM,0)
 P=matrix(NA,n,k-1)
 for(j in 1:(k-1)) P[,j] = X %*% (PARAM[j,])+rnorm(n)
 P=cbind(P,0)
S=apply(exp(P),1,sum)
Pb = exp(P)/S
tirage = function(i){
      sample(1:10,size=1,prob = Pb[i,])
}
Y = LETTERS[Vectorize(tirage)(1:n)]
dbase = data.frame(Y=as.factor(Y),X1,X2,X3,X4,X5)

In the paragraph previously, I suggested to take the most likely one. Being wrong means that it was not the first choice. But perhaps being the second or the third choice is not that bad, actually. So in my simulations, I look at the proportion of predictions were our prediction is the good one (top 1), if the true one is either the most likely or the second most likely (top 2), or in the top 3. That will be on my x-axis. I draw some line, but we simply have three points (top 1, top 2 and top 3). I compute the proportion of good prediction, using cross-validation techniques (10-fold). The black lines are one of the methods described above. The red one is the standard multinomial model (with a logistic link function). For the Siamese model, I tried several models. I tried a logistic regression, and some smooth version (GAM) on top

and a classification tree, on the left, as well as some random forest on the right, below.

It looks like the multinomial approach performs always better than any Siamese one… and to be honest, I am disappointed.

Here is the code I did use when I considered a logistic regression on the Siamese dataset,

set.seed(1)
kfold = sample(rep(1:10,n/10))
 
KFOLDglm = function(i){
i_test=which(kfold==i)
i_calibration=which(kfold!=i)
y=credit[i_calibration,"Y"]
tirage = function(){
v=c(sample(i_calibration[y==levels(y)[1]],size=1),
    sample(i_calibration[y==levels(y)[2]],size=1),
    sample(i_calibration[y==levels(y)[3]],size=1),
    sample(i_calibration[y==levels(y)[4]],size=1),
    sample(i_calibration[y==levels(y)[5]],size=1),
    sample(i_calibration[y==levels(y)[6]],size=1),
    sample(i_calibration[y==levels(y)[7]],size=1),
    sample(i_calibration[y==levels(y)[8]],size=1),
    sample(i_calibration[y==levels(y)[9]],size=1),
    sample(i_calibration[y==levels(y)[10]],size=1))
names(v)=levels(y)
return(v)
}
 
LogisticModel <- multinom(Y ~ ., data = credit[i_calibration,], trace=FALSE)
 
comparaisonx = function(base,x=base[1,]){
  mix_base = base
  for(j in 1:ncol(base)){
    xj = as.numeric(x[j])-base[,j]
    mix_base[,j] = (xj)
  }
  mix_base
}
comparaisony = function(base,y=base[1]){
  as.factor(base == y)
}
creditx = credit[,-which(names(credit) == "Y")]
nc=length(i_calibration)
B=comparaisonx(base = creditx[i_calibration[2:nc],],x=creditx[i_calibration[1],])
B$Y=comparaisony(base = credit[i_calibration[2:nc],"Y"],y=credit[i_calibration[1],"Y"])
for(i in 2:(nc-1)){
  B0=comparaisonx(base = creditx[i_calibration[(i+1):nc],],x=creditx[i_calibration[i],])
  B0$Y=comparaisony(base = credit[i_calibration[(i+1):nc],"Y"],y=credit[i_calibration[i],"Y"])
  B=rbind(B,B0)
}
credit_mix = B
 
OneShotLogisticModel <- glm(Y ~ ., data = credit_mix, family=binomial)
A_ref = table(credit[i_calibration,"Y"])/length(i_calibration)
 
vect_oneshot = function(i){
  B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,])
  predict(OneShotLogisticModel,type="response",newdata=B2)
}
 
prediction_oneshot = function(i,type=1){
B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,])
p=predict(OneShotLogisticModel,type="response",newdata=B2)
y=credit[i_calibration,"Y"]
base = data.frame(p,y)
base = base[rev(order(base$p)),]
if(type==1){T = table(base$y[1:11])
return(names(which.max(T)))}
if(type==2){return(base$y[1])}
if(type==3){A=table(base$y[1:10])/10
T=A/A_ref
return(names(which.max(T)))}
if(type==4){
  costheta = rep(NA,length(i_calibration))
  for(j in 1:length(i_calibration)){
    vecteur_proba = vect_oneshot(i_calibration[j])
    costheta[j] = sum(vecteur_proba*p)/(sqrt(sum(vecteur_proba^2))*sqrt(sum(p^2)))
  }
  return(y[which.max(costheta)])}
if(type==5){
  kl = rep(NA,length(i_calibration))
  for(j in 1:length(i_calibration)){
    vecteur_proba = vect_oneshot(i_calibration[j])
    kl[j] = sum(p*log(vecteur_proba/p))
  }
  return(y[which.max(as.vector(kl))])}
if(type==6){ ## one shot : tirer au hasard un de chaque, et dire lequel est plus credible !
 
y=credit[i_calibration,"Y"]
tirage = function(){
    v=c(sample(i_calibration[y==levels(y)[1]],size=1),
        sample(i_calibration[y==levels(y)[2]],size=1),
        sample(i_calibration[y==levels(y)[3]],size=1),
        sample(i_calibration[y==levels(y)[4]],size=1),
        sample(i_calibration[y==levels(y)[5]],size=1),
        sample(i_calibration[y==levels(y)[6]],size=1),
        sample(i_calibration[y==levels(y)[7]],size=1),
        sample(i_calibration[y==levels(y)[8]],size=1),
        sample(i_calibration[y==levels(y)[9]],size=1),
        sample(i_calibration[y==levels(y)[10]],size=1))
    names(v)=levels(y)
    return(v)
  }
  pd=rep(NA,101)
for(ix in 1:101){
ids = tirage()
B2=comparaisonx(base = creditx[ids,],x=creditx[i,])
p=predict(OneShotLogisticModel,type="response",newdata=B2)
pd[ix]=levels(y)[which.max(p)]
}
levels(y)[which.max(table(pd))]
}
}
PRED0=as.character(credit[i_test,"Y"])
PRED1=as.character(predict(LogisticModel,type = "class",
                           newdata=credit[i_test,]))
PRED21=as.character(Vectorize(function(i) prediction_oneshot(i,type=1))
                    (i_test))
PRED22=as.character(Vectorize(function(i) prediction_oneshot(i,type=2))(i_test))
PRED23=as.character(Vectorize(function(i) prediction_oneshot(i,type=3))(i_test))
PRED24=as.character(Vectorize(function(i) prediction_oneshot(i,type=4))(i_test))
PRED25=as.character(Vectorize(function(i) prediction_oneshot(i,type=5))(i_test))
PRED26=as.character(Vectorize(function(i) prediction_oneshot(i,type=6))(i_test))
B=data.frame(PRED0,PRED1,PRED21,PRED22,PRED23,PRED24,PRED25,PRED26)
B
}
 
s=1/100;setTxtProgressBar(pb, s*2)
for(i in 2:10){
  PREDICTION = rbind(PREDICTION,KFOLDglm(i))
s=s+1/100;setTxtProgressBar(pb, s*2)}
for(j in 1:5) PREDICTION[,j]=as.character(PREDICTION[,j])
L=list()
v=mean(PREDICTION[,1]!=PREDICTION[,2])
names(v)="logistic"
L[["logistic"]]=v
v=c(mean(PREDICTION[,1]!=PREDICTION[,3]),
  mean(PREDICTION[,1]!=PREDICTION[,4]),
  mean(PREDICTION[,1]!=PREDICTION[,5]),
  mean(PREDICTION[,1]!=PREDICTION[,6]),
  mean(PREDICTION[,1]!=PREDICTION[,7]),
  mean(PREDICTION[,1]!=PREDICTION[,8]))
names(v)=c("top10","max","10norm","cos","KL","OS")
L[["glm"]]=v

Autocalibration for Insurance Pricing with Machine Learning

With Michel Denuit and Julien Trufin, we recently uploaded a joint paper on ArXiv, entitled Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning.

Boosting techniques and neural networks are particularly effective machine learning methods for insurance pricing.
Often in practice, there are nevertheless endless debates about the choice of the right loss function to be used to train the machine learning model, as well as about the appropriate metric to assess the performances of competing models. Also,
the sum of fitted values can depart from the observed totals to a large extent and this often confuses actuarial analysts.
The lack of balance inherent to training models by minimizing deviance outside the familiar GLM with canonical
link setting has been empirically documented in Wüthrich (2019, 2020) who attributes it to the early stopping rule in gradient descent methods for model fitting.
The present paper aims to further study this phenomenon when learning proceeds by minimizing Tweedie deviance.
It is shown that minimizing deviance involves a trade-off between the integral of weighted differences of lower partial moments and the bias measured on a specific scale.
Autocalibration is then proposed as a remedy. This new method to correct for bias adds an extra local GLM step to the analysis.
Theoretically, it is shown that it implements the autocalibration concept in pure premium calculation and ensures that balance also holds on a local scale, not only at portfolio level as with existing bias-correction techniques. The convex order appears to be the natural tool to compare competing models, putting a new light on the diagnostic graphs and associated metrics proposed by Denuit et al. (2019).

In this paper, we started with a simple observation : with GLM (like a Poisson regression) we have unbiased predictions, in the sense that \widehat{y}_1+...+\widehat{y}_n=y_1+...+y_n (on the training dataset, we should expect to have \widehat{y}_1+...+\widehat{y}_n\approx y_1+...+y_n on the validation dataset). But with any machin learning algorithm, this is no longer the case. For instance, with a boosting algorithm, we can end up with a significant bias. This is an application on an insurance dataset from the CASdataset package, with a GLM (Poisson regression) on the left, and additive smooth version, GAM in the middle and some boosting algorithm on the right. Here is the dispersion of the predictions, \widehat{y}_i‘s,

Let \widehat{\pi} denote a fited model. If \widehat{\pi} was close to the true parameter of the Poisson model, \mu where \mu(\boldsymbol{x})=\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}], then function  s\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=s] should be very close to the first diagonal (since \mathbb{E}[Y|\mu(\boldsymbol{X})=s]=s). And this is was we observe below, for the GLM. But the boosting algorithm has a significant bias.

Furthermore, the bias is not the same everywhere. We can see it more precisely on a quantile version of the x-axis

that is u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)]. Here is a multiplicative version of that bias u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)]/u

The idea of autocalibration is to use that multiplicative factor to correct from the bias. Thus, if \widehat{\pi}(\boldsymbol{x}) is close to the 75% upper quantile, then the true value should be 20% larger for the bootsting algorithm. Here is the new distribution of the predictions, with that correction

and we can observe that u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)] is now much closer to the first diagonal (below is the multiplicative bias)

We discuss in the paper this correction for local bias, so that u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)] becomes as close as possible to the first diagonal… see https://arxiv.org/abs/2103.03635 for a complete version of the paper, and https://github.com/freakonometrics/autocalibration/ for the complete version of the R codes

Principal Component Analysis: A Generalized Gini Approach

Our paper, with Stéphane Mussard and Téa Ouraga, entitle Principal Component Analysis: A Generalized Gini Approach is finally out in the European Journal of Operations Research.

A principal component analysis based on the generalized Gini correlation index is proposed (Gini PCA). The Gini PCA generalizes the standard PCA based on the variance. It is shown, in the Gaussian case, that the standard PCA is equivalent to the Gini PCA. It is also proven that the dimensionality reduction based on the generalized Gini correlation matrix, that relies on city-block distances, is robust to outliers. Monte Carlo simulations and an application on cars data (with outliers) show the robustness of the Gini PCA and provide different interpretations of the results compared with the variance PCA.

Some general thoughts on Partial Dependence Plots with correlated covariates

The partial dependence plot is a nice tool to analyse the impact of some explanatory variables when using nonlinear models, such as a random forest, or some gradient boosting.The idea (in dimension 2), given a model m(x_1,x_2) for \mathbb{E}[Y|X_1=x_1,X_2=x_2]. The partial dependence plot for variable x_1 is model m is function p_1 defined as x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)]. This can be approximated, using some dataset using \widehat{p}_1(x_1)=\frac{1}{n}\sum_{i=1}^n m(x_1,x_{2,i})My concern here what the interpretation of that plot when there are some (strongly) correlated covariates. Let us generate some dataset to start with

n=1000
library(mnormt)
r=.7
set.seed(1234)
X = rmnorm(n,mean = c(0,0),varcov = matrix(c(1,r,r,1),2,2))
Y = 1+X[,1]-2*X[,2]+rnorm(n)/2
df = data.frame(Y=Y,X1=X[,1],X2=X[,2])

As we can see, the true model is here is y_i=\beta_0+\beta_1 x_{1,i}+\beta_2x_{2,i}+\varepsilon_i where \beta_1 =1 but the two variables are positively correlated, and the second one has a strong negative impact. Note that here

reg = lm(Y~.,data=df)
summary(reg)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.01414    0.01601   63.35   <2e-16 ***
X1           1.02268    0.02305   44.37   <2e-16 ***
X2          -2.03248    0.02342  -86.80   <2e-16 ***

If we estimate a wrongly specified model y_i=b_0+b_1 x_{1,i}+\eta_i, we would get

reg1 = lm(Y~X1,data=df)
summary(reg1)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.03522    0.04680  22.121   <2e-16 ***
X1          -0.44148    0.04591  -9.616   <2e-16 ***

Thus, on the proper model, \widehat{\beta}_1\sim+1.02 while \widehat{b}_1\sim-0.44  on the mispecified model.

Now, let us look at the parial dependence plot of the good model, using standard R dedicated packages,

library(pdp) 
pdp::partial(reg, pred.var = "X1", plot = TRUE,
              plot.engine = "ggplot2")

which is the linear line y=1+x, that corresponds to y=\beta_0+\beta_1x.

library(DALEX)
plot(DALEX::single_variable(DALEX::explain(reg,
data=df),variable = "X1",type = "pdp"))

which corresponds to the previous graph. Here, it is also possible to creaste our own function to compute that partial dependence plot,

pdp1 = function(x1){
  nd = data.frame(X1=x1,X2=df$X2)
  mean(predict(reg,newdata=nd))
}

that will be the straight line below (the dotted line is the theoretical one y=1+x,

vx=seq(-3.5,3.5,length=101)
vpdp1 = Vectorize(pdp1)(vx)
plot(vx,vpdp1,type="l")
abline(a=1,b=1,lty=2)

which is very different from the univariate regression on x_1

abline(reg1,col="red")

Actually, the later is very consistent with a local regression, only on x_1

library(locfit)
lines(locfit(Y~X1,data=df),col="blue")

Now, to get back to the definition of the partial dependence plot, x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)], in the context of correlated variable, I was wondering if it would not make more sense to consider some local version actually, something like x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2|X_1}}[m(x_1,X_2)]. My intuition was that, somehow, it did not make any sense to consider any X_2 while X_1 was fixed (and equal to x_1). But it would make more sense actually to look at more valid X_2‘s given the value of X_1. And a natural estimate could be some k neareast-neighbors, i.e. \tilde{p}_1(x_1)=\frac{1}{k}\sum_{i\in\mathcal{V}_k(x)}^n m(x_1,x_{2,i}) where \mathcal{V}_k(x) is the set of indices of the k x_i‘s that are the closest to x, i.e.

lpdp1 = function(x1){
  nd = data.frame(X1=x1,X2=df$X2)
  idx = rank(abs(df$X1-x1))
  mean(predict(reg,newdata=nd[idx<50,]))
}
vlpdp1 = Vectorize(lpdp1)(vx)
lines(vx,vlpdp1,col="darkgreen",lwd=2)

Surprisingly (?), this local partial dependence plot gives a curve that corresponds to the simple regression…

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 = 2017
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…