# Classification from scratch, bagging and forests 10/8

Tenth post of our series on classification from scratch. Today, we’ll see the heuristics of the algorithm inside bagging techniques.

Often, bagging is associated with trees, to generate forests. But actually, it is possible using bagging for any kind of model. Recall that bagging means “boostrap aggregation”. So, consider a model $m:\mathcal{X}\rightarrow \mathcal{Y}$. Let $\widehat{m}_{S}$ denote the estimator of $m$ obtained from sample $S=\{y_i,\mathbf{x}_i\}$ with $i=\{1,\cdots,n\}$.

Consider now some boostrap sample, $S_b=\{y_i,\mathbf{x}_i\}$ with $i$ is randomly drawn from $\{1,\cdots,n\}$ (with replacement). Based on that sample, estimate $\widehat{m}_{S_b}$. Then draw many samples, and consider the agregation of the estimators obtained, using either a majority rule, or using the average of probabilities (if a probabilist model was considered). Hence$$\widehat{m}^{bag}(\mathbf{x})=\frac{1}{B}\sum_{b=1}^B \widehat{m}_{S_b}(\mathbf{x})$$

## Bagging logistic regression #1

Consider the case of the logistic regression. To generate a bootstrap sample, it is natural to use the technique describe above. I.e. draw pairs $(y_i,\mathbf{x}_i)$ randomly, uniformly (with probability $1/n$) with replacement. Consider here the small dataset, just to visualize. For the b part of bagging, use the following code

L_logit = list() n = nrow(df) for(s in 1:1000){ df_s = df[sample(1:n,size=n,replace=TRUE),] L_logit[[s]] = glm(y~., df_s, family=binomial)}

Then we should aggregate over the 1000 models, to get the agg part of bagging,

p = function(x){ nd=data.frame(x1=x[1], x2=x[2]) unlist(lapply(1:1000,function(z) predict(L_logit[[z]],newdata=nd,type="response")))}

We now have a prediction for any new observation

vu = seq(0,1,length=101) vv = outer(vu,vu,Vectorize(function(x,y) mean(p(c(x,y))))) image(vu,vu,vv,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(vu,vu,vv,levels = .5,add=TRUE) ## Bagging logistic regression #2 Another technique that can be used to generate a bootstrap sample is to keep all $\mathbf{x}_i$‘s, but for each of them, to draw (randomly) a value for $y$, with$$Y_{i,b}\sim\mathcal{B}(\widehat{m}_{S}(\mathbf{x}_i))$$since$$\widehat{m}(\mathbf{x})=\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}].$$Thus, the code for the b part of bagging algorithm is now L_logit = list() n = nrow(df) reg = glm(y~x1+x2, df, family=binomial) for(s in 1:100){ df_s = df df_s$y = factor(rbinom(n,size=1,prob=predict(reg,type="response")),labels=0:1) L_logit[[s]] = glm(y~., df_s, family=binomial) }

The agg part of bagging algorithm remains unchanged. Here we obtain

vu = seq(0,1,length=101) vv = outer(vu,vu,Vectorize(function(x,y) mean(p(c(x,y))))) image(vu,vu,vv,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(vu,vu,vv,levels = .5,add=TRUE) Of course, we can use that code we check the prediction obtain on the observations we have in our sample. Just to change, consider here the myocarde data. The entiere code is here L_logit = list() reg = glm(as.factor(PRONO)~., myocarde, family=binomial) for(s in 1:1000){ myocarde_s = myocarde myocarde_s$PRONO = 1*rbinom(n,size=1,prob=predict(reg,type="response")) L_logit[[s]] = glm(as.factor(PRONO)~., myocarde_s, family=binomial) } p = function(x){ nd=data.frame(FRCAR=x[1], INCAR=x[2], INSYS=x[3], PRDIA=x[4], PAPUL=x[4], PVENT=x[5], REPUL=x[6]) unlist(lapply(1:1000,function(z) predict(L_logit[[z]],newdata=nd,type="response")))}

For the first observation, with our 1000 simulated datasets, and our 1000 models, we obtained the following estimation for the probability to die.

## Fronm bags to forest

Here, we grew a lot of trees, but it is not stricto sensus a random forest algorithm, as introduced in 1995, in Random decision forests. Actually, the difference is in the creation of decision trees. To understand what happens, get back to the previous post on classification trees. As we’ve seen, when we have a node, we look at possible splits : we consider all possible variable, and all possible threshold. The startegy here will be to draw randomly $k$ variables out of $p$ (with of course $k, for instance $k=\sqrt{p}$). That's interesting in high dimension, because at each split, we should look for all variables, all cutoffs, and that can take quite some time (especially with the bootstrap procedure, where the goal will be to grow 1000 trees).

To be continued…

# Classification from scratch, penalized Lasso logistic 5/8

Fifth post of our series on classification from scratch, following the previous post on penalization using the $\ell_2$ norm (so-called Ridge regression), this time, we will discuss penalization based on the $\ell_1$ norm (the so-called Lasso regression).

First of all, one should admit that if the name stands for least absolute shrinkage and selection operator, that’s actually a very cool name… Funny story, a few years before, Leo Breiman introduce a concept of garrote technique… “The garrote eliminates some variables, shrinks others, and is relatively stable”.

I guess that somehow, the lasso is the extension of the garotte technique

## Normalization of the covariates

As previously, the first step will be to consider linear transformations of all covariates $x_j$ to get centered and scaled variables (with unit variance)

y = myocarde$PRONO X = myocarde[,1:7] for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j]) X = as.matrix(X) ## Ridge Regression (from scratch) The heuristics about Lasso regression is the following graph. In the background, we can visualize the (two-dimensional) log-likelihood of the logistic regression, and the blue square is the constraint we have, if we rewite the optimization problem as a contrained optimization problem, LogLik = function(bbeta){ b0=bbeta[1] beta=bbeta[-1] sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*log(1 + exp(b0+X%*%beta)))} u = seq(-4,4,length=251) v = outer(u,u,function(x,y) LogLik(c(1,x,y))) image(u,u,v,col=rev(heat.colors(25))) contour(u,u,v,add=TRUE) polygon(c(-1,0,1,0),c(0,1,0,-1),border="blue") The nice thing here is that is works as a variable selection tool, since some components can be null here. That’s the idea behind the following (popular) graph (with lasso on the left, and ridge on the right). Heuristically, the maths explanation is the following. Consider a simple regression $y_i=x_i\beta+\varepsilon$, with $\ell_1$-penality and a $\ell_2$-loss fuction. The optimization problem becomes$$\min\big\{\mathbf{y}^T\mathbf{y}-2\mathbf{y}^T\mathbf{x}\beta+\beta\mathbf{x}^T\mathbf{x}\beta+2\lambda{\color{red}{|}}\beta{\color{red}{|}}\big\}$$The first order condition can be written$$-2\mathbf{y}^T\mathbf{x}+2\mathbf{x}^T\mathbf{x}\widehat{\beta}{\color{red}{\pm} }2\lambda=0$$(the sign in ${\color{red}{\pm}}$ being the sign of $\widehat{\beta}$). Assume that $\mathbf{y}^T\mathbf{x}>0$, then solution is $$\widehat{\beta}_{\lambda}^{lasso}=\max\left\lbrace\frac{\mathbf{y}^T\mathbf{x}-\lambda}{\mathbf{x}^T\mathbf{x}},0\right\rbrace$$(we get a corner solution when $\lambda$ is large). ## Optimization routine As in our previous post, let us start with standard (R) optimization routines, such as BFGS PennegLogLik = function(bbeta,lambda=0){ b0=bbeta[1] beta=bbeta[-1] -sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*log(1 + exp(b0+X%*%beta)))+lambda*sum(abs(beta)) } opt_lasso = function(lambda){ beta_init = lm(PRONO~.,data=myocarde)$coefficients logistic_opt = optim(par = beta_init*0, function(x) PennegLogLik(x,lambda), hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) logistic_opt$par[-1] } v_lambda=c(exp(seq(-4,2,length=61))) est_lasso=Vectorize(opt_lasso)(v_lambda) library("RColorBrewer") colrs=brewer.pal(7,"Set1") plot(v_lambda,est_lasso[1,],col=colrs[1],type="l") for(i in 2:7) lines(v_lambda,est_lasso[i,],col=colrs[i],lwd=2) But it is very heratic… or non stable. ## Using glmnet Just to compare, with R routines dedicated to lasso, we get the following library(glmnet) glm_lasso = glmnet(X, y, alpha=1) plot(glm_lasso,xvar="lambda",col=colrs,lwd=2) plot(glm_lasso,col=colrs,lwd=2) If we look carefully what’s in the ouput, we can see that there is variable selection, in the sense that some $\widehat{\beta}_{j,\lambda}=0$, in the sense “really null” glmnet(X, y, alpha=1,lambda=exp(-4))$beta 7x1 sparse Matrix of class "dgCMatrix" s0 FRCAR . INCAR 0.11005070 INSYS 0.03231929 PRDIA . PAPUL . PVENT -0.03138089 REPUL -0.20962611

Of course, with out optimization routine, we cannot expect to have null values

opt_lasso(.2) FRCAR INCAR INSYS PRDIA 0.4810999782 0.0002813658 1.9117847987 -0.3873926427 PAPUL PVENT REPUL -0.0863050787 -0.4144139379 -1.3849264055

So clearly, it will be necessary to spend more time today, to understand how it works…

## Orthogonal covariates

Before getting into the maths, observe that when covariates are orthogonal, there is some very clear “variable” selection process,

library(factoextra) pca = princomp(X) pca_X = get_pca_ind(pca)$coord glm_lasso = glmnet(pca_X, y, alpha=1) plot(glm_lasso,xvar="lambda",col=colrs) plot(glm_lasso,col=colrs) ## Interior Point approach The penalty is now expressed using the $\ell_1$ so intuitively, it should be possible to consider algorithms related to linear programming. That was actually suggested in Koh, Kim & Boyd (2007), with some implementation in matlab, see http://web.stanford.edu/~boyd/l1_logreg/. If I can find some time, later one, maybe I will try to recode it. But actually, it is not the technique used in most R functions. Now, o be honest, we face a double challenge today: the first one is to understand how lasso works for the “standard” (least square) problem, the second one is to see how to adapt it to the logistic case. ## Standard lasso (with weights) If we get back to the original Lasso approach, the goal was to solve$$\min\left\lbrace\frac{1}{2n}\sum_{i=1}^n [y_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2+\lambda \sum_j |\beta_j|\right\rbrace$$(with standard notions, as in wikipedia or Jocelyn Chi’s post – most of the code in this section is inspired by Jocelyn’s great post). Observe that the intercept is not subject to the penalty. The first order condition is then$$\frac{\partial}{\partial\beta_0}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}-\beta_0\mathbf{1}\|^2=(\mathbf{X}\mathbf{\beta}-\mathbf{y})^T\mathbf{1}+\beta_0\|\mathbf{1}\|^2=0$$i.e.$$\beta_0=\frac{1}{n^2}(\mathbf{X}\mathbf{\beta}-\mathbf{y})^T\mathbf{1}$$Assume now that KKT conditions are satisfied, since we cannot differentiate (to find points where the gradient is $\mathbf{0}$), we can check if $\mathbf{0}$ contains the subdifferential at the minimum. Namely$$\mathbf{0}\in\partial \left(\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2+\lambda\|\mathbf{\beta}\|_{\ell_1}\right)=\frac{1}{2}\nabla\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2+\partial(\lambda\|\mathbf{\beta}\|_{\ell_1})$$ For the term on the left, we recognize $$\frac{1}{2}\nabla\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2=-\mathbf{X}^T(\mathbf{y}-\mathbf{X}\mathbf{\beta})=-\mathbf{g}$$so that the previous equation can be writen$$g_k\in\partial(\lambda|\beta_k|)=\begin{cases}\{+\lambda\}\text{ if }\beta_k>0 \\ \{-\lambda\}\text{ if }\beta_k<0 \\ (-\lambda,+\lambda)\text{ if }\beta_k=0\end{cases}$$i.e. if $\beta_k\neq 0$, then $g_k = \text{sign}(\beta_k)\cdot\lambda$. Then we write the KKT conditions for this formulation and simplify them to produce a set of rules for checking our solution We can split $\beta_j$ into a sum of its positive and negative parts by replacing $\beta_j$ with $\beta_j^+-\beta_j^-$ where $\beta_j^+,\beta_j^-\geq0$. Then the Lasso problem becomes$$-\log\mathcal{L}(\mathbf{\beta})+\lambda\sum_j(\beta_j^+-\beta_j^-)$$with constraints $\beta_j^+-\beta_j^-$. Let $\alpha_j^+,\alpha_j^-$ denote the Lagrange multipliers for $\beta_j^+,\beta_j^-$, respectively. $$L({\mathbf{\beta}}) + \lambda \sum_{j} (\beta_{j}^{+} - \beta_{j}^{-}) - \sum_{j}\alpha_{j}^{+}\beta_{j}^{+} - \sum_{j} \alpha_{j}^{-}\beta_{j}^{-}.$$To satisfy the stationarity condition, we take the gradient of the Lagrangian with respect to $\beta_{j}^{+}$ and set it to zero to obtain$$\nabla L({\mathbf{\beta}})_{j} + \lambda - \alpha_{j}^{+} = 0$$We do the same with respect to $\beta_{j}^{-}$ to obtain$$-\nabla L({\mathbf{\beta}})_{j}+\lambda-\alpha_{j}^{-} = 0$$ As discussed in Jocelyn Chi’s post, primal feasibility requires that the primal constraints be satisfied so this gives us $\beta_{j}^{+} \ge 0$ and $\beta_{j}^{-} \ge 0$. Then dual feasibility requires non-negativity of the Lagrange multipliers so we get $\alpha_{j}^{+} \ge 0$ and $\alpha_{j}^{-} \ge 0$. And finally, complementary slackness requires that $\alpha_{j}^{+}\beta_{j}^{+} = 0$ and $\alpha_{j}^{-}\beta_{j}^{-} = 0$. We can simplify these conditions to obtain a simple set of rules for checking whether or not our solution is a minimum. The following is inspired by Jocelyn Chi’s post. From $\nabla L(\beta)_{j} + \lambda - \alpha_{j}^{+} = 0$, we have $\nabla L(\beta)_{j} + \lambda= \alpha_{j}^{+} \ge 0$. This gives us $\nabla L(\beta)_{j} \ge -\lambda$. From $-\nabla L(\beta)_{j} + \lambda - \alpha_{j}^{-} = 0$, we have $-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} \ge 0$. This gives us $-\nabla L(\beta)_{j} \ge -\lambda$, which gives us $\nabla L(\beta)_{j} \le \lambda$. Hence, $\lvert \nabla L(\beta)_{j} \rvert \le \lambda \; \forall j$ When $\beta_{j}^{+} > 0, \lambda > 0$, complementary slackness requires $\alpha_{j}^{+} = 0$. So $\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{+} = 0$. Hence, $\nabla L(\beta)_{j} = -\lambda < 0$ since $\lambda > 0$. At the same time, $-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} \ge 0$ so $2 \lambda = \alpha_{j}^{-} > 0$ since $\lambda > 0$. Then complementary slackness requires $\beta_{j}^{-} = 0$. Hence, when $\beta_{j}^{+} > 0$, we have $\beta_{j}^{-}=0$ and $\nabla L(\beta)_{j} = -\lambda$ Similarly, when $\beta_{j}^{-} > 0, \lambda > 0$, complementary slackness requires $\alpha_{j}^{-}=0$. So $-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} = 0$ and $\nabla L(\beta)_{j}=\lambda>0$ since $\lambda > 0$. Then from $\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{+} \ge 0$ and the above, we get $2 \lambda = \alpha_{j}^{+} > 0$. Then complementary slackness requires $\beta_{j}^{+} = 0$. Hence, when $\beta_{j}^{-} > 0$, we have $\beta_{j}^{+}=0$ and $\nabla L(\beta)_{j} = \lambda$. Since $\beta_{j} = \beta_{j}^{+} - \beta_{j}^{-}$, this means that when $\beta_{j} > 0$, $\nabla L(\beta)_{j} = -\lambda$. And when $\beta_{j} <0$, $\nabla L(\beta)_{j} = \lambda$. Combining this with $\lvert \nabla L(\beta)_{j} \rvert \le \lambda \; \forall j$, we arrive at the same convergence requirements that we obtained before using subdifferential calculus. For conveniency, introduce the soft-thresholding function$$S(z,\gamma)=\text{sign}(z)\cdot(|z|-\gamma)_+=\begin{cases}z-\gamma&\text{ if }\gamma>|z|\text{ and }z<0\\z+\gamma&\text{ if }\gamma<|z|\text{ and }z<0 \\0&\text{ if }\gamma\geq|z|\end{cases}$$ Noticing that the optimization problem $$\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|_{\ell_2}^2+\lambda\|\mathbf{\beta}\|_{\ell_1}$$can also be written $$\min\left\lbrace\sum_{j=1}^p -\widehat{\beta}_j^{ols}\cdot\beta_j+\frac{1}{2}\beta_j^2+\lambda|\beta_j|\right\rbrace$$observe that$$\widehat{\beta}_{j,\lambda}=S(\widehat{\beta}_j^{ols},\lambda)$$which is a coordinate-wise update. Now, if we consider a (slightly) more general problem, with weights in the first part$$\min\left\lbrace\frac{1}{2n}\sum_{i=1}^n{\color{red}{\omega_i}} [y_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2+\lambda \sum_j |\beta_j|\right\rbrace$$the coordinate-wise update becomes $$\widehat{\beta}_{j,\lambda,{\color{red}{\omega}}}=S(\widehat{\beta}_j^{{\color{red}{\omega-}}ols},\lambda)$$ An alternative is to set$\mathbf{r}_j=\mathbf{y} - \left(\beta_0\mathbf{1}+\sum_{k\neq j}\beta_k\mathbf{x}_k\right)=\mathbf{y}-\widehat{\mathbf{y}}^{(j)}$ so that the optimization problem can be written, equivalently $$\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p [\mathbf{r}_j-\beta_j\mathbf{x}_j]^2+\lambda |\beta_j|\right\rbrace$$ hence$$\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p \beta_j^2\|\mathbf{x}_j\|-2\beta_j\mathbf{r}_j^T\mathbf{x}_j+\lambda |\beta_j|\right\rbrace$$ and one gets $$\beta_{j,\lambda} = \frac{1}{\|\mathbf{x}_j\|^2}S(\mathbf{r}_j^T\mathbf{x}_j,n\lambda)$$ or, if we develop $$\beta_{j,\lambda} = \frac{1}{\sum_i x_{ij}^2}S\left(\sum_ix_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)$$ Again, if there are weights $\mathbf{\omega}=(\omega_i)$, the coordinate-wise update becomes $$\beta_{j,\lambda,{\color{red}{\omega}}} = \frac{1}{\sum_i {\color{red}{\omega_i}}x_{ij}^2}S\left(\sum_i{\color{red}{\omega_i}}x_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)$$ The code to compute this componentwise descent is soft_thresholding = function(x,a){ result = numeric(length(x)) result[which(x &gt; a)] a)] - a result[which(x &lt; -a)] &lt;- x[which(x &lt; -a)] + a return(result) } and the code lasso_coord_desc = function(X,y,beta,lambda,tol=1e-6,maxiter=1000){ beta = as.matrix(beta) X = as.matrix(X) omega = rep(1/length(y),length(y)) obj = numeric(length=(maxiter+1)) betalist = list(length(maxiter+1)) betalist[[1]] = beta beta0list = numeric(length(maxiter+1)) beta0 = sum(y-X%*%beta)/(length(y)) beta0list[1] = beta0 for (j in 1:maxiter){ for (k in 1:length(beta)){ r = y - X[,-k]%*%beta[-k] - beta0*rep(1,length(y)) beta[k] = (1/sum(omega*X[,k]^2))*soft_thresholding(t(omega*r)%*%X[,k],length(y)*lambda) } beta0 = sum(y-X%*%beta)/(length(y)) beta0list[j+1] = beta0 betalist[[j+1]] = beta obj[j] = (1/2)*(1/length(y))*norm(omega*(y - X%*%beta - beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta)) if (norm(rbind(beta0list[j],betalist[[j]]) - rbind(beta0,beta),'F') &lt; tol) { break } } return(list(obj=obj[1:j],beta=beta,intercept=beta0)) } Let’s keep that one warm, and let’s get back to our initial problem. ## The lasso logistic regression The trick here is that the logistic problem can be formulated as a quadratic programming problem. Recall that the log-likelihood is here $$\log\mathcal{L}=\frac{1}{n}\sum_{i=1}^n y_i\cdot(\beta_0+\mathbf{x}_i^T\mathbf{\beta})-\log[1+\exp(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]$$ which is a concave function of the parameters. Hence, one can use a quadratic approximation of the log-likelihood – using Taylor expansion,$$\log\mathcal{L}\approx\log\mathcal{L}'=\frac{1}{n}\sum_{i=1}^n \omega_i\cdot[z_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2$$ where $z_i$ is the working response $$z_i=(\beta_0+\mathbf{x}_i^T\mathbf{\beta})+\frac{y_i-p_i}{p_i[1-p_i]}$$ $p_i$ is the prediction$$p_i = \frac{\exp[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}{1+\exp[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}$$and $\omega_i$ are weights $\omega_i = p_i[1-p_i]$. Thus, we obtain a penalized least-square problem. And we can use what was done previously lasso_coord_desc = function(X,y,beta,lambda,tol=1e-6,maxiter=1000){ beta = as.matrix(beta) X = as.matrix(X) obj = numeric(length=(maxiter+1)) betalist = list(length(maxiter+1)) betalist[[1]] = beta beta0 = sum(y-X%*%beta)/(length(y)) p = exp(beta0*rep(1,length(y)) + X%*%beta)/(1+exp(beta0*rep(1,length(y)) + X%*%beta)) z = beta0*rep(1,length(y)) + X%*%beta + (y-p)/(p*(1-p)) omega = p*(1-p)/(sum((p*(1-p)))) beta0list = numeric(length(maxiter+1)) beta0 = sum(y-X%*%beta)/(length(y)) beta0list[1] = beta0 for (j in 1:maxiter){ for (k in 1:length(beta)){ r = z - X[,-k]%*%beta[-k] - beta0*rep(1,length(y)) beta[k] = (1/sum(omega*X[,k]^2))*soft_thresholding(t(omega*r)%*%X[,k],length(y)*lambda) } beta0 = sum(y-X%*%beta)/(length(y)) beta0list[j+1] = beta0 betalist[[j+1]] = beta obj[j] = (1/2)*(1/length(y))*norm(omega*(z - X%*%beta - beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta)) p = exp(beta0*rep(1,length(y)) + X%*%beta)/(1+exp(beta0*rep(1,length(y)) + X%*%beta)) z = beta0*rep(1,length(y)) + X%*%beta + (y-p)/(p*(1-p)) omega = p*(1-p)/(sum((p*(1-p)))) if (norm(rbind(beta0list[j],betalist[[j]]) - rbind(beta0,beta),'F') &lt; tol) { break } } return(list(obj=obj[1:j],beta=beta,intercept=beta0)) } It looks like what can get when calling glmnet… and here, we do have null components for some $\lambda$ large enough ! Really null… and that’s cool actually. ## Application on our second dataset Consider now the second dataset, with two covariates. The code to get lasso estimates is df0 = df df0$y = as.numeric(df$y)-1 plot_lambda = function(lambda){ m = apply(df0,2,mean) s = apply(df0,2,sd) for(j in 1:2) df0[,j] &lt;- (df0[,j]-m[j])/s[j] reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=1,lambda=lambda) u = seq(0,1,length=101) p = function(x,y){ xt = (x-m[1])/s[1] yt = (y-m[2])/s[2] predict(reg,newx=cbind(x1=xt,x2=yt),type="response")} v = outer(u,u,p) image(u,u,v,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+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)}

Consider some small values, for [\lambda], so that we only have some sort of shrinkage of parameters,

reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=1) par(mfrow=c(1,2)) plot(reg,xvar="lambda",col=c("blue","red"),lwd=2) abline(v=exp(-2.8)) plot_lambda(exp(-2.8)) But with a larger $\lambda$, there is variable selection: here $\widehat{\beta}_{1,\lambda}=0$ reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=1) par(mfrow=c(1,2)) plot(reg,xvar="lambda",col=c("blue","red"),lwd=2) abline(v=exp(-2.1)) plot_lambda(exp(-2.1))

# Classification from scratch, logistic with kernels 3/8

Third post of our series on classification from scratch, following the previous post introducing smoothing techniques, with (b)-splines. Consider here kernel based techniques. Note that here, we do not use the “logistic” model… it is purely non-parametric.

## kernel based estimated, from scratch

I like kernels because they are somehow very intuitive. With GLMs, the goal is to estimate $\hat{m}(\mathbf{x})=\mathbb{E}(Y|\mathbf{X}=\mathbf{x})$. Heuritically, we want to compute the (conditional) expected value on the neighborhood of $\mathbf{x}$. If we consider some spatial model, where $\mathbf{x}$ is the location, we want the expected value of some variable $Y$, “on the neighborhood” of $\mathbf{x}$. A natural approach is to use some administrative region (county, departement, region, etc). This means that we have a partition of $\mathcal{X}$ (the space with the variable(s) lies). This will yield the regressogram, introduced in Tukey (1961). For convenience, assume some interval / rectangle / box type of partition. In the univariate case, consider $$\hat{m}_{\mathbf{a}}(x)=\frac{\sum_{i=1}^n \mathbf{1}(x_i\in[a_j,a_{j+1}))y_i}{\sum_{i=1}^n \mathbf{1}(x_i\in[a_j,a_{j+1}))}$$or the moving regressogram $$\hat{m}(x)=\frac{\sum_{i=1}^n \mathbf{1}(x_i\in[x\pm h])y_i}{\sum_{i=1}^n \mathbf{1}(x_i\in[x\pm h])}$$In that case, the neighborhood is defined as the interval $(x\pm h)$. That’s nice, but clearly very simplistic. If $\mathbf{x}_i=\mathbf{x}$ and $\mathbf{x}_j=\mathbf{x}-h+\varepsilon$ (with $\varepsilon>0$), both observations are used to compute the conditional expected value. But if $\mathbf{x}_{j'}=\mathbf{x}-h-\varepsilon$, only $\mathbf{x}_i$ is considered. Even if the distance between $\mathbf{x}_{j}$ and $\mathbf{x}_{j'}$ is extremely extremely small. Thus, a natural idea is to use weights that are function of the distance between $\mathbf{x}_{i}$‘s and $\mathbf{x}$.Use$$\tilde{m}(x)=\frac{\sum_{i=1}^ny_i\cdot k_h\left({x-x_i}\right)}{\sum_{i=1}^nk_h\left({x-x_i}\right)}$$where (classically)$$k_h(x)=k\left(\frac{x}{h}\right)$$for some kernel $k$ (a non-negative function that integrates to one) and some bandwidth $h$. Usually, kernels are denoted with capital letter $K$, but I prefer to use $k$, because it can be interpreted as the density of some random noise we add to all observations (independently).

Actually, one can derive that estimate by using kernel-based estimators of densities. Recall that$$\tilde{f}(\mathbf{y})=\frac{1}{n|\mathbf{H}|^{1/2}}\sum_{i=1}^n k\left(\mathbf{H}^{-1/2}(\mathbf{y}-\mathbf{y}_i)\right)$$
Now, use the fact that the expected value can be defined as$$m(x)=\int yf(y|x)dy=\frac{\int y f(y,x)dy}{\int f(y,x)dy}$$Consider now a bivariate (product) kernel to estimate the joint density. The numerator is estimated by$$\frac{1}{nh}\sum_{i=1}^n\int y_i k\left(t,\frac{x-x_i}{h}\right)dt=\frac{1}{nh}\sum_{i=1}^ny_i \kappa\left(\frac{x-x_i}{h}\right)$$while the denominator is estimated by$$\frac{1}{nh^2}\sum_{i=1}^n \int k\left(\frac{y-y_i}{h},\frac{x-x_i}{h}\right)=\frac{1}{nh}\sum_{i=1}^n\kappa\left(\frac{x-x_i}{h}\right)$$In a general setting, we still use product kernels between $Y$ and $\mathbf{X}$ and write $$\widehat{m}_{\mathbf{H}}(\mathbf{x})=\displaystyle{\frac{\sum_{i=1}^ny_i\cdot k_{\mathbf{H}}(\mathbf{x}_i-\mathbf{x})}{\sum_{i=1}^n k_{\mathbf{H}}(\mathbf{x}_i-\mathbf{x})}}$$for some symmetric positive definite bandwidth matrix $\mathbf{H}$, and $$k_{\mathbf{H}}(\mathbf{x})=\det[\mathbf{H}]^{-1}k(\mathbf{H}^{-1}\mathbf{x})$$

Now that we know what kernel estimates are, let us use them. For instance, assume that $k$ is the density of the $\mathcal{N}(0,1)$ distribution. At point $x$, with a bandwidth $h$ we get the following code

mean_x = function(x,bw){ w = dnorm((myocarde$INSYS-x)/bw, mean=0,sd=1) weighted.mean(myocarde$PRONO,w)} u = seq(5,55,length=201) v = Vectorize(function(x) mean_x(x,3))(u) plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19)

and of course, we can change the bandwidth.

v = Vectorize(function(x) mean_x(x,2))(u) plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19)

We observe what we can read in any textbook : with a smaller bandwidth, we get more variance, less bias. “More variance” means here more variability (since the neighborhood is smaller, there are less points to compute the average, and the estimate is more volatile), and “less bias” in the sense that the expected value is supposed to be compute at point $x$, so the smaller the neighborhood, the better.

## Using ksmooth R function

Actually, there is a function in R to compute this kernel regression.

reg = ksmooth(myocarde$INSYS,myocarde$PRONO,"normal",bandwidth = 2*exp(1)) plot(reg$x,reg$y,ylim=0:1,type="l",col="red",lwd=2,xlab="INSYS",ylab="") points(myocarde$INSYS,myocarde$PRONO,pch=19)

We can replicate our previous estimate. Nevertheless, the output is not a function, but two series of vectors. That’s nice to get a graph, but that’s all we get. Furthermore, as we can see, the bandwidth is not exactly the same as the one we used before. I did not find any information online, so I tried to replicate the function we wrote before

g=function(bk=3){ reg = ksmooth(myocarde$INSYS,myocarde$PRONO,"normal",bandwidth = bk) f=function(bm){ v = Vectorize(function(x) mean_x(x,bm))(reg$x) z=reg$y-v sum((z[!is.na(z)])^2)} optim(bk,f)$par} x=seq(1,10,by=.1) y=Vectorize(g)(x) plot(x,y) abline(0,exp(-1),col="red") abline(0,.37,col="blue") There is a slope of $0.37$, which is actually $e^{-1}$. Coincidence ? I don’t know to be honest… ## Application in higher dimension Consider now our bivariate dataset, and consider some product of univariate (Gaussian) kernels u = seq(0,1,length=101) p = function(x,y){ bw1 = .2; bw2 = .2 w = dnorm((df$x1-x)/bw1, mean=0,sd=1)* dnorm((df$x2-y)/bw2, mean=0,sd=1) weighted.mean(df$y=="1",w) } v = outer(u,u,Vectorize(p)) image(u,u,v,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) We get the following prediction Here, the different colors are probabilities. ## k-nearest neighbors An alternative is to consider a neighborhood not defined using a distance to point $\mathbf{x}$ but the $k$-neighbors, with the $n$ observations we got.$$\tilde{m}_k(\mathbf{x})=\frac{1}{n}\sum_{i=1}^n\omega_{i,k}(\mathbf{x})y_i$$ where $\omega_{i,k}(\mathbf{x})=n/k$ if $i\in\mathcal{I}_{\mathbf{x}}^k$ with $$\mathcal{I}_{\mathbf{x}}^k=\{i:\mathbf{x}_i\text{ one of the }k\text{ nearest observations to }\mathbf{x}\}$$ The difficult part here is that we need a valid distance. If units are very different on each component, using the Euclidean distance will be meaningless. So, quite naturally, let us consider here the Mahalanobis distance Sigma = var(myocarde[,1:7]) Sigma_Inv = solve(Sigma) d2_mahalanobis = function(x,y,Sinv){as.numeric(x-y)%*%Sinv%*%t(x-y)} k_closest = function(i,k){ vect_dist = function(j) d2_mahalanobis(myocarde[i,1:7],myocarde[j,1:7],Sigma_Inv) vect = Vectorize(vect_dist)((1:nrow(myocarde))) which((rank(vect)))} Here we have a function to find the $k$ closest neighbor for some observation. Then two things can be done to get a prediction. The goal is to predict a class, so we can think of using a majority rule : the prediction for $y_i$ is the same as the one the majority of the neighbors. k_majority = function(k){ Y=rep(NA,nrow(myocarde)) for(i in 1:length(Y)) Y[i] = sort(myocarde$PRONO[k_closest(i,k)])[(k+1)/2] return(Y)}

But we can also compute the proportion of black points among the closest neighbors. It can actually be interpreted as the probability to be black (that’s actually what was said at the beginning of this post, with kernels),

k_mean = function(k){ Y=rep(NA,nrow(myocarde)) for(i in 1:length(Y)) Y[i] = mean(myocarde$PRONO[k_closest(i,k)]) return(Y)} We can see on our dataset the observation, the prediction based on the majority rule, and the proportion of dead individuals among the 7 closest neighbors cbind(OBSERVED=myocarde$PRONO, MAJORITY=k_majority(7),PROPORTION=k_mean(7)) OBSERVED MAJORITY PROPORTION [1,] 1 1 0.7142857 [2,] 0 1 0.5714286 [3,] 0 0 0.1428571 [4,] 1 1 0.5714286 [5,] 0 1 0.7142857 [6,] 0 0 0.2857143 [7,] 1 1 0.7142857 [8,] 1 0 0.4285714 [9,] 1 1 0.7142857 [10,] 1 1 0.8571429 [11,] 1 1 1.0000000 [12,] 1 1 1.0000000

Here, we got a prediction for an observed point, located at $\boldsymbol{x}_i$, but actually, it is possible to seek the $k$ closest neighbors of any point $\boldsymbol{x}$. Back on our univariate example (to get a graph), we have

mean_x = function(x,k=9){ w = rank(abs(myocarde$INSYS-x),ties.method ="random") mean(myocarde$PRONO[which(w&lt;=9)])} u=seq(5,55,length=201) v=Vectorize(function(x) mean_x(x,3))(u) plot(u,v,ylim=0:1,type="l",col="red",lwd=2,xlab="INSYS",ylab="") points(myocarde$INSYS,myocarde$PRONO,pch=19)

That’s not very smooth, but we do not have a lot of points either.

If we use that technique on our two-dimensional dataset, we obtain the following

Sigma_Inv = solve(var(df[,c("x1","x2")])) u = seq(0,1,length=51) p = function(x,y){ k = 6 vect_dist = function(j) d2_mahalanobis(c(x,y),df[j,c("x1","x2")],Sigma_Inv) vect = Vectorize(vect_dist)(1:nrow(df)) idx = which(rank(vect)&lt;=k) return(mean((df$y==1)[idx]))} v = outer(u,u,Vectorize(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+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) This is the idea of local inference, using either kernel on a neighborhood of $\mathbf{x}$ or simply using the $k$ nearest neighbors. Next time, we will investigate penalized logistic regressions, to be continued # Classification from scratch, logistic regression 1/8 Let us start today our series on classification from scratch The logistic regression is based on the assumption that given covariates $\mathbf{x}$, $Y$ has a Bernoulli distribution,$$Y|\mathbf{X}=\mathbf{x}\sim\mathcal{B}(p_{\mathbf{x}}),~~~~p_\mathbf{x}=\frac{\exp[\mathbf{x}^T\mathbf{\beta}]}{1+\exp[\mathbf{x}^T\mathbf{\beta}]}$$The goal is to estimate parameter $\mathbf{\beta}$. Recall that the heuristics for the use of that function for the probability is that$$\log[\text{odds}(Y=1)]=\log\frac{\mathbb{P}[Y=1]}{\mathbb{P}[Y=0]}=\mathbf{x}^T\mathbf{\beta}$$ ## Maximimum of the (log)-likelihood function The log-likelihood is here$$\log\mathcal{L} = \sum_{i=1}^n y_i\log p_i+(1-y_i)\log (1-p_i)$$ where $p_{i}=(1+\exp[-\mathbf{x}_i^T\mathbf{\beta}])^{-1}$. Numerical techniques are based on (numerical) gradient descent to compute the maximum of the likelihood function. The (negative) log-likelihood is the following function y = myocarde$PRONO X = cbind(1,as.matrix(myocarde[,1:7])) negLogLik = function(beta){ -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta))) }

We use the minus sign since standard optimization routines compute minima, not maxima. Now, to find the minimum of that function, we need a starting point to initiate the algorithm

beta_init = lm(PRONO~.,data=myocarde)$coefficients Why not start with the parameter of the OLS. Somehow, we might think that at least, sign should be ok for instance. Anyway, we need a starting point, and let us use that one. logistic_opt = optim(par = beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) Here, we obtain  logistic_opt$par (Intercept) FRCAR INCAR INSYS 1.656926397 0.045234029 -2.119441743 0.204023835 PRDIA PAPUL PVENT REPUL -0.102420095 0.165823647 -0.081047525 -0.005992238

Let us verify here that this output is valid. For instance, what if we change the value of the starting point (randomly)

simu = function(i){ logistic_opt_i = optim(par = rnorm(8,0,3)*beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) logistic_opt_i$par[2:3] } v_beta = t(Vectorize(simu)(1:1000)) plot(v_beta) par(mfrow=c(1,2)) hist(v_beta[,1],xlab=names(myocarde)[1]) hist(v_beta[,2],xlab=names(myocarde)[2]) Ooops. There is a problem here. Clearly, we cannot rely on numerical optimization here. We can think about using another optimization routine library(optimx) logit = function(mX, vBeta) { exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) } logLikelihoodLogitStable = function(vBeta, mX, vY) { -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta))) + (1-vY)*(-log(1 + exp(mX %*% vBeta)))) } likelihoodScore = function(vBeta, mX, vY) { return(t(mX) %*% (logit(mX, vBeta) - vY) ) } optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable, method = 'L-BFGS-B', gr = likelihoodScore, mX = X, vY = y, hessian=TRUE) The optimum is here attr(optimLogitLBFGS, "details")[[2]] [,1] 0.066680272 FRCAR 0.003080542 INCAR 0.079031364 INSYS -0.001586194 PRDIA 0.040500697 PAPUL -0.041870705 PVENT -0.014162756 REPUL 0.195632244 Let’s be honest here, I do not feel confortable with those techniques. So, what happened here ? Here, the technique we use is based on the following idea,$$\mathbf{\beta}_{new}=\mathbf{\beta}_{old} -\left(\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}\right)^{-1}\cdot \frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}$$The problem is that my computer does not know this first and second derivatives. So it will compute them using approximation techniques. Actually, it is possible to use functions dedicated to such computation library(numDeriv) library(MASS) logit = function(x){1/(1+exp(-x))} logLik = function(beta, X, y){ -sum(y*log(logit(X%*%beta)) + (1-y)*log(1-logit(X%*%beta))) } optim_second = function(beta, num_iter){ LL = vector() for(i in 1:num_iter){ grad = (t(X)%*%(logit(X%*%beta) - y)) H = hessian(logLik, beta, method = "complex", X = X, y = y) beta = beta - ginv(H)%*%grad LL[i] = logLik(beta, X, y) } result = list(beta, H) return(result) } With our OLS starting point, we obtain opt0 = optim_second(beta_init,500) opt0[[1]] [,1] [1,] 0.951074420 [2,] 0.018860280 [3,] 0.275428978 [4,] 0.144803636 [5,] -0.058535606 [6,] 0.001182178 [7,] -0.108651776 [8,] -0.002940315 But if we try with another starting point opt1 = optim_second(beta_init*runif(8),500) opt1[[1]] [,1] [1,] 0.052894794 [2,] 0.024718435 [3,] 0.167953661 [4,] 0.171662947 [5,] -0.057458066 [6,] -0.011361034 [7,] -0.107532114 [8,] -0.002679064 Clearly, some coefficients are rather close. But other aren’t. From my point of viezw, that is a major problem (keep in mind that we do not deal here with massive data ! There are only 7 explanatory variables, and only 71 observations). Why not try to be clever, and use the analytical values of those derivatives ? Even if some people claim the oppositive, sometimes, it can actually be usefull to do the maths, instead of considering only numerical values. ## Newton (or Fisher) Algorithm If you open any Econometrics textbooks (one can also try to derive it), you will get $$\frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}=\mathbf{X}^T(\mathbf{y}-\mathbf{p}_{old})$$ while$$\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=-\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}$$ Y=myocarde$PRONO X=cbind(1,as.matrix(myocarde[,1:7])) colnames(X)=c("Inter",names(myocarde[,1:7])) beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} Observe that here, I use only ten iterations of the algorithm !  beta[,8:10] [,1] [,2] [,3] XInter -10.187641685 -10.187641696 -10.187641696 XFRCAR 0.138178119 0.138178119 0.138178119 XINCAR -5.862429035 -5.862429037 -5.862429037 XINSYS 0.717084018 0.717084018 0.717084018 XPRDIA -0.073668171 -0.073668171 -0.073668171 XPAPUL 0.016756506 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187 The thing is that is seems to converge extremely fast. And it is rather robust ! Look at what we get if we change our starting point beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} beta[,8:10] [,1] [,2] [,3] XInter -10.187641586 -10.187641696 -10.187641696 XFRCAR 0.138178118 0.138178119 0.138178119 XINCAR -5.862429017 -5.862429037 -5.862429037 XINSYS 0.717084013 0.717084018 0.717084018 XPRDIA -0.073668172 -0.073668171 -0.073668171 XPAPUL 0.016756508 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187

Nice, isn’t it? Looks like we got our winner, don’t we? And one can use the inverse of the Hessian matrix to get standard deviations.

## Weighted Least-Squares

Let us go one step further. We’ve seen that we want to compute something like$$\mathbf{\beta}_{new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}$$(if we do substitute matrices in the analytical expressions) where $\mathbf{z}=\mathbf{X}\mathbf{\beta}_{old}+\mathbf{\Delta}_{old}^{-1}[\mathbf{y}-\mathbf{p}_{old}]$. But actually, that’s simply a standard least-square problem$$\mathbf{\beta}_{new} = \text{argmin}\left\lbrace(\mathbf{z}-\mathbf{X}\mathbf{\beta})^T\mathbf{\Delta}_{old}^{-1}(\mathbf{z}-\mathbf{X}\mathbf{\beta})\right\rbrace$$The only problem here is that weights $\mathbf{\Delta}_{old}$ are functions of unknown $\mathbf{\beta}_{old}$. But actually, if we keep iterating, we should be able to solve it : given the $\mathbf{\beta}$ we got the weights, and with the weights, we can use weighted OLS to get an updated $\mathbf{\beta}$. That’s the idea of iteratively reweighted least squares.

The algorithm will be

df = myocarde beta_init = lm(PRONO~.,data=df)$coefficients X = cbind(1,as.matrix(myocarde[,1:7])) beta = beta_init for(s in 1:1000){ p = exp(X %*% beta) / (1+exp(X %*% beta)) omega = diag(nrow(df)) diag(omega) = (p*(1-p)) df$Z = X %*% beta + solve(omega) %*% (df$PRONO - p) beta = lm(Z~.,data=df[,-8], weights=diag(omega))$coefficients }

and the output is here

 beta (Intercept) FRCAR INCAR INSYS PRDIA -10.187641696 0.138178119 -5.862429037 0.717084018 -0.073668171 PAPUL PVENT REPUL 0.016756506 -0.106776012 -0.003154187

which is almost what we’ve obtained before. Nice isn’t it ? Actually, here we also have standard deviations of estimators

summary( lm(Z~.,data=df[,-8], weights=diag(omega)))   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) -10.187642 10.668138 -0.955 0.343 FRCAR 0.138178 0.102340 1.350 0.182 INCAR -5.862429 6.052560 -0.969 0.336 INSYS 0.717084 0.503527 1.424 0.159 PRDIA -0.073668 0.261549 -0.282 0.779 PAPUL 0.016757 0.306666 0.055 0.957 PVENT -0.106776 0.099145 -1.077 0.286 REPUL -0.003154 0.004386 -0.719 0.475

## The standard glm function

Of course, it is possible to use an R built-in function to get our estimate

summary(glm(PRONO~.,data=myocarde,family=binomial(link = "logit")))   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -10.187642 11.895227 -0.856 0.392 FRCAR 0.138178 0.114112 1.211 0.226 INCAR -5.862429 6.748785 -0.869 0.385 INSYS 0.717084 0.561445 1.277 0.202 PRDIA -0.073668 0.291636 -0.253 0.801 PAPUL 0.016757 0.341942 0.049 0.961 PVENT -0.106776 0.110550 -0.966 0.334 REPUL -0.003154 0.004891 -0.645 0.519

## Application and visualisation

Let us visualize the prediction obtained from the logistic regression, on our second dataset

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x,x2=y,y=as.factor(z)) reg = glm(y~x1+x2,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(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Here level curves – or iso-probabilities – are linear, so the space is divided in two (0 and 1, survival and death, white and black) by a straight line (or an hyperplane in higher dimension). Furthermore, since we have a linear model, if we change the cutoff (the threshold used to create the two classes), we obtain another straight line (or hyperplane) parallel to the first one.

Next time, we will introduce splines to smooth those continuous covariates… to be continued.

# Regression on factors

Most of our intuitions about regression models come from the Gaussian standard linear model. One interesting feature is that, when we have a factor explanatory variable, the sum of predictions per class is the sum of observations of the endogeneous variable, per class. To be more specific, consider some factor variable $x_1\in\{0,1\}$, and a regression model

$y_i=\beta_0+\beta_1 \boldsymbol{1}(x_1=1)+\beta_2 x_2+\varepsilon_i$

Use ordinary least squares to fit that model

$\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 \boldsymbol{1}(x_1=1)+\widehat{\beta}_2 x_2$

Then for all $x\in\{0,1\}$

$\sum_{i:x_i=x} y_i = \sum_{i:x_i=x} \widehat{y}_i$

> n=200 > X1=rep(0:1,each=n/2) > set.seed(1) > X2=runif(2*n) > L=X1-X2 > B=data.frame(Y=rnorm(n,L),X1=as.factor(X1),X2=X2) > pd=aggregate(x=B$Y,by=list(B$X1),mean)$x > pd [1] -0.4881735 0.5341301 > fit=lm(Y~X1+X2,data=B) > B2=data.frame(x=B$X1,y=predict(fit)) > aggregate(x=B2$y,by=list(B2$x),mean)$x [1] -0.4881735 0.5341301 # Classification on the German Credit Database In our data science course, this morning, we’ve use random forrest to improve prediction on the German Credit Dataset. The dataset is > url="http://freakonometrics.free.fr/german_credit.csv" > credit=read.csv(url, header = TRUE, sep = ",") Almost all variables are treated a numeric, but actually, most of them are factors, > str(credit) 'data.frame': 1000 obs. of 21 variables:$ Creditability   : int  1 1 1 1 1 1 1 1 1 1 ...
$Account.Balance : int 1 1 2 1 1 1 1 1 4 2 ...$ Duration        : int  18 9 12 12 12 10 8  ...
$Purpose : int 2 0 9 0 0 0 0 0 3 3 ... (etc). Let us convert categorical variables as factors, > F=c(1,2,4,5,7,8,9,10,11,12,13,15,16,17,18,19,20) > for(i in F) credit[,i]=as.factor(credit[,i]) Let us now create our training/calibration and validation/testing datasets, with proportion 1/3-2/3 > i_test=sample(1:nrow(credit),size=333) > i_calibration=(1:nrow(credit))[-i_test] The first model we can fit is a logistic regression, on selected covariates > LogisticModel <- glm(Creditability ~ Account.Balance + Payment.Status.of.Previous.Credit + Purpose + Length.of.current.employment + Sex...Marital.Status, family=binomial, data = credit[i_calibration,]) Based on that model, it is possible to draw the ROC curve, and to compute the AUC (on ne validation dataset) > fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) > library(ROCR) > pred = prediction( fitLog, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCLog1=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCLog1,"\n")
AUC:  0.7340997

An alternative is to consider a logistic regression on all explanatory variables

> LogisticModel <- glm(Creditability ~ .,
+  family=binomial,
+  data = credit[i_calibration,])

We might overfit, here, and we should observe that on the ROC curve

> fitLog <- predict(LogisticModel,type="response",
+                   newdata=credit[i_test,])
> pred = prediction( fitLog, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCLog2=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCLog2,"\n") AUC: 0.7609792 There is a slight improvement here, compared with the previous model, where only five explanatory variables were considered. Consider now some regression tree (on all covariates) > library(rpart) > ArbreModel <- rpart(Creditability ~ ., + data = credit[i_calibration,]) We can visualize the tree using > library(rpart.plot) > prp(ArbreModel,type=2,extra=1) The ROC curve for that model is > fitArbre <- predict(ArbreModel, + newdata=credit[i_test,], + type="prob")[,2] > pred = prediction( fitArbre, credit$Creditability[i_test])
> perf <- performance(pred, "tpr", "fpr")
> plot(perf)
> AUCArbre=performance(pred, measure = "auc")@y.values[[1]]
> cat("AUC: ",AUCArbre,"\n")
AUC:  0.7100323

As expected, a single has a lower performance, compared with a logistic regression. And a natural idea is to grow several trees using some boostrap procedure, and then to agregate those predictions.

> library(randomForest)
> RF <- randomForest(Creditability ~ .,
+ data = credit[i_calibration,])
> fitForet <- predict(RF,
+                     newdata=credit[i_test,],
+                     type="prob")[,2]
> pred = prediction( fitForet, credit$Creditability[i_test]) > perf <- performance(pred, "tpr", "fpr") > plot(perf) > AUCRF=performance(pred, measure = "auc")@y.values[[1]] > cat("AUC: ",AUCRF,"\n") AUC: 0.7682367 Here this model is (slightly) better than the logistic regression. Actually, if we create many training/validation samples, and compare the AUC, we can observe that – on average – random forests perform better than logistic regressions, > AUC=function(i){ + set.seed(i) + i_test=sample(1:nrow(credit),size=333) + i_calibration=(1:nrow(credit))[-i_test] + LogisticModel <- glm(Creditability ~ ., + family=binomial, + data = credit[i_calibration,]) + summary(LogisticModel) + fitLog <- predict(LogisticModel,type="response", + newdata=credit[i_test,]) + library(ROCR) + pred = prediction( fitLog, credit$Creditability[i_test])
+   AUCLog2=performance(pred, measure = "auc")@y.values[[1]]
+   RF <- randomForest(Creditability ~ .,
+   data = credit[i_calibration,])
+   fitForet <- predict(RF,
+                       newdata=credit[i_test,],
+                       type="prob")[,2]

# Variable Selection using Cross-Validation (and Other Techniques)

A natural technique to select variables in the context of generalized linear models is to use a stepŵise procedure. It is natural, but contreversial, as discussed by Frank Harrell  in a great post, clearly worth reading. Frank mentioned about 10 points against a stepwise procedure.

• It yields R-squared values that are badly biased to be high.
• The F and chi-squared tests quoted next to each variable on the printout do not have the claimed distribution.
• The method yields confidence intervals for effects and predicted values that are falsely narrow (see Altman and Andersen (1989)).
• It yields p-values that do not have the proper meaning, and the proper correction for them is a difficult problem.
• It gives biased regression coefficients that need shrinkage (the coefficients for remaining variables are too large (see Tibshirani (1996)).
• It has severe problems in the presence of collinearity.
• It is based on methods (e.g., F tests for nested models) that were intended to be used to test prespecified hypotheses.
• Increasing the sample size does not help very much (see Derksen and Keselman (1992)).
• It allows us to not think about the problem.
• It uses a lot of paper.

# Visualising a Classification in High Dimension, part 2

A few weeks ago, I published a post on Visualising a Classification in High Dimension, based on the use of a principal component analysis, to get a projection on the first two components. Following that post, I was wondering what could be done in the context of a classification on categorical covariates. A natural idea would be to consider a correspondance analysis, and to run a similar code.

Consider here the dataset used in a recent post,

> source("http://freakonometrics.free.fr/import_data_credit.R")

If we consider a correspondance analysis, we get

> library(FactoMineR)
> acm=MCA(train.db,quali.sup =
+ which(names(train.db,)=="class"),ncp=10)

For the covariates (including also the variable we want to model, considered here as some supplementary variable), the visualisation – on the first two components – is

and for the individuals

# Visualising a Classification in High Dimension

So far, when discussing classification, we’ve been playing on my toy-dataset (actually, I should no claim it’s mine, it is inspired by the one used in the introduction of Boosting, by Robert Schapire and Yoav Freund). But in ral life, there are more observations, and more explanatory variables.With more than two explanatory variables, it starts to be more complicated to visualise. For instance, consider

MYOCARDE=read.table(
"http://freakonometrics.free.fr/saporta.csv",
head=TRUE,sep=";")

where we have observations from people in E.R., for infarctus, and we want to understand who did survive, to get a predictive model. But before running some classifier, let us visualise our data. Since we have seven explanatory variables and our class (survival or death), we can go for a PCA.

library(FactoMineR) # ACP (sur les var continues)
X=MYOCARDE[,1:7]
acp=PCA(X)

To add the death/survival variable, treat it as numerical 0/1 variable (at least to get a direction)

MYOCARDE2=MYOCARDE
MYOCARDE2$PRONO=(MYOCARDE2$PRONO=="SURVIE")*1
acp=PCA(MYOCARDE2,quanti.sup=8,graph=TRUE)

The nice thing is that we see here where variables are colinear with that one. It is also possible to visualise individuals, and classes, too

acp=PCA(MYOCARDE,quali.sup=8,graph=TRUE)
plot(acp, habillage = 8,col.hab=c("red","blue"))

# Supervised Classification, beyond the logistic

In our data-science class, after discussing limitations of the logistic regression, e.g. the fact that the decision boundary line was a straight line, we’ve mentioned possible natural extensions. Let us consider our (now) standard dataset

 clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1))
clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2))
x <- c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
y <- c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
z <- c(1,1,1,1,1,0,0,1,0,0)
df <- data.frame(x,y,z)
plot(x,y,pch=19,cex=2,col=clr1[z+1])

One can consider a quadratic function of the covariates (instead of a linear one)

 reg=glm(z~x+y+I(x^2)+I(y^2)+I(x*y),
data=df,family=binomial)
summary(reg)

pred_1 <- function(x,y){
predict(reg,newdata=data.frame(x=x,
y=y),type="response")>.5 }

x_grid<-seq(0,1,length=101)
y_grid<-seq(0,1,length=101)
z_grid <- outer(x_grid,y_grid,pred_1)
image(x_grid,y_grid,z_grid,col=clr2)
points(x,y,pch=19,cex=2,col=clr1[z+1])

# Supervised Classification, Logistic and Multinomial

We will start, in our Data Science course,  to discuss classification techniques (in the context of supervised models). Consider the following case, with 10 points, and two classes (red and blue)

> clr1 <- c(rgb(1,0,0,1),rgb(0,0,1,1))
> clr2 <- c(rgb(1,0,0,.2),rgb(0,0,1,.2))
> x <- c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
> y <- c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
> z <- c(1,1,1,1,1,0,0,1,0,0)
> df <- data.frame(x,y,z)
> plot(x,y,pch=19,cex=2,col=clr1[z+1])

To get a prediction, i.e. a partition of the space in two parts, consider some logistic regression

> reg=glm(z~x+y,data=df,family=binomial)
> summary(reg)

Call:
glm(formula = z ~ x + y, family = binomial, data = df)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.6593  -0.4400   0.2564   0.5830   1.5374

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.706      1.999  -0.854    0.393
x             -5.489      5.360  -1.024    0.306
y              8.568      5.515   1.554    0.120

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 13.4602  on 9  degrees of freedom
Residual deviance:  8.1445  on 7  degrees of freedom
AIC: 14.144

Number of Fisher Scoring iterations: 5

Given some point, the predicted class is obtained using

> pred_1 <- function(x,y){
+ predict(reg,newdata=data.frame(x=x,
+ y=y),type="response")>.5
+ }

(here, the predicted class is simply the one that is the most likely). To visualize it use

> x_grid<-seq(0,1,length=101)
> y_grid<-seq(0,1,length=101)
> z_grid <- outer(x_grid,y_grid,pred_1)
> image(x_grid,y_grid,z_grid,col=clr2)
> points(x,y,pch=19,cex=2,col=clr1[z+1])


Since the logistic regression is a (generalized) linear model, the line that separate the two regions is a straight line.

# Regression on variables, or on categories?

I admit it, the title sounds weird. The problem I want to address this evening is related to the use of the stepwise procedure on a regression model, and to discuss the use of categorical variables (and possible misinterpreations). Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db2.txt",header=TRUE,sep=";")

First, let us change the reference in our categorical variable  (just to get an easier interpretation later on)

> db$X3=relevel(as.factor(db$X3),ref="E")

If we run a logistic regression on the three variables (two continuous, one categorical), we get

> reg=glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-3.0758   0.1226   0.2805   0.4798   2.0345

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.39528    0.86649  -6.227 4.77e-10 ***
X1           0.51618    0.09163   5.633 1.77e-08 ***
X2           0.24665    0.05911   4.173 3.01e-05 ***
X3A         -0.09142    0.32970  -0.277   0.7816
X3B         -0.10558    0.32526  -0.325   0.7455
X3C          0.63829    0.37838   1.687   0.0916 .
X3D         -0.02776    0.33070  -0.084   0.9331
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 806.29  on 999  degrees of freedom
Residual deviance: 582.29  on 993  degrees of freedom
AIC: 596.29

Number of Fisher Scoring iterations: 6

Now, if we use a stepwise procedure, to select variables in the model, we get

> step(reg)
Start:  AIC=596.29
Y ~ X1 + X2 + X3

Df Deviance    AIC
- X3    4   587.81 593.81
<none>      582.29 596.29
- X2    1   600.56 612.56
- X1    1   617.25 629.25

Step:  AIC=593.81
Y ~ X1 + X2

Df Deviance    AIC
<none>      587.81 593.81
- X2    1   606.90 610.90
- X1    1   622.44 626.44

So clearly, we should remove the categorical variable if our starting point was the regression on the three variables.

Now, what if we consider the same model, but slightly different: on the five categories,

> X3complete = model.matrix(~0+X3,data=db)
> db2 = data.frame(db,X3complete)
> head(db2)
Y       X1       X2 X3 X3A X3B X3C X3D X3E
1 1 3.297569 16.25411  B   0   1   0   0   0
2 1 6.418031 18.45130  D   0   0   0   1   0
3 1 5.279068 16.61806  B   0   1   0   0   0
4 1 5.539834 19.72158  C   0   0   1   0   0
5 1 4.123464 18.38634  C   0   0   1   0   0
6 1 7.778443 19.58338  C   0   0   1   0   0

From a technical point of view, it is exactly the same as before, if we look at the regression,

> reg = glm(Y~X1+X2+X3A+X3B+X3C+X3D+X3E,family=binomial,data=db2)
> summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3A + X3B + X3C + X3D + X3E, family = binomial,
data = db2)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-3.0758   0.1226   0.2805   0.4798   2.0345

Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.39528    0.86649  -6.227 4.77e-10 ***
X1           0.51618    0.09163   5.633 1.77e-08 ***
X2           0.24665    0.05911   4.173 3.01e-05 ***
X3A         -0.09142    0.32970  -0.277   0.7816
X3B         -0.10558    0.32526  -0.325   0.7455
X3C          0.63829    0.37838   1.687   0.0916 .
X3D         -0.02776    0.33070  -0.084   0.9331
X3E               NA         NA      NA       NA
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 806.29  on 999  degrees of freedom
Residual deviance: 582.29  on 993  degrees of freedom
AIC: 596.29

Number of Fisher Scoring iterations: 6

Both regressions are equivalent. Now, what about a stepwise selection on this new model?

> step(reg)
Start:  AIC=596.29
Y ~ X1 + X2 + X3A + X3B + X3C + X3D + X3E

Step:  AIC=596.29
Y ~ X1 + X2 + X3A + X3B + X3C + X3D

Df Deviance    AIC
- X3D   1   582.30 594.30
- X3A   1   582.37 594.37
- X3B   1   582.40 594.40
<none>      582.29 596.29
- X3C   1   585.21 597.21
- X2    1   600.56 612.56
- X1    1   617.25 629.25

Step:  AIC=594.3
Y ~ X1 + X2 + X3A + X3B + X3C

Df Deviance    AIC
- X3A   1   582.38 592.38
- X3B   1   582.41 592.41
<none>      582.30 594.30
- X3C   1   586.30 596.30
- X2    1   600.58 610.58
- X1    1   617.27 627.27

Step:  AIC=592.38
Y ~ X1 + X2 + X3B + X3C

Df Deviance    AIC
- X3B   1   582.44 590.44
<none>      582.38 592.38
- X3C   1   587.20 595.20
- X2    1   600.59 608.59
- X1    1   617.64 625.64

Step:  AIC=590.44
Y ~ X1 + X2 + X3C

Df Deviance    AIC
<none>      582.44 590.44
- X3C   1   587.81 593.81
- X2    1   600.73 606.73
- X1    1   617.66 623.66

What do we get now? This time, the stepwise procedure recommends that we keep one category (namely C). So my point is simple: when running a stepwise procedure with factors, either we keep the factor as it is, or we drop it. If it is necessary to change the design, by pooling together some categories, and we forgot to do it, then it will be suggested to remove that variable, because having 4 categories meaning the same thing will cost us too much if we use the Akaike criteria. Because this is exactly what happens here

> library(car)
> reg = glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)
> linearHypothesis(reg,c("X3A=X3B","X3A=X3D","X3A=0"))
Linear hypothesis test

Hypothesis:
X3A - X3B = 0
X3A - X3D = 0
X3A = 0

Model 1: restricted model
Model 2: Y ~ X1 + X2 + X3

Res.Df Df  Chisq Pr(>Chisq)
1    996
2    993  3 0.1446      0.986

So here, we should pool together categories A, B, D and E (which was here the reference). As mentioned in a previous post, it is necessary to pool together categories that should be pulled together as soon as possible. If not, the stepwise procedure might yield to some misinterpretations.

# Logistic regression and categorical covariates

A short post to get back – for my nonlife insurance course – on the interpretation of the output of a regression when there is a categorical covariate. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)
> tail(db)
Y       X1       X2 X3
995  1 4.801836 20.82947  A
996  1 9.867854 24.39920  C
997  1 5.390730 21.25119  D
998  1 6.556160 20.79811  D
999  1 4.710276 21.15373  A
1000 1 6.631786 19.38083  A

Let us run a logistic regression on that dataset

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.45885    1.04646  -4.261 2.04e-05 ***
X1           0.51664    0.11178   4.622 3.80e-06 ***
X2           0.21008    0.07247   2.899 0.003745 **
X3B          1.74496    0.49952   3.493 0.000477 ***
X3C         -0.03470    0.35691  -0.097 0.922543
X3D          0.08004    0.34916   0.229 0.818672
X3E          2.21966    0.56475   3.930 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 397.69  on 993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

Here, the reference is modality $A$. Which means that for someone with characteristics $(X_1,X_2,X_3=A)$, we predict the following probability

$p=H(\widehat\beta_0+\widehat\beta_1 X_1+\widehat\beta_2 X_2)$

where $H(\cdot)$ denotes the cumulative distribution function of the logistic distribution

$H(x)=\frac{e^x}{1+e^x}$

For someone with characteristics $(X_1,X_2,X_3=B)$, we predict the following probability

$p=H(\widehat\beta_0+\widehat\beta_1 X_1+\widehat\beta_2 X_2+\widehat\beta_3^{\ (B)})$

For someone with characteristics $(X_1,X_2,X_3=C)$, we predict the following probability

$p=H(\widehat\beta_0+\widehat\beta_1 X_1+\widehat\beta_2 X_2+\widehat\beta_3^{\ (C)})$

(etc.) Here, if we accept $H_0:\beta_3^{\ (C)}=0$ (against $H_1:\beta_3^{\ (C)}\neq0$), it means that modality $C$ cannot be considerd as different from $A$.

A natural idea can be to change the reference modality, and to look at the $p$-values. If we consider the following loop, we get

> M = matrix(NA,5,5)
> rownames(M)=colnames(M)=LETTERS[1:5]
> for(k in 1:5){
+ db$X3 = relevel(X3,LETTERS[k]) + reg = glm(Y~X1+X2+X3,family=binomial,data=db) + M[levels(db$X3)[-1],k] = summary(reg)$coefficients[4:7,4] + } > M A B C D E A NA 0.0004771853 9.225428e-01 0.8186723647 8.482647e-05 B 4.771853e-04 NA 4.841204e-04 0.0009474491 4.743636e-01 C 9.225428e-01 0.0004841204 NA 0.7506242347 9.194193e-05 D 8.186724e-01 0.0009474491 7.506242e-01 NA 1.730589e-04 E 8.482647e-05 0.4743636442 9.194193e-05 0.0001730589 NA and if we simply want to know if the $p$-value exceeds – or not – 5%, we get the following, > M.TF = M>.05 > M.TF A B C D E A NA FALSE TRUE TRUE FALSE B FALSE NA FALSE FALSE TRUE C TRUE FALSE NA TRUE FALSE D TRUE FALSE TRUE NA FALSE E FALSE TRUE FALSE FALSE NA The first column is obtained when $A$ is the reference, and then, we see which parameter should be considered as null. The interpretation is the following: • $C$ and $D$ are not different from $A$ • $E$ is not different from $B$ • $A$ and $D$ are not different from $C$ • $A$ and $C$ are not different from $D$ • $B$ is not different from $E$ Note that we only have, here, some kind of intuition. So, let us run a more formal test. Let us consider the following regression (we remove the intercept to get a model easier to understand) > library(car) > db$X3=relevel(X3,"A")
> reg=glm(Y~0+X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
X1   0.51664    0.11178   4.622 3.80e-06 ***
X2   0.21008    0.07247   2.899  0.00374 **
X3A -4.45885    1.04646  -4.261 2.04e-05 ***
X3E -2.23919    1.06666  -2.099  0.03580 *
X3D -4.37881    1.04887  -4.175 2.98e-05 ***
X3C -4.49355    1.06266  -4.229 2.35e-05 ***
X3B -2.71389    1.07274  -2.530  0.01141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1386.29  on 1000  degrees of freedom
Residual deviance:  397.69  on  993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

It is possible to use Fisher test to test if some coefficients are equal, or not (more generally if some linear constraints are satisfied)

> linearHypothesis(reg,c("X3A=X3C","X3A=X3D","X3B=X3E"))
Linear hypothesis test

Hypothesis:
X3A - X3C = 0
X3A - X3D = 0
- X3E  + X3B = 0

Model 1: restricted model
Model 2: Y ~ 0 + X1 + X2 + X3

Res.Df Df  Chisq Pr(>Chisq)
1    996
2    993  3 0.6191      0.892

Here, we clearly accept the assumption that the first three factors are equal, as well as the last two. What is the next step? Well, if we believe that there are mainly two categories, $\{A,C,D\}$ and $\{B,E\}$, let us create that factor,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis) > reg=glm(Y~X1+X2+X3bis,family=binomial,data=db) > summary(reg) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.39439 1.02791 -4.275 1.91e-05 *** X1 0.51378 0.11138 4.613 3.97e-06 *** X2 0.20807 0.07234 2.876 0.00402 ** X3bisBE 1.94905 0.36852 5.289 1.23e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 552.64 on 999 degrees of freedom Residual deviance: 398.31 on 996 degrees of freedom AIC: 406.31 Number of Fisher Scoring iterations: 7 Here, all the categories are significant. So we do have a proper model. # Large claims, and ratemaking During the course, we have seen that it is natural to assume that not only the individual claims frequency can be explained by some covariates, but individual costs too. Of course, appropriate families should be considered to model the distribution of the cost $Y$, given some covariates $\boldsymbol{X}$.Here is the dataset we’ll use, > sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt", + header=TRUE,sep=";") > sinistres=sinistre[sinistre$garantie=="1RC",]
>  sinistres=sinistres[sinistres$cout>0,] > contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt", + header=TRUE,sep=";") > couts=merge(sinistres,contrat) > tail(couts) nocontrat no garantie cout exposition zone puissance agevehicule 1919 6104006 11933 1RC 5376.04 0.37 E 6 1 1920 6107355 12349 1RC 51.63 0.74 E 4 1 1921 6108364 13229 1RC 1320.00 0.74 B 9 1 1922 6109171 11567 1RC 1320.00 0.74 B 13 1 1923 6111208 14161 1RC 970.20 0.49 E 10 5 1924 6111650 14476 1RC 1940.40 0.48 E 4 0 ageconducteur bonus marque carburant densite region 1919 32 57 12 E 93 10 1920 45 57 12 E 72 10 1921 32 100 12 E 83 0 1922 56 50 12 E 93 13 1923 30 90 12 E 53 2 1924 69 50 12 E 93 13 Here, each line is a claim. Usual families to model the cost are the Gamma distribution, or the inverse Gaussian. Or the lognormal distribution (which is not in the exponential family, but one can assume that the logarithm of the cost can be modeled with a Gaussian distribution). Consider here only one covariate, e.g. the age of the car, and two different models: a Gamma one, and a lognormal one. > age=0:20 > reggamma.sp <- glm(cout~agevehicule,family=Gamma(link="log"), + data=couts) > Pgamma <- predict(reggamma.sp,newdata=data.frame(agevehicule=age),type="response") For the Gamma regression, it is a simple GLM, so it is not difficult. For a lognormal distribution, one should remember that the expected value of a lognormal distribution is not the exponential of the underlying Gaussian distribution. A correction should be made, here to get an unbiased estimator for the average cost, > reglm.sp <- lm(log(cout)~agevehicule,data=baseCOUT) > sigma <- summary(reglm.sp)$sigma
> mu <- predict(reglm.sp,newdata=data.frame(agevehicule=age))
> Pln <- exp(mu+sigma^2/2)

We can plot those two predictions on a single graph,

> plot(age,Pgamma,xlab="",ylab="",col="red",type="b",pch=4)
> lines(age,Pln,col="blue",type="b")

Here it is,

Observe that it is also possible to use splines, since there might be no reason for the age to appear here in a multiplicative way,

Here, the two models are rather close. Nevertheless, one should remember that the Gamma model can be extremely sensitive to large claims (I mean here really large claims). On the other hand, with the log-transformation for the lognormal model, it seams that this model is less sensitive to large events. Actually, if I use the complete dataset, the regressions are the following,

i.e. with a lognormal distribution, the average cost is decreasing with the age of the car, while it is increasing with a Gamma model. The main reason here is that there is one large (not to say huge) claim in the dataset,

> couts[which.max(couts$cout),] cout exposition zone puissance agevehicule ageconducteur 7842 4024601 0.22 B 9 13 19 marque carburant densite region 7842 2 E 93 24 One young driver got a$ 4 million claim, with a 13 year old car. This is an outliers for the Gamma regression, that clearly influences the estimation (the second largest if only one third of this one). Since there is a clear influence of large claims on the estimation of the average cost, a natural idea might be to remove those large claims. Or perhaps to see them as different from normal claims: normal claims can be explained by some covariates, but perhaps that those large claims should be shared not only within its own class, but within all the insured on the portfolio. To formalize this idea, observe that we can write

$\mathbb{E}(Y|\boldsymbol{X}) = {\color{Blue} {\underbrace{\mathbb{E}(Y|\boldsymbol{X},Y\leq s)}_{A} \cdot {\underbrace{\mathbb{P}(Y\leq s|\boldsymbol{X})}_{B}}}}+{\color{Red} {{\underbrace{\mathbb{E}(Y|Y> s, \boldsymbol{X}) }_{C}}\cdot {\underbrace{\mathbb{P}(Y> s| \boldsymbol{X})}_{B}}}}$

where the blue part is associated to normal-sized claims, while large ones correspond to the red part. It is then possible to run three regressions: one on normal sized claims, one on large claims, and one on the indicator of having a large claims, given that a claim occurred. The code here is something like that: a large claim – here – is above $10,000 (one has a fix it) > s= 10000 > couts$normal=(couts$cout<=s) > mean(couts$normal)
[1] 0.9818087

which represent 2% of the claims in our dataset.We can run 3 sets of regressions, with smoothed regression on the age of the car. The first one to model large claims individual costs,

> indice = which(couts$cout>s) > mean(couts$cout[indice])
[1] 34471.59
> library(splines)
> regB=glm(cout~bs(agevehicule),data=couts,
+ subset=indice,family=Gamma(link="log"))
> ypB=predict(regB,newdata=data.frame(agevehicule=age),type="response")
> ypB2=mean(couts$cout[indice]) the second one to model normal claims individual costs, > indice = which(couts$cout<=s)
> mean(couts$cout[indice]) [1] 1335.878 > regA=glm(cout~bs(agevehicule),data=couts, + subset=indice,family=Gamma(link="log")) > ypA=predict(regA,newdata=data.frame(agevehicule=age),type="response") > ypA2=mean(couts$cout[indice])

And finally, a third one, on the probability of having a normal sized claim, given that a claim occurred

> regC=glm(normal~bs(agevehicule),data=couts,family=binomial)
> ypC=predict(regC,newdata=data.frame(agevehicule=age),type="response")
> regC2=glm(normal~1,data=couts,family=binomial)
> ypC2=predict(regC2,newdata=data.frame(agevehicule=age),type="response")

Note that we to have, each time something that can be interpreted either as $\mathbb{E}(Y|\boldsymbol{X},Y\gtrless s)$, or $\mathbb{E}(Y|Y\gtrless s)$ – i.e. no covariate is considered on the later. On the graph below, we did plot

$\mathbb{E}(Y|\boldsymbol{X}) = {\color{Blue} {\underbrace{\mathbb{E}(Y|\boldsymbol{X},Y\leq s)}_{A} \cdot {\underbrace{\mathbb{P}(Y\leq s|\boldsymbol{X})}_{B}}}}+{\color{Red} {{\underbrace{\mathbb{E}(Y|Y> s, \boldsymbol{X}) }_{C}}\cdot {\underbrace{\mathbb{P}(Y> s| \boldsymbol{X})}_{B}}}}$

where Gamma regressions – with splines – are considered for the average costs, while logistic regressions – again with splines – are considered to model probabilities.

(but careful with splines: on borders, since we do not have a lot of observations, the behavior can be… odd. And adjustments should be made to obtain an adequate level of premium).  If it is legitimate to assume that normal-sized claims can be explained by some covariates, perhaps large claims (or extremely large ones) are just purely random, i.e. not function of any covariate, at all. I.e.

$\mathbb{E}(Y|\boldsymbol{X}) = {\color{Blue} {\underbrace{\mathbb{E}(Y|\boldsymbol{X},Y\leq s)}_{A} \cdot {\underbrace{\mathbb{P}(Y\leq s|\boldsymbol{X})}_{B}}}}+{\color{Red} {{\underbrace{\mathbb{E}(Y|Y> s) }_{C'}}\cdot {\underbrace{\mathbb{P}(Y> s| \boldsymbol{X})}_{B}}}}$

To go one step further, it might also be possible to assume that not only the size of the claim (given that it is a large one) is not a function of any covariate, but perhaps neither is the probability of having an extremely large claim, too

$\mathbb{E}(Y|\boldsymbol{X}) = {\color{Blue} {\underbrace{\mathbb{E}(Y|\boldsymbol{X},Y\leq s)}_{A} \cdot {\underbrace{\mathbb{P}(Y\leq s)}_{B'}}}}+{\color{Red} {{\underbrace{\mathbb{E}(Y|Y> s) }_{C'}}\cdot {\underbrace{\mathbb{P}(Y> s)}_{B'}}}}$

From the first part, we’ve seen that the distribution considered had an impact on the prediction, and in the second, we’ve seen that the definition of large claims (and how to deal with them) also has an impact. So clearly, actuaries have some leverage when working on ratemaking…