Partial Dependency Plots and GBM

My favorite "go-to-first" modeling algorithm for a classification or regression task is Gradient Boosted Regression Trees (Friedman, 2001 and 2002), especially as coded in the GBM package of R. Gradient boosted regression is one of the focuses of my masters thesis (along with random forests) and has gotten more and more attention as more and more data mining contests have been won, at least in part, by employing this method. Some of the reasons I pick Gradient Boosted Regression Trees as the best off the shelf predictive modeling algorithm available are: 

• High predictive accuracy across many domains and target data types
• Ability to specify various loss functions (Gaussian, Bernoulli,  Poisson etc.) as well as run survival analysis via Cox proportional hazards, quantile regression etc.   
• Handles mixed data types (continuous and nominal)
• Seamlessly deals with missing values
• Contains out-of-bag (OOB) estimates of error
• Contains variable importance measures
• Contains variable interaction measures / detection
• Allows estimates of marginal effects of a predictor (s) via Partial Dependency Plots.

This latter point is a nice feature coded into the GBM package that gives the analyst the ability to produce univariate and bivariate partial dependency plots. These plots enable the researcher to understand the effect of a predictor variable (or interaction between two predictors) on the target outcome, given the other predictors (partialling them out - after accounting for their average effects). 

The technique is not unique to Gradient Boosted Regression Trees and in fact is a general method for understanding any black box modeling algorithm. When we use ensembles for instance, this tends to be the only way to understand how changing the value of a predictor, say a binary predictor  from 0 to 1 or a continuous predictor from within it's observed range, effects the outcome, given the model (i.e. accounting for the impact of other predictors in the model).

The idea of these “marginal” plots is to displays the modeled outcome, over the range of a given predictor. Hastie (2009) describes this general technique as considering the full model function $f$, depending on a small subset of predictors we are interested in, $X_{s}$, where this subset is typically one or two predictors in practice, and the other predictors in the model , $X_{c}$. This full model function is then written as $f(X)=f(X_{s},X_{c})$ where $X=X_{s} \cup X_{c}$ . A partial dependence plot displays the marginal expected value of $f(X_{s})$, by averaging over the values of $X_{c}$. Practically, this would be given by $f(X_{s}=x_{s})=Avg(f(X_{s}=x_{s}, X_{ci}))$ where the expected value of the function for a given value (vector) of the subset of predictors is the average value of the model output setting the subset predictors to the value of interest and averaging the result with the values of $X_{c}$ as they exist in the data set. 

A brute force method would be to create a new data set from the training data, but repeating it for every pattern in $X_{s}$ of interest plugged in and then using the model to output the average value for each distinct value of  $X_{s}$.

Need for Repeated Hold Outs in Predictive Models

Many models are built and deployed using a training/validation partitioning approach. Under this construct a data set is randomly split in two parts - one part that is used to train the model and the other part which is held out to validate how well the model works. The split is sometimes done 50/50 and other times more data is left for training than validating- especially when one is dealing with a rare event as the target variable (which almost always is the case in database marketing). Personally, I use 65/35 train/validate as the standard. Sometimes, a third partition might be sampled out to further tune the model before the final validation on a "test" set. The only time I have used this approach is with calibrating scores into probabilities, but it is used frequently by others and is the default in SAS Enterprise Miner when partitioning.

There are many types of validation, from this split sample approach to bootstrapping, k-fold cross validation, Leave One Out Cross Validation (LOOCV) and many others. In R, the Caret package is a fantastic work bench for tuning and validating predictive models.

In database marketing, it seems that many (majority?) of models are validated with a single hold-out partition.

Here is quick example on the need to understand the variability that is inherent with the random sampling split between train and validation. If it is done one time, the results may be quite misleading.

Below displays the error curves for a model built on a data set with 152,000 records. 65% was sampled to train the model and 35% was held back to validate. In the full data set there were about 3% positive events (this is a binary target). The model is a form of Generalized Naive Bayes.

The process followed to create this plot was simple:

• Split the data set into train and validate with a random seed.
• Do all processing to the data (e.g. feature selection) that involves the target variable or could be considered a tuning parameter. [It is vital to do this independently each time in any cross validation process].
• Train the model.
• Predict the outcome on the validation set.
• Rank the validation data into 10 equally sized groups (deciles) by the predicted score.
• For each decile, calculate the average actual "response" (i.e. proportion of 1's), the average predicted response and the difference between these.

Repeat the above multiple times - here it was done 10 times. To make this feasible, the pre-processing, feature selection and modeling needs to be something that is rather automated (programmatically). This is not always possible.

There is a large degree of variability between deciles for each run, that a decision maker likely needs to understand and that the modeler needs to be aware of - that could well be lost with a single hold out partition.

I would be interested to know what others do in these cases. I lean towards presenting something like below- from the same data as above - where a smooth is fit (loess) to the average actual prediction by decile, along with this same error data.

Delta Method in Logistic Regression

In the last post we looked at how to construct and evaluate a simple linear hypothesis test within the frame work of logistic regression (or any generalized linear model fit via maximum likelihood). While the example was simple, the idea can be quite powerful for quickly testing hypothesis concerning regression coefficients.

For this method to work however you need to be able to write your (null) hypothesis in terms of a linear combination of the model parameters or coefficients. In our simple example, this was achieved as we wished to test the null hypothesis $\beta_{2}-\beta_{3}=0$. What happens when we can not construct the weights ($L$) to place on the coefficients with a vector (or a matrix when interested in multiple simultaneous contrasts such as  Ho:  $\beta_{2}=\beta_{3}=0$ )? We need to apply other methods to the estimation and inference task. An often used classical procedure is called the Delta Method.

Delta Method

The main idea of the delta method is in those cases where you are not facing a simple linear sum of observations (random variables), create a linear approximation to that more complicated function using a Taylor expansion and derive a variance of that approximation instead. This method will allow you to derive an approximate variance in large samples and thus construct a confidence interval for theses non-linear estimators.

There are many sources for information on the theory of the delta method, including:

• Wikipedia
• Xu and Long
• Greene Page 70
• Cooch and White

but the main result of interest here is that in large samples under "standard regularity conditions" a                                                      
function $G(b)$ of the coefficients of our generalized linear model ($b_{0},b_{1}....b_{p}$) will be distributed as $N(G(\beta),\frac{\partial G(\beta)}{\partial \beta'}VAR(\hat{\beta})\frac{\partial G(\beta)}{\partial \beta}$ ) which we evaluate at the ML estimates:

Point Estimate: $G(\hat{\beta})$
Variance: $\frac{\partial G(\hat{\beta})}{\partial \hat{\beta}'}VAR(\hat{\beta})\frac{\partial G(\hat{\beta})}{\partial \hat{\beta}}$

For example, returning to the last post, lets say we are interested in the ratio quantity $\frac{\beta_{2}}{\beta_{3}}$.

We have the vector of maximum likelihood estimates :

$[-2.8211865  -0.5794947  -0.4589188  -0.3883266 ]$

and the estimated variance-covariance matrix:

                       (Intercept)                 Gifts1                Gifts2                 Gifts3

(Intercept)  0.0002674919  -0.0002674919 -0.0002674919  -0.0002674919
Gifts1          -0.0002674919  0.0034667138   0.0002674919   0.0002674919
Gifts2          -0.0002674919  0.0002674919   0.0039732911   0.0002674919
Gifts3          -0.0002674919  0.0002674919   0.0002674919   0.0039964383

From here, we know the point estimate of this quantity is just the expression evaluated at the maximum likelihood estimates:

$\frac{ -0.4589188}{ -0.3883266}=1.181778$

Next, we need to get the first derivative of the function:

$\frac{\partial G(\beta)}{\partial \beta}$ which is $[0,0,\frac{1}{\beta_{3}}, \frac{-\beta_{2}}{(\beta_{3})^{2}}]^{t}$ and evaluated at the ML estimates $[0,0, \frac{1}{-0.3883266}, \frac{-0.4589188}{-0.3883266^{2}}]^{t}$

Thus the variance of $\frac{\beta_{2}}{\beta_{3}}$ is thus 0.05916905 and a 95% confidence interval could be constructed as $1.181778\pm1.96\sqrt{0.0591690}=(0.7050221, 1.65855)$

The car package in R provides a very nice short-cut to all of this, supplying symbolic differentiation to save you from calculating the derivative. The deltaMethod function accepts a variety of different model objects or you can supply a mean vector and vcov matrix. Here is a replication of the simple analysis above:



Another common usage for this method is to compute confidence intervals for predicted values out of generalized linear models on the response scale - for example in logistic regression we may be interested in the effect of a predictor variable on the probability p and not on the link scale of the logit. Xu and Long show the derivation for several predicted probabilities in the link above.

• SAS NLMIXED uses the delta method for estimating the variance with its "estimate" statement - for example when you are interested in the variance of the difference in estimated probabilities between two different predictor values (or vector of values).  

 When you have time and computing power, a bootstrap analysis will typically provide better coverage. In a future post we will compare bootstrapping, simulation and the delta method.