Tag Archives: splines

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 !

Infidelity and econometrics

On http://www.bakadesuyo.com, there was recently an interesting discussion about infidelity, the key question being “at what ages are men and women most likely to have affairs?” The discussion is based on some graphs, e.g.

The source is a paper by Donald Cox. Based on a sample of 36 men and 22 women 3,432 respondent (NHSLS dataset) . And to be honest, I have been surprised by the shape of the curves. Especially for men… In order to compare, it is possible to use another dataset that can be found in R,

> library(Ecdat)
> data(Fair)
> tail(Fair)
sex age   ym child religious education occupation
596   male  47 15.0   yes         3        16          4
597   male  22  1.5   yes         1        12          2
598 female  32 10.0   yes         2        18          5
599   male  32 10.0   yes         2        17          6
600   male  22  7.0   yes         3        18          6
601 female  32 15.0   yes         3        14          1
rate nbaffairs
596    2         7
597    5         1
598    4         7
599    5         2
600    2         2
601    5         1

with 601 observations (from Fair (1977)). It is possible to run a Poisson regression to describe the number of affairs in the past year. E.g for men

> library(splines)
> regM=glm(nbaffairs~bs(age),family=poisson,
+ data=Fair[Fair$sex=="male",]) > a=seq(20,60) > N=predict(regM,newdata=data.frame(age=a),type="response") > plot(a,N,type="l",lwd=2,col="red") or for women, > regF=glm(nbaffairs~bs(age),family=poisson, + data=Fair[Fair$sex=="female",])
> N=predict(regF,newdata=data.frame(age=a),type="response")
> plot(a,N,type="l",lwd=2,col="red",lty=2)

On that (larger) dataset, we obtain curves that are more intuitive… But maybe the Poisson distribution is not an appropriate model. For instance, having no affairs do not mean that the person did not want to… So perhaps, a more interesting model would be a Poisson model with a zero-inflation, i.e. some people are honest and do not want to have affairs (and appear as 0), while some do want to have some affairs, and the number of affairs is Poisson distributed (and can take the value 0). If we focus on people wo do not want to have affairs, the model (and the prediction) is the following, where we plot the probability of not being interested in having an affair,

> library(pscl)
> regM0=zeroinfl(nbaffairs~bs(age)|bs(age),family=poisson,
> N0=predict(regM0,newdata=data.frame(age=a),type="zero")
> plot(a,N0,type="l",lwd=2,col="blue")

For those willing to have an affair, here is the parameter of the Poisson distribution of the number of affairs,

> Nc=predict(regM0,newdata=data.frame(age=a),type="count")
> plot(a,Nc,type="l",lwd=2,col="purple")

The same can be done for women, with the probability of no-willing to have an affair,

and to Poisson rate for women willing to have an affair,

If we focus on people willing to have an affair, the curves are the following,

i.e. men below 40 have more interested, but after 40, the probability drops, while women are still more and more likely to be willing to have an affair. On the other hand, young women having affairs might be less, but they usually have much more affairs than men…