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.
42,824 characters · 10 sections · 26 citation commands
Machine Collaboration
In recent decades, various learning methods, including deep neural networks (DNN), decision trees (DT, breiman1984classification), support vector machines cortes1995support, and k-nearest neighbors (kNN) have been developed for regression in supervised learning. As argued by hastie2009elements, each of these methods may have advantages over the others in some respects, but not in others. For example, while DNN is effective at approximating complicated nonlinear functions, the problems of overfitting and vanishing gradients could harm their performances , especially when the sample size is small. As another popular learning method, DT is robust to irrelevant predictor variables and outliers and insensitive to monotone transformations of the input data. However, the lack of smoothness of the prediction surface is one limitation of DT. Another drawback of DT is that a slight change in data can result in quite different splits of the DT, leading to a potentially large variance of prediction. Ridge regression is also a popular learning method, being robust against the multicollinearity of the predictors. However, it can only capture the linear relationship between the response and predictors. There is no single learning method that will dominate all others in all scenarios. In this article, we propose a new ensemble learning framework to combine the strengths of these different learning methods through collaboration.
Ensemble learning has emerged and been extensively studied by many in the past few decades (e.g., dasarathy1979composite, Schapire_1990, ho1995random, and breiman1996bagging), with its popularity recently skyrocketing (e.g., NIPS2017_49ad23d1, yu2018semi, NEURIPS2019_3f7bcd0b, and tian2021rasea). mendes2012ensemble, sagi2018ensemble, and dong2020survey are some recent comprehensive surveys. The general idea of ensemble learning is to combine the predictions obtained from different learning methods (hereafter, base machines), or predictions based on different subsamples or different feature spaces, in order to improve prediction performance. Bagging, stacking, and boosting are three prominent examples. In bagging breiman1996bagging and stacking, base machines are first run in parallel and independently, and then the final prediction is constructed as a simple/weighted average of the predictions from these base machines. In boosting Schapire1998boosting, the base machines work jointly in a top-down manner. In all three algorithms, the output from each base machine is fixed after being calculated. Like human collaboration, an idea that may yield potential improvement is to let different kinds of base machines communicate with each other and update their outputs after observing the predictions of the other base machines. Based on this idea, we propose the Machine Collaboration (MaC) ensemble learning framework with heterogeneous base machines, where the word heterogeneous stands for that the base machines are of different types (e.g. DT, DNN, Ridge Regression).
Compared with bagging, stacking, and boosting, MaC has the following desirable features. Figure (ref) provides the schematic for bagging, stacking, boosting, and MaC. As illustrated, bagging and stacking are parallel & independent, boosting is sequential & top-down, while MaC is circular & recursive. In the framework of MaC, base machines work in a circular manner. Further, the circulation goes multiple rounds. Valuable information is passed recursively through base machines around a “round table", but not top-down. In this process, the base machines update their structures and/or parameters once in each round, according to the information received from the other machines. We demonstrate that MaC can deliver competitive performance when compared with the base machines or other ensemble methods.
Three existing ensemble learning methods are closely related to MaC. One is the so-called super learner (hereafter, SL) proposed in vanderLaanPolleyHubbard2007, which can be viewed as a stacking method. SL receives predictions from different kinds of base machines, then outputs a weighted average of these predictions as the final prediction, where the weights are obtained using cross-validation with some specified loss function. The second one is the LS-Boost proposed in friedman2001greedy. LS-Boost works as a forward regression, using a particular base machine (e.g. tree) each time to fit the current residuals, and ensemble all machines in the end. LS-Boost is then a sequential & top-down method. While both ensemble methods work well in certain situations, we show that our newly proposed MaC performs better in our extensive simulation studies and real data analysis, partly due to its distinctive circular & recursive learning structure. The third one is the so-called additive groves which were proposed to train an additive ensemble of regression trees (See sorokina2007additive and sorokina2008detecting). The main difference is that the additive groves are meant to build an ensemble of homogeneous models (regression trees), more in the spirit of LS-boost or bagging, while MaC focuses on heterogeneous ensembles in the spirit of stacking.
The main contributions of this work are fourfold. First, we propose a new type of ensemble learning framework, MaC, which is circular & recursive. The circular & recursive aspect could be a potential direction for exploring new methods of ensemble learning. Second, we present some desirable finite statistical properties of MaC. Third, we demonstrate via extensive simulations that MaC performs better than all individual base machines and the ensemble methods SL and LS-Boost. Lastly, in the real data analysis, we compare MaC with the competing methods on $119$ benchmark datasets in the Penn Machine Learning Benchmarks (PMLB) Olson2017PMLB, which demonstrate notable advantages of MaC over competing methods for most datasets.
We now introduce the details of the new ensemble learning framework, machine collaboration (MaC). In MaC, we consider a collection of base machines and allow them to collaborate to improve the prediction performance. The base machines could contain both hyperparameters and non-hyperparameters. Here, a hyperparameter is a parameter whose value controls the structure of the base machine and could be set using domain knowledge or selected by cross-validation or data-splitting. The remaining parameters, whose values are estimated by fitting the base machine with fixed hyperparameters on the training data, are non-hyperparameters. For example, in DNN, the learning rate, number of nodes, number of layers, and the activation function are hyperparameters, while for DT, the maximum depth is a hyperparameter. For DNN, the weights of each layer are non-hyperparameters, and for DT, the split variable and value for each split are non-hyperparameters. Hereafter, for the sake of simplicity, we use “parameter” instead of “non-hyperparameter” unless confusion arises. First, we provide a simple sketch of MaC with two base machines in the next subsection. Then the general MaC framework is provided in subsection (ref).
Suppose we have a regression task and two base machines, $M_{A}$ and $M_{B}$ (e.g., a DT and a DNN). Denote the dataset as $\mathbf{D}\ensuremath{=}\left\{ \text{feature }\ensuremath{X},\text{ target }Y\right\} $. For two real-valued $n_{z}$-dimensional vectors $Z=\left(z_{1},\cdots,z_{n_{z}}\right)$ and $\hat{Z}=\left(\hat{z}_{1},\cdots,\hat{z}_{n_{z}}\right)$, denote the loss function as $L\left(z,\hat{z}\right)$, the empirical risk function $R(Z,\hat{Z})\equiv\frac{1}{n_{z}}\sum_{j=1}^{n_{z}}L\left(z_{j},\hat{z}_{j}\right)$. Let $\hat{Y}_{A}$ and $\hat{Y}_{B}$ denote the currently fitted values from machines $M_{A}$ and $M_{B}$, respectively. Note that the key idea of MaC is to update $\hat{Y}_{A}$ and $\hat{Y}_{B}$ alternatively throughout the collaboration process. Figure (ref) is a schematic of MaC with two base machines with the detailed steps as follows.
Note that we adopt a hold-out method to evaluate different machines. In theory, a cross-validation method (hear after, CV), e.g. $k$-fold CV scheme could be adopted, for our experiments, however, Optuna, the hyperparameter optimization software framework taken by us prevents us to adopt a CV method. Taking the $k$-fold CV as an example, Optuna will not give us the evaluation values for the same hyperparameters for all the $k$ folds so that we can't take the average of the evaluation values across the $k$-fold. Developing a new optimization framework for $k$-fold is an interesting future work.
We can easily extend the idea in the MaC for two machines to a situation with more machines. To describe the general algorithm (Algorithm 1) of MaC with more than $2$ base machines, we need the following setup and notations. Suppose we have an independent and identically distributed sample of size $n$ \[ \mathbb{\mathbf{D}}=\left(D_{1},\cdots,D_{n}\right)=\left(\left(X_{1},Y_{1}\right),\cdots,\left(X_{n},Y_{n}\right)\right) \] generated from the true distribution $P_{0}$, which is an element of a statistical model $\mathbb{P}$. The support of $P_{0}$ is $\mathcal{D}\equiv\mathbb{\mathcal{X}}\times\mathcal{Y}=\left\{ d\in\mathcal{\mathbb{R}}^{l}\times\mathbb{R}|P_{0}\left(d\right)\neq0\right\} $, where $l$ is a positive integer denoting the dimension of $\mathbb{\mathcal{X}}$. MaC is constructed based on $K_{n}$ different base machines $\{m_{k,\lambda_{k},\theta_{k}}\}_{k=1,\cdots,K_{n}}$, where $\lambda_{k}$ and $\theta_{k}$ are the vectors of hyperparameters and parameters of the $k$th machine, respectively. Here, we use $K_{n}$ to allow the number of base machines $K$ to grow alongside sample size $n$, as demonstrated in Section (ref). We use the same loss function and risk function for tuning the parameters and hyperparameters of the base machines and the MaC throughout Algorithm 1. Moreover, assume the tuning and estimation algorithm for each base machine is given. Let $\hat{Y}_{k}^{t}$ and $\hat{Y}_{k}^{v}$ denote the predictions of the training sample and validation sample based on the $k$th machine, respectively. For machine $k=1,\cdots,K_{n}$, define $\hat{Y}_{k}\equiv(\hat{Y}_{k}^{t},\hat{Y}_{k}^{v})$, $\hat{Y}_{-k}\equiv\sum_{j\neq k}^{K_{n}}\hat{Y}_{j}$ and $\hat{Y}^{v}\equiv\sum_{j=1}^{K_{n}}\hat{Y}_{j}^{v}$. Define the index pair of the outer and inner loops when the loop stops as $\left\{ i^{*},k^{*}\right\} $.
When the sample size is not large, if we do not use the full data but only $\mathbf{D}^{t}$ to estimate the non-hyperparameter vector, there could be a large information loss. To avoid this, in our experiments in Section (ref), we follow common practice in machine learning (See Listing 4.1 in p.98 of chollet2017deep), namely to use the full data to estimate the non-hyperparameter vector in step 6 of Algorithm 1. Note that, there is a trade-off between computation cost and the accuracy of prediction. To get a higher accuracy of prediction, you would need bigger $T$ or $K_{n}$, leading to higher computation cost.
Under the setting of subsection (ref), denote the cumulative distribution function of $X_{1}$ as $F_{0}$. Let $B_{n}\in\left\{ 0,1\right\} ^{n}$ be a random binary $n$-vector whose observed value defines a split of the data $\mathbf{D}$ into a training sample $\mathbf{D}^{t}$ and a validation sample $\mathbf{D}^{v}$, with $1$ for validation and $0$ for training. Let $p$ denote the proportion of observations in the validation sample, and $P_{n}$, $P$$_{n,B_{n}}^{t}$ and $P_{n,B_{n}}^{v}$ denote the empirical distributions of $\mathbf{D}$, $\mathbf{D}^{t}$, and $\mathbf{D}^{v}$, respectively.
Suppose we have a set of $K_{n}$ base machines $\left\{ M_{1},\cdots,M_{K_{n}}\right\} $. Assume the space of the hyperparameter $\lambda_{j}$ and the parameter vector $\theta_{j}$ of the $j$th base machine are $\Lambda_{j}\subseteq\mathbb{R}^{d_{\lambda,j}}$ and $\Theta_{j}\subseteq\mathbb{R}^{d_{\theta,j}}$, respectively. Then, each base machine $M_{j}:\mathbb{P}\rightarrow\mathcal{S}_{j}\left(\mathcal{X}|\Lambda_{j}\times\Theta_{j}\right)$ is a mapping from $\mathbb{P}$ into $\mathcal{S}_{j}\left(\mathcal{X}\right)\equiv\mathcal{S}_{j}\left(\mathcal{X}|\Lambda_{j}\times\Theta_{j}\right)$. $\mathcal{S}_{j}\left(\mathcal{X}|\Lambda_{j}\times\Theta_{j}\right)$ is a space of real-valued parametric functions from $\mathcal{X}$ to $\mathbb{R}$, taking the vectors of the hyperparameters and parameters in the space $\Lambda_{j}\times\Theta_{j}$. For $j=1,\cdots,K_{n}$, denote the realization of $M_{j}$ as $m_{j,\lambda_{j},\theta_{j}}:\mathcal{X}\rightarrow\mathbb{R}$, which is a function with hyperparameter vector $\lambda_{j}\in\Lambda_{j}$ and parameter vector $\theta_{j}\in\Theta_{j}$.
We consider candidate MaCs constructed by the sum of $K_{n}$ base machines with certain hyperparameters and parameters. In particular, we have $\mathcal{M}\equiv M_{1}(\cdot)+\cdots+M_{K_{n}}(\cdot):\mathbb{P}\rightarrow\mathcal{S}(\mathcal{X})\equiv\mathcal{S}\left(\mathcal{X}|\Lambda(n)\times\Theta(n)\right)$, where $\Lambda(n)=\prod_{j=1}^{K_{n}}\Lambda_{j}$ and $\Theta(n)=\prod_{k=1}^{K_{n}}\Theta_{k}$ denotes the space of the hyperparameter vector $\lambda\equiv\left\{ \lambda_{1},\lambda_{2}\cdots,\lambda_{K_{n}}\right\} $ and parameter vector $\theta\equiv\left\{ \theta_{1},\theta_{2},\cdots,\theta_{K_{n}}\right\} $ of MaC, respectively.
Define the space of all candidate MaCs as $\mathbb{M}\equiv\left\{ \mathcal{M}\left(P\right):P\in\mathbb{P}\right\} \subseteq\mathcal{S}(\mathcal{X})$, then each realization of $\mathcal{M}$, $\mathfrak{m}\equiv\mathfrak{m}_{\lambda,\theta}:\mathcal{X}\rightarrow\mathbb{R}$ is an element of $\mathbb{M}$ with hyperparameter vector $\lambda\in\Lambda(n)$ and parameter vector $\theta\in\Theta(n)$. Denote the loss function of $m$ by $L\left(D,m\right)$ and the risk of $m$ by $E_{P_{0}}L\left(D,m\right)$ for $D\in\mathcal{D}$, and $m$ denoting a base machine or a MaC. Endow the space $\mathcal{S}(\mathcal{X})$ a dissimilarity function $\delta:\mathcal{S}(\mathcal{X})\times\mathcal{S}(\mathcal{X})\rightarrow\mathbb{R}$, with the dissimilarity between $\mathfrak{m}_{1}$ and $\mathfrak{m}_{2}$ in $\mathbb{M}$ defined as \[ \delta\left(\mathfrak{m}_{1},\mathfrak{m}_{2}\right)\equiv\int\left|L(d,\mathfrak{m}_{1})-L\left(d,\mathfrak{m}_{2}\right)\right|\mathrm{d}P_{0}(d). \] Define the pseudo-true MaC as
and the risk difference of $\mathfrak{m}$ as $\delta\left(\mathfrak{m},\mathfrak{m}_{0}\right)$.
Denote the space of the hyperparameter vector in which the MaC algorithm searches by $\tilde{\Lambda}(n)\subseteq\Lambda(n)$ with cardinality $\left|\tilde{\Lambda}(n)\right|=J_{n}$. Indexing the elements of $\tilde{\Lambda}(n)$ as $\left\{ \lambda_{\left(1\right)},\cdots,\lambda_{\left(J_{n}\right)}\right\} $, we construct $J_{n}$ subspace $\mathbb{M}_{k}\equiv\left\{ \mathfrak{m}_{\lambda_{\left(k\right)},\theta}:\theta\in\Theta(n)\right\} \subseteq\mathbb{M}$ according to the selection of hyperparameter vector for $k\in\kappa_{n}\equiv\left\{ 1,\cdots,J_{n}\right\} $. For a particular subspace $\mathbb{M}_{k}$ with fixed hyperparameter vector $\lambda_{\left(k\right)}$, the MaC algorithm searches for the optimal parameter vector on the parameter space $\tilde{\Theta}_{k}(n)\subseteq\Theta(n)$. After collecting MaCs whose hyperparameter vector is $\lambda_{\left(k\right)}$ and parameter vector $\theta\in\tilde{\Theta}_{k}(n)$, we create the space $\mathbb{M}_{k,\tilde{\Theta}_{k}}\equiv\left\{ \mathfrak{m}_{k,\theta}\equiv\mathfrak{m}_{\lambda_{\left(k\right)},\theta}:\theta\in\tilde{\Theta}(n)\right\} \subseteq\mathbb{M}_{k}$. Define the risk approximation error of $\mathbb{M}_{k,\tilde{\Theta}_{k}}$ as $B_{0}\left(k\right)\equiv\min_{\mathfrak{m}\in\mathbb{M}_{k,\tilde{\Theta}_{k}}}\delta\left(\mathfrak{m},\mathfrak{m}_{0}\right)$. Note that, in van der Laan (2006), $\mathbb{M}_{k,\tilde{\Theta}_{k}}$ is called epsilon-net, which is a form of sieving net for the parameters.
For any empirical distribution $P_{n}$, define the estimated model with fixed $k$ based on our algorithm as \[ \hat{\mathcal{M}}_{k}\left(P_{n}\right)\equiv\underset{\mathfrak{m}\in\mathbb{M}_{k,\tilde{\Theta}_{k}(n)}}{\arg\!\min}\int L(d,\mathfrak{m})\mathrm{d}P_{n}(d). \] Then, the estimated MaC is $\hat{\mathcal{M}}\left(P_{n}\right)\equiv\hat{\mathcal{M}}_{k\left(P_{n}\right)}\left(P_{n}\right),$ where $k\left(P_{n}\right)\equiv\arg\!\min_{k\in\kappa_{n}}\int L(d,\hat{\mathcal{M}}_{k}\left(P_{n}^{t}\right))\mathrm{d}P_{n}^{v}(d)$.
To introduce the following theorem, we need the following definition. \theoremstyle{definition}
Note that in the definition above, the space $\tilde{\Theta}_{k}(n)$ is first fixed. If we first fix $\varepsilon_{k}$, and then select a space $\tilde{\Theta}_{k}(n)$ with minimum $N_{k}$ so that $\cup_{j=1}^{N_{k}}B\left(\mathfrak{m}_{\left(k,j\right)},\varepsilon\right)$ can cover $\mathbb{M}_{k}$, then $N_{k}$ is just what we would call a covering number.
Let us refer to the value of the parameter vector that results in the minimum of the empirical loss of estimation as the optimal parameter vector. {In the theoretical analysis, following vanderLaanPolleyHubbard2007, we only consider a MaC algorithm that searches the optimal parameter vector on a discrete set.\footnote{{[}}1{]}{{As $N_{k}$ goes to infinity with $n$, Theorem (ref) can be regarded as an approximate result for the MaC algorithm that searches the optimal parameter vector on the whole domain.}} For $\tilde{\Theta}_{k}(n)=\left\{ \theta_{\left(k,1\right)},\cdots,\theta_{\left(k,N_{k}\right)}\right\} $ with $N_{k}<\infty$, we have the following finite sample result.}
The proof of Theorem (ref) is provided in the supplement. Theorem (ref) shows that the expectation of the risk difference of $\hat{\mathcal{M}}\left(P_{n,B_{n}}^{t}\right)(x)$ and $\mathfrak{m}_{0}(x)$ depends on the searching number $N_{k}$, and the cardinality $J_{n}$ of the searching space of the hyperparameters $\tilde{\Lambda}(n)$. Suppose we search for the optimal MaC on $\mathbb{M}_{k}$. If the searching number $N_{k}$ of $\mathbb{M}_{k}$ is small, the searching resolution $\varepsilon_{k}\left(N_{k}\right)$ may be large, and then $\tilde{B}_{0}\left(k\right)$ could be large too. As a result, a small $N_{k}$ could result in the risk bound increasing. Moreover, if $J_{n}$ is small, the value of $\min_{k\in\kappa_{n}}\tilde{B}_{0}\left(k\right)$ could be large. To reduce the prediction risk, we need to not only adjust $N_{k}$ to strike a balance between the term associated with $\log\left(N_{k}\right)$ and that associated with risk in the approximation error $\tilde{B}_{0}\left(k\right)$ but also adjust the balance between $\min_{k\in\kappa_{n}}\tilde{B}_{0}\left(k\right)$ and the term with $\log\left(J_{n}\right)$. Note that increasing the number of base machines, $K_{n}$, may reduce the first term of eq.((ref)), but cause the term with $\log\left(J_{n}\right)$ to increase, hence a proper selection of the base machines and a moderate $K_{n}$ would also be helpful to reduce the risk.
Theorem (ref) also illustrates a merit of the circular & recursive feature of MaC, that is, MaC only needs to add up the final $K_{n}$ machines during the loop of the algorithm without the need to record all single machines in the loop. Note that $K_{n}$ is ordinarily much smaller than the total number of all single machines of the algorithm. Using a circular & recursive type algorithm can reduce $J_{n}$ compared with a sequential & top-down type algorithm because a \textit{sequential }&\textit{ top-down} type algorithm needs to add all single machines into its final estimate.
We first generate data for $10$ independent variables, $x_{1},x_{2},\cdots,x_{10}$, following a multivariate normal distribution $N(0,\Sigma)$, with a variance--covariance matrix $\Sigma=[\Sigma_{ij}]_{10\times10}$, where $\Sigma_{ij}=\rho^{|i-j|}$ with $\rho=0.1$. Then, generate the error term $\varepsilon\sim N\left(0,1\right)$. We consider the following two data-generating processes (DGPs).
where $I(\cdot)$ stands for indicator function and the constants $c_{0}$ through $c_{5}$ are chosen to standardize each term (mean 0 and variance 1). It is clear to see that DGP 2 has stronger nonlinearity than DGP 1. Note that the variables $x_{6},\cdots,x_{10}$ are not included in either of the DGPs, but are included as predictors for all methods. The sample size $n$ and the number of replications for both is $1000$. We set $p=0.25$. The sample is randomly split into training data of $600$ observations, validation data of $200$ observations, and test data of $200$ observations. We choose the mean squared prediction error $\textrm{MSPE}_{i}=\frac{1}{n}\sum_{j=1}^{n}\left(\hat{y}_{i,j}-y_{j}\right)^{2}$ as the loss for both the artificial and real data experiments, where $\hat{y}_{i,j}$ is the prediction of $y_{j}$ for the $i$th replication.
We use DT (CART, breiman1984classification), DNN, and Ridge regression as base machines. The DNN used here has $5$ dense layers. Each of the first $4$ dense layers is associated with a dropout layer. We treat the number of nodes of the latent layers, the types of activation functions, the dropout ratio, the learning rate of optimization, the batch size, and the number of epochs as hyperparameters. The maximum depth of the DT is fixed as $10$. The DT is pruned according to the cost complexity parameter, which is a hyperparameter. The tuning parameter that controls the strength of the penalty of Ridge is also a hyperparameter. The order of the base machines is also treated as a hyperparameter. As mentioned in Section (ref), all the hyperparameters are tuned with respect to the MSPE in the validation data. Using these three base machines for MaC, LS-boost, and SL, we compare the prediction performances on the test data between MaC, SL, and LS-Boost for both DGPs. For MaC, we set $T=50$ and $\tau=10$.
We conduct LS-boost with a similar setup of MaC. The three base machines are repeatedly used with a fixed order for at most $T=50$ times. The order is selected according to the empirical risk of the validation data. The processes will be stopped if the empirical risk of the validation data does not decrease for $\tau=10$ times. For LS-boost and SL, the hyperparameters of each base machine are tuned according to the empirical risk of the validation data, when we fit each base machine.
We use the Greene HPC Cluster of New York University to carry out all the experiments. All the experiments are carried out by CPU with a 2x Intel Xeon Platinum 8268 24C 205W 2.9GHz Processor. For each simulation replication, the computing time of the artificial simulation is about $64$ minutes, while that of the real data experiments with datasets of different sizes ranges from $42$ minutes to $50$ hours.
The simulation results for DGP 1 are in Figures (ref) and (ref), with that of DGP 2 in Figures (ref) and (ref). We count the number of replications for which MaC wins, i.e., it produces a smaller MSPE for the test data than its competitors, namely a particular ensemble method or a base machine. As depicted in Figures (ref) and (ref), for both DGPs, MaC generally outperforms all the other methods. The boxplots in Figure (ref) and (ref) show that MaC has a smaller mean and median of prediction errors on the test data than all the other methods, whose means and medians are marked by the blue rectangles and the orange lines, respectively. For the simulation and the real data experiment (in the next subsection), we calculate the mean and median of MSPE for each method, the paired t-statistic, and Cohen's $d$ cohen1988statistical of the means of MSPEs for each pair, which consists of an alternative method and MaC. The results of the simulations are tabulated in the first two panels of Table (ref). For both DGPs, all the paired t-statistics are large with values of at least $4.6$, and most of Cohen's $d$s are greater than $0.2$, implying that the mean and median of MSPE for the MaC are significantly smaller than those of the other methods.
We use the same base machines and the same setting as the simulations for our real data experiment. The real datasets are from PMLB, which is an open-source dataset collection for benchmarking machine learning methods. All the datasets do not contain personally identifiable information or offensive content. There are 122 datasets for regression in PMLB in total, including data about automobile prices, faculty salaries, pollution, and crime. The sample sizes of the datasets range from $47$ to $1025010$. We dropped three datasets with sample sizes greater than or equal to $1$ million because of our limits in computing resources. As a result, the datasets used have sample sizes ranging from $47$ to $177147$ with a mean of $5476.69$, while the numbers of features are from $2$ to $1000$ with a mean of $26.05$. To make the comparison between the different methods simpler, we standardize all variables in the datasets so they all have a mean of $0$ and a variance of $1$. We apply MaC, SL, LS-Boost, and the three individual base machines to predict
the target variable in the test data for each dataset. We perform the experiment $20$ times, in which for each time, each dataset is randomly split into training ($64\%$), validation ($16\%$), and test ($20\%$). All the prediction results are in Figures (ref), (ref), and (ref). We count the number of datasets for which MaC outperforms its competitors in terms of MSPE on the test data. The results are in Figure (ref), which shows that MaC wins for more than $60\%$ datasets against all competitors. Moreover, as illustrated in Figure (ref), MaC has a smaller mean and median of MSPE than all the competing methods. To obtain a closer look at the settings where MaC has an advantage, we calculate the MSPE differences by subtracting the MSPEs of the other methods from that of MaC. The results of the MSPE differences are plotted in Figure (ref) for all the methods. As shown, for each subfigure, there are more points with positive MSPE difference than with negative MSPE difference. It means that MaC has smaller MSPE than corresponding alternative methods. In Figure (ref), we draw an orange cross when MaC outperforms others and a blue cross otherwise. MaC has superior performance across the entire range of sample sizes and the number of variables. Similar to the simulation, the paired t-statistics and Cohen's $d$ for the real data in the third panel of Table (ref) show that the mean and median of MSPE of the MaC are significantly smaller than those of the other methods. \footnote{{[}}2{]}{Note that the MSPE of SL for one dataset is extremely large (42.48), which we ignore in some figures for better visibility (we calculate the paired t-statistic and Cohen's d after deleting the pair corresponding to this dataset).}
In this paper, we propose a new ensemble learning framework, MaC, for regression problems. The key feature of MaC being circular & recursive helps it to communicate among the base machines, yielding better performance in various scenarios. However, the framework of making an ensemble algorithm circular & recursive is not limited to supervised regression learning tasks. Some interesting extensions would be to apply this same MaC framework to other types of tasks such as classification and semisupervised learning. In our theoretical analysis, we derive the risk bound of MaC using a quadratic risk function. Deriving risk bounds for more general risk functions and comparing the risks of MaC with individual base machines, super learner, and boosting could also lead to a deeper understanding of MaC. We intend to address these extensions in future work. Similar to boosting methods, MaC has a limitation on the burden of computation, since the current version of the algorithm of MaC can't be executed in parallel. It is another important future work to reduce the burden of computation.
Partially supported by NIH Grant 1R21AG074205-01 and NSF Grant DMS-2324489 (Feng) and JSPS KAKENHI (grant number JP19K01582 and JP22H00833) (Liu).