<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/boost-series-i-advantage-in-smooth-functions/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://insightr.wordpress.com/2018/08/20/boost-series-i-advantage-in-smooth-functions/"> R – insightR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<div><strong>By Gabriel Vasconcelos and Yuri Fonseca</strong></div>
</div>
<h2>Introduction</h2>
<p>This is the first of a series of post on the BooST (Boosting Smooth Trees). If you missed the first post introducing the model click <a rel="nofollow" href="https://insightr.wordpress.com/2018/08/14/boost-boosting-smooth-trees-a-new-machine-learning-model-for-partial-effect-estimation-in-nonlinear-regressions/" target="_blank">here</a> and if you want to see the full article click <a rel="nofollow" href="https://arxiv.org/abs/1808.03698" target="_blank">here</a>. The BooST is a model that uses Smooth Trees as base learners, which makes it possible to approximate the derivative of the underlying model. In this post, we will show some examples on generated data of how the BooST approximates the derivatives and we also will discuss how the BooST may be a good choice when dealing with smooth functions if compared to the usual discrete Regression Trees.</p>
<p><span></span></p>
<h2>Installation</h2>
<p>To install the BooST package in R you should run the following code:</p>
<pre>library(devtools)
install_github(&quot;gabrielrvsc/BooST&quot;)
</pre>
<p> </p>
<p>Note that this R implementation is suited just for small problems. If you want to use the BooST on larges instances with more speed we recommend the Julia implementation until the C++ is ready. The Julia package can be installed with:</p>
<p><code>Pkg.clone(&quot;<a rel="nofollow" href="https://github.com/gabrielrvsc/BooSTjl.jl%22" target="_blank">https://github.com/gabrielrvsc/BooSTjl.jl&quot;</a>)</code></p>
<p>and then <code>using BooSTjl</code> in your Julia terminal. Both packages have documentation for all exported functions.</p>
<h2>Example 1: Cosine</h2>
<p>The first example is the one we briefly discussed in the previous post. We are going to generate data from:</p>
<p><img alt="displaystyle y_i = cos( pi [x_{i,1}+x_{i,2}])+ varepsilon_i" title="displaystyle y_i = cos( pi [x_{i,1}+x_{i,2}])+ varepsilon_i" src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+y_i+%3D+%5Ccos%28+%5Cpi+%5Bx_%7Bi%2C1%7D%2Bx_%7Bi%2C2%7D%5D%29%2B+%5Cvarepsilon_i&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>where <img alt="x_{i,1} sim N(0,1)" title="x_{i,1} sim N(0,1)" src="https://s0.wp.com/latex.php?latex=x_%7Bi%2C1%7D+%5Csim+N%280%2C1%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, <img alt="x_{i,2}" title="x_{i,2}" src="https://s0.wp.com/latex.php?latex=x_%7Bi%2C2%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is a Bernoulli with <img alt="p = 0.5" title="p = 0.5" src="https://s0.wp.com/latex.php?latex=p+%3D+0.5&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="varepsilon_i sim N(0,1)" title="varepsilon_i sim N(0,1)" src="https://s0.wp.com/latex.php?latex=%5Cvarepsilon_i+%5Csim+N%280%2C1%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. The function used to generate the data is:</p>
<pre>dgp = function(N,r2){
  X = matrix(rnorm(N*2,0,1),N,2)
  X[,ncol(X)] = base::sample(c(0,1),N,replace=TRUE)
  aux = X
  yaux = cos(pi*(rowSums(X)))
  vyaux = var(yaux)
  ve = vyaux*(1-r2)/r2
  e = rnorm(N,0,sqrt(ve))
  y = yaux+e
  var(yaux)/var(y)
  return(list(y = y, X = X))
}
</pre>
<p>We are going to generate 1000 observations with an R2 of 0.3 (a lot of noise). The code below generates the data and runs the BooST and the Boosting using the xgboost package. We estimated 300 trees in each model with a step of 0.2. The last lines in the code just organize the results in a data.frame and generate values for the real function we want to recover.</p>
<pre>library(BooST)
library(tidyverse)
library(xgboost)
library(reshape2)

set.seed(1)
data = dgp(N = 1000, r2 = 0.3)
y = data$y
x = data$X

set.seed(1)
BooST_Model = BooST(x,y, v = 0.2, M = 300 ,display = TRUE)
xgboost_Model = xgboost(x,label = y, nrounds = 300, params = list(eta = 0.2, max_depth = 3))

x1r = rep(seq(-4,4,length.out = 1000), 2)
x2r = c(rep(0,1000), rep(1,1000))
yr = cos(pi*(x1r+x2r))
real_function = data.frame(x1 = x1r, x2 = as.factor(x2r), y = yr)
fitted = data.frame(x1 = x[,1],x2 = as.factor(x[,2]), BooST = fitted(BooST_Model),
                    xgboost = predict(xgboost_Model,x), y = y)
</pre>
<p>Before going into the results let’s have a look at the data in the figure below. The two black lines are the real cosine function we used and the dots are the data we generated. In this first look it seems hard to recover the real function from this data.</p>
<pre>ggplot() + geom_point(data = fitted, aes(x = x1, y = y, color = x2)) + geom_line(data = real_function, aes(x = x1, y = y, linetype = x2))
</pre>
<p><img alt="plot of chunk unnamed-chunk-7" src="https://i2.wp.com/i.imgur.com/kzF3lYz.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The next figure shows the result with regular Boosting using the xgboost package. The points in blue and red are what we fitted and the points in gray are the data from the previous plot. All the structure in the model comes from the cosine function represented by the two black lines. The main conclusion here is that we over-fitted the data.</p>
<pre>ggplot() + geom_point(data = fitted, aes(x = x1, y = y), color = &quot;gray&quot;) + geom_point(data = fitted, aes(x = x1, y = xgboost, color = x2)) + geom_line(data = real_function, aes(x = x1, y = y, linetype = x2))
</pre>
<p><img alt="plot of chunk unnamed-chunk-8" src="https://i2.wp.com/i.imgur.com/aQxZt5X.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The next plot shows what we obtained with the BooST. Again, red and blue points are fitted values and the data is in gray. The model fits the function very well with a few exceptions on extreme points where we have much less data.</p>
<pre>ggplot() + geom_point(data = fitted, aes(x = x1, y = y), color = &quot;gray&quot;) + geom_point(data = fitted, aes(x = x1, y = BooST, color = x2)) + geom_line(data = real_function, aes(x = x1, y = y, linetype = x2))
</pre>
<p><img alt="plot of chunk unnamed-chunk-9" src="https://i1.wp.com/i.imgur.com/IWX0vXX.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Next, let’s have a look at the derivatives. The code below estimates them and organizes the results for the plot.</p>
<pre>BooST_derivative = estimate_derivatives(BooST_Model, x, 1)
derivative = data.frame(x1 = x[,1],x2 = as.factor(x[,2]), derivative = BooST_derivative)
dr = -1*sin(pi*(x1r+x2r))*pi
real_function$derivative = dr
</pre>
<p>The results are in the next plot. As we can see, the model also estimates the derivatives very well. However, the performance deteriorates as we go to the boarders where we have less data. This is a natural feature of many nonparametric models, which are more precise where the data are more dense.</p>
<pre>ggplot() + geom_point(data = fitted, aes(x = x1, y = BooST_derivative, color = x2)) + geom_line(data = real_function, aes(x = x1, y = derivative, linetype = x2))
</pre>
<p><img alt="plot of chunk unnamed-chunk-11" src="https://i1.wp.com/i.imgur.com/rAj9rTJ.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Finally, let’s generate new data from the same dgp and see how the BooST and the Boosting perform. The output of the following code is the BooST RMSE divided by the Boosting RMSE.</p>
<pre>set.seed(2)
dataout=dgp(N=1000,r2=0.3)
yout = dataout$y
xout = dataout$X

p_BooST = predict(BooST_Model, xout)
p_xgboost = predict(xgboost_Model, xout)

sqrt(mean((p_BooST - yout)^2))/sqrt(mean((p_xgboost - yout)^2))
</pre>
<pre>## [1] 0.9165293
</pre>
<h2>Example 2: More variables interacting</h2>
<p>In the previous example we only had two variables. What happens if we add more? In the next example we are going to generate data from a dgp with 10 variables, where the function we are interested in comes from all second order interactions with all variables:</p>
<p><img alt="displaystyle y_i = f(x_i) + varepsilon_i," title="displaystyle y_i = f(x_i) + varepsilon_i," src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+y_i+%3D+f%28x_i%29+%2B+%5Cvarepsilon_i%2C&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"><br>
where,</p>
<p><img alt="displaystyle f(x_i) = sum_{j=1}^{10}sum_{k=i}^{10}x_{ij}x_{ik}," title="displaystyle f(x_i) = sum_{j=1}^{10}sum_{k=i}^{10}x_{ij}x_{ik}," src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+f%28x_i%29+%3D+%5Csum_%7Bj%3D1%7D%5E%7B10%7D%5Csum_%7Bk%3Di%7D%5E%7B10%7Dx_%7Bij%7Dx_%7Bik%7D%2C&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>with all variables and <img alt="varepsilon_i" title="varepsilon_i" src="https://s0.wp.com/latex.php?latex=%5Cvarepsilon_i&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> generated from a standard Normal distribution. The difficulty comes from the fact that we have to many interactions. Nevertheless, let’s see how the BooST recovers derivatives in this setup. Since we now have 10 variables interacting a plot like the one we did with the cosine is no longer possible. A good way to make this example more visual is to keep all variables fixed except one and see the derivative as we move on this one variable. However, this strategy needs a lot of data to work because we may be calculating derivatives in parts of the space that were poorly mapped. First, the dgp function:</p>
<pre>dgp2 = function(N,k,r2){
  yaux = rep(0,N)
  x = matrix(rnorm(N*k),N,k)
  for(i in 1:k){
    for(j in i:k)
      yaux = yaux  +  x[,i]*x[,j]
  }
  vyaux=var(yaux)
  ve=vyaux*(1-r2)/r2
  e=rnorm(N,0,sqrt(ve))
  y=yaux+e
  var(yaux)/var(y)
  return(list(y=y,x=x))
}
</pre>
<p>Now let’s estimate the model. Given the complexity of this model, we adopted a more conservative strategy with a step of 0.1 and smaller values for gamma (controls the transition on the logistic function). These adjustment may require more trees to converge. Therefore, we used M=1000 trees.</p>
<pre>set.seed(1)
data = dgp2(1000,10,0.7)
x = data$x
y = data$y

set.seed(1)
BooSTmv = BooST(x,y, v = 0.1, M = 1000, display = TRUE, gamma = seq(0.5,1.5,0.01))
</pre>
<p>The next step is to put the data in the way we need to calculate the derivative. We are going to look at the derivative of <img alt="y_i" title="y_i" src="https://s0.wp.com/latex.php?latex=y_i&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> with respect to the first variable <img alt="x_{i1}" title="x_{i1}" src="https://s0.wp.com/latex.php?latex=x_%7Bi1%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> keeping all other variables in the mean, which is 0. The solution in this case will be <img alt="2x_{i1}" title="2x_{i1}" src="https://s0.wp.com/latex.php?latex=2x_%7Bi1%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> for the derivative.</p>
<pre>xrep = x
xrep[,2:ncol(x)] = 0
derivative_mv = estimate_derivatives(BooSTmv,xrep,1)

df = data.frame(x1 = x[,1], derivative = derivative_mv)
xr1 = seq(-4,4,0.01)
dfr = data.frame(x1 = xr1, derivative = 2*xr1)
</pre>
<p>Finally, the results. The black line is the real derivative and the blue dots are what we estimated. Given the complexity of the problem the results are very good. The blue dots are always close to the black line except by some extreme values where we have less data.</p>
<pre>ggplot() + geom_point(data = df, aes(x = x1, y = derivative),color = &quot;blue&quot;) + geom_line(data = dfr, aes(x = x1, y = derivative)) + xlim(-4,4)
</pre>
<p><img alt="plot of chunk unnamed-chunk-16" src="https://i1.wp.com/i.imgur.com/k2oEs1p.png?w=456&ssl=1" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://insightr.wordpress.com/2018/08/20/boost-series-i-advantage-in-smooth-functions/"> R – insightR</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/jsRo62pXUh8" height="1"></a></div>