Efficient Estimation and Validation of Shrinkage Estimators in Big Data Analytics

Abstract

Shrinkage estimators are often used to mitigate the consequences of multicollinearity in linear regression models. Despite the ease with which these techniques can be applied to small- or moderate-size datasets, they encounter significant challenges in the big data domain. Some of these challenges are that the volume of data often exceeds the storage capacity of a single computer and that the time required to obtain results becomes infeasible due to the computational burden of a high volume of data. We propose an algorithm for the efficient model estimation and validation of various well-known shrinkage estimators to be used in scenarios where the volume of the data is large. Our proposed algorithm utilises sufficient statistics that can be computed and updated at the row level, thus minimizing access to the entire dataset. A simulation study, as well as an application on a real-world dataset, illustrates the efficiency of the proposed approach.

1. Introduction

The ability to collect data for which the volume, velocity, and variety exceeds the capacity of standard analytical tools has been made possible by advancements in computer technology employed in practically every field of science and daily life. Access to these data types encourages innovation across a wide range of scientific fields but poses challenges to data storage, computational efficiency, and current statistical and computational methods [1]. Although the high variety and velocity aspects of big data require novel statistical methodology, the data considered in this paper speak to the high-volume aspect of big data. A dataset with a size exceeding $20\%$ of a computer’s RAM is considered large, whereas a dataset with a size exceeding $50\%$ of the computer’s RAM is considered massive since simple calculations would consume all remaining RAM [2]. Access to large datasets enables us to gain a better understanding of the relationship between a response variable and a set of predictor variables, but we often face major issues during the analysis of such data. Multicollinearity is one of these issues [3].

Multicollinearity refers to the existence of linear relationships between the predictor variables in a regression problem. Two major consequences of multicollinearity are unstable regression coefficients due to inflated standard errors and poor out-of-sample prediction performance. The consequences of multicollinearity, such as unstable regression coefficients, are often addressed with shrinkage estimators. Shrinkage estimators restrict the parameter space to avoid parameter explosion due to inflated standard errors. However, these estimators are often nonlinear with respect to the shrinkage parameter, which complicates their optimisation. For an overview of ongoing research relating to shrinkage estimators for multicollinearity and optimal tuning parameter selection, see [3,4].

Despite the ease with which these techniques can be applied to small- or moderate-size datasets, they encounter significant challenges in the big data domain. Some of these challenges are that the volume of the data frequently exceeds the storage capacity of a single computer and that the time required to obtain results becomes infeasible due to the computational burden of a high volume of data. Many strategies for the efficient computation of shrinkage estimators in ultrahigh dimensions exist, but the body of works relating to shrinkage estimators for datasets with many observations is limited. In the context of a large number of observations, Ref. [5] proposed an algorithm that utilises sufficient statistics to obtain the ridge estimator of the multiple linear regression model to address the above-mentioned big data challenges. Their proposal relied on a grid search strategy for finding the optimal shrinkage parameter.

In this paper, we extend the work of [5] by proposing a method for the efficient model estimation of various well-known ridge- and Liu type estimators with closed-from solutions for the optimal shrinkage parameter. Our approach also allows efficient model validation through K-fold cross-validation to evaluate the generalisability of the model. We utilise an array of sufficient statistics and can overcome the memory and computational efficiency barrier since we do not need access to the entire dataset simultaneously to obtain sufficient statistics. Furthermore, the dataset can be removed from memory once the sufficient statistic has been computed, and hence, minimal access to the dataset is required. For illustrative purposes, we consider five shrinkage estimators that aim to address multicollinearity in the linear regression model. These estimators are the ridge estimator [6], the modified ridge type estimator [7], the Liu estimator [8], the modified one-parameter Liu estimator [9] and the Kibria–Lukman estimator [10]. This paper is structured as follows: Section 2 introduces the linear regression model and the shrinkage estimators under consideration. An algorithm for efficient model estimation and validation is proposed in Section 3. Section 4 illustrates the efficiency of our proposed approach through a simulation study and application on a real-world dataset, and lastly, Section 5 concludes the paper.

2. Statistical Methodology

Consider the multiple linear regression model given by

$$\begin{array}{c}\hfill \mathbf{y}=\mathbf{X}\mathit{\beta }+\mathbf{e},\end{array}$$

(1)

with $\mathbf{y}={\left(y_1,\dots ,y_n\right)}^{\top }$ the n-dimensional response variable, $\mathbf{X}={\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)}^{\top }$ the $n\times p$ matrix of nonstochastic explanatory variables, $\mathit{\beta }={\left({\beta }_1,{\beta }_2,\dots ,{\beta }_p\right)}^{\top }$ the p-dimensional vector of regression coefficients, and $\mathbf{e}={\left(e_1,\dots ,e_n\right)}^{\top }\sim N\left(0,{\sigma }^2\mathbf{I}\right)$ the n-dimensional error vector. If rank(${\mathbf{X}}^{\top }\mathbf{X}$) = p, then let $\mathbf{Z}=\mathbf{X}\mathbf{Q}$ and $\alpha ={\mathbf{Q}}^{\top }\mathit{\beta }$, where $\mathbf{Q}$ is the orthogonal matrix whose columns constitute the eigenvectors of ${\mathbf{X}}^{\top }\mathbf{X}$. Thus, the ordinary least squares estimator of $\mathit{\alpha }$ is $\widehat{\mathit{\alpha }}={\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^{\top }\mathbf{y}$. Accordingly, ${\widehat{\sigma }}^2={\mathbf{y}}^{\top }\left[\mathbf{I}-\mathbf{Z}{\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^{\top }\right]\mathbf{y}/\left(n-p\right)$ and $V\left(\widehat{\mathit{\alpha }}\right)={\widehat{\sigma }}^2{\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}$.

Now, consider any estimator of $\mathit{\alpha }$ with the form

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\alpha }}}_g& =& g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right){\mathbf{Z}}^{\top }\mathbf{y}\hfill \end{array}$$

(2)

where $g\in {\mathbf{R}}^{n\times p}\to \mathbf{S}\left(p\right)$, with $\mathbf{S}\left(p\right)$ the space of all positive definite, symmetric, $p\times p$ matrices. The variance of ${\widehat{\mathit{\alpha }}}_g$ is given by

$$\begin{array}{ccc}\hfill V\left({\widehat{\mathit{\alpha }}}_g\right)& =& g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right){\mathbf{Z}}^{\top }COV\left(\mathbf{y}\right)\mathbf{Z}g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)\hfill \\ & =& {\sigma }^{2}g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right).\hfill \end{array}$$

In the presence of multicollinearity, ${\mathbf{Z}}^{\top }\mathbf{Z}$ is ill-conditioned, and as a result, the OLS estimator of $\mathit{\alpha }$ has a large variance. Various shrinkage estimators have been proposed to mitigate the consequences of multicollinearity in the linear regression model. The shrinkage estimators considered in this paper are listed below. These can all be written in terms of the form given in (2), which is used in the algorithm proposed in Section 3.

• The ridge estimator, with $k\ge 0$, proposed by [6] as a possible solution to multicollinearity is given by

  $$\begin{array}{ccc}\hfill {\widehat{\mathit{\alpha }}}_{RRE}\left(k\right)& =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+k{\mathbf{I}}_p\right)}^{-1}{\mathbf{Z}}^{\top }\mathbf{y}.\hfill \end{array}$$
  In this case, $g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)={\left({\mathbf{Z}}^{\top }\mathbf{Z}+k{\mathbf{I}}_p\right)}^{-1}.$ The shrinkage parameter, k, should be estimated using the data, but it is unclear which value of k produces the best estimator. A large number of efforts have been aimed at accurately estimating the shrinkage parameter. Some proposed techniques for estimating k were proposed by refs. [11,12,13,14] and recently, [15]. We consider the estimator of the shrinkage parameter proposed by [6]. The estimator is $\widehat{k}={\widehat{\sigma }}^2/{\widehat{\alpha }}_{max}^2$, where ${\widehat{\alpha }}_{max}$ is the maximum element of $\widehat{\alpha }$.
• The modified ridge-type estimator with $g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)={\left({\mathbf{Z}}^{\top }\mathbf{Z}+k\left(1+d\right){\mathbf{I}}_p\right)}^{-1}$ is

  $$\begin{array}{ccc}\hfill {\widehat{\mathit{\alpha }}}_{MRT}\left(k,d\right)& =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+k\left(1+d\right){\mathbf{I}}_p\right)}^{-1}{\mathbf{Z}}^{\top }\mathbf{y}.\hfill \end{array}$$
  This estimator includes the ridge and OLS estimators as special cases. The authors of Ref. [7] suggest an iterative approach to estimating the shrinkage parameters. Their approach is as follows:

  (a)

  Let the initial estimate of d be calculated as $\widehat{d}=min\left({\widehat{\sigma }}^2/{\widehat{\alpha }}_i^2\right)$ where ${\widehat{\alpha }}_i$ is the ith element of $\widehat{\alpha }$.

  (b)

  Using $\widehat{d}$ from 1, estimate k as $\widehat{k}=p{\widehat{\sigma }}^2/\left({\sum }_{i=1}^p\left(1+\widehat{d}\right){\widehat{\alpha }}_i^2\right)$.

  (c)

  Let $\widehat{d}={\widehat{\sigma }}^2/\left(k{\widehat{{\alpha }_i}}^2\right)-1$. Estimate $d_{mrt}$ as ${\widehat{d}}_{mrt}=p/{\sum }_{i=1}^p{\widehat{d}}^{-1}$.

  (d)

  Use ${\widehat{d}}_{mrt}=\widehat{d}$ if ${\widehat{d}}_{mrt}$ is not between 0 and 1.

• The Liu estimator with the combined benefits of the ridge estimator and the Stein type estimator [16] was proposed by [8]. Given that rank(${\mathbf{X}}^{\top }\mathbf{X}$) = p, the Liu estimator is given by

  $$\begin{array}{ccc}\hfill {\widehat{\mathit{\alpha }}}_{LIU}\left(d\right)& =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{y}+d\widehat{\mathit{\alpha }}\right)\hfill \\ & =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{Z}+d{\mathbf{I}}_p\right){\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^{\top }\mathbf{y}\hfill \end{array}$$

  with $g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)={\left({\mathbf{Z}}^{\top }\mathbf{Z}+{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{Z}+d{\mathbf{I}}_p\right){\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}$. The estimator of the shrinkage parameter is

  $$\begin{array}{ccc}\hfill \widehat{d}& =& 1-{\widehat{\sigma }}^2\left[\frac{{\sum }_{i=1}^p\left(1/\left({\lambda }_i\left({\lambda }_i+1\right)\right)\right)}{{\sum }_{i=1}^p\left({\widehat{{\alpha }_i}}^2/{\left({\lambda }_i+1\right)}^2\right)}\right]\hfill \end{array}$$

  where ${\widehat{\alpha }}_i$ and ${\lambda }_i$ are the ith element of $\widehat{\mathit{\alpha }}$ and the vector of eigenvalues of ${\mathbf{X}}^{\top }\mathbf{X}$, respectively.
• A limitation of the estimator of the shrinkage parameter of the Liu estimator is that, in some instances, it has a negative value that affects the estimator’s performance [17]. The modified one-parameter Liu estimator proposed by [9] yields a positive value of $\widehat{d}$ and provides a significant improvement in the performance of the estimator. Given that rank(${\mathbf{X}}^{\top }\mathbf{X}$) = p, the modified one-parameter Liu estimator is given by

  $$\begin{array}{ccc}\hfill {\widehat{\mathit{\alpha }}}_{MLIU}\left(d\right)& =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{Z}-d{\mathbf{I}}_p\right){\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^{\top }\mathbf{y}\hfill \end{array}$$

  with $g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)={\left({\mathbf{Z}}^{\top }\mathbf{Z}+{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{Z}-d{\mathbf{I}}_p\right){\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}$. The optimal value of the shrinkage parameter is

  $$\begin{array}{ccc}\hfill \widehat{d}& =& \min \left(\frac{{\lambda }_i\left({\widehat{\sigma }}^2+{\widehat{{\alpha }_i}}^2\right)}{{\widehat{\sigma }}^2+{\lambda }_i{\widehat{{\alpha }_i}}^2}\right)\hfill \end{array}$$

  where ${\widehat{\alpha }}_i$ and ${\lambda }_i$ are the ith element of $\widehat{\mathit{\alpha }}$ and the vector of eigenvalues of ${\mathbf{X}}^{\top }\mathbf{X}$, respectively.
• Lastly, the Kibria–Lukman estimator, proposed by [10], is a one-parameter estimator that combines the characteristics of both the ridge and Liu estimators. Given that rank(${\mathbf{X}}^{\top }\mathbf{X}$) = p, the estimator is given by

  $$\begin{array}{ccc}\hfill {\widehat{\mathit{\alpha }}}_{KL}\left(k\right)& =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+k{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{y}-k\widehat{\mathit{\alpha }}\right)\hfill \\ & =& {\left({\mathbf{Z}}^{\top }\mathbf{Z}+k{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{Z}-k{\mathbf{I}}_p\right){\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^{\top }\mathbf{y}\hfill \end{array}$$

  with $g\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)={\left({\mathbf{Z}}^{\top }\mathbf{Z}+k{\mathbf{I}}_p\right)}^{-1}\left({\mathbf{Z}}^{\top }\mathbf{Z}-k{\mathbf{I}}_p\right){\left({\mathbf{Z}}^{\top }\mathbf{Z}\right)}^{-1}.$ The estimator of the shrinkage parameter proposed by [10], is

  $$\begin{array}{ccc}\hfill \widehat{k}& =& \frac{{\widehat{\sigma }}^2}{2{\widehat{{\alpha }_i}}^2+\left({\widehat{\sigma }}^2/{\lambda }_i\right)}\hfill \end{array}$$

  with ${\widehat{\alpha }}_i$ the ith element of $\widehat{\mathit{\alpha }}$ and ${\lambda }_i$ the ith eigenvalue of ${\mathbf{X}}^{\top }\mathbf{X}$. Following [18], the harmonic mean version is given by

  $$\begin{array}{ccc}\hfill \widehat{k}& =& \frac{{\widehat{\sigma }}^2}{{\sum }_{i=1}^p\left[2{\widehat{{\alpha }_i}}^2+\left({\widehat{\sigma }}^2/{\lambda }_i\right)\right]}.\hfill \end{array}$$

3. Efficient Model Estimation and Validation in Big Data Analytics

The consequences of multicollinearity in linear regression models are often mitigated using shrinkage methods. Despite the ease with which these techniques can be applied to small or moderate datasets, they encounter significant challenges in the big data domain. Some of these challenges are that the high volume of data frequently exceeds the storage capacity of a single computer and that the time required to obtain results becomes infeasible due to the computational burden of such a high volume of data. To address these big data challenges, Ref. [5] proposed an algorithm that utilises an array of sufficient statistics to obtain the ridge estimator of the multiple linear regression model. The sufficient statistics array can be computed and updated at the row or mini-batch level. However, the ridge estimator is nonlinear with respect to the shrinkage parameter, which complicates its optimisation. The shrinkage estimators given in Section 2 are frequently utilised to address the adverse effects of multicollinearity in linear regression models and has the added advantage of having closed-form solutions. In this section, we extend the work of [5] to include various other well-known shrinkage estimators and also adapt their algorithm to allow an efficient model validation by employing K-fold cross-validation. We utilise the array of sufficient statistics given in (3) in our proposed approach. The sufficient statistic array is used to determine the estimator of the regression coefficients as well as the estimator of the covariance matrix of the regression coefficients in the multiple linear regression model when using the estimators in Section 2. For estimation purposes, access to the entire dataset is not needed once the sufficient statistics array has been computed. The sufficient statistics array for shrinkage estimators with the form given in (2) is

$$\begin{array}{c}\hfill \mathbf{A}:={\displaystyle \sum _{i=1}^n}{\mathbf{A}}_i=\left[\begin{array}{cc}{S}_{yy}& {\mathbf{S}}_{zy}^{\top }\\ {\mathbf{S}}_{zy}& {\mathbf{S}}_{zz}\end{array}\right]=\left[\begin{array}{cc}\sum _{i=1}^n{y}_i^2& \sum _{i=1}^n{y}_i{\mathbf{z}}_i^{\top }\\ \sum _{i=1}^n{\mathbf{z}}_i{y}_i& \sum _{i=1}^n{\mathbf{z}}_i{\mathbf{z}}_i^{\top }\end{array}\right].\end{array}$$

(3)

Considering the sufficient statistics array, the form of the estimators in Section 2 is

$$\begin{array}{c}\hfill {\widehat{\mathit{\alpha }}}_g=g\left({\mathbf{S}}_{zz}\right){\mathbf{S}}_{zy}.\end{array}$$

The storage capacity of a single computer may be much less than the size of a large dataset, and hence the computation of (3) has to be distributed across multiple processors, in which case, algorithms based on multiple processors should be utilised. Parallel computation is extremely useful in this regard, and hence the concept of MapReduce is important in our approach since it interleaves parallel and sequential computation. The MapReduce concept enables us to perform K-fold cross-validation with ease, and hence we base our algorithm on a strategy suitable for multiple processors. Model validation often relies on resampling techniques that utilise two datasets, one for training the model and one for evaluating the prediction performance of the model on unseen data. The testing dataset must be accessed after the regression coefficients are estimated, irrespective of the method used to obtain the regression coefficients, since the predicted response is calculated using only the observed testing data and the estimated regression coefficients. The advantage of using K-fold cross-validation as a model validation technique lies in the fact that (3) can be obtained for each of the K training datasets without any more access to the data than would be required to calculate (3) for the complete dataset. Suppose we partition the dataset into K blocks containing $n_k$ observations. The sufficient statistics array of each block is then given by ${\sum }_{i=1}^{n_k}{\mathbf{A}}_i$, such that the sufficient statistics array of the complete dataset is

$$\begin{array}{c}\hfill \mathbf{A}:={\displaystyle \sum _{k=1}^K}{\mathbf{A}}_{\left(k\right)}:={\displaystyle \sum _{k=1}^K}{\displaystyle \sum _{i=1}^{n_k}}{\mathbf{A}}_i.\end{array}$$

We can easily obtain the array of sufficient statistics for any training dataset since each of the K training datasets contains $(K-1)$ of the K blocks. Algorithm 1 is used for the model estimation and validation of the shrinkage estimators given in Section 2. Furthermore, Algorithm 1 is partitioned into four stages. The eigenvectors and eigenvalues of ${\mathbf{X}}^{\top }\mathbf{X}$ corresponding to the K training datasets are calculated in the first stage. The second step involves calculating the sufficient statistics array for each of the K blocks. In the third stage, we obtain the estimators proposed in Section 2, and finally, the prediction performance is evaluated in the fourth stage.

It is important to note that the estimators we considered have a general mathematical form, and hence the approach is not limited to only those estimators presented in this paper. Should feature selection be of interest, our approach can be applied to various other penalties, such as the elastic net, by using the local quadratic approximation of [19]. In these cases, it is not necessary to use the canonical form of the multiple linear regression model since the optimal turning parameter needs to be obtained through K-fold cross-validation. Furthermore, with minor modifications to the algorithm, data issues such as outliers could be addressed with the robust ridge estimator [20]. However, our main concern is multicollinearity, and hence the simulation and application considered the five estimators given in Section 2.

Algorithm 1 Model estimation and validation based on the sufficient statistics array

Input: Subsets of the data Output: ${\widehat{\mathit{\alpha }}}_g$, ${\widehat{\mathit{\sigma }}}^2$, $V\left({\widehat{\alpha }}_g\right)$ and measures of prediction performance based on the testing dataset. 1. Divide the data into K blocks. MAP tasks 2. For each of the K blocks: 2.1 Calculate ${\mathbf{S}}_{xx\left(k\right)}={\sum }_{i=1}^{n_k}{\mathbf{x}}_i{\mathbf{x}}_i^{\top }$ Reduce tasks 3. Let j be the index of the block to be excluded from each of the K training datasets such that $-j$ represent the remaining $K-1$ blocks. For the K training datasets: 3.1 Calculate ${\sum }_{k\ne j}{\mathbf{S}}_{xx\left(k\right)}$ 3.2 Calculate ${\mathbf{Q}}_{\left(-j\right)}$, the eigenvectors of ${\sum }_{k\ne j}{\mathbf{S}}_{xx\left(k\right)}$. 3.3 Calculate ${\left({\lambda }_1,\dots ,{\lambda }_p\right)}_{\left(-j\right)}$, the eigenvalues of ${\sum }_{k\ne j}{\mathbf{S}}_{xx\left(k\right)}$. MAP tasks 4. For each of the K blocks: 4.1 Calculate $S_{yy\left(k\right)}$, ${\mathbf{S}}_{zy\left(k\right)}$, ${\mathbf{S}}_{zz\left(k\right)}$ and ${\mathbf{A}}_{\left(k\right)}$ for each ${\mathbf{Q}}_{\left(-j\right)}$ with $j\in \left\{1,\dots ,K\right\}$. Reduce tasks 5. For each of the K training datasets with $j\in \left\{1,\dots ,K\right\}$: 5.1 Calculate ${\sum }_{k\ne j}{\mathbf{A}}_{\left(k\right)}$. 5.2 Calculate ${\widehat{\mathit{\alpha }}}_{\left(-j\right)}$, and ${\widehat{\sigma }}_{\left(-j\right)}^2$. 5.3 Calculate the shrinkage parameter of the estimator under consideration. 5.4 Calculate ${\widehat{\mathit{\alpha }}}_g$ and $V\left({\widehat{\mathit{\alpha }}}_g\right)$. 6. For each of the K testing datasets: 6.1 Calculate the predicted response. 6.2 Calculate prediction performance measures based on the observed and predicted responses.

4. Numerical Analysis

In this section, the performance of the proposed approach is evaluated through an extensive simulation study and an application on a real-world dataset. A high-performance computer with 60 cores and 256GB RAM was used for all computations.

4.1. Simulation Study

The simulation study considered two experiments. The first experiment aimed to compare the execution time of our proposed approach with the execution time of the traditional approach, where we performed the relevant matrix algebra from Section 2 using the entire dataset at once. The second experiment was to verify that our algorithm yielded the same estimators and performance measures as the traditional approach. Theoretically, our proposed approach should yield the same answers as the traditional approach, which was verified by the results of the second experiment. The predictor variables were generated using

$$\begin{array}{ccc}\hfill X_{ij}& =& {\left(1-{\gamma }^2\right)}^{1/2}{z}_{ij}+\gamma z_{i\left(p+1\right)}\hfill \end{array}$$

with $i\in \left\{1,\dots ,n\right\}$, $j\in \left\{1,\dots ,p\right\}$, $z_{ij}$ independent standard normal random numbers and ${\gamma }^2$ the correlation between any two predictor variables. The response was generated by

$$\begin{array}{ccc}\hfill y_i& =& {\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\dots +{\beta }_p{x}_{ip}+e_i\hfill \end{array}$$

with $i\in \left\{1,\dots ,n\right\}$ and $e_i\sim N\left(0,{\sigma }^2\right)$. We considered $n=10,000,000$, $p=100$, $\gamma =0.9,0.99$, and ${\sigma }^2=1$. The regression coefficients were generated from a uniform($-1,1$) distribution. Furthermore, we included both 5-fold and 10-fold cross-validation and considered three performance measures, namely the mean squared error (MSE) of the estimators, the mean prediction error (MPE) of the predicted response on the testing dataset, and the symmetric mean absolute percentage error (SMAPE) of the predicted response on the testing dataset. For an extended performance evaluation, the measures used by [21] can also be considered. The MSE, MPE, and SMAPE were calculated using

$$\begin{array}{ccc}\hfill \mathrm{MSE}\left(\widehat{\mathit{\beta }}\right)& =& \frac{1}{K}{\displaystyle \sum _{k=1}^K}{\left({\widehat{\mathit{\beta }}}_k-\mathit{\beta }\right)}^{\top }\left({\widehat{\mathit{\beta }}}_k-\mathit{\beta }\right)\hfill \\ \hfill \mathrm{MPE}\left(\widehat{\mathit{y}}\right)& =& \frac{1}{n}{\displaystyle \sum _{k=1}^K}{\left({\widehat{\mathit{y}}}_k-\mathit{y}\right)}^{\top }\left({\widehat{\mathit{y}}}_k-\mathit{y}\right).\hfill \\ \hfill \mathrm{SMAPE}\left(\widehat{\mathit{y}}\right)& =& \frac{100\%}{n}{\displaystyle \sum _{k=1}^K\sum _{i=1}^{\frac{n}{K}}}\frac{|{\widehat{y}}_{ik}-y_{ik}|}{\left(\left|y_{ik}\right|+\left|{\widehat{y}}_{ik}\right|\right)/2}.\hfill \end{array}$$

Figure 1 and

Table 1. Difference in performance measures calculated for the approach proposed in Section 3 and the traditional approach given in Section 2 when $K=5$ and $\mathbf{K}=\mathbf{10}$.

$\mathit{\gamma }$	Estimator	$\mathbf{MSE}\left(\widehat{\mathit{\beta }}\right)$	$\mathbf{MSE}\left(\widehat{\mathit{y}}\right)$	$\mathbf{SMAPE}\left(\widehat{\mathit{y}}\right)$
0.9	Ridge	$3.866\times {10}^{-15}$ $\mathbf{4}.\mathbf{387}\times {\mathbf{10}}^{-\mathbf{11}}$	−2.776 $\times {10}^{-17}$ $\mathbf{1}.\mathbf{206}\times {\mathbf{10}}^{-\mathbf{13}}$	6.501 $\times {10}^{-13}$ $\mathbf{3}.\mathbf{270}\times {\mathbf{10}}^{-\mathbf{9}}$
0.9	Modified ridge	−1.784 $\times {10}^{-9}$ $-\mathbf{9}.\mathbf{465}\times {\mathbf{10}}^{-\mathbf{9}}$	3.679 $\times {10}^{-10}$ $\mathbf{3}.\mathbf{539}\times {\mathbf{10}}^{-\mathbf{10}}$	−4.567 $\times {10}^{-7}$ $-\mathbf{6}.\mathbf{069}\times {\mathbf{10}}^{-\mathbf{7}}$
0.9	Liu	2.568 $\times {10}^{-8}$ $-\mathbf{5}.\mathbf{866}\times {\mathbf{10}}^{-\mathbf{8}}$	1.825 $\times {10}^{-9}$ $\mathbf{1}.\mathbf{263}\times {\mathbf{10}}^{-\mathbf{9}}$	1.760 $\times {10}^{-6}$ $-\mathbf{4}.\mathbf{399}\times {\mathbf{10}}^{-\mathbf{6}}$
0.9	Modified Liu	3.077 $\times {10}^{-9}$ $\mathbf{5}.\mathbf{432}\times {\mathbf{10}}^{-\mathbf{9}}$	5.115 $\times {10}^{-12}$ $-\mathbf{1}.\mathbf{852}\times {\mathbf{10}}^{-\mathbf{12}}$	3.470 $\times {10}^{-7}$ $\mathbf{2}.\mathbf{999}\times {\mathbf{10}}^{-\mathbf{7}}$
0.9	Kibria–Lukman	−1.865 $\times {10}^{-15}$ $-\mathbf{3}.\mathbf{193}\times {\mathbf{10}}^{-\mathbf{16}}$	−1.388 $\times {10}^{-17}$ $-\mathbf{6}.\mathbf{939}\times {\mathbf{10}}^{-\mathbf{17}}$	1.421 $\times {10}^{-13}$ $\mathbf{4}.\mathbf{974}\times {\mathbf{10}}^{-\mathbf{14}}$
0.999	Ridge	4.184 $\times {10}^{-12}$ $-\mathbf{4}.\mathbf{048}\times {\mathbf{10}}^{-\mathbf{8}}$	2.637 $\times {10}^{-16}$ $-\mathbf{2}.\mathbf{637}\times {\mathbf{10}}^{-\mathbf{11}}$	2.540 $\times {10}^{-12}$ $-\mathbf{1}.\mathbf{835}\times {\mathbf{10}}^{-\mathbf{8}}$
0.999	Modified ridge	−6.462 $\times {10}^{-6}$ $-\mathbf{4}.\mathbf{155}\times {\mathbf{10}}^{-\mathbf{6}}$	−6.894 $\times {10}^{-9}$ $\mathbf{5}.\mathbf{182}\times {\mathbf{10}}^{-\mathbf{9}}$	−7.709 $\times {10}^{-6}$ $-\mathbf{1}.\mathbf{219}\times {\mathbf{10}}^{-\mathbf{7}}$
0.999	Liu	8.086 $\times {10}^{-5}$ $-\mathbf{4}.\mathbf{591}\times {\mathbf{10}}^{-\mathbf{6}}$	9.276 $\times {10}^{-8}$ $\mathbf{8}.\mathbf{133}\times {\mathbf{10}}^{-\mathbf{9}}$	−4.793 $\times {10}^{-6}$ $-\mathbf{9}.\mathbf{801}\times {\mathbf{10}}^{-\mathbf{6}}$
0.999	Modified Liu	3.448 $\times {10}^{-6}$ $\mathbf{5}.\mathbf{979}\times {\mathbf{10}}^{-\mathbf{6}}$	2.908 $\times {10}^{-10}$ $\mathbf{1}.\mathbf{749}\times {\mathbf{10}}^{-\mathbf{10}}$	1.503 $\times {10}^{-6}$ $\mathbf{1}.\mathbf{051}\times {\mathbf{10}}^{-\mathbf{6}}$
0.999	Kibria–Lukman	2.920 $\times {10}^{-13}$ $\mathbf{4}.\mathbf{710}\times {\mathbf{10}}^{-\mathbf{13}}$	2.776 $\times {10}^{-17}$ $-\mathbf{1}.\mathbf{388}\times {\mathbf{10}}^{-\mathbf{17}}$	1.386 $\times {10}^{-13}$ $\mathbf{1}.\mathbf{847}\times {\mathbf{10}}^{-\mathbf{13}}$

contain the results obtained in the simulation study.

The computational time required to obtain the estimated model using the proposed approach was ±eight times faster than the traditional approach for the 5-fold cross-validation and ±five times faster for the 10-fold cross-validation. The time required to obtain the predicted responses and the performance measures on the testing dataset was less for the traditional approach. When using the traditional approach, the complete dataset was imported before the model estimation was performed. The dataset was never discarded, and hence it was readily available when prediction on the testing set commenced. This came at a cost since the memory used to store the complete dataset was unavailable at the model estimation stage. Although our proposed approach required the testing datasets to be imported again, the total time required to obtain the models, predict, and calculate the relevant performance measures was significantly less when using the approach proposed in Section 3.

Furthermore, the difference between the performance measures calculated for the proposed approach and the traditional approach, where we performed the relevant matrix algebra from Section 2 using the entire dataset at once, was extremely close to zero. That is, the proposed approach was as accurate as the traditional approach where the mean squared error of the estimators, the mean prediction error, and the symmetric mean absolute percentage error of the predicted response on the testing dataset were concerned. The two experiments were also conducted on datasets where the level of multicollinearity was low, e.g., cases where $\gamma =0.5,0.6$. The results obtained were very similar to the results shown in Figure 1 and Table 1, and the conclusions arrived at also held for lower levels of multicollinearity.

4.2. Application

The dataset used in this application was obtained from the Bureau of Transportation Statistics, which provides information on U.S. transportation systems. The dataset contains information such as the origin and destination airport of various flights, the time at which flights are scheduled to arrive and depart, the delay in departing and arriving flights, the time and distance to destination airports, and the date and time of these flights, for major air carriers. We selected air traffic data from 2010 to 2020, resulting in a dataset with $\pm 6,000,000$ observations and $\pm 90$ features. The response variable was the arrival delay in minutes. As in the simulation study, we considered two experiments. The first experiment aimed to compare the execution time of our proposed approach with the execution time of the traditional approach, where we performed the relevant matrix algebra from Section 2 using the entire dataset at once. The second experiment was to show that our algorithm was as accurate as the traditional approach. The performance measures used were the mean prediction error (MPE) and the symmetric mean absolute percentage error (SMAPE) of the predicted response on the testing dataset. The results obtained in the application are given in Figure 2 and

Table 2. Difference in performance measures calculated for the proposed approach and the traditional approach given in Section 2 when $K=5$ and $\mathbf{K}=\mathbf{10}$.

Estimator	$\mathbf{MPE}\left(\widehat{\mathit{y}}\right)$	$\mathbf{SMAPE}\left(\widehat{\mathit{y}}\right)$
Ridge	$4.025\times {10}^{-5}$ $-\mathbf{2}.\mathbf{948}\times {\mathbf{10}}^{-\mathbf{4}}$	6.719$\times {10}^{-7}$ $\mathbf{1}.\mathbf{781}\times {\mathbf{10}}^{-\mathbf{5}}$
Modified ridge	1.803$\times {10}^{-3}$ $\mathbf{1}.\mathbf{766}\times {\mathbf{10}}^{-\mathbf{3}}$	2.871$\times {10}^{-4}$ $\mathbf{1}.\mathbf{673}\times {\mathbf{10}}^{-\mathbf{3}}$
Liu	5.253$\times {10}^{-4}$ $\mathbf{4}.\mathbf{835}\times {\mathbf{10}}^{-\mathbf{4}}$	1.496$\times {10}^{-5}$ $\mathbf{3}.\mathbf{623}\times {\mathbf{10}}^{-\mathbf{5}}$
Modified Liu	3.131$\times {10}^{-5}$ $\mathbf{9}.\mathbf{429}\times {\mathbf{10}}^{-\mathbf{5}}$	9.764$\times {10}^{-7}$ $\mathbf{1}.\mathbf{000}\times {\mathbf{10}}^{-\mathbf{5}}$
Kibria–Lukman	7.488$\times {10}^{-6}$ $\mathbf{7}.\mathbf{397}\times {\mathbf{10}}^{-\mathbf{5}}$	5.760$\times {10}^{-7}$ $\mathbf{9}.\mathbf{347}\times {\mathbf{10}}^{-\mathbf{6}}$

.

From the results given in Figure 2 and Table 2, it is clear that the total computational time required to obtain the models, predict, and calculate the relevant performance measures was significantly less when using the approach proposed in Section 3. Furthermore, the differences between the performance measures calculated for the proposed approach and the traditional approach were close to zero.

5. Conclusions and Future Work

Shrinkage methods are often used to mitigate the consequences of multicollinearity in linear regression models. These methods can easily be applied to small- or moderate-size datasets but encounter significant challenges, such as the memory and time required to obtain an optimal model in the big data domain. The goal of this paper was to propose a method for the efficient model estimation and validation of various well-known shrinkage estimators—mainly used to address multicollinearity—when considering large datasets. We proposed an algorithm that utilised an array of sufficient statistics for model estimation of various ridge- and Liu type estimators with closed-form solutions for the tuning parameter. Our algorithm also enabled model validation through K-fold cross-validation. A simulation study and an application on a real-world dataset illustrated the efficiency and accuracy of our approach. The estimators that we considered have a general mathematical form, and hence our approach is not limited to only the estimators presented in this paper. Furthermore, some ideas for future work include:

• Extending the methodology to generalised linear models.
• Extending the methodology to streaming data.
• Extending the methodology to robust estimators.
• Extending the methodology to estimators for feature selection, for which one may consider measures of validation and evaluation similar to those in [21].