Fitting a Boosted-Tree Model With gbm
The gbm package provides an R interface for fitting gradient-boosted tree models. A useful fit is not just a call to a function: the response distribution, tree structure, learning rate, sampling settings, and number of trees all need to match the modeling task. This tutorial focuses on specifying that fit and obtaining predictions at candidate iteration counts. The earlier material on gradient boosting explains how successive trees update the model; here, the goal is to translate those choices into a careful R workflow.
The examples use hypothetical data. A fitted model should be trained using only the data designated for training. Any validation or test data used to assess predictions should remain separate from that fit.
Arguments That Define the Fit
A common call has the form gbm(formula, distribution, data, n.trees, interaction.depth, shrinkage, bag.fraction, n.minobsinnode, train.fraction). The formula names the response and predictors, and data supplies the data frame containing them. The remaining arguments control the model fit. They should be stated explicitly in an analysis so that the modeling choices can be read and checked.
| Argument | Role in the fit | Check before fitting |
|---|---|---|
| distribution | Chooses the response loss, such as gaussian or bernoulli | Does the choice match the response type and coding? |
| n.trees | Sets the maximum number of boosting iterations | Is it large enough to include all candidate counts to be examined? |
| interaction.depth | Limits the depth of each tree and hence the complexity of interactions it can represent | Does the chosen depth match the intended tree complexity? |
| shrinkage | Scales each tree’s contribution to the current model | Is the learning rate recorded alongside the iteration count? |
| bag.fraction | Sets the fraction of training observations sampled for each tree | Is the intended stochastic or full-row fitting behavior clear? |
| n.minobsinnode | Sets the minimum number of observations in a terminal node | Can the tree make the desired splits without creating overly small leaves? |
| train.fraction | Uses a fraction of the supplied data for the package’s training portion | Is this consistent with the intended data partition? |
For a deliberately specified training data frame, train.fraction = 1 uses all rows passed to gbm for fitting. This does not mean that a separate validation set should be passed into the fit: keep that set outside the data argument, then supply it as newdata when predicting. If using a fraction below 1 or package cross-validation, be clear about what observations each operation uses. These controls do not replace a sound data-partitioning plan.
The basic syntax is illustrated by gbm(y ~ x1 + x2, distribution = "gaussian", data = train_data, n.trees = 300, interaction.depth = 2, shrinkage = 0.03, bag.fraction = 0.7, n.minobsinnode = 5, train.fraction = 1). For binary classification, change the response to a 0/1 variable and choose distribution = "bernoulli". The syntax alone cannot determine appropriate tuning values: the choices depend on the modeling task and should be evaluated using a validation design that does not reuse the final test data for tuning.
A Justification for Comparing Iteration Counts
The prediction formula follows directly from the stagewise construction. At each iteration, a new tree is fitted to improve the current model, and shrinkage scales the size of that update. After the first iteration, the model is \(F_1(x)=F_0(x)+\eta h_1(x)\). If after \(m-1\) iterations the model is \(F_{m-1}(x)=F_0(x)+\eta\sum_{j=1}^{m-1}h_j(x)\), adding the next scaled tree gives
This is why gbm can predict from a fit using, for example, 20 trees or 80 trees, provided the fitted object contains at least 80. Each request uses a different prefix of the same fitted sequence. It does not refit the model using a new learning rate or a different tree depth.
Worked Example: A One-Tree Gaussian Fit
A hypothetical response has four training observations. Fit a Gaussian model with one tree, depth one, shrinkage 0.5, full row sampling, and a minimum terminal-node size of one. The observations are \((x,y)=(1,2),(2,2),(3,8),(4,8)\). Work out the first tree’s predictions and check the resulting training mean squared error.
Method. A depth-one tree can make one split. For this simple squared-error fit, the starting prediction is the response mean, and the tree’s terminal-node adjustments are the mean residuals within the leaves. Full row sampling avoids changing the four-row data for this hand calculation.
The response mean is \((2+2+8+8)/4=20/4=5\). The initial residuals, in order of increasing \(x\), are \(2-5=-3\), \(2-5=-3\), \(8-5=3\), and \(8-5=3\). A split between \(x=2\) and \(x=3\) gives residuals \((-3,-3)\) in the lower-\(x\) leaf and \((3,3)\) in the upper-\(x\) leaf. Their mean residuals are \((-3-3)/2=-3\) and \((3+3)/2=3\), respectively.
Shrinkage scales each leaf adjustment by 0.5. Thus the lower-leaf prediction is \(5+0.5(-3)=3.5\), and the upper-leaf prediction is \(5+0.5(3)=6.5\). The R call is gbm(y ~ x, distribution = "gaussian", data = train_data, n.trees = 1, interaction.depth = 1, shrinkage = 0.5, bag.fraction = 1, n.minobsinnode = 1, train.fraction = 1).
The fitted values for the four training rows are \(3.5,3.5,6.5,6.5\). Their squared errors are \((2-3.5)^2=2.25\), \(2.25\), \((8-6.5)^2=2.25\), and \(2.25\). Therefore, the mean squared error is \((2.25+2.25+2.25+2.25)/4=9/4=2.25\). As a check, the initial constant prediction of 5 has squared error 9 on each row and mean squared error 9, so this fitted tree reduces the training error. The calculation avoids the mistake of adding the full leaf residual rather than its shrinkage-scaled value.
Fitting a Bernoulli Model and Requesting Predictions
For binary classification, code the outcome as 0 or 1 and use distribution = "bernoulli". The number of trees in the fit must be at least as large as the largest candidate iteration count you intend to request. To compare counts from one fitted model, pass each count as n.trees in a prediction call. When probabilities are required, request type = "response"; this makes the intended output scale explicit.
Worked Example: Fitting a Binary Classifier for Candidate Counts
A hypothetical screening data set has eight training observations and a binary response y. Four responses are 1 and four are 0. A single numeric predictor x is available. The analyst wants probability predictions at 20, 50, and 100 iterations, using depth-two trees, shrinkage 0.03, row sampling fraction 0.7, and minimum terminal-node size 2.
Method. Fit one Bernoulli model with a maximum of 100 trees, then request predictions at each candidate count. The observed prevalence is \(4/8=0.5\). As a baseline check, the intercept-only probability is 0.5; its log odds are \(\log(0.5/(1-0.5))=\log(1)=0\). This checks the response coding and gives a simple reference point, but it does not determine the fitted model’s later predictions.
The fit is gbm(y ~ x, distribution = "bernoulli", data = train_data, n.trees = 100, interaction.depth = 2, shrinkage = 0.03, bag.fraction = 0.7, n.minobsinnode = 2, train.fraction = 1). For a separate validation data frame, the three calls are predict(fit, newdata = validation_data, n.trees = 20, type = "response"), predict(fit, newdata = validation_data, n.trees = 50, type = "response"), and predict(fit, newdata = validation_data, n.trees = 100, type = "response").
Each call returns one predicted probability per validation row. The 50-tree prediction uses the first 50 trees, not an independently fitted 50-tree model. If the intended candidate is 120 trees, this object cannot provide that prediction because it was fitted with only 100; the model must first be fitted with a maximum of at least 120 trees. The common mistake this example avoids is fitting only 20 trees and then attempting to request 50 or 100.
Comparing Candidate Predictions on Validation Data
Obtaining predictions is not the same as selecting an iteration count. The metric must suit the task, and candidate counts should be compared on data not used to fit the model. For a binary response, one possible metric is mean log loss. Given observed binary outcomes \(y_i\) and predicted probabilities \(p_i\), the mean log loss is
Lower values indicate better probabilistic predictions on the evaluated observations. A probability of exactly 0 or 1 for an outcome that did not occur makes the log loss infinite, so predictions and metric calculations should be checked carefully. This metric is only one possible criterion; the choice should reflect the modeling objective.
Worked Example: Checking Two Candidate Counts With Log Loss
For four hypothetical validation cases, the observed outcomes are \(0,1,1,0\). Suppose predictions from the fitted model are, respectively, \((0.30,0.60,0.55,0.45)\) at 10 trees and \((0.20,0.70,0.60,0.40)\) at 30 trees. Compute mean log loss for both counts and choose the lower-loss candidate for this validation set.
Method. Apply the binary log-loss formula separately to each case, using \(-\log(p_i)\) when \(y_i=1\) and \(-\log(1-p_i)\) when \(y_i=0\). At 10 trees, the four losses are \(-\log(0.70)=0.3567\), \(-\log(0.60)=0.5108\), \(-\log(0.55)=0.5978\), and \(-\log(0.55)=0.5978\). Their sum is \(2.0631\), so the mean is \(2.0631/4=0.5158\), to four decimal places.
At 30 trees, the four losses are \(-\log(0.80)=0.2231\), \(-\log(0.70)=0.3567\), \(-\log(0.60)=0.5108\), and \(-\log(0.60)=0.5108\). Their sum is \(1.6014\), so the mean is \(1.6014/4=0.4004\), to four decimal places. The 30-tree candidate has lower validation log loss for these four cases. The calculation avoids choosing a count from training performance alone; it does not establish that the same count is best on every sample.
Practical Checks and Common Traps
Before interpreting a gbm fit, check that the response matches the selected distribution, the model was trained on the intended data, and each prediction request uses a count no greater than the fitted maximum. Also check that newdata contains the predictor variables required by the formula in compatible forms. Record the chosen settings with the predictions so that results from different fits are not accidentally compared as if they differed only in iteration count.
- Wrong distribution: A continuous response and a 0/1 outcome are not interchangeable. Use the loss appropriate to the target.
- Too few fitted trees: n.trees in predict cannot request iterations that were not fitted.
- Unclear prediction scale: When the desired output is a probability, specify type = "response" in a Bernoulli prediction call.
- Data leakage: Do not include validation observations in the fitting data and then describe their predictions as held out.
- Unfair candidate comparison: Predictions from different fits may differ in more than iteration count. State which settings were held fixed and how the candidates were evaluated.
Check Your Understanding
For each problem, identify the relevant fit or prediction choices and check the final answer against the stated data.
- A continuous response is to be modeled with squared-error loss. The analyst wants predictions at 25 and 60 trees, using depth-one trees and shrinkage 0.05. What distribution should be specified, and what is the minimum n.trees for the fit? (Answer: gaussian; at least 60 trees.)
- A Bernoulli model is fitted with n.trees = 80, but a prediction call requests n.trees = 100. Is that request supported by the fitted object? What change is needed? (Answer: No; fit at least 100 trees before requesting 100-tree predictions.)
- A binary training sample has 30 positive outcomes among 120 observations. Calculate its prevalence and the intercept-only log odds. (Answer: Prevalence \(=30/120=0.25\); log odds \(=\log(0.25/0.75)=\log(1/3)\approx-1.0986\).)
- For two validation outcomes \(y=(0,1)\), the predicted probabilities are \(p=(0.25,0.80)\). Calculate mean log loss to four decimal places. (Answer: \([-\log(0.75)-\log(0.80)]/2=(0.2877+0.2231)/2=0.2554\).)
- A fit uses train.fraction = 1 on a training data frame, and a separate validation data frame is available. Where should the validation rows be supplied when requesting predictions? (Answer: Keep them outside the fitting data and supply them as newdata in predict, with the desired n.trees and prediction type.)