Big Data/Analytics Zone is brought to you in partnership with:

Arthur Charpentier, ENSAE, PhD in Mathematics (KU Leuven), Fellow of the French Institute of Actuaries, professor at UQàM in Actuarial Science. Former professor-assistant at ENSAE Paritech, associate professor at Ecole Polytechnique and professor assistant in economics at Université de Rennes 1. Arthur is a DZone MVB and is not an employee of DZone and has posted 160 posts at DZone. You can read more from them at their website. View Full User Profile

# Exposure with Binomial Responses

03.04.2013
| 1365 views |

Last week, we’ve seen how to take into account the exposure to compute nonparametric estimators of several quantities (empirical means, and empirical variances) incorporating exposure. Let us see what can be done if we want to model a binomial response. The model here is the following: ,

• the number of claims $N_i$ on the period $[0,1]$
• the number of claims $Y_i$ on $[0,E_i]$)

that can be visualize below

Consider the case where the variable of interest is not the number of claims, but simply the indicator of the occurrence of a claim. Then we wish to model the event $\{N=0\}$ versus $\{N>0\}$, interpreted as non-occurrence and occurrence. Given the fact that we can only observe $\{Y=0\}$ versus $\{Y>0\}$. Having an inclusion is not enough to derive a model. Actually, with a Poisson process model, we can get easily that

$\mathbb{P}(Y=0) = \mathbb{P}(N=0)^E$

With words, it means that the probability of not having a claim in the first six months of the year is the square root of not have a claim over a year. Which makes sense. Assume that the probability of not having a claim can be explained by some covariates, denoted $\boldsymbol{X}$, through some link function (using the GLM terminology),

$\mathbb{P}(N=0|\boldsymbol{X})=h(\boldsymbol{X}'\boldymbol{\beta})$

Now, since we do observe $Y$ – and not $N$ – we have
$\mathbb{P}(Y=0|\boldsymbol{X},E)=h(\boldsymbol{X}'\boldymbol{\beta})^E$

The dataset we will use is always the same

> sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
> sinistres=sinistre[sinistre$garantie=="1RC",] > sinistres=sinistres[sinistres$cout>0,]
> T=table(sinistres$nocontrat) > T1=as.numeric(names(T)) > T2=as.numeric(T) > nombre1 = data.frame(nocontrat=T1,nbre=T2) > I = contrat$nocontrat%in%T1
> T1= contrat$nocontrat[I==FALSE] > nombre2 = data.frame(nocontrat=T1,nbre=0) > nombre=rbind(nombre1,nombre2) > sinistres = merge(contrat,nombre) > sinistres$nonsin = (sinistres$nbre==0) The first model we can consider is based on the standard logistic approach, i.e. $\mathbb{P}(Y=0|\boldsymbol{X},E)=\left(\frac{\exp(\boldsymbol{X}'\boldymbol{\beta})}{1+\exp(\boldsymbol{X}'\boldymbol{\beta})}\right)^E$ That’s nice, but difficult to handle with standard functions. Nevertheless, it is always possible to compute numerically the maximum likelihood estimator of $\boldymbol{\beta}$ given $(Y_i,\boldsymbol{X}_i,E_i)$. > Y=sinistres$nonsin
> X=cbind(1,sinistres$ageconducteur) > E=sinistres$exposition
> logL = function(beta){
+ 	pi=(exp(X%*%beta)/(1+exp(X%*%beta)))^E
+ 	-sum(log(dbinom(Y,size=1,prob=pi)))
+ }
> optim(fn=logL,par=c(-0.0001,-.001),
+ method="BFGS")
$par [1] 2.14420560 0.01040707$value
[1] 7604.073
$counts function gradient 42 10$convergence
[1] 0
$message NULL > parametres=optim(fn=logL,par=c(-0.0001,-.001), + method="BFGS")$par

Now, let us look at alternatives, based on standard regression models. For instance a binomial-log model. Because the exposure appears as a power of the annual probability, everything would be fine if $h$ was the exponential function (or $h^{-1}$ was the log link function), since

$\mathbb{P}(Y=0|\boldsymbol{X},E)=\exp(E+\boldsymbol{X}'\boldymbol{\beta})$
Now, if we try to code it, it starts quickly to be problematic,

> reg=glm(nonsin~ageconducteur+offset(exposition),
Error: no valid set of coefficients has been found: please supply starting values

I tried (almost) everything I could, but I could not get rid of that error message,

> startglm=c(0,-.001)
> names(startglm)=c("(Intercept)","ageconducteur")
> etaglm=rep(-.01,nrow(sinistresI))
> etaglm[sinistresI$nonsin==0]=-10 > muglm=exp(etaglm) > reg=glm(nonsin~ageconducteur+offset(exposition), + data=sinistresI,family=binomial(link="log"), + control = glm.control(epsilon=1e-5,trace=TRUE,maxit=50), + start=startglm, + etastart=etaglm,mustart=muglm) Deviance = NaN Iterations - 1 Error: no valid set of coefficients has been found: please supply starting values So I decided to give up. Almost. Actually, the problem comes from the fact that $\mathbb{P}(Y=0)$ is closed to 1. I guess everything would be nicer if we could work with probability close to 0. Which is possible, since $\mathbb{P}(Y>0)=1-\mathbb{P}(Y=0) = 1-[1-\mathbb{P}(N>0)]^E$ where $\mathbb{P}(N>0)$ is close to 0. So we can use Taylor’s expansion, $\mathbb{P}(Y>0)\sim1-1+E\cdot \mathbb{P}(N>0)]=E\cdot \mathbb{P}(N>0)]$ Here, the exposure does no longer appears as a power of the probability, but appears multiplicatively. Of course, there are higher order terms. But let us forget them (so far). If – one more time – we consider a log link function, then we can incorporate the exposure, or to be more specific, the logarithm of the exposure. > regopp=glm((1-nonsin)~ageconducteur+offset(log(exposition)), + data=sinistresI,family=binomial(link="log")) which now works perfectly. Now, to see a final model, perhaps we should get back to our Poisson regression model since we do have a model for the probability that $\mathbb{P}(Y=\cdot)$. > regpois=glm(nbre~ageconducteur+offset(log(exposition)), + data=sinistres,family=poisson(link="log")) We can now compare those three models. Perhaps, we should also include the prediction without any explanatory variable. For the second model (actually, it does run without any explanatory variable), we run > regreff=glm((1-nonsin)~1+offset(log(exposition)), + data=sinistres,family=binomial(link="log")) so that the prediction is here > exp(coefficients(regreff)) (Intercept) 0.06776376 This value is comparable with the logistic regression, > logL2 = function(beta){ + pi=(exp(beta)/(1+exp(beta)))^E + -sum(log(dbinom(Y,size=1,prob=pi)))} > param=optim(fn=logL2,par=.01,method="BFGS")$par
> 1-exp(param)/(1+exp(param))
[1] 0.06747777

But is quite different from the Poisson model,

> exp(coefficients(glm(nbre~1+offset(log(exposition)),
(Intercept)
0.07279295

Let us produce a graph, to compare those models,

> age=18:100
> yml1=exp(parametres[1]+parametres[2]*age)/(1+exp(parametres[1]+parametres[2]*age))
> plot(age,1-yml1,type="l",col="purple")
> yp=predict(regpois,newdata=data.frame(ageconducteur=age,
+ exposition=1),type="response")
> yp1=1-exp(-yp)
> ydl=predict(regopp,newdata=data.frame(ageconducteur=age,
+ exposition=1),type="response")
> plot(age,ydl,type="l",col="red")
> lines(age,yp1,type="l",col="blue")
> lines(age,1-yml1,type="l",col="purple")
> abline(h=exp(coefficients(regreff)),lty=2)

Observe here that the three models are quite different. Actually, with two models, it is possible to run more complex regression, e.g. with splines, to visualize the impact of the age on the probability of having – or not – a car accident. If we compare the Poisson regression (still in red) and the log-binomial model, with Taylor’s expansion, we get

The next step is to see how to incorporate the exposure in a tree. But that’s another story…

Published at DZone with permission of Arthur Charpentier, author and DZone MVB. (source)

(Note: Opinions expressed in this article and its replies are the opinions of their respective authors and not those of DZone, Inc.)

Tags:
"Starting from scratch" is seductive but disease ridden