# IME2020 in Montréal, canceled

The conference, that  we were organizing in July, in Montréal, is officially cancelled.

# Reprise des cours (à venir)

Lundi, les cours de la session d’hiver reprennent… à distance. Je vais mettre en ligne une série de capsules vidéos pour finir le cours, sur les GLM. J’ai mis en ligne une première vidéo (slides 0) pour annoncer le plan. J’ai fait des slides rapidement (ça changera des cours que je faisais au tableau) et je fais des enregistrement unique, sans montage, histoire de mettre en ligne rapidement le cours en ligne. Le pdf des slides est aussi en ligne (slides 0). J’essayerais de mettre les vidéos en ligne au fur et à mesure. Je ne suis pas particulièrement fier de moi, mais quitte à perdre du temps à faire le guignol devant la caméra, autant que ça serve au plus grand nombre.

Ah oui, il y a probablement des coquilles, voire des erreurs dans les slides… les commentaires sont ouverts pour me faire part de toute suggestion !

# Il était prévisible que Léonard DiCaprio ne survive pas au naufrage du Titanic

Vendredi dernier, j’ai mis en ligne un petit quiz pour prédire qui a survécu au naufrage du Titanic, en 1912. J’avais mis en ligne une sous-base avec quelques passagers

loc_fichier = "http://freakonometrics.free.fr/titanic.RData" download.file(loc_fichier, "titanic.RData") load("titanic.RData") str(base)

et plus précisément

• Survived: Passenger survival indicator (1 if survived)
• Pclass: Passenger class
• Sex: Sex of the passenger
• Age: Age of the passenger
• SibSp: Number of siblings/spouses aboard
• Parch: Number of parents/children aboard
• Embarked: Port of Embarkation
• Name: Name of passenger

La dernière variable ne sert pas à grand chose, donc on la vire,

base = base[,1:7]

On va maintenant pouvoir répondre aux questions.

# Reinforcement Learning in Economics and Finance

With Romuald Elie and Carl Remlinger we recently uploaded on ArXiv a paper on Reinforcement Learning in Economics and Finance

Reinforcement learning algorithms describe how an agent can learn an optimal action policy in a sequential decision process, through repeated experience. In a given environment, the agent policy provides him some running and terminal rewards. As in online learning, the agent learns sequentially. As in multi-armed bandit problems, when an agent picks an action, he can not infer ex-post the rewards induced by other action choices. In reinforcement learning, his actions have consequences: they influence not only rewards, but also future states of the world. The goal of reinforcement learning is to find an optimal policy — a mapping from the states of the world to the set of actions, in order to maximize cumulative reward, which is a long term strategy. Exploring might be sub-optimal on a short-term horizon but could lead to optimal long-term ones. Many problems of optimal control, popular in economics for more than forty years, can be expressed in the reinforcement learning framework, and recent advances in computational science, provided in particular by deep learning algorithms, can be used by economists in order to solve complex behavioral problems. In this article, we propose a state-of-the-art of reinforcement learning techniques, and present applications in economics, game theory, operation research and finance.

# Qui a survécu au naufrage du Titanic?

also known as quiz numéro 4 du cours STT5100 (de la session d’hiver). Vendredi dernier, alors que nous terminions le cours vers midi, François Legault a décrété l’état d’urgence sanitaire pour une quinzaine de jours. Le cours est donc sur la glace depuis une semaine, et (pour l’instant) encore une bonne semaine.

Les deux derniers cours ont été l’occasion de voir comment faire, et comment analyser, une régression logistique. Histoire de briser un peu le confinement que l’on connaît depuis une semaine, j’ai envoyé ce matin un quiz aux étudiants, et j’en profite pour le mettre en ligne (la correction arrivera dans le milieu de la semaine prochaine) pour tous ceux qui voudraient se pratiquer un peu sur les modèles de classification…

J’ai mis en ligne une base de données avec des informations sur les passagers du Titanic, qui est un paquebot qui a fait naufrage il y a plus d’un siècle (certains ont peut être vu le documentaire)

loc_fichier = "http://freakonometrics.free.fr/titanic.RData" download.file(loc_fichier, "titanic.RData") load("titanic.RData") str(base)

On a dans la base

• Survived: Passenger survival indicator (1 if survived)
• Pclass: Passenger class
• Sex: Sex of the passenger
• Age: Age of the passenger
• SibSp: Number of siblings/spouses aboard
• Parch: Number of parents/children aboard
• Embarked: Port of Embarkation
• Name: Name of passenger

La dernière variable ne sert pas à grand chose, donc on la vire,

base = base[,1:7]

J’ai ensuite posé une série de questions,

1. quelle proportion de passagers ont survécu ?
2. quelle proportion de passagers de première classe ont survécu ?
3. quelle est la statistique de test pour un test d’indépendance du chi-deux [vu en cours vendredi dernier] entre le fait de survivre (ou pas) et la classe ?
4. ensuite il faut faire une régression logistique, et prévoir la probabilité de survie pour deux passagers (fictifs) du paquebot
newbase = data.frame( Pclass = as.factor(c(1,3)), Sex = as.factor(c("female","male")), Age = c(17,20), SibSp = c(1,0), Parch = c(2,0), Embarked = as.factor(c("S","S")), Name = as.factor(c("Winslet, Miss. Kate","DiCaprio, Mr. Leonardo")))

Je ne peux mettre le formulaire de réponse en ligne ici, mais les commentaires sont ouverts pour ceux qui voudraient d’entraîner un peu pendant cette période un peu longue de confinement…

# 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}$$and$$I^\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)])&lt;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 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="") # De la qualité d’un classifieur On va profiter de la quarantaine pour mettre en ligne un billet sur la courbe ROC, la receiver operating characteristic. Considérons une petite base de données avec $n=10$ observations, deux variables continues, $x_1$ et $x_2$, et la variable d’intérêt binaire $y\in\{0,1\}$. On peut représenter les points dans le plan $(x_1,x_2)$, et on utilise une couleur différente pour $y\in\{0,1\}$. x1 = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) x2 = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) y = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x1,x2=x2,y=as.factor(y)) plot(x1,x2,col=c("red","blue")[1+y],pch=19,cex=1.5) On peut alors faire une régression logistique, de telle sorte que $$\mathbb{P}(Y=1|x_1,x_2)=\frac{e^{\beta_0+\beta_1x_1+\beta_2x_2}}{1+e^{\beta_0+\beta_1x_1+\beta_2x_2}}$$On peut visualiser l’ensemble des points $(x_1,x_2)$ pour lesquels $$\mathbb{P}(Y=0|x_1,x_2)=\mathbb{P}(Y=1|x_1,x_2)$$(on a alors autant de chance d’être rouge et bleu) soit $\beta_0+\beta_1x_1+\beta_2x_2=0$ qui correspond à une droite, reg = glm(y~x1+x2,data=df,family=binomial(link = "logit")) b = coefficients(reg) abline(a=-b[1]/b[3],b=-b[2]/b[3]) On peut alors représenter $y_i$ en fonction du score, i.e. l’estimation de $\mathbb{P}(Y=1|x_{1,i},x_{2,i})$, Y = dfy S = predict(reg,type="response") plot(S,y,xlab="probabilité prédite",ylab="y")

On va alors se donner un seuil (par exemple $50\%$) : si la probabilité que $Y$ prenne la valeur $1$ excède le seuil, on prédit $1$ (et sinon $0$). Sur la figure ci-dessus, on a alors 4 sortes de points : ceux à gauche du seuil (et qui sont prédits $0$), qui sont bien classés s’ils sont en bas, et mal classés en haut; à droite du seuil (et qui sont prédits $1$), ils sont bien classés s’ils sont en haut, et mal classés en bas [dans le code ci-dessous, le symbole &gt désigne l’opérateur “supérieur“, qui malheureusement ne passe pas dans cet éditeur]

seuil = .5 Yhat = (S&gt;seuil)*1 plot(S,y,xlab="probabilité prédite",ylab="y",pch=19, col=c("red","blue")[1+(y==Yhat)]) abline(v=seuil,lty=2)

Les couleurs reflètent le bon ou mauvais classement : les points rouges correspondent à des erreurs de classement. On peut retrouver tout ça dans le tableau de contingence ci-dessous,  qui correspond au tableau standard d’un test d’hypothèse

table(Yhat,Y) Y Yhat 0 1 0 3 1 1 1 5

Ce qui va nous intéresser ici à deux grandeurs particulières : le taux de faux positifs et le taux de vrais positifs,

 FP=sum((Ps==1)*(Y==0))/sum(Y==0) TP=sum((Ps==1)*(Y==1))/sum(Y==1)

On a obtenu ce tableau à un seuil donné (ici $50\%$) mais on peut regarder ce qui se passe lors que le seuil change, comme sur l’animation ci-dessous, où on trace, à droite, le taux de vrais positifs (sur l’axe $y$) en fonction du taux de faux positifs (sur l’axe $x$

L’ensemble des points donne la courbe ROC.

roc.curve=function(s,print=FALSE){ Ps=(S&gt;s)*1 FP=sum((Ps==1)*(Y==0))/sum(Y==0) TP=sum((Ps==1)*(Y==1))/sum(Y==1) if(print==TRUE){ print(table(Observed=Y,Predicted=Ps)) } vect=c(FP,TP) names(vect)=c("FPR","TPR") return(vect)} u = seq(0,1,length=251) V = Vectorize(roc.curve)(u) plot(t(V),type="s",xlab="Faux Positifs",ylab="Vrais Positifs") segments(0,0,1,1,col="light blue")

On peut vérifier que le point qu’on avait obtenu avec un seuil de $50\%$ est bien sur la courbe

table(Yhat,Y) Y Yhat 0 1 0 3 1 1 1 5 (FP = sum((Yhat)*(Y==0))/sum(Y==0)) [1] 0.25 (TP = sum((Yhat==1)*(Y==1))/sum(Y==1)) [1] 0.8333333 abline(v=FP,lty=2,col="blue") abline(h=TP,lty=2,col="blue") points(FP,TP,pch=19,cex=1.5)

Bien entendu, il y a (beaucoup) de packages R qui permettent d’avoir cette courbe,

library(ROCR) pred = prediction(S,Y) plot(performance(pred,"tpr","fpr"))

Une grandeur intéressante est appelée aire sous la courbe (ou AUC) qu’on peut calculer ici à la main (on a une simple fonction en escalier)

p1 = roc.curve(1/3) p2 = roc.curve(.7) p2[1]*p2[2]+(p1[1]-p2[1])*p1[2]+(1-p1[1]) [1] 0.875

mais qu’on peut avoir automatiquement

auc.perf = performance(pred, measure = "auc") auc.perf@y.values[[1]] [1] 0.875

Allez, tentons un autre classifieur : toujours une régression logistique, mais sur un facteur obtenu en coupant la seconde variable en deux, $\boldsymbol{1}_{[s,\infty)}(x_2)$

reg = glm(y~I(x2&gt;.525),data=df,family=binomial(link = "logit")) abline(h=.525)

La droite horizontale n’est plus la droite qui donne autant de chance d’être rouge que bleu, mais qui coupe la variable $x_2$. Ici, on prédit juste deux valeurs : $40\%$ de chance d’être bleu en bas, et $80\%$ de chance d’être bleu, en haut. Si on représente observations $y_i$ en fonction des probabilités prédites, on obtient

Y = df$y S = predict(reg,type="response") plot(S,y,xlab="probabilité prédite",ylab="y",xlim=0:1) Avec un seuil à $50\%$, on obtient le tableau de contingence suivant (avec 3 erreurs, contre 2 auparavant) seuil = .5 Yhat = (S&gt;seuil)*1 table(Yhat,Y) Y Yhat 0 1 0 3 2 1 1 4 Si on trace la courbe ROC, on obtient roc.curve=function(s,print=FALSE){ Ps=(S&gt;s)*1 FP=sum((Ps==1)*(Y==0))/sum(Y==0) TP=sum((Ps==1)*(Y==1))/sum(Y==1) if(print==TRUE){ print(table(Observed=Y,Predicted=Ps)) } vect=c(FP,TP) names(vect)=c("FPR","TPR") return(vect)} u = seq(0,1,length=251) V = Vectorize(roc.curve)(u) plot(t(V),type="l",xlab="Faux Positifs",ylab="Vrais Positifs") segments(0,0,1,1,col="light blue") Cette fois, la courbe n’est plus constante par morceaux, mais linéaire par morceaux, et continue… L’interprétation est un peu plus subtile: cette fois, on a deux régions de l’espace, et dans chaque région, on ne sait pas trop comment distinguer (la probabilité est plate partout ici, contraitement à la régression précédante). Autrement dit, dans cette région, on a une proba constante, par exemple $40\%$ (en bas) : quand on doit prévoir pour un individu dans cette région, on lui attribue les valeurs $\{0,1\}$ respectivement avec les probabilités $\{40\%,60\%\}$. Quand on a une probabilité constante, on parle de classifieur aléatoire… La diagonale bleue sur la figure ci-dessus est justement un classifieur aléatoire… C’est ce qu’on obtient si on prédit au hasard pred = prediction(S,Y) plot(performance(pred,"tpr","fpr")) Le point est obtenu avec un seuil de $50\%$ (ou en fait, n’importe quelle valeur entre $40\%$ et $80\%$). On peut là encore calculer l’aire sous la courbe, cette fois à l’aide de trapèzes (ou de triangles) p1 = roc.curve(.5) p2[1]*p2[2]/2+(1-p1[1])*p1[2]+(1-p1[1])*(1-p1[2])/2 [1] 0.7083333 auc.perf = performance(pred, measure = "auc") auc.perf@y.values[[1]] [1] 0.7083333 # 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+\varepsilon$$but 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(&gt;|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(&gt;|t|) (Intercept) 42.98 2.15 19.988 &lt; 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 that$$y=\beta_0+\beta_1x+\beta_2(x-s_1)_++\cdots+\beta_k(x-s_{k-1})_++\varepsilon$$where $(x-s)_+$ equals $0$ if $x 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_++\varepsilon$$with 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_++\varepsilon$$for 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…

# Gini index, poverty and top shares

Consider some ordered income $\{y_1,y_2,\dots,y_n\}$, with $y_1\leq y_2\leq\dots\leq y_n$. A classical tool to visualize inequality is Lorenz curve: define the proportion of people $F_{i}=i/n$ (with the convention $F_{0}=0$); then the cumulated wealth $S_{i}=\sum_{j=1}^{i}y_{j}$ and the fraction of cumulated wealth $L_{i}=S_{i}/S_{n}$ (with again ${\displaystyle L_{0}=0}$). Then Lorenz curve is simply the plot $\{F_i,L_i\}$ : it plots the proportion of the total income of the population ($y$ axis) that is cumulatively earned by the bottom $x$\% of the population. And Gini index is the ratio of the area that lies between the line of equality (the first diagonal, $(0,0)-(1,1)$) and the Lorenz curve over the total area under the line of equality. A simple formula would be $${\displaystyle G={\frac {2\sum _{i=1}^{n}iy_{i}}{n\sum _{i=1}^{n}y_{i}}}-{\frac {n+1}{n}}}$$but let us keep in mind simply the fact that it is simply the area below the first diagonal. Note further that the Lorenz curve is increasing, and convex. So actually, for a given Gini index – say $G=60\%$, we can have the two following situations below : on the left, 60% of the poor people get absolutely nothing, and the top 40% shares equally the remaining wealth; on the right, one person gets 60% of the wealth, and everyone else shares equally the remaining wealth.

The two areas are equals (the triangles are the same – up to some symmetrys and rotations) so the two Lorenz curve exhibit the same Gini index. On the left, the 10% the poorest own 0% of the wealth (in green) while the 10% of the richest own 25% of the wealth (in red). On the right the 10% the poorest own 4% of the wealth (in green) while the 10% of the richest own 64% of the wealth (in red). Which can be seen as some sort of paradox : the two cases exhibit the same over inequality, but the one where the poorest get more is also the one where the richest get more.

# 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…

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

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

# Nice Thematic Semester in Montréal: The Mathematics of Decision Making

We should have a lot a great talks this semester, on a very interesting topic,

# 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 !

OpenEdition