Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
54,964 characters · 18 sections · 24 citation commands
pystacked: Stacking generalization and machine learning in Stata
\graphicspath{{fig/}}
\makeatletter \makeatother
\inserttype[st0001]{article} { Achim Ahrens\\ETH Z\"urich \\[email removed] \and Christian B. Hansen\\ University of Chicago\\[email removed] \and Mark E. Schaffer\\ Heriot-Watt University\\Edinburgh, United Kingdom \\ [email removed] }
When faced with a new prediction or classification task, it is a priori rarely obvious which machine learning algorithm is best suited. A common approach is to evaluate the performance of a set of machine learners on a hold-out partition of the data or via cross-validation and then select the machine learner that minimizes a chosen loss metric. However, this approach is incomplete as combining multiple learners into one final prediction might lead to superior performance compared to each individual learner. This possibility motivates stacked generalization, or simply stacking, due to Wolpert1992 and Breiman1996a. Stacking is a form of model averaging. Theoretical results in Laan2007 support the use of stacking as it performs asymptotically at least as well as the best-performing individual learner as long as the number of base learners is not too large.
In this article, we introduce \stcmd{pystacked}\xspace for stacking regression and binary classification in Stata. \stcmd{pystacked}\xspace allows users to fit multiple machine learning algorithms via Python's {\it scikit-learn}\xspace scikit-learn,sklearn_api and combine these into one final prediction as a weighted average of individual predictions. \stcmd{pystacked}\xspace adds to the growing number of programs for machine learning in Stata, including, among others, lassopack for regularized regression Ahrens2019, rforest for random forests Schonlau2020 and svm for support vector machines Schonlau2016,Guenther2018. Similarly to \stcmd{pystacked}\xspace, cerulli2022machine and Droste2020 provide an interface to {\it scikit-learn}\xspace in Stata. mlrtime allows Stata users to make use of R's parsnip machine learning library Huntington2021. \stcmd{pystacked}\xspace differs from these in that it is, to our knowledge, the first to make stacking available to Stata users. Furthermore, \stcmd{pystacked}\xspace can also be used to fit a single machine learner and thus provides an easy-to-use and versatile API to {\it scikit-learn}\xspace's machine learning algorithms.
Stacking is widely used in applied predictive modeling in many disciplines, e.g., for predicting mortality hwangbo2022, bankruptcy filings liang2020,fedorova2022, or temperatures hooker2018. The use of stacking as a method and \stcmd{pystacked}\xspace as a program is, however, not only restricted to pure prediction or classification tasks. A growing literature exploits machine learning to facilitate causal inference Athey2019c. Indeed, a motivation for writing \stcmd{pystacked}\xspace is that it can be used in combination with \stcmd{ddml}\xspace ddml_package,ddml_papers, which implements the Double-Debiased Machine Learning (DDML) methodology of Chernozhukov2018. DDML utilizes cross-fitting, a form of iterative sample splitting, which allows leveraging a wide class of supervised machine learners, including stacking, for the estimation of causal parameters. For instance, in the context of DDML, stacking can be used to estimate the conditional expectation of an outcome with respect to confounders or propensity scores.
We stress that \stcmd{pystacked}\xspace relies on Python's {\it scikit-learn}\xspace (version 0.24 or higher) and the ongoing work of the {\it scikit-learn}\xspace contributors. Thus, \stcmd{pystacked}\xspace relies on Stata's Python integration which was introduced in Stata 16.0. We kindly ask users to cite {\it scikit-learn}\xspace along with this article when using \stcmd{pystacked}\xspace. Throughout we refer to version 0.7 of \stcmd{pystacked}\xspace, which is the latest version at the time of writing. Please check for updates to \stcmd{pystacked}\xspace on a regular basis and consult the help file to be informed about new features. The \stcmd{pystacked}\xspace help file includes information on how to install a recent Python version and set up Stata's Python integration.
Section (ref) introduces stacking. Section (ref) presents the main features of the \stcmd{pystacked}\xspace program. Section (ref) demonstrates the use of the program using examples.
In this section, we briefly summarize the stacking approach for regression and binary classification tasks. van2011targeted provides a book-length treatment; for a concise textbook treatment, see Hastie2009.
We first focus on stacking for regression problems where the aim is to predict the continuous outcome $y_i$ using predictors $\bm{x}_i$. The idea of stacking is to combine a set of “base” (or “level-0”) learners using a “final” (or “level-1”) estimator. It is advisable to include a relatively large and diverse set of base learners to capture different types of patterns in the data. The same algorithm can also be included more than once using different tuning or hyper-tuning parameters. Typical choices for base learners are regularized regression or ensemble methods, such as random forests or gradient boosting.
In the first step of stacking, we obtain cross-validated predicted values $\hat{y}^{-k(i)}_{(j),i}$ for each base learner $j$ and observation $i$. The super-script “$-k(i)$” indicates that we form the cross-fitted predicted value for observation $i$ by fitting the learner to all folds except fold $k(i)$, which is the fold that includes observation $i$. The use of cross-validation is necessary as stacking would otherwise give more weight to base learners that suffer from over-fitting. The second step is to fit a final learner using the observed $y_i$ as the outcome and the cross-validated predicted values $\hat{y}^{-k(i)}_{(1),i},\ldots,\hat{y}^{-k(i)}_{(J),i}$ as predictors. A typical choice for the final learner is constrained least squares, which enforces the stacking weights to be non-negative and sum to one. This restriction facilitates the interpretation of stacking as a weighted average of base learners and may lead to better performance Breiman1996a,Hastie2009. Algorithm 1 summarizes the stacking algorithm for regression problems as it is implemented in \stcmd{pystacked}\xspace.
It is instructive to compare Step 2 with the classical `winner-takes-all' approach that selects one base learner as the one which exhibits the lowest cross-validated loss. Stacking, in contrast, may assign non-zero weights to multiple base learners, thus combining their strengths to produce a better overall predictor than any of the individual base learners. Other choices for the final learner are possible. In addition to the default final learner, \stcmd{pystacked}\xspace supports, among others, non-negative least squares without the constraint $\sum w_j=1$, ridge regression, the aforementioned `winner-takes-all' approach that selects the base learner with the smallest cross-validated mean-squared error, and unconstrained least squares.
Stacking can be applied in a similar way to classification problems. \stcmd{pystacked}\xspace supports stacking for binary classification problems where the outcome $y_i$ takes the values 0 or 1. The main difference to stacking regression is that $\hat{y}^{-k}_{(j),i}$ represent cross-validated predicted probabilities.
In the following paragraphs, we briefly describe the base learners supported by \stcmd{pystacked}\xspace and highlight central tuning parameters. We repeat that each of the machine learners discussed below can be fit using \stcmd{pystacked}\xspace as a regular stand-alone machine learner without the stacking layer.
Note that it goes beyond the scope of the article to describe each learner in detail. Familiarity with linear regression, logistic regression, and classification and regression trees is assumed. We recommend consulting machine learning textbooks, e.g.\ Hastie2009, for more detailed discussion.
\paragraph{Regularized regression} imposes a penalty on the size of coefficients to control over-fitting. The lasso penalizes the absolute size of coefficients, whereas the ridge penalizes the sum of squared coefficients. Both methods shrink coefficients toward zero, but only the lasso yields sparse solutions where some coefficient estimates are set to exactly zero. The elastic net combines lasso and ridge-type penalties. For classification tasks with a binary outcome, logistic versions of lasso, ridge and elastic net are available.
The severity of the penalty is most commonly chosen by cross-validation. For lasso only, \stcmd{pystacked}\xspace also supports selecting the penalty by AIC or BIC. The use of AIC or BIC has the advantage that it is computationally less intensive than cross-validation. Ahrens2019 compare the two approaches.
\paragraph{Random forests} rely on fitting a large number of regression or decision trees on bootstrap samples of the data. The random forest prediction is obtained as the average across individual trees. A crucial aspect of random forests is that, at each split when growing a tree, one may consider only a random subset of predictors. This restriction aims at de-correlating the individual trees. Central tuning parameters are the number of trees (\verb+n_estimators()+), the maximum depth of individual trees (\verb+max_depth()+), the minimum number of observations per leaf (\verb+min_samples_leaf()+), the number of features to be considered at each split (\verb+max_features()+) and the size of the bootstrap samples (\verb+max_samples()+).
\paragraph{Gradient boosted trees} also rely on fitting a large number of trees. In contrast to random forests, these trees are fit sequentially to the residuals from the current model. The learning rate determines how much the latest tree contributes to the overall model. Individual trees are usually fit to the whole sample, although sub-sampling is possible. In addition to tuning parameters relating to the trees, the learning rate (\verb+learnings_rate()+) and the number of trees (\verb+n_estimators()+) are the most important tuning parameters.
\paragraph{Support vector machines (SVM).} Support vector classifiers span a hyperplane that separates observations by their outcome class. The hyperplane is chosen to maximize the distance (margin) to correctly classified observations while allowing for some classification errors. The tuning parameter $C$ (\verb+C()+) controls the frequency and degree of classification mistakes. The hyperplane can be either linear or fitted using kernels. The SVM algorithm can also be adapted for regression tasks. To this end, the hyperplane is constructed to include as many observations as possible in a tube of size $2\epsilon$ around the hyperplane. Central tuning parameters for regression are $\epsilon$ (\verb+epsilon()+) and $C$ (\verb+C()+), which determines the cost of observations outside of the tube.
\paragraph{Feed-forward neural networks} consist of hidden layers that link the predictors (referred to as input layers) to the outcome. Each hidden layer is composed of multiple units (nodes) which pass signals to the next layer using an activation function. Central tuning parameters are the choice of the activation function (\verb+activation()+), and the number and size of hidden layers (\verb+hidden_layer_sizes()+). Further tuning choices relate to stochastic gradient descent algorithms which are typically used to fit neural networks. The default solver is Adam kingma2014adam. The option \verb+early_stopping+ can be used to set aside a random fraction of the data for validation. The optimization algorithm stops if there is no improvement in performance over a pre-specified number of iterations (see related options \verb+n_iter_no_change()+ and \verb+tol()+).
This section introduces the program \stcmd{pystacked}\xspace and its main features. \stcmd{pystacked}\xspace offers two alternative syntaxes between which the user can choose (see Section (ref) and (ref)). The two syntaxes offer the same functionality and are included to accommodate different user preferences. Section (ref) and (ref) list post-estimation commands and general options, respectively. Section (ref) discusses supported base learners. Section (ref) is a note on learner-specific predictors. Section (ref) explains the pipeline feature.
The first syntax uses {\tt\underbar{m}ethod}(\ststring) to select base learners, where {\it string} is a list of base learners. Options are passed on to base learners via cmdopt1(string), cmdopt2(string), etc. That is, base learners can be specified and options are passed on in the order in which they appear in {\tt\underbar{m}ethod}(\ststring) (see Section (ref)). Likewise, the \texttt{pipe*(}\emph{string}\texttt{)} option can be used for pre-processing predictors within Python on the fly, where `\texttt{*}' is a placeholder for `\texttt{1}', `\texttt{2}', etc. (see Section (ref)). Finally, \texttt{xvars*(}\emph{predictors}\texttt{)} allows specifying a learner-specific variable lists of predictors.
where {\it general_options} are discussed in Section (ref) below.
In the second syntax, base learners are added before the comma using {\tt\underbar{m}ethod}(\ststring) along with further learner-specific settings and separated by `\verb+||+'.
\paragraph{Predicted values.} To get predicted values:
To get fitted values for each base learner:
Predicted values (in- and out-of-sample) are calculated when \stcmd{pystacked}\xspace is run and stored in Python memory. predict pulls the predicted values from Python memory and saves them in Stata memory. This storage structure means that no changes on the data in Stata memory should be made between the \stcmd{pystacked}\xspace call and the predict call. If changes to the data set are made, predict will return an error.
The option {\tt basexb} returns predicted values for each base learner. By default, the predicted values from re-fitting base learners on the full estimation sample are returned. If combined with {\tt cvalid}, the cross-fitted predicted values are returned for each base learner.
\paragraph{Tables.} After estimation, \stcmd{pystacked}\xspace can report a table of in-sample (both cross-validated and full-sample) and, optionally, out-of-sample (or holdout sample) performance for both the stacking regression and the base learners. For regression problems, the table reports the root mean-squared prediction error (RMSPE). For classification problems, a confusion matrix is reported. The default holdout sample used for out-of-sample performance with the holdout option is all observations not included in the estimation. Alternatively, the user can specify the holdout sample explicitly using the syntax holdout(varname). The table can be requested after estimation as a replay command or as part of the \stcmd{pystacked}\xspace estimation.
\paragraph{Graphs.} \stcmd{pystacked}\xspace can also create graphs of in-sample and, optionally, out-of-sample performance for both the stacking regression and the base learners. For regression problems, the graphs compare predicted and actual values of depvar. For classification problems, the default is to generate receiver operator characteristic (ROC) curves. Optionally, histograms of predicted probabilities are reported. As with the {\tt table} option, the default holdout sample used for out-of-sample performance is all observations not included in the estimation, but the user can instead specify the holdout sample explicitly. The table can be requested after estimation or as part of the \stcmd{pystacked}\xspace estimation command. The {\tt graph} option on its own reports the graphs using \stcmd{pystacked}\xspace's default settings. Because graphs are produced using Stata's twoway, roctab and histogram commands, the user can control either the combined graph (graph(options\texttt{)}) or the individual learner graphs (\texttt{lgraph(}\emph{options}\texttt{)}) by passing options to these commands.
A full list of general options is provided in the \stcmd{pystacked}\xspace help file. We list only the most important general options here:
The base learners are chosen using the option \uline{m}ethod(string) in combination with type(string). The latter can take the value \emph{reg(ress)} for regression and \emph{class} for classification problems. Table (ref) provides an overview of supported base learners and their underlying {\it scikit-learn}\xspace routines.
cmdopt*(string) (Syntax 1) and opt(string) (Syntax 2) are used to pass options to the base learners. Due to the large number of options, we do not list all options here. We instead provide a tool that lists options for each base learner. For example, to get the default options for lasso with cross-validated penalty, type
The naming of the options follows {\it scikit-learn}\xspace. Allowed settings for each option can be inferred from the {\it scikit-learn}\xspace documentation. We strongly recommend that the user reads the {\it scikit-learn}\xspace documentation carefully.\footnote{The {\it scikit-learn}\xspace documentation is available at \url{https://scikit-learn.org/stable/index.html} (last accessed on August 18, 2022).}
By default, \stcmd{pystacked}\xspace uses the same set of predictors for each base learner. Using the same predictors for each method is often not desirable as the optimal set of predictors may vary across base learners. For example, when using linear machine learners such as the lasso, adding polynomials, interactions and other transformations of the base set of predictors might greatly improve out-of-sample prediction performance. The inclusion of transformations of base predictors is especially worth considering if the base set of observed predictors is small (relative to the sample size) and the relationship between outcome and predictors is likely non-linear. Tree-based methods (e.g., random forests and boosted trees), on the other hand, can detect certain types of non-linear patterns automatically. While adding transformations of the base predictors may still lead to performance gains, the added benefit is less striking relative to linear learners and might not justify the additional costs in terms of computational complexity.
There are two approaches to implement learner-specific sets of predictors: Pipelines, discussed in the next section, can be used to create some transformations on the fly for specific base learners. A more flexible approach is the xvars*(predictors) option, which allows specifying predictors for a particular learner. xvars*() supports standard Stata factor variable notation.
{\it scikit-learn}\xspace uses pipelines to pre-preprocess input data on the fly. In \stcmd{pystacked}\xspace, pipelines can be used to impute missing values, create polynomials and interactions, and to standardize predictors. Table (ref) lists the pipelines currently supported by \stcmd{pystacked}\xspace.
\paragraph{Remarks.} First, regularized regressors (i.e., the methods {\it lassoic}, {\it lassocv}, {\it ridgecv} and {\it elasticcv}) use the stdscaler pipeline by default. pipe*(nostdscaler) disables this behavior. Second, the stdscaler0 pipeline is useful in combination with sparse, which transforms the predictor matrix into a sparse matrix. \emph{stdscaler0} does not center predictors so that the predictor matrix retains its sparsity property.
This section demonstrates how to apply \stcmd{pystacked}\xspace for regression and classification tasks. Before we discuss stacking, we first show how to use \stcmd{pystacked}\xspace as a `regular' machine learning program for fitting a single supervised machine learner (see next sub-section). We then illustrate stacking regression and stacking classification in Section (ref) and (ref), respectively.
We import the Kelley1997 California house price data, and split the sample randomly into training and validation partitions using a 75/25 split. The aim of the prediction task is to predict median house prices ({\tt medhousevalue}) using a set of house price characteristics. \footnote{The following predictors are included in the data set in this order: district longitude (longitude), latitude (latitude), median house age (houseage), average number of rooms per household (rooms), average number of bedrooms per household (bedrooms), block group population (population), average number of household members (\texttt{households}), median income in block group (\texttt{medinc}). The data set was retrieved from the StatLib repository at \verb|https://www.dcc.fc.up.pt/ ltorgo/Regression/cal_housing.html| (last accessed on February 19, 2022). There are 20,640 observations in the data set. } We prepare the data for analysis as follows.
\paragraph{Gradient-boosted trees.} As a first example, we use \stcmd{pystacked}\xspace to fit gradient-boosted regression trees and save the out-of-sample predicted values. Since we consider a regression task rather than a classification task, we specify type(reg) (which is also the default). The option method(gradboost) selects gradient boosting. We will later see that we can specify more than one learner in methods(), and that we can also fit gradient-boosted classification trees.
The output shows the stacking weights associated with each base learner. Since we only consider one method, the output is not particularly informative and simply shows a weight of one for gradient boosting. Yet, \stcmd{pystacked}\xspace has fitted 100 boosted trees (the default) in the background using {\it scikit-learn}\xspace's ensemble.\-Gradient\-Boosted\-Regressor.
Before we tune our gradient boosting learner, we retrieve a list of available options. The default options for gradient boosting can be listed in the console, and the {\it scikit-learn}\xspace documentation provides more detail on the allowed parameters of each option.
We consider two additional specifications. Note that we restrict ourselves to a few selected specifications for illustrative purposes. We stress that careful parameter tuning should consider a grid of values across multiple learner parameters. The first specification reduces the learning rate from 0.1 (the default) to 0.01. The second specification reduces the learning rate and increases the number of trees from 100 (the default) to 1000. We use cmdopt1() since gradient boosting is the first (and only) method listed in methods().
We can then compare the performance across the three models using the out-of-sample mean-squared prediction error (MSPE):
The initial gradient booster achieves an out-of-sample MSPE of 29.87. The second gradient booster uses a reduced learning rate of 0.01 and performs much worse, with an MSPE of 71.10. The third gradient booster performs only slightly worse than the first, illustrating the trade-off between the learning rate and the number of trees.
\paragraph{Pipelines.} We can make use of pipelines to pre-process our predictors. This is especially useful in the context of stacking when we want to, for example, use second-order polynomials of predictors as inputs for one method, but only use elementary predictors for another method. Here, we compare lasso with and without the poly2 pipeline:
We could replace { pipe1(\tt poly2)} with xvars1(c.(medh longi-medi)\#\#c.(medh longi -medi)). In fact, the latter is more flexible and allows, for example, to create interactions for some predictors and not for others.
We again calculate the out-of-sample MSPE:
The poly2 pipeline improves the performance of the lasso, indicating that squared and interaction terms constitute important predictors. However, the lasso does not perform as well as gradient boosting in this application.
We now consider a stacking regression application with five base learners: (1) linear regression, (2) lasso with penalty chosen by cross-validation, (3) lasso with second-order polynomials and interactions, (4) random forest with default settings, and (5) gradient boosting with a learning rate of 0.01 and 1000 trees. That is, we use the lasso twice---once with and once without the poly2 pipeline. Indeed, nothing keeps us from using the same algorithm multiple times. This way, we can combine the same algorithm with different settings.
Note the numbering of the pipe*() and cmdopt*() options below. We apply the \emph{poly2} pipe to the first and third methods (\emph{ols} and \emph{lassoic}). We also change the default learning rate and the number of estimators for gradient boosting (the 5th estimator).
The above syntax becomes a bit difficult to read with many methods and many options. We offer an alternative syntax for easier use with many base learners.
The stacking weights shown in the output determine how much each method contributes to the final stacking predictor. In this example, OLS and lasso based on both sets of predictors all receive a weight of zero. Random forest receives a weight of 83.8%, and gradient boosting contributes the remaining 16.2% of the weight to the final predictor.
\paragraph{Predicted values.} In addition to the stacking predicted values, we can also get the predicted values of each base learner using the {\tt basexb} option:
\paragraph{Plotting.} \stcmd{pystacked}\xspace also comes with plotting features. The graph option creates a scatter plot of predicted values on the vertical and observed values on the horizontal axis for stacking and each base learner, see Figure (ref). The black line is a 45-degree line and shown for reference. Since \stcmd{pystacked}\xspace with graph can be used as a post-estimation command, there is no need to re-run the stacking estimation.
Figure (ref) shows the out-of-sample predicted values. To see the in-sample predicted values, simply omit the holdout option. Note that the holdout option will not work if the estimation was run on the whole sample.
\stcmd{pystacked}\xspace can be applied to binary classification problems. For demonstration, we consider the Spambase Data Set of cranor1998spam, which we retrieve from the UCI Machine Learning Repository. We load the data and split the data into training (75%) and validation sample (%25).
The example below is more complicated. We go through it step-by-step:
As in the previous regression example, gradient boosting receives the largest stacking weight and thus contributes most to the final stacking prediction.
\paragraph{Confusion matrix.} Confusion matrices allow comparing actual and predicted outcomes in a $2\times 2$ matrix. \stcmd{pystacked}\xspace provides a compact table format that combines confusion matrices for each base learner and for the final stacking classifier, both for the training and validation partition.
For example, the table shows 29 false positives for stacking and 294 for logistic regression in the validation partition, while the number of false negatives is 29 and 32, respectively. The accuracy of stacking and logistic regression is thus given by $(678+397)/1133=94.9$% and $(413+394)/1133=71.2$%.
\paragraph{Plotting.} \stcmd{pystacked}\xspace supports ROC curves which allow assessing the classification performance for varying discrimination thresholds. The $y$-axis in a ROC plot corresponds to sensitivity (true positive rate) and the $x$-axis corresponds to 1--specificity (false positive rate). The Area Under the Curve (AUC) displayed below each ROC plot is a common evaluation metric for classification problems.
In this article, we introduce the Stata program pystacked. The program not only makes a range of popular supervised machine learners available, such as regularized regression, random forests, and gradient-boosted trees, but also facilitates combining multiple learners for stacking regression or classification. pystacked comes with a number of practically relevant features: for example, sparse matrix support, learner-specific predictors, pipelines for predictor transformations, and two alternative syntaxes. Nevertheless, we also see the potential for improvement in at least two directions. First, pystacked offers few diagnostics for individual learners. Features such as variable importance plots and coefficient estimates for parametric learners could be added in later versions. Secondly, the deep learning features of pystacked are currently relatively limited and could be improved by adding support for other deep learning algorithms.
We thank Jan Ditzen, Ben Jann, Blaise Melly, Alessandro Oliveira, Matthias Schonlau, and Thomas Wiemann for their helpful feedback. All remaining errors are our own.