The NN-Stacking: Feature weighted linear stacking through neural networks

TL;DR

The paper introduces NN-Stacking, a neural network-based method that dynamically adjusts feature weights in linear stacking models, enhancing predictive accuracy by capturing regional variations in model performance.

Contribution

It generalizes traditional linear stacking by allowing coefficients to vary with input features using neural networks, improving prediction in diverse data regions.

Findings

Improved predictive accuracy over traditional stacking methods.

Effective in large sample datasets with complex feature interactions.

Maintains interpretability at a local level.

Abstract

Stacking methods improve the prediction performance of regression models. A simple way to stack base regressions estimators is by combining them linearly, as done by \citet{breiman1996stacked}. Even though this approach is useful from an interpretative perspective, it often does not lead to high predictive power. We propose the NN-Stacking method (NNS), which generalizes Breiman's method by allowing the linear parameters to vary with input features. This improvement enables NNS to take advantage of the fact that distinct base models often perform better at different regions of the feature space. Our method uses neural networks to estimate the stacking coefficients. We show that while our approach keeps the interpretative features of Breiman's method at a local level, it leads to better predictive power, especially in datasets with large sample sizes.

Equations

$$G_{\theta }(\mathbf{x}):=\sum _{i=1}^k{\theta }_i{g}_i(\mathbf{x})={\theta }^{\prime }{g}_{\mathbf{x}}$$

$$\arg \underset{\theta }{\min }\sum _{i=1}^n(y_i-G_{\theta }^{(-i)}({\mathbf{x}}_i){)}^2=\arg \underset{\theta }{\min }\sum _{i=1}^n(y_i-{\theta }^{\prime }{g}_{{\mathbf{x}}_i}^{(-i)}{)}^2,$$

$$G_{\theta }(\mathbf{x}):=\sum _{i=1}^k{\theta }_i(\mathbf{x})g_i(\mathbf{x})={\theta }_x^{\prime }{g}_{\mathbf{x}},$$

$$R(G_{\theta })=\text{E}\left[{\left(Y-G_{\theta }(\mathbf{X})\right)}^2\right],$$

$$\frac{1}{n}\sum _{k=1}^n(G_{\theta }({\mathbf{x}}_k)-y_k{)}^2.$$

$$G_{\theta }^{\prime }(\mathbf{x})={\theta }_{\mathbf{x}}^{\prime }{g}_{\mathbf{x}}+\varphi (\mathbf{x}).$$

$$\arg \underset{{\theta }_x}{\min }R(G_{\theta })$$

$${\theta }_{\mathbf{x}}=\frac{{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}}{{\mathbf{e}}^{\prime }{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}},$$

$${\mathbb{M}}_{\mathbf{x}}$$

$$=\text{E}[Y^2|\mathbf{x}]-\text{E}[Y|\mathbf{x}](g_i(\mathbf{x})-g_j(\mathbf{x})))+g_i(\mathbf{x})g_j(\mathbf{x}).$$

$$R(G_{\theta })=\text{E}\left[\text{E}\left[{\left(Y-G_{\theta }(\mathbf{X})\right)}^2\right|\mathbf{X}]\right].$$

$$\text{E}\left[{\left(Y-G_{\theta }(\mathbf{x})\right)}^2|\mathbf{X}=\mathbf{x}\right]$$

$$=\sum _{i,j}{\theta }_i(\mathbf{x}){\theta }_j(\mathbf{x})\text{E}\left[(Y-g_i(\mathbf{x}))(Y-g_j(\mathbf{x}))|\mathbf{X}=\mathbf{x}\right]$$

$$={\theta }_{\mathbf{x}}^t{\mathbb{M}}_{\mathbf{x}}{\theta }_{\mathbf{x}},$$

$$f({\theta }_{\mathbf{x}},\lambda ):={\theta }_{\mathbf{x}}^t{\mathbb{M}}_{\mathbf{x}}{\theta }_{\mathbf{x}}-\lambda ({\mathbf{e}}^t{\theta }_{\mathbf{x}}-1).$$

$$\frac{\partial f({\theta }_{\mathbf{x}},\mathbf{x})}{\partial {\theta }_{\mathbf{x}}}=2{\mathbb{M}}_{\mathbf{x}}{\theta }_{\mathbf{x}}-\lambda \mathbf{e},$$

$$f({\theta }_{\mathbf{x}}^{\ast },\lambda )=-\frac{{\lambda }^2}{4}{\mathbf{e}}^t{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}+\lambda ,$$

$$\frac{\partial f({\theta }_{\mathbf{x}}^{\ast },\lambda )}{\partial \lambda }=0⟺\lambda =\frac{2}{{\mathbf{e}}^t{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}},$$

$${\theta }_{\mathbf{x}}^{\ast }=\frac{\lambda }{2}{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}=\frac{2}{{\mathbf{e}}^t{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}}\frac{1}{2}{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}=\frac{{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}}{{\mathbf{e}}^t{\mathbb{M}}_{\mathbf{x}}^{-1}\mathbf{e}}.$$

Full text

1 Introduction

The standard procedures for model selection in prediction problems is cross-validation and data splitting. However, such an approach is known to be sub-optimal [2, 3, 4]. The reason is that one might achieve more accurate predictions by combining different regression estimators rather then by selecting the best one. Stacking methods [5] are a way of overcoming such a drawback from standard model selection.

A well known stacking method was introduced by Breiman [1]. This approach consists in taking a linear combination of base regression estimators. That is, the stacked regression has the shape $\sum_{i=1}^{k}\theta_{i}g_{i}({\mathbf{x}})$, where $g_{i}$’s are the individual regression estimators (such as random forests, linear regression or support vector regression), $\theta_{i}$ are weights that are estimated from data and ${\mathbf{x}}$ represents the features.

Even though this linear stacking method leads to combined estimators that are easy to interpret, it may be sub-optimal in cases where models have different local accuracy, i.e., situations where the performance of these estimators vary over the feature space. Example 1.1 illustrates this situation.

Example 1.1**.**

Consider predicting $Y$ based on a single feature, $x$, using the data in Figure 1. We fit two least squares estimators: $g_{1}({\mathbf{x}})=\theta_{01}+\theta_{11}x$ and $g_{2}({\mathbf{x}})=\theta_{02}+\theta_{12}x^{2}$. None of the models is uniformly better; for example, the linear fit has better performance when $x\leq 5$, but the quadratic fit yields better performance for $x\in(-2.5,2.5)$. One may take this into account when creating the stacked estimator by assigning different weights for each regression according to $x$: while one can assign a larger weight to the linear fit on the regime $x\leq 5$, a lower weight should be assigned to it if $x\in(-2.5,2.5)$.

It is well known that different regression methods may perform better on different regions of the feature space. For instance, because local estimators do not suffer from boundary effects, they achieve good performance closer to the edges of the feature space [6]. Random forests, on the other hand, implicitly perform feature selection, and thus may have better performance in regions where some features are not relevant [7].

In this work we improve Breiman’s approach so that it can take local accuracy into account. That is, we develop a meta-learner that is able to learn which models have higher importance on each region of the feature space. We achieve this goal by allowing each parameter $\theta_{i}$ to vary as a function of the features ${\mathbf{x}}$. In this way, the meta-learner can adapt to each region of the feature space, which yields higher predictive power. Our approach keeps the local interpretability of the linear stacking model.

The remaining of the work is organized as follows. Section 2 introduces the notation used in the paper, as well as our method. Section 3 shows details on its implementation. Section 4 shows applications of our method to a variety of datasets to evaluate its performance. Section 5 concludes the paper.

2 Notation and Motivation

The stacking method proposed by Breiman [1] is a linear combination of $k$ regression functions for a label $Y\in\mathbb{R}$. More precisely, let $g_{\mathbf{x}}=(g_{1}({\mathbf{x}}),g_{2}({\mathbf{x}}),\ldots,g_{k}({\mathbf{x}}))^{\prime}$ be a vector of regression estimators, that is, $g_{i}({\mathbf{x}})$ is an estimate of ${\textbf{E}}[Y|{\mathbf{x}}],\forall\,i=1,2,\ldots,k$. The linear stacked regression is defined as

[TABLE]

where $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{k})^{\prime}$ are meta-parameters. One way to estimate the meta-parameters using data $({\mathbf{x}}_{1},y_{1}),\ldots,({\mathbf{x}}_{n},y_{n})$ is through the least squares method, computed using a leave-one-out setup:

[TABLE]

where $g_{j}^{(-i)}({\mathbf{x}}_{i})$ is the prediction for ${\mathbf{x}}_{i}$ made by the $j$-th regression fitted without the $i$-th instance. Note that it is important to use this hold-out approach because if the base regression functions $g_{1}({\mathbf{x}}),g_{2}({\mathbf{x}}),\ldots,g_{k}({\mathbf{x}})$ are constructed using the same data as $\theta_{1},\ldots,\theta_{k}$, this can cause $G_{\theta}({\mathbf{x}})$ to over-fit the training data.

In order for the stacked estimator to be easier to interpret, Breiman [1] also requires $\theta_{i}$’s to be weights, that is $\theta_{i}\geq 0\,\forall\,i=1,2,\ldots,k$ and that $\sum_{i=1}^{k}\theta_{i}=1$.

Even though Breiman’s solution works on a variety of settings, it does not take into account that each regression method may perform better in distinct regions of the feature space. In order to overcome this limitation, we propose the Neural Network Stacking (NNS) which generalizes Breiman’s approach by allowing $\theta$ on Equation 2.1 to vary with ${\mathbf{x}}$. That is, our meta-learner has the shape

[TABLE]

where $\theta_{x}:=(\theta_{1}({\mathbf{x}}),\theta_{2}({\mathbf{x}}),\ldots,\theta_{k}({\mathbf{x}}))^{\prime}$. In other words, the NNS is a local linear meta-learner. Example 2.1 shows that NNS can substantially decrease the prediction error of Breiman’s approach.

Example 2.1**.**

We fit both Breiman’s linear meta-learner and our NNS local linear meta-learner to the models fitted in Example 1.1. Figure 2 shows that Breiman’s meta-learner is not able to fit the true regression satisfactorily because both estimators have poor performance on specific regions of the data. On the other hand, feature-varying weights yield a better fit.

3 Methodology

Our goal is to find $\theta_{\mathbf{x}}=(\theta_{1}({\mathbf{x}}),\ldots,\theta_{k}({\mathbf{x}}))^{\prime}$, $\theta_{i}:\mathcal{X}\longrightarrow\mathbb{R}$, that minimizes the mean squared risk,

[TABLE]

where $G_{\theta}({\mathbf{x}})$ is defined as in Equation in 2.3.

We estimate $\theta_{\mathbf{x}}$ via an artificial neural network. This network takes ${\mathbf{x}}$ as input and produces an output $\theta_{\mathbf{x}}$, which is then used to obtain $G_{\theta}({\mathbf{x}})$. To estimate the weights of the networks, we introduce an appropriate loss function that captures the goal of having a small $R(G_{\theta})$. This is done by using the loss function

[TABLE]

Notice that the base regression estimators are used only when evaluating the loss function; they are not the inputs of the network. With this approach, we allow each $\theta_{i}({\mathbf{x}})$ to be a complex function of the data. We call this method Unconstrained Neural Network Stacking (UNNS). Figure 3 illustrates a UNNS that stacks 2 base estimators in a regression problem with four features.

In addition to the linear stacking, this approach allows the user to easily take advantage of the neural network architecture by directly adding a network output node, $\phi({\mathbf{x}})$, to the stacking. That is, we also consider a variation of UNNS which takes the shape

[TABLE]

This has some similarity to adding a single neural network estimator to the stacking. However, we use the same architecture to create the additional term, mitigating computation time. Algorithm 1 shows how this method is implemented. In order to avoid over-fitting, $\theta_{i}$’s and $g_{i}$’s are estimated using different folds of the training set.

In order to achieve an interpretable stacked solution, we follow Breiman’s suggestion and consider a second approach to estimate $\theta_{i}$’s which consists in minimizing $R(G_{\theta})$ under the constrain that $\theta_{i}$’s are weights, that is, $\theta_{i}({\mathbf{x}})\geq 0$ and $\sum_{i=1}^{k}\theta_{i}({\mathbf{x}})=1$. Unfortunately, it is challenging to directly impose this restriction to the solution of the neural network. Instead, we use a different parametrization of the problem, which is motivated by Theorem 3.1.

Theorem 3.1**.**

The solution of

[TABLE]

under the constrain that $\theta_{i}({\mathbf{x}})\geq 0$ and $\sum_{i=1}^{k}\theta_{i}({\mathbf{x}})=1$ is given by

[TABLE]

where $\mathbf{e}$ is a $k$-dimensional vector of ones and

[TABLE]

with $(i,j)\in\{1,...,k\}^{2}$.

Theorem 3.1 shows that, under the given constrains, $\theta({\mathbf{x}})$ is uniquely defined by $\mathbb{M}^{-1}_{\mathbf{x}}$. Now, because $\mathbb{M}_{\mathbf{x}}$ is a covariance matrix, then $\mathbb{M}^{-1}_{\mathbf{x}}$ is positive definite, and thus Cholesky decomposition can be applied to it. It follows that $\mathbb{M}^{-1}_{\mathbf{x}}=L_{\mathbf{x}}L_{\mathbf{x}}^{\prime}$, where $L_{\mathbf{x}}$ is a lower triangular matrix. This suggests that we estimate $\theta_{\mathbf{x}}$ by first estimating $L_{\mathbf{x}}$ and then plugging the estimate back into Equation 3.1. That is, in order to obtain a good estimator under the above mentioned restrictions, the output of the network is set to be $L_{\mathbf{x}}$ rather than the weights themselves111Since the gradients for all matrix operations are implemented for Pytorch tensor classes, the additional operations of the CNNS method will be automatically backpropagated once Pytorch’s backward method is called on the loss evaluation.. We name this method Constrained Neural Network Stacking (CNNS). Figure 4 illustrates a CNNS that stacks 2 base regressors (that is, $L_{\mathbf{x}}=[l_{ij}]$ is a $2x2$ triangular matrix) in a 4 feature regression problem.

Algorithm 2 shows the implementation of this method. As with UNNS, we also explore a variation which adds an extra network output $\phi({\mathbf{x}})$ to $G_{\theta}$.

Figure 5 illustrates the full training process. For simplicity, the neural network early stopping patience criterion is set to a single epoch and the additional parameter $\phi_{\mathbf{x}}$ is not used.

3.1 Comparison with standard stacking methods

Most stacking methods create a meta-regression model by applying a regression method directly on the outputs of individual predictions. In particular, a meta-regression method can be a neural network. Such procedure differs from NN-Stacking by the shape of both the input and of the output of the network. While standard stacking uses base regression estimates ($g_{x}$) as input and $Y$ as output, NN-Stacking uses the features as input and either the weights $\theta_{x}$ (for UNSS) or $L_{\mathbf{x}}$ (for CNSS) as outputs. The base regression estimates are used only on the loss function. Thus, the NN-Stacking method leads to more interpretable models. Section 4 compares these methods in terms of their predictive power. We also point out that our approach has some similarity to Sill et al. [4], which allows each $\theta_{i}$ to depend on meta-features computed from ${\mathbf{x}}$ using a specific parametric form. Neural networks, on the other hand, provide a richer family of functions to model such dependencies (in fact, they are universal approximators; Csáji [8]).

3.2 Selecting base regressors

Consider the extreme case where $g_{i}({\mathbf{x}})=g_{j}({\mathbf{x}})\,\forall\,{\mathbf{x}}\in\mathcal{X}$ for some $i\neq j$, that is, the case in which two base regressors generate the same prediction over all feature space. Now, suppose that one fits a NNS (either CNNS or UNNS) for this case. Then $\theta_{i}g_{i}({\mathbf{x}})+\theta_{j}g_{j}({\mathbf{x}})=(\theta_{i}+\theta_{j})g_{i}({\mathbf{x}}).$ Thus, one of the regressions can be dropped from the stacking with no loss in predictive power.

In practice, our experiments (Section 4) show that regression estimators that have strongly correlated results do not contribute to the meta-learner. This suggests that one should choose base regressors with considerably distinct nature.

3.3 Implementation details

A Python package that implements the methods proposed in this paper is available at github.com/randommm/nnstacking. The scripts for the experiments in Section 4 are availiable at github.com/vcoscrato/NNStacking. We work with the following specifications for the artificial neural networks:

• •

Optimizer: we use the Adam algorithm [9] and decrease its learning rate after the validation loss stops improving for a user-defined number of epochs.

• •

Initialization: we use the Xavier Gaussian method proposed by Glorot and Bengio [10] to sample the initial parameters of the neural network.

• •

Layer activation and regularization: we use ELU [11] as the activation function, and do not use regularization.

• •

Normalization: we use batch normalization [12] to speed-up the training process.

• •

Stopping criterion: in order to address the risk of having strong over-fit on the neural networks, we worked with a 90%/10% split early stopping for small datasets and a higher split factor for larger datasets (increasing the proportion of training instances) and a patience of 10 epochs without improvement on the validation set.

• •

Dropout: We use dropout (with a rate of 50%) to address the problem of over-fitting [13].

• •

Software: we use PyTorch [14].

• •

Architecture: as default values we use a 3 layer depth network with hidden layer size set to 100; these values have been experimentally found to be suitable in our experiments (Section 4).

4 Experiments

We compare stacking methods for the following UCI datasets:

• •

The GPU kernel performance dataset (241600 instances, 13 features) [15],

• •

The music year prediction dataset [16] (515345 instances and 90 features),

• •

The blog feedback dataset [17] (60021 instances, 280 features),

• •

The superconductivity dataset [18] (21263 instances, 80 features).

First, we fit the following regression estimators (that will be stacked):

• •

Three linear models: with L1, L2, and no penalization [19],

• •

Two tree based models: bagging and random forests [19],

• •

A gradient boosting method (GBR) [20].

The tuning parameters of these estimators are chosen by cross-validation using scikit-learn [21].

Using these base estimators, we then fit four variations of NNS (both CNNS and UNNS with and without the additional $\phi_{x}$) using the following specifications:

• •

Tuning: three different architectures were tested for each neural network approach. The layer size was fixed at 100 and the number of hidden layers were set to 1, 3, and 10. We choose the architecture with the lowest validation mean-squared error.

• •

Train/validation/test split: for all datasets, we use 75% of the instances to fit the models, among which 10% are used for performing early stop. The remaining 25% of the instances are used as a test set to compare the performance of the various models. The train/test split is performed at random. The cross-validated predictions (the matrix $P$ denoted on Algorithm 1) are obtained using a 10-fold cross-validation on the training data (i.e., $F=10$).

• •

Total fitting time: we compute the total fitting time (in seconds; including the time for cross-validating the network architecture) of each method on two cores of an AMD Ryzen 7 1800X processor running at 3.6Gz.

We compare our methods with Breiman’s linear stacking and the usual neural net stacking model described in Section 3.1. In addition to these, we also include a comparison with a direct neural network that has ${\mathbf{x}}$ as its input and $Y$ as its output.

The comparisons are made by evaluating the mean squared error (MSE, $n^{-1}\sum_{i=1}^{n}(y_{i}-g({\mathbf{x}}_{i}))^{2}$) and the mean absolute error (MAE, $n^{-1}\sum_{i=1}^{n}|y_{i}-g({\mathbf{x}}_{i})|$) of each model $g$ on a test set. We also compute the standard error for each of these metrics, which enables one to compute confidence intervals for the errors of each method.

4.1 GPU kernel performance dataset

Table 1 shows the results that were obtained for the GPU kernel performance dataset. Our UNNS methods outperforms both Breiman’s stacking and the usual meta-regression stacking approaches in terms of MSE. Moreover, the UNNS model is also the best one in terms of MAE, even though the gap between the models is lower in this case. Our stacking methods also perform better than all base estimators. This suggests that each base model performs better on a distinct region of the feature space.

Figure 6 shows a boxplot with the distribution of the fitted $\theta_{i}$’s for UNNS. Many fitted values fall out of the range $[0,1]$, which explains why UNNS gives better results than Breiman’s and CNNS (which have the restriction that $\theta_{i}$’s must be proper weights).

Table 2 shows the correlation between the prediction errors for base estimators. The linear estimators had an almost perfect pairwise correlation, which indicates that removing up to 2 of them from the stacking would not affect predictions. Indeed, after refitting UNNS without using ridge regression and lasso, we obtain exactly the same results. We also refit the best UNNS removing all of the linear estimators to check if poor performing estimators are making stacking results worse. In this setting, we obtain an MSE of $11074.13(\pm 227.57)$, and a MAE of $45.76(\pm 0.39)$. Note that although the point estimates of the errors are lower than those obtained in Table 1, the confidence intervals have an intersection, which leads to the conclusion that the poor performance of linear estimators is not damaging the stacked estimator.

4.2 Music year dataset

Table 3 shows the accuracy metrics results for the music year dataset. In this case, the CNNS gave the best results, both in terms of MSE and MAE. For this dataset, Breiman’s stacking was worse than using gradient boosting, one of the base regressors. The same happens with the usual meta-regression neural network approach. On the other hand, NNS could find a representation that combines the already powerful GBR estimator with less powerful ones in a way that leverages their individual performance.

All base estimators had high prediction error correlations (Table 4). In particular, two of the linear estimators could be removed from the stacking without affecting its performance. However, when removing all three linear estimators the MSE for the best NNS increased to $83.92(\pm 0.57)$ and its MAE increased to $6.44(\pm 0.02)$.

Figure 7 shows that the fitted NNS weights have a large dispersion. This illustrates the flexibility added by our method. Models with very distinctive nature (e.g., ridge regression - which imposes a linear shape on the regression function, and random forests - which is fully non-parametric) can add to each other, getting weights of different magnitudes depending on the region of the feature space that the new instance lies on.

4.3 Blog feedback dataset

Table 5 shows the results for the blog feedback dataset. All stacked estimators had similar performance in terms of MSE. However, UNNS had slightly worse performance with respect to MAE. This may happen because the NNS is designed to minimize the MSE and not the MAE. Overall, for this small dataset, the NNS shows no improvement over Breiman’s stacking or the usual meta-regression neural network.

GBR had the lowest MSE for the base estimators, while bagging and random forests had the lowest MAE. This explains why these models have larger fitted weights (Figure 8). Moreover, the linear models prediction errors had an almost perfect error correlation (Table 6). This suggests that removing up to 2 of them from the NNS would not impact its performance. Also, the linear estimators has a poor performance when compared to the other base regressors. We thus refit the best NNS for this data after removing these estimators, and achieve an MSE of 531.88 ($\pm 62.67$) and a MAE of 5.31 ($\pm 0.20$). We conclude that the linear estimators did not damage nor improved the NNS.

4.4 Superconductivity dataset

The results for the superconductivity dataset (Table 7) were similar to those obtained for the blog feedback data: the NNS methods perform slightly better than Breiman’s in terms of MSE, and worse in terms of MAE. Moreover, both tree-based models had the best MSE among base estimators, competing with the GBR in terms of MAE. Hence, they got larger fitted weights (Figure 9).

Table 8 shows that GBR did not have a high correlation error to the tree-based estimators (0.72 in both cases). This is another reason why although having higher MSE, the GBR has high fitted weight for some instances. One can also note that bagging and random forest had an almost perfect error correlation. This implies that removing one of them would lead to no changes in the NNS. Finally, removing the linear models did not change the MSE and the MAE for the stacking methods.

5 Conclusion and future extensions

NN-Stacking is a stacking tool with good predictive power that keeps the simplicity in interpretation of Breiman’s method. The key idea of the method is to take advantage of the fact that distinct base models often perform better at different regions of the feature space, and thus it allows the weight associated to each model to vary with ${\mathbf{x}}$.

Our experiments show that both CNNS and UNNS can be suitable in different settings: in cases where the base estimators do not capture the complexity from the whole data, the freedom adopted by UNNS can lead to a larger improvement in performance. On the other hand, when base estimators already have high performance, UNNS the CNNS have similar predictive power, but the restrictions imposed by CNNS guarantee a more interpretable solution. Both CNNS and UNNS have comparable computational cost.

In our experiments, we observe that NNS improves over standard stacking approaches especially on large datasets. This can be explained by the fact that NNS methods have a higher complexity (i.e., larger number of parameters) than the other approaches. Thus, a larger sample size is needed to satisfactorily estimate them. The experiments also show that including weak regression methods (such as linear methods) might decrease the errors of NNS. In a few cases, however, adding such weak regressors slightly increases the prediction errors of the stacked estimators This suggests that adding a penalization to the loss function that encourages $\theta_{i}$’s to be zero may lead to improved results.

Future work includes extending these ideas to classification problems, as well as developing a leave-one-out version based on super learners [22]. Also, we desire to develop a method of regularization on population moments estimation to avoid over-fitting, as well as to study asymptotic properties for the estimator of $L_{\mathbf{x}}$.

Acknowledgments

Victor Coscrato and Marco Inácio are grateful for the financial support of CAPES: this study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. Rafael Izbicki is grateful for the financial support of FAPESP (grant 2017/03363-8) and CNPq (grant 306943/2017-4). The authors are also grateful for the suggestions given by Rafael Bassi Stern.