Generalized Boosted Models: A guide to the gbm package

1 Generalized Boosted Models: a guide to the gbm packageGreg RidgewayAugust 3, 2007 Boosting takes on various forms with different programs using different lossfunctions, different base models, and different optimization schemes. The gbmpackage takes the approach described in [2] and [3]. Some of the terminologydiffers, mostly due to an effort to cast boosting terms into more standard sta-tistical terminology ( deviance). In addition, the gbm package implementsboosting for models commonly used in statistics but not commonly associatedwith boosting. The Cox proportional hazard model , for example, is an incred-ibly useful model and the boosting framework applies quite readily with onlyslight modification [5]. Also some algorithms implemented in the gbm packagediffer from the standard implementation. The AdaBoost algorithm [1] has aparticular loss function and a particular optimization algorithm associated withit.

2 The gbm implementation of AdaBoost adopts AdaBoost s exponential lossfunction (its bound on misclassification rate) but uses Friedman s gradient de-scent algorithm rather than the original one proposed. So the main purposes ofthis document is to spell out in detail what the gbm package Gradient boostingThis section essentially presents the derivation of boosting described in [2]. Thegbm package also adopts the stochastic gradient boosting strategy, a small butimportant tweak on the basic algorithm, described in [3]. Friedman s gradient boosting machineFriedman (2001) and the companion paper Friedman (2002) extended the workof Friedman, Hastie, and Tibshirani (2000) and laid the ground work for a newgeneration of boosting algorithms. Using the connection between boosting andoptimization, this new work proposes the Gradient Boosting any function estimation problem we wish to find a regression function, f(x), that minimizes the expectation of some loss function, (y, f), as shownin (4).

3 F(x) = arg minf(x)Ey,x (y, f(x))1 Initialize f(x) to be a constant, f(x) = arg min Ni=1 (yi, ).Fortin 1, .. , Tdo1. Compute the negative gradient as the working responsezi= f(xi) (yi, f(xi))|f(xi)= f(xi)(1)2. Fit a regression model ,g(x), predictingzifrom the Choose a gradient descent step size as = arg min N i=1 (yi, f(xi) + g(xi))(2)4. Update the estimate off(x) as f(x) f(x) + g(x)(3)Figure 1: Friedman s Gradient Boost algorithm2= arg minf(x)Ex[Ey|x (y, f(x)) x](4)We will focus on finding estimates off(x) such that f(x) = arg minf(x)Ey|x[ (y, f(x))|x](5)Parametric regression models assume thatf(x) is a function with a finite numberof parameters, , and estimates them by selecting those values that minimize aloss function ( squared error loss) over a training sample ofNobservationson (y,x) pairs as in (6). = arg min N i=1 (yi, f(xi; ))(6)When we wish to estimatef(x) non-parametrically the task becomes more dif-ficult.

4 Again we can proceed similarly to [4] and modify our current estimateoff(x) by adding a new functionf(x) in a greedy fashion. Lettingfi=f(xi),we see that we want to decrease theNdimensional functionJ(f) =N i=1 (yi, f(xi))=N i=1 (yi, Fi).(7)The negative gradient ofJ(f) indicates the direction of the locally greatestdecrease inJ(f). Gradient descent would then have us modifyfas f f J(f)(8)where is the size of the step along the direction of greatest descent. Clearly,this step alone is far from our desired goal. First, it only fitsfat values ofxfor which we have observations. Second, it does not take into account thatobservations with similarxare likely to have similar values off(x). Both theseproblems would have disastrous effects on generalization error. However, Fried-man suggests selecting a class of functions that use the covariate informationto approximate the gradient, usually a regression tree. This line of reasoningproduces his Gradient Boosting algorithm shown in Figure 1.

5 At each itera-tion the algorithm determines the direction, the gradient, in which it needs toimprove the fit to the data and selects a particular model from the allowableclass of functions that is in most agreement with the direction. In the case ofsquared-error loss, (yi, f(xi)) = Ni=1(yi f(xi))2, this algorithm correspondsexactly to residual are various ways to extend and improve upon the basic frameworksuggested in Figure 1. For example, Friedman (2001) substituted several choices3in for to develop new boosting algorithms for robust regression with leastabsolute deviation and Huber loss functions. Friedman (2002) showed thata simple subsampling trick can greatly improve predictive performance whilesimultaneously reduce computation time. Section 2 discusses some of Improving boosting methods using control ofthe learning rate, sub-sampling, and a decom-position for interpretationThis section explores the variations of the previous algorithms that have thepotential to improve their predictive performance and interpretability.

6 In par-ticular, by controlling the optimization speed or learning rate, introducing low-variance regression methods, and applying ideas from robust regression we canproduce non-parametric regression procedures with many desirable a by-product some of these modifications lead directly into implementationsfor learning from massive datasets. All these methods take advantage of thegeneral form of boosting f(x) f(x) + E(z(y, f(x))|x).(9)So far we have taken advantage of this form only by substituting in our favoriteregression procedure for Ew(z|x). I will discuss some modifications to estimatingEw(z|x) that have the potential to improve our Decreasing the learning rateAs several authors have phrased slightly differently, ..boosting, whatever flavor,seldom seems to overfit, no matter how many terms are included in the additiveexpansion . This is not true as the discussion to [4] points the update step of any boosting algorithm we can introduce a learningrate to dampen the proposed move.

7 F(x) f(x) + E(z(y, f(x))|x).(10)By multiplying the gradient step by as in equation 10 we have control on therate at which the boosting algorithm descends the error surface (or ascends thelikelihood surface). When = 1 we return to performing full gradient (2001) relates the learning rate to regularization through optimal number of iterations,T, and the learning rate, , depend oneach other. In practice I set to be as small as possible and then selectTbycross-validation. Performance is best when is as small as possible performancewith decreasing marginal utility for smaller and smaller . Slower learning ratesdo not necessarily scale the number of optimal iterations. That is, if when = and the optimalTis 100 iterations, doesnotnecessarily imply thatwhen = the optimalTis 1000 Variance reduction using subsamplingFriedman (2002) proposed the stochastic gradient boosting algorithm that sim-ply samples uniformly without replacement from the dataset before estimatingthe next gradient step.

8 He found that this additional step greatly improved per-formance. We estimate the regression E(z(y, f(x))|x) using a random subsampleof the ANOVA decompositionCertain function approximation methods are decomposable in terms of a func-tional ANOVA decomposition . That is a function is decomposable asf(x) = jfj(xj) + jkfjk(xj, xk) + jk`fjk`(xj, xk, x`) + .(11)This applies to Boosted trees. Regression stumps (one split decision trees) de-pend on only one variable and fall into the first term of 11. Trees with two splitsfall into the second term of 11 and so on. By restricting the depth of the treesproduced on each boosting iteration we can control the order of additive components are sufficient to approximate a multivariate functionwell, Generalized additive models, the na ve Bayes classifier, and Boosted stumpsfor example. When the approximation is restricted to a first order we can alsoproduce plots ofxjversusfj(xj) to demonstrate how changes inxjmight affectchanges in the response Relative influenceFriedman (2001) also develops an extension of a variable s relative influence forboosted estimates.

9 For tree based methods the approximate relative influenceof a variablexjis J2j= splits onxjI2t(12)whereI2tis the empirical improvement by splitting onxjat that point. Fried-man s extension to Boosted models is to average the relative influence of variablexjacross all the trees generated by the boosting Common user optionsThis section discusses the options to gbm that most users will need to changeor Loss functionThe first and foremost choice isdistribution. This should be easily dictatedby the application. For most classification problems eitherbernoulliorad-aboostwill be appropriate, the former being recommended. For continuous5 Select a loss function (distribution) the number of iterations,T( ) the depth of each tree,K( ) the shrinkage (or learning rate) parameter, (shrinkage) the subsampling rate,p( )Initialize f(x) to be a constant, f(x) = arg min Ni=1 (yi, )Fortin 1, .. , Tdo1. Compute the negative gradient as the working responsezi= f(xi) (yi, f(xi))|f(xi)= f(xi)(13)2.

10 Randomly selectp Ncases from the dataset3. Fit a regression tree withKterminal nodes,g(x) = E(z|x). This tree isfit using only those randomly selected observations4. Compute the optimal terminal node predictions, 1, .. , K, as k= arg min xi Sk (yi, f(xi) + )(14)whereSkis the set ofxs that define terminal Update f(x) as f(x) f(x) + k(x)(15)wherek(x) indicates the index of the terminal node into which an obser-vation with featuresxwould fall. Again this step uses only the randomlyselected observationsFigure 2: Boosting as implemented ingbm()6outcomes the choices aregaussian(for minimizing squared error),laplace(forminimizing absolute error), and quantile regression (for estimating percentilesof the conditional distribution of the outcome). Censored survival outcomesshould requirecoxph. Count outcomes may usepoissonalthough one mightalso considergaussianorlaplacedepending on the analytical The relationship between shrinkage and number of it-erationsThe issues that most new users of gbm struggle with are the choice It is important to know that smaller values ofshrinkage(al-most) always give improved predictive performance.