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.
72,241 characters · 12 sections · 42 citation commands
On the Subbagging Estimation for Massive Data
\def\spacingset#1{ {#1}} \spacingset{1}
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem} \newtheorem{proposition}{Proposition} \newtheorem{corollary}{Corollary} \newtheorem{assumption}{Assumption} \newtheorem{lemma}{Lemma} \newtheorem{example}{Example}
\if00 {\spacingset{1.5}
} \fi
\if10 {\spacingset{2}
} \fi
{\bf \sl Keywords:} Big Data; Incomplete $U$-Statistics; Memory Constraint; {\color{black} Parameter Estimation; Subbagging}.
\spacingset{1.45}
{\color{black}With the rapid development of information and technology, more and more datasets with massive data sizes are unprecedentedly available for econometric and statistical analysis. Such enormous sizes of data bring in both a blessing and a curse. On one hand, the ultra large data size can lead to high estimation efficiency, which means much more accurate estimation output. On the other hand, analyzing the immense size of data is possibly beyond the limits of our computers. For example, the massive datasets are sometimes too large to be read into a computer's memory. A personal computer (PC) typically encounters this bottleneck, while a supercomputer with petabytes of data storage also has limited node memory in the process of computation (e.g., the supercomputer Gadi at \url{https://opus.nci.org.au/display/Help/Queue+Limits}).
To alleviate the bottleneck of memory constraints, one possible approach is to downsize the data volumes by constructing a subset of data randomly, or in other words, drawing a subsample, instead of analyzing the whole dataset. Though the economic wisdom of diminishing returns also suggests that we may not need the entire large dataset to estimate parameters or perform hypothesis testing, constructing subsamples randomly for massive data analysis have not been sufficiently studied in the econometrics literature except for a brief exposition by ng2017opportunities, a seminal review article by lee2020econometric, and lee2020sketching about two-stage least squares estimation for big data, where these references refer to this approach as the term “sketching”. Analyzing a sketch/subsample of big data has also drawn the attention of statisticians and machine learning researchers in recent years; see, e.g., dhillon2013new, wang2018optimal and yu2020optimal. A more complete review can be found in lee2020econometric.
Except for the memory constraint, lee2020econometric also provides other motivations for working with sketches/subsamples of big data, and we summarise those motivations below. First, the massive data are of limited practical use if they are too expensive to store. In this case, a researcher may need to get a sense from a small subsample of whether the expensive full sample analysis would be worth pursuing. Second, even if a computer's memory is enough to cache the big dataset, drawing and loading subsamples are much faster than loading the whole big dataset into the memory. A subsample can be adequate when a researcher is still attempting to specify an appropriate model. Third, the computation for complex estimating problems based on a subsample can be much easier, and debugging in the process is much faster as well. Last, it may only be allowed to circulate a subsample of dataset in real practice due to the reasons of confidentiality.
Though sketching/subsampling can alleviate the computer bottlenecks, the statistical view for efficient estimation and inference inevitably favors using as many observations as possible. Specifically, the estimation based on one subsample that can satisfy the memory constraint, is far less efficient than the full sample estimation. In order to improve the efficiency, the aggregation methods (e.g., bagging), have been used as a general-purpose procedure for reducing the variance in statistical learning hastie2009elements. Similar ideas can be adopted via aggregating more than one subsample, which leads to the subsample aggregating or subbagging approaches (suggested by, e.g., andonova2002simple, buhlmann2003bagging and zaman2009effect, for other machine learning problems). lee2020econometric refers to this approach as the term “combining sketches” and they find that pooling over subsamples always provides more efficient estimates. By the principle of divide and conquer, drawing several subsamples and computing their estimates can be realized in parallel, and hence are still computationally inexpensive. Instead of efficiency gains, combining several subsamples which are drawn uniformly at random can help avoid the possible sample bias problem caused by drawing only one subsample.
Different from the existing literature, the focus of this article is on establishing the asymptotic theory for the subbagging estimator/combining sketches, such that the efficiency gain from the subbagging approach can be clearly revealed in the asymptotic variance of the estimator. We specifically target on a general estimation approach of solving estimating equations. This general estimation framework comprises a wide class of estimators, including but not limited to the maximum likelihood, least squares, method of moments and generalized estimating equation estimators as special cases. Since econometric models usually specify a set of moment conditions for the true parameter, the approach of solving estimating equations is also widely used in the econometrics literature (see, e.g., rilstone1996second,kim2016higher).
Specifically, for the whole dataset with massive data size $N$, we propose the subbagging estimator which is the aggregation of the subsample estimators computed based on $m_N$ subsamples, respectively. The size of each subsample $k_N$ is set to be far less than $N$ to meet the memory constraint. Based on the subbagging literature (see, e.g., buhlmann2003bagging and mentch2016quantifying), the $m_N$ subsamples are independently and identically drawn from $S_N^{k_N}$ with equal probability $1/{N\choose k_N}$, where $S_N^{k_N}$ is a set consisting of all the $N\choose k_N$ subsamples of data with size $k_N$. This sampling procedure allows us to directly draw subsamples from the hard drive, which is important in the sense that it is almost impossible to load the terabytes or petabytes of full data into memory for random sampling. Specifically, drawing each subsample from $S_N^{k_N}$ with equal probability $1/{N\choose k_N}$ can be accomplished via the simple random sampling without replacement, which only requires $O(k_N)$ of memory and $O(k_N\log_2 k_N)$ of sampling time (see, e.g., Gupta1984algorithm). Hence, drawing $m_N$ subsamples independently and identically from the hard drive is feasible and computationally efficient in real practice.
Based on the sampling procedure, the major challenge that we encounter in establishing the asymptotic theory for the subbagging estimator is that any one subsample (one subset of data) from the drawn $m_N$ subsamples is possibly overlapping with the other subsample (the other subset). The standard independent asymptotics cannot deal with the aggregation of subsamples that possibly have overlap with each other. To overcome this challenge, we extend the theory of the incomplete $U$-statistics with an infinite order kernel, which was originally developed for the random forests (see, e.g., mentch2016quantifying,peng2019asymptotic), and then we use it to show the theoretical properties of the subbagging estimator.
This novel asymptotic framework leads to three important findings for the subbagging estimation approach. The first is that subbagging can indeed reduce the asymptotic variance of estimators as the number of drawn subsamples $m_N$ increases. In particular, the subbagging estimator can achieve $\sqrt{N}$-consistency and asymptotic normality, but has an inflation rate of $1/\alpha$ in its asymptotic variance compared to the full sample estimator, where $\alpha=\lim (k_Nm_N/N)\in (0,\infty]$. The second is that the subbagging approach cannot improve the order of the asymptotic bias, which relies on the subsample size $k_N$. As a consequence, in order to achieve $\sqrt{N}$-consistency for the subbagging estimation, the asymptotic bias needs to be taken into the consideration. To this end, we propose two subbagging algorithms with proper hyperparameter selection of $k_N$ and $m_N$, to balance the computation cost and achieve $\sqrt{N}$-consistency at the same time. The third is that a subbagging variance estimation approach is proposed in this article for inference of parameters under the big data setting. This can be treated as extensions of the $m$-out-of-$n$ bootstrap (see, e.g., bickel1988richardson,bickel2012resampling), or the subsampling methods developed for frequentist inference (see, e.g., hong2006fast).
The rest of this paper is organized as follows. Section (ref) introduces the subbagging estimator based on a simple average ensemble and develops its theoretical properties. Section (ref) proposes bias correction before the simple average ensemble and demonstrates its advantage. Section (ref) provides a consistent subbagging variance estimation approach for inference of parameters. Simulation studies and an application of the American airline dataset are elaborated in Section (ref) to illustrate the usefulness of the proposed subbagging methods. Section (ref) concludes the paper with a discussion. All theoretical proofs are relegated to Appendix. Their extensive derivations are presented in the supplementary material.
}
{\color{black}
Consider that we have observations of $p$-dimensional data vectors $Z_1,\cdots,Z_N\stackrel{iid}\sim Z$ and $Z$ follows a distribution $F_{Z}$. Let $\theta_0$ be the $d$-dimensional true parameter vector of interest. The aim of this article is to consider the estimation and inference problems for parameter $\theta_0$ under the massive data setting, i.e., the sample size $N$ is extremely large. In the econometrics and statistics literature, various estimation approaches have been developed to estimate the parameter $\theta_0$ based on data $Z_1,\cdots,Z_N$. Without loss of generality, we can denote any one of those estimators by $\hat\theta_N(Z_1,\cdots,Z_N)$, which is a function of the full sample $Z_1,\cdots,Z_N$.
The ultra large data size $N$ can lead to high estimation efficiency of the full sample estimator $\hat\theta_N(Z_1,\cdots,Z_N)$, for instance, the usual $\sqrt{N}$-consistency for the estimator. However, computing $\hat\theta_N(Z_1,\cdots,Z_N)$ requires loading the whole big dataset in the memory, which is not plausible under the memory constraint. In order to satisfy the memory constraint of a certain computer, we possibly can only analyze a smaller size $k_N$ of data, where $k_N< N$. Let $s=\{s_1,\cdots,s_{k_N}\}\subset\{1,\cdots, N\}$ be a subsample (subset) of data $\{1,\cdots,N\}$. In total, we have $N \choose k_N$ subsamples (subsets). We denote $S_{N}^{k_N}$ as a set consisting of all the subsamples with size $k_N$, so that each $s\in S_N^{k_N}$. Therefore, the size of $S_{N}^{k_N}$, denoted by $|S_{N}^{k_N}|$, is equal to $N \choose k_N$. A subsample estimator can be obtained by $\hat\theta_{k_N}(Z_{s_1},\cdots,Z_{s_{k_N}})$, where the estimation function $\hat\theta_{k_N}()$ is the subsample version of $\hat\theta_N()$. If the full sample estimator $\hat\theta_N(Z_1,\cdots,Z_N)$ is $\sqrt{N}$-consistent, then the subsample estimator is $\sqrt{k_N}$-consistent.
Clearly, the rate of the subsample estimator reduces to $\sqrt{k_N}$. Accordingly, even though the subsample estimator requires less memory to compute, the bias and variance of the estimator both increase. In order to reduce the variance, we consider the approach of bootstrap aggregating, or bagging, which is a general-purpose procedure for reducing the variance in statistical learning hastie2009elements. Here, we adopt similar ideas and consider the subsample aggregating, which can be called subbagging. Ideally we can take all $N\choose {k_N}$ subsamples in $S_{N}^{k_N}$, and build a subbagging estimator
by averaging the $N\choose {k_N}$ subsample estimators.
Of course, building $N\choose {k_N}$ subsample estimators is computationally infeasible even for moderate size $N$, and a substantial improvement in computation efficiency can be made by building and taking average over only $m_N<{N\choose {k_N}}$ subsample estimators. In this case, the subbagging estimator becomes
where $|\mathcal{S}|=m_N$. Based on the subbagging literature (see, e.g., buhlmann2003bagging and mentch2016quantifying), we consider that these $m_N$ subsamples in $\mathcal{S}$ are independently and identically drawn from $S_N^{k_N}$ with equal probability $1/{N\choose k_N}$. In the following sections, we develop a novel theoretical framework to derive the theoretical properties of the subbagging estimator, and explicitly show how subbagging improves the efficiency of the subsample estimator.
It is worth noting that in the above subbaging approach, one subsample $s^{(1)}=\{s_1^{(1)},\cdots,$ $s_{k_N}^{(1)}\}\in\mathcal{S}$ is possibly overlapping with the other subsample $s^{(2)}=\{s_1^{(2)},\cdots,s_{k_N}^{(2)}\}\in\mathcal{S}$. It is also even possible that $s^{(2)}=s^{(1)}$ due to the sampling strategy of ${\mathcal{S}}$. Therefore, the standard independent asymptotic techniques cannot be directly applied for the statistic ((ref)).
In fact, the subbagging estimators ((ref)) and ((ref)), are a complete $U$-statistic and an incomplete $U$-statistic, respectively, with the kernel function being $\hat\theta_{k_N}()$ and the order of kernel being $k_N$. When $k_N$ is fixed, the resulting incomplete $U$-statistic remains asymptotically normal; see janson1984asymptotic or lee2019u. However, we require $k_N\to\infty$ to guarantee $\sqrt{k_N}$-consistency for the convergence of each subsample estimator. frees1989infinite studied the complete $U$-statistic with infinite $k_N$; mentch2016quantifying introduced a central limit theorem (CLT) for the incomplete $U$-statistics with the order of kernel $k_N\to\infty$; peng2019asymptotic further extended the CLT. Nevertheless, mentch2016quantifying focused more on solving the inference problem for the random forests, and hence their conditions in the CLT, including the conditions for $k_N$ and the kernel function, can be improved for the subbagging estimation problem given here. For instance, they require the subsample size $k_N=o(N^{1/2})$ which will cause a large bias for the subsample estimator. peng2019asymptotic studied the CLT for the random forests under condition $k_N=o(N)$ and some restrictions imposed on the kernel function, but the CLT was established based on a different resampling setting compared to this paper. To this end, we show an improved version of CLT for the general incomplete $U$-statistics with infinite $k_N$ and more general conditions on the kernel function in Lemma (ref) of Appendix, where its extensive proof is relegated to the supplementary material. Using Lemma (ref), we obtain the theoretical properties of the subbagging estimator in the following subsection.
Throughout the following sections, we denote the subsample and full sample estimators by $\hat\theta_{k_N,s}=\hat\theta_{k_N}(Z_{s_1},\cdots,Z_{s_{k_N}})$ and $\hat\theta_N=\hat\theta_N(Z_1,$ $\cdots,Z_N)$, respectively, to simplify the notation. In this section, we target on two general estimators, the $Z$-estimator and $M$-estimator (named by van2000asymptotic), which are widely used in the literature. Without loss of generality, we introduce the $Z$-estimator or the $M$-estimator for a subsample $s=\{s_1,\cdots,s_{k_N}\}$. The full sample estimator can be naturally obtained by letting $k_N=N$ and $s=\{1,\cdots, N\}$.
Specifically, the subsample $Z$-estimator $\hat\theta_{k_N,s}$ satisfies a system of estimating equations
with respect to $\theta$, where $\psi_\theta(z)$ is a $d$-dimensional measurable function of $z$. In contrast, the subsample $M$-estimator can be obtained by minimizing an objective function $ \sum_{i=1}^{k_N} \mathcal{M}_{\theta}(Z_{s_{i}}) $ with respect to $\theta$, where $\mathcal{M}_\theta(z)$ is a $1$-dimensional measurable function of $z$. To unify these two estimation approaches, if we consider $\mathcal{M}_\theta(z)$ to be a convex smooth function with respective to $\theta$, and denote its partial derivative $\partial \mathcal{M}_\theta(z)/\partial \theta$ by $\psi_{\theta}(Z_{s_{i}})$, then the $M$-estimator can also be obtained via solving estimating equations ((ref)). Accordingly, we consider that the subsample estimator $\hat\theta_{k_N,s}$ solves the estimating equations ((ref)) in the rest of this paper. It is worth noting that this estimator comprises the maximum likelihood, least squares, method of moments and generalized estimating equation estimators as special cases, and is widely used in the econometrics and statistics literature (see, e.g., rilstone1996second,kim2016higher), since econometric models usually specify a set of moment conditions ${\rm E} \psi_{\theta_0}(Z_i)=0$ for the true parameter $\theta_0$, which is the population version of ((ref)).
Based on the subsample estimator $\hat\theta_{k_N,s}$, the subbagging estimator $\hat \theta_{k_N,m_N}$ can be obtained via the subsample aggregation in ((ref)). In this subsection, we make use of the CLT for the general incomplete $U$-statistics with an infinite order kernel in Lemma (ref) of Appendix to derive the asymptotic property of a general subbagging estimator $\hat \theta_{k_N,m_N}$. For the illustration purpose, we first provide the asymptotic theory for the subbagging mean estimator in the following example.
Let $\|\cdot\|_2$ be the Euclidean norm for any generic vector. We then obtain the following theorem for the subbagging mean estimator.
Compared to the variance of the subsample estimator ${\rm Var}(\hat\mu_{k_N,s})=\Sigma_0/k_N$, Theorem (ref) shows that the variance of the subbagging estimator is indeed reduced for all scenarios (1) -- (3). In contrast to the variance of the full sample estimator ${\rm Var}(\hat\mu_{N})=\Sigma_0/N$, the variance inflation $ \Sigma_{0}\{ 1-{1}/{{N\choose k_N}}\}/(k_Nm_N)$ in ${\rm Var}( \mu_{k_N,{m_N}})$ is due to the extra randomness in sampling $s$ from $S_N^{k_N}$. However, this variance inflation can tend to be smaller than $O(N^{-1})$ as $m_N$ gets larger, which means the subbagging estimator can achieve $\sqrt{N}$-consistency for large $m_N$ as shown in scenarios (1) and (3).
Based on Theorem (ref), in order to result in the asymptotic normality for the subbagging mean estimator, we require the condition $k_N/N\to 0$. If $m_N$ is selected too small such that $m_Nk_N/N\to \alpha\in[0,\infty)$, then we require a further moment condition ${\rm E}\|Z\|_2^4<\infty$ to obtain the asymptotic normality. Similar theoretical results can be established for a general subbagging estimator $\hat\theta_{k_N,m_N}$ using the CLT for the incomplete $U$-statistics in Lemma (ref) of Appendix. We summarize the asymptotic property of a general subbagging estimator $\hat\theta_{k_N,m_N}$ based on estimating equations ((ref)) in the theorem below.
To present the theorem, let $ \Sigma_{\theta_0}={\rm Var}\big( \psi_{\theta_0}(Z)\big)$ be the variance-covariance matrix of $\psi_{\theta_0}(Z)$, $V_{\theta_0}={\rm E}\{{\partial \psi_{\theta_0}(Z)}/{\partial \theta^\top}\}$, $V^{(2)}_{\theta_0}={\rm E}\{{\partial^2 \psi_{\theta_0}(Z)}/{(\partial \theta^\top\otimes \partial \theta^\top)}\}$, $H_{\theta_0}={\rm E}\big[ \psi_{\theta_0}(Z)^\top \otimes \{{\partial \psi_{\theta_0}(Z) }/{\partial \theta^\top}\} \big]$, $J_{k_N,s}=k_N^{-1/2} \sum_{i=1}^{k_N}\psi_{\theta_0}(Z_{s_i})$, $J^{(1)}_{k_N,s}=k_N^{-1/2}\sum_{i=1}^{k_N}\{\partial \psi_{\theta_0}(Z_{s_i})/$ $\partial \theta^\top-V_{\theta_0}\}$ and
where $\otimes$ denotes the Kronecker product, and if we denote $(i,j)$-th entry of ${\partial \psi_{\theta}(Z)}/{\partial \theta^\top}$ by $\big({\partial \psi_{\theta}(Z)}/{\partial \theta^\top}\big)_{ij}$, then we define \[ \frac{\partial^2 \psi_{\theta}(Z)}{\partial \theta^\top\otimes \partial \theta^\top}=\left( \frac{\big({\partial \psi_{\theta}(Z)}/{\partial \theta^\top}\big)_{ij}}{\partial \theta^\top} \right). \] Also denote ${\rm vec}(\cdot)$ as the vectorization for any generic matrix.
It is worth noting that the asymptotic normality is from the leading term $-k_N^{-1/2}V_{\theta_0}^{-1}{m_N^{-1}}$ $\sum_{s\in{\mathcal{S}}}{J_{k_N,s}}$ in the asymptotic expansion, which is an incomplete $U$-statistic. Consequently, the CLT for the incomplete $U$-statistics in Lemma (ref) of Appendix is used to establish the asymptotic normality in this theorem; see details of the proof in Appendix (ref). Based on setting (i) in this theorem, we find that $1/\alpha$ is the variance inflation rate for the subbagging estimator $\hat\theta_{k_N,m_N}$ compared to the full sample {\color{black}estimator} $\hat\theta_{N}$. If $\alpha=\infty$, define $1/\alpha=0$. In this case, there is no variance inflation for the subbagging estimator, and its asymptotic variance $N^{-1}V_{\theta_0}^{-1} \Sigma_{\theta_0}( V_{\theta_0}^{-1})^\top$ is exactly equal to that of the full sample estimator $\hat\theta_N$.
Different from Example (ref) and Theorem (ref), the general subbagging estimator $\hat\theta_{k_N,m_N}$ in Theorem (ref) is a consistent but biased estimator, whose asymptotic bias can be obtained by the mean of $k_N^{-1}{m_N^{-1}}\sum_{s\in{\mathcal{S}}}{\mathcal{B}_{k_N,s}}$ in the asymptotic expansion. Using the results of Theorem (ref), we find that this asymptotic bias is $B_{\theta_0}/k_N=O(1/k_N)$, which is the same as the asymptotic bias of the subsample estimator $\hat\theta_{k_N,s}$; e.g., see Lemma (ref) of Appendix. The reason behind this fact is that solving estimating equations ((ref)) generally results in a biased estimator. In sum, the subbagging approach can improve the efficiency compared to the subsample estimator as shown in Theorem (ref) -- (ref), but cannot improve the bias.
Under setting (i) of Theorem (ref), the asymptotic variance of the subbagging estimator is $N^{-1}(1+1/\alpha)V_{\theta_0}^{-1} \Sigma_{\theta_0}( V_{\theta_0}^{-1})^\top=O(N^{-1})$, and hence the asymptotic mean squared error is MSE = bias$^2$ + variance = $O(k_N^{-2}+N^{-1})$. In order to make the square of the bias much smaller than the variance, we further require $k_N/\sqrt{N}\to \infty$, which leads to MSE = $O(N^{-1})$, namely the $\sqrt{N}$-consistency for the subbagging estimator. This gives an intuitive explanation about the reason why the $\sqrt{N}$-consistency requires a technical condition $k_N/\sqrt{N}\to \infty$ in setting (i). Setting (ii) of Theorem (ref) seems to be computationally most efficient because $m_N$ only requires a small setting $m_N/(N/k_N)\to 0$. However, the asymptotic variance of $\hat\theta_{k_N,m_N}$ is of order $(k_Nm_N)^{-1}$ which is far less than $N^{-1}$ compared to the asymptotic variance of the full sample estimator $\hat\theta_{N}$. Consequently, we recommend setting (i) in real practice to balance the computational feasibility and the estimation efficiency.
To illustrate the usefulness of Theorem (ref), we provide the asymptotic theory for the subbagging covariance matrix estimator in the following example.
Let $K^{(3)}={\rm E}\left[ \left\{ (Z-\mu_0)\otimes (Z-\mu_0)\right\} (Z-\mu_0)^\top\right]$, $K^{(4)}={\rm E}[\{ (Z-\mu_0)(Z-\mu_0)^\top \} \otimes \{(Z-\mu_0)(Z-\mu_0)^\top \}]$, $L_p$ be the $p(p+1)/2\times p^2$ elimination matrix defined in lutkepohl2005new, and $K^{(4)}_L=L_p \{K^{(4)}-{\rm vec}(\Sigma_0)$ ${\rm vec}^\top(\Sigma_0)\} L_p^\top$. Then we have the following corollary for the subbagging covariance matrix estimator.
Corollary (ref) is a direct result of Theorem (ref) based on the setting of Example (ref). It is worth noting that the subbagging covariance matrix estimator is consistent but biased, because the subsample estimator $\hat\Sigma_{k_N,s}=k_N^{-1}\sum_{i=1}^{k_N} (Z_{s_i}-\hat\mu_{k_N,s})(Z_{s_i}-\hat\mu_{k_N,s})^\top$ is a biased estimator. Similar to the general subbagging estimator in Theorem (ref), this biasedness is the reason why we require a technical condition $k_N/\sqrt{N}\to \infty$ to achieve $\sqrt{N}$-consistency in setting (i) of Corollary (ref).
It is worth noting that the memory constraint of a computer can be treated as a constraint on the subsample size $k_N$. In order to both satisfy the memory constraint of the computer and achieve $\sqrt{N}$-consistency for the subbagging estimation, the hyperparameters of $k_N$ and $m_N$ need to be chosen properly based on the theoretical results of Theorems (ref) -- (ref). We summarize the $\sqrt{N}$-consistent subbagging estimation algorithm with hyperparameter selection in the following algorithm.
If the subbagging estimator is biased (e.g., the estimator in Theorem (ref)), $k_N/\sqrt{N}\to\infty$ is required to achieve $\sqrt{N}$-consistency, and the asymptotic bias is of order $O(1/k_N)$. As a consequence, $k_N$ needs to be set as large as possible to reduce the bias. If the subbagging estimator is unbiased (e.g., the subbagging mean estimator in Theorem (ref), and subbagging the ordinary least squares estimation for linear regression), achieving $\sqrt{N}$-consistency does not require $k_N/\sqrt{N}\to\infty$.
}
From the perspectives of the computation efficiency and estimation efficiency, a small $\delta_m>0$ {\color{black}in the above algorithm} is preferred, because it results in not only a relatively small $m_N$ and but also a small asymptotic variance of $\hat\theta_{k_N,m_N}$, which is exactly equal to that of the full sample {\color{black}estimator} $\hat\theta_N$. Setting $\delta_m=0$ can lead to an even smaller $m_N$ and thus a faster computation speed, but with a sacrifice of the variance inflation rate $1/\alpha$ for the subbagging estimator.
{\color{black}
Based on the subbagging literature (see, e.g., buhlmann2003bagging and mentch2016quantifying), we draw $m_N$ subsamples independently and identically from $S_{N}^{k_N}$ with equal probability $1/{N\choose k_N}$ in Algorithm (ref). Each drawn subsample $s=\{s_1,\cdots, s_{k_N}\}$ is a subset of $\{1,\cdots,N\}$ with size $k_N$, and thus $\{s_1,\cdots, s_{k_N}\}$ should be sampled via the simple random sampling (SRS) without replacement. However, due to the memory constraint, subsamples with size $k_N$ need to be drawn directly from the data stored on the hard drive. In order to tackle this task, SRS algorithms without replacement from the hard drive need to be investigated. It is worth noting that those algorithms have been studied in the literature. For example, Gupta1984algorithm proposes an algorithm of SRS without replacement called GSEL, which has the time complexity of $O(k_N\log_2 k_N)$ and requires $O(k_N)$ storage locations. This algorithm meets the memory constraint and hence can realize SRS without replacement from the hard drive for massive data analysis.
}
Algorithm (ref) can {\color{black}result in a $\sqrt{N}$-consistent subbagging estimator} if the memory constraint allows $k_N=O(N^{1/2+\delta_k})$. {\color{black}However i}n real practice, $k_N$ is possibly restricted to be smaller by the memory constraint of a computer, e.g., $k_N=O(N^\delta)$ and $\delta\leq 1/2$. In this case, $k_N/\sqrt{N}\to \infty$ in Theorem (ref) and Algorithm (ref) is no longer satisfied{\color{black}, and hence it is possible that Algorithm (ref) cannot lead to a $\sqrt{N}$-consistent result. It is worth noting that the requirement $k_N/\sqrt{N}\to \infty$ in Theorem (ref) and Algorithm (ref) is more to make the square of the bias {\color{black}much} smaller than the variance in the subbagging procedure. Thus}, if we can reduce the bias {\color{black}when we solve the estimating equations ((ref)), the requirement for} large $k_N$ can be released. {\color{black}For example, in the following we reduce the bias of the covariance matrix estimator in Example (ref), and propose a $\sqrt{N}$-consistent subbagging estimator without requiring $k_N/\sqrt{N}\to \infty$.}
{\color{black}
Theorem (ref) indicates that the bias correction before the subbagging approach can lead to a $\sqrt{N}$-consistency estimator as $k_N$ grows slower than $\sqrt{N}$. This motivates us to consider a similar idea to deal with the general subbagging estimators. Specifically,} we {\color{black}first} propose three bias-corrected subsample {\color{black}estimator}s $\hat\theta_{k_N,s}^{(bc1)}$, $\hat\theta_{k_N,s}^{(bc2)}$ and $\hat\theta_{k_N,s}^{(bc3)}${\color{black}, relative to the estimator $\hat\theta_{k_N,s}$ that directly solves the estimating equations ((ref)).}
The first bias-corrected subsample {\color{black}estimator} is defined by $\hat\theta_{k_N,s}^{(bc1)}=\hat\theta_{k_N,s}-k_N^{-1}B_{\hat\theta_{k_N,s}},$ which is simply the subsample {\color{black}estimator} minus its bias estimation, where $B_{\theta}$ is given {\color{black}in Theorem (ref)}. {\color{black}Based on different model settings}, it is possible that $B_{\theta}$ is hard to derive or does not have a closed-form. In this case, we can replace $B_{\theta}$ by its sample version and let \[ \widehat V_{k_N,s}(\theta)=k_N^{-1}\sum_{i=1}^{k_N}\frac{\partial \psi_{\theta}(Z_{s_i})}{\partial \theta^\top},\widehat V^{(2)}_{k_N,s}(\theta)=k_N^{-1}\sum_{i=1}^{k_N}\frac{\partial^2 \psi_{\theta_0}(Z_{s_i})}{\partial \theta^\top\otimes \partial \theta^\top},\textrm{ and } \]
rilstone1996second and kim2016higher have shown $\widehat B_{k_N,s}(\hat\theta_{k_N,s})$ is a consistent estimator of $B_{\theta_0}$. Accordingly, we have the second bias-corrected subsample {\color{black}estimator} defined by $\hat\theta_{k_N,s}^{(bc2)}=\hat\theta_{k_N,s}-k_N^{-1}\widehat B_{k_N,s}(\hat\theta_{k_N,s}).$ Making use of $\widehat B_{k_N,s}(\theta)$ defined above, firth1993bias and kim2016higher further proposed a new {\color{black}estimator} $\hat\theta_{k_N,s}^{(bc3)}$ which solves the estimating equations \[\color{black}\sum_{i=1}^{k_N}\psi_\theta(Z_{s_i}) +\widehat V_{k_N,s}(\theta)\widehat B_{k_N,s}(\theta)=0, \] and $\hat\theta_{k_N,s}^{(bc3)}$ is the proposed third bias-corrected subsample {\color{black}estimator}.
Let $\hat\theta_{k_N,s}^{(bc)}$ be any one of the proposed three bias-corrected subsample {\color{black}estimator}s. Using these estimators, we obtain the bias-corrected subbagging estimator {\color{black}by} \[\color{black} \hat \theta_{k_N,m_N}^{(bc)}=\frac{1}{m_N}\sum_{s\in\mathcal{S}}\hat\theta_{k_N,s}^{(bc)}. \] {\color{black}The following theorem derives the asymptotic property of a general bias-corrected subbagging estimator $\hat\theta_{k_N,m_N}^{(bc)}$.
}
{\color{black} Different from Example (ref) and Theorem (ref), the general bias-corrected subbagging estimator $\hat\theta_{k_N,m_N}$ in Theorem (ref) is consistent, but still a biased estimator. However, the asymptotic bias is no longer from the mean of $k_N^{-1}{m_N^{-1}}\sum_{s\in{\mathcal{S}}}{(\mathcal{B}_{k_N,s}-B_{\theta_0})}$ in the asymptotic expansion, because using the results of Theorem (ref), we obtain ${\rm E}\{k_N^{-1}{m_N^{-1}}\sum_{s\in{\mathcal{S}}}{(\mathcal{B}_{k_N,s}-B_{\theta_0})}\}$ $=0$. As a consequence of Theorems (ref) -- (ref), we recommend the following algorithm for hyperparameter selection and obtaining a $\sqrt{N}$-consistent subbagging estimator under the more restrict memory constraint $k_N=O(N^{1/3+\delta_k})$ for some $0<\delta_k\leq 1/6$.
}
{\color{black} Compared to Algorithm (ref), the subbagging approach with bias correction in Algorithm (ref) can achieve $\sqrt{N}$-consistency under a smaller $k_N=O(N^{1/3+\delta_k})$. For example, the 12GiB American airline dataset analyzed in Section (ref) of this article has $N=$ 118,914,459 observations in total. It suffices to set $k_N=\lfloor N^{1/2+1/1000}\rfloor=$11,109 if we apply Algorithm (ref). However, $k_N$ can be reduced to $\lfloor N^{1/3+1/1000}\rfloor=500$ if we apply Algorithm (ref). It is worth noting that the reduction of $k_N$ causes the rise of $m_N$ in Algorithm (ref), and hence the computation cost increases.
}
Based on the idea {\color{black}of Algorithm (ref)}, one can expect that a second-order bias correction or higher-order bias-correction is needed for more strict memory {\color{black}constraints}, but similar bias correction analyses can be made to resolve this issue.
{\color{black}Analogously}, other bias correction methods can be {\color{black}implemented to reduce the bias}, e.g., jackknife and {\color{black}bootstrap}. However, the jackknife or bootstrap bias correction requires a further resampling procedure and solving estimating equations of the resampled data, which is computationally more expensive than the aforementioned bias correction methods.
In order to realize the {\color{black}inference} on parameter $\theta$, we need to estimate $\Xi_{\theta_0}\triangleq V_{\theta_0}^{-1} \Sigma_{\theta_0}( V_{\theta_0}^{-1})^\top$ which is the {\color{black}variance in the asymptotic normal distribution of the full sample estimator and subbagging estimators}. For the full sample {\color{black}estimator} $\hat\theta_N$, {\color{black}$\Xi_{\theta_0}$ can either be estimated} by $ \Xi_{\hat\theta_N}=V_{\hat\theta_N}^{-1} \Sigma_{\hat\theta_N}( V_{\hat\theta_N}^{-1})^\top$ or
The former variance estimator $ \Xi_{\hat\theta_N}$ is possibly not easy to obtain because $\Xi_{\theta}$ can be difficult to derive or does not have a closed-form. The latter variance estimator $\widehat\Xi_N(\hat\theta_N)$ requires the use of full data set $\{1,\cdots,N\}$ even if we have already computed $\hat\theta_N$, and hence does not satisfy the memory constraints $k_N=O(N^{1/2+\delta_k})$ and $k_N=O(N^{1/3+\delta_k})$ in Algorithms (ref) -- (ref).
Similar to bootstrap/bagging, subbagging also provides the variance estimation as $m_N\to \infty$. Specifically, we can consider the subbagging variance estimator
Using the theoretical results in {\color{black}Theorem (ref) and Theorem (ref), along with} the incomplete $U$-statistics theory of Lemma (ref) in Appendix, we obtain the consistency of the subbagging variance estimators in the following theorem.
Compared to the full sample variance estimators $\Xi_{\hat\theta_N}$ and $\widehat\Xi_N(\hat\theta_N)$, there are three advantages of using the subbagging variance estimators $k_N \widehat\Omega_{k_N,m_N}$ and $k_N \widehat\Omega^{(bc)}_{k_N,m_N}$. First, similar to bootstrap, obtaining the subbagging variance estimators does not require the derivation of the variance $\color{black}\Xi_{\theta_0}= V_{\theta_0}^{-1} \Sigma_{\theta_0}( V_{\theta_0}^{-1})^\top$. Second, the subbagging variance estimators can satisfy the memory constraint $k_N=O(N^{1/2+\delta_k})$ or $k_N=O(N^{1/3+\delta_k})$. Third, if we have already computed the subsample {\color{black}estimates as in Algorithms (ref) -- (ref)}, the computation of the subbagging variance estimator $k_N \widehat\Omega_{k_N,m_N}$ or $k_N \widehat\Omega^{(bc)}_{k_N,m_N}$ is more efficient than that of the full sample variance estimator $\widehat\Xi_N(\hat\theta_N)$, because computing $k_N \widehat\Omega_{k_N,m_N}$ or $k_N \widehat\Omega^{(bc)}_{k_N,m_N}$ only requires the aggregation of $m_N$ subsample {\color{black}estimate}s but not $N$ data points, and in both Algorithms (ref) -- (ref), $m_N=\lfloor \alpha N^{1+\delta_m}/k_N\rfloor\ll N$. This computation advantage is demonstrated in the real data analysis of Section (ref).
{\color{black}To assess the finite sample performance of our proposed algorithms and theoretical results, we conduct numerical studies} by simulating data from {\color{black}linear regression and logistic regression. The ordinary least squares estimation of linear regression and the maximum likelihood estimation of logistic regression are both under the estimation framework that we have discussed in Section (ref). In this section, we only present the results for logistic regression to save space. The results for linear regression and their related discussions are provided in Section (ref) of the supplementary material.
Specifically, for each observation $i\in\{1,\cdots,N\}$ from the logistic regression model}, the data are simulated by $Z_i=(Y_i,X_i)$ with a binary response $Y_i$ and $d=2$ exogenous covariates $X_{i}=(1, X_{i 1})^{\top},$ where $X_{i 1}$ is independently generated from the standard normal distribution. The corresponding repression coefficients of $X_{i}$ are set to be $\theta_0=(\theta_{10},\theta_{20})^\top=(0,1)^{\top} .$ The binary response $Y_{i}$ is generated from a Bernoulli distribution with the probability given by \[ P\left(Y_{i}=1 | X_{i}\right)=\frac{\exp \left(X_{i}^{\top} \theta_0\right)}{1+\exp \left(X_{i}^{\top} \theta_0\right)}. \] We consider three sample sizes $N=$ 2,000, 10,000, and 50,000. {\color{black}Since the maximum likelihood estimator of logistic regression is a biased estimator, to} satisfy the requirements of Algorithm (ref) and Algorithm (ref), we consider four subsample sizes $k_N=\lfloor N^{5/12}\rfloor$, $\lfloor N^{6/12}\rfloor$, $\lfloor N^{7/12}\rfloor$ and $\lfloor N^{8/12}\rfloor$, and two settings of $m_N=\lfloor \alpha N/k_N\rfloor$ and $\lfloor \alpha N^{4/3}/k_N\rfloor$ with $\alpha=1$ and $1/3$.
For each setting, all simulation results are obtained via 1,000 realizations. To evaluate the performance of parameter estimates, we define $\hat\theta^{(r)}=(\hat\theta_1^{(r)},\hat\theta_2^{(r)})^\top$ as the vector estimate of $\theta$ obtained via the full sample {\color{black}estimation} or Algorithm (ref) or Algorithm (ref) in the $r$-th realization. For each component of $\theta$, say $\theta_j$, the averaged bias of $\hat\theta_j^{(r)}$, $r=1,\cdots,1000$, is BIAS $=1000^{-1}\sum_{r}( \hat{\theta}_{j}^{(r)}-{\theta}_{j0})$, and the standard deviation of $\hat\theta_j^{(r)}$ is SD $=\{1000^{-1}\sum_{r}(\hat{\theta}_{j}^{(r)}-1000^{-1}{\sum_{r}\hat{\theta}_{j}^{(r)}})^2\}^{1/2}$. So the root mean squared error is RMSE $=\sqrt{\textrm{BIAS}^2+\textrm{SD}^2}$. For any generic square matrix $G$, let $(G)_{jj}$ be the $j$-th diagonal element of $G$. The asymptotic standard deviation (ASD) of $\hat \theta_j^{(r)}$ is approximated by $1000^{-1}\sum_r\sqrt{\big(N^{-1}\widehat \Xi_N^{(r)} (\theta_0)\big)_{jj}}$ and the averaged subbagging standard error (SSE) is $ 1000^{-1}\sum_r\sqrt{\big(N^{-1}k_N\widehat\Omega_{k_N,m_N}^{(r)}\big)_{jj}}$, where $\widehat \Xi_N^{(r)} (\cdot)$ and $\widehat\Omega_{k_N,m_N}^{(r)}$ are those defined in ((ref)), ((ref)) and ((ref)) but for the $r$-th realization. Due to the variance inflation of the subbagging estimators, we further define the $\alpha$-adjusted ASD and $\alpha$-adjusted SSE by $\sqrt{1+1/\alpha}$ times ASD and SSE, respectively. The SSE and $\alpha$-adjusted SSE, together with the asymptotic normality of the subbagging estimators, lead to 95% confidence intervals and $\alpha$-adjusted confidence intervals of ${\theta}_{j}$, respectively. If we denote either of the confidence intervals as CI$_j^{(r)}$, the empirical coverage probability for CI$_j^{(r)}$ is $\mathrm{CP}=1000^{-1} \sum_{r} I_{\{\theta_{j0} \in \mathrm{CI}_j^{(r)}\}}$, where $I_{\{\cdot\}}$ is an indicator function.
Table (ref) presents the BIAS, SD and RMSE of the full sample {\color{black}estimate}s with three sample sizes. We also report their memory usage (MEMORY) for programming in Python by utilizing an Intel Core i7 CPU (3.2GHz) with 32GB of 2667MHz DDR4 memory of a PC as a benchmark. To compare with the proposed subbagging estimates, Tables (ref) -- (ref) report the aforementioned performance measures for $\hat \theta_{k_N, m_N}$ in Algorithm (ref) and $\hat \theta_{k_N, m_N}^{(bc2)}$, $\hat \theta_{k_N, m_N}^{(bc3)}$ in Algorithm (ref), respectively, by setting $\alpha=1$. As the logistic regression does not have the closed-form for the first bias-corrected subsample {\color{black}estimator} $\hat\theta_{k_N,s}^{(bc1)}$ introduced in Section (ref), we only report the performance measures of subbagging estimates established based on $\hat\theta_{k_N,s}^{(bc2)}$ and $\hat\theta_{k_N,s}^{(bc3)}$. The results for $\alpha=1/3$ yield similar findings and are presented in Tables (ref) -- (ref) of the supplementary material to save space.
We obtain four interesting findings if the settings of $k_N$ satisfy the requirements of Algorithms (ref) -- (ref) (settings of $k_N$ in bold in Tables (ref) -- (ref){\color{black}, which are due to the biasedness of the maximum likelihood estimation for logistic regression}), respectively. The first is that the SD and the RMSE of the subbagging estimates are close to those of the full sample estimates when $N$ is large. This finding indicates the $\sqrt{N}$-consistency of the subbagging estimators. The second is that the subbagging estimates have larger BIAS compared to the full sample estimates, and larger $k_N$ can reduce the BIAS. The third is that the larger setting of $m_N=\lfloor \alpha N^{4/3}/k_N\rfloor$ leads to a better efficiency but not a smaller BIAS, and under this setting, the SD is almost identical to the ASD and SSE. However, the SD for the smaller setting of $m_N=\lfloor \alpha N/k_N\rfloor$ is much closer to the $\alpha$-adjusted ASD and $\alpha$-adjusted SSE, because the variance inflation rate $1/\alpha$ should be taken into the consideration as $k_Nm_N/N\to \alpha\in(0,\infty)$. The corresponding CP performances in Tables (ref) -- (ref) further confirm this result. These second and third findings are consistent with {\color{black}Theorem (ref) and Theorem (ref)}, which have shown that the bias of the subbagging estimators is determined by $k_N$ only but not $N$ and $m_N$, and the efficiency of the subbagging estimators can be improved by increasing $m_N$. Lastly, the memory usage of the subbagging estimates is significantly smaller than that of the full sample estimates, and smaller $k_N$ indeed leads to a smaller memory usage.
We also investigate the performance of the subbagging estimates if the settings of $k_N$ do not satisfy the requirements of Algorithms (ref) -- (ref). Table (ref) obtained by Algorithm (ref) shows that the ASD, SSE, CP and their adjusted versions all break down if a smaller $k_N=O(N^{1/2})$ is set. In contrast, Tables (ref) -- (ref) obtained by Algorithm (ref) have a significant reduce in BIAS under the same $k_N=O(N^{1/2})$ setting as $N$ is large. The corresponding estimates, their standard errors, and empirical coverages also perform well. Consequently, we recommend using Algorithm (ref) in practice when the memory of a computer is not large enough for Algorithm (ref).
For {\color{black}the illustration purpose}, we demonstrate the application of our proposed subbagging methods on {\color{black}a 12GiB} American airline dataset (\url{https://doi.org/10.7910/DVN/HG7NV7}){\color{black}, although our approach can definitely be applied to analyze terabytes of data in real practice. The dataset} contains the flight arrival and departure information for all commercial flights from 1987 to 2008 in US. {\color{black}The full sample size of data is $N=$ 118,914,459.} Each sample in the dataset corresponds to one flight record.
The research target of this dataset is to investigate factors that can influence the delayed status of a flight. To this end, we obtain the variable “Delayed", indicating whether or not the flight is delayed for arrival, as the response variable, and consider three variables as covariates, including the year of the flight (Year), the scheduled departure time (CRSDepTime, in HHmm), and the actual elapsed time (ActualElapsedTime, in minutes), to predict the response variable via logistic regression. It is worth noting that the {\color{black}12GiB} data contain more variables of flight records. However, for the sake of exogeneity and interpretation purposes, we select these three covariates. We further scale the covariates by $1/1000$ before the model fitting. The intercept and the three variables lead to the parameter estimates of their regression coefficients $\theta=(\theta_1,\cdots,\theta_4)^\top$.
To apply our subbagging algorithms to estimate the {\color{black}logistic} regression coefficients, we consider $k_N=\lfloor N^{1/2+1/1000}\rfloor=$11,109 and $k_N=\lfloor N^{1/3+1/1000}\rfloor=500$ along with $m_N=\lfloor \alpha N/k_N\rfloor$ to balance the computation efficiency and satisfy the minimum requirements of Algorithms (ref) -- (ref), respectively. Similar to simulation in Section (ref), we only report the results of the two bias-corrected subbagging {\color{black}estimate}s $\hat\theta_{k_N,m_N}^{(bc2)}$ and $\hat\theta_{k_N,m_N}^{(bc3)}$ for Algorithm (ref). Further, various settings of $\alpha=0.01$, 0.02, 0.04, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8 and 1 are investigated. To save space, we only report the parameter estimates under settings $\alpha=0.01$ and 0.2 in Table (ref), respectively. However, the $\alpha$-adjusted subbagging standard errors ($\alpha$-adjusted SSE, namely $\sqrt{1+1/\alpha}\times\sqrt{\big(N^{-1}k_N\widehat\Omega_{k_N,m_N}\big)_{jj}}$ used in Section (ref)) for all $\alpha$s are reported in Figure (ref).
For the illustration purpose, we program in Python via using an Intel Core i7 CPU (3.2GHz) with 32GB of 2667MHz DDR4 memory of a PC, which allows us to load the {\color{black}12GiB} data in the memory, obtain the response and three covariates, and compute the full sample {\color{black}estimate only for the purpose of} comparison. The loading of the {\color{black}full} dataset takes {\color{black}27.7 minutes} shown in Table (ref), while our proposed subbagging algorithms only require less than or equal to {\color{black}6.5} minutes for sampling from data on the hard {\color{black}drive as} $\alpha=0.01$. {\color{black} This fact demonstrates the advantage of the proposed subbagging approaches under smaller $\alpha$ even if the memory of a computer is enough to load the whole dataset in. It is worth noting that the reported loading time for the whole dataset refers to the time that it takes to read the CSV file into a “NumPy” array in Python. To accomplish this task, we apply the commonly used function “pandas.read_csv(CSV file).values” from the “Pandas” package. Note that compared to “read.csv” function in R, the loading time of using “Pandas” package in Python has already been improved (around 31 minutes if loading in R), not to mention that using other packages or functions in Python (e.g., “NumPy” function “numpy.genfromtxt”) can cause a loading failure for the 12 GiB full dataset.}
We also compare the execution times for computing estimates (Estimation Time) and computing SEs (SE Time) in Table (ref). The SE is the standard error for the full sample estimate or the $\alpha$-adjusted SSE for the subbagging estimate, and the SE Time is measured after the estimation is completed and saved in the memory. The results in Table (ref) show that the full sample SE takes 1.7 seconds to compute, while Algorithms (ref) -- (ref) compute SEs within 0.002 seconds for all settings. Further, the Estimation Times of $\hat \theta_{k_N,m_N}$ in Algorithm (ref) and $\hat \theta_{k_N,m_N}^{(bc2)}$, $\hat \theta_{k_N,m_N}^{(bc3)}$ in Algorithm (ref) under $\alpha=0.01$ are {\color{black}0.4, 2.0 and 3.3} minutes, respectively, much smaller than the {\color{black}20.9} minutes used for computing the full sample estimate. Note that it takes longer time in accomplishing Algorithm (ref) under $\alpha=0.2$. This is because that the sacrifice of the memory usage (MEMORY in Table (ref)) for Algorithm (ref) causes the increase of $m_N=\lfloor \alpha N/k_N\rfloor$ (from $m_N=$ 2,140 in Algorithm (ref) to $m_N=47,565$ in Algorithm (ref) under $\alpha=0.2$), which corresponds to larger computational complexity. {\color{black}But if parallel computing is considered}, since both Algorithms (ref) -- (ref) can {\color{black} actually be implemented in parallel}, the computation time for larger $\alpha$ can be reduced to Table (ref)'s time divided by the number of parallel {\color{black}processes}. It is also noted that all the memory usage of Algorithms (ref) -- (ref) is much smaller, i.e., only {\color{black}2666.1KiB $=$ 2.6MiB} and {\color{black}120.0KiB}, respectively, compared to the memory usage of the full sample estimation, {\color{black}12686071.1KiB $=$ 12.1GiB}. Therefore, our proposed algorithms can be used to meet various memory constraints of computers in real practice.
Table (ref) further shows two interesting findings in parameter estimates and SEs. Firstly, the estimates obtained by Algorithm (ref) -- (ref) are both numerically close to the full sample estimates. Secondly, the SEs for the full sample estimate and the subbagging estimates have the same order of magnitude under $\alpha=0.2$, and the SEs for the subbagging estimates under $\alpha=0.01$ are slightly larger. This finding is not surprising since all of the estimators considered in Table (ref) are $\sqrt{N}$-consistent, but due to the variance inflation rate $1/\alpha$ for the subbagging estimator, smaller $\alpha$ can result in a larger variance. Figure (ref) further compares the SEs for the subbagging estimates under larger $\alpha$s and the SE for the full sample estimate. It is found in Figure (ref) that they are almost identical as $\alpha=1$.
To investigate the trade-off between the estimation efficiency and the computation efficiency of the proposed subbagging estimator, we report the $\alpha$-adjusted SSE versus $\alpha$ in Figure (ref) for parameter estimates $\hat\theta=(\hat\theta_1,\cdots,\hat\theta_4)^\top$ obtained by Algorithm (ref). Note that larger $\alpha$ is equivalent to larger $m_N=\lfloor \alpha N/k_N\rfloor$, which is less computationally efficient, but more efficient in estimation. Based on Figure (ref), if one would like to obtain a subbagging estimate whose SE is as close as that of the full sample estimate, then a larger $\alpha$ is preferred. However, in order to balance the computation cost, we can find in Figure (ref) that when $\alpha>0.2$, the reduction of SE is quite limited. {\color{black}As a consequence, we recommend using $\alpha=0.2$}. Similar analyses can be made in practice for selecting an optimal $\alpha$ via a “scree plot” of Figure (ref).
It is worth noting that the $\alpha$-adjusted SSE for the estimate of $\theta_j$ ($j=1,\cdots,4$) in Figure (ref), namely $\sqrt{1+1/\alpha}\times\sqrt{\big(N^{-1}k_N\widehat\Omega_{k_N,m_N}\big)_{jj}}$, is computed based on a given $m_N=\lfloor \alpha N/k_N\rfloor$, which requires a longer time to compute for larger $\alpha$ (e.g., see Table (ref) given $\alpha=0.2$). So Figure (ref) is not convenient enough for users to fast select an optimal $\alpha$ in real practice. Note that using {\color{black}Theorem (ref)}, we have that $\sqrt{1+1/\alpha}\times\sqrt{\big(N^{-1}k_N\widehat\Omega_{k_N,\lfloor 0.01 N/k_N\rfloor}\big)_{jj}}$ is also a consistent SE for the subbagging estimate, and $\sqrt{\big(N^{-1}k_N\widehat\Omega_{k_N,\lfloor 0.01 N/k_N\rfloor}\big)_{jj}}$ is a consistent SE for the full sample {\color{black}estimate}. We call these standard errors anticipated SEs in Figure (ref). The advantage of reporting the anticipated SEs is that we only need to draw a small number of subsamples ($m_N=\lfloor 0.01 N/k_N\rfloor$) to compute those SEs, and hence they can be obtained within a short period of time ({\color{black}0.9} minutes only to {\color{black}generate} Figure (ref)).
{\color{black}In addition}, we are capable of anticipating the computation time of using a PC for any subbagging estimation algorithms given $\alpha$ ($m_N=\lfloor \alpha N/k_N\rfloor$), as long as we know the actual computation time given $\alpha=0.01$ ($m_N=\lfloor 0.01 N/k_N\rfloor$). The idea behind it is that the number of drawn subsamples $m_N$, which determines the computational complexity of the algorithms, increases linearly with respect to $\alpha$. Accordingly, in Figure (ref) we report {\color{black}the} anticipated computation time (the sum of Loading Time, Estimation Time and SE Time) for Algorithm (ref) given $\alpha$, which is equal to $\alpha/0.01$ times the actual computation time given $\alpha=0.01$. {\color{black}Accordingly}, Figure (ref) reveals how much time it will take to reach some certain level of estimation efficiency for the proposed subbagging estimation.
In summary, Table (ref) given $\alpha=0.01$ provides a fast {\color{black}sketch on} the big dataset with ensuring $\sqrt{N}$-consistency for the estimation, while Figure (ref) can help users make decisions on whether or not it is worth taking more time or drawing more subsamples in parallel to achieve higher estimation efficiency. All of the above findings demonstrate the usefulness of the subbagging methods for {\color{black}econometric and statistical analysis} on massive data.
In this article, in order to meet different memory constraints under the massive data setting, we introduce two computationally efficient subbagging {\color{black}Algorithms (ref) -- (ref) }. Both subbagging algorithms lead to $\sqrt{N}$-consistent estimators that have the same order of estimation efficiency as the full sample estimator. {\color{black}By the principle of divide and conquer, both Algorithms (ref) -- (ref) can be realized in parallel, and hence are still computationally inexpensive.} The first algorithm uses the simple average subbagging but {\color{black}generally} requires a larger memory of $O(N^{1/2+\delta_k})$ for computation, where $\delta_k>0$. The second algorithm adopts the bias-corrected subbagging method and can satisfy a more strict memory constraint $O(N^{1/3+\delta_k})$. Obtaining the limiting distributions of the resulting estimators based on these two algorithms is quite challenging due to the fact that one subsample can have overlap with the other subsample in ensembles; see detailed derivations via utilizing the incomplete $U$-statistics theory in the supplementary material. We subsequently obtain that both algorithms provide asymptotically normal estimators but with a variance inflation compared to the full sample estimator. These properties allow us to {\color{black}perform} hypothesis testing and {\color{black}construct confidence intervals.} The performances of the proposed algorithms and the corresponding {\color{black}theoretical} properties are supported by both simulation studies and a real American airline dataset analysis.
We conclude this work by identifying two potential avenues for future research. First, the first-order bias-correction method proposed in this paper can be extended to second-order or higher-order bias-correction, which is expected to satisfy more strict memory constraints. Second, {\color{black}though the proposed subbagging approach has been analyzed under a general estimation framework of solving estimating equations in this article, other complex estimating problems can also be discussed, especially when the data have a probabilistic structure. In this case, having more data (considering a larger subsample size or drawing more subsamples) is not necessarily more desirable boivin2006more}. We believe these extensions would further strengthen the usefulness of the subbagging methods for big data analysis.