Next week, with , I will be at the Summer School of the Swiss Association of Actuaries, in Lausanne, with Jean-Philippe Boucher (UQAM) and Ewen Gallic (AMSE).
There will be some hands-on applications, on R. I will share some codes in the slides.
Next week, with , I will be at the Summer School of the Swiss Association of Actuaries, in Lausanne, with Jean-Philippe Boucher (UQAM) and Ewen Gallic (AMSE).
There will be some hands-on applications, on R. I will share some codes in the slides.
This article presents a set of tools for the modeling of a spatial allocation problem in a large geographic market and gives examples of applications. In our settings, the market is described by a network that maps the cost of travel between each pair of adjacent locations. Two types of agents are located at the nodes of this network. The buyers choose the most competitive sellers depending on their prices and the cost to reach them. Their utility is assumed additive in both these quantities. Each seller, taking as given other sellers prices, sets her own price to have a demand equal to the one we observed. We give a linear programming formulation for the equilibrium conditions. After formally introducing our model we apply it on two examples: prices offered by petrol stations and quality of services provided by maternity wards (only the later is described here for privacy issues). These examples illustrate the applicability of our model to aggregate demand, rank prices and estimate cost structure over the network. We insist on the possibility of applications to large scale data sets using modern linear programming solvers such as Gurobi.
Demand for gas in gas stations in Britanny, and demand for maternity in France (with border correction)
In addition to this paper we released a R toolbox to implement our results and an online tutorial, optimalnetwork.github.io.
Before moving to cross-validation, it was natural to say “I will burn 50% (say) of my data to train a model, and then use the remaining to fit the model”. For instance, we can use training data for variable selection (e.g. using some stepwise procedure in a logistic regression), and then, once variable have been selected, fit the model on the remaining set of observations. A natural question is usually “does it really matter ?”.
In order to visualize this problem, consider my (simple) dataset
MYOCARDE=read.table( "http://freakonometrics.free.fr/saporta.csv", head=TRUE,sep=";") |
Let us generate 100 training samples (where we keep about 50% of the observations). On each of them, we use a stepwise procedure, and we keep the estimates of the remaining variables (and their standard deviation actually)
n=nrow(MYOCARDE) M=matrix(NA,100,ncol(MYOCARDE)) colnames(M)=c("(Intercept)",names(MYOCARDE)[1:7]) S1=S2=M1=M2=M for(i in 1:100){ idx = which(sample(0:1,size=n, replace=TRUE)==1) reg=step(glm(PRONO=="DECES"~.,data=MYOCARDE[idx,])) nm=names(reg$coefficients) M1[i,nm]=reg$coefficients S1[i,nm]=summary(reg)$coefficients[,2] f=paste("PRONO=='DECES'~",paste(nm[-1],collapse="+"),sep="") reg=glm(f,data=MYOCARDE[-idx,]) M2[i,nm]=reg$coefficients S2[i,nm]=summary(reg)$coefficients[,2] } |
Then, for the 7 covariates (and the constant) we can look at the value of the coefficient in the model fitted on the training sample, and the value on the model fitted on the validation sample (of course, only when they were remaining)
for(j in 1:8){ idx=which(!is.na(M1[,j])) plot(M1[idx,j],M2[idx,j]) abline(a=0,b=1,lty=2,col="gray") segments(M1[idx,j]-2*S1[idx,j],M2[idx,j],M1[idx,j]+2*S1[idx,j],M2[idx,j]) segments(M1[idx,j],M2[idx,j]-2*S2[idx,j],M1[idx,j],M2[idx,j]+2*S2[idx,j]) } |
For instance, with the intercept, we have the following

where horizontal segments are confidence intervals of the parameter on the model fitted on the training sample, the vertical on the validation sample. The green part means some sort of consistency, while the red one means that actually, the coefficient was negative with one model, positive with the other one. Which is odd (but in that case, observe that coefficients are rarely significant).
We can also visualize the joint distribution of the two estimators,
for(j in 1:8){ library(ks) idx = which(!is.na(M1[,j])) Z = cbind(M1[idx,j],M2[idx,j]) H = Hpi(x=Z) fhat = kde(x=Z, H=H) image(fhat$eval.points[[1]], fhat$eval.points[[2]],fhat$estimate) abline(a=0,b=1,lty=2,col="gray") abline(v=0,lty=2) abline(h=0,lty=2) } |
which are here, almost on the diagonal,

meaning that the intercept on the two samples is (more or less) the same. We can then look at other parameters (which is actually more interesting).


On that variable, it seems that it is significant on the training dataset (somehow, it is consistent with the fact that it is remaining in the model after the stepwise procedure) but not on the validation sample (or hardly significant).
Others are much more consistent (with some possible outliers)






On the next one, we have again significance on the training sample, but not on the validation sample,




and probably more interesting


where the two are very consistent.
In my previous post on GLMs, I discussed power link functions. But there are much more links that can be used :

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.
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>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>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(>|z|) (Intercept) -0.002274 0.048413 -0.047 0.963 X1 1.019814 0.026186 38.945 <2e-16 *** X2 1.004814 0.024172 41.570 <2e-16 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|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…
Yesterday, I published a post on extracting information from a picture, but it did not work as expected. I claimed that it was because of the original graph I had. More precisely, the was based on some weird projection, and I could not reconcile. So I decide to cheat a little bit, by creating my own map,

Colors are ugly, I know. But I got them using
u = seq(0,1,length=30) couleurs = rgb(u,rev(u),0,1) |
The picture is
url = "https://freakonometrics.hypotheses.org/files/2018/12/chomage3.png" library(pixmap) library(png) IMG = readPNG(url) |
I used those colors because it would make things easy when extracting reds and greens…
ROUGE=t(IMG[,,1])[x1:x2,] ROUGE=ROUGE[,y2:y1] library(scales) image(x1:x2,y1:y2,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) VERT=t(IMG[,,2])[x1:x2,] VERT=VERT[,y2:y1] image(x1:x2,y1:y2,VERT,col=alpha(colour=rgb(0,1,0,1), alpha = seq(0,1,by=.01))) |

Let us see if the contour of France can be overlaid
library(maptools) library(PBSmapping) download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds") FR=readRDS("FRA_adm0.rds") library(maptools) PP = SpatialPolygons2PolySet(FR) par(mfrow=c(1,1)) PP=PP[(PP$X<=8.25)&(PP$Y>=42.2),] u=(x1:x2)-x1 v=(y1:y2)-y1 ax=min(PP$X) bx=max(PP$X)-min(PP$X) ay=min(PP$Y) by=max(PP$Y)-min(PP$Y) PP$X=(PP$X-ax)/bx*max(u) PP$Y=(PP$Y-ay)/by*max(v) image(u,v,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) points(PP$X,PP$Y) |
We have a perfect match, don’t we…?

Let us now use a shapefile based on départements,
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm2.rds","FRA_adm2.rds") FR2=readRDS("FRA_adm2.rds") library(maptools) PP = SpatialPolygons2PolySet(FR2) image(u,v,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) k=35 pX=(PP$X[PP$PID==k]-ax)/bx*max(u) pY=(PP$Y[PP$PID==k]-ay)/by*max(v) points(pX,pY)nge(pX) |
For instance, the thirty-fifth polygon is the following

Let us extract the color inside that polygon
u=1:nrow(ROUGE) v=1:ncol(ROUGE) |
The code would be
pX=(PP$X[PP$PID==k]-ax)/bx*max(u) pY=(PP$Y[PP$PID==k]-ay)/by*max(v) E=expand.grid(u,v) M=matrix(point.in.polygon(E[,1],E[,2],pX,pY)>0,length(u),length(v)) image(u,v,ROUGE*M,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) points(pX,pY) |

Now, for each département, I extract the average value of red, and the average value of green,
extract_info = function(k){ pX=(PP$X[PP$PID==k]-ax)/bx*max(u) pY=(PP$Y[PP$PID==k]-ay)/by*max(v) E=expand.grid(u,v) M=matrix(point.in.polygon(E[,1],E[,2],pX,pY)>0,length(u),length(v)) nom=FR2[FR2$OBJECTID ==k,c("NAME_2","CCA_2")] return(c(as.numeric(nom$CCA_2),sum(ROUGE[M==1])/sum(M),sum(VERT[M==1])/sum(M))) } donnees = Vectorize(extract_info)(1:95) x2=donnees[1,] y2=donnees[2,]/(donnees[2,]+donnees[3,]) df2=data.frame(dpt=x2,extract=y2) x1=as.numeric(as.character(baseChomage$no)) y1=baseChomage$chomagePremierTrimestre2017 df1=data.frame(dpt=x1,obs=y1) df=merge(df1,df2) plot(df$obs,df$extract) |
On the graph below, we have the original values on the x-axis (unemployement, in percent) and the “average value of red”. Note that points are almost perfectly correlated… The accumulation can be explained because on the original map, different values could have the same color

So far, I can claim that we’ve been able to extract useful information from the original picture.
Consider the case now that the original map was the following one

The picture can be downloaded using the following code
url = "https://freakonometrics.hypotheses.org/files/2018/12/chomage5.png" library(pixmap) library(png) IMG = readPNG(url) |
Here, the colors are obtained from a standard palette,
library(pals) couleurs = rev(brewer.rdylgn(30)) |
Here again, we use our previous code to extract reds and greens

And if we use our function
extract_info = function(k){ pX=(PP$X[PP$PID==k]-ax)/bx*max(u) pY=(PP$Y[PP$PID==k]-ay)/by*max(v) E=expand.grid(u,v) M=matrix(point.in.polygon(E[,1],E[,2],pX,pY)>0,length(u),length(v)) nom=FR2[FR2$OBJECTID ==k,c("NAME_2","CCA_2")] return(c(as.numeric(nom$CCA_2),sum(ROUGE[M==1])/sum(M),sum(VERT[M==1])/sum(M))) } donnees = Vectorize(extract_info)(1:95) x2=donnees[1,] y2=donnees[2,]/(donnees[2,]+donnees[3,]) df2=data.frame(dpt=x2,extract=y2) x1=as.numeric(as.character(baseChomage$no)) y1=baseChomage$chomagePremierTrimestre2017 df1=data.frame(dpt=x1,obs=y1) df=merge(df1,df2) plot(df$obs,df$extract) |
we obtain the following graph

Here again, we have a strong correlation, not to say comonotonic variables (in the sense that ranks are identical). Nice, isn’t it ?
This week, I wanted to get information I found on the nice map, below. I could not get access to the original dataset, per zip code… and I was wondering, if (assuming that the map was with high resolution) it was actually possible to extract information, using a simple R function…

As we can see, there is red, and green on the map, and I would love to know which are the green and the red cities, in France. One important issue is actually the background. Here it’s nice, it white… but white is a strange color, achromatic and very light. More specifically, if I search red areas, the background is very red. And very green, too. So, to avoid those issues, I did use gimp to change the background, into black. On the opposite, where it’s black, it’s neither red, nor green !

Let us get the map, and extract information from the file
url="https://freakonometrics.hypotheses.org/files/2018/12/inondation3.png" download.file(url,"inondation3.png") image="inondation3.png" library(pixmap) library(png) IMG=readPNG(image)
Information is stored in several matrices – or in arrays. Dimension 1 is the height of the picture (in pixels), dimension 2 is the width, and the third one is either 1 (red), 2 (green) or 3 (blue), based on the rgb decomposition of each pixel. Then, I try to find the border of the map
nl=dim(IMG)[1] nc=dim(IMG)[2] MAT=(IMG[,,1]+IMG[,,2])/2 x=apply(MAT,2,max) plot(x,type="l")

When it’s null, it means no color on the line of the matrix, i.e. completly black (initially, I used the mean function, but the maximum really behaves like a step function)
y=apply(MAT,1,max) plot(y,type="l")

Let us find cutoff values, on the left and on the right, on top and on the bottom
image(1:nc,1:nl,t(MAT)) abline(v=min(which(x>.2)),col="blue") abline(v=max(which(x>.2)),col="blue") abline(h=min(which(y>.2)),col="blue") abline(h=max(which(y>.2)),col="blue")
We obtain the following (forget about the fact that – somehow – France is upside-down)

We can zoom-in, just to make sure that our border are fine
par(mfrow=c(1,2)) image(min(which(x>.2))+(-5):5,1:nl,t(MAT)[min(which(x>.2))+(-5):5,]) abline(v=min(which(x>.2))+(-5):5,col="white") abline(v=min(which(x>.2)),col="blue") x1=min(which(x>.2))-1

and on the vertical range
image(max(which(x>.2))+(-5):5,1:nl,t(MAT)[max(which(x>.2))+(-5):5,]) abline(v=max(which(x>.2))+(-5):5,col="white") abline(v=max(which(x>.2)),col="blue") x2=max(which(x>.2))+1

So far so good. Let us keep the subpart of the picture,
image(x1:x2,y1:y2,t(MAT)[x1:x2,y1:y2])

Now, let us focus on the red part / component of that picture
ROUGE=t(IMG[,,1])[x1:x2,] ROUGE=ROUGE[,y2:y1] library(scales) image(x1:x2,y1:y2,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))

That’s not bad, isn’t it ? And get can have a similar graph for the green part
VERT=t(IMG[,,2])[x1:x2,] VERT=VERT[,y2:y1] image(x1:x2,y1:y2,VERT,col=alpha(colour=rgb(0,1,0,1), alpha = seq(0,1,by=.01)))

Now, I wanted to ajust a map of France on that one. Using shapefiles of administrative regions, it would be possible to get the proportion of red and green parts (départements, cantons, etc). As a starting point (before going to ‘départements’), let us use a standard shapefile for France
library(maptools)
library(PBSmapping)
url="http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds"
download.file(url,"FRA_adm0.rds")
FR=readRDS("FRA_adm0.rds")
library(maptools)
PP = SpatialPolygons2PolySet(FR)
PP=PP[(PP$X<=8.25)&(PP$Y>=42.2),]
u=(x1:x2)-x1
v=(y1:y2)-y1
ax=min(PP$X)
bx=max(PP$X)-min(PP$X)
ay=min(PP$Y)
by=max(PP$Y)-min(PP$Y)
PP$X=(PP$X-ax)/bx*max(u)
PP$Y=(PP$Y-ay)/by*max(v)
image(u,v,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01)))
points(PP$X,PP$Y)
We try here to rescale it. The left part should be on the left part of the picture as well as the right part. And the same holds for the top, and the bottom,

Unfortunately, even if we change the projection technique, I could not match perfectly the contour of France. I am quite sure that it’s a projection problem ! But I did try a dozen popular ones, with no success… so if anyone has a clever idea…
Usually, when I give a course on GLMs, I try to insist on the fact that the link function is probably more important than the distribution. In order to illustrate, consider the following dataset, with 5 observations
x = c(1,2,3,4,5) y = c(1,2,4,2,6) base = data.frame(x,y) |
Then consider several model, with various distributions, and either an identity link (and in that case \mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbf{x}^T\mathbf{\beta}) or a log link function (so that \mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=e^{\mathbf{x}^T\mathbf{\beta}})
regNId = glm(y~x,family=gaussian(link="identity"),data=base) regNlog = glm(y~x,family=gaussian(link="log"),data=base) regPId = glm(y~x,family=poisson(link="identity"),data=base) regPlog = glm(y~x,family=poisson(link="log"),data=base) regGId = glm(y~x,family=Gamma(link="identity"),data=base) regGlog = glm(y~x,family=Gamma(link="log"),data=base) regIGId = glm(y~x,family=inverse.gaussian(link="identity"),data=base) regIGlog = glm(y~x,family=inverse.gaussian(link="log"),data=base |
One can also consider some Tweedie distribution, to be even more general
library(statmod) regTwId = glm(y~x,family=tweedie(var.power=1.5,link.power=1),data=base) regTwlog = glm(y~x,family=tweedie(var.power=1.5,link.power=0),data=base) |
Consider the prediction obtained in the first case, with the linear link function
library(RColorBrewer) darkcols = brewer.pal(8, "Dark2") plot(x,y,pch=19) abline(regNId,col=darkcols[1]) abline(regPId,col=darkcols[2]) abline(regGId,col=darkcols[3]) abline(regIGId,col=darkcols[4]) abline(regTwId,lty=2) |

The predictions are very very close, aren’t they ? In the case of the exponential prediction, we obtain
plot(x,y,pch=19) u=seq(.8,5.2,by=.01) lines(u,predict(regNlog,newdata=data.frame(x=u),type="response"),col=darkcols[1]) lines(u,predict(regPlog,newdata=data.frame(x=u),type="response"),col=darkcols[2]) lines(u,predict(regGlog,newdata=data.frame(x=u),type="response"),col=darkcols[3]) lines(u,predict(regIGlog,newdata=data.frame(x=u),type="response"),col=darkcols[4]) lines(u,predict(regTwlog,newdata=data.frame(x=u),type="response"),lty=2) |

We can actually look closer. For instance, in the linear case, consider the slope obtained with a Tweedie model (that will include all the parametric familes mentioned here, actually)
pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base))$coefficients[2,1:2] Vgamma = seq(-.5,3.5,by=.05) Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3,ylim=c(.965,1.03),xlab="power",ylab="slope") |

The slope here is always very very close to one ! Even more if we add a confidence interval
plot(Vgamma,Vpente[1,]) lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2) |

Heuristically, for the Gamma regression, or the Inverse Gaussian one, because the variance is a power of the prediction, if the prediction is small (here on the left), the variance should be small. So, on the left of the graph, the error should be small with a higher power for the variance function. And that’s indeed what we observe here
erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(-.1,.04),xlab="power",ylab="error") abline(h=0,lty=2) |

Of course, we can do the same with the exponential models
pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base))$coefficients[2,1:2] Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3) |

or, if we add the confidence bands, we obtain
plot(Vgamma,Vpente[1,],ylim=c(0,.8),type="l",lwd=3,xlab="power",ylab="slope") lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2) |

So here also, the “slope” is rather similar… And if we look at the error we make on the left part of the graph, we obtain
erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(.001,.32),xlab="power",ylab="error") |

So my point is that the distribution is usually not the most important point on GLMs, even if chapters of books on GLMs are distribution based… But as mentioned in an another post, if you consider a nonlinear transformation, like we have with GAMs, the story is more complicated…
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) [1] 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_jFrom 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 1 B [1] 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,]) [1] 26.48313 sum(B*E[2,]) [1] 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]) [1] 107 sum(A*E[,2]) [1] 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(>|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) [1] 355 sum(df$y) [1] 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)[1]+c(0,coefficients(reg)[2:5])) a/a[1] religion2 religion3 religion4 religion5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 A/A[1] 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)[1]+c(0,coefficients(reg)[6:11])) b/b[1] occupation2 occupation3 occupation4 occupation5 occupation6 occupation7 1.0000000 0.5612551 1.8043769 1.5463210 1.0435773 1.6584203 3.5756477 B/B[1] 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).
Real-time computation (or estimation) of the “probability to win” is difficult. We’ve seem that in soccer games, in elections… but actually, as a professor, I see that frequently when I grade my students.
Consider a classical multiple choice exam. After each question, imagine that you try to compute the probability that the student will pass. Consider here the case where we have 50 questions. Students pass when they have 25 correct answers, or more. Just for simulations, I will assume that students just flip a coin at each question… I have n students, and 50 questions
set.seed(1) n=10 M=matrix(sample(0:1,size=n*50,replace=TRUE),50,n) |
Let X_{i,j} denote the score of student i at question j. Let S_{i,j} denote the cumulated score, i.e. S_{i,j}=X_{i,1}+\cdots+X_{i,j}. At step j, I can get some sort of prediction of the final score, using \hat{T}_{i,j}=50\times S_{i,j}/j. Here is the code
SM=apply(M,2,cumsum) NB=SM*50/(1:50) |
We can actually plot it
plot(NB[,1],type="s",ylim=c(0,50)) abline(h=25,col="blue") for(i in 2:n) lines(NB[,i],type="s",col="light blue") lines(NB[,3],type="s",col="red") |

But that’s simply the prediction of the final score, at each step. That’s not the computation of the probability to pass !
Let’s try to see how we can do it… If after j questions, the students has 25 correct answer, the probability should be 1 – i.e. if S_{i,j}\geq 25 – since he cannot fail. Another simple case is the following : if after j questions, the number of points he can get with all correct answers until the end is not sufficient, he will fail. That means if S_{i,j}+(50-i+1)< 25 the probability should be 0. Otherwise, to compute the probability to sucess, it is quite straightforward. It is the probability to obtain at least 25-S_{i,j} correct answers, out of 50-j questions, when the probability of success is actually S_{i,j}/j. We recognize the survival probability of a binomial distribution. The code is then simply
PB=NB*NA for(i in 1:50){ for(j in 1:n){ if(SM[i,j]>=25) PB[i,j]=1 if(SM[i,j]+(50-i+1)<25) PB[i,j]=0 if((SM[i,j]<25)&(SM[i,j]+(50-i+1)>=25)) PB[i,j]=1-pbinom(25-SM[i,j],size=(50-i),prob=SM[i,j]/i) }} |
So if we plot it, we get
plot(PB[,1],type="s",ylim=c(0,1)) abline(h=25,col="red") for(i in 2:n) lines(PB[,i],type="s",col="light blue") lines(PB[,3],type="s",col="red") |

which is much more volatile than the previous curves we obtained ! So yes, computing the “probability to win” is a complicated exercice ! Don’t blame those who try to find it hard to do !
Of course, things are slightly different if my students don’t flip a coin… this is what we obtain if half of the students are good (2/3 probability to get a question correct) and half is not good (1/3 chance),

If we look at the probability to pass, we usually do not have to wait until the end (the 50 questions) to know who passed and who failed

PS : I guess a less volatile solution can be obtained with a Bayesian approach… if I find some spare time this week, I will try to code it…
In 2012, Danielle Herbert, Adrian Barnett, Philip Clarke and Nicholas Graves published an article entitled “on the time spent preparing grant proposals: an observational study of Australian researchers“, whose conclusions had been included in Nature under a more explicit title, “Australia’s grant system wastes time” ! In this study, they included 3700 grant applications sent to the National Health and Medical Research Council, and showed that each application represented 37 working days: “Extrapolating this to all 3,727 submitted proposals gives an estimated 550 working years of researchers’ time (95% confidence interval, 513-589)“. But in these times when I have to write my funding application, I find that losing 37 days of work is huge. Because it’s become the norm! And somehow, it’s sad.

Forget about the crazy idea that I would rather, in fact, spend more time doing my research. In fact, the thought I had this morning was that it is rather sad that in the Faculty of Science, mathematicians are asked to spend a considerable amount of time, comparable to that required of physicists or chemists, for often smaller amounts of funding… And I thought it could be easily verified. We start by retrieving the discipline codes
url="http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSC-ResultatsCSS_eng.asp" download.file(url,destfile = "GSC.html") library(XML) tables=readHTMLTable("GSC.html") GSC=tables[[1]]$V1 GSC=as.character(GSC[-(1:2)]) namesGSC=tables[[1]]$V2 namesGSC=as.character(namesGSC[-(1:2)]) |
We’re going to need a small function, to remove the $ and other symbols that pollute the data (and prevent them from being treated as numbers)
library(stringr) Correction = function(x) as.numeric(gsub('[$,]', '', x)) |
We will now read the 12 pages, and harvest (we will just take the 2017 data, but we could go back a few years before)
grants= function(gsc){ url=paste("http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSCDetail-ResultatsCSSDetails_eng.asp?Year=2017&GSC=",gsc,sep="") download.file(url,destfile = "GSC.html") library(XML) tables=readHTMLTable("GSC.html") X=as.character(tables[[1]]$"Awarded Amount") A=as.numeric(Vectorize(Correction)(X)) return(c(median(A),mean(A),as.numeric(quantile(A,(1:99)/100)))) } M=Vectorize(grants)(GSC[1:12]) |
The average amounts of individual grants can be compared,
barplot(M[2,]) |

In mathematics, the average grant amount is $24400. If we normalize by this quantity, we obtain
barplot(M[2,]/M[2,8]) |

In other words, the average amount of a (individual) grant in chemistry (to pay for students, conferences, etc.) is twice that in mathematics, 60% higher in physics than in maths…
We can also look at the median values (rather than the averages)
barplot(M[1,]) |

Here again, it is in mathematics that it is the weakest….
barplot(M[1,]/M[1,8]) |

in comparable proportions. If we think that the time spent writing should be proportional to the amount allocated, we should spend half as much time in math as in chemistry.
Cumulative functions can also be ploted,
plot(M[3:101,8],(1:99)/100,type="s",xlim=range(M)) lines(M[3:101,5],(1:99)/100,type="s",col="red") lines(M[3:101,4],(1:99)/100,type="s",col="blue") |

with math in black, physics in red, and chemistry in blue. What is surprising is the bottom part: a “bad” researcher in chemistry or physics will earn more than the median researcher in mathematics…
Now that my intuition is confirmed, I have to go back, writing my proposal… and explain to my coauthors that I have to postpone some research projects because, well, you know…
Each time we face real applications in an applied econometrics course, we have to deal with categorial variables. And the same question arise, from students : how can we combine automatically factor levels ? Is there a simple R function ?
I did upload a few blog posts, over the pas years. But so far, nothing satistfying. Let me write down a few lines about what could be done. And if some wants to write a nice R function, that would be awesome. To illustrate the idea, consider the following (simulated dataset)
n=200 set.seed(1) x1=runif(n) x2=runif(n) y=1+2*x1-x2+rnorm(n,0,.2) LB=sample(LETTERS[1:10]) b=data.frame(y=y,x1=x1, x2=cut(x2,breaks= c(-1,.05,.1,.2,.35,.4,.55,.65,.8,.9,2), labels=LB)) str(b) 'data.frame': 200 obs. of 3 variables: $ y : num 1.345 1.863 1.946 2.481 0.765 ... $ x1: num 0.266 0.372 0.573 0.908 0.202 ... $ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ... table(b$x2)[LETTERS[1:10]] A B C D E F G H I J 11 12 23 34 23 36 12 32 3 14 |
There is one (continuous) dependent variable y, one continuous covariable x_1 and one categorical variable x_2, with here ten levels. We can plot the data using
plot(b$x1,y,col="white",xlim=c(0,1.1)) text(b$x1,y,as.character(b$x2),cex=.5) |

The output of a linear regression yield the following predictions
for(i in 1:10){ p=function(x) predict(lm(y~x1+x2,data=b),newdata=data.frame(x1=x,x2=LETTERS[i])) u=seq(-1,1.065,by=.01) v=Vectorize(p)(u) lines(u,v)} |

the slope for x_1 is the same, we simply add a different constant for each level. As we can see, some levels are very very close, so it seems legitimate to combine them into one single category. Here is the output of the linear regression,
summary(lm(y~x1+x2,data=b)) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.843802 0.119655 7.052 3.23e-11 *** x1 1.992878 0.053838 37.016 < 2e-16 *** x2A 0.055500 0.131173 0.423 0.6727 x2H 0.009293 0.121626 0.076 0.9392 x2F -0.177002 0.121020 -1.463 0.1452 x2B -0.218152 0.130192 -1.676 0.0955 . x2D -0.206970 0.121294 -1.706 0.0896 . x2G -0.407417 0.129999 -3.134 0.0020 ** x2C -0.526708 0.123690 -4.258 3.24e-05 *** x2J -0.664281 0.128126 -5.185 5.54e-07 *** x2E -0.816454 0.123625 -6.604 3.94e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2014 on 189 degrees of freedom Multiple R-squared: 0.8995, Adjusted R-squared: 0.8942 F-statistic: 169.1 on 10 and 189 DF, p-value: < 2.2e-16 AIC(lm(y~x1+x2,data=b)) [1] -60.74443 BIC(lm(y~x1+x2,data=b)) [1] -21.16463 |
Here the reference category is “I”. And it looks like we could actually combine that category with several others. One strategy here would be to select all categories that seem to be not significantly different, and to run a (multiple) test
library(car) linearHypothesis(lm(y~x1+x2,data=b), c("x2A = 0", "x2H = 0", "x2F = 0")) Hypothesis: x2A = 0 x2H = 0 x2F = 0 Model 1: restricted model Model 2: y ~ x1 + x2 Res.Df RSS Df Sum of Sq F Pr(>F) 1 192 8.4651 2 189 7.6654 3 0.79971 6.5726 3e-04 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
It seems that we can combine those four categories together.
Here, we can see what’s going on when we change the reference category (actually, loop on all categories)
P=matrix(NA,nlevels(b$x2),nlevels(b$x2)) colnames(P)=rownames(P)=LETTERS[1:10] plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5), ylim=c(0,10.5)) text(1:10,0,LETTERS[1:10]) text(0,1:10,LETTERS[1:10]) for(i in 1:nlevels(b$x2)){ #levels(b$x2)=LETTERS[1:10] b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,3) P[LETTERS[i],names(p)]=p p=P[LETTERS[i],] idx=which(p>.05) points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2) idx=which(p>.1) points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)} |

We are glad to see that it is symmetric : if “H” should be combined with “I”, “I” should also be combined with “H”.
Here black points are related with the 10% p-value, and white points the 5% p-value. This graph is actually hard to read… And actually, this reminds us of Bertin (1967).

Here, we can predefine manually some ordering (we will see below how it might be automatised)
LETTERSord=c("I","A","H","F","B","D","G","C","J","E") P=matrix(NA,nlevels(b$x2),nlevels(b$x2)) colnames(P)=rownames(P)=LETTERSord plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5), ylim=c(0,10.5)) ct=c(3,3,2,1,1) abline(v=.5+c(0,cumsum(ct)),lty=2) abline(h=.5+c(0,cumsum(ct)),lty=2) text(1:10,0,LETTERSord) text(0,1:10,LETTERSord) for(i in 1:nlevels(b$x2)){ #levels(b$x2)=LETTERS[1:10] b$x2=relevel(b$x2,LETTERSord[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,3) P[LETTERSord[i],names(p)]=p p=P[LETTERSord[i],] idx=which(p>.05) points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2) idx=which(p>.1) points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2) } |
Here we get the following

It looks like we have our combined categories…
Actually, it is possible to use another strategy. We start from some level, say “A”. Then, we merge it with all non-significantly different levels. If “B” is not one of them, we use it as the new reference. Etc.
for(i in 1:nlevels(b$x2)){ if(LETTERS[i]%in%levels(b$x2)){ b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,nchar(p)) idx=which(p>.05) mix=c(LETTERS[i],names(p)[idx]) b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep="")) }} |
The final categories are
table(b$x2) A+I+H B+D+F C+G E J 46 82 35 23 14 |
with the following regression output
summary(lm(y~x1+x2,data=b)) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.86407 0.03950 21.877 < 2e-16 *** x1 1.99180 0.05323 37.417 < 2e-16 *** x2B+D+F -0.21517 0.03699 -5.817 2.44e-08 *** x2C+G -0.50545 0.04528 -11.164 < 2e-16 *** x2E -0.83617 0.05128 -16.305 < 2e-16 *** x2J -0.68398 0.06131 -11.156 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2008 on 194 degrees of freedom Multiple R-squared: 0.8975, Adjusted R-squared: 0.8948 F-statistic: 339.6 on 5 and 194 DF, p-value: < 2.2e-16 AIC(lm(y~x1+x2,data=b)) [1] -66.76939 BIC(lm(y~x1+x2,data=b)) [1] -43.68117 |
Which is consistent with the group we got before. But actually, if we change the order, we can get different combinations. For instance, if we go from “J” to “A”, instead of “A” to “J”, we obtain
for(i in nlevels(b$x2):1){ #levels(b$x2)=LETTERS[1:10] if(LETTERS[i]%in%levels(b$x2)){ b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,nchar(p)) idx=which(p>.05) mix=c(LETTERS[i],names(p)[idx]) b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep="")) }} table(b$x2) E G+C I+A+B+D+F+H J 23 35 128 14 |
with different information criteria here
AIC(lm(y~x1+x2,data=b)) [1] -36.61665 BIC(lm(y~x1+x2,data=b)) [1] -16.82675 |
I guess it would be necessary to run randomly the order we go through the levels. Last, but not least, one can use regression trees (even if it not per se in the syllabus of the course). The problem is that there is another explanatory variable that might interphere. So I would suggest (1) to fit a linear model y=\beta_0+\beta_1x_1+u_i, to calculate the residuals, \widehat{u}_i (2) to run a regression tree, to explain \widehat{u}_i with categorical variable x_2 (I did explain how trees are build when the explanatory variable is a categorical one in a previous post)
library(rpart) library(rpart.plot) b$e=residuals(lm(y~x1,data=b)) arbre=rpart(e~x2,data=b) prp(arbre,type=2,extra=1) |

Observe that the leaves have the same groups as the one we got.
arbre n= 200 node), split, n, deviance, yval * denotes terminal node 1) root 200 22.563500 7.771561e-18 2) x2=G,C,J,E 72 4.441495 -3.232525e-01 4) x2=J,E 37 1.553520 -4.578492e-01 * 5) x2=G,C 35 1.509068 -1.809646e-01 * 3) x2=I,A,H,F,B,D 128 6.366628 1.818295e-01 6) x2=F,B,D 82 2.983381 1.048246e-01 * 7) x2=I,A,H 46 2.030229 3.190993e-01 * |
I guess that it should be possible to put all that in an R function, to suggest combinations of level that might improve the regression.
This morning during the lecture on nonlinear regression, I mentioned (very) briefly the case of convex regression. Since I forgot to mention the codes in R, I will publish them here. Assume that y_i=m(\mathbf{x}_i)+\varepsilon_i where m:\mathbb{R}^d\rightarrow \mathbb{R} is some convex function.
Then m is convex if and only if \forall\mathbf{x}_1,\mathbf{x}_2\in\mathbb{R}^d, \forall t\in[0,1], m(t\mathbf{x}_1+[1-t]\mathbf{x}_2) \leq tm(\mathbf{x}_1)+[1-t]m(\mathbf{x}_2)Hidreth (1954) proved that if m^\star=\underset{m \text{ convex}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-m(\mathbf{x_i})\big)^2\right\rbracethen \mathbf{\theta}^\star=(m^\star(\mathbf{x_1}),\cdots,m^\star(\mathbf{x_n})) is unique.
Let \mathbf{y}=\mathbf{\theta}+\mathbf{\varepsilon}, then \mathbf{\theta}^\star=\underset{\mathbf{\theta}\in \mathcal{K}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-\theta_i)\big)^2\right\rbracewhere\mathcal{K}=\{\mathbf{\theta}\in\mathbb{R}^n:\exists m\text{ convex },m(\mathbf{x}_i)=\theta_i\}. I.e. \mathbf{\theta}^\star is the projection of \mathbf{y} onto the (closed) convex cone \mathcal{K}. The projection theorem gives existence and unicity.
For convenience, in the application, we will consider the real-valued case, m:\mathbb{R}\rightarrow \mathbb{R}, i.e. y_i=m(x_i)+\varepsilon_i. Assume that observations are ordered x_1\leq x_2\leq\cdots \leq x_n. Here \mathcal{K}=\left\lbrace\mathbf{\theta}\in\mathbb{R}^n:\frac{\theta_2-\theta_1}{x_2-x_1}\leq \frac{\theta_3-\theta_2}{x_3-x_2}\leq \cdots \leq \frac{\theta_n-\theta_{n-1}}{x_n-x_{n-1}}\right\rbrace
Hence, quadratic program with n-2 linear constraints.
m^\star is a piecewise linear function (interpolation of consecutive pairs (x_i,\theta_i^\star)).
If m is differentiable, m is convex if m(\mathbf{x})+ \nabla m(\mathbf{x})^{\text{T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})
More generally, if m is convex, then there exists \xi_{\mathbf{x}}\in\mathbb{R}^n such that m(\mathbf{x})+ \xi_{\mathbf{x}}^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})
\xi_{\mathbf{x}} is a subgradient of m at {\mathbf{x}}. And then \partial m(\mathbf{x})=\big\lbrace m(\mathbf{x})+ \xi^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y}),\forall \mathbf{y}\in\mathbb{R}^n\big\rbrace
Hence, \mathbf{\theta}^\star is solution of \text{argmin}\big\lbrace\|\mathbf{y}-\mathbf{\theta}\|^2\big\rbrace\text{subject to }\theta_i+\xi_i^{\text{ T}}[\mathbf{x}_j-\mathbf{x}_i]\leq\mathbf{\theta}_j,~\forall i,j and \xi_1,\cdots,\xi_n\in\mathbb{R}^n. Now, to do it for real, use cobs package for constrained (b)splines regression,
library(cobs) |
To get a convex regression, use
plot(cars) x = cars$speed y = cars$dist rc = conreg(x,y,convex=TRUE) lines(rc, col = 2) |

Here we can get the values of the knots
rc Call: conreg(x = x, y = y, convex = TRUE) Convex regression: From 19 separated x-values, using 5 inner knots, 7, 8, 9, 20, 23. RSS = 1356; R^2 = 0.8766; needed (5,0) iterations |
and actually, if we use them in a linear-spline regression, we get the same output here
reg = lm(dist~bs(speed,degree=1,knots=c(4,7,8,9,,20,23,25)),data=cars) u = seq(4,25,by=.1) v = predict(reg,newdata=data.frame(speed=u)) lines(u,v,col="green") |
Let us add vertical lines for the knots
abline(v=c(4,7,8,9,20,23,25),col="grey",lty=2) |

Sixth post of our series on classification from scratch. The latest one was on the lasso regression, which was still based on a logistic regression model, assuming that the variable of interest Y has a Bernoulli distribution. From now on, we will discuss technique that did not originate from those probabilistic models, even if they might still have a probabilistic interpretation. Somehow. Today, we will start with neural nets.
Maybe I should start with a disclaimer. The goal is not to replicate well designed R functions, used for predictive modeling. It is simply to get a basic understanding of what’s going on.
First of all, neurals nets are nets, or networks. I will skip the parallel with “neural” stuff because it does not help me understanding what is happening (all apologies for my poor knowledge on biology, and cells)

So, it’s about some network. Networks have nodes, and edges (possibly connected) that connect nodes,

or maybe, to more specific (at least it helped me understanding what’s going on), some sort of flow network,

In such a network, we usually have sources (here multiple) sources (here \color{red}\{s_1,s_2,s_3\}), on the left, on a sink (here \{\color{blue}t\}), on the right. To continue with this metaphorical introduction, information from the sources should reach the sink. An usually, sources are explanatory variables, \{\mathbf{x}_1,\cdots,\mathbf{x}_p\}, and the sink is our variable of interest \mathbf{y}. And we want to create a graph, from the sources to the sink. We will have directed edges, with only one (unique) direction, where we will put weights. It is not a flow, the parallel with flow will stop here. For instance, the most simple network will be the following one, with no layer (i.e no node between the source and the sink)

The output here is a binary variable y\in\{0,1\} (it can also be y\in\{-1,+1\} but here, it’s not a big deal). In our network, our output will be y\in(0,1), because it is more easy to handly. For instance, consider y=f(something), for some function f taking values in (0,1). One can consider the sigmoid functionf(x)=\frac{1}{1+e^{-x}}=\frac{e^{x}}{e^{x}+1}which is actually the logistic function (so we should not be surprised to have results somehow close the logistic regression…). This function f is called the activation function, and there are thousands of such functions. If y\in\{-1,+1\}, people consider the hyperbolic tangentf(x)=\tanh(x)={\frac {(e^{x}-e^{-x})}{(e^{x}+e^{-x})}}or the inverse tangent function
f(x)=\tan ^{-1}(x)And as input for such function, we consider a weighted sum of incoming nodes. So herey_i=f\left(\sum_{j=1}^p\omega_j x_{j,i}\right)We can also add a constant actuallyy_i=f\left(\omega_0+\sum_{j=1}^p\omega_j x_{j,i}\right)So far, we are not far away from the logistic regression. Except that our starting point was a probabilistic model, in the sense that the later was interpreted as a probability (the probability that Y=1) and we wanted the model with the highest likelihood. But we’ll talk about selection of weights later on. First, let us construct our first (very simple) neural network. First, we have the sigmoid function
sigmoid = function(x) 1 / (1 + exp(-x)) |
The consider some weights. In our model with seven explanatory variables, with need 7 weights. Or 8 if we include the constant term. Let us consider \mathbf{\omega}=\mathbf{1},
weights_0 = rep(1,8) X = as.matrix(cbind(1,myocarde[,1:7])) y_5_1 = sigmoid(X %*% weights_0) |
that’s kind of stupid because all our predictions are 1, here. Let us try something else. Like \mathbf{\omega}=\widehat{\mathbf{\beta}}^{ols}. It is optimized, somehow, but we needed something to visualize what’s going on
weights_0 = lm(PRONO~.,data=myocarde)$coefficients |
then use
y_5_1 = sigmoid(X %*% weights_0) |
In order to see if we get a “good” prediction, let use plot the ROC curve, and compare it with the one we got with a (simple) logistic regression
library(ROCR) pred = ROCR::prediction(y_5_1,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) reg = glm(PRONO~.,data=myocarde,family=binomial(link = "logit")) y_0 = predict(reg,type="response") pred0 = ROCR::prediction(y_0,myocarde$PRONO) perf0 = ROCR::performance(pred0,"tpr", "fpr") plot(perf0,add=TRUE,col="red") |

That’s not bad for a very first attempt. Except that we’ve been cheating here, since we did use \mathbf{\omega}=\widehat{\mathbf{\beta}}^{ols}. How, for real, should we choose those weights?
Well, if we want an “optimal” set of weights, we need to “optimize” an objective function. So we need to quantify the loss of a mistake, between the prediction, and the observation. Consider here a quadratic loss function
loss = function(weights){ mean( (myocarde$PRONO-sigmoid(X %*% weights))^2) } |
It might be stupid to use a quadratic loss function for a classification, but here, it’s not the point. We just want to understand what is the algorithm we use, and the loss function \ell is just one parameter. Then we want to solve\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n\ell\left(y_i,f(\omega_0+\mathbf{x}_i^T\mathbf{\omega})\right)\right\rbraceThus, consider
weights_1 = optim(weights_0,loss)$par |
(where the starting point is the OLS estimate). Again, to see what’s going on, let us visualize the ROC curve
y_5_2 = sigmoid(X %*% weights_1) pred = ROCR::prediction(y_5_2,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf0,add=TRUE,col="red") |

That’s not amazing, but again, that’s only a first step.
Let us add a single layer in our network.

Those nodes are connected to the sources (incoming from sources) from the left, and then connected to the sink, on the right. Those nodes are not inter-connected. And again, for that network, we need edges (i.e series of weights). For instance, on the network above, we did add one single layer, with (only) three nodes.
For such a network, the prediction formula is \mathbf{y}=f\left( \omega_0+ \sum_{h=1}^3\omega_h f_h\left(\omega_{h,0}+ \sum_{j=1}^p \omega_{h,j} x_j\right)\right)or more synthetically\mathbf{y}=f\left( \omega_0+ \sum_{h=1}^3 \omega_hf_h\left(\omega_{h,0}+ \mathbf{x}^T\mathbf{\omega}_h\right)\right)Usually, we consider the same activation function everywhere. Don’t ask me why, I find that weird.
Now, we have a lot of weights to choose. Let us use again OLS estimates
weights_1 <- lm(PRONO~1+FRCAR+INCAR+INSYS+PAPUL+PVENT,data=myocarde)$coefficients X1 = as.matrix(cbind(1,myocarde[,c("FRCAR","INCAR","INSYS","PAPUL","PVENT")])) weights_2 <- lm(PRONO~1+INSYS+PRDIA,data=myocarde)$coefficients X2=as.matrix(cbind(1,myocarde[,c("INSYS","PRDIA")])) weights_3 <- lm(PRONO~1+PAPUL+PVENT+REPUL,data=myocarde)$coefficients X3=as.matrix(cbind(1,myocarde[,c("PAPUL","PVENT","REPUL")])) |
In that case, we did specify edges, and which sources (explanatory variables) should be used for each additional node. Actually, here, other techniques could be have been used, like using a PCA. Each node will then be one of the components. But we’ll use that idea later on…
X = cbind(sigmoid(X1 %*% weights_1), sigmoid(X2 %*% weights_2), sigmoid(X3 %*% weights_3)) |
But we’re not done here. Those were weights from the source to the know nodes, in the layer. We still need the weights from the nodes to the sink. Here, let use use a simple average
weights = c(1/3,1/3,1/3) y_5_3 <- sigmoid(X %*% weights) |
Again, we can plot the ROC curve to see what we’ve done…
pred = ROCR::prediction(y_5_3,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf0,add=TRUE,col="red") |

Now, we need some optimal selection of those weights. Observe that with only 3 nodes, there are already (7+1)\times3+3=27 parameters in that model! Clearly, parcimony is not the major issue when you start using neural nets! If p(\mathbf{x})=f\left( \omega_0+ \sum_{h=1}^3 \omega_hf_h\left(\omega_{h,0}+ \mathbf{x}^T\mathbf{\omega}_h\right)\right)we want to solve\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n\ell\left(y_i,p(\mathbf{x}_i)\right)\right\rbracefor some loss function, which is\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n (y_i-p(\mathbf{x}_i))^2 \right\rbracefor the quadratic norm, or\mathbf{\omega}^\star=\text{argmin}\left\lbrace\frac{1}{n}\sum_{i=1}^n (y_i\log p(\mathbf{x}_i)+[1-y_i]\log [1-p(\mathbf{x}_i)]) \right\rbraceif we want to use cross-entropy.
For convenience, let us center all the variable we create, otherwise, we get numerical problems.
center = function(z) (z-mean(z))/sd(z) loss = function(weights){ weights_1 = weights[0+(1:7)] weights_2 = weights[7+(1:7)] weights_3 = weights[14+(1:7)] weights_ = weights[21+1:4] X1=X2=X3=as.matrix(myocarde[,1:7]) Z1 = center(X1 %*% weights_1) Z2 = center(X2 %*% weights_2) Z3 = center(X3 %*% weights_3) X = cbind(1,sigmoid(Z1), sigmoid(Z2), sigmoid(Z3)) mean( (myocarde$PRONO-sigmoid(X %*% weights_))^2)} |
Now that we have our objective function, consider some starting points. We can consider weights from a PCA, and then use a gradient descent algorithm,
pca = princomp(myocarde[,1:7]) W = get_pca_var(pca)$contrib weights_0 = c(W[,1],W[,2],W[,3],c(-1,rep(1,3)/3)) weights_opt = optim(weights_0,loss)$par |
The prediction is then obtained using
weights_1 = weights_opt[0+(1:7)] weights_2 = weights_opt[7+(1:7)] weights_3 = weights_opt[14+(1:7)] weights_ = weights_opt[21+1:4] X1=X2=X3=as.matrix(myocarde[,1:7]) Z1 = center(X1 %*% weights_1) Z2 = center(X2 %*% weights_2) Z3 = center(X3 %*% weights_3) X = cbind(1,sigmoid(Z1), sigmoid(Z2), sigmoid(Z3)) y_5_4 = sigmoid(X %*% weights_) |
And as previously, why not plot the ROC curve of that model
pred = ROCR::prediction(y_5_4,myocarde$PRONO) perf = ROCR::performance(pred,"tpr", "fpr") plot(perf,col="blue",lwd=2) plot(perf,add=TRUE,col="red") |

That’s not too bad. But with 27 coefficients, that’s what we would expect, no?
That’s more or less what is done in neural nets functions. Let us now have a look at some dedicated R functions.
library(nnet) myocarde_minmax = myocarde minmax = function(z) (z-min(z))/(max(z)-min(z)) for(j in 1:7) myocarde_minmax[,j] = minmax(myocarde_minmax[,j]) |
Here, variables are linearly transformed, to take values in (0,1). Then we can construct a neural network with one single layer, and three nodes,
model_nnet = nnet(PRONO~.,data=myocarde_minmax,size=3) summary(model_nnet) a 7-3-1 network with 28 weights options were - b->h1 i1->h1 i2->h1 i3->h1 i4->h1 i5->h1 i6->h1 i7->h1 -9.60 -1.79 21.00 14.72 -20.45 -5.05 14.37 -17.37 b->h2 i1->h2 i2->h2 i3->h2 i4->h2 i5->h2 i6->h2 i7->h2 4.72 2.83 -3.37 -1.64 1.49 2.12 2.31 4.00 b->h3 i1->h3 i2->h3 i3->h3 i4->h3 i5->h3 i6->h3 i7->h3 -0.58 -6.03 25.14 18.03 -1.19 7.52 -19.47 -12.95 b->o h1->o h2->o h3->o -1.32 29.00 -10.32 26.27 |
Here, it is the complete full network. And actually, there are (online) some functions that can he used to visualize that network
library(devtools) source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r') plot.nnet(model_nnet) |

Nice, isn’t it? We clearly see the intermediary layer, with three nodes, and on top the constants. Edges are the plain lines, the darker, the heavier (in terms of weights).
Other R functions can actually be considered.
library(neuralnet) model_nnet = neuralnet(formula(glm(PRONO~.,data=myocarde_minmax)), myocarde_minmax,hidden=3, act.fct = sigmoid) plot(model_nnet) |

Again, for the same network structure, with one (hidden) layer, and three nodes in it.
The good thing is that it’s not possible to add more layers. Like two layers. Nodes from the first layer are no longuer connected with the sink, but with nodes in the second layer. And those nodes will then be connected to the sink. We now have something like
p(\mathbf{x})=f\left( \omega_0+ \sum_{h=1}^3 \omega_h f_h\left(\omega_{h,0}+ \mathbf{z}_h^T\mathbf{\omega}_h\right)\right)where\mathbf{z}_h=f\left( \omega_{h,0}+ \sum_{j=1}^{k_h} \omega_{h,j} f_{h,j}\left(\omega_{h,j,0}+ \mathbf{x}^T\mathbf{\omega}_{h,j}\right)\right)I may be rambling here (a little bit) but that’s a lot of parameters. Here is the visualization of such a network,
library(neuralnet) model_nnet = neuralnet(formula(glm(PRONO~.,data=myocarde_minmax)), myocarde_minmax,hidden=3, act.fct = sigmoid) plot(model_nnet) |

Let us get back on our simple dataset, with only two covariates.
library(neuralnet) df_minmax =df df_minmax$y=(df_minmax$y=="1")*1 minmax = function(z) (z-min(z))/(max(z)-min(z)) for(j in 1:2) df_minmax[,j] = minmax(df[,j]) X = as.matrix(cbind(1,df_minmax[,1:2])) |
Consider only one layer, with two nodes
model_nnet = neuralnet(formula(lm(y~.,data=df_minmax)), df_minmax,hidden=c(2)) plot(model_nnet) |

Here, we did not specify it, but the activation function is the sigmoid (actually, it is called logistic here)
model_nnet$act.fct function (x) { 1/(1 + exp(-x)) } attr(,"type") [1] "logistic" f=model_nnet$act.fct |
The weights (on the figure) can be obtained using
w0 = model_nnet$weights[[1]][[2]][,1] w1 = model_nnet$weights[[1]][[1]][,1] w2 = model_nnet$weights[[1]][[1]][,2] |
Now, to get our prediction,
we should usep(\mathbf{x})=f\left( \omega_0+ \omega_1 f(\omega_{1,0}+ \mathbf{x}_h^T\mathbf{\omega}_{1,1:2})+\omega_1 f(\omega_{2,0}+ \mathbf{x}_h^T\mathbf{\omega}_{2,1:2})\right)which can be obtained using
f(cbind(1,f(X%*%w1),f(X%*%w2))%*%w0) [,1] [1,] 0.7336477343 [2,] 0.7317999050 [3,] 0.7185803540 [4,] 0.7404005280 [5,] 0.7518482779 [6,] 0.4939774149 [7,] 0.4965876378 [8,] 0.7101714888 [9,] 0.5050760026 [10,] 0.5049877644 |
Unfortunately, it is not the output of the model here,
neuralnet::prediction(model_nnet) Data Error: 0; $rep1 x1 x2 y 1 0.1250 0.0000000000 0.02030470787 2 0.0625 0.1176470588 0.89621706711 3 0.9375 0.2352941176 0.01995171956 4 0.0000 0.4705882353 1.10849420363 5 0.5000 0.4705882353 -0.01364966058 6 0.3125 0.5294117647 -0.02409150561 7 0.6875 0.8235294118 0.93743057765 8 0.3750 0.8823529412 1.01320924782 9 1.0000 0.9058823529 1.04805134309 10 0.5625 1.0000000000 1.00377379767 |
If anyone has a clue, I’d be glad to know what went wrong here… I find that odd to have outputs outside the (0,1) interval, but the output is neitherp(\mathbf{x})=\omega_{0,0}+ \omega_{0,1} f(\omega_{1,0}+ \mathbf{x}_h^T\mathbf{\omega}_{1,1:2})+\omega_{0,2} f(\omega_{2,0}+ \mathbf{x}_h^T\mathbf{\omega}_{2,1:2})
cbind(1,f(X%*%w1),f(X%*%w2))%*%w0 [,1] [1,] 1.01320924782 [2,] 1.00377379767 [3,] 0.93743057765 [4,] 1.04805134309 [5,] 1.10849420363 [6,] -0.02409150561 [7,] -0.01364966058 [8,] 0.89621706711 [9,] 0.02030470787 [10,] 0.01995171956 |
Fourth post of our series on classification from scratch, following the previous post which was some sort of detour on kernels. But today, we’ll get back on the logistic model.
We’ve seen before that the classical estimation technique used to estimate the parameters of a parametric model was to use the maximum likelihood approach. More specifically, \widehat{\mathbf{\beta}}=\text{argmax}\lbrace \log\mathcal{L}(\mathbf{\beta}|\mathbf{x},\mathbf{y})\rbraceThe objective function here focuses (only) on the goodness of fit. But usually, in econometrics, we believe something like non sunt multiplicanda entia sine necessitate (“entities are not to be multiplied without necessity”), the parsimony principle, simpler theories are preferable to more complex ones. So we want to penalize for too complex models.
This is not a bad idea. It is mentioned here and there in econometrics textbooks, but usually, for model choice, not about the inference. Usually, we estimate parameters using maximum likelihood techniques, and them we use AIC or BIC to compare two models. Recall that Akaike (AIC) criteria is based on-2\log\mathcal{L}(\widehat{\mathbf{\beta}}|\mathbf{x},\mathbf{y})+2\text{dim}(\widehat{\mathbf{\beta}})We have on the left a measure for the goodness of fit, and on the right, a penalty increasing with the “complexity” of the model.
Very quickly, here, the complexity is the number of variates used. I will not enter into details about the concept of sparsity (and the true dimension of the problem), I will recommend to read the book by Martin Wainwright, Robert Tibshirani and Trevor Hastie on that issue. But assume that we do not make and variable selection, we consider the regression on all covariates. Define\Vert\mathbf{a} \Vert_{\ell_0}=\sum_{i=1}^d \mathbf{1}(a_i\neq 0), ~~\Vert\mathbf{a} \Vert_{\ell_1}=\sum_{i=1}^d |a_i|,~~\Vert\mathbf{a} \Vert_{\ell_2}=\left(\sum_{i=1}^d a_i^2\right)^{1/2}for any \mathbf{a}\in\mathbb{R}^d. One might say that the AIC could be written-2\log\mathcal{L}(\widehat{\mathbf{\beta}}|\mathbf{x},\mathbf{y})+2\|\widehat{\mathbf{\beta}}\|_{\ell_0}And actually, this will be our objective function. More specifically, we will consider
\widehat{\mathbf{\beta}}_{\lambda}=\text{argmin}\lbrace -\log\mathcal{L}(\mathbf{\beta}|\mathbf{x},\mathbf{y})+\lambda\|\mathbf{\beta}\|\rbracefor some norm \|\cdot\|. I will not get back here on the motivation and the (theoretical) properties of those estimates (that will actually be discussed in the Summer School in Barcelona, in July), but in this post, I want to discuss the numerical algorithm to solve such optimization problem, for \|\cdot\|_{\ell_2} (the Ridge regression) and for \|\cdot\|_{\ell_1} (the LASSO regression).
The problem of \|\mathbf{\beta}\| is that the norm should make sense, somehow. A small \mathbf{\beta}_j is with respect to the “dimension” of x_j‘s. So, 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) |
Before running some codes, recall that we want to solve something like\widehat{\mathbf{\beta}}_{\lambda}=\text{argmin}\lbrace -\log\mathcal{L}(\mathbf{\beta}|\mathbf{x},\mathbf{y})+\lambda\|\mathbf{\beta}\|_{\ell_2}^2\rbrace In the case where we consider the log-likelihood of some Gaussian variable, we get the sum of the square of the residuals, and we can obtain an explicit solution. But not in the context of a logistic regression.
The heuristics about Ridge regression is the following graph. In the background, we can visualize the (two-dimensional) log-likelihood of the logistic regression, and the blue circle is the constraint we have, if we rewite the optimization problem as a contrained optimization problem : \min_{\mathbf{\beta}:\|\mathbf{\beta}\|^2_{\ell_2}\leq s} \lbrace \sum_{i=1}^n -\log\mathcal{L}(y_i,\beta_0+\mathbf{x}^T\mathbf{\beta}) \rbracecan be written equivalently (it is a strictly convex problem)\min_{\mathbf{\beta},\lambda} \lbrace -\sum_{i=1}^n \log\mathcal{L}(y_i,\beta_0+\mathbf{x}^T\mathbf{\beta}) +\lambda \|\mathbf{\beta}\|_{\ell_2}^2 \rbraceThus, the constrained maximum should lie in the blue disk
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) u = seq(-1,1,length=251) lines(u,sqrt(1-u^2),type="l",lwd=2,col="blue") lines(u,-sqrt(1-u^2),type="l",lwd=2,col="blue") |

Let us consider the objective function, with the following code
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(beta^2) } |
Why not try a standard optimisation routine ? In the very first post on that series, we did mention that using optimization routines were not clever, since they were strongly relying on the starting point. But here, it is not the case
lambda = 1 beta_init = lm(PRONO~.,data=myocarde)$coefficients vpar = matrix(NA,1000,8) for(i in 1:1000){ vpar[i,] = optim(par = beta_init*rnorm(8,1,2), function(x) PennegLogLik(x,lambda), method = "BFGS", control = list(abstol=1e-9))$par} par(mfrow=c(1,2)) plot(density(vpar[,2]),ylab="",xlab=names(myocarde)[1]) plot(density(vpar[,3]),ylab="",xlab=names(myocarde)[2]) |

Clearly, even if we change the starting point, it looks like we converge towards the same value. That could be considered as the optimum.
The code to compute \widehat{\mathbf{\beta}}_{\lambda} would then be
opt_ridge = function(lambda){ beta_init = lm(PRONO~.,data=myocarde)$coefficients logistic_opt = optim(par = beta_init*0, function(x) PennegLogLik(x,lambda), method = "BFGS", control=list(abstol=1e-9)) logistic_opt$par[-1]} |
and we can visualize the evolution of \widehat{\mathbf{\beta}}_{\lambda} as a function of {\lambda}
v_lambda = c(exp(seq(-2,5,length=61))) est_ridge = Vectorize(opt_ridge)(v_lambda) library("RColorBrewer") colrs = brewer.pal(7,"Set1") plot(v_lambda,est_ridge[1,],col=colrs[1]) for(i in 2:7) lines(v_lambda,est_ridge[i,],col=colrs[i]) |

At least it seems to make sense: we can observe the shrinkage as \lambda increases (we’ll get back to that later on).
We’ve seen that we can also use Newton Raphson to solve this problem. Without the penalty term, the algorithm was\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}}where
\frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}=\mathbf{X}^T(\mathbf{y}-\mathbf{p}_{old})and\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=-\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}where \mathbf{\Delta}_{old} is the diagonal matrix with terms \mathbf{p}_{old}(1-\mathbf{p}_{old}) on the diagonal.
Thus\mathbf{\beta}_{new} = \mathbf{\beta}_{old} + (\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T[\mathbf{y}-\mathbf{p}_{old}]that we can also write\mathbf{\beta}_{new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}where \mathbf{z}=\mathbf{X}\mathbf{\beta}_{old}+\mathbf{\Delta}_{old}^{-1}[\mathbf{y}-\mathbf{p}_{old}]. Here, on the penalized problem, we can easily prove that\frac{\partial\log\mathcal{L}_p(\mathbf{\beta}_{\lambda,old})}{\partial\mathbf{\beta}}=\frac{\partial\log\mathcal{L}(\mathbf{\beta}_{\lambda,old})}{\partial\mathbf{\beta}}-2\lambda\mathbf{\beta}_{old}while\frac{\partial^2\log\mathcal{L}_p(\mathbf{\beta}_{\lambda,old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{\lambda,old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}-2\lambda\mathbb{I}Hence\mathbf{\beta}_{\lambda,new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}+2\lambda\mathbb{I})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}
The code is then
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) X = cbind(1,X) 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])) Delta = matrix(0,nrow(X),nrow(X));diag(Delta)=(pi*(1-pi)) z = X%*%beta[,s] + solve(Delta)%*%(Y-pi) B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%Delta%*%z) beta = cbind(beta,B)} beta[,8:10] [,1] [,2] [,3] XInter 0.59619654 0.59619654 0.59619654 XFRCAR 0.09217848 0.09217848 0.09217848 XINCAR 0.77165707 0.77165707 0.77165707 XINSYS 0.69678521 0.69678521 0.69678521 XPRDIA -0.29575642 -0.29575642 -0.29575642 XPAPUL -0.23921101 -0.23921101 -0.23921101 XPVENT -0.33120792 -0.33120792 -0.33120792 XREPUL -0.84308972 -0.84308972 -0.84308972 |
Again, it seems that convergence is very fast.
And interestingly, with that algorithm, we can also derive the variance of the estimator\text{Var}[\widehat{\mathbf{\beta}}_{\lambda}]=[\mathbf{X}^T\mathbf{\Delta}\mathbf{X}+2\lambda\mathbb{I}]^{-1}\mathbf{X}^T\mathbf{\Delta}\text{Var}[\mathbf{z}]\mathbf{\Delta}\mathbf{X}[\mathbf{X}^T\mathbf{\Delta}\mathbf{X}+2\lambda\mathbb{I}]^{-1}where\text{Var}[\mathbf{z}]=\mathbf{\Delta}^{-1}
The code to compute \widehat{\mathbf{\beta}}_{\lambda} as a function of \lambda is then
newton_ridge = function(lambda=1){ beta = as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8) for(s in 1:20){ pi = exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) Delta = matrix(0,nrow(X),nrow(X));diag(Delta)=(pi*(1-pi)) z = X%*%beta[,s] + solve(Delta)%*%(Y-pi) B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%Delta%*%z) beta = cbind(beta,B)} Varz = solve(Delta) Varb = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% t(X)%*% Delta %*% Varz %*% Delta %*% X %*% solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) return(list(beta=beta[,ncol(beta)],sd=sqrt(diag(Varb))))} |
We can visualize the evolution of \widehat{\mathbf{\beta}}_{\lambda} (as a function of \lambda)
v_lambda=c(exp(seq(-2,5,length=61))) est_ridge=Vectorize(function(x) newton_ridge(x)$beta)(v_lambda) library("RColorBrewer") colrs=brewer.pal(7,"Set1") plot(v_lambda,est_ridge[1,],col=colrs[1],type="l") for(i in 2:7) lines(v_lambda,est_ridge[i,],col=colrs[i]) |

and to get the evolution of the variance
v_lambda=c(exp(seq(-2,5,length=61))) est_ridge=Vectorize(function(x) newton_ridge(x)$sd)(v_lambda) library("RColorBrewer") colrs=brewer.pal(7,"Set1") plot(v_lambda,est_ridge[1,],col=colrs[1],type="l") for(i in 2:7) lines(v_lambda,est_ridge[i,],col=colrs[i],lwd=2) |

Recall that when \lambda=0 (on the left of the graphs), \widehat{\mathbf{\beta}}_{0}=\widehat{\mathbf{\beta}}^{mco} (no penalty). Thus as \lambda increase (i) the bias increase (estimates tend to 0) (ii) the variances deacrease.
As always, there are R functions availble to run a ridge regression. Let us use the glmnet function, with \alpha=0
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) library(glmnet) glm_ridge = glmnet(X, y, alpha=0) plot(glm_ridge,xvar="lambda",col=colrs,lwd=2) |

as a function of the norm

the \ell_1 norm here, I don’t know why. I don’t know either why all graphs obtained with different optimisation routines are so different… Maybe that will be for another post…
An interesting case is obtained when covariates are orthogonal. This can be obtained using a PCA of the covariates.
library(factoextra) pca = princomp(X) pca_X = get_pca_ind(pca)$coord |
Let us run a ridge regression on those (orthogonal) covariates
library(glmnet) glm_ridge = glmnet(pca_X, y, alpha=0) plot(glm_ridge,xvar="lambda",col=colrs,lwd=2) |

plot(glm_ridge,col=colrs,lwd=2) |

We clearly observe the shrinkage of the parameters, in the sense that \widehat{\mathbf{\beta}}_{\lambda}^{\perp}=\frac{\widehat{\mathbf{\beta}}^{mco}}{1+\lambda}
Let us try with our second set of data
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] = (df0[,j]-m[j])/s[j] reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=0,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=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) } |
We can try various values of \lambda
reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=0) par(mfrow=c(1,2)) plot(reg,xvar="lambda",col=c("blue","red"),lwd=2) abline(v=log(.2)) plot_lambda(.2) |

or
reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=0) par(mfrow=c(1,2)) plot(reg,xvar="lambda",col=c("blue","red"),lwd=2) abline(v=log(1.2)) plot_lambda(1.2) |

Next step is to change the norm of the penality, with the \ell_1 norm (to be continued…)
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.
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 asm(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.
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…
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.
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<=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)<=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…