When you start discussing with people in machine learning, you quickly hear something like “forget your econometric models, your GLMs, I can easily find a machine learning ‘model’ that can beat yours”. I am usually very sceptical, especially when I hear “easily” or “always“. I have no problem about the fact that I use old econometric models, but I had the feeling that things aren’t that easy. I can understand that we might have problems when we do have a lot of features (I am still working on that, I’ll get back to this point soon), but I have the feeling that I can still capture interactions, and non-linearities with standard econometric models as well as any machine learning algorithm.

Just to illustrate, consider the following ‘model‘

where is (just to illustrate)

> n <- 5000
> rtf <- function(x1, x2) { sin(x1+x2)/(x1+x2) }
> xgrid <- seq(1,6,length=31)
> ygrid <- seq(1,6,length=31)
> zgrid <- outer(xgrid,ygrid,rtf)
> persp(xgrid,ygrid,zgrid,theta=30, phi=30, 
+ col="green", ticktype="detailed",shade=TRUE)

Here

> df <- data.frame(x1=(runif(n, min=1, max=6)),
+                  x2=(runif(n, min=1, max=6)))
> df$m=rtf(df$x1, df$x2)
> df$y=df$m+rnorm(n,sd=.07)

To visualise those observations on the graph above, use

> m=persp(xgrid,ygrid,zgrid*0+100,theta=30,
+ phi=30, col="white")
> pt=trans3d(df$x1,df$x2,df$y,m)
> points(pt$x,pt$y,cex=.7)

or, if we play with the color (blue on top, when is large, and red on the bottom, when is small)

How could we estimate , based on that sample ? What alternative to regression models can we get ? (that is more or less what I plan to discuss today).

I am an econometrician (that’s what I usually claim). What about a linear regression?

> reg=lm(y~x1+x2,data=df)
> p=function(x1,x2) predict(reg,newdata=data.frame(x1=x1,x2=x2))
> xgrid=seq(1,6,length=31)
> ygrid=seq(1,6,length=31)
> zgrid=outer(xgrid,ygrid,p)
> persp(xgrid,ygrid,zgrid,theta=30, phi=30,
+ col="green",shade=TRUE)

It is not great, let us be honest. But we shouldn’t be suprised. It is a nonlinear model, and we try to approximate that nonlinear function by some linear one. We should expect that nonlinearity cannot be captured be a linear regression. Especially when the model is not monotonous. Why not an additive model?

> library(splines)
> reg=lm(y~bs(x1)+bs(x2),data=df)
> p=function(x1,x2) predict(reg,newdata=data.frame(x1=x1,x2=x2))

It is slightly better. I said ‘slightly’. We need something else, here. Why not try a multivariate adaptive regression splines model (so called MARS)? It is claimed that those models “can be seen as an extension of linear models that automatically models non-linearities and interactions between variables“. That what we’re looking for, right?

> library(earth)
> reg=earth(y~x1+x2,data=df)
> p=function(x1,x2) predict(reg,newdata=data.frame(x1=x1,x2=x2))

We’re still flat, here. Why not use some bivariate splines, directly?

> library(mgcv)
> reg=gam(y~s(x1,x2),data=df)
> p=function(x1,x2) predict(reg,newdata=data.frame(x1=x1,x2=x2))

Here we’re good! Aren’t we? So with econometric models, it is possible to capture nonlinearities and interractions.

Now, what about alternatives? What are those machine learning algorithm everyone is talking about? A popular one is the regression tree,

> library(rpart)
> reg=rpart(y~x1+x2,data=df,method="anova")
> p=function(x1,x2) predict(reg,newdata=data.frame(x1=x1,x2=x2),type="vector")
> zgrid=outer(xgrid,ygrid,p)
> persp(xgrid,ygrid,zgrid,theta=30, phi=30, 
+ col="green",shade=TRUE)

It is not that bad. Of course, it is simplistic, because with a tree, predictions are constant on rectangular areas. But still.

Using bagging, or random forests, we should get something nice

> library(randomForest)
> reg=randomForest(y~x1+x2,data=df)
> p=function(x1,x2) as.numeric(predict(reg,newdata=              data.frame(x1=x1,x2=x2),type="response"))
> zgrid=outer(xgrid,ygrid,p)
> persp(xgrid,ygrid,zgrid,theta=30, phi=30,
+ col="green",shade=TRUE)

Not bad. What else can be consider? -nearest neighbours can be also an interesting alternative.

> k=10
> p=function(x1,x2){
+   d=(df$x1-x1)^2+(df$x2-x2)^2
+   return(mean(df$y[which(rank(d)<=k)]))
+ }
> zgrid=outer(xgrid,ygrid,Vectorize(p))
> persp(xgrid,ygrid,zgrid,theta=30, phi=30, 
+ col="green",shade=TRUE)

Here it is working well, too.  Last, but not least, what about boosting techniques? To be honest, it still sounds like magic to me. I can get the heuristics about it (the fact that a set of weak learners can create a strong learner, somehow). I think I also understand the theoretical foundations of gradient boosting techniques – but to be honest, it took me some time, since a gradient descent is usually based on Taylor’s expansion

but here, we have something like

where the part on the right is interpreted as a gradient (but in a much larger space, since it is a functional space). Weak learners – here decision stumps – are used to approximate that complicated function. Why not…

Anyway, I think I understand what this boosting idea is, but so far, I have not been able to write my own function (where I loop, and where I grow a tree on the residuals, and then create new residuals, and another tree, etc). I am still working on it, but if anyone is willing to share some R codes (DIY boosting regression tree), I am interested.

If we use some R functions, we can get

> library(dismo)
> reg=gbm.step(data=df, gbm.x = 1:2, gbm.y = 4,
+ family = "gaussian", tree.complexity = 5,
+ learning.rate = 0.01, bag.fraction = 0.5)
 
 
 GBM STEP - version 2.9 
 
Performing cross-validation optimisation of a boosted regression tree model 
for y and using a family of gaussian 
Using 5000 observations and 2 predictors 
creating 10 initial models of 50 trees 
 
 folds are unstratified 
total mean deviance =  0.0182 
tolerance is fixed at  0 
ntrees resid. dev. 
50    0.0136 
now adding trees... 
100   0.0109 
150   0.0093 
200   0.0082 
250   0.0075 
300   0.0069 
350   0.0065

(etc)

1800   0.005 
1850   0.005 
1900   0.005 
1950   0.005 
fitting final gbm model with a fixed number of  1450  trees for  y 
 
mean total deviance = 0.018 
mean residual deviance = 0.003 
 
estimated cv deviance = 0.005 ; se = 0 
 
training data correlation = 0.924 
cv correlation =  0.852 ; se = 0.013 
 
elapsed time -  0.09 minutes

To visualize the output, use

> p=function(x1,x2) predict(reg,newdata=data.frame(x1=x1,x2=x2),
+ n.trees=1400)
> zgrid=outer(xgrid,ygrid,p)
> persp(xgrid,ygrid,zgrid,theta=30, phi=30, 
+ col="green",shade=TRUE)

Nice isn’t it? So clearly, there are alternatives to standard regression models. But econometric models aren’t that bad, when we take into account nonlinearities. And so far, I still feel more confortable with my good old econometric models (but I am still learning about those alternatives).

---

---