Tag Archives: R

Interprétabilité et explicabilité (formalisé) des modèles prédictifs

Dans Confessiones, Saint Augustin écrivait

quid est ergo tempus? si nemo ex me quaerat, scio; si quaerenti explicare velim, nescio

que l’on traduit

Qu’est-ce-que le temps ? Si personne ne me le demande, je le sais. Si je veux l’expliquer à qui me le demande, je ne le sais plus.

Pour aller un peu plus loin (car souvent, si on nous demande d’expliquer, on a quelques idées), dans Une étude en rouge de Sir Arthur Conan Doyle, paru en 1887, on a l’échange suivant, entre Sherlock Holmes et le docteur Watson,

– Je me demande ce que cherche ce type là-bas, demandai-je,
désignant un grand individu habillé simplement qui suivait
l’autre côté de la rue, en examinant anxieusement les numéros.
Il tenait à la main une grande enveloppe bleue et, de toute
évidence, portait un message.
– Vous parlez de ce sergent d’infanterie de marine ? dit Sherlock Holmes.

puis, comme il s’avère que la personne est effectivement sergent dans la marine (tout comme un autre personnage de l’histoire, un certain Arthur Charpentier), le docteur Holmes lui demande une explication, il veut savoir comment il est arrivé à cette conclusion

« Comment diable avez-vous pu deviner cela ? demandai-je.
– Deviner quoi ? fit-il sans aménité.
– Eh bien, qu’il était un sergent de marine en retraite ?
– Je n’ai pas de temps à perdre en bagatelles ! répondit-il
avec brusquerie avant d’ajouter dans un sourire : excusez ma rudesse ! Vous avez rompu le fil de mes pensées. Mais c’est peut-être aussi bien. Ainsi donc vous ne voyiez pas que cet homme était un sergent de marine ?
Non, certainement pas !
– Décidément, l’explication de ma méthode me coûte plus que son application ! Si l’on vous demandait de prouver que deux et deux font quatre, vous seriez peut-être embarrassé ; et cependant, vous êtes sûr qu’il en est ainsi. Malgré la largeur de la rue, j’avais pu voir une grosse ancre bleue tatouée sur le dos de la main du gaillard. Cela sentait la mer. Il avait la démarche militaire et les favoris réglementaires ; c’était, à n’en pas douter, un marin. Il avait un certain air de commandement et d’importance. Rappelez-vous son port de tête et le balancement de sa canne ! En outre, son visage annonçait un homme d’âge moyen, sérieux, respectable. Tous ces détails m’ont amené à penser qu’il était sergent.
– C’est merveilleux ! m’écriai-je.
– Peuh ! L’enfance de l’art ! dit Holmes, mais d’un air qui me
parut trahir sa satisfaction devant ma surprise et mon admiration manifestes.

(pour être honnête, c’est Liu Cixin qui en parle dans Le problème à trois corps). Pour l’anecdote, c’est la première histoire du couple Holmes-Watson, qui introduit la méthode de travail de Sherlock Holmes. Pour ceux qui sont familier avec les nouvelles, cette approche narrative sera largement reprise par la suite: Sherlock Holmes énonce un fait, le docteur Watson est étonné et demande une explication, et Sherlock Holmes explique, point par point, comment il est arrivé à cette conclusion. C’est un peu cette approche qu’on tente de mettre en place quand on va construire un modèle prédictif : sur la base des données du Titanic, si on prédit que telle personne va mourir, et que telle autre va survivre, on veut comprendre pourquoi le modèle arrive à cette conclusion.

Continue reading Interprétabilité et explicabilité (formalisé) des modèles prédictifs

Interpretability and explainability of predictive models

In 400 AD, in his Confessiones, Augustine wrote

quid est ergo tempus? si nemo ex me quaerat, scio; si quaerenti explicare velim, nescio

that can be translated as

What then is time? If no one asks me, I know what it is. If I wish to explain it to him who asks, I do not know.

To go a little further (because often, if we are asked to explain, we have some ideas), in A Study in Scarlet by Sir Arthur Conan Doyle, published in 1887, we have the following exchange, between Sherlock Holmes and Doctor Watson

– “I wonder what that fellow is looking for?” I asked, pointing to a stalwart, plainly-dressed individual who was walking slowly down the other side of the street, looking anxiously at the numbers. He had a large blue envelope in his hand, and was evidently the bearer of a message.
– “You mean the retired sergeant of Marines,” said Sherlock Holmes.

then, as it turns out that the person is indeed a sergeant in the navy (as is another character in the story, someone named Arthur Charpentier), Dr. Holmes asks him for an explanation, he wants to know how he arrived at this conclusion

– “How in the world did you deduce that?” I asked.
“Deduce what?” said he, petulantly.
“Why, that he was a retired sergeant of Marines.”
“I have no time for trifles,” he answered, brusquely; then with a smile, “Excuse my rudeness. You broke the thread of my thoughts; but perhaps it is as well. So you actually were not able to see that that man was a sergeant of Marines?”
“No, indeed.”
– “It was easier to know it than to explain why I knew it. If you were asked to prove that two and two made four, you might find some difficulty, and yet you are quite sure of the fact. Even across the street I could see a great blue anchor tattooed on the back of the fellow’s hand. That smacked of the sea. He had a military carriage, however, and regulation side whiskers. There we have the marine. He was a man with some amount of self-importance and a certain air of command. You must have observed the way in which he held his head and swung his cane. A steady, respectable, middle-aged man, too, on the face of him – all facts which led me to believe that he had been a sergeant.”

(to be honest, it is Liu Cixin who talks about it in The Three-Body Problem). For the record, this is the first story of the Holmes-Watson couple, which introduces Sherlock Holmes’ working method. For those who are familiar with the short stories, this narrative approach will be widely used thereafter: Sherlock Holmes states a fact, Dr. Watson is astonished and asks for an explanation, and Sherlock Holmes explains, point by point, how he arrived at this conclusion. This is a bit like the approach we try to implement when we build a predictive model: on the basis of the Titanic data, if we predict that such and such a person will die, and that such and such a person will survive, we want to understand why the model arrives at this conclusion.
Continue reading Interpretability and explainability of predictive models

Report des voix entre les deux tours aux élections présidentielles en France

Allez, on va profiter des élections présidentielles, à deux tours, en France, pour jouer un peu, maintenant que les données commencent à être disponibles. On peut avoir les données par département (une centaine) ou par commune (36000)

# département
loc_dpt_1 = "https://www.data.gouv.fr/fr/datasets/r/312b99d0-27a7-429e-be8a-b930975545e9"
loc_dpt_2 = "https://www.data.gouv.fr/fr/datasets/r/c7b72e78-5e78-4976-995f-30533f0a35d3"
# commune
loc_com_1 = "https://www.data.gouv.fr/fr/datasets/r/54782507-e795-4f9d-aa70-ed06feba22e3"
loc_com_2 = "https://www.data.gouv.fr/fr/datasets/r/1601b027-48d1-42b1-96f1-b0fb49f55aa3"
download.file(loc_dpt_1,"dpt1.csv")
download.file(loc_dpt_2,"dpt2.csv")
d1 = read.csv("dpt1.csv")
d2 = read.csv("dpt2.csv")
download.file(loc_com_1,"com1.csv")
download.file(loc_com_2,"com2.csv")
c1 = read.csv("com1.csv")
c2 = read.csv("com2.csv")

Commençons par les départements. On peut récupérer les pourcentages de voix obtenus par les différents candidats, au premier tour, et regarder les pourcentages obtenus par Marine Le Pen. Par exemple, on peut regarder, département par département, le score de Marine Le Pen au premier et au second tour, pour vérifier que le report des voix s’opère bien

D1 = d1[d1$cand_nom=="LE PEN",c("dep_code","cand_rapport_inscrits","cand_rapport_exprim")]
D2 = d2[,c("CodeDépartement","LE.PEN.ins","LE.PEN.exp")]
names(D1) = c("dep_code","ins1","exp1")
names(D2) = c("dep_code","ins2","exp2")
D = merge(D1,D2)
par(mfrow=c(1,2))
plot(D$ins1,D$ins2,xlab="inscrits (%), premier tour",ylab="inscrits (%), second tour",col=colr2[1])
plot(D$exp1,D$exp2,xlab="exprimés (%), premier tour",ylab="exprimés (%), second tour",col=colr2[2],main="Le Pen / Le Pen",cex.main=1.1)

Comme on pouvait s’y attendre, on a une belle droite de régression, avec les nombres de voix obtenus, relativement aux inscrits, en bleu, et aux exprimés, à droite. Pour Emmanuel Macron au premier et Marine Le Pen au second tour, on a une belle corrélation négative, mais moins prononcée qu’auparavant. Autrement dit, un département qui votait peu Macron au premier tour a beaucoup voté Le Pen au second

Mais là, on parlait des deux candidats qui se sont maintenus. Si on regarde l’autre candidat de l’extrême droite, les reports sont assez étranges… après, les taux étaient “relativement” faibles au premier tour pour Erik Zemmour.

Pour le candidat en tête à gauche, on a une corrélation négative, mais sur les votes exprimées, on observera des outliers en haut à droite, avec manifestement un vote d’opposition au président sortant dans certains départements

Pour le candidat écologiste, la corrélation est elle encore plus clair

autrement dit, quand on a voté écologiste au premier tour, on avait peu tendance à voter Le Pen (ou disons les départements où on votait fortement écologiste n’ont pas voté Le Pen au second #ecologicalfallacy, etc).

Mais on peut faire le même exercice par commune. Mon seul soucis a été que les codes communes dans les deux bases étaient différents, donc j’ai utilisé les noms des communes, ce que je n’aime pas faire (avec les accents, les tirets, etc, la fusion peut être assez mauvaises… mais là ça marche bien). Je me suis ensuite limité aux grosses communes, de plus de 10000 inscrits,

C1 = c1[c1$cand_nom=="LE PEN",c("commune_name","cand_rapport_inscrits","cand_rapport_exprim","votants_nb")]
C2 = c2[,c("Commune","LE.PEN.ins","LE.PEN.exp")]
names(C1) = c("com_code","ins1","exp1","taille")
names(C2) = c("com_code","ins2","exp2")
D = merge(C1,C2)
D = D[D$taille>5000,]
par(mfrow=c(1,2))
plot(D$ins1,D$ins2,xlab="inscrits (%), premier tour",ylab="inscrits (%), second tour",col=colr2[1])
plot(D$exp1,D$exp2,xlab="exprimés (%), premier tour",ylab="exprimés (%), second tour",col=colr2[2],main="Le Pen / Le Pen",cex.main=1.1)

On peut alors faire le même exercice, avec Le Pen au premier et second tour

ou Macron au premier, Le Pen au second

On observe des beaux nuages, sur la gauche, respectivement en haut et en bas… Beaucoup de communes qui avaient peu voté pour Le Pen semblent avoir eu peur, et on massivement voté au second pour elle. Alors que ce n’est pas trop le cas pour celles où Le Pen avait obtenu beaucoup de voix. Et pour Macron, c’est en bas à gauche, autrement dit, beaucoup de communes qui n’ont pas voté Macron au premier tour n’ont pas pour autant voté Le Pen au second. Pour ceux que ça amuse, les codes sont disponibles pour regarder qui sont ces communes dans les coins.

Pour Zemmour, c’est assez confus, j’avoue, le personnage est déroutant (ou dégoutant ? comme le suggérait le correcteur automatique)

ou alors Jean Luc Mélanchon

ou la relativement franche corrélation des écologistes, qui semblent résolument aller en sens inverse avec le vote d’extrême droite…

 

Insurance Pricing Game

Would you like to put your data science skills to the test?

Imperial College London, Universite du Quebec à Montreal (UQAM), and actuarial institutes in Singapore, the UK, including the IFoA, and Australia, ASTIN, the Casualty Actuarial Society are co-organising a global data science competition.

Would you like to accurately predict the cost of insurance by putting your data science skills to the test? We are hosting two competitions with separate datasets, a loss prediction competition on Kaggle with synthetic workers’ compensation data, and a pricing competition in a simulated market hosted on AI Crowd with real-world motor insurance contracts. Codes can be either in R or python. The competition is being sponsored by a number of different organisations, with a total of US$12,000 in cash prizes to be won. For more information about how to take part please visit www.pricing-game.com

Hidding values in the output of the summary function for a (linear) regression

Since our Fall 2020 session will be 100% online (and off-site), I have to work hard this summer to prepare online quizz and exams. I started intensively to play with Achim’s awesome r-exams package. But there are still a few things that I wanted to add, so I will post a series of posts on my blogs to keep tracks of updated functions I will write. Most of them a modification of R internal functions, so the code is hard to read. Here is the file, and I will update it frequently

url = "http://freakonometrics.free.fr/onlineExams.R"
source(url)

I have updated the summary function (more precisely the summary.lm function). To see how it works, consider a simple regression

library(car)
reg = lm(prestige ~ women, data=Prestige)
my_summary(reg)
 
Call:
lm(formula = prestige ~ women, data = Prestige)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-33.444 -12.391  -4.126  13.034  39.185 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 48.69300    2.30760  21.101    2e-16 ***
women       -0.06417    0.05385  -1.192    0.236    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 17.17 on 100 degrees of freedom
Multiple R-squared:  0.014,	Adjusted R-squared:  0.004143
F-statistic:  1.42 on 1 and 100 DF,  p-value: 0.2362

A classical question I ask in my quizz is to hide the p-value of the F-test, and ask what it is (to make sure that students understand the equivalence between the F and the t test, in a simple regression). To hide the p-value, use

my_summary(reg, Fisher=TRUE)
 
Call:
lm(formula = prestige ~ women, data = Prestige)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-33.444 -12.391  -4.126  13.034  39.185 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 48.69300    2.30760  21.101    2e-16 ***
women       -0.06417    0.05385  -1.192    0.236    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 17.17 on 100 degrees of freedom
Multiple R-squared:  0.014,	Adjusted R-squared:  0.004143
F-statistic:  1.42 on 1 and 100 DF,  p-value: ■■■■■

(and then, in a multiple choice exam, I ask if it is 1%, 5%, 12%, 23%, 47%, for example). That one was easy, since all those lines are based on the cat function, so I just modify it, if necessary

if(Fisher) cat("\nF-statistic:", formatC(x$fstatistic[1L], 
    digits = digits), "on", x$fstatistic[2L], "and", 
        x$fstatistic[3L], "DF,  p-value:", "■■■■■")
    if(!Fisher) cat("\nF-statistic:", formatC(x$fstatistic[1L], 
    digits = digits), "on", x$fstatistic[2L], "and", 
                   x$fstatistic[3L], "DF,  p-value:", format.pval(pf(x$fstatistic[1L], 
                   x$fstatistic[2L], x$fstatistic[3L], lower.tail = FALSE), 
                   digits = digits))

(here I use the unicode ‘black square‘ symbol to hide numbers). Of course, I can hide the value of \sigma, or the (adjusted or not) R ^2, etc.

Now, a little bit more tricky: what if we want to change the regression table, with the coefficients, their standard deviation, etc.  It is tricky since those values are numeric, with an appropriate format (not too many digits), but it can be done easily since that formating is done through the printCoefmat function.  So in my code, I have my internal function, where I ask to put some ‘black square‘ (and the good number to keep a readable format) at some specific locations. Consider a more complex regression

reg = lm(prestige ~ ., data=Prestige)

and assume that we want to hide the value of the intercet, \widehat{\beta}_0 (i.e. located at (1,1) in the matrix) and the p-value of the t-test for the fourth one (i.e. located at (4,4) in the matrix – since the first colum is \widehat{\beta}_3, the second one its standard deviation, the thirst one the t value, and then, the fourth one, the p-value of the test). I use the following two vectors

vligne = c(1,4),
vcolonne = c(1,4)

with rows and columns in the matrix (of course, the two should have the same length). Then, the good thing is that the printCoefmat function convert numerical values into characters (to have things that look like columns actually). So we simply have to remove numerical digits, and use squares instead

Cf2=Cf
  if(length(vligne)>0){  
    for(i in 1:length(vligne)){
      long = nchar(Cf[vligne[i],vcolonne[i]])
      Cf2[vligne[i],vcolonne[i]] = paste(rep("■",long),collapse = "")
    }}

Then, we print the updated version of the table

print.default(Cf2, quote = quote, right = right, na.print = na.print,...)

For example, here, it would be

my_summary(reg, vligne=c(1,4), vcolonne=c(1,4))
 
Call:
lm(formula = prestige ~ ., data = Prestige)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-12.9863  -4.9813   0.6983   4.8690  19.2402 
 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) ■■■■■■■■■■  8.018e+00  -1.513  0.13380    
education    3.933e+00  6.535e-01   6.019 3.64e-08 ***
income       9.946e-04  2.601e-04   3.824  0.00024 ***
women        1.310e-02  3.019e-02   0.434  ■■■■■■■    
census       1.156e-03  6.183e-04   1.870  0.06471 .  
typeprof     1.077e+01  4.676e+00   2.303  0.02354 *  
typewc       2.877e-01  3.139e+00   0.092  0.92718    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 7.037 on 91 degrees of freedom
  (4 observations deleted due to missingness)
Multiple R-squared:  0.841,	Adjusted R-squared:  0.8306
F-statistic: 80.25 on 6 and 91 DF,  p-value: < 2.2e-16

Of course, it is hand made, I do not check for typos (like you should not ask to put squares in the seventh column), but that works well enough to generate random regressions in a quizz (or identical regressions on subsamples of a large dataset), and to hide values, in a quizz.

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.

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.

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.

Function basis and regression

In the first part of the course on linear models, we’ve seen how to construct a linear model when the vector of covariates \boldsymbol{x} is given, so that \mathbb{E}(Y|\boldsymbol{X}=\boldsymbol{x}) is either simply \boldsymbol{x}^\top\boldsymbol{\beta} (for standard linear models) or a functional of \boldsymbol{x}^\top\boldsymbol{\beta} (in GLMs). But more generally, we can consider transformations of the covariates, so that a linear model can be used. In a very general setting, consider \sum_{j=1}^m\beta_j h_j(\boldsymbol{x})with h_j:\mathbb{R}^p\rightarrow\mathbb{R}. The standard linear model is obtained when m=p and h_j(\boldsymbol{x})=x_j , but of course, much more general models can be obtained, for instance with h_k(\boldsymbol{x})=x_j^2 or h_k(\boldsymbol{x})=x_{j}x_{j'}, that could be used to achieve high-order Taylor expansions. In that case, we will obtain the polynomial regression, that we will discuss first. We might also think of piecewise constant functions, h_k(\boldsymbol{x})=\boldsymbol{1}(x_j\in [a,b]) , that could be related to regression trees (but that is not in the scope in the STT5100 course). And if we go on step futher, we might think of piecewise linear or piecewise polynomial function, possibly with additional continuity constraints, that will lead us to spline basis.

  • Polynomial regression

For pedagogical purpose, when I talk about polynomial regression, I always have in mind (in the univariate case) y=\beta_0+\beta_1x+\beta_2x^2+\cdots+\beta_kx^k+\varepsilonbut if we use

lm(y~poly(x,k))

in R, the output is not the \beta_j‘s.

As discussed in Kennedy & Gentle (1980) Statistical Computing,

Recall that orthogonal polynomials are defined with respect to the classical inner-product (on the finite interval (a,b)){\displaystyle \langle f,g\rangle =\int _{a}^{b}f(x)g(x)~\mathrm {d} x} And a sequence of orthogonal polynomials is (P_n) where P_n is a polynomial of degree n, for all n, and such that P_m\perp P_n for all m\neq n. Note that those polyomials are orthogonal with respect to the inner product defined above, i.e. given some finite interval (a,b). But if (a,b) changes, the polynomials will be different.

A popular family of orthogonal polynomial, on finite interval (-1,+1) is the family of Legendre polynomials, satisfying{\displaystyle \int _{-1}^{1}P_{m}(x)P_{n}(x)~\mathrm {d} x=0}as soon as m\neq n. Those polynomials satisfy Bonnet’s recursion formula{\displaystyle (n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x)} or Rodrigues’ formula {\displaystyle P_{n}(x)={\frac {1}{2^{n}n!}}{\frac {d^{n}}{dx^{n}}}(x^{2}-1)^{n}}The first values are here{\displaystyle P_{0}(x)=1} {\displaystyle P_{1}(x)=x}{\displaystyle P_{2}(x)={\frac {3x^{2}-1}{2}}}{\displaystyle P_{3}(x)={\frac {5x^{3}-3x}{2}}} {\displaystyle P_{4}(x)={\frac {35x^{4}-30x^{2}+3}{8}}}

Interestingly, we can get those polynomial functions using

library(orthopolynom)
(leg4coef = legendre.polynomials(n=4))
[[1]]
1 
 
[[2]]
x 
 
[[3]]
-0.5 + 1.5*x^2 
 
[[4]]
-1.5*x + 2.5*x^3 
 
[[5]]
0.375 - 3.75*x^2 + 4.375*x^4

Of course, there are many families of orthogonal polynomials (Jacobi polynomials, Laguerre polynomials, Hermite polynomials, etc). Now, in R, there is the standard poly function, that we use in polynomial regression.

x = seq(-1,1,length=101)
y = poly(x,4)
y
                   1            2             3            4
  [1,] -1.706475e-01  0.215984813 -2.480753e-01  0.270362873
  [2,] -1.672345e-01  0.203025724 -2.183063e-01  0.216290298
...
[100,]  1.672345e-01  0.203025724  2.183063e-01  0.216290298
[101,]  1.706475e-01  0.215984813  2.480753e-01  0.270362873
attr(,"coefs")
attr(,"coefs")$alpha
[1] 3.157229e-17 2.655145e-16 9.799244e-17 5.368224e-16
 
attr(,"coefs")$norm2
[1]   1.0000000 101.0000000  34.3400000   9.3377328   2.4472330   0.6330176
 
attr(,"degree")
[1] 1 2 3 4
attr(,"class")
[1] "poly"   "matrix"

But these are not Legendre polynomials… As explained in 李哲源‘s post on stackoverflow, the idea is to start with P_{-1}(x)=0, P_{0}(x)=1 and P_{1}(x)=x, and then define \ell_n=\langle P_n,P_n\rangle  as well as \alpha_n=\langle P_nP_1,P_1\rangle/\ell_n=\langle P_n^2,P_1\rangle/\ell_i= and \beta_n=\ell_n/\ell_{n-1}. Finally, define recursively{\displaystyle P_{n}(x)=(x-\alpha_{n-1})P_{n-1}(x)-\beta_{i-1}P_{i-2}(x)}and its normalized version, \tilde{P}_{n}=P_n/\sqrt{\ell_n}. That is what poly computes.

So, for pedagogical purpose, I said that I like to use y=\boldsymbol{x}^\top\boldsymbol{\beta}+\varepsilon where\boldsymbol{x}=(1,x,x^2,\cdots,xˆ{k-1},x^k)And actually, when using poly, we use the QR decomposition of that matrix. As discussed in in 李哲源‘s post, we can almost reproduce the poly function using

my_poly - function (x, degree = 1) {
    xbar = mean(x)
    x = x - xbar
    QR = qr(outer(x, 0:degree, "^"))
    X = qr.qy(QR, diag(diag(QR$qr), length(x), degree + 1))[, -1, drop = FALSE]
    X2 = X * X
    norm2 = colSums(X * X)   
    alpha = drop(crossprod(X2, x)) / norm2
    beta = norm2 / (c(length(x), norm2[-degree]))
    colnames(X) = 1:degree
    scale = sqrt(norm2)
    X = X * rep(1 / scale, each = length(x))
    X}

Nevertheless, the two models are equivalent. More precisely,

plot(cars)
reg1 = lm(dist~speed+I(speed^2)+I(speed^3),data=cars)
reg2 = lm(dist~poly(speed,3),data=cars)
u = seq(3,26,by=.1)
v1 = predict(reg1,newdata=data.frame(speed=u))
v2 = predict(reg2,newdata=data.frame(speed=u))
lines(u,v1,col="blue")
lines(u,v2,col="red",lty=2)

We have exactly the same prediction here

v1[u==15]
     121 
38.43919 
v2[u==15]
     121 
38.43919

And probably also quite interesting : the coefficients do not have the same interpretation (since we do not have the same basis), but the p-value for the highest degree is exactly the same here ! Here the two models reject, with the same confidence, the polynomial of degree three,

summary(reg1)
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -19.50505   28.40530  -0.687    0.496
speed         6.80111    6.80113   1.000    0.323
I(speed^2)   -0.34966    0.49988  -0.699    0.488
I(speed^3)    0.01025    0.01130   0.907    0.369
 
Residual standard error: 15.2 on 46 degrees of freedom
Multiple R-squared:  0.6732,	Adjusted R-squared:  0.6519 
F-statistic: 31.58 on 3 and 46 DF,  p-value: 3.074e-11
 
summary(reg2)
 
Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)        42.98       2.15  19.988  < 2e-16 ***
poly(speed, 3)1   145.55      15.21   9.573  1.6e-12 ***
poly(speed, 3)2    23.00      15.21   1.512    0.137    
poly(speed, 3)3    13.80      15.21   0.907    0.369    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.2 on 46 degrees of freedom
Multiple R-squared:  0.6732,	Adjusted R-squared:  0.6519 
F-statistic: 31.58 on 3 and 46 DF,  p-value: 3.074e-11
  • B-splines regression (and GAMs)

Splines are also important in regression models, especially when we start talking about Generalized Additive Models. See Perperoglou, Sauerbrei, Abrahamowicz & Schmid (2019) for a review. In the univariate case, I introduce (linear) splines through positive parts, in the sense thaty=\beta_0+\beta_1x+\beta_2(x-s_1)_++\cdots+\beta_k(x-s_{k-1})_++\varepsilonwhere (x-s)_+ equals 0 if x<s and x-s if x>s. Those functions are nice since they are continuous, so the model is continuous (the weighted sum of continuous functions is continuous). And we can go one step further, with y=\beta_0+\beta_1x+\beta_2x^2+\beta_3(x-s_1)^2_++\cdots+\beta_k(x-s_{k-2})^2_++\varepsilonwith quadratic splines, or y=\beta_0+\beta_1x+\beta_2x^2+\beta_3x^3+\beta_4(x-s_1)^3_++\cdots+\beta_k(x-s_{k-3})^3_++\varepsilonfor cubic splines. Interestingly, quadratic splines are not only continuous, but their first derivative is also continuous (and the second one for cubic splines). So the knot discontinuity is s_1,s_2,\cdots is now invisible…

I like those models since they are easy to interprete. For example, the simple model \beta_1 x+\beta_2(x-s)_+ is the following piecewise linear function, continuous, with a “rupture” at knot s.

Observe also the following interpretation: for small values of x, there is a linear increase, with slope \beta_1, and for lager values of x, there is a linear decrease, with slope \beta_1+\beta_2. Hence, \beta_2 is interpreted as a change of the slope.

Unfortunately, it is now what R is using when using the bs function in R, which are the standard B-splines. Just to visualize (I will skip the maths here), with R, we have

library(splines)
clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02")
x = seq(5,25,by=.25)
B = bs(x,knots=c(10,20),Boundary.knots=c(5,55),degre=1)
matplot(x,B,type="l",lty=1,lwd=2,col=clr6)
B=bs(x,knots=c(10,20),Boundary.knots=c(5,55),degre=2)
matplot(x,B,type="l",col=clr6,lty=1,lwd=2)

while the functions I mentioned were (more or less) the following

pos = function(x,s) (x-s)*(x&gt;s)
par(mfrow=c(1,2))
clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02")
x = seq(5,25,by=.25)
B = cbind(pos(x,5),pos(x,10),pos(x,20))
matplot(x,B,type="l",lty=1,lwd=2,col=clr6)
pos2 = function(x,s) (x-s)^2*(x&gt;s)
B = cbind(pos(x,5)*20,pos2(x,5),pos2(x,10),pos2(x,20))
matplot(x,B,type="l",col=clr6,lty=1,lwd=2)

And as for the polynomial regression, the two models are equivalent. For example

plot(cars)
reg1 = lm(dist~speed+pos(speed,10)+pos(speed,20),data=cars)
reg2 = lm(dist~bs(speed,degree=1,knots=c(10,20)),data=cars)
v1 = predict(reg1,newdata=data.frame(speed=u))
v2 = predict(reg2,newdata=data.frame(speed=u))
lines(u,v1,col="blue")
lines(u,v2,col="red",lty=2)

or more specifically

v1[u==15]
     121 
39.35747 
v2[u==15]
     121 
39.35747

So one more time, the two models are equivalent, but I still find the approach with the positive part more intuitive, and easy to understand. As well as the interpretation of coefficients,

summary(reg1)
 
Coefficients:
               Estimate Std. Error t value Pr(&gt;|t|)  
(Intercept)     -7.6305    16.2941  -0.468   0.6418  
speed            3.0630     1.8238   1.679   0.0998 .
pos(speed, 10)   0.2087     2.2453   0.093   0.9263  
pos(speed, 20)   4.2812     2.2843   1.874   0.0673 .
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15 on 46 degrees of freedom
Multiple R-squared:  0.6821,	Adjusted R-squared:  0.6613 
F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11
 
summary(reg2)
 
Coefficients:
                                          Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)                                  4.621      9.344   0.495   0.6233    
bs(speed, degree = 1, knots = c(10, 20))1   18.378     10.943   1.679   0.0998 .  
bs(speed, degree = 1, knots = c(10, 20))2   51.094     10.040   5.089 6.51e-06 ***
bs(speed, degree = 1, knots = c(10, 20))3   88.859     12.047   7.376 2.49e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15 on 46 degrees of freedom
Multiple R-squared:  0.6821,	Adjusted R-squared:  0.6613 
F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11

Here we can see directly that the first knot was not interesting (the slope did not change significantly) while the second one was…

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 !

Combiner les modalités d’une variable factorielle

Un billet rapide pour reprendre un point qu’on a vu ce matin en cours STT5100 pour illustrer le test de Fisher. On va utiliser les données de prix d’appartements en Pologne (données pas mal utilisées dans mon ébauche de notes de cours)

library(DALEX)
data(apartments)
with(data = apartments, boxplot(m2.price ~ district))

On a envie de faire ici des regroupements de modalités (c’est d’ailleurs suggéré par la régression simple, 5 variables explicatives étant ici non significatives). Pour mieux voir, on peut réordonner les modalités

A = with(data = apartments, aggregate(m2.price,by=list(district),FUN=mean))
A = A[order(A$x),]
L = as.character(A$Group.1)
apartments$district = factor(apartments$district, level=L)
with(data = apartments, boxplot(m2.price ~ district))

On va prendre ici le district le moins cher comme référence,

reg=lm(m2.price ~ district, data=apartments)
&gt; summary(reg)
 
Coefficients:
                    Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)          2968.36      58.02  51.160   &lt;2e-16 ***
districtBielany        17.38      84.16   0.207    0.836    
districtPraga          26.45      85.12   0.311    0.756    
districtUrsynow        42.01      82.65   0.508    0.611    
districtBemowo         80.10      83.71   0.957    0.339    
districtUrsus         102.01      82.25   1.240    0.215    
districtZoliborz      829.59      83.94   9.884   &lt;2e-16 ***
districtMokotow       887.10      81.86  10.837   &lt;2e-16 ***
districtOchota        987.93      84.16  11.738   &lt;2e-16 ***
districtSrodmiescie  2214.39      83.28  26.591   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 597.4 on 990 degrees of freedom
Multiple R-squared:  0.5698,	Adjusted R-squared:  0.5659 
F-statistic: 145.7 on 9 and 990 DF,  p-value: &lt; 2.2e-16

On peut tester si les 5 premières modalités sont nulles, ce qui est un test multiple, et on va utiliser le test de Ficher :

library(car)
linearHypothesis(reg, c("districtBielany = 0", 
                        "districtPraga = 0",
                        "districtUrsynow = 0",
                        "districtBemowo = 0",
                        "districtUrsus = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F Pr(&gt;F)
1    995 354051715                           
2    990 353269202  5    782513 0.4386 0.8217

La statistique de Fisher est faible, et avec une p-value de 82%. On peut tenter le diable, et rajouter encore une modalité

library(car)
linearHypothesis(reg, c("districtBielany = 0", 
                        "districtPraga = 0",
                        "districtUrsynow = 0",
                        "districtBemowo = 0",
                        "districtUrsus = 0",
                        "districtZoliborz = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F    Pr(&gt;F)    
1    996 405455409                                  
2    990 353269202  6  52186207 24.374 &lt; 2.2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Mais on a peut être été trop gourmand cette fois. On va regrouper les 6 premières modalités (et appeler A le regroupement de districts). Si on regarde les prix moyens, par districts, on obtient

levels(apartments$district) = c(rep("A",6),levels(apartments$district)[7:11])
with(data = apartments, boxplot(m2.price ~ district))

apartments$district = relevel(apartments$district,"Zoliborz")

On recommence, en mettant le district le moins cher comme référence, et on veut tester si les deux suivants ont des coefficients nuls dans la régression linéaire.

reg=lm(m2.price ~ district, data=apartments)
linearHypothesis(reg, c("districtMokotow = 0",
                        "districtOchota = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F Pr(&gt;F)
1    997 355292524                           
2    995 354051715  2   1240809 1.7435 0.1754

Avec une p-value de 17%, on peut accepter de regrouper les trois modalités ensemble. On a alors trois groupes de districts, dont les noms sont A, B et C. On obtient les boîtes à moustaches suivantes

levels(apartments$district) = c("B","A",rep("B",2),"C")
apartments$district = relevel(apartments$district,"A")
with(data = apartments, boxplot(m2.price ~ district))

Je laisse les plus courageux vérifier, mais on a trois districts effectivement différents, et si le but est de prévoir le prix des logements, inutiles d’utiliser un découpage avec 10 modalités, un découpage avec 3 suffit !