(This article was originally published at Statistics – Freakonometrics, and syndicated at StatsBlogs.)
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).
Please comment on the article here: Statistics – Freakonometrics
