# Exotic link functions for GLMs

In my previous post on GLMs, I discussed power link functions. But there are much more links that can be used : • The square root link (for the Poisson model)

Consider some random variable $Y$ with mean $\mu$ and variance $\sigma^2$. Using Taylor’s expansion,$$g(Y)\sim g(\mu)+(Y-\mu)g'(\mu)+\frac{1}{2}(Y-\mu)^2g''(\mu)$$we can write$$\mathbb{E}[g(Y)]\sim g(\mu)+\frac{\sigma^2}{2}g''(\mu)$$ $$\text{Var}[g(Y)]\sim [g'(\mu)]^2\sigma^2$$

Assume that $Y\sim\mathcal{P}(\lambda)$, a consider a square root transformation, $g(y)=\sqrt{y}$, then the second equality becomes $$\text{Var}[\sqrt{Y}]\sim \left[\frac{1}{2\sqrt{\mathbb{E}[Y]}}\right]^2\text{Var}[Y]=\frac{1}{4}$$

So, somehow, with a square-root transformation, we have variance stability, which might be interpreted as some homoscedasticity.

• The complementary log-log function for the Bernoulli model

Assume that the true variable of interest is a Poisson one, $N|\mathbf{X}=\mathbf{x}\sim\mathcal{P}(\lambda_{\mathbf{x}})$ where $\lambda_{\mathbf{x}}=\exp[\mathbf{x}^T\mathbf{\beta}]$Thus,$$\mathbb{P}[N=0|\mathbf{X}=\mathbf{x}]=\exp[-\lambda_{\mathbf{x}}]=\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]$$while$$\mathbb{P}[N>0|\mathbf{X}=\mathbf{x}]=1-\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]=H(\mathbf{x}^T\mathbf{\beta})$$where $H(s)=1-\exp[-\exp(s)]$. Let $Y=\mathbf{1}(N>0)$. The previous model seems like a Bernoulli regression with $H$ as link function,$$\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=H(\mathbf{x}^T\mathbf{\beta})$$

So, assume now that instead of observing $N$ we observe $Y=\boldsymbol{1}(N>0)$. In that case, running a Bernoulli regression with a complementary log-log link function would be the same (?) as running first a Poisson regression on the original data, and then use it on our binary variable, zero vs. non-zero. Let us generate some data, and see what’s going on. Let us compare $e^{\lambda_{\mathbf{x}}}$ and $p_{\mathbf{x}}$ obtained from a standard logistic regression

n=563 set.seed(1) base=data.frame(X1=rnorm(n),X2=rnorm(n)) lambda=base$X1+base$X2 base$Y=rpois(n,exp(lambda)) regPois = glm(Y~.,data=base,family=poisson(link="log")) lambda = predict(regPois,type="response") regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit")) prob = predict(regBinom, type="response") plot(prob,exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") What if $p_{\mathbf{x}}$ was obtained from a Bernoulli regression, with a cloglog link function ? regBinom = glm((Y&gt;0)~.,data=base,family=binomial(link="cloglog")) prob = predict(regBinom, type="response") plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") It looks like the fit is very good here ! Now, what if we have real data, like the dataset from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy (with 563 observations, and nine variables) base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE) str(base) x=base$SEX base$SEX="M" base$SEX[x=="0"]="F" x=base$CHILDREN base$CHILDREN="YES" base$CHILDREN[x==0]="NO" regPois = glm(Y~.,data=base,family=poisson(link="log")) lambda = predict(regPois,type="response") regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit")) prob = predict(regBinom, type="response") plot(prob,exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") In that case the two models are very different. But actually, so is the second one regBinom = glm((Y&gt;0)~.,data=base,family=binomial(link="cloglog")) prob = predict(regBinom, type="response") plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") How can we interpret that ? Could it be because the Poisson model is not good ? Actually, if we run a zero-inflated model here, library(pscl) regZIP = zeroinfl(Y ~ . | ., data = base) summary(regZIP) Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.002274 0.048413 -0.047 0.963 X1 1.019814 0.026186 38.945 &lt;2e-16 *** X2 1.004814 0.024172 41.570 &lt;2e-16 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -4.90190 2.07846 -2.358 0.0184 * X1 -2.00227 0.86897 -2.304 0.0212 * X2 -0.01545 0.96121 -0.016 0.9872 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Hence, we reject here the Poisson distribution assumption, because of the inflation of zeros… It looks like the cloglog link can be used to check if the Poisson distribution is a good model, or not… # Bailey (1963) and Poisson regression on two factors Consider the following dataset, from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy, with 563 observations, and nine variables : eight covariates, and the variable of interest, the number of extramarital affairs, over a year, base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE) str(base) 'data.frame': 563 obs. of 9 variables:$ SEX : int 1 0 0 1 1 0 0 1 0 1 ... $AGE : num 37 27 32 57 22 32 22 57 32 22 ...$ YEARMARRIAGE: num 10 4 15 15 0.75 1.5 0.75 15 15 1.5 ... $CHILDREN : int 0 0 1 1 0 0 0 1 1 0 ...$ RELIGIOUS : int 3 4 1 5 2 2 2 2 4 4 ... $EDUCATION : int 18 14 12 18 17 17 12 14 16 14 ...$ OCCUPATION : int 7 6 1 6 6 5 1 4 1 4 ... $SATISFACTION: int 4 4 4 5 3 5 3 4 2 5 ...$ Y : int 0 0 0 0 0 0 0 0 0 0 ...

Let us focus on two categorical covariates, related to the importance of religion, and the occupation

df=data.frame(y=base$Y, religion=as.factor(base$RELIGIOUS), occupation=as.factor(base$OCCUPATION), expo = 1) (E=xtabs(expo~religion+occupation,data=df)) occupation religion 1 2 3 4 5 6 7 1 4 1 8 4 16 9 0 2 23 3 11 17 56 36 6 3 29 1 10 12 39 25 2 4 38 7 12 21 59 44 2 5 13 1 3 10 19 19 3 (N=xtabs(y~religion+occupation,data=df)) occupation religion 1 2 3 4 5 6 7 1 4 1 13 3 13 7 0 2 1 1 13 10 25 43 10 3 15 0 12 11 34 35 1 4 24 1 3 15 11 9 10 5 6 0 0 6 11 7 0 The two tables above are the exposure (number of observations) and the number of extramarital affairs, here as contingency tables. Without any covariate, one can assume that $N\sim\mathcal{P}(\lambda\cdot E)$, where $\lambda$ would be sum(N)/sum(E)  0.6305506 The idea with the margin method is to assume that $N_{i,j}=E_{i,j}\cdot\lambda_{i,j}$ where $\lambda_{i,j}=A_i\cdot B_j$. Bailey (1963) added two series of constraints : per row, $$\sum_j N_{i,j}=\sum_j E_{i,j}\cdot A_i\cdot B_j$$ for any $i$ and similarly, for any $j$ $$\sum_i N_{i,j}=\sum_i E_{i,j}\cdot A_i\cdot B_j$$From the first series of constraints, write $$A_i=\frac{\sum_j N_{i,j}}{\sum_j E_{i,j}\cdot B_j}$$ and use the second series to write $$B_j=\frac{\sum_i N_{i,j}}{\sum_i E_{i,j}\cdot A_i}$$Because we need $A_i$‘s to compute $B_j$‘s, and conversely, it is natural to consider some iterative procedure to solve it. Observe that we do not have unicity… Consider here some starting values for $A_i$‘s and $B_j$‘s A=rep(1,length(levels(df$religion))) B=rep(1,length(levels(df$occupation)))*sum(N)/sum(E) A  1 1 1 1 1 B  0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 The predicted number of extramarital affairs would be $\hat N_{i,j}=E_{i,j}\cdot\hat A_i\cdot \hat B_j$ E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.5222025 0.6305506 5.0444050 2.5222025 10.0888099 5.6749556 0.0000000 2 14.5026643 1.8916519 6.9360568 10.7193606 35.3108348 22.6998224 3.7833037 3 18.2859680 0.6305506 6.3055062 7.5666075 24.5914742 15.7637655 1.2611012 4 23.9609236 4.4138544 7.5666075 13.2415631 37.2024867 27.7442274 1.2611012 5 8.1971581 0.6305506 1.8916519 6.3055062 11.9804618 11.9804618 1.8916519 sum(B*E[1,])  26.48313 sum(B*E[2,])  95.84369 apply(t(B*t(E)),1,sum) 1 2 3 4 5 26.48313 95.84369 74.40497 115.39076 42.87744 sum(A*E[,1])  107 sum(A*E[,2])  13 apply(A*E,2,sum) 1 2 3 4 5 6 7 107 13 44 64 189 133 13 From expressions above, observe that one can very easily write expressions of $A_i$‘s and $B_j$‘s as functions of $B_j$‘s and $A_i$‘s respectively A=apply(N,1,sum)/apply(t(B*t(E)),1,sum) B=apply(N,2,sum)/apply(A*E,2,sum) Let it iterate one thousand times for(i in 1:1000){ A=apply(N,1,sum)/apply(t(B*t(E)),1,sum) B=apply(N,2,sum)/apply(A*E,2,sum) } We obtain here A 1 2 3 4 5 1.5404346 1.0447195 1.4825650 0.6553159 0.6634763 B 1 2 3 4 5 6 7 0.4685515 0.2629769 0.8454435 0.7245310 0.4889697 0.7770553 1.6753750 E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586111 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450811 0.3898804 12.5342484 12.8899708 28.2722423 28.8008726 4.9677044 4 11.6678702 1.2063307 6.6483904 9.9707299 18.9053460 22.4055332 2.1957997 5 4.0413463 0.1744790 1.6827951 4.8070914 6.1639760 9.7955975 3.3347148 That is our prediction, per category, of the number of affairs. Observe that here, sums per row are equal to observed numbers, apply(N,1,sum) 1 2 3 4 5 41 103 108 73 30 apply(E * A%*%t(B),1,sum) 1 2 3 4 5 41 103 108 73 30 as well as sums per colums apply(N,2,sum) 1 2 3 4 5 6 7 50 3 41 45 94 101 21 apply(E * A%*%t(B),2,sum) 1 2 3 4 5 6 7 50 3 41 45 94 101 21 Now, why should I mention that here, in the section on the Poisson regression in our course ? Because actually, this is exactly what we get if we run a Poisson regression on those two covariates reg=glm(y~religion+occupation,data=df,family=poisson) summary(reg) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.32604 0.21325 -1.529 0.126285 religion2 -0.38832 0.18791 -2.066 0.038783 * religion3 -0.03829 0.18585 -0.206 0.836771 religion4 -0.85470 0.19757 -4.326 1.52e-05 *** religion5 -0.84233 0.24416 -3.450 0.000561 *** occupation2 -0.57758 0.59549 -0.970 0.332083 occupation3 0.59022 0.21349 2.765 0.005699 ** occupation4 0.43588 0.20603 2.116 0.034381 * occupation5 0.04265 0.17590 0.242 0.808399 occupation6 0.50587 0.17360 2.914 0.003569 ** occupation7 1.27415 0.26298 4.845 1.27e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 First of all, observe that the total sum of predictions equals the total sum of observations yp = predict(reg,type="response") sum(yp)  355 sum(df$y)  355

But actually, the predicted number of affairs, for our 35 classes, is exactly what we got using Bailey’s technique

xtabs(yp~df$religion+df$occupation) df$occupation df$religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586112 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450813 0.3898804 12.5342484 12.8899708 28.2722424 28.8008726 4.9677044 4 11.6678703 1.2063307 6.6483904 9.9707300 18.9053460 22.4055332 2.1957997 5 4.0413464 0.1744790 1.6827951 4.8070914 6.1639761 9.7955975 3.3347148 E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586111 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450811 0.3898804 12.5342484 12.8899708 28.2722423 28.8008726 4.9677044 4 11.6678702 1.2063307 6.6483904 9.9707299 18.9053460 22.4055332 2.1957997 5 4.0413463 0.1744790 1.6827951 4.8070914 6.1639760 9.7955975 3.3347148

To be more specific, up to a multiplicate constant, series of coefficients are equal here, e.g. for $A_i$‘s

a=exp(coefficients(reg)+c(0,coefficients(reg)[2:5])) a/a religion2 religion3 religion4 religion5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 A/A 1 2 3 4 5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072

but also for $B_j$‘s

b=exp(coefficients(reg)+c(0,coefficients(reg)[6:11])) b/b occupation2 occupation3 occupation4 occupation5 occupation6 occupation7 1.0000000 0.5612551 1.8043769 1.5463210 1.0435773 1.6584203 3.5756477 B/B 1 2 3 4 5 6 7 1.0000000 0.5612551 1.8043770 1.5463210 1.0435773 1.6584203 3.5756478

This will have major implications in non-life insurance models (for claims reserving).

# Actuariat de l’Assurance Non-Vie #5

Petit retour sur les modèles de comptage, avec la digression sur les GLM, avec les modèles permettant d’avoir de la surdispersion (théoriquement impossible dans un modèle de Poisson, mais plus réaliste compte tenu des particularités des données d’assurance). Les slides sont en ligne.

# Actuariat de l’Assurance Non-Vie #3

Mardi, seconde partie du cours d’actuariat, avec les modèles de classification le matin, mais l’après midi, on devrait commencer les modèles de comptage. Les slides sont en ligne.

Le matin, je dois intervenir 15 minutes à l’IHP vers 9 heures, dans un colloque sur Artificial Intelligence for Fintech and Insurtech: si le RER traîne un peu trop, entre Luxembourg et Lozère, j’aurais peut-être 5 minutes de retard….

# Actuariat de l’Assurance non-Vie #3

Après plusieurs jours d’interruption, le cours d’actuariat non-vie reprendra ce lundi. Les slides sont en ligne. Il s’agit d’une introduction aux modèles de comptage.

# Econometrics vs. Machine Learning with Temporal Patterns

A few months ago, I did publish a (long) post entitled ‘some thoughts on economics, mathematics, econometrics, machine learning, etc‘. In that post, I was discussing possible differences between foundations of econometrics, and machine learning. I wanted to get back today on an important point, related to training/sampling datasets, when we have temporal data.

I was discussing this morning, with a student of the Data Science for Actuaries program, an interesting point related to claim frequency models, for insurance ratemaking. Since the goal is to predict claims frequency (to assess the level of the insurance premium), he suggested to use old data to train the model, and more recent one to test it. The problem is that the model did not incorporate any temporal pattern, and we got surprising results.

Consider here a simple dataset,

> set.seed(1)
> n=50000
> X1=runif(n)
> T=sample(2000:2015,size=n,replace=TRUE)
> L=exp(-3+X1-(T-2000)/20)
> E=rbeta(n,5,1)
> Y=rpois(n,L*E)
> B=data.frame(Y,X1,L,T,E)

Claims frequency is driven by a Poisson process, with one covariate, X1, and we assume that the intensity decreases (with an exponential rate). Consider here a standard linear regression, without any time effect

> reg=glm(Y~X1+offset(log(E)),data=B,
+ family=poisson)

We can also compute the empirical annualized claims frequency

> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B[abs(B$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp=Vectorize(p)(seq(.05,.95,by=.1)) and plot the two curves on the same graph, > plot(seq(.05,.95,by=.1),vp,type="b") > lines(u,exp(v),lty=2,col="red") This is what we usually do in econometrics. In machine learning, and more specifically to assess the quality of the model, and for model selection, it is common to split the dataset in two parts. A training sample, and a validation sample. Consider some randomized training/validation samples, then fit a model on the training sample, and finally use it to get a prediction, > idx=sample(1:nrow(B),size=nrow(B)*7/8) > B_a=B[idx,] > B_t=B[-idx,] > reg=glm(Y~X1+offset(log(E)),data=B_a, + family=poisson) > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B_a[abs(B_a$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_t=Vectorize(p)(seq(.05,.95,by=.1)) > lines(seq(.05,.95,by=.1),vp_t,col="red") The blue curve is the prediction on the training sample (as we usually do in econometrics), but then the red curve is the prediction on the testing sample. Here, volatility probably comes from the small size of the testing sample (1 observation out of 8). Now, what if we use the year as a splitting criteria : we fit a model on old years to fit a model, and we test it on recent years, > B_a=subset(B,T<2014) > B_t=subset(B,T>=2014) > reg=glm(Y~X1+offset(log(E)),data=B_a,family=poisson) > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B_a[abs(B_a$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_t=Vectorize(p)(seq(.05,.95,by=.1)) > lines(seq(.05,.95,by=.1),vp_t,col="red") Clearly, we miss something here… We were looking at such a graph this morning, and it took me some time to understand how training and validation samples were designed, and that there was a possible temporal effect (actually, this morning, it was based on a 3 year training sample, and a 1 year validation sample). Since there is a temporal pattern, let us capture it. As an econometrician, let me use a regression model > reg=glm(Y~X1+T+offset(log(E)),data=B, + family=poisson) > C=coefficients(reg) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C+C*(2000:2015))) > lines(u,v,lty=2,col="red") (I focus only on the evolution of the temporal variate on that graph). Here, we use a linear model, but there are usually no reason to assume linearity. So we might consider splines > library(splines) > reg=glm(Y~X1+bs(T)+offset(log(E)), + data=B,family=poisson) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > v2=predict(reg,newdata=data.frame(X1=0, + T=2000:2015,E=1)) > plot(2000:2015,exp(v2),type="b") > lines(u,v,lty=2,col="red") But here again, why should we assume that there is an underlying smooth function? There might be some ruptures… So let us consider a regression on factors > reg=glm(Y~0+X1+as.factor(T)+offset(log(E)), + data=B,family=poisson) > C=coefficients(reg) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[2:17]),type="b") > lines(u,v,lty=2,col="red") An alternative might be to consider some more general model, like a regression tree > library(rpart) > reg=rpart(Y~X1+T+offset(log(E)),data=B, + method="poisson",cp=1e-4) > p=function(t){ + B=B[B$T==t,]
+   B$E=1 + mean(predict(reg,newdata=B)) + } > y_m=Vectorize(function(t) p(t))(2000:2015) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3+.5) > plot(2000:2015,y_m,ylim=c(.02,.085),type="b") > lines(u,v,lty=2,col="red") Here, it seems that something went wrong. I guess it’s coming from the exposure. So consider a simplier model, on the annualized frequency, and with weights that are related to the exposure > reg=rpart(Y/E~X1+T,data=B,weights=B$E,cp=1e-4)
> p=function(t){
+   B=B[B$T==t,] + B$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red") That was for the econometrician perspective. With a machine learning perspective, consider a training sample (here based on old data) and a validation sample (based on more recent ones)

> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)

If we consider a model, it is easy to get a prediction on recent years, even if the model was designed to model older ones,

> reg_a=glm(Y~X1+T+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C+C*c(2000:2013,
+ NA,NA)),type="b")
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(C+C*2014:2015),
+ pch=19,col="blue") But if we use years as factors, things are more complicated.

> reg_a=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> RMSE=function(A){
+   L=exp(C*B_t$X1+ A*(B_t$T==2014) + A*(B_t$T==2015)) + Y_t=L*B_t$E
+   sum( (Y_t - B_t$Y )^2)} > i=optim(c(.4,.4),RMSE)$par
> plot(2000:2015,c(exp(C[2:15]),NA,NA),)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(i),pch=19,col="blue") becase we need to get a prediction on levels that were not in our training sample. Here, we minimize the RMSE to quantify factor levels for recent years. And the output is not that bad.

So yes, it is possible to get a training dataset on older data, and test it on recent years. But one should be careful, and take into account, properly, temporal patterns.

# Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that  has a similar distribution… so there is no reason to test for a Gamma model for  before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if  has a joint Gaussien distribution, then both marginals are Gaussian, but also . So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for . The econometric interpretation is that with a standard Gaussian linear model, if is normally distributed, not only the conditional distribution  is Gaussian but also the non-conditional distribution of .

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20)) Indeed, here the distribution of  is also Gaussian

> library(nortest)

Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually) Here  is continuous. What if we consider a finite mixture here, i.e. takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for . But in practice, would we really reject the Gaussian assumption, for ? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
+ ad.test(Y)$p.value} > plot(2:100,Vectorize(pv)(2:100),type="l") > abline(h=.05,col="red") So here, it could be possible to have also a Gaussian distribution, for . As least to accept that assumption, statistically. In the context of a Poisson regression, it is well know that it’s not possible to have at the same time that is Poisson distributed (that’s a Poisson regression) and also that is Poisson distributed. That simply comes from the fact that while and because of the conditional Poisson distribution, then Thus, So cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for . More generally, it is very difficult to have a distribution family for that is also the distribution of the non-conditional variable . In the context of a finite mixture ( takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference). Hence, as a keep saying, before running a regression model on with some given family, it is never a good idea to check if the non-conditional distribution has the same distribution. Because there is no reason, usually, to remain in the same family. # De la fréquence des ouragans aux Etats-Unis Pour la fin du cours de statistique, on joue enfin avec des données. Après les taux de réussite au brevet des collèges, on a travaillé un peu sur la survenance des ouragans aux Etats-Unis. La base utilisée est la suivante > base=read.table( "http://freakonometrics.free.fr/ouragan.csv", sep=";",header=TRUE,dec=",") > base=base[1:207,1:7] Le nombre annuel d’ouragan est la série suivante > TB <- table(base$Year)
>  years <- as.numeric(names(TB))
>  counts <- as.numeric(TB)
>  years0=(1900:2005)[which(!(1900:2005)%in%years)]
>  db <- data.frame(years=c(years,years0),
counts=c(counts,rep(0,length(years0))))
>  X=db$counts > plot(db,type="h") On peut comparer les fréquences observées avec ce que donnerait une loi de Poisson (paramètre obtenu par maximum de vraisemblance) > lambda=mean(X) > n=length(X) > cbind(X,dpois(0:6,lambda)*n) [,1] [,2] 0 16 15.038430 1 34 29.367500 2 27 28.674870 3 13 18.665717 4 5 9.112744 5 6 3.559128 6 5 1.158396 On peut alors faire un test d’ajustement, > library(vcd) > fit=goodfit(X,"poisson") > plot(fit) > summary(fit) Goodness-of-fit test for poisson distribution X^2 df P(> X^2) Likelihood Ratio 14.17686 5 0.01452415 On voit que le test d’ajustement d’une loi Poisson n’est pas très concluant ici. Notons que ni la loi binomiale > fit=goodfit(X,"binomial") > plot(fit) > summary(fit) Goodness-of-fit test for binomial distribution X^2 df P(> X^2) Likelihood Ratio 46.21196 5 8.222959e-09 ni la loi binomiale négative ne semble marcher. > fit=goodfit(X,"nbinomial") > plot(fit) > summary(fit) Goodness-of-fit test for nbinomial distribution X^2 df P(> X^2) Likelihood Ratio 10.91507 4 0.02753523 On peut faire un test de stabilité temporelle du paramètre d’une loi de Poisson, à l’aide du test de Przyborowski & Wilenski (1940). Un ancien package permettait de mettre en place ce test. Mais on peut heureusement utiliser une autre fonction, avec une version différente du test (on pourra aussi consulter Krishnamoorthy & Thomson (2004) pour aller plus loin). L’idée du premier test est assez simple: notons que si sur la première période (de durée ), on a un processus de Poisson de paramètre et sur la seconde période (de durée ), alors et Un test de s’écrit Notons que la loi conditionnelle de sachant suit une loi binomiale, donc le paramètre de probabilité est fonction de ce ratio. Plus précisément, conditionnellement à, suit un loi binomiale avec A partir de là, il est facile de mettre en place une stratégie de test. Une autre stratégie est d’utiliser > seuil=1960 > X1=db[db$years<seuil,"counts"]
>  X2=db[db$years>=seuil,"counts"] > N1=length(db$years<seuil)
>  N2=length(db$years>=seuil) > poisson.test(c(sum(X1),sum(X2)),c(N1,N2)) Comparison of Poisson rates data: c(sum(X1), sum(X2)) time base: c(N1, N2) count1 = 85, expected count1 = 103.5, p-value = 0.01216 alternative hypothesis: true rate ratio is not equal to 1 95 percent confidence interval: 0.5218662 0.9265666 sample estimates: rate ratio 0.6967213 Bref, on peut accepter que le nombre annuel d’ouragan n’est pas un processus de Poisson homogène dans le temps… # Autres Données pour le cours de Statistiques Encore quelques données pour faire un TP de statistique. Ce sont des données de coûts d’ouragans, aux Etats-Unis, > library(gdata) > db=read.xls( + "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls", + sheet=1) > stupidcomma = function(x){ + x=as.character(x) + for(i in 1:10){x=sub(",","",as.character(x))} + return(as.numeric(x))} > base=db[,1:4] > base$Base.Economic.Damage=
Vectorize(stupidcomma)(db$Base.Economic.Damage) > base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05) > base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05) Lire la base est un peu pénible, à cause de virgules qui traînent comme séparateur de milliers. Mais peu importe, le petit code permet de nettoyer la base. • Fréquence d’ouragans La première étape pourra être de travailler sur la fréquence. On pourra ajuster des lois de Poisson, ou binomiales, ou binomiales négatives, et faire des tests de stabilité du paramètre, avec le temps (on pourra couper en 4 ou 5 blocs). > TB <- table(base$Year)
> years <- as.numeric(names(TB))
> counts <- as.numeric(TB)
> years0=(1900:2005)[which(!(1900:2005)%in%years)]
> db <- data.frame(years=c(years,years0),
+       counts=c(counts,rep(0,length(years0))))
> db[88:93,]
years counts
88  2003      3
89  2004      6
90  2005      6
91  1902      0
92  1905      0
93  1907      0
> plot(years,counts,type='h') • Coût des ouragans

Dans un second temps, on pourra travailler sur le coût des ouragans

> plot(base$Normalized.PL05/1e9,type="h") On pourra tenter d’ajuster des lois, log normale, gamma, Pareto, etc. Faire des tests d’ajustement. Et tester la stabilité des paramètres, dans le temps (en coupant en 2 blocs, par exemple, ou davantage). # Actuariat de l’Assurance Non-Vie #5 Pour le quatrième cours d’actuariat de l’assurance non-vie à l’ENSAE, la semaine prochaine, on abordera la notion de sur-dispersion dans les modèles de comptage. Les slides sont en ligne (la version pdf téléchargeable est plus complète que celle sur slideshare) # Actuariat de l’Assurance Non-Vie #3 Pour le troisème cours d’actuariat de l’assurance non-vie à l’ENSAE (oui, j’ai pris un peu d’avance), on abordera la modélisation de la fréquence de sinistres. Les slides sont en ligne, # 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, family=poisson(link="log")))  438.6314 > AIC(glm(dist ~ speed, data=train_cars, family=poisson(link="identity")))  436.3997 > AIC(glm(dist ~ bs(speed), data=train_cars, family=poisson(link="log")))  425.6434 > AIC(glm(dist ~ bs(speed), data=train_cars, family=poisson(link="identity")))  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. # 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… # Reinterpreting Lee-Carter Mortality Model Last week, while I was giving my crash course on R for insurance, we’ve been discussing possible extensions of Lee & Carter (1992) model. If we look at the seminal paper, the model is defined as follows # The odds of a cluster of airplane accidents Recently, there have been a lot of airplane accidents. • July, 17th 2014, Hrabove, Ukraine, Malaysia Airlines, Boeing 777, fatalities 298 (/298) • July, 23rd 2014, Magong, Taiwan, TransAsia Airways, ATR 72-500, fatalities 47 (/58) • July, 24th 2014, Aguelhok, Mali, Air Algerie, Mc Donnell Douglas MD-83, fatalities 116 (/116) It is simple to find a lot of datasets about airplane crashes. For instance on http://ntsb.gov/aviationquery. The dataset is nice, with a lot of information, > planes=read.table( + "cbad3ca6-6b8f-4c98-9ee0-601faAviationData.txt", + sep="|",header=TRUE) for instance the exact location of the crashes, > library(maps) > map("world", interior = FALSE) > points(planes$Longitude,planes$Latitude, + pch=19,cex=planes$Total.Fatal.Injuries/50,
+ col="red") 