Response-Aided Score-Matching Representative Approaches for Big Data Analysis and Model Selection under Generalized Linear Models ^†

Abstract

In this paper, we propose an efficient method called the response-aided score-matching representative (RASMR) approach to facilitate massive data model selection and data analysis with generalized linear models (GLMs) and a predetermined data partition due to data localization. Similar to the original score-matching representative (SMR) approach, RASMR constructs an artificial data point, called the representative, for each data block. It then fits a GLM on the representative dataset, which provides not only an efficient approach for massive data analysis but also an ideal solution in response to privacy concerns by avoiding the transfer of sensitive data. By further splitting the data blocks according to the values of the response variables, RASMR can obtain more accurate parameter estimates than SMR. Furthermore, by theoretical justifications and simulation studies, we show that RASMR can be more efficiently utilized for model selection and variable selection for a massive dataset by approximating the Akaike information criterion (AIC) and the aggregated prediction errors for cross-validation, which are commonly used for choosing the most appropriate statistical model and drawing reliable conclusions. We also apply the proposed RASMR approach to the airline on-time performance data, which consists of 371 data files labeled by month, and show that RASMR can be successfully used for selecting the most appropriate model for real massive data analysis.

1. Introduction

The numerous innovations in data analysis in the last decade have dramatically affected the technologies used in our daily lives. Datasets with unprecedented sizes and complexities are rapidly generated and collected from a great variety of resources [1]. While big datasets bring us incredible opportunities for new discoveries, many traditional methods that perform well for moderate sample sizes are no longer realistic for analyzing massive amounts of data [2].

To address the big data challenges, many traditional statistical tools have been reinvented or adapted to deal with the gigantic volume or size of big data, which is a major goal of our work. Comprehensive reviews, such as [3,4], have been provided for algorithmic solutions to big data problems, including divide-and-conquer, subsampling-based approaches, stochastic gradient descent, and online updating.

Divide-and-conquer approaches split the big data into manageable blocks, extract local summaries from each data block, and then generate overall insight. Various types of algorithms have been proposed to adapt classic statistical tools to big data problems [5,6,7,8,9,10,11,12].

Subsampling-based approaches draw subsamples from the original data by elaborately designed sampling mechanisms and then perform downstream inference and prediction based on the subsamples. For example, leveraging for big data regression, [13] constructed nonuniform sampling probabilities so that influential data points were sampled with high probabilities, and [14] proposed a novel method called information-based optimal sub-data selection (IBOSS), which selects samples from a big dataset based on the D-optimality criterion; it has been extended to various statistical models, such as logistic regression [15].

When multiple subsamples of the original data are available, different approaches have been proposed to aggregate the estimates obtained from different subsamples, including bagging or bootstrap aggregating [16], stacking [17,18], magging or maximin aggregating [19], neagging or normalized entropy aggregating [20,21]. These aggregation approaches are especially useful for inhomogeneous large-scale data under regression analysis.

Stochastic gradient descent (SGD) algorithms [22,23,24] update in a sequential manner based on a noisy gradient. Local SGD is the key building block of different federated learning algorithms [25,26], which were discussed for both homogeneous nodes [27,28,29,30,31,32] and heterogeneous nodes [29,33,34,35,36,37,38,39,40,41,42].

As stated in [43], new challenges in the big data era involve data security and localization legislation. According to [44,45], legislation on data localization has emerged as a global trend, with more nations developing laws that prohibit the transfer of sensitive data. Data localization may not kill global cooperation in data science, but it will undoubtedly create unpleasant barriers to data communication. A novel approach called the score-matching representative (SMR) was proposed by [43] for analyzing big data under communication restraints. Given an existing data partition, such as data blocks labeled by countries, regions, or sources, SMR constructs model-specified data representative(s) for each block and performs downstream analysis on the constructed representatives. Unlike the divide-and-conquer or subsampling approaches, SMR does not require communicating the actual data but rather their representatives, which are not part of the original data and, thus, avoid their transfer. Comprehensive studies show that the accuracy of the estimated model parameters based on SMR can be comparable to the full data estimates [43]. According to [43], the computational complexity of SMR is $O\left(Np\right)$ given the sample size N and the number p of model parameters.

Nevertheless, there are three hidden traps in [43]’s SMR method (see Section 2.1 for more details), as follows: (1) there can be multiple solutions to the score-matching equations; (2) the constructed representatives may be quite far away from their data blocks; and (3) the SMR algorithm may fail to converge. Due to these issues, SMR may not be accurate enough for more elaborate tasks, such as model selection and variable selection.

In this paper, we develop a new representative approach, called the response-aided score-matching representative (RASMR) approach, for big data analysis under generalized linear models (GLMs). It improves SMR significantly by splitting the data blocks further with the aid of the response variable (see Figure 1 for a graphical display). Since the refined data blocks yield only one solution to the score-matching equation, RASMR ensures the uniqueness of the data representatives. Compared with SMR, RASMR can not only provide much more accurate estimates for model parameters but also be used for model selection and variable selection efficiently, which are critical for data scientists to draw reliable conclusions from data analysis.

2. RASMR for Big Data Analysis under GLM

2.1. SMR Approach for GLM

In this section, we review the key components of the original SMR approach proposed by [43], along with its potential issues.

Following the notation of [43], the original data $\{({\mathbf{x}}_i,y_i),i=1,\dots ,N\}$ are provided with an index partition $\{I_1,\dots ,I_K\}$ of $I=\{1,\dots ,N\}$, where ${\mathbf{x}}_i\in {\mathbb{R}}^d$ represents the ith covariate vector, $y_i\in \mathbb{R}$ is the corresponding response, and $I={\cup }_{j=1}^K{I}_j$, $I_j\ne \varnothing$ for each j, and $I_i\cap I_j=\varnothing$ for each $i\ne j$.

Under a generalized linear model (GLM) [46,47], there is a link function g, p known predictor functions $h_1,\dots ,h_p$, and p unknown regression parameters ${\beta }_1,\dots ,{\beta }_p$, such that the expectation of the response variable $Y_i$ given ${\mathbf{x}}_i$ satisfies the following:

$$E(Y_i\mid {\mathbf{x}}_i)={\mu }_i \mathrm{a}\mathrm{n}\mathrm{d} {\eta }_i=g\left({\mu }_i\right)={\mathbf{X}}_i^T\mathit{\beta },$$

(1)

where ${\mathbf{X}}_i={(h_1\left({\mathbf{x}}_i\right),\dots ,h_p\left({\mathbf{x}}_i\right))}^T$, $i=1,\dots ,N$, and $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_p)}^T$.

Given the GLM (1), we denote the kth data block by ${\mathcal{D}}_k=\{({\mathbf{X}}_i,y_i),i\in I_k\}$ and let $s_k\left(\mathit{\beta }\right)={\sum }_{i\in I_k}[y_i-G\left({\eta }_i\right)]\nu \left({\eta }_i\right){\mathbf{X}}_i$ be the contribution made by the kth data block to the score function $s\left(\mathit{\beta }\right)={\sum }_{k=1}^K{s}_k\left(\mathit{\beta }\right)$ [43,46], where $G\left(\eta \right)=g^{-1}\left(\eta \right)$, $\nu \left(\eta \right)=G^{\prime }\left(\eta \right)/h\left(\eta \right)$, and $h\left({\eta }_i\right)=\mathrm{Var}\left(Y_i\right)$ are functions of $\eta$ or ${\eta }_i$. According to Section 2.5 of [46], the maximum likelihood estimate (MLE) of $\mathit{\beta }$ solves the score equation $s\left(\mathit{\beta }\right)=0$. The SMR algorithm [43] was designed to find ${\tilde{y}}_k\in \mathbb{R}$ and ${\tilde{\mathbf{X}}}_k\in {\mathbb{R}}^p$ solving $s_k\left(\mathit{\beta }\right)=n_k[{\tilde{y}}_k-G\left({\tilde{\eta }}_k\right)]\nu \left({\tilde{\eta }}_k\right){\tilde{\mathbf{X}}}_k\stackrel{▵}{=}{\tilde{s}}_k\left(\mathit{\beta }\right)$, where $n_k=\left|I_k\right|$ is the size of $I_k$, and ${\tilde{\eta }}_k={\tilde{\mathbf{X}}}_k^T\mathit{\beta }$. More specifically, the SMR algorithm first chooses the following:

$${\tilde{y}}_k={\left[\sum _{i\in I_k}\nu \left({\eta }_i\right){\eta }_i\right]}^{-1}\sum _{i\in I_k}\nu \left({\eta }_i\right){\eta }_i{y}_i,$$

(2)

then solves the score-matching equation as follows:

$$n_k\nu \left({\tilde{\eta }}_k\right)[{\tilde{y}}_k-G\left({\tilde{\eta }}_k\right)]{\tilde{\eta }}_k=\sum _{i\in I_k}\nu \left({\eta }_i\right)[y_i-G\left({\eta }_i\right)]{\eta }_i$$

(3)

for ${\tilde{\eta }}_k\in \mathbb{R}$, and then constructs the kth representative $({\tilde{\mathbf{X}}}_k,{\tilde{y}}_k)$ by calculating the following:

$${\tilde{\mathbf{X}}}_k={\left\{n_k\nu \left({\tilde{\eta }}_k\right)[{\tilde{y}}_k-G\left({\tilde{\eta }}_k\right)]\right\}}^{-1}\sum _{i\in I_k}\nu \left({\eta }_i\right)[y_i-G\left({\eta }_i\right)]{\mathbf{X}}_i.$$

(4)

Then the score function $s\left(\mathit{\beta }\right)={\sum }_{i=1}^N[y_i-G\left({\eta }_i\right)]\nu \left({\eta }_i\right){\mathbf{X}}_i={\sum }_{k=1}^K$ $s_k\left(\mathit{\beta }\right)={\sum }_{k=1}^K{\tilde{s}}_k\left(\mathit{\beta }\right)\stackrel{▵}{=}\tilde{s}\left(\mathit{\beta }\right)$. The MLE $\widehat{\mathit{\beta }}$ of $\mathit{\beta }$ based on the full data, which solves the score equation $s\left(\mathit{\beta }\right)=0$, is expected to be the same as the SMR estimate $\tilde{\mathit{\beta }}$ based on the weighted representative data $\{(n_k,{\tilde{\mathbf{X}}}_k,{\tilde{y}}_k),k=1,\dots ,K\}$, which solves the matched score equation $\tilde{s}\left(\mathit{\beta }\right)=0$. In practice, since $\mathit{\beta }$ is unknown, the SMR algorithm first solves the matched score function $\tilde{s}\left(\mathit{\beta }\right)=0$ with some initial representatives for $\tilde{\mathit{\beta }}$, and then solves the score-matching Equation (3) for the representatives. The procedure may continue iteratively until reaching a predetermined accuracy level. Three iterations were suggested by [43] for general applications, whose computational complexity is $O\left(Np\right)$.

Given the successful applications of SMR in estimating parameter values for big data analysis using a GLM, its accuracy level may not be high enough for more delicate applications, including model selection and variable selection, due to the following three issues. Firstly, there can be more than one solution for Equation (3) (see Section 2.2). SMR chooses the solution whose representative is closest to the mean representative (MR, that is, ${\overline{\mathbf{X}}}_k=n_k^{-1}{\sum }_{i\in I_k}{\mathbf{X}}_i,{\overline{y}}_k=n_k^{-1}{\sum }_{i\in I_k}{y}_i$), which, however, may not match the likelihood well. Secondly, the predictor representative ${\tilde{\mathbf{X}}}_k$ obtained from Equation (4) may be far away from its corresponding data block (see Section S2 in the Supplementary Materials), which may not represent its data block well. Thirdly, the SMR algorithm may not converge well in practice, which may occur along with misspecified link functions or highly skewed predictor distributions (see Section 4).

2.2. Solving the Score-Matching Equation with Splitting Points

In this section, we investigate the number of solutions for the score-matching Equation (3) and explain why more stable and accurate solutions can be obtained for Equation (3) by further splitting the data blocks, which is our motivation to propose the RASMR algorithms in later sections.

According to Equations (2) and (4), the kth representative $({\tilde{\mathbf{X}}}_k,{\tilde{y}}_k)$ can be obtained as weighted averages of the kth data block ${\mathcal{D}}_k=\{({\mathbf{X}}_i,y_i),i\in I_k\}$. The weights of $y_i$’s and ${\mathbf{X}}_i$’s are $\nu \left({\eta }_i\right){\eta }_i$ and $\nu \left({\eta }_i\right)[y_i-G\left({\eta }_i\right)]$, respectively. To stabilize the constructed representatives, we propose to keep the weights with the same sign (that is, all positive or all negative) in each individual data block. That is, besides splitting the data block according to $\mathrm{sgn}\left({\eta }_i\right)$, as suggested in Remark 3.1 of [43], we suggest further splitting the data blocks by $\mathrm{sgn}(y_i-G\left({\eta }_i\right))$.

Now, we investigate the number of solutions to Equation (3) and how to further split the data block to ensure a unique solution. Following [43], we denote $S\left(\eta \right)=\nu \left(\eta \right)[{\tilde{y}}_k-G\left(\eta \right)]\eta$. According to the proof for Theorem 3.1 in [43], Equation (3) can be rewritten as follows:

$$S\left({\tilde{\eta }}_k\right)={\displaystyle \frac{1}{n_k}}\sum _{i\in I_k}S\left({\eta }_i\right).$$

(5)

We further denote ${\eta }_k^{\wedge }={min}_{i\in I_k}\left\{{\eta }_i\right\}$ and ${\eta }_k^{\vee }={max}_{i\in I_k}\left\{{\eta }_i\right\}$. If both $\nu \left(\eta \right)$ and $G\left(\eta \right)$ are continuous, then there exists a solution ${\eta }_{*}\in [{\eta }_k^{\wedge },{\eta }_k^{\vee }]$ that solves Equation (5) (see Theorem 3.1 in [43]).

Examples of $\nu \left(\eta \right)$ and $G\left(\eta \right)$ for commonly used GLMs can be found in Table 1 of [43], which are all continuous. If $S\left(\eta \right)$ is strictly monotone on $[{\eta }_k^{\wedge },{\eta }_k^{\vee }]$, then there exists a unique ${\eta }_{*}\in [{\eta }_k^{\wedge },{\eta }_k^{\vee }]$ that solves (3) (see Theorem A1 in Appendix A).

If $\nu \left(\eta \right)$ is a constant, such as for the normal model with an identity link, the Bernoulli model with a logit link, the Poisson model with a log link, the gamma model with a reciprocal link, and the inverse Gaussian model with an inverse-square link, without any loss of generality, we rewrite $S\left(\eta \right)=[{\tilde{y}}_k-G\left(\eta \right)]\eta$. Then, its first derivative $S^{\prime }\left(\eta \right)={\tilde{y}}_k-\left[G\left(\eta \right)+G^{\prime }\left(\eta \right)\eta \right]$, with its key component $T\left(\eta \right)\stackrel{▵}{=}G\left(\eta \right)+G^{\prime }\left(\eta \right)\eta$. In this case, we have the following: ${\tilde{y}}_k={\displaystyle \frac{{\sum }_{i\in I_k}{\eta }_i{y}_i}{{\sum }_{i\in I_k}{\eta }_i}}={\displaystyle \frac{n_k^{-1}{\sum }_{i\in I_k}{\eta }_i{y}_i}{{\overline{\eta }}_k}}$, where ${\overline{\eta }}_k=n_k^{-1}{\sum }_{i\in I_k}{\eta }_i$ corresponds to the mean representative (MR). We denote the following: ${\tilde{G}}_k={\displaystyle \frac{{\sum }_{i\in I_k}{\eta }_{i}G\left({\eta }_i\right)}{{\sum }_{i\in I_k}{\eta }_i}}={\displaystyle \frac{n_k^{-1}{\sum }_{i\in I_k}{\eta }_{i}G\left({\eta }_i\right)}{{\overline{\eta }}_k}}$. Then, we have the following:

$$\overline{S}\stackrel{▵}{=}{\displaystyle \frac{1}{n_k}}\sum _{i\in I_k}S\left({\eta }_i\right)={\displaystyle \frac{1}{n_k}}\sum _{i\in I_k}[{\tilde{y}}_k-G\left({\eta }_i\right)]{\eta }_i={\overline{\eta }}_k({\tilde{y}}_k-{\tilde{G}}_k).$$

(6)

For the normal model with an identity link, that is, the usual linear model, there can be up to two solutions to (3) (see Theorem A2 in Appendix A).

For Bernoulli models (see Table 1 of [43]), $y_i\in \{0,1\}$ and $G\left({\eta }_i\right)\in (0,1)$ for all i. Then, $y_i<G\left({\eta }_i\right)$ always implies $y_i=0$, and $y_i>G\left({\eta }_i\right)$ always implies $y_i=1$. Suppose (i) either ${\eta }_i>0$ for all $i\in I_k$ or ${\eta }_i<0$ for all $i\in I_k$; and (ii) either $y_i>G\left({\eta }_i\right)$ for all $i\in I_k$ or $y_i<G\left({\eta }_i\right)$ for all $i\in I_k$ . Then, ${\tilde{y}}_k$ is either 0 or 1. Depending on ${\tilde{y}}_k=0$ or 1, we denote $S\left(\eta \right)$ as $S_0\left(\eta \right)$ or $S_1\left(\eta \right)$, and the potential splitting point as ${\eta }_l$ or ${\eta }_r$ (see

Table 1. Splitting points for Bernoulli models with different links.

Link Function	${\mathit{\eta }}_{\mathit{l}}$ for ${\mathit{S}}_0\left(\mathit{\eta }\right)$	${\mathit{\eta }}_{\mathit{r}}$ for ${\mathit{S}}_1\left(\mathit{\eta }\right)$
logit	−1.278464542761	1.278464542761
probit	−0.839923675692	0.839923675692
cloglog	−1	0.729114174900
loglog	−0.729114174900	1
cauchit	−0.801916425045	0.801916425045

), respectively. For Bernoulli models with logit, probit, cloglog, loglog, or cauchit links, a data block $I_k$ after splitting at ${\eta }_l$ or ${\eta }_r$ yields a unique solution solving (3) (see Theorem A3 in Appendix A and Lemmas S1–S5 in Section S6 of the Supplementary Materials). We summarize in Table 1 the potential splitting points ${\eta }_l$ for blocks with ${\tilde{y}}_k=0$ and ${\eta }_r$ for blocks with ${\tilde{y}}_k=1$.

For the Poisson model with a log link, if either ${\eta }_i>0$ for all $i\in I_k$ or ${\eta }_i<0$ for all $i\in I_k$, then there are up to two solutions solving (3) (see Theorem A4 in Appendix A).

For the gamma model with a reciprocal link, ${\overline{\eta }}_k=n_k^{-1}{\sum }_{i\in I_k}{\eta }_i$ is the unique solution solving (3) (see Theorem A5 in Appendix A).

For the inverse Gaussian model with an inverse-square link, in general, there are up to two solutions solving (3) (see Theorem A6 in Appendix A).

In conclusion, the score-matching Equation (3) often yields two solutions for commonly used GLMs, confirming the first hidden trap mentioned in the Introduction section for the original SMR algorithm. On the other hand, there are, at most, two solutions for Equation (3) (see Theorems A2–A6 in Appendix A), which motivates us to further split the corresponding data block in a way that each sub-block yields a unique solution.

2.3. Response-Aided Score-Matching Representative Approach

In this section, we propose a new algorithm (see Algorithm 1 for its pseudocode) with further splits based on the response variable $y_i$’s, and call it the response-aided score-matching representative (RASMR) approach. This suggests that further splits may at most quadruple the original number of data blocks. Since the time complexity of the original SMR is $O\left(Np\right)$ [43], which does not depend on the number of data blocks, the RASMR algorithm consumes no significantly more time than the SMR algorithm (see Section 4.2). Its time complexity is $O\left(Np\right)$ as well.

Algorithm 1: RASMR.

Data: $\mathcal{D}=\{({\mathbf{X}}_i,y_i),i=1,\dots ,N\}$ with a partition of K blocks indexed by $\{I_1,\dots ,I_K\}$. Denote the kth data block ${\mathcal{D}}_k=\{({\mathbf{X}}_i,y_i),i\in I_k\}$. Result: Parameter estimate $\tilde{\mathit{\beta }}$ of a given GLM with a predetermined number T of iterations. Calculate the initial weighted representative data: ${\tilde{\mathcal{D}}}^{\left(0\right)}=\left\{\right(n_k,{\tilde{\mathbf{X}}}_k^{\left(0\right)}={\overline{\mathbf{X}}}_k,$ ${\tilde{y}}_k^{\left(0\right)}={\overline{y}}_k{\left)\right\}}_{k=1}^K$ with ${\overline{\mathbf{X}}}_k=n_k^{-1}{\sum }_{i\in I_k}{\mathbf{X}}_i$ and ${\overline{y}}_k=n_k^{-1}{\sum }_{i\in I_k}{y}_i$, that is, the mean representatives Implement the iteratively reweighted least squares (IRLS) procedure [48] on ${\tilde{\mathcal{D}}}^{\left(0\right)}$ to obtain the initial estimate ${\tilde{\mathit{\beta }}}^{\left(0\right)}$

Remark 1.

The iteratively reweighted least squares (IRLS or IWLS) procedure has been widely used for finding the MLE of a GLM [48,49,50]. Nevertheless, some more robust variants of IRLS have been proposed to make the estimates less sensitive to outliers (see [51] and references therein). One may choose a more robust procedure than IRLS in Algorithm 1, given that the targeted parameter estimate $\widehat{\mathit{\beta }}$ is under the same criterion.

Remark 2.

Following the splits described in (2) of Algorithm 1, according to Theorems A2–A6 in Appendix A, many sub-blocks yielded a unique solution to (3). There are still leftover cases for normal/linear, Bernoulli, Poisson, and inverse Gaussian models, under which, some sub-blocks may yield up to two solutions to (3). For those cases, according to the proofs of Theorems A2–A4 and A6 (see Section S6 in the Supplementary Materials), we may further split such a block into two sub-blocks according to the peak value η of $S\left(\eta \right)$, which is ${\displaystyle \frac{1}{2}}{\tilde{y}}_{kl}$ for normal/linear models, ${\eta }_l$ or ${\eta }_r$ for Bernoulli models, $u\left({\tilde{y}}_{kl}\right)$ for the Poisson models, or ${\left(4{\tilde{y}}_{kl}^2\right)}^{-1}$ for the inverse Gaussian models. Based on our experience, such a split often leads to data blocks with a unique solution to (3).

2.4. RASMR with the Delta Ratio Split

As we will demonstrate in Section 4.2, RASMR significantly improves upon SMR in estimating model parameters. However, it may not perform sufficiently well when approximating likelihood with non-Gaussian covariates ${\mathbf{x}}_i$ or predictors ${\mathbf{X}}_i$ for model selection purposes.

To investigate this, we start with the mean representative (MR) of the predictors in block $I_k$, ${\overline{\mathbf{X}}}_k=n_k^{-1}{\sum }_{i\in I_k}{\mathbf{X}}_i$, which is a natural choice for the block center. Then, the predictor radius of $I_k$ can be defined as ${\Delta }_k=\underset{i\in I_k}{max}\Vert {\mathbf{X}}_i-{\overline{\mathbf{X}}}_k\Vert$, where ‖·‖ is the Euclidean norm. We call the relative distance of ${\tilde{\mathbf{X}}}_k$ from ${\overline{\mathbf{X}}}_k$, the delta ratio, defined as ${\tilde{\delta }}_k=\Vert {\tilde{\mathbf{X}}}_k-{\overline{\mathbf{X}}}_k\Vert /{\Delta }_k$ .

By definition, ${\tilde{\delta }}_k=0$ if $n_k=1$. If ${\tilde{\delta }}_k>1$ for some $I_k$, then its predictor representative ${\tilde{\mathbf{X}}}_k$ is outside the corresponding data block, which implies that it may not be a good representative for calculating the likelihood. Compared with SMR, RASMR’s delta ratios are significantly smaller (see Figure S1 in the Supplementary Materials), which partly explains why RASMR improves SMR significantly. Nevertheless, the delta ratios of RASMR may tend to inflate when the distribution of covariates is too extreme or complicated (see Section S2 in the Supplementary Materials for more details).

If ${\tilde{\delta }}_k$ is greater than a predefined threshold ${\delta }_0$, a fourth layer of further splitting at the mean ${\overline{\eta }}_k=n_k^{-1}{\sum }_{i\in I_k}{\eta }_i$ of ${\eta }_i$’s may be applied. Since the sample mean is sensitive to outliers, using ${\overline{\eta }}_k$ as the cutoff point of the split may enable us to separate the outliers from the majority of the data block, which will make the leftover members in the data blocks closer to each other. In other words, we may identify the minorities of data points in the original data block and separate them into their own data block.

In practice, if the likelihood approximation is not a serious concern, setting the threshold for delta ratios to be one is a conservative choice, since it only requires the representatives to stay inside its data block, not necessarily very close to the center (see Section 4.3 for further discussion on choosing ${\delta }_0$). The pseudocode of the RASMR algorithm with the delta ratio split is described in Algorithm 2.

Algorithm 2: RASMR algorithm with the delta ratio split.

2.5. Learning Rate Scheduling

As an iterative method to maximize the log-likelihood function, the RASMR approach may suffer from a convergence issue, since the solution at each iteration may be far away from the global maximizer, or jump back and forth at the final stages (see Figure S13 in the Supplementary Materials for such an example).

The convergence issue of RASMR comes from various reasons. It could be caused by a bad initial estimate via the mean representative (MR). A small number of data blocks may also make the convergence poor (see Figure S12 in the Supplementary Materials for such an example; see also Section S1 for generating data blocks), especially when the number of parameters is moderately large (see Figure S13). The inherent complexities of the original data, such as ultra-high-dimensional, highly skewed, or heavily tailed distributions of covariates or predictors, may cause similar problems as well.

To make the RASMR approach more robust and practically useful for big data analysis, we adopt a learning rate scheduling strategy, which is a popular idea in the modern development of machine learning and artificial intelligence to resolve convergence issues (see, e.g., [52,53], and references therein). More specifically, we update the RASMR parameter estimates with an exponential learning rate decay as follows:

$${\mathit{\beta }}^{\left(t\right)}={\mathit{\beta }}^{(t-1)}+e^{-\theta t}({\tilde{\mathit{\beta }}}^{\left(t\right)}-{\mathit{\beta }}^{(t-1)}),$$

(7)

where ${\tilde{\mathit{\beta }}}^{\left(t\right)}$ is the original RASMR parameter estimate after the tth iteration, $\theta >0$ is a hyperparameter to control the learning rate $e^{-\theta t}$. Apparently, if $\theta =0$ or with the learning rate 1, (7) yields the original RASMR estimates. A larger value of $\theta$ posts a quicker rate of decay, and will force the convergence of iterative parameter estimates. On the other hand, if the learning rate quickly shrinks to zero, the estimate may be updated very slowly, making the improvement of later iterations negligible. To overcome the issue that the learning rate decays too fast, we suggest truncating the learning rate at, say, the 10th iteration and making the learning rate constant afterward. The truncation strategy for the exponential learning rate decay works well when the number of variables is as large as 100 (see Figure S13).

3. Model Selection and Variable Selection Using RASMR

3.1. Information-Based Criteria and Model Selection

During the process of RASMR parameter estimation by Algorithms 1 and 2, two very useful side products are also generated, namely the response representative ${\tilde{y}}_k$ and the predictor representative ${\tilde{\mathbf{X}}}_k$ . The quadruple $(n_k,{\tilde{y}}_k,{\tilde{\mathbf{X}}}_k,\tilde{\mathit{\beta }})$ provides information for approximating the maximum likelihood, which is a key component for information-based criteria, such as the Akaike information criterion (AIC, [54,55]) and Bayesian information criterion (BIC, [55]).

Recall that the log-likelihood function $l(\mathit{\beta };\mathbf{y},\mathbf{X})$ through the full dataset $\mathbf{y}=(y_1,\dots ,$ $y_N{)}^T$ and $\mathbf{X}={({\mathbf{X}}_1,\dots ,{\mathbf{X}}_N)}^T$ attains its maximum at the MLE $\widehat{\mathit{\beta }}$. We denote the log-likelihood approximated through the RASMR representative quadruple set $\{(n_k,{\tilde{y}}_k,{\tilde{\mathbf{X}}}_k,\tilde{\mathit{\beta }}),$ $k=1,\cdots ,K\}$ by $l(\tilde{\mathit{\beta }};\tilde{\mathbf{y}},\tilde{\mathbf{X}})$, where $\tilde{\mathbf{y}}={({\tilde{y}}_1,\dots ,{\tilde{y}}_K)}^T$, and $\tilde{\mathbf{X}}={({\tilde{\mathbf{X}}}_1,\dots ,{\tilde{\mathbf{X}}}_K)}^T$.

Denoting $\tilde{\Delta }={max}_k{max}_{i\in I_k}\Vert {\mathbf{X}}_i-{\tilde{\mathbf{X}}}_k\Vert$, the maximum Euclidean distance of predictor vectors away from their corresponding predictor representatives across all data blocks, the following theorem guarantees the consistency of the estimated log-likelihood $l(\tilde{\mathit{\beta }};\mathbf{y},\mathbf{X})$, as $\tilde{\Delta }$ goes to zero.

Theorem 1.

Suppose the log-likelihood function $l(\mathit{\beta };\mathbf{y},\mathbf{X})$ is twice differentiable and strictly concave on a compact set $C\subset {\mathcal{R}}^P$, with its maximum not located on the boundary of C. Suppose $\tilde{\mathit{\beta }}$ is the estimate obtained by Algorithms 1 or 2, and ${\tilde{y}}_k=n_k^{-1}{\sum }_{i\in I_k}{y}_i+O\left(\tilde{\Delta }\right)$ for each k, then $l(\tilde{\mathit{\beta }};\mathbf{y},\mathbf{X})-l(\widehat{\mathit{\beta }};\mathbf{y},\mathbf{X})=O\left({\tilde{\Delta }}^{1/2}\right)$ as $\tilde{\Delta }\to 0$.

The proof of Theorem 1, as well as all other proofs in this paper, has been relegated to Section S6 of the Supplementary Materials.

We denote the AIC and BIC values calculated based on $l(\tilde{\mathit{\beta }};\mathbf{y},\mathbf{X})$ as $\widetilde{AIC}$ and $\widetilde{BIC}$, respectively. As a direct conclusion of Theorem 1, we have the following corollary:

Corollary 1.

For a generalized linear model under the same technical conditions of Theorem 1, the values of $\widetilde{AIC}$ and $\widetilde{BIC}$ converge to the original values of AIC and BIC, respectively, as $\tilde{\Delta }\to 0$. Moreover, $\widetilde{AIC}=AIC+O\left({\tilde{\Delta }}^{1/2}\right)$, and $\widetilde{BIC}=BIC+O\left({\tilde{\Delta }}^{1/2}\right)$.

The computation of $l(\tilde{\mathit{\beta }};\tilde{\mathbf{y}},\tilde{\mathbf{X}})$ only requires the K representatives rather than the whole dataset. Its accuracy not only depends on the estimate $\tilde{\mathit{\beta }}$ of the parameters but also on the representatives $\tilde{\mathbf{y}}$ and $\tilde{\mathbf{X}}$. The next theorem confirms that the estimated log-likelihood $l(\tilde{\mathit{\beta }};\tilde{\mathbf{y}},\tilde{\mathbf{X}})$ based on RASMR algorithms is consistent as $\tilde{\Delta }$ goes to zero.

Theorem 2.

Under the same technical conditions of Theorem 1, $l(\tilde{\mathit{\beta }};\tilde{\mathbf{y}},\tilde{\mathbf{X}})-l(\widehat{\mathit{\beta }};\mathbf{y},\mathbf{X})=O\left({\tilde{\Delta }}^{1/2}\right)$, as $\tilde{\Delta }\to 0$.

We denote the AIC and BIC calculated with $l(\tilde{\mathit{\beta }};\tilde{\mathbf{y}},\tilde{\mathbf{X}})$ as $\widetilde{\tilde{AIC}}$ and $\widetilde{\tilde{BIC}}$, respectively. As a direct conclusion of Theorem 2, we have the following corollary:

Corollary 2.

For a generalized linear model under the same technical conditions of Theorem 1, the values of $\widetilde{\tilde{AIC}}$ and $\widetilde{\tilde{BIC}}$ converge to the original values of AIC and BIC, respectively, as $\tilde{\Delta }\to 0$. Moreover, $\widetilde{\tilde{AIC}}=AIC+O\left({\tilde{\Delta }}^{1/2}\right)$, and $\widetilde{\tilde{BIC}}=BIC+O\left({\tilde{\Delta }}^{1/2}\right)$.

Corollary 2 ensures that $\widetilde{\tilde{AIC}}$ and $\widetilde{\tilde{BIC}}$ can be used to approximate the original value of $AIC$ and $BIC$, and the error vanishes along with the maximum size of data blocks. Both SMR and RASMR algorithms can be used to generate $\widetilde{\tilde{AIC}}$ and $\widetilde{\tilde{BIC}}$, but numerical experiments overwhelmingly favor RASMR as it improves both the parameter estimates and the representatives.

3.2. Link Function Selection

For GLM (1), the link function g plays a key role in modeling the relationship between a linear combination of predictors and the expectation of the response variable. It is critical for data analysis to choose the most appropriate link function for a given dataset. An information-based approach, such as $AIC$ or $BIC$, requires precise estimation of the maximum likelihood. Theorem 2 and Corollary 2 provide theoretical justifications for using $\widetilde{\tilde{AIC}}$ or $\widetilde{\tilde{BIC}}$ to select the most appropriate link function.

In our simulation studies (see Section S3 in the Supplementary Materials), SMR fails to choose the logit link function from its competitors (see Table S1), due to the low quality of its representatives. The RASMR with the delta ratio split (Algorithm 2) outperforms SMR with representatives of better quality, making the corresponding $\widetilde{\tilde{AIC}}$ and $\widetilde{\tilde{BIC}}$ more reliable (see Table S2).

3.3. Variable Selection

Variable selection is an essential step in statistical data analysis, including subset selection and stepwise selection (see, for example, [55]).

RASMR can be directly extended for subset selection in big data variable selection problems with a moderate number p of predictors. In the steps involving model fitting and evaluating information criteria ($AIC$ or $BIC$), RASMR can be implemented readily on the dataset to draw accurate results (see Section S5.1 in the Supplementary Materials).

To balance the processing time and performance of variable selection, we introduce a quick variable screening process using the mean representatives (MR, see Section S4 in the Supplementary Materials). The information criteria $\overline{AIC}$ and $\overline{BIC}$ based on MR perform sufficiently well (see Table S1).

For forward stepwise variable selection, we recommend using RASMR to perform a finer selection until the stopping criterion is met (see Section S4 in the Supplementary Materials).

3.4. Cross-Validation

Cross-validation is a widely used data-driven technique for evaluating model performance (see, for example, [55]). For big data analysis under GLMs, we extend RASMR for V-fold cross-validation (VFCV) in the big data analysis as follows:

• Data $\{({\mathbf{X}}_i,y_i),i=1,\dots ,N\}$ are given with a partition $I_1,\dots ,I_K$ of $\{1,\dots ,N\}$.
• A random partition $A_1,\dots ,A_V$ of $\{1,\dots ,N\}$ is given for V-fold cross-validation.
• For $j=1,\dots ,V$, fit the target model on the training set $\{({\mathbf{X}}_i,y_i),i\notin A_j\}$ using RASMR with blocks $I_1\setminus A_j,\dots ,I_K\setminus A_j$ after removing empty ones, and then calculate the aggregated prediction errors ${\widehat{\mathcal{R}}}_j^{VFCV}$ when applying the fitted model on the testing set $\{({\mathbf{X}}_i,y_i),i\in A_j\}$.
• Report $N^{-1}{\sum }_{j=1}^V{\widehat{\mathcal{R}}}_j^{VFCV}$ as the estimated average predictor error.

Note that we do not need to run the same number of iterations on each model fitting. Since RASMR is an iterative approach that can benefit from good initial values of parameter estimates, the estimate from the previous model fitting can be used as the initial values for the next model fitting procedure. Due to this reason, we may set a larger number of iterations for the first RASMR model fitting and then use a smaller number of iterations for later model fitting (see Section S5.2 in the Supplementary Materials).

4. Simulation Studies and Numerical Justifications

4.1. Simulation Setup and Evaluation Method

To justify the performance improvement of the RASMR algorithm on data of massive sizes, we conduct extensive simulations on a main-effects GLM as follows:

$$E\left(Y_i\right)={\mu }_i\mathrm{and}{\eta }_i=g\left({\mu }_i\right)={\mathbf{X}}_i^T\mathit{\beta }={\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_d{x}_{id},$$

(8)

where $i=1,\dots ,N$. Following [43,56], we choose $d=7$, ${\beta }_0=0$, and ${\beta }_1=\cdots ={\beta }_7=0.5$. The feature variables ${\mathbf{x}}_i={(x_{i1},x_{i2},\cdots ,x_{id})}^T$ are randomly generated from one of seven different distributions, namely mzNormal, nzNormal, ueNormal, mixNormal, $T_3$, EXP, and BETA (see Section S5 in the Supplementary Materials for more details).

For Bernoulli regression models (see Example 1 in Section 4.2), we consider four commonly used link functions, namely logit, cloglog, probit, and cauchit. We first simulate $N={10}^6$ data points from each of the seven distributions, respectively, under the logit link. We then obtain the data blocks $I_1,\dots ,I_K$ using K-means clustering with $K=1000$. By assuming each of the four link functions, we obtain the corresponding parameter estimates using the MR, SMR, and RASMR algorithms, respectively. We then compare the estimate ${\tilde{\beta }}_i$ with the full data estimate ${\widehat{\beta }}_i$ under the corresponding link function. The performance of estimation is measured by the root mean square error (RMSE, that is, ${[{\sum }_{i=1}^7{({\tilde{\beta }}_i-{\widehat{\beta }}_i)}^2/7]}^{1/2}$). A smaller RMSE indicates a better approach.

Similarly, we consider a Poisson regression model with a log link (see Example 2 in Section 4.2) and a Gamma regression model with a reciprocal link (see Example 3 in Section 4.2). For each example in Section 4.2 and Section 4.3, the corresponding simulation is repeated 100 times. The summarized results, including the average RMSE and the sample standard deviation (std in parenthesis) of 100 RMSEs, are listed in the corresponding tables.

Note that the K-means clustering algorithm performed for each simulation study is used for illustration purposes. It generates a partition with blocks whose data points are homogeneous. Other clustering methods, such as k-medoids (see, for example, [57] and references therein), may be used for the same purpose as well. Our goal here is to evaluate the performance of MR, SMR, and RASMR given the same partition, not to compare different clustering algorithms.

For SMR, each simulation study recommends a number of iterations, T, set to 3, as suggested by [43], as improvement beyond the third iteration is often negligible (see Figure S12 in the Supplementary Materials). In contrast, RASMR typically continues to improve its parameter estimates as T increases, provided the number of data blocks is sufficiently large (see Figure S12). Overall, we recommend setting T to 10 for RASMR algorithms to balance estimation accuracy with computational time. It should be noted that RASMR still outperforms SMR even with the same number of iterations (see Table 6).

4.2. Performance of RASMR, Algorithm 1

In this section, we evaluate the performance of Algorithm 1 on Bernoulli, Poisson, and gamma regression models, respectively.

Example 1.

In this example, we consider Bernoulli regression models (see (8)) for binary classifications, with one of four commonly used link functions, namely logit, cloglog, probit, and cauchit. The goal is to check not only the estimation performance of RASMR when the link function is correctly specified but also its performance when the link function is misspecified, which is crucial when selecting the link function.

The results of the corresponding simulation study as described in Section 4.1 are summarized in

Table 2. Average (std) of RMSEs (${10}^{-3}$) over 100 simulations for binary classification.

Simulation	Binary Classification, K-Means ($\mathit{K}=1000$), True Link Function = Logit
Representatives	MR				Original SMR				RASMR
Setup	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit
mzNormal	17.9	9.2	9.2	33.4	1.8	1.6	1.0	2.4	$3.7\times {10}^{-5}$	$2.3\times {10}^{-1}$	$6.4\times {10}^{-5}$	$1.0\times {10}^{-2}$
	(0.3)	(0.3)	(0.2)	(1.2)	(0.6)	(0.5)	(0.3)	(2.4)	$(1.0\times {10}^{-5})$	$(1.8\times {10}^{-1})$	$(2.6\times {10}^{-5})$	$(4.2\times {10}^{-3})$
nzNormal	14.3	5.0	7.0	135	4.0	0.8	1.5	51.6	$1.5\times {10}^{-3}$	$2.5\times {10}^{-2}$	$2.5\times {10}^{-4}$	69.6
	(0.7)	(0.3)	(0.3)	(45.8)	(1.1)	(0.4)	(0.5)	(34.6)	$(5.0\times {10}^{-4})$	$(3.0\times {10}^{-2})$	$(7.4\times {10}^{-5})$	(20.0)
ueNormal	211	114	110	455	3.3	11	2.1	19.4	$7.2\times {10}^{-5}$	$4.3$	$4.7\times {10}^{-4}$	$10.4$
	(1.4)	(1.5)	(0.7)	(5.2)	(1.4)	(3.4)	(1.1)	(10.1)	$(3.3\times {10}^{-6})$	$\left(1.2\right)$	$(9.9\times {10}^{-5})$	$\left(10.2\right)$
mixNormal	17.5	8.6	8.6	48.1	3.0	1.8	1.2	3.2	$3.6\times {10}^{-4}$	$5.0\times {10}^{-1}$	$1.6\times {10}^{-4}$	$1.9\times {10}^{-1}$
	(0.4)	(0.6)	(0.2)	(3.1)	(0.9)	(0.7)	(0.3)	(2.5)	$(1.6\times {10}^{-4})$	$(1.8\times {10}^{-1})$	$(1.8\times {10}^{-4})$	$(5.5\times {10}^{-2})$
$T_3$	12.2	10.1	7.8	10.0	10.7	15.0	6.7	8.6	$6.3\times {10}^{-2}$	$1.2\times {10}^{-1}$	$5.2\times {10}^{-2}$	$2.8\times {10}^{-2}$
	(3.1)	(2.6)	(1.9)	(2.7)	(3.1)	(39.5)	(1.9)	(2.7)	$(3.5\times {10}^{-2})$	$(7.2\times {10}^{-2})$	$(5.2\times {10}^{-2})$	$(2.8\times {10}^{-2})$
EXP	12.4	3.9	6.2	10.4	5.8	1.4	2.8	4.5	$1.4\times {10}^{-6}$	$3.3\times {10}^{-4}$	$5.3\times {10}^{-6}$	$9.6\times {10}^{-3}$
	(0.9)	(0.5)	(0.5)	(1.4)	(1.0)	(0.3)	(0.5)	(1.4)	$(3.3\times {10}^{-7})$	$(8.1\times {10}^{-5})$	$(6.9\times {10}^{-7})$	$(4.4\times {10}^{-3})$
BETA	3.1	1.3	1.7	7.6	2.0	0.9	1.0	11.3	$9.2\times {10}^{-7}$	$4.7\times {10}^{-5}$	$2.7\times {10}^{-5}$	$5.6\times {10}^{-3}$
	(0.8)	(0.4)	(0.5)	(1.7)	(0.7)	(0.3)	(0.3)	(2.3)	$(6.2\times {10}^{-7})$	$(2.4\times {10}^{-5})$	$(2.6\times {10}^{-6})$	$(6.7\times {10}^{-3})$

. In almost all scenarios, the RMSEs of RASMR are significantly lower than those of the mean representative (MR) approach or the original SMR approach, which implies that the RASMR algorithm (Algorithm 1) outperforms MR and the original SMR algorithm in this example. The estimates based on RASMR are not only much more accurate than the estimates using MR or SMR but also more robust.

Nevertheless, for the case with the cauchit link and nzNormal distribution, the RASMR algorithm performs as poorly as other methods. In Section 4.3, we will show further improvement by using Algorithm 2.

Example 2.

In this example, we consider Poisson regression with its canonical link function, namely, log link. The corresponding simulation is repeated for 100 times and the results are summarized in

Table 3. Average (std) of RMSEs (${10}^{-3}$) over 100 simulations for the Poisson regression.

Simulation	Poisson Regression, K-Means ($\mathit{K}=1000$)
Representatives	MR	Percentage of NAs	Original SMR	Percentage of NAs	RASMR	Percentage of NAs
mzNormal	37.5 (12.9)	0%	13.2 (10.2)	0%	2.0 (0.5)	0%
nzNormal	37.5 (12.9)	0%	23.6 (15.7)	3%	0.2 $(3.6\times {10}^{-2})$	0%
ueNormal	82.5 (25.9)	0%	9.5 (7.2)	69%	$8.5\times {10}^{-2}$$(2.3\times {10}^{-2})$	0%
mixNormal	49.9 (20.7)	0%	29.5 (17.7)	1%	0.7 (0.1)	0%
$T_3$	298 (458)	0%	97.2 (144)	10%	31.3 (79.1)	0%
EXP	31.2 (0.7)	0%	6.2 (1.0)	0%	$1.9\times {10}^{-4}$$(5.8\times {10}^{-5})$	0%
BETA	0.5 (0.2)	0%	2.3 (0.3)	0%	$1.2\times {10}^{-7}$$(1.8\times {10}^{-9})$	0%

.

From Table 3 we can see that (1) In terms of the average RMSE (the smaller, the better), the original SMR outperforms MR, and RASMR further improves the accuracy of estimates impressively, although not as accurate as the full data estimate in this case (see Table 3 in [43]); (2) In terms of the percentage of NAs (the smaller, the more robust) out of 100 simulations, both MR and RASMR achieve $0\%$ and outperform the original SMR.

Example 3.

In the example, we consider Gamma regression with its canonical link function, namely, reciprocal link. Since Gamma regression with reciprocal is relatively fragile on the distribution of features, we only consider features generated from Beta distribution. Again, the data clusters are generated from K-means clustering with $K=1000$. Similar to Examples 1 and 2, we obtain the RMSE of the MR, SMR, and RASMR estimates from the full data estimate and summarize the results in

Table 4. Average (std) of RMSEs (${10}^{-3}$) over 100 Simulations for Gamma Regression.

Simulation	Gamma Regression, K-Means ($\mathit{K}=1000$)
Representatives	MR	Percentage of NAs	Original SMR	Percentage of NAs	RASMR	Percentage of NAs
Beta	7.4 (0.7)	0%	5.8 (1.2)	11%	$2.0\times {10}^{-2}$ (0.2)	0%

. The pattern of the results is similar to the ones in Table 3, which confirms the significant improvement of RASMR against both MR and SMR under Gamma regression models.

In [43], extensive comparisons between SMR and other big data approaches have been made in terms of computing speed. According to [43], given a partition of data, SMR is noticeably faster than a divide-and-conquer (DC) approach proposed by [5], but slower than the subsampling approach under A-optimality [56] for logistic regression. Using their results of computing time comparison as a reference, we compare the computing speed of RASMR with SMR for binary classification with four link functions: logit, cloglog, probit, and cauchit. The covariates are generated from the seven distributions described previously and the computation is conducted on a PC running Windows 10 Home (Version 2004) with 2.80 GHz 4-core Intel i7-7700HQ and 16 GB of memory. The simulation is repeated 100 times and the average CPU time for one iteration of SMR and RASMR is recorded in

Table 5. Average CPU time (seconds) over 100 simulations for binary classification.

Simulation	Binary Classification, K-Means (K = 1000), True Link Function = Logit
Representatives	SMR ($\mathit{T}$ = 3)				RASMR ($\mathit{T}$ = 3)
Setup	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit
mzNormal	6.04	8.95	11.09	6.51	7.32	9.75	12.95	7.63
nzNormal	5.99	8.90	12.82	6.22	6.60	9.82	13.65	6.70
ueNormal	6.75	10.33	13.55	7.14	7.68	10.38	15.82	7.79
mixNormal	5.90	8.52	11.26	6.22	6.63	9.52	12.64	7.02
$T_3$	6.32	8.43	8.10	6.36	6.77	9.06	8.50	7.00
EXP	5.66	7.98	10.51	6.24	6.55	9.20	12.11	6.90
BETA	5.66	7.96	10.73	6.36	6.51	9.10	12.57	7.01

.

According to Table 5, RASMR is slightly slower than SMR. It is because additional splits are made on a proportion of data blocks, which result in a larger number of data blocks in RASMR, given the same partition of data. Nevertheless, based on our experience, most of the RASMR split requirements are not fulfilled, thus the resulting final number of data blocks of RASMR is only slightly larger than that of SMR.

Similar to [43], we calculate the corresponding RMSE from the full data estimate ${[{\sum }_{i=1}^7{({\tilde{\beta }}_i-{\widehat{\beta }}_i)}^2/7]}^{1/2}$ and the RMSE from the true parameter values ${[{\sum }_{i=1}^7{({\tilde{\beta }}_i-{\beta }_i)}^2/7]}^{1/2}$, respectively. The simulation follows the same settings as in Example 1, and the comparisons are among MR, SMR with $T=3$, RASMR with $T=3$, and DC with $K=1000$. The results are listed in

Table 6. Average (std) of RMSEs (${10}^{-3}$) over 100 simulations from full data estimates or true parameter value.

Simulation		Binary Classification, K-Means (K = 1000), True Link Function = Logit
Benchmark		RMSE from Full Data Estimate				RMSE from True Parameter Value
	Methods	Representative Approaches			DC	Full Data	Representative Approaches			DC
Setup		MR	SMR ($\mathit{T}$ = 3)	RASMR ($\mathit{T}$ = 3)			MR	SMR ($\mathit{T}$ = 3)	RASMR ($\mathit{T}$ = 3)
mzNormal		17.98(0.31)	1.92(0.62)	0.054(0.018)	6.93(0.12)	3.72(1.08)	18.40(1.02)	4.17(1.16)	3.72(1.08)	7.90(1.04)
nzNormal		14.30(0.70)	4.53(1.18)	0.29(0.085)	20.20(0.35)	7.23(2.05)	16.04(1.74)	8.58(2.02)	7.24(2.07)	21.49(1.64)
ueNormal		211.17(1.44)	3.83(1.57)	0.72(0.018)	13.12(0.26)	2.11(0.82)	211.17(1.35)	4.41(1.77)	2.22(0.84)	13.24(1.24)
mixNormal		17.37(0.33)	2.79(0.83)	0.14(0.043)	11.20(0.20)	4.96(1.33)	17.94(1.05)	5.91(1.48)	4.96(1.33)	12.09(1.15)
$T_3$		12.23(3.13)	11.3(3.23)	1.89(0.74)	12.06(0.34)	16.00(4.43)	20.51(5.44)	20.03(5.54)	16.03(4.51)	19.63(3.76)
EXP		12.4(9.15)	6.58(0.82)	0.014(0.0013)	16.88(0.31)	6.18(1.67)	14.5(2.18)	9.25(2.02)	6.18(1.67)	18.25(2.30)
BETA		3.03(0.80)	2.34 0.68)	0.00031(0.000090)	5.92(0.20)	7.49(2.38)	7.89(2.30)	7.73(2.34)	7.49(2.38)	9.31(2.54)

.

According to Table 6, with the same number of iterations ($T=3$), RASMR still outperforms MR, SMR, and DC. When the comparison is made with respect to the true parameter values, RASMR performs almost as well as the full data estimate with negligible differences. Recall that the data contain the entire information for estimating the parameters and the RASMR approach achieves its theoretical extreme in terms of estimation accuracy in this case.

4.3. Performance of RASMR with the Delta Ratio Split, Algorithm 2

In Section 2.4, we develop Algorithm 2 for estimating the maximum likelihood better, especially for non-Gaussian covariates or predictors. To justify the improvement of Algorithm 2 against RASMR, in this section, we employ the same simulation setting as described in Section 4.1 (see Example 1). For illustration purposes, the threshold ${\delta }_0$ for delta ratio split is chosen to be $0.05$, $0.1$, $0.5$, and 1, respectively. We fit the logistic regression model using SMR, RASMR, and RASMR with the delta ratio split and various thresholds. The estimated maximum log-likelihood approaches the full data value as the threshold of the delta ratio decreases, except in the $T_3$ case, where the improvement is negligible (see Figure S11 in the Supplementary Materials). One reason is that in the case of $T_3$, the delta ratios are much smaller than those in other cases (see Figure S1 in the Supplementary Materials), and even $0.05$ does not capture a good amount of blocks.

To evaluate how the delta ratio split impacts the accuracy of parameter estimates, we choose $K=100,500$, and 1000 for K-means clustering, respectively. With moderately large K, such as $K=500$ or 1000, RASMR with a delta ratio split performs roughly the same as the original RASMR (see Figure S12 in the Supplementary Materials), and both RASMR and RASMR with a delta ratio split outperform SMR. Nevertheless, when K is as small as 100, the RASMR algorithms may fail for some data structures, such as nzNormal (where the data have imbalanced responses), $T_3$ (where the data distribute with heavy tails), or mixNormal (where the data come from a mixture of two distributions).

In

Table 7. Average (std) of RMSEs (${10}^{-3}$) over 100 simulations using RASMR with the delta ratio split.

Simulation	K-Means (K = 1000), True Link Function = Logit
Approach	RASMR with the Delta Ratio Split, Threshold = $\mathbf{0.1}$
Setup	Logit	Cloglog	Probit	Cauchit
mzNormal	$4.0\times {10}^{-5}$	$3.2\times {10}^{-2}$	$6.1\times {10}^{-5}$	$8.2\times {10}^{-3}$
	$(9.9\times {10}^{-6})$	$(3.7\times {10}^{-2})$	$(2.4\times {10}^{-5})$	$(3.2\times {10}^{-3})$
nzNormal	$1.6\times {10}^{-3}$	$1.5\times {10}^{-2}$	$2.7\times {10}^{-4}$	$8.8$
	$(5.6\times {10}^{-4})$	$(1.4\times {10}^{-2})$	$(8.5\times {10}^{-5})$	$\left(4.5\right)$
ueNormal	$2.0\times {10}^{-4}$	$9.9\times {10}^{-2}$	$1.2\times {10}^{-3}$	$2.3\times {10}^{-2}$
	$3.4\times {10}^{-3}$	$9.6\times {10}^{-1}$	$2.4\times {10}^{-3}$	$3.3\times {10}^{-2}$
mixNormal	$3.8\times {10}^{-4}$	$1.9\times {10}^{-1}$	$1.8\times {10}^{-4}$	$7.3\times {10}^{-2}$
	$(1.5\times {10}^{-4})$	$(1.4\times {10}^{-1})$	$(1.8\times {10}^{-4})$	$(3.0\times {10}^{-2})$
$T_3$	$6.3\times {10}^{-2}$	$1.0\times {10}^{-1}$	$5.2\times {10}^{-2}$	$2.7\times {10}^{-2}$
	$3.4(\times {10}^{-2})$	$(5.6\times {10}^{-2})$	$(2.5\times {10}^{-2})$	$(1.5\times {10}^{-2})$
EXP	$2.3\times {10}^{-6}$	$4.2\times {10}^{-4}$	$5.8\times {10}^{-6}$	$1.4\times {10}^{-2}$
	$(1.8\times {10}^{-6})$	$(9.0\times {10}^{-5})$	$(7.7\times {10}^{-7})$	$(9.7\times {10}^{-3})$
BETA	$4.6\times {10}^{-6}$	$7.3\times {10}^{-5}$	$2.6\times {10}^{-5}$	$1.0\times {10}^{-2}$
	$(3.9\times {10}^{-7})$	$(3.4\times {10}^{-5})$	$(3.4\times {10}^{-6})$	$(3.8\times {10}^{-3})$

, we list the performance of RASMR with the delta ratio split and ${\delta }_0=0.1$ when the link function is misspecified. Comparing Table 7 with Table 2, we can see that RASMR with the delta ratio split (Algorithm 2) is, in general, better than or comparable with Algorithm 1 in terms of estimation accuracy. The improvements are significant in the unbalanced or heavy-tailed cases, such as nzNormal with cauchit link, and ueNormal with cloglog or cauchit links (see Section S5 in the Supplementary Materials for more details).

5. Real Data Analysis

In this section, we use a real example—the airline on-time performance data—collected from the Bureau of Transportation Statistics (https://www.transtats.bts.gov/, accessed on 31 August 2018), to illustrate how RASMR works for analyzing real data with a massive size.

5.1. The Airline On-Time Performance Data

The airline on-time performance data contain detailed information on US domestic flights since October 1987. The original data are saved in individual files labeled by month. For illustration purposes, in this study, we use data from between October 1987 and August 2018, 371 files/months in total, which contain 182,751,882 flights. After removing the records with missing or invalid inputs, N = 179,528,198 is actually considered in this study.

Following [43], we formulate the data into a binary classification problem aiming to model the flight delay status, a binary response ArrDelayLabel (arrival delay label), which is generated by cutting the continuous variable ARRIVAL DELAY at the 15-min point. For simplicity, the departure time block DepTimeBlk, originally a factor of 17 levels, is regrouped into a factor with 4 levels (1 for 12:00 a.m.–05:59 a.m., 2 for 06:00 a.m.–11:59 a.m., 3 for 12:00 p.m.–05:59 p.m., and 4 for 06:00 p.m.–11:59 p.m.). Moreover, three categorical variables (QUARTER, DayOfWeek, and DepTimeBlk) and three continuous variables (DISTANCE, DepDelay, and CRSTimeElapsed) are used as illustrations to explain the status of flight arrival delays (see Table S8 in the Supplementary Materials for a list of the variables).

5.2. Model Selection

In this section, we select the most appropriate model for the airline data, which consists of 371 months and 179,528,198 flight records.

As described in Section 3.2, we first perform the link function selection using $\widetilde{\tilde{AIC}}$ and five-fold cross-validation based on RASMR with the delta ratio split (Algorithm 2). In this step, all variables listed in Table S8 are included, and the candidate link functions include logit, cloglog, probit, and cauchit.

According to the link function selection results listed in

Table 8. Link function selection using $\widetilde{\tilde{AIC}}$ and 5-fold cross-validation with RASMR.

Link Function	Logit	Cloglog	Probit	Cauchit
$\widetilde{\tilde{AIC}}$	90284022	97172503	95844750	113052242
5-fold CV with Cross-entropy Loss	1.38441	1.7491	1.8228	2.1535

, we select the logit link for our model (see Figure S14 in the Supplementary Materials for the fitted model based on “glm2” function in [58] and the RASMR representative data).

As described in Section 3.3, we further perform the variable selection for the airline data by implementing the stepwise selection strategy with an initial screening based on MR and fine selection using RASMR (see Section S4 in the Supplementary Materials). The selection criterion is $\widetilde{\tilde{AIC}}$ again. When using “glm2”, we estimate the dispersion parameter by Pearson’s Chi-square statistic divided by its degree of freedom (see Figure S15 in the Supplementary Materials). If n, N, and p represent the number of representatives, the full data size, and the number of parameters in the model, respectively, the degrees of freedom for Pearson’s Chi-square test are $n-p$ for RASMR and $N-p$ for the full data. To estimate the dispersion parameter of full data, the RASMR dispersion parameter estimate needs to be adjusted by the degrees of freedom, i.e., $\widehat{\varphi }={\widehat{\varphi }}_{RASMR}·(n-p)/(N-p)$, where $\widehat{\varphi }$ is the dispersion parameter estimated from the full data and ${\widehat{\varphi }}_{RASMR}$ is the estimate from the RASMR representative dataset. After such an adjustment, the estimated dispersion parameters are 3291.14 and 3341.67 for the model with all variables and the model with the selected variables, respectively, which indicates strong over-dispersion.

In conclusion, the GLM after variable selection contains only one variable, namely, DepDelay (departure delay in minutes), which is the most informative variable for predicting the arrival delay of a flight.

5.3. Comparison Analysis

For performance evaluation and comparison purposes, we create an oracle model by fitting a logistic regression model using data from between March 2012 and February 2017, 60 months in total, utilizing all explanatory variables listed in Table S8. Then, the oracle response of ArrDelayLabel for the entire dataset is generated accordingly based on the oracle model.

Assuming that we do not know the true link (i.e., logit), we fit the Bernoulli model with one of four different links (namely, logit, cloglog, probit, and cauchit) using four approaches: (1) iterative reweighted least square (IRLS) on the combined data; (2) the mean representative (MR) approach; (3) the score-matching representative (SMR) approach; and (4) RASMR with the delta ratio split and exponential learning rate decay.

The data blocks for MR, SMR, and RASMR are created by the following strategy. The original airline data are saved in monthly files, which are regarded as natural data blocks. Before applying MR, SMR, or RASMR, we further split each of the monthly files according to the distinct values of three categorical variables, namely, QUARTER, DayOfWeek, and DepTimeBlk. Then, each monthly file is divided into $4\times 7\times 4=112$ sub-blocks. We then split each sub-block further according to the three continuous variables via one of the following two schemes. The first scheme is a correlation-based quantile split. That is, we split each data block at the 25%, 50%, and 75% quantiles of the three continuous variables as the cutting points, and obtain $4\times 4\times 4=64$ sub-blocks. The second scheme is to apply K-means clustering to the three continuous variables with $K=64$. As a result, each monthly file is divided into $112\times 64=7168$ sub-blocks (see Section S1 in the Supplementary Materials for more discussion).

We compare MR, SMR, and RASMR with four different data sizes, namely, 60 months, 120 months, 240 months, and 371 months. The IRLS estimates on the combined data are referred to as the full data estimate, but only available for 60-month and 120-month datasets. As a result, in

Table 9. Average (std) of RMSEs (${10}^{-3}$) over 10 simulations of the airline on-time performance data using K-means clustering.

Simulation	Binary Classification, K-Means (K = 64), True Link Function = Logit
Representatives	MR				SMR				RASMR
Setup	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit
60 months	28.4	44.3	12.1	135.7	29.3	212.3 (NA removed)	12.6	115.5	1.6	3.0	0.3	34.3
	$(4.5\times {10}^{-4})$	$(4.6\times {10}^{-5})$	$(1.2\times {10}^{-4})$	$(7.0\times {10}^{-4})$	$(2.2\times {10}^{-3})$	8.7 (NA removed)	$(7.6\times {10}^{-4})$	$(1.6\times {10}^{-2})$	$(1.0\times {10}^{-5})$	$(8.5\times {10}^{-6})$	$(2.0\times {10}^{-6})$	$(2.5\times {10}^{-5})$
120 months	28.3	44.3	12.2	135.4	29.3	212.1 (NA removed)	12.6	115.3	1.6	3.0	0.3	34.3
	$(4.5\times {10}^{-4})$	$(4.5\times {10}^{-5})$	$(1.2\times {10}^{-4})$	$(7.0\times {10}^{-4})$	$(2.2\times {10}^{-3})$	8.6 (NA removed)	$(7.6\times {10}^{-4})$	$(1.6\times {10}^{-2})$	$(1.0\times {10}^{-5})$	$(8.5\times {10}^{-6})$	$(2.0\times {10}^{-6})$	$(2.5\times {10}^{-5})$
240 months	24.7	41.8	9.0	111.4	33.2	219.8 (NA removed)	12.7	142.0	1.6	3.2	0.3	31.1
	$(3.3\times {10}^{-4})$	$(2.6\times {10}^{-5})$	$(8.9\times {10}^{-5})$	$(7.4\times {10}^{-4})$	$(4.3\times {10}^{-3})$	9.4 (NA removed)	$(7.6\times {10}^{-4})$	$(2.2\times {10}^{-2})$	$(9.7\times {10}^{-6})$	$(1.0\times {10}^{-5})$	$(2.8\times {10}^{-6})$	$(1.7\times {10}^{-5})$
371 months	24.2	41.4	9.0	111.2	35.9	216.1 (NA removed)	12.7	140.3	1.6	3.2	0.3	31.2
	$(3.1\times {10}^{-4})$	$(2.6\times {10}^{-5})$	$(8.9\times {10}^{-5})$	$(7.4\times {10}^{-4})$	$(6.2\times {10}^{-3})$	9.4 (NA removed)	$(7.6\times {10}^{-4})$	$(2.0\times {10}^{-2})$	$(9.7\times {10}^{-6})$	$(1.0\times {10}^{-5})$	$(2.8\times {10}^{-6})$	$(1.6\times {10}^{-5})$

Notes: RMSEs of 60 months and 120 months are calculated with respect to the IRLS estimates; RMSEs of 240 months and 371 months are calculated with respect to the oracle parameter values.

and

Table 10. Average (std) of RMSEs (${10}^{-3}$) over 10 simulations of the airline on-time performance data using correlation-based quantile split.

Simulation	Binary Classification, Correlation-Based Quantile Split, True Link Function = Logit
Representatives	MR				SMR				RASMR
Setup	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit	Logit	Cloglog	Probit	Cauchit
60 months	62.7	182.8	30.1	376.4	150.1	239.2	49.1	384.0	21.3	20.6	2.4	37.2
	$(1.9\times {10}^{-3})$	$(7.9\times {10}^{-2})$	$(1.5\times {10}^{-2})$	$(4.2\times {10}^{-2})$	0.7	0.4	$(2.6\times {10}^{-2})$	1.5	$(1.5\times {10}^{-4})$	$(1.1\times {10}^{-3})$	$(1.2\times {10}^{-4})$	$(6.6\times {10}^{-3})$
120 months	62.5	182.7	30.1	374.2	150.0	239.5	49.4	384.7	21.3	20.6	2.4	37.2
	$(1.9\times {10}^{-3})$	$(7.9\times {10}^{-2})$	$(1.5\times {10}^{-2})$	$(4.2\times {10}^{-2})$	0.7	0.4	$(2.6\times {10}^{-2})$	1.5	$(1.5\times {10}^{-4})$	$(1.1\times {10}^{-3})$	$(1.2\times {10}^{-4})$	$(6.6\times {10}^{-3})$
240 months	60.4	180.4	27.7	385.4	147.3	242.1	48.4	373.5	20.2	20.8	2.4	36.9
	$(1.1\times {10}^{-3})$	$(8.1\times {10}^{-2})$	$(1.0\times {10}^{-2})$	$(4.4\times {10}^{-2})$	0.7	0.4	$(2.7\times {10}^{-2})$	1.4	$(1.4\times {10}^{-4})$	$(1.1\times {10}^{-3})$	$(1.2\times {10}^{-4})$	$(6.4\times {10}^{-3})$
371 months	60.3	180.1	27.7	385.2	147.3	242.1	48.4	373.4	20.2	20.84	2.4	36.9
	$(1.1\times {10}^{-3})$	$(8.1\times {10}^{-2})$	$(1.0\times {10}^{-2})$	$(4.4\times {10}^{-2})$	0.7	0.4	$(2.6\times {10}^{-2})$	1.4	$(1.4\times {10}^{-4})$	$(1.1\times {10}^{-3})$	$(1.2\times {10}^{-4})$	$(6.4\times {10}^{-3})$

Notes: RMSEs of 60 months and 120 months are calculated with respect to the IRLS estimates; RMSEs of 240 months and 371 months are calculated with respect to the oracle parameter values.

, we list RMSE of MR, SMR, and RASMR (with the delta ratio split) from the IRLS estimates for 60 months and 120 months, and their RMSE from the oracle parameter values for 240 months and 371 months.

According to Table 9 and Table 10, RASMR produces consistently more accurate estimates, regardless of the sample size or clustering scheme.

To further compare the convergence properties of SMR and RASMR under different link functions, we plot log (RMSE) based on SMR or RASMR with K-means clustering in Figure 2 and quantile split clustering in Figure 3, respectively, against the number of iterations. Both Figure 2 and Figure 3 show that RASMR has a better convergence property than the original SMR in this case, which makes RASMR benefit considerably from the increasing number of iterations.

6. Conclusions

The proposed RASMR approach manages to generate more accurate estimates of parameters due to better-quality representatives after additional model-specific splitting. The estimation accuracy improves along with the increasing number of iterations. For data with a relatively simple structure, typically a small number of iterations can provide accurate enough estimates for further data analysis, including model selection. For data with a more complicated structure or a larger number of parameters, we recommend an increased number of iterations and an exponential learning rate decay to ensure better convergence behavior.

For a practical implementation of RASMR, K-means clustering may consume a considerable amount of time. When a natural partition is not available, we recommend small K-means clustering or correlation-based quantile split for a faster clustering process (see Section S1 in the Supplementary Materials for more details). When the number of covariates or predictors is large, we recommend the correlation-based quantile split strategy due to its fast speed.

When applying RASMR to the model selection, we recommend RASMR with the delta ratio split. A default threshold of delta ratio is 1 since the delta ratios rarely go beyond 1 for typical applications. When the competing models have very close performance, we recommend a smaller threshold, such as $0.1$, to ensure better-quality representatives.

A key step with RASMR is to construct a good set of representatives and find an accurate approximation $\tilde{\mathit{\beta }}$ for the MLE $\widehat{\mathit{\beta }}$ based on the full data. When some collinearity exists among the covariates or predictors, $\widehat{\mathit{\beta }}$ may not be unique and the RASMR algorithm may encounter an identifiability issue. It is worth exploring how to adjust RASMR for variable selection under the presence of collinearity.

The proposed RASMR approach is especially useful for analyzing big data with binary responses, such as intrusion detection for cyber security [59] and fraud detection for insurance companies [60]. Nevertheless, when the responses have three or more categories, the data may be modeled by multinomial logistic models [61,62,63] instead of GLMs. In this case, it is challenging to construct representatives for given data blocks that are more efficient than the corresponding mean representatives.