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

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

# Testing for a causal effect (with 2 time series)

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

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

and I tried to explain that it was just stupid : if we consider the regression of the temperature on day $t+1$ against the number of cyclist on day $t$, the $R^2$ exceeds 80%… but it is hard to claim that the number of cyclists on specific day will actually cause the temperature on the next day… Nevertheless, that was frustrating, and I was wondering if there was a clever way to test for causality in that case. A popular one is Granger causality (I can mention a paper we published a few years ago where we use such a test, Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media). To explain that test, consider a bivariate time series (just like the one we have here), $\boldsymbol{z}_t=(x_t,y_t)$, and consider some bivariate autoregressive model
$${\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) $\mathbb{E} (\boldsymbol{\varepsilon}_{t})=\boldsymbol{0}$ (the noise is centered) (ii) $\mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t}^\top)=\Omega }$, so the variance matrix is constant, but possibly non-diagonal (iii) $\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$${\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  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 ! # Lasso Regression (home made) Again, this post is related to my MAT7381 course, where we will see that it is actually possible to write our own code to compute Lasso regression, $$\min\left\lbrace\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|_{\ell_2}^2+\lambda\|\mathbf{\beta}\|_{\ell_1}\right\rbrace$$We have to define the soft-thresholding function$$S(z,\gamma)=\text{sign}(z)\cdot(|z|-\gamma)_+=\begin{cases}z-\gamma&\text{ if }\gamma>|z|\text{ and }z<0\\z+\gamma&\text{ if }\gamma<|z|\text{ and }z<0 \\0&\text{ if }\gamma\geq|z|\end{cases}$$The R function would be soft_thresholding = function(x,a){ sign(x) * pmax(abs(x)-a,0) } To solve our optimization problem, set$$\mathbf{r}_j=\mathbf{y} - \left(\beta_0\mathbf{1}+\sum_{k\neq j}\beta_k\mathbf{x}_k\right)=\mathbf{y}-\widehat{\mathbf{y}}^{(j)}$$ so that the optimization problem can be written, equivalently $$\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p [\mathbf{r}_j-\beta_j\mathbf{x}_j]^2+\lambda |\beta_j|\right\rbrace$$ hence$$\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p \beta_j^2\|\mathbf{x}_j\|-2\beta_j\mathbf{r}_j^T\mathbf{x}_j+\lambda |\beta_j|\right\rbrace$$ and one gets $$\beta_{j,\lambda} = \frac{1}{\|\mathbf{x}_j\|^2}S(\mathbf{r}_j^T\mathbf{x}_j,n\lambda)$$ or, if we develop $$\beta_{j,\lambda} = \frac{1}{\sum_i x_{ij}^2}S\left(\sum_ix_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)$$ Again, if there are weights $\mathbf{\omega}=(\omega_i)$, the coordinate-wise update becomes $$\beta_{j,\lambda,{\color{red}{\omega}}} = \frac{1}{\sum_i {\color{red}{\omega_i}}x_{ij}^2}S\left(\sum_i{\color{red}{\omega_i}}x_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)$$ The code to compute this componentwise descent is lasso_coord_desc = function(X,y,beta,lambda,tol=1e-6,maxiter=1000){ beta = as.matrix(beta) X = as.matrix(X) omega = rep(1/length(y),length(y)) obj = numeric(length=(maxiter+1)) betalist = list(length(maxiter+1)) betalist[] = beta beta0list = numeric(length(maxiter+1)) beta0 = sum(y-X%*%beta)/(length(y)) beta0list = beta0 for (j in 1:maxiter){ for (k in 1:length(beta)){ r = y - X[,-k]%*%beta[-k] - beta0*rep(1,length(y)) beta[k] = (1/sum(omega*X[,k]^2))* soft_thresholding(t(omega*r)%*%X[,k],length(y)*lambda) } beta0 = sum(y-X%*%beta)/(length(y)) beta0list[j+1] = beta0 betalist[[j+1]] = beta obj[j] = (1/2)*(1/length(y))*norm(omega*(y - X%*%beta - beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta)) if (norm(rbind(beta0list[j],betalist[[j]]) - rbind(beta0,beta),'F') &lt; tol) { break } } return(list(obj=obj[1:j],beta=beta,intercept=beta0)) } For instance, consider the following (simple) dataset, with three covariates chicago = read.table("http://freakonometrics.free.fr/chicago.txt",header=TRUE,sep=";") that we can “normalize” (or “standardize“) X = model.matrix(lm(Fire~.,data=chicago))[,2:4] for(j in 1:3) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j]) y = chicago$Fire y = (y-mean(y))/sd(y)

To initialize the algorithm, use the OLS estimate

beta_init = lm(Fire~0+.,data=chicago)$coef For instance lasso_coord_desc(X,y,beta_init,lambda=.001)$obj  0.001014426 0.001008009 0.001009558 0.001011094 0.001011119 0.001011119   $beta [,1] X_1 0.0000000 X_2 0.3836087 X_3 -0.5026137$intercept  2.060999e-16

and we can get the standard Lasso plot by looping, # Quantile Regression (home made, part 2)

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

## Median

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

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

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

To illustrate, consider a sample from a lognormal distribution,

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

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

So far so good…

## Quantile Regression

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

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

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

require(lpSolve) tau = .3 n=nrow(base) X = cbind( 1, base$area) y = base$rent_euro K = ncol(X) N = nrow(X) A = cbind(X,-X,diag(N),-diag(N)) c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N)) b = base$rent_euro const_type = rep("=",N) r = lp("min",c,A,const_type,b) beta = r$sol[1:K] - r$sol[(1:K+K)] beta  148.946864 3.289674 Of course, we can use R function to fit that model library(quantreg) rq(rent_euro~area, tau=tau, data=base) Coefficients: (Intercept) area 148.946864 3.289674 Here again, it seems to work quite well. We can use a different probability level, of course, and get a plot plot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")), ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5) sf=0:250 yr=r$solution[2*n+1]+r$solution[2*n+2]*sf lines(sf,yr,lwd=2,col="blue") tau = .9 r = lp("min",c,A,const_type,b) tail(r$solution,2)  121.815505 7.865536 yr=r$solution[2*n+1]+r$solution[2*n+2]*sf lines(sf,yr,lwd=2,col="blue") And we can adapt the later to multiple regressions, of course,

X = cbind(1,base$area,base$yearc) K = ncol(X) N = nrow(X) A = cbind(X,-X,diag(N),-diag(N)) c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N)) b = base$rent_euro const_type = rep("=",N) r = lp("min",c,A,const_type,b) beta = r$sol[1:K] - r$sol[(1:K+K)] beta  -5542.503252 3.978135 2.887234 to be compared with library(quantreg) rq(rent_euro~ area + yearc, tau=tau, data=base) Coefficients: (Intercept) area yearc -5542.503252 3.978135 2.887234 Degrees of freedom: 4571 total; 4568 residual # De la pratique de la régression Depuis le début de la session, j’ai imposé une petite innovation, en donnant, environ une semaine sur deux, un petit exercice (obligatoire mais non noté) avant le cours, en vue de forcer à réfléchir (et de donner des éléments de réponse). Par exemple pour le premier cours, il fallait “prévoir” une valeur manquante, et le but était de montrer que, naturellement, on choisit la valeur moyenne. Pour demain, j’avais posé un exercice un peu plus compliqué, sachant qu’on avait vu, lors du dernier cours, comme faire une régression linéaire, et qu’on avait fini en discutant les tests simples (en lien avec la significativité) et les tests multiples. Pour l’exercice, j’avais mis en ligne une petite base de données, download.file("http://freakonometrics.free.fr/data3.RData","data3.RData") load("data3.RData") str(df) 'data.frame': 147 obs. of 10 variables:$ Y : num 11.72 15.91 14.19 11.15 8.31 ... $X1: num 1.33 3.18 0.28 2.08 0.11 1.67 1.97 1.27 4.38 0.52 ...$ X2: num 3.66 3.75 3.32 2.68 4.97 2.98 4.56 1.78 2.83 6.36 ... $X3: num 1.41 3.01 0.34 2.19 0.25 1.69 2.01 1.25 4.41 0.43 ...$ X4: num -3.53 -4.46 -3.35 -7.54 -7.02 -2.53 -6.1 -5.99 -3.92 -5.84 ... $X5: num 0.57 0.01 -0.7 1.62 -0.95 -1.37 1.18 -0.72 2.63 -1.63 ...$ X6: num -0.82 1.2 3.03 -0.91 -1.6 1.77 1 -1.33 1.31 -0.7 ... $X7: num 1.01 0.06 2.02 3.63 2.66 2.53 1.29 3.5 1.17 1.8 ...$ X8: num 8.31 8.52 9.78 7.34 7.26 ... \$ X9: num 6.04 6.53 7.52 5.61 4.52 6.06 6.2 5.99 6.93 4.38 ...

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

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

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

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

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

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

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

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

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

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

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

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

Mais ce n’était pas la vraie question…

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# L’assurance aujourd’hui

Très joli numéro des annales des mines, coordonnée par Pierre-Charles Pradier sur l’assurance aujourd’hui Pierre-Charles a eu la gentillesse de me demander deux courtes contribution dont une avec Laurence Barry et Ewen Gallic sur les probabilités prédictives # Quelle responsabilité pour les algorithmes ?

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

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

# On my way to London, UK Tonight, I will be on my way to Imperial College, in London, to finalize our 2020 actuarial pricing game… More to come soon, on the blog !

(and yes, it’s London, UK, not London, Ontario… where the 6th Workshop on Insurance Mathematics will be organized at the end of the week)