# 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… # Convex Regression Model This morning during the lecture on nonlinear regression, I mentioned (very) briefly the case of convex regression. Since I forgot to mention the codes in R, I will publish them here. Assume that $y_i=m(\mathbf{x}_i)+\varepsilon_i$ where $m:\mathbb{R}^d\rightarrow \mathbb{R}$ is some convex function. Then $m$ is convex if and only if $\forall\mathbf{x}_1,\mathbf{x}_2\in\mathbb{R}^d$, $\forall t\in[0,1]$, $$m(t\mathbf{x}_1+[1-t]\mathbf{x}_2) \leq tm(\mathbf{x}_1)+[1-t]m(\mathbf{x}_2)$$Hidreth (1954) proved that if$$m^\star=\underset{m \text{ convex}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-m(\mathbf{x_i})\big)^2\right\rbrace$$then $\mathbf{\theta}^\star=(m^\star(\mathbf{x_1}),\cdots,m^\star(\mathbf{x_n}))$ is unique. Let $\mathbf{y}=\mathbf{\theta}+\mathbf{\varepsilon}$, then $$\mathbf{\theta}^\star=\underset{\mathbf{\theta}\in \mathcal{K}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-\theta_i)\big)^2\right\rbrace$$where$$\mathcal{K}=\{\mathbf{\theta}\in\mathbb{R}^n:\exists m\text{ convex },m(\mathbf{x}_i)=\theta_i\}$$. I.e. $\mathbf{\theta}^\star$ is the projection of $\mathbf{y}$ onto the (closed) convex cone $\mathcal{K}$. The projection theorem gives existence and unicity. For convenience, in the application, we will consider the real-valued case, $m:\mathbb{R}\rightarrow \mathbb{R}$, i.e. $y_i=m(x_i)+\varepsilon_i$. Assume that observations are ordered $x_1\leq x_2\leq\cdots \leq x_n$. Here $$\mathcal{K}=\left\lbrace\mathbf{\theta}\in\mathbb{R}^n:\frac{\theta_2-\theta_1}{x_2-x_1}\leq \frac{\theta_3-\theta_2}{x_3-x_2}\leq \cdots \leq \frac{\theta_n-\theta_{n-1}}{x_n-x_{n-1}}\right\rbrace$$ Hence, quadratic program with $n-2$ linear constraints. $m^\star$ is a piecewise linear function (interpolation of consecutive pairs $(x_i,\theta_i^\star)$). If $m$ is differentiable, $m$ is convex if $$m(\mathbf{x})+ \nabla m(\mathbf{x})^{\text{T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})$$ More generally, if $m$ is convex, then there exists $\xi_{\mathbf{x}}\in\mathbb{R}^n$ such that $$m(\mathbf{x})+ \xi_{\mathbf{x}}^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})$$ $\xi_{\mathbf{x}}$ is a subgradient of $m$ at ${\mathbf{x}}$. And then $$\partial m(\mathbf{x})=\big\lbrace m(\mathbf{x})+ \xi^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y}),\forall \mathbf{y}\in\mathbb{R}^n\big\rbrace$$ Hence, $\mathbf{\theta}^\star$ is solution of $$\text{argmin}\big\lbrace\|\mathbf{y}-\mathbf{\theta}\|^2\big\rbrace$$$$\text{subject to }\theta_i+\xi_i^{\text{ T}}[\mathbf{x}_j-\mathbf{x}_i]\leq\mathbf{\theta}_j,~\forall i,j$$ and $\xi_1,\cdots,\xi_n\in\mathbb{R}^n$. Now, to do it for real, use cobs package for constrained (b)splines regression, library(cobs) To get a convex regression, use plot(cars) x = cars$speed y = cars$dist rc = conreg(x,y,convex=TRUE) lines(rc, col = 2) Here we can get the values of the knots rc Call: conreg(x = x, y = y, convex = TRUE) Convex regression: From 19 separated x-values, using 5 inner knots, 7, 8, 9, 20, 23. RSS = 1356; R^2 = 0.8766; needed (5,0) iterations and actually, if we use them in a linear-spline regression, we get the same output here reg = lm(dist~bs(speed,degree=1,knots=c(4,7,8,9,,20,23,25)),data=cars) u = seq(4,25,by=.1) v = predict(reg,newdata=data.frame(speed=u)) lines(u,v,col="green") Let us add vertical lines for the knots abline(v=c(4,7,8,9,20,23,25),col="grey",lty=2) # Classification from scratch, boosting 11/8 Eleventh post of our series on classification from scratch. Today, that should be the last one… unless I forgot something important. So today, we discuss boosting. ## An econometrician perspective I might start with a non-conventional introduction. But that’s actually how I understood what boosting was about. And I am quite sure it has to do with my background in econometrics. The goal here is to solve something which looks like$$m^\star=\underset{m\in\mathcal{M}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,m(\mathbf{x}_i))\right\rbrace$$for some loss function $\ell$, and for some set of predictors $\mathcal{M}$. This is an optimization problem. Well, optimization is here in a function space, but still, that’s simply an optimization problem. And from a numerical perspective, optimization is solve using gradient descent (this is why this technique is also called gradient boosting). And the gradient descent can be visualized like below Again, the optimum is not some some real value $x^\star$, but some function $m^\star$. Thus, here we will have something like$$m^{(k)}=m^{(k-1)}+\underset{h\in\mathcal{H}}{\text{argmin}}\left\lbrace \sum_{i=1}^n \ell(y_i,m^{(k-1)}(\mathbf{x}_i)+h(\mathbf{x}_i))\right\rbrace$$(as they write it is serious articles) where the term on the right can also be written$$m^{(k)}=m^{(k-1)}+\underset{h\in\mathcal{H}}{\text{argmin}}\left\lbrace \sum_{i=1}^n \ell(\underbrace{y_i-m^{(k-1)}(\mathbf{x}_i)}_{\varepsilon_{k,i}},h(\mathbf{x}_i))\right\rbrace$$I prefer the later, because we see clearly that $f$ is some model we fit on the remaining residuals. We can rewrite it like that: define$$r_{i,k}=-\left.\frac{\partial \ell(y_i,m(\mathbf{x}_i))}{\partial m(\mathbf{x}_i)}\right\vert_{m(\mathbf{x}_i)=m^{(k-1)}(\mathbf{x}_i)}$$for all $i=1,\cdots,n$. The goal is to fit a model so that $r_{i,k}=h^\star(\mathbf{x}_i)$, and when we have that optimal function, set $m_k(\mathbf{x})=m_{k-1}(\mathbf{x})+\gamma_k h^\star(\mathbf{x})$ (yes, we can include some shrinkage here). Two important comments here. First of all, the idea should be weird to any econometrician. First, we fit a model to explain $y$ by some covariates $\mathbf{x}$. Then consider the residuals $\widehat{\varepsilon}$, and to explain them with the same covariate $\mathbf{x}$. If you try that with a linear regression, you’d done at the end of step 1, since residuals $\widehat{\varepsilon}$ are orthogonal to covariates $\mathbf{x}$: no way that we can learn from them. Here it works because we consider simple non linear model. And actually, something that can be used is to add a shrinkage parameter. Do not consider $\widehat{\varepsilon}=y-\widehat{m}(\mathbf{x})$ but $\widehat{\varepsilon}=y-\gamma\widehat{m}(\mathbf{x})$. The idea of weak learners is extremely important here. The more we shrink, the longer it will take, but that’s not (too) important. I should also mention that it’s nice to keep learning from our mistakes. But somehow, we should stop, someday. I said that I will not mention this part in this series of posts, maybe later on. But heuristically, we should stop when we start to overfit. And this can be observed either using a split training/validation of the initial dataset or to use cross validation. I will get back on that issue later one in this post, but again, those ideas should probably be dedicated to another series of posts. ## Learning with splines Just to make sure we get it, let’s try to learn with splines. Because standard splines have fixed knots, actually, we do not really “learn” here (and after a few iterations we get to what we would have with a standard spline regression). So here, we will (somehow) optimize knots locations. There is a package to do so. And just to illustrate, use a Gaussian regression here, not a classification (we will do that later on). Consider the following dataset (with only one covariate) n=300 set.seed(1) u=sort(runif(n)*2*pi) y=sin(u)+rnorm(n)/4 df=data.frame(x=u,y=y) For an optimal choice of knot locations, we can use library(freeknotsplines) xy.freekt=freelsgen(df$x, df$y, degree = 1, numknot = 2, 555) With 5% shrinkage, the code it simply the following v=.05 library(splines) xy.freekt=freelsgen(df$x, df$y, degree = 1, numknot = 2, 555) fit=lm(y~bs(x,degree=1,knots=xy.freekt@optknot),data=df) yp=predict(fit,newdata=df) df$yr=df$y - v*yp YP=v*yp for(t in 1:200){ xy.freekt=freelsgen(df$x, df$yr, degree = 1, numknot = 2, 555) fit=lm(yr~bs(x,degree=1,knots=xy.freekt@optknot),data=df) yp=predict(fit,newdata=df) df$yr=df$yr - v*yp YP=cbind(YP,v*yp)} nd=data.frame(x=seq(0,2*pi,by=.01)) viz=function(M){ if(M==1) y=YP[,1] if(M&gt;1) y=apply(YP[,1:M],1,sum) plot(df$x,df$y,ylab="",xlab="") lines(df$x,y,type="l",col="red",lwd=3) fit=lm(y~bs(x,degree=1,df=3),data=df) yp=predict(fit,newdata=nd) lines(nd$x,yp,type="l",col="blue",lwd=3) lines(nd$x,sin(nd$x),lty=2)} To visualize the ouput after 100 iterations, use viz(100) Clearly, we see that we learn from the data here… Cool, isn’t it? ## Learning with stumps (and trees) Let us try something else. What if we consider at each step a regression tree, instead of a linear-by-parts regression (that was considered with linear splines). library(rpart) v=.1 fit=rpart(y~x,data=df) yp=predict(fit) df$yr=df$y - v*yp YP=v*yp for(t in 1:100){ fit=rpart(yr~x,data=df) yp=predict(fit,newdata=df) df$yr=df$yr - v*yp YP=cbind(YP,v*yp)} Again, to visualise the learning process, use viz=function(M){ y=apply(YP[,1:M],1,sum) plot(df$x,df$y,ylab="",xlab="") lines(df$x,y,type="s",col="red",lwd=3) fit=rpart(y~x,data=df) yp=predict(fit,newdata=nd) lines(nd$x,yp,type="s",col="blue",lwd=3) lines(nd$x,sin(nd$x),lty=2)} This time, with those trees, it looks like not only we have a good model, but also a different model from the one we can get using a single regression tree. What if we change the shrinkage parameter? viz=function(v=0.05){ fit=rpart(y~x,data=df) yp=predict(fit) df$yr=df$y - v*yp YP=v*yp for(t in 1:100){ fit=rpart(yr~x,data=df) yp=predict(fit,newdata=df) df$yr=df$yr - v*yp YP=cbind(YP,v*yp)} y=apply(YP,1,sum) plot(df$x,df$y,xlab="",ylab="") lines(df$x,y,type="s",col="red",lwd=3) fit=rpart(y~x,data=df) yp=predict(fit,newdata=nd) lines(nd$x,yp,type="s",col="blue",lwd=3) lines(nd$x,sin(nd$x),lty=2)} There is clearly an impact of that shrinkage parameter. It has to be small to get a good model. This is the idea of using weak learners to get a good prediction. ## Classification and Adaboost Now that we understand how bootsting works, let’s try to adapt it to classification. It will be more complicated because residuals are usually not very informative in a classification. And it will be hard to shrink. So let’s try something slightly different, to introduce the adaboost algorithm. In our initial discussion, the goal was to minimize a convex loss function. Here, if we express classes as $\{-1,+1\}$, the loss function we consider is $e^{-y\cdot m(\mathbf{x})}$ (this product $y\cdot m(\mathbf{x})$) was already discussed when we’ve seen the SVM algorithm. Note that the loss function related to the logistic model would be $\log(1+e^{-y\cdot m(\mathbf{x})})$. What we do here is related to gradient descent (or Newton algorithm). Previously, we were learning from our errors. At each iteration, the residuals are computed and a (weak) model is fitted to these residuals. The the contribution of this weak model is used in a gradient descent optimization process. Here things will be different, because (from my understanding) it is more difficult to play with residuals, because null residuals never exist in classifications. So we will add weights. Initially, all the observations will have the same weights. But iteratively, we ill change them. We will increase the weights of the wrongly predicted individuals and decrease the ones of the correctly predicted individuals. Somehow, we want to focus more on the difficult predictions. That’s the trick. And I guess that’s why it performs so well. This algorithm is well described in wikipedia, so we will use it. We start with $\mathbf{\omega}_0=\mathbf{1}/n$, then at each step fit a model (a classification tree) with weights $\mathbf{\omega}_k$(we did not discuss weights in the algorithms of trees, but it is straigtforward in the formula actually). Let $\widehat{h}_{\mathbf{\omega}_k}$ denote that model (i.e. the probability in each leaves). Then consider the classifier $2~\mathbf{1}[\widehat{h}_{\mathbf{\omega}_k}(\cdot)>0.5]-1$ which returns a value in $\{-1,+1\}$. Then set $$\varepsilon_k=\sum_{i\in\mathcal{I}_k}\omega_i$$where $\mathcal{I}_k$ is the set of misclassified individuals,$$\mathcal{I}_k=\big\lbrace i:2~\mathbf{1}[\widehat{h}_{\mathbf{\omega}_k}(\mathbf{x}_i)>0.5]-1\neq y_i\big\rbrace$$Then set $$\alpha_k = \frac{1}{2} \ln \left(\frac{1-\epsilon_k}{\epsilon_k}\right)$$and update finally the model using$$m_{k=1}=m_k+\alpha_k\widehat{h}_{\mathbf{\omega}_k}$$as well as the weights$$\mathbf{\omega}_{k+1}=\mathbf{\omega}_k e^{-\mathbf{y} \alpha_k \widehat{h}_{\mathbf{\omega}_k}(\mathbf{x}_i)}$$(of course, devide by the sum to insure that the total sum is then 1). And as previously, one can include some shrinkage. To visualize the convergence of the process, we will plot the total error on our dataset. n_iter = 100 y = (myocarde[,"PRONO"]==1)*2-1 x = myocarde[,1:7] error = rep(0,n_iter) f = rep(0,length(y)) w = rep(1,length(y)) # alpha = 1 library(rpart) for(i in 1:n_iter){ w = exp(-alpha*y*f) *w w = w/sum(w) rfit = rpart(y~., x, w, method="class") g = -1 + 2*(predict(rfit,x)[,2]&gt;.5) e = sum(w*(y*g&lt;0)) alpha = .5*log ( (1-e) / e ) alpha = 0.1*alpha f = f + alpha*g error[i] = mean(1*f*y&lt;0) } plot(seq(1,n_iter),error,type=&quot;l&quot;, ylim=c(0,.25),col=&quot;blue&quot;, ylab=&quot;Error Rate&quot;,xlab=&quot;Iterations&quot;,lwd=2) Here we face a classical problem in machine learning: we have a perfect model. With zero error. That is nice, but not interesting. It is also possible in econometrics, with polynomial fits: with 10 observations, and a polynomial of degree 9, we have a perfect fit. But a poor model. Here it is the same. So the trick is to split our dataset in two, a training dataset, and a validation one set.seed(123) id_train = sample(1:nrow(myocarde), size=45, replace=FALSE) train_myocarde = myocarde[id_train,] test_myocarde = myocarde[-id_train,] We construct the model on the first one, and we check on the second one that it’s not that bad… y_train = (train_myocarde[,"PRONO"]==1)*2-1 x_train = train_myocarde[,1:7] y_test = (test_myocarde[,"PRONO"]==1)*2-1 x_test = test_myocarde[,1:7] train_error = rep(0,n_iter) test_error = rep(0,n_iter) f_train = rep(0,length(y_train)) f_test = rep(0,length(y_test)) w_train = rep(1,length(y_train)) alpha = 1 for(i in 1:n_iter){ w_train = w_train*exp(-alpha*y_train*f_train) w_train = w_train/sum(w_train) rfit = rpart(y_train~., x_train, w_train, method="class") g_train = -1 + 2*(predict(rfit,x_train)[,2]&gt;.5) g_test = -1 + 2*(predict(rfit,x_test)[,2]&gt;.5) e_train = sum(w_train*(y_train*g_train&lt;0)) alpha = .5*log ( (1-e_train) / e_train ) alpha = 0.1*alpha f_train = f_train + alpha*g_train f_test = f_test + alpha*g_test train_error[i] = mean(1*f_train*y_train&lt;0) test_error[i] = mean(1*f_test*y_test&lt;0)} plot(seq(1,n_iter),test_error,col='red') lines(train_error,lwd=2,col='blue') Here, as previously, after 80 iterations, we have a perfect model on the training dataset, but it behaves badly on the validation dataset. But with 20 iterations, it seems to be ok… ## R function Of course, it’s possible to use R functions, library(gbm) gbmWithCrossValidation = gbm(PRONO ~ .,distribution = "bernoulli", data = myocarde,n.trees = 2000,shrinkage = .01,cv.folds = 5,n.cores = 1) bestTreeForPrediction = gbm.perf(gbmWithCrossValidation) Here cross-validation is considered, and not training/validation, as well as forests instead of single trees, but overall, the idea is the same… Off course, the output is much nicer (here the shrinkage is a very small parameter, and learning is extremely slow) # Classification from scratch, logistic with splines 2/8 Today, second post of our series on classification from scratch, following the brief introduction on the logistic regression. ## Piecewise linear splines To illustrate what’s going on, let us start with a “simple” regression (with only one explanatory variable). The underlying idea is natura non facit saltus, for “nature does not make jumps”, i.e. process governing equations for natural things are continuous. That seems to be a rather strong assumption, because we can assume that there is a fixed threshold to explain death. For instance, if patients die (for sure) if the “stroke index” exceeds a threshold, we might expect some discontinuity. Exceept that if that threshold is an heterogeneous (non-observable continuous) variable, then we get back to the continuity assumption. The most simple model we can think of to extend the linear model we’ve seen in the previous post is to consider a piecewise linear function, with two parts : small values of $x$, and larger values of $x$. The most convenient way to do so is to use the positive part function $(x-s)_+$ which is the difference between $x$ and $s$ if that difference is positive, and $0$ otherwise. For instance $$\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. And of course, it is possible to consider more than one knot. The function to get the positive value is the following pos = function(x,s) (x-s)*(x&gt;=s) then we can use it direcly in our regression model reg = glm(PRONO~INSYS+pos(INSYS,15)+ pos(INSYS,25),data=myocarde,family=binomial) The output of the regression is here summary(reg) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.1109 3.2783 -0.034 0.9730 INSYS -0.1751 0.2526 -0.693 0.4883 pos(INSYS, 15) 0.7900 0.3745 2.109 0.0349 * pos(INSYS, 25) -0.5797 0.2903 -1.997 0.0458 * Hence, the original slope, for very small values is not significant, but then, above 15, it become significantly positive. And above 25, there is a significant change again. We can plot it to see what’s going on u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,type="l") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2) ## Using bs() linear splines Using the GAM function, things are slightly different. We will use here so called b-splines, library(splines) We can define spline functions with support $(5,55)$ and with knots $\{15,25\}$ clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02") x = seq(0,60,by=.25) B = bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=1) matplot(x,B,type="l",lty=1,lwd=2,col=clr6) as we can see, the functions defined here are different from the one before, but we still have (piecewise) linear functions on each segment $(5,15)$, $(15,25)$ and $(25,55)$. But linear combinations of those functions (the two sets of functions) will generate the same space. Said differently, if the interpretation of the output will be different, predictions should be the same reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=1), data=myocarde,family=binomial) summary(reg) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.9863 2.0555 -0.480 0.6314 bs(INSYS,..)1 -1.7507 2.5262 -0.693 0.4883 bs(INSYS,..)2 4.3989 2.0619 2.133 0.0329 * bs(INSYS,..)3 5.4572 5.4146 1.008 0.3135 Observe that there are three coefficients, as before, but again, the interpretation is here more complicated… v=predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2) Nevertheless, the prediction is the same… and that’s nice. ## Piecewise quadratic splines Let us go one step further… Can we have also the continuity of the derivative ? Yes, and that’s easy actually, considering parabolic functions. Instead of using a decomposition on $x,(x-s_1)_+$ and $(x-s_2)_+$ consider now a decomposition on $x,x^{\color{red}{2}},(x-s_1)^{\color{red}{2}}_+$ and $(x-s_2)^{\color{red}{2}}_+$.  pos2 = function(x,s) (x-s)^2*(x&gt;=s) reg = glm(PRONO~poly(INSYS,2)+pos2(INSYS,15)+pos2(INSYS,25), data=myocarde,family=binomial) summary(reg) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 29.9842 15.2368 1.968 0.0491 * poly(INSYS, 2)1 408.7851 202.4194 2.019 0.0434 * poly(INSYS, 2)2 199.1628 101.5892 1.960 0.0499 * pos2(INSYS, 15) -0.2281 0.1264 -1.805 0.0712 . pos2(INSYS, 25) 0.0439 0.0805 0.545 0.5855 As expected, there are here five coefficients: the intercept and two for the part on the left (three parameters for the parabolic function), and then two additional terms for the part in the center – here $(15,25)$ – and for the part on the right. Of course, for each portion, there is only one degree of freedom since we have a parabolic function (three coefficients) but two constraints (continuity, and continuity of the first order derivative). On a graph, we get the following v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2,xlab="INSYS",ylab="") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2) ## Using bs() quadratic splines Of course, we can do the same with our R function. But as before, the basis of function is expressed here differently  x = seq(0,60,by=.25) B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=2) matplot(x,B,type="l",xlab="INSYS",col=clr6) If we run R code, we get reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=2),data=myocarde, family=binomial) summary(reg) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 7.186 5.261 1.366 0.1720 bs(INSYS, ..)1 -14.656 7.923 -1.850 0.0643 . bs(INSYS, ..)2 -5.692 4.638 -1.227 0.2198 bs(INSYS, ..)3 -2.454 8.780 -0.279 0.7799 bs(INSYS, ..)4 6.429 41.675 0.154 0.8774 But that’s not really a big deal since the prediction is exactly the same v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2) ## Cubic splines Last, but not least, we can reach the cubic splines. With our previous notions, we would consider a decomposition on (guess what) $x,x^2,x^{\color{red}{3}},(x-s_1)^{\color{red}{3}}_+,(x-s_2)^{\color{red}{3}}_+$, to get this time continuity, as well as continuity of the first two derivatives (and to get a very smooth function, since even variations will be smooth). If we use the bs function, the basis is the followin B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=3) matplot(x,B,type="l",lwd=2,col=clr6,lty=1,ylim=c(-.2,1.2)) abline(v=c(5,15,25,55),lty=2) and the prediction will now be reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=3), data=myocarde,family=binomial) u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2) Two last things before concluding (for today), the location of the knots, and the extension to additive models. ## Location of knots In many applications, we do not want to specify the location of the knots. We just want – say – three (intermediary) knots. This can be done using reg = glm(PRONO~1+bs(INSYS,degree=1,df=4),data=myocarde,family=binomial) We can actually get the locations of the knots by looking at attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 1L, knots = c(15.8, 21.4, 27.15), Boundary.knots = c(8.7, 54), intercept = FALSE)

which provides us with the location of the boundary knots (the minumun and the maximum from from our sample) but also the three intermediary knots. Observe that actually, those five values are just (empirical) quantiles

quantile(myocarde$INSYS,(0:4)/4) 0% 25% 50% 75% 100% 8.70 15.80 21.40 27.15 54.00 If we plot the prediction, we get v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2)

If we get back on what was computed before the logit transformation, we clealy see ruptures are the different quantiles

B = bs(x,degree=1,df=4) B = cbind(1,B) y = B%*%coefficients(reg) plot(x,y,type="l",col="red",lwd=2) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2) Note that if we do specify anything about knots (number or location), we get no knots… reg = glm(PRONO~1+bs(INSYS,degree=2),data=myocarde,family=binomial) attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 2L, knots = numeric(0), Boundary.knots = c(8.7,54), intercept = FALSE)

and if we look at the prediction

u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19)

actually, it is the same as a quadratic regression (as expected actually)

reg = glm(PRONO~1+poly(INSYS,degree=2),data=myocarde,family=binomial) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19)

Consider now the second dataset, with two variables. Consider here a model like
$$\mathbb{P}[Y|X_1=x_1,X_2=x_2]=\frac{\exp[\eta(x_1,x_2)]}{1+\exp[\eta(x_1,x_2)]}$$
where
$$\exp[\eta(x_1,x_2)]=\beta_0+\color{red}{s_1(x_1)}+\color{blue}{s_2(x_2)}$$
$$\color{red}{s_1(x_1)}=\beta_{1,0}x_1+\beta_{1,1}(x_1-s_{11})_++\beta_{1,2}(x_1-s_{12})_+$$
and
$$\color{blue}{s_2(x_2)}=\beta_{2,0}x_2+\beta_{2,1}(x_2-s_{21})_++\beta_{2,2}(x_2-s_{22})_+$$
It might seem a little bit restrictive, but that’s actually the idea of additive models.

reg = glm(y~bs(x1,degree=1,df=3)+bs(x2,degree=1,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) Now, if think about is, we’ve been able to get a “perfect” model, so, somehow, it seems no longer continuous… persp(u,u,v,theta=20,phi=40,col="green" Of course, it is… it is piecewise linear, with hyperplane, some being almost vertical. And one can also consider piecewise quadratic functions reg = glm(y~bs(x1,degree=2,df=3)+bs(x2,degree=2,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Funny thing, we now have two “perfect” models, with different areas for the white and the black dots… Don’t ask me how to choose on that one.

In R, it is possible to use the mgcv package to run a gam regression. It is used for generalized additive models, but here, we have only one variable, so it is difficult to see the “additive” part, actually. And to be more specific, mgcv is using penalized quasi-likelihood from the nlme package (but we’ll get back on penalized routines later on).

But maybe I should also mention another smoothing tool before, kernels (and maybe also $k$-nearest neighbors). To be continued

On Thursday, March 2nd, I will give the first lecture of the PhD course on advanced tools for econometrics, on nonlinearities. Slides are available online.

# Data Science for Actuaries, Regression Models with R

After an introduction to Advanced R, we will discuss for the last part of our crash course visualization and graphs (from the previous set of slides), and I just uploaded additional slides on regression models (including some pdf version)

# Regression with Splines: Should we care about Non-Significant Components?

Following the course of this morning, I got a very interesting question from a student of mine. The question was about having non-significant components in a splineregression.  Should we consider a model with a small number of knots and all components significant, or one with a (much) larger number of knots, and a lot of knots non-significant?

My initial intuition was to prefer the second alternative, like in autoregressive models in R. When we fit an AR(6) model, it’s not really a big deal if most coefficients are not significant (but the last one). It’s won’t affect much the forecast. So here, it might be the same. With a larger number of knots, we should be able to capture small bumps that we’ll never capture with a smaller number.

Here is what a have with a small number of knots, and cubic splines

and with a larger number of knots

In order to understand what’s going on, consider a simple model, with the two splines above, in red

> set.seed(1)
> library(splines)
> x=seq(0,1,by=.01)
> v=bs(x,10)
> x2=v[,2]
> x10=v[,10]
> set.seed(1)
> y=1+3*x2+5*x10+rnorm(length(x))/4
> y_test=1+3*x2+5*x10+rnorm(length(x))/4

Note that here I have generated two sets of data, one to train a model, and one to test it.  Here, the data looks like that

> plot(x,y)

It is based on two splines,

> lines(df$x,1+3*x2+5*x10) If we use a spline model with 10 degrees of freedom, we get > df=data.frame(x,y) > reg=lm(y~bs(x,10),data=df) > summary(reg) Coefficients: Estimate Std. Er t value Pr(>|t|) (Intercept) 0.91671 0.17068 5.371 6.08e-07 *** bs(x, 10)1 0.20485 0.32696 0.627 0.533 bs(x, 10)2 3.15593 0.22534 14.005 < 2e-16 *** bs(x, 10)3 0.04847 0.25075 0.193 0.847 bs(x, 10)4 0.09373 0.21597 0.434 0.665 bs(x, 10)5 0.11624 0.22939 0.507 0.614 bs(x, 10)6 0.24829 0.22293 1.114 0.268 bs(x, 10)7 -0.06825 0.23498 -0.290 0.772 bs(x, 10)8 0.19633 0.26241 0.748 0.456 bs(x, 10)9 0.27557 0.26976 1.022 0.310 bs(x, 10)10 4.78134 0.24116 19.826 < 2e-16 *** which makes sense, from what we have generated. Indeed, most of the components are not significant, but the second and the tenth. We can actually test that all those components are null (at the same time) > A=matrix(0,8,11) > colnames(A)=names(coefficients(reg)) > A[1,2]=A[2,4]=A[3,5]=A[4,6]=A[5,7]= + A[6,8]=A[7,9]=A[8,10]=1 > b=rep(0,8) > linearHypothesis(reg, A,b) Linear hypothesis test Hypothesis: bs(x, 10)1 = 0 bs(x, 10)3 = 0 bs(x, 10)4 = 0 bs(x, 10)5 = 0 bs(x, 10)6 = 0 bs(x, 10)7 = 0 bs(x, 10)8 = 0 bs(x, 10)9 = 0 Model 1: restricted model Model 2: y ~ bs(x, 10) Res.Df RSS Df Sum of Sq F Pr(>F) 1 98 4.8766 2 90 4.6196 8 0.25701 0.6259 0.754 and yes, those coefficients are not significant. > yp10=predict(reg) > lines(df$x,yp10,col="red")

# Computing AIC on a Validation Sample

This afternoon, we’ve seen in the training on data science that it was possible to use AIC criteria for model selection.

> library(splines)
> AIC(glm(dist ~ speed, data=train_cars,
[1] 438.6314
> AIC(glm(dist ~ speed, data=train_cars,
[1] 436.3997
> AIC(glm(dist ~ bs(speed), data=train_cars,
[1] 425.6434
> AIC(glm(dist ~ bs(speed), data=train_cars,
[1] 428.7195

And I’ve been asked why we don’t use a training sample to fit a model, and then use a validation sample to compare predictive properties of those models, penalizing by the complexity of the model.    But it turns out that it is difficult to compute the AIC of those models on a different dataset. I mean, it is possible to write down the likelihood (since we have a Poisson model) but I want a code that could work for any model, any distribution….

Hopefully, Heather suggested a very clever idea, using her package

And actually, it works well.

# Choosing a Classifier

In order to illustrate the problem of chosing a classification model consider some simulated data,

> n = 500
> set.seed(1)
> X = rnorm(n)
> ma = 10-(X+1.5)^2*2
> mb = -10+(X-1.5)^2*2
> M = cbind(ma,mb)
> set.seed(1)
> Z = sample(1:2,size=n,replace=TRUE)
> Y = ma*(Z==1)+mb*(Z==2)+rnorm(n)*5
> df = data.frame(Z=as.factor(Z),X,Y)

A first strategy is to split the dataset in two parts, a training dataset, and a testing dataset.

> df1 = training = df[1:300,]
> df2 = testing  = df[301:500,]
• The Holdout Method: Training and Testing Datasets

The two datasets can be visualised below, with the training dataset on top, and the testing dataset below

> plot(df1$X,df1$Y,pch=19,col=c(rgb(1,0,0,.4),
+ rgb(0,0,1,.4))[df1$Z]) # An Update on Boosting with Splines In my previous post, An Attempt to Understand Boosting Algorithm(s), I was puzzled by the boosting convergence when I was using some spline functions (more specifically linear by parts and continuous regression functions). I was using > library(splines) > fit=lm(y~bs(x,degree=1,df=3),data=df) The problem with that spline function is that knots seem to be fixed. The iterative boosting algorithm is • start with some regression model $\boldsymbol{y}_1=h_1(\boldsymbol{x})$ • compute the residuals, including some shrinkage parameter,$\boldsymbol{\varepsilon}_{1}=\boldsymbol{y}-\nu_1 h_1(\boldsymbol{x})$ then the strategy is to model those residuals • at step $j$, consider regression $\boldsymbol{\varepsilon}_j=h_j(\boldsymbol{x})$ • update the residuals $\boldsymbol{\varepsilon}_{j+1}=\boldsymbol{\varepsilon}_j-\nu_j h_j(\boldsymbol{x})$ and to loop. Then set $\widehat{\boldsymbol{y}}=\sum_{j=1}^M \nu_j\boldsymbol{\varepsilon}_{j}=\sum_{j=1}^M \nu_jh_j(\boldsymbol{x})$ I thought that boosting would work well if at step $j$, it was possible to change the knots. But the output was quite disappointing: boosting does not improve the prediction here. And it looks like knots don’t change. Actually, if we select the ‘best‘ knots, the output is much better. The dataset is still > n=300 > set.seed(1) > u=sort(runif(n)*2*pi) > y=sin(u)+rnorm(n)/4 > df=data.frame(x=u,y=y) For an optimal choice of knot locations, we can use > library(freeknotsplines) > xy.freekt=freelsgen(df$x, df$y, degree = 1, + numknot = 2, 555) The code of the previous post can simply be updated > v=.05 > library(splines) > xy.freekt=freelsgen(df$x, df$y, degree = 1, + numknot = 2, 555) > fit=lm(y~bs(x,degree=1,knots= + xy.freekt@optknot),data=df) > yp=predict(fit,newdata=df) > df$yr=df$y - v*yp > YP=v*yp > for(t in 1:200){ + xy.freekt=freelsgen(df$x, df$yr, degree = 1, + numknot = 2, 555) + fit=lm(yr~bs(x,degree=1,knots= + xy.freekt@optknot),data=df) + yp=predict(fit,newdata=df) + df$yr=df$yr - v*yp + YP=cbind(YP,v*yp) + } > nd=data.frame(x=seq(0,2*pi,by=.01)) > viz=function(M){ + if(M==1) y=YP[,1] + if(M>1) y=apply(YP[,1:M],1,sum) + plot(df$x,df$y,ylab="",xlab="") + lines(df$x,y,type="l",col="red",lwd=3)
+    fit=lm(y~bs(x,degree=1,df=3),data=df)
+    yp=predict(fit,newdata=nd)
+    lines(nd$x,yp,type="l",col="blue",lwd=3) + lines(nd$x,sin(nd$x),lty=2)} > viz(100) I like that graph. I had the intuition that using (simple) splines would be possible, and indeed, we get a very smooth prediction. # I Fought the (distribution) Law (and the Law did not win) A few days ago, I was asked if we should spend a lot of time to choose the distribution we use, in GLMs, for (actuarial) ratemaking. On that topic, I usually claim that the family is not the most important parameter in the regression model. Consider the following dataset > db <- data.frame(x=c(1,2,3,4,5),y=c(1,2,4,2,6)) > plot(db,xlim=c(0,6),ylim=c(-1,8),pch=19) To visualize a regression model, use the following code > nd=data.frame(x=seq(0,6,by=.1)) > add_predict = function(reg){ + prd1=predict(reg,newdata=nd,se.fit = TRUE,type="response") + y1=prd1$fit
+ y1_upp=prd1$fit+prd1$residual.scale*1.96*
prd1$se.fit + y1_low=prd1$fit-prd1$residual.scale*1.96* prd1$se.fit
+ polygon(c(nd$x,rev(nd$x)),c(y1_upp,
rev(y1_low)),col="light green",angle=90,
density=40,border=NA)
+ lines(nd$x,y1,col="red",lwd=2) + } For instance, with a Poisson regression (with a log link function) we get > plot(db) > reg1=glm(y~x,family=poisson(link="log"), + data=db) > add_predict(reg1)  while, with a Gaussian regresion (but still with a log link function), we get > plot(db) > reg2=glm(y~x,family=gaussian(link="log"), + data=db) > add_predict(reg2) If we just care about the expected value of our prediction, the output is more or less the same > plot(db) > lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd, + type="response"),col="blue",lwd=1.5) So, indeed, forget about the (distribution) law when running a GLM. Not convinced? Consider – on the same dataset – a Poisson regression (with an identity link function this time) > plot(db) > reg1=glm(y~x,family=poisson(link="identity"), + data=db) > add_predict(reg1)  while, with a Gaussian regresion (but still with an identity link function), we get > plot(db) > reg2=glm(y~x,family=gaussian(link="identity"), + data=db) > add_predict(reg2) Again, if we just plot the expected value of our prediction, the output is more or less the same > plot(db) > lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd, + type="response"),col="blue",lwd=1.5) So clearly, the simplistic message you should not care too much about the (distribution) law seems to be valid… # Some heuristics about spline smoothing Let us continue our discussion on smoothing techniques in regression. Assume that .$\mathbb{E}(Y\vert X=x)=h(x)$ where $h(\cdot)$ is some unkown function, but assumed to be sufficently smooth. For instance, assume that $h(\cdot)$ is continuous, that $h'(\cdot)$ exists, and is continuous, that $h''(\cdot)$ exists and is also continuous, etc. If $h(\cdot)$ is smooth enough, Taylor’s expansion can be used. Hence, for $x\in(\alpha,\beta)$ $h(x)=h(\alpha)+\sum_{k=1}^ d \frac{(x-\alpha)^k}{k!}h^{(k)}(x_0)+\frac{1}{d!}\int_{\alpha}^x [x-t]^d h^{(d+1)}(t)dt$ which can also be writen as $h(x)=\sum_{k=0}^ d a_k (x-\alpha)^k +\frac{1}{d!}\int_{\alpha}^x [x-t]^d h^{(d+1)}(t)dt$ for some $a_k$‘s. The first part is simply a polynomial. The second part, is some integral. Using Riemann integral, observe that $\frac{1}{d!}\int_{\alpha}^x [x-t]^d h^{(d+1)}(t)dt\sim \sum_{i=1}^ j b_i (x-x_i)_+^d$ for some $b_i$‘s, and some $\alpha < x_1< x_2< \cdots < x_{j-1} < x_j < \beta$ Thus, $h(x) \sim \sum_{k=0}^ d a_k (x-\alpha)^k +\sum_{i=1}^ j b_i (x-x_i)_+^d$ Nice! We have our linear regression model. A natural idea is then to consider a regression of $Y$ on $\boldsymbol{X}$ where $\boldsymbol{X} = (1,X,X^2,\cdots,X^d,(X-x_1)_+^d,\cdots,(X-x_k)_+^d )$ given some knots $\{x_1,\cdots,x_k\}$. To make things easier to understand, let us work with our previous dataset, plot(db) If we consider one knot, and an expansion of order 1, attach(db) library(splines) B=bs(xr,knots=c(3),Boundary.knots=c(0,10),degre=1) reg=lm(yr~B) lines(xr[xr<=3],predict(reg)[xr<=3],col="red") lines(xr[xr>=3],predict(reg)[xr>=3],col="blue") The prediction obtained with this spline can be compared with regressions on subsets (the doted lines) reg=lm(yr~xr,subset=xr<=3) lines(xr[xr<=3],predict(reg)[xr<=3],col="red",lty=2) reg=lm(yr~xr,subset=xr>=3) lines(xr[xr>=3],predict(reg),col="blue",lty=2) It is different, since we have here three parameters (and not four, as for the regressions on the two subsets). One degree of freedom is lost, when asking for a continuous model. Observe that it is possible to write, equivalently reg=lm(yr~bs(xr,knots=c(3),Boundary.knots=c(0,10),degre=1),data=db) So, what happened here? B=bs(xr,knots=c(2,5),Boundary.knots=c(0,10),degre=1) matplot(xr,B,type="l") abline(v=c(0,2,5,10),lty=2) Here, the functions that appear in the regression are the following Now, if we run the regression on those two components, we get B=bs(xr,knots=c(2,5),Boundary.knots=c(0,10),degre=1) matplot(xr,B,type="l") abline(v=c(0,2,5,10),lty=2) If we add one knot, we get the prediction is reg=lm(yr~B) lines(xr,predict(reg),col="red") Of course, we can choose much more knots, B=bs(xr,knots=1:9,Boundary.knots=c(0,10),degre=1) reg=lm(yr~B) lines(xr,predict(reg),col="red") We can even get a confidence interval reg=lm(yr~B) P=predict(reg,interval="confidence") plot(db,col="white") polygon(c(xr,rev(xr)),c(P[,2],rev(P[,3])),col="light blue",border=NA) points(db) reg=lm(yr~B) lines(xr,P[,1],col="red") abline(v=c(0,2,5,10),lty=2) And if we keep the two knots we chose previously, but consider Taylor’s expansion of order 2, we get B=bs(xr,knots=c(2,5),Boundary.knots=c(0,10),degre=2) matplot(xr,B,type="l") abline(v=c(0,2,5,10),lty=2) So, what’s going on? If we consider the constant, and the first component of the spline based matrix, we get k=2 plot(db) B=cbind(1,B) lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1) If we add the constant term, the first term and the second term, we get the part on the left, before the first knot, k=3 lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1) and with three terms from the spline based matrix, we can get the part between the two knots, k=4 lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1) and finallty, when we sum all the terms, we get this time the part on the right, after the last knot, k=5 lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1) This is what we get using a spline regression, quadratic, with two (fixed) knots. And can can even get confidence intervals, as before reg=lm(yr~B) P=predict(reg,interval="confidence") plot(db,col="white") polygon(c(xr,rev(xr)),c(P[,2],rev(P[,3])),col="light blue",border=NA) points(db) reg=lm(yr~B) lines(xr,P[,1],col="red") abline(v=c(0,2,5,10),lty=2) The great idea here is to use functions $(x-x_i)_+$, that will insure continuity at point $x_i$. Of course, we can use those splines on our Dexter application, Here again, using linear spline function, it is possible to impose a continuity constraint, plot(data$no,data$mu,ylim=c(6,10)) abline(v=12*(0:8)+.5,lty=2) reg=lm(mu~bs(no,knots=c(12*(1:7)+.5),Boundary.knots=c(0,97), degre=1),data=db) lines(c(1:94,96),predict(reg),col="red") But we can also consider some quadratic splines, plot(data$no,data$mu,ylim=c(6,10)) abline(v=12*(0:8)+.5,lty=2) reg=lm(mu~bs(no,knots=c(12*(1:7)+.5),Boundary.knots=c(0,97), degre=2),data=db) lines(c(1:94,96),predict(reg),col="red") # Multiple (smoothed) regression and portfolio exposure Wednesday, in class, we’ve seen how to visualize a multiple regression model (with two continuous explanatory variables). Here, the goal is to predict the average cost of an insurance claim, using some covariates, e.g. the age of the driver, and the age of the car (recall that losses here are liability losses). The prediction obtained from a (standard) generalized linear model, with a log-link > reg1=glm(cout~ageconducteur+agevehicule,data=base,family=Gamma(link="log")) The code to visualize the predicted average cost is the following: first, we have to compute predictions for specific values, > pred=function(x,y){ + predict(reg,newdata=data.frame(ageconducteur=x, + agevehicule=y),type="response") Then, we use this function to compute values on a grid, > X=seq(20,80,by=5) > Y=0:20 > Z=outer(X,Y,p) > image(X,Y,Z,col=rev(heat.colors(101))) > contour(X,Y,Z,add=TRUE, + levels=c(1400,1800,2000,2200,2400,2600,2800,3000,3200,4000,5000)) If we use factors, and not continuous variates (cut versions of those two variates), > reg2=glm(cout~cut(ageconducteur,breaks=c(0,22,35,55,80,100))* + cut(agevehicule,breaks=c(-1,1,3,5,10,100)), + data=base,family=Gamma(link="log")) (note that we consider the Cartesian product, so values are computed for each product of factors, age of the driver and age of the car) we obtain Obviously, we’re missing something here: the most expensive class with one model is the cheapeast for the other one! Of course, it might come from our classes (that were chosen a bit randomly), but it might be interesting to use nonlinear functions of the ages. So, let us use splines to smooth those two variables, > reg3=glm(cout~bs(ageconducteur)+bs(agevehicule),data=base, + family=Gamma(link="log")) With additive smoothed functions, we obtained a symmetric graph (due to the additive property) while with a bivariate spline > library(mgcv) + reg4=gam(cout~s(ageconducteur,agevehicule),data=base, + family=Gamma(link="log")) (for some odd reasons, I could not use – easily – a bivariate spline in the Generalized Linear Model, but it did work considering a Generalized Additive Model – which is, by no means additive now). We can identify here some regions where the average cost can be extremely expensive… But, as mentioned wednesday, one should keep in mind that some parts of the square above are not reached. More precisely, the distribution of the portfolio, as a function of those two covariates is the following Thus, the proportion of young drivers driving a brand new car, and the proportion of old drivers driving a very old car is rather small… If the goal is to find niches, one should look at the prediction more carefully, but if the goal is to make that everyone gets an insurance cover, maybe we should allow that some drivers are under-priced (especially when are rare in the portfolio). And one should keep in mind that average costs are extremely sensitive to large losses, as discussed previously http://freakonometrics.hypotheses.org/3490 (and in class) In the univariate case, I have migrated an old post, we I tried to reproduce (in R and in French) some standard graphs in the insurance industry: it is always interesting to visualize not only the prediction obtained from our models, but also the size of each class in the portfolio, The post is online here http://freakonometrics.hypotheses.org/1224 # Natura non facit saltus (see John Wilkins’ article on the – interesting – history of that phrase http://scienceblogs.com/evolvingthoughts/…). We will see, this week in class, several smoothing techniques, for insurance ratemaking. As a starting point, assume that we do not want to use segmentation techniques: everyone will pay exactly the same price. • no segmentation of the premium And that price should be related to the pure premium, which is proportional to the frequency (or the annualized frequency, as discussed previously), since $\mathbb{E}_{\mathbb{P}}\left(\sum_{i=1}^N Y_i\right)=\mathbb{E}_{\mathbb{P}}(N) \cdot \mathbb{E}_{\mathbb{P}}(Y_i)$ The probability measure is mentioned here just to recall that we can use any measure. Even $\mathbb{P}_{\boldsymbol{X}}$ (based on some covariates). Without any covariate, the expected frequency should be > regglm0=glm(nbre~1+offset(log(exposition)),data=sinistres,family=poisson) > summary(regglm0) Call: glm(formula = nbre ~ 1 + offset(log(exposition)), family = poisson, data = sinistres) Deviance Residuals: Min 1Q Median 3Q Max -0.5033 -0.3719 -0.2588 -0.1376 13.2700 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.6201 0.0228 -114.9 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 12680 on 49999 degrees of freedom Residual deviance: 12680 on 49999 degrees of freedom AIC: 16353 Number of Fisher Scoring iterations: 6 > exp(coefficients(regglm0)) (Intercept) 0.07279295 Thus, if we do not want to take into account potential heterogeneity, we should assume that $N\sim\mathcal{P}(\lambda)$ where $\lambda$ is closed to 7.28%. Yes, as mentioned in class, it is rather common to see $\lambda$ as a percentage, i.e. a probability, since $\mathbb{P}(N\neq 0)=1-e^{-\lambda}\approx \lambda$ i.e. $\lambda$ can be interpreted as the probability of not have a claim (see also the law of small numbers). Let us visualize this as a function of the age of the driver, > a=18:100 > yp=predict(regglm0,newdata=data.frame(ageconducteur=a,exposition=1),type="response",se.fit=TRUE) > yp0=yp$fit
> yp1=yp$fit+2*yp$se.fit
> yp2=yp$fit-2*yp$se.fit
> plot(a,yp0,type="l",ylim=c(.03,.12))
> abline(v=40,col="grey")
> lines(a,yp1,lty=2)
> lines(a,yp2,lty=2)
> k=23
> points(a[k],yp0[k],pch=3,lwd=3,col="red")
> segments(a[k],yp1[k],a[k],yp2[k],col="red",lwd=3)

We do predict the same frequency for all drivers, e.g. for some drive aged 40,

> cat("Frequency =",yp0[k]," confidence interval",yp1[k],yp2[k])
Frequency = 0.07279295  confidence interval 0.07611196 0.06947393

Let us now consider the case where we try to take into account heterogeneity, e.g. by age,

• The (standard) Poisson regression

The idea of the (log-)Poisson regression is to assume that instead of having $N\sim\mathcal{P}(\lambda)$, we should have $N|\boldsymbol{X}\sim\mathcal{P}(\lambda_{\boldsymbol{X}})$, where

$\lambda_{\boldsymbol{X}}=\exp(\beta_0+\beta_1 \boldsymbol{X}_1+\cdots+\beta_k\boldsymbol{X}_k)$

in a very general setting. Here, let us consider only one explanatory variable, i.e.

$\lambda_{X}=\exp(\beta_0+\beta_1 {X})$

Here, we have

> yp=predict(regglm1,newdata=data.frame(ageconducteur=a,exposition=1),
+ type="response",se.fit=TRUE)
> yp0=yp$fit > yp1=yp$fit+2*yp$se.fit > yp2=yp$fit-2*yp$se.fit > plot(a,yp0,type="l",ylim=c(.03,.12)) > abline(v=40,col="grey") > lines(a,yp1,lty=2) > lines(a,yp2,lty=2) > points(a[k],yp0[k],pch=3,lwd=3,col="red") > segments(a[k],yp1[k],a[k],yp2[k],col="red",lwd=3) i.e. the prediction for the annualized claim frequency for our 40 year old driver is now 7.74% (which is slightly higher than what we had before, 7.28%) > cat("Frequency =",yp0[k]," confidence interval",yp1[k],yp2[k]) Frequency = 0.07740574 confidence interval 0.08117512 0.07363636 It is possible to compute not the expected frequency , but the ratio $\mathbb{E}(N|X)/\mathbb{E}(N)$. Above the horizontal blue line, the premium will be higher than the one obtained without segmentation, and (of course) lower below. Here, drivers younger than 44 year old will pay more, while driver older than 44 year old will be less. We have discussed, in the introduction, the necessity of segmentation. If we consider two companies, one segmenting, while the other one has a flat rate, then older drivers will go to the first company (since insurance is cheaper) while younger ones will go to the second one (again, it is cheaper). The problem is that the second company implicitly hopes that older drivers will compensate the risk. But since they’re gone, insurance will be too cheap, and the company will loose money (if not goes bankrupt). So companies have to use segmentation techniques to survive. Now, the problem is that we cannot be sure that this exponential decay of the premium is the proper way the premium should evolve as a function of the age. An alternative can be to use nonparametric techniques to visualize to true influence of the age on claims frequency. • A pure nonparametric model A first model can be to consider a premium, per age. This can be done considering the age of the driver as a factor in the regression, > regglm2=glm(nbre~as.factor(ageconducteur)+offset(log(exposition)), + data=sinistres,family=poisson) > yp=predict(regglm2,newdata=data.frame(ageconducteur=a0,exposition=1), + type="response",se.fit=TRUE) > yp0=yp$fit
> yp1=yp$fit+2*yp$se.fit
> yp2=yp$fit-2*yp$se.fit
> plot(a0,yp0,type="l",ylim=c(.03,.12))
> abline(v=40,col="grey")

Here, the forecast for our 40 year old driver is slightly lower than be previous one, but the confidence interval is much larger (since we focus on a very small subclass of the portfolio: drivers aged exactly 40)

Frequency = 0.06686658  confidence interval 0.08750205 0.0462311

Here, we consider too small classes, and the premium is too erratic: the premium will decrease of 20% from age 40 to 41, and then increase of 50% from age 41 to 42,

> diff(log(yp0[23:25]))
24         25
-0.2330241  0.5223478

There is no chance that the company will keep the insured with this strategy. This discontinuity of the premium is clearly an important issue here.

• Using age classes

An alternative can be to consider age classes, from very young drivers to senior drivers.

> level1=seq(15,105,by=5)
> regglmc1=glm(nbre~cut(ageconducteur,level1)+offset(log(exposition)),
+ data=sinistres,family=poisson)
> summary(regglmc1)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                         -1.6036     0.1741  -9.212  < 2e-16 ***
cut(ageconducteur, level1)(20,25]   -0.4200     0.1948  -2.157   0.0310 *
cut(ageconducteur, level1)(25,30]   -0.9378     0.1903  -4.927 8.33e-07 ***
cut(ageconducteur, level1)(30,35]   -1.0030     0.1869  -5.367 8.02e-08 ***
cut(ageconducteur, level1)(35,40]   -1.0779     0.1866  -5.776 7.65e-09 ***
cut(ageconducteur, level1)(40,45]   -1.0264     0.1858  -5.526 3.28e-08 ***
cut(ageconducteur, level1)(45,50]   -0.9978     0.1856  -5.377 7.58e-08 ***
cut(ageconducteur, level1)(50,55]   -1.0137     0.1855  -5.464 4.65e-08 ***
cut(ageconducteur, level1)(55,60]   -1.2036     0.1939  -6.207 5.40e-10 ***
cut(ageconducteur, level1)(60,65]   -1.1411     0.2008  -5.684 1.31e-08 ***
cut(ageconducteur, level1)(65,70]   -1.2114     0.2085  -5.811 6.22e-09 ***
cut(ageconducteur, level1)(70,75]   -1.3285     0.2210  -6.012 1.83e-09 ***
cut(ageconducteur, level1)(75,80]   -0.9814     0.2271  -4.321 1.55e-05 ***
cut(ageconducteur, level1)(80,85]   -1.4782     0.3371  -4.385 1.16e-05 ***
cut(ageconducteur, level1)(85,90]   -1.2120     0.5294  -2.289   0.0221 *
cut(ageconducteur, level1)(90,95]   -0.9728     1.0150  -0.958   0.3379
cut(ageconducteur, level1)(95,100] -11.4694   144.2817  -0.079   0.9366
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> yp=predict(regglmc1,newdata=data.frame(ageconducteur=a,exposition=1),
+ type="response",se.fit=TRUE)
> yp0=yp$fit > yp1=yp$fit+2*yp$se.fit > yp2=yp$fit-2*yp$se.fit > plot(a,yp0,ylim=c(.03,.12),type="s") > abline(v=40,col="grey") > lines(a,yp1,lty=2,type="s") > lines(a,yp2,lty=2,type="s") Here we obtain the following predictions, and for our 40 year old driver, the frequency is now 6.84%. Frequency = 0.0684573 confidence interval 0.07766717 0.05924742 But our classes were defined arbitrarily here. Perhaps should we consider other classes, to see if the prediction is sensitive to the cutting values, > level2=level1-2 > regglmc2=glm(nbre~cut(ageconducteur,level2)+offset(log(exposition)), + data=sinistres,family=poisson) which yields the following values for our 40 year old driver, Frequency = 0.07050614 confidence interval 0.07980422 0.06120807 So here, we did not remove the discontinuity problem. An idea here can be to consider moving regions: if the goal is to predict the frequency for a 40 year old driver, perhaps the class should be (somehow) centered around 40. And center the interval around 35 for drivers aged 35. Etc. • Moving average Thus, it is natural to consider some local regressions, where only drivers aged almost 40 should be considered. This almost concept is related to the bandwidth. For instance, drivers between 35 and 45 can be considered as being almost40. In practice we can either consider a subset function, or we can use weights in the regressions > value=40 > h=5 > sinistres$omega=(abs(sinistres$ageconducteur-value)<=h)*1 > regglmomega=glm(nbre~ageconducteur+offset(log(exposition)), + data=sinistres,family=poisson,weights=omega) To see what’s going on, let us consider an animated plot, where the age of interest is changing, Here, for our 40 year old drive, we get Frequency = 0.06913391 confidence interval 0.07535564 0.06291218 We do obtain a curve that can be interpreted as a local regression. But here, we do not take into account that 35 is not as close to 40 as 39 could be. An here, 34 is assumed to be very far away from 40. Clearly, we could improve that technique: kernel functions can considered, i.e. the closer to 40, the larger the weight. > value=40 > h=5 > sinistres$omega=dnorm(abs(sinistres$ageconducteur-value)/h) > regglmomega=glm(nbre~ageconducteur+offset(log(exposition)), + data=sinistres,family=poisson,weights=omega) which can be plotted below Here, our prediction for our 40 year old drive is Frequency = 0.07040464 confidence interval 0.07981521 0.06099408 This is the idea of kernel regression techniques. But as explained in the slides, other non parametric techniques can be considered, like spline functions. • Smoothing with splines In R, it is simple to use spline function (somehow much more simple than kernel smoothers) > library(splines) > regglmbs=glm(nbre~bs(ageconducteur)+offset(log(exposition)), + data=sinistres,family=poisson) The prediction for our 40 year old driver is now Frequency = 0.06928169 confidence interval 0.07397124 0.06459215 Note that this techniques is related to another class of models, the so-called Generalized Additive Models, i.e. GAMs. > library(mgcv) > reggam=gam(nbre~s(ageconducteur)+offset(log(exposition)), + data=sinistres,family=poisson) The prediction is extremely close to the one we obtained above (the main differences being observed for very old drivers) Frequency = 0.06912683 confidence interval 0.07501663 0.06323702 • Comparison of the different models Somehow, one way or another, all those models are valid. So perhaps we should compare them, On the graph above, we can visualize the upper and the lower bound of the prediction, for the 9 models. The horizontal line is the predicted value without taking into account heterogeneity. It is possible to consider relative values, with respect to this value, # Sondages et prévisions de séries temporelles Ce soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant sondage=read.table("http://freakonometrics.blog.free.fr/ public/data/sondage1ertour2012.csv", header=TRUE,sep=";",dec=",") sondage[36,1]="11/03/12" sondeur=unique(sondage$SONDEUR)
sondage$date=as.Date(as.character(sondage$DATE),
"%d/%m/%y")

(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
couleur=c(CL[1],CL[6],"orange","grey")
datefin=as.Date("2012/04/12")
k=1
plot(sondage$date,sondage[,k+2],col=couleur[k], pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40), xlim=c(min(sondage$date),datefin))
rl=lowess(sondage$date,sondage[,k+2]) lines(rl,lwd=2,col=couleur[k]) for(k in 2:4){ points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
rl=lowess(sondage$date,sondage[,k+2]) lines(rl,lwd=2,col=couleur[k]) } Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste rl=lowess(sondage$date,sondage[,k+2])

Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser),

D=data.frame(date=min(sondage$date)+0:148) library(splines) k=1 plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin)) rs=lm(sondage[,k+2]~bs(date),data=sondage) prl=predict(rs,newdata=D) prlse=sqrt(prl/100*(1-prl/100)/1000)*100 polygon(c(D$date,rev(D$date)),c(prl+2*prlse, rev(prl-2*prlse)),col=CL[3],border=NA) lines(D$date,prl,lwd=2,col=CL[1])
points(sondage$date,sondage[,k+2], col=couleur[k],pch=19,cex=.7,) for(k in 2:4){ points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
}

avec, là encore, la commande suivante pour lisser

rs=lm(sondage[,k+2]~bs(date),data=sondage)

Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets),

On a ici la courbe suivante

avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

k=2
plot(sondage$date,sondage[,k+2], col=couleur[k],pch=19,cex=.7, xlab="",ylab="",ylim=c(0,40),xlim= c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[4],border=NA)
lines(D$date,prl,lwd=2,col=CL[6]) points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7,)
for(k in c(1,3:4)){
points(sondage$date,sondage[,k+2],col= couleur[k],pch=19,cex=.7) } Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête, Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“. Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait… Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via D=data.frame(date=min(sondage$date)+0:148)
prl=predict(rs,newdata=D)

On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections)

On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

library(mnormt)
simulation=function(S){
proba=c(S,100-sum(S))/100
variance= -proba%*%t(proba)
diag(variance)=proba*(1-proba)
variance=variance/1000
return(rmnorm(1,proba[1:4],variance[1:4,1:4]))}
simulsondages=function(M){
Z=rep(NA,ncol(M))
for(i in 1:nrow(M)){
Z=rbind(Z,simulation(M[i,]))
}
return(Z[-1,])}
prediction4=matrix(NA,nrow(D),4)
for(k in 1:4){
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prediction4[,k]=predict(rs,newdata=D)
}

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

set.seed(1)
S4=simulsondages(prediction4)
S100=100*S4
k=1
plot(D$date,S100[,k],col=couleur[k], pch=19,cex=.7, xlab="",ylab="",ylim=c(0,40),xlim= c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D), lwd=2,col=couleur[k]) for(k in 2:4){ points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D), lwd=2,col=couleur[k]) } (les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé), library(forecast) k=1 X=S100[,k] ETS=ets(X) F=forecast(ETS,h=60) fdate=max(D$date)+1:60
k=1
plot(D$date,S100[,k],col=couleur[k], pch=19,cex=.7, xlab="",ylab="",ylim=c(0,40),xlim= c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D), lwd=2,col=couleur[k]) for(k in 2:4){ points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k]) } polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[3],border=NA) polygon(c(fdate,rev(fdate)), c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))), col=CL[2],border=NA) lines(fdate,as.numeric(F$mean),lwd=2,col=CL[1])

Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes

ETS=ets(X)
F=forecast(ETS,h=60)

Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

k=2
X=S100[,k]
ETS = ets(X)
F=forecast(ETS,h=60)
k=1
plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7, xlab="",ylab="",ylim=c(0,40),xlim= c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D), lwd=2,col=couleur[k]) for(k in 2:4){ points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D), lwd=2,col=couleur[k]) } polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[4],border=NA) polygon(c(fdate,rev(fdate)), c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))), col=CL[5],border=NA) lines(fdate,as.numeric(F$mean),lwd=2,col=CL[6])

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

set.seed(1)
for(sim in 1:20){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
lines(fdate,as.numeric(F1$mean),col=CL[1]) lines(fdate,as.numeric(F2$mean),col=CL[6])}

Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour

set.seed(1)
VICTOIRE=rep(NA,1000)
for(sim in 1:1000){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
VICTOIRE[sim]=(F1$mean[30]>F2$mean[30])}
mean(VICTOIRE)

Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes !