Shrinkage Approaches for Ridge-Type Estimators Under Multicollinearity

1. Introduction

Regression analysis provides answers to queries concerning the dependence of a variable, known as the response, on one or a set of independent variables, known as predictors. These include problems like predicting the response value for a given collection of predictors or identifying the most significant group of predictors with a plausible impact on the response variable. A commonly used method to estimate the functional relationship between the response variable and the set of predictors is the ordinary least squares method (OLS). This method provides the best linear unbiased estimator (BLUE) of the regression vector of coefficients based on the Gauss-Markov theorem. This result holds, assuming the columns of the design matrix, which are represented by the independent variables, are not correlated. However, the OLS estimator becomes less efficient if a strong or near-to-strong linear relationship exists among the columns of design matrix, known as multicollinearity.

The multicollinearity issue lessens the accuracy of the OLS estimated coefficients, also the coefficient estimates can exhibit significant fluctuations depending on the inclusion of different independent variables in the model, beside to a high degree of sensitivity to even minor change in the regression model. There are several estimation methods developed to improve the OLS estimator when multicollinearity exists. For instance, Hoerl and Kennard [1] proposed the ridge regression estimator analytically and numerically the superiority of the new estimator over the OLS estimator. Later, Liu K [2] introduced a novel biased estimator and verified analytically and numerically the superiority of the new estimator over the OLS estimator known as Liu type estimator . Recently, Kibria and Lukman [3] developed a new ridge type estimator for the regression parameter vector. They showed the new estimator performed better than the ridge and Liu estimators in terms of the MSE criterion.

The OLS technique uses the sample information obtained from the data set to draw inference about the unknown population regression vector of parameters, which refers to classical inference method. However, in Bayesian context, the sample information is combined with non-sample information (NSI) to make an inference about the regression parameters. Such NSI may not be available at all times. However, model selection and building procedures as AIC, BIC, penalization methods, or machine learning algorithms could still be utilized to yield NSIs. Bancroft [4] was one of the first to try estimating regression coefficients by merging NSIs with sample information and produced what is known as the pretest estimator. The Pretest estimator relies on evaluating the statistical significance of certain regression coefficients. After the determination is made, the pretest selects either the estimate of the full model or the estimator of the revised model, known as sub-model, which has fewer number of coefficients. The pertest algorithm chooses either the full or sub model estimators using binary weights. A modified version of the pretest estimator, known as the shrinkage estimator, has been formulated by Stein [5]. This estimator utilizes smooth weights to merge the estimations from both the overall and sub-model estimators. The regression coefficients are modified to converge towards a desired value that is influenced by the NSI. Nevertheless, the enhanced shrinkage estimator sometimes encounters an occurrence of excessive reduction. Following that, a more improved version of this estimator is proposed by Stein [6] to efficiently address the issue of excessive shrinkage, known as the positive shrinkage estimator.

Many researches have been drawn to the idea of using shrinkage and pretest estimation approaches. For instance, [7] proposed a novel approach utilizing pretest and shrinkage approaches to accurately estimate the regression coefficients vector of the marginal model for multinomial responses. [8] presented the Liu-type pretest, shrinkage, and positive shrinkages estimators for the conditional autoregressive regression model’s large-scale effect parameter vector, and demonstrated that these estimators are more efficient than Liu-type estimators. By employing the concept of shrinkage, [9] proposed an enhanced version of the Liu-type estimator. The proposed method’s superiority was demonstrated using analytical and numerical data. Subsequently, [10] proposed the utilization of the ridge estimator as an appropriate method for managing high-dimensional multicollinear data. [11] proposed the use of pretest and shrinkage ridge estimation techniques for the linear regression model. Demonstrated the advantages of employing the suggested estimators alongside specific penalty estimators. [12] introduced the pretest and shrinkage approaches that use generalized ridge regression estimation to address issues related to multicollinearity and high-dimensional situations.  [13] introduced an effective shrinkage and penalty estimators for regression coefficients in spatial error models. It showcases the efficiency improvements using asymptotic and numerical studies. In a recent study, [14] shown that the ridge-type pretest and shrinkage estimators outperformed the maximum likelihood estimator in the presence of multicollinearity in a spatial error regression model. To obtain further information on the shrinkage estimators, please consult with [15,16,17], among other relevant sources.

The primary objectives of this research are threefold. The initial goal of this study is to improve the ridge estimator (RE) proposed by [3] of the regression vector by integrating pretest and shrinkage approaches in the presence of multicollinearity among the regressor variables. This enhancement primarily focuses on situations where particular regression coefficients are considered insignificant. The second purpose is to examine the analytical findings on the bias and quadratic risks of the estimators for the newly proposed set. The third step involves evaluating the performance of the estimators by numerical simulation, specifically by measuring the mean squared error. Additionally, the estimators will be assessed using real-world data to determine the prediction error.

The subsequent sections of the paper are structured in accordance with our objectives. In Section (2), we provide the regression model and the KL estimator. The strategy improvements for the RE estimator are outlined in Section (3). In Section (4), we studied some analytical properties of the array of estimators, specifically the bias and quadratic risks. Some penalizing techniques are presented in Section (5), and we presented three machine learning algorithms in Section (6). In Section (7), we conducted a comparison of the array of estimators using both Monte Carlo simulation and a real data example. Section (8) provides concluding remarks, and a supplementary section is appended at the conclusion of the manuscript, which includes the proofs of the analytical findings.

2. Statistical Model and New Ridge Estimation

To explain the problem, let us consider the linear regression model given below

$$\begin{array}{ccc}\hfill \mathbf{Y}& =& \mathbf{X}\mathsf{\theta }+\mathbf{\upsilon },\hfill \end{array}$$

(1)

where $\mathbf{Y}$ is an $(n\times 1)$ vector of random responses, $\mathbf{X}={\left({\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_p\right)}^T$ is an $(n\times p)$ full rank design matrix with ${\mathbf{X}}_j={(X_{j1},X_{j2},\dots ,X_{jp})}^T$, $\mathsf{\theta }=({\theta }_1,{\theta }_2,\dots ,{\theta }_p)$ is a vector of unknown but fixed regression parameters, and $\mathbf{\upsilon }={({\upsilon }_1,{\upsilon }_2,\dots ,{\upsilon }_n)}^T$ is an $(n\times 1)$ random error vector such that the expected value and the variance of ${\upsilon }_i$ are, respectively $E\left({\upsilon }_i\right)=0$, and $Var\left({\upsilon }_i\right)={\tau }^2>0$. The ordinary least squares (OLS) estimator of $\mathsf{\theta }$ is given by:

$$\begin{array}{ccc}\hfill \widehat{\mathsf{\theta }}& =& {\mathbf{G}}^{-1}{\mathbf{X}}^T\mathbf{Y},\hfill \end{array}$$

where $\mathbf{G}={\mathbf{X}}^T\mathbf{X}$. Based on Gauss-Markov theorem, $\widehat{\mathsf{\theta }}$ is the best linear unbiased estimator (BLUE) of $\mathsf{\theta }$, which is the case when the columns of the design matrix $\mathbf{X}$ are not correlated, with $Var\left(\widehat{\mathsf{\theta }}\right)={\tau }^2{\mathbf{G}}^{-1}$. However, $\widehat{\mathsf{\theta }}$ becomes less efficient if a strong or near to strong linear relationship exists among the columns of $\mathbf{X}$, which is known as multicollinearity. The ridge-type estimator proposed by [3] estimator of the regression coefficient $\mathsf{\theta }$ is denoted by ${\widehat{\mathsf{\theta }}}^{RE}$ and given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{RE}& =& {(\mathbf{G}+k{\mathbf{I}}_p)}^{-1}(\mathbf{G}-k{\mathbf{I}}_p)\widehat{\mathsf{\theta }},\hfill \end{array}$$

(2)

where ${\mathbf{I}}_p$ is a $(p\times p)$ identity matrix and $k\ge 0$ known as a biasing parameter to be estimated using the data.

The incorporation of NSI in a model often involves the inclusion of a hypothesized restriction on the model parameters, resulting in the emergence of potential sub-models. Bayesian statistical approaches have been developed in response to the need of incorporating NSI into models fitted to objective sample data. This allows for the consideration of the uncertainty introduced by both sources of information. If the NSI asserts that some of the regression coefficients are irrelevant, then it is possible to integrate this information into the estimation process via testing a linear hypothesis of the form:

$$\begin{array}{ccc}\hfill H_0:\mathbf{R}\mathsf{\theta }& =& {\mathbf{0}}_q,\hfill \end{array}$$

(3)

where $\mathbf{R}$ is a $(q\times p)$ known matrix of rank $(q\le p)$, and ${\mathbf{0}}_q$ is a $(q\times 1)$ vector of zeros.

The new sub-model often incorporates a reduced number of regression variables, facilitating interpretation and mitigating the complexity associated with a large number of irrelevant variables. A Candidate sub-model may be obtained by using some known variable selection methods, such as AIC , BIC, or by penalization techniques as ridge estimation, LASSO, Elastic net, SCAD, or adaptive LASSO, among others. Under the restriction given by the hypothesis in (3), and using Lagrange multipliers the sub-model estimator of $\mathsf{\theta }$, denoted by ${\widehat{\mathsf{\theta }}}^R$ is given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^R& =& \widehat{\mathsf{\theta }}-{\mathbf{G}}^{-1}{\mathbf{R}}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}\mathbf{R}\widehat{\mathsf{\theta }}\hfill \end{array}$$

(4)

Obviously, if the restriction in (3) is correct, then ${\widehat{\mathsf{\theta }}}^R=\widehat{\mathsf{\theta }}$, which is unbiased estimator for $\mathsf{\theta }$. On the other hand, when the restriction is not true, ${\widehat{\mathsf{\theta }}}^R$ becomes less efficient than $\widehat{\mathsf{\theta }}$, particularly in cases where there is ongoing multicollinearity among the columns of the $\mathbf{X}$ columns. This can be diminished by employing the Kibria and Lukman estimation method to calculate the coefficients of sub-model regression, which will be denoted by ${\widehat{\mathsf{\theta }}}^{RRE}$ and given as follows:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{RRE}& =& {(\mathbf{G}+k{\mathbf{I}}_p)}^{-1}(\mathbf{G}-k{\mathbf{I}}_p){\widehat{\mathsf{\theta }}}^R\hfill \end{array}$$

(5)

3. Strategy Improvements for the RE Estimator

Given that a sub-model has been acquired via the use of penalization methods or machine learning algorithms, we will explore various shrinkage and pretest procedures to combine the estimations from both the full and sub-models. First, we establish the definition of a test statistic, designated as ${\varphi }_n$, which is used to test the hypothesis stated in equation (3) as follows:

$$\begin{array}{ccc}\hfill {\varphi }_n& =& {\displaystyle \frac{{\left(\mathbf{R}\widehat{\mathsf{\theta }}\right)}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}\left(\mathbf{R}\widehat{\mathsf{\theta }}\right)}{q{s}^2}},\hfill \end{array}$$

where $s^2={\displaystyle \frac{{(\mathbf{Y}-\mathbf{X}\widehat{\mathsf{\theta }})}^T(\mathbf{Y}-\mathbf{X}\widehat{\mathsf{\theta }})}{n-p}}$ is an estimator of $Var\left({\nu }_i\right)={\tau }^2$. The test statistic ${\varphi }_n$, assuming the null hypothesis (3), is asymptotically distributed according to a $F-$distribution with $(q,r),r=n-p$ degrees of freedom. However, if the alternative hypothesis $\mathbf{R}\mathsf{\theta }\ne {\mathbf{0}}_q$ holds, the test statistic follows a non-central $F-$distribution with a non-centrality parameter ${\Delta }^2={\displaystyle \frac{{\left(\mathbf{R}\mathsf{\theta }\right)}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}\left(\mathbf{R}\mathsf{\theta }\right)}{{\tau }^2}}$. An improved set of shrinkage estimators using [3] technique can be formulated as follows:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{RE-Shrinkage}& =& {\widehat{\mathsf{\theta }}}^{RE}-\left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right)f\left({\varphi }_n\right),\hfill \end{array}$$

(6)

where $f\left({\varphi }_n\right)$ is a Borel measurable function of ${\varphi }_n$ and certain selections of $f\left({\varphi }_n\right)$ may provide shrinkage estimators that are both plausible and helpful.

3.1. Preliminary and Shrinkage Estimators

The preliminary test estimator can be obtained by when the function $f\left({\varphi }_n\right)=I\left({\varphi }_n\le F_{q,r}\left(\alpha \right)\right)$, it is denoted by ${\widehat{\mathsf{\theta }}}^{REPT}$, and given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{REPT}& =& {\widehat{\mathsf{\theta }}}^{RE}-\left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right)I\left({\varphi }_n\le F_{q,r}\left(\alpha \right)\right),\hfill \end{array}$$

(7)

where, $I(.)$ is the indicator function, $\alpha$ is the level of significance, and $F_{q,r}\left(\alpha \right)$ represents the upper critical $\alpha -$value of the $F-$distribution with $(q,r)$ degrees of freedom. Clearly, if the indicator function is zero, then ${\widehat{\mathsf{\theta }}}^{REPT}={\widehat{\mathsf{\theta }}}^{RE}$, otherwise ${\widehat{\mathsf{\theta }}}^{REPT}={\widehat{\mathsf{\theta }}}^{RER}$. One limitation of this estimator is its reliance on a binary choice of the two estimators that is influenced by the level of significance $\alpha$. A more refined method of assigning weights may be attained by setting $f\left({\varphi }_n\right)=\pi {\varphi }_n^{-1},\pi ={\displaystyle \frac{(q-2)(n-p)}{q(n-p+2)}}$ the James-Stein shrinkage estimator, which is formally stated as follows.

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{RES}& =& {\widehat{\mathsf{\theta }}}^{RE}-\left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right)\pi {\varphi }_n^{-1}.\hfill \end{array}$$

(8)

The shrinkage estimator ${\widehat{\mathsf{\theta }}}^{RES}$ experiences over-shrinkage, resulting in negative $\mathsf{\theta }$ coordinates when $({\varphi }_n<\pi )$. The positive James–Stein estimator solves this issue. It is denoted by ${\widehat{\mathsf{\theta }}}^{REPS}$, and can by obtained by setting $f\left({\varphi }_n\right)=\pi {\varphi }_n^{-1}+(1-\pi {\varphi }_n^{-1})I({\varphi }_n\le \pi )$

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{REPS}& =& {\widehat{\mathsf{\theta }}}^{RE}\hfill \\ \hfill & -& \left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right)\left[\pi {\varphi }_n^{-1}+(1-\pi {\varphi }_n^{-1})I({\varphi }_n\le \pi )\right],\hfill \end{array}$$

(9)

which can be simplifies as:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{REPS}& =& {\widehat{\mathsf{\theta }}}^{RES}-\left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right)\hfill \\ \hfill & & \left(1-\pi {\varphi }_n^{-1}\left)I\right({\varphi }_n\le \pi \right).\hfill \end{array}$$

(10)

3.2. Modified Preliminary and Shrinkage Estimators

Another way to enhance the regression estimate when the null hypothesis in (3) is true, which could any type of NSI, is to consider the linear shrinkage estimator of the sub and full models RE estimator, [25] provided evidence that the estimator is competitive when the NSI is correct. The estimator, denoted by ${\widehat{\mathsf{\theta }}}^{LSRE}$, and given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{LSRE}& =& d{\widehat{\mathsf{\theta }}}^{RRE}+(1-d){\widehat{\mathsf{\theta }}}^{RE}\hfill \\ & =& {\widehat{\mathsf{\theta }}}^{RE}-d\left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right),\hfill \end{array}$$

(11)

where $0\le d\le 1$ serves as a tuning parameter which is selected to minimize the estimator’s mean-squared error. Obviously, when $d=0$, ${\widehat{\mathsf{\theta }}}^{LSRE}={\widehat{\mathsf{\theta }}}^{RE}$, and when $d=1$, ${\widehat{\mathsf{\theta }}}^{LSRE}={\widehat{\mathsf{\theta }}}^{RRE}$. Additionally, we can employ the idea proposed by [26] to produces a new RE estimator that is often referred to as the shrinkage pretest estimator. It is denoted by ${\widehat{\mathsf{\theta }}}^{SPTRE}$ and given as follows:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{SPTRE}& =& {\widehat{\mathsf{\theta }}}^{RE}-d\left({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}\right)\hfill \\ & & I\left({\varphi }_n\le F_{q,r}\left(\alpha \right)\right)\hfill \end{array}$$

(12)

4. Analytical Properties

This section presents the bias and quadratic risk functions of the proposed RE estimators. For this purpose, assume ${\widehat{\mathsf{\theta }}}^{*RE}$ is any of these RE-type estimators; consequently the bias of ${\widehat{\mathsf{\theta }}}^{*}$ will be $B\left({\widehat{\mathsf{\theta }}}^{*}\right)=E\left({\widehat{\mathsf{\theta }}}^{*}-\mathsf{\theta }\right)$. The bias expressions are given in the following theorem.

Theorem 1

 

Let ${\mathbf{A}}_k=(\mathbf{G}+k{\mathbf{I}}_p)(\mathbf{G}-k{\mathbf{I}}_p)$, and $\mathbf{B}={\mathbf{G}}^{-1}{\mathbf{R}}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}{\mathbf{R}}^T$, then:

• $B\left({\widehat{\mathsf{\theta }}}^{RE}\right)=[{\mathbf{A}}_k-{\mathbf{I}}_p]\mathsf{\theta }$
• $B\left({\widehat{\mathsf{\theta }}}^{RRE}\right)=B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }$
• $B\left({\widehat{\mathsf{\theta }}}^{REPT}\right)=B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }{H}_{q+2,r}(c_{\alpha },{\Delta }^2)$
• $B\left({\widehat{\mathsf{\theta }}}^{RES}\right)=B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-\pi {\mathbf{A}}_k\mathbf{B}\mathsf{\theta }E\left(F_{q+2,r}^{-1}\left({\Delta }^2\right)\right)$
• $B\left({\widehat{\mathsf{\theta }}}^{REPS}\right)=B\left({\widehat{\mathsf{\theta }}}^{RES}\right)-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }\left\{\pi E\left(F_{q+r,2}^{-1}\left({\Delta }^2\right)\right)+E\left((1-\pi F_{q+2,r}^{-1}(c_{\alpha };{\Delta }^2))I(F_{q+2,r}^{-1}\left({\Delta }^2\right)\le \pi )\right)\right\}$
• $B\left({\widehat{\mathsf{\theta }}}^{LSRE}\right)=B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-d{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }$
• $B\left({\widehat{\mathsf{\theta }}}^{SPTRE}\right)=B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-d{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }{H}_{q+2,r}(c_{\alpha },{\Delta }^2)$

where $H_{m,r}(.)$ is the non-central cumulative distribution function of F random variable with $(m,r)$ degrees of freedom, $c_{\alpha }$ is the $\alpha -$critical value form the $F-$distribution, and ${\Delta }^2={\displaystyle \frac{{\mathsf{\theta }}^T\mathsf{\theta }}{2}}$ as a non-centrality parameter. For the proof of the above theorem see the Appendix A.

In order to obtain the quadratic risk expressions, we use the quadratic loss function of any RE-type estimator ${\widehat{\mathsf{\theta }}}^{*}$, which is defined for any positive definite $p\times p$ matrix $\mathbf{M}$ as follows:

$$\begin{array}{ccc}\hfill L({\widehat{\mathsf{\theta }}}^{*},\mathsf{\theta })& =& {({\widehat{\mathsf{\theta }}}^{*}-\mathsf{\theta })}^T\mathbf{M}({\widehat{\mathsf{\theta }}}^{*}-\mathsf{\theta })\hfill \end{array}$$

The quadratic risk function of any estimator, denoted by $\Gamma \left({\widehat{\mathsf{\theta }}}^{*}\right)$, and defined as $\Gamma \left({\widehat{\mathsf{\theta }}}^{*}\right)=E\left(L({\widehat{\mathsf{\theta }}}^{*},\mathsf{\theta })\right)$. The quadratic risk expressions are given in the following theorem.

Theorem 2

 

Let $\mathbf{M}$ be a $(p\times p)$ positive definite matrix, then:

$$\begin{array}{c}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)={\tau }^{2}tr\left({\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{A}}_k^T\mathbf{M}\right)+{\mathsf{\theta }}^T({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p)\mathsf{\theta }\end{array}$$

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RRE}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2{\tau }^{2}tr\left({\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{B}}^T{\mathbf{A}}_k^T\mathbf{M}\right)\hfill \\ & -& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p){\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }\hfill \\ & +& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}\right)\mathsf{\theta }\hfill \\ & +& {\tau }^{2}tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{G}}^{-1}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\hfill \\ & +& {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p){\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{REPT}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)H_{q+2,r}(c_{\alpha },{\Delta }^2)\hfill \\ & -& 2\mathsf{\theta }tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }{H}_{q+4,r}(c_{\alpha },{\Delta }^2)\hfill \\ & -& 2\mathsf{\theta }tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}\right)\mathsf{\theta }{H}_{q+2,r}(c_{\alpha },{\Delta }^2)\hfill \\ & +& tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)H_{q+2,r}(c_{\alpha },{\Delta }^2)\hfill \\ & +& {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }{H}_{q+4,r}(c_{\alpha },{\Delta }^2)\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RES}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2\pi tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)E\left(F_{q+2,r}^{-1}\left({\Delta }^2\right)\right)\hfill \\ & -& 2\pi {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }E\left(F_{q+4,r}^{-1}\left({\Delta }^2\right)\right)\hfill \\ & -& 2\pi {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}\right)\mathsf{\theta }E\left(F_{q+2,r}^{-1}\left({\Delta }^2\right)\right)\hfill \\ & +& {\pi }^{2}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k\right)E\left(F_{q+2,r}^{-2}\left({\Delta }^2\right)\right)\hfill \\ & +& {\pi }^2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k\right)\mathsf{\theta }E\left(F_{q+4,r}^{-2}\left({\Delta }^2\right)\right)\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{REPS}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\right)\hfill \\ & & E[\pi F_{q+2,r}^{-1}\left({\Delta }^2\right)+(1-\pi F_{q+2,r}^{-1}\left({\Delta }^2\right))\hfill \\ & & I\left(F_{q+2,r}^{-1}\left({\Delta }^2\right)\le \pi \right)]\hfill \\ & -& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }\hfill \\ & & E[\pi F_{q+4,r}^{-1}\left({\Delta }^2\right)+(1-\pi F_{q+4,r}^{-1}\left({\Delta }^2\right))\hfill \\ & & I\left(F_{q+4,r}^{-1}\left({\Delta }^2\right)\le \pi \right)]\hfill \\ & -& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }\hfill \\ & & E[\pi F_{q+2,r}^{-1}\left({\Delta }^2\right)+(1-\pi F_{q+2,r}^{-1}\left({\Delta }^2\right))\hfill \\ & & I\left(F_{q+2,r}^{-1}\left({\Delta }^2\right)\le \pi \right)]\hfill \\ & +& tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\hfill \\ & & E[\pi F_{q+2}^{-1}\left({\Delta }^2\right)+(1-\pi F_{q+2,r}^{-1}\left({\Delta }^2\right))\hfill \\ & & I\left(F_{q+2,r}^{-1}\left({\Delta }^2\right)\le \pi \right)]\hfill \\ & +& {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }\hfill \\ & & E[\pi F_{q+4,r}^{-1}\left({\Delta }^2\right)+(1-\pi F_{q+4,r}^{-1}\left({\Delta }^2\right))\hfill \\ & & I\left(F_{q+4,r}^{-1}\left({\Delta }^2\right)\le \pi \right)]\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{LSRE}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2d{\tau }^{2}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{G}}^{-1}{\mathbf{A}}^T\right)\hfill \\ & -& 2d{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p){\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }\hfill \\ & +& d^2{\tau }^{2}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{G}}^{-1}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)+2d{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}\right)\mathsf{\theta }\hfill \\ & +& d^2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p){\mathbf{B}}^T{\mathbf{A}}_K^T\right)\mathsf{\theta }\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{SPTRE}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2dtr\left({\mathbf{A}}_k{\mathbf{B}}^T{\mathbf{A}}_k^T\mathbf{M}\right)H_{q+2,r}(c_{\alpha },{\Delta }^2)\hfill \\ & -& 2d{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k{\mathbf{B}}^T{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }{H}_{q+4,r}(c_{\alpha },{\Delta }^2)\hfill \\ & -& 2d{\mathsf{\theta }}^{T}tr\left(\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }{H}_{q+2,r}(c_{\alpha },{\Delta }^2)\hfill \\ & +& d^{2}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)H_{q+2,r}(c_{\alpha },{\Delta }^2)\hfill \\ & +& d^2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }{H}_{q+4,r}(c_{\alpha },{\Delta }^2),\hfill \end{array}$$

For the proof of the above theorem see the Appendix A.

In the following section, we shall enumerate some extant penalization methods to produce a new sub-model from the literature.

5. Some Penalizing Techniques

Penalty estimators are produced as a consequence of simultaneous model selection and parameter estimation processes by applying a penalty to the least squares equation. As a result, model selection and estimating processes are included in penalty techniques. Some of the penalty techniques are provided below.

5.1. Ridge Estimator

The ridge estimator proposed by [1], efficiently addresses the issue of multicollinearity by using a technique known as coefficient shrinkage, which reduces the magnitudes of the coefficients associated with strongly correlated variables. This approach aids in achieving model stability and mitigating the influence of multicollinearity on the estimation of coefficients. It is denoted by ${\widehat{\mathsf{\theta }}}^{Ridge}$, and can obtained by minimizing the penalized residual sum of squares, and it is given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{Ridge}& =& \underset{\mathsf{\theta }}{argmin}\left[{(\mathbf{Y}-\mathbf{X}\mathsf{\theta })}^T(\mathbf{Y}-\mathbf{X}\mathsf{\theta })+{\lambda }_R({\mathsf{\theta }}^T\mathsf{\theta }-c)\right],\hfill \end{array}$$

(13)

where ${\lambda }_R\ge 0$ is a constant known as a biasing parameter.

5.2. LASSO Estimator

The LASSO or the least absolute selection and shrinkage operator proposed by [18]. It is a method for selecting variables and estimating parameters in linear models. In the usual least squares estimation of regression coefficients, the LASSO algorithm employs the $L1$ norm of the vector to define a penalty term. The LASSO estimator, denoted by ${\widehat{\mathsf{\theta }}}^{lasso}$, and given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{lasso}& =& \underset{\mathsf{\theta }}{argmin}\left[{(\mathbf{Y}-\mathbf{X}\mathsf{\theta })}^T(\mathbf{Y}-\mathbf{X}\mathsf{\theta })+{\lambda }_L\sum _{i=1}^p\left|{\theta }_i\right|\right],\hfill \end{array}$$

(14)

where ${\lambda }_L\ge 0$ is tuning parameter to be estimated from the data. However, it is known that the LASSO technique may not be the most ideal approach when dealing with a set of columns in a design matrix that exhibit significant levels of correlation. As a solution, [19] introduced the elastic net (ELNT) method, which combines an $L1$ and $L2$ penalty terms in a linear manner.

5.3. Elastic Net Estimator

It is denoted by ${\widehat{\mathsf{\theta }}}^{ELNT}$ and obtained as follows:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{ELNT}& =& \underset{\mathsf{\theta }}{argmin}\left[{(\mathbf{Y}-\mathbf{X}\mathsf{\theta })}^T(\mathbf{Y}-\mathbf{X}\mathsf{\theta })+{\lambda }_L\sum _{i=1}^p\left|{\theta }_i\right|+{\lambda }_R\sum _{i=1}^p{\theta }_i^2\right],\hfill \end{array}$$

(15)

where ${\lambda }_L$, ${\lambda }_R$ are respectively, the LASSO and the ridge tuning parameters. A variable selection process is said to possess the oracle property when it successfully identifies the correct subset of zero coefficients inside the regression model being examined. Additionally, the estimators of the remaining non-zero coefficients demonstrate consistency and asymptotic normality. The LASSO estimator does not enjoy this property. However, [20] and [21] developed two approaches that possess the oracle property . In the following two subsections, we establish the definitions for these methods.

5.4. SCAD Estimator

The SCAD or smoothly clipped absolute deviation estimator, denoted by ${\widehat{\mathsf{\theta }}}^{SCAD}$ and obtained by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{SCAD}& =& {\displaystyle \underset{\mathsf{\theta }}{argmin}\left\{{(\mathbf{Y}-\mathbf{X}\mathsf{\theta })}^T(\mathbf{Y}-\mathbf{X}\mathsf{\theta })+{\lambda }_S\sum _{i=1}^p{P}_{a,{\lambda }_S}\left({\theta }_i\right)\right\},}\hfill \end{array}$$

(16)

where $P_{a,{\lambda }_S}\left(t\right)$ is a continuous function of t and its derivative is given by:

$$\begin{array}{ccc}\hfill P_{a,{\lambda }_S}^{\prime }\left(t\right)& =& {\lambda }_S\left[I\left(\right|t|\le {\lambda }_S)+{\displaystyle \frac{(a{\lambda }_S-{\left|t\right|)}_{+}}{(a-1){\lambda }_S}}I\left(\right|t|>{\lambda }_S)\right],\hfill \\ & & t\ge 0,\hfill \end{array}$$

where ${\lambda }_S>0$, $a>2$ are the known as tuning parameters. Note that when $a=\infty$, the function $P_{a,{\lambda }_S}\left(t\right)$ is equivalent to $L1$ penalty.

5.5. Adaptive LASSO Estimator

Using adaptive weights on $L1$ penalties on regression coefficients, the adaptive LASSO modifies the LASSO penalty. Theoretically, it has been demonstrated that the adaptive LASSO estimator enjoys the oracle property. It is denoted by ${\widehat{\mathsf{\theta }}}^{Alasso}$, and obtained by:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{Alasso}& =& \underset{\mathsf{\theta }}{argmin}\left[{(\mathbf{Y}-\mathbf{X}\mathsf{\theta })}^T(\mathbf{Y}-\mathbf{X}\mathsf{\theta })+{\lambda }_A\sum _{i=1}^p{\widehat{h}}_i\left|{\theta }_i\right|\right],\hfill \end{array}$$

(17)

where ${\widehat{h}}_i$ is a weight function defined as ${\widehat{h}}_i={\displaystyle \frac{1}{|{\theta }_i^{★}{|}^{\gamma }}}\gamma >0$. The estimator $\widehat{{\theta }_i^{★}}$ is considered to be root-n consistent for ${\theta }_i$. The minimization process for the adaptive LASSO solution does not provide any computing challenges and can be easily solved. One possible choice for $\widehat{{\theta }_i^{★}}$ is OLS estimate of ${\theta }_i$. Now, we will describe some machine learning algorithms that will be used in this manuscript.

6. Machine Learning

The area of regression analysis has been significantly transformed by the advent of machine learning, providing robust methodologies and strategies to extract meaningful insights and achieve precise predictions from intricate data sets. Yet, the emergence of machine learning has expanded the scope of regression analysis, allowing us to use the predictive capacities of algorithms in order to unveil latent patterns and make informed judgments based on data. Three machine learning algorithms are briefly described: Random forest, K-nearest neighbors, and neural network.

6.1. Random Forest

 [27] suggested a technique for creating tree-based classifiers that may be expanded as needed to improve accuracy on both training and hidden data, which is known as random forest. Random Forest is one of the most effective and flexible algorithms that can be used for both classification and regression. It’s a kind of ensemble learning that takes the results of many different decision trees and uses them together to come up with more reliable results. It build several decision trees on distinct subsets of the training data, and let them generate their own predictions. Randomization is introduced via bootstrapping using random data points and training set replacements and selecting a random subset of features for each tree. Later on, [28] discussed the problem of over fitting and getting the most accurate results by building a decision tree-based predictor that stays as accurate as possible on training data and gets more accurate as it gets more complicated. In a Random Forest, the final prediction is derived by aggregating the predictions of all the individual trees, naturally through regression averaging. This ensemble method reduces over-fitting and increases the algorithm’s precision and stability. For more details about the extensions and developments of the algorithm can be found in [29] and [30], among others.

6.2. K-Nearest Neighbors

The K-nearest neighbor (KNN)is an effective machine-learning algorithm used for classification and regression as the Random Forest algorithm. It was initially proposed by  [31], and later on, [32] made modifications to it. The KNN method employs a supervised learning in a non-parametric fashion. In the context of economic forecasting, [33] examined the application of the KNN method, and showed it exhibits greater efficacy compared to alternative methodologies. [34] examined the implementation of the Euclidean distance formula in KNN in comparison to the normalized Euclidean distance, Manhattan, and normalized Manhattan. [35] provided a comprehensive examination of several methodologies employed in Nearest Neighbour classification.

The KNN uses instance-based, non-parametric learning. The method does not rely heavily on assumptions about the underlying data distribution. However, it produces predictions via evaluating similarities between the data points. In regression context, The method computes the average or weighted average of the response value of the K nearest neighbours to predict the response value for the new data point. The KNN algorithm is often used for initial machine learning applications due to its simplicity in comprehension and implementation. Yet, the performance of the system may be influenced by the selection of K and the specific distance measure used. Moreover, the computational cost associated with this approach is high, since it necessitates the calculation of distances for every data point in the data set during the prediction stage.

6.3. Neural Network

A neural network is a computer model that draws inspiration from the anatomical and functional characteristics of the human brain. One of the firsts attempts about this method proposed by [36] who established the foundation for comprehending the ability of basic computational units, resembling neurons, to execute intricate logical operations. Subsequently, other improvements and enhancements have been incorporated into this approach. [37] introduced a detection approach based on a neural networks method for identifying changes in the properties of systems with unknown structures. [38] presented a case study showcasing the effective operation of neural networks in the evaluation of credit risk. [39] provided a thorough examination of clustering neural networks approaches based on competitive learning. [40] provided an overview of academic articles on neural networks that have been published in archival research journals from 1989 to 2000.

The neural network purpose is to effectively process information and acquire knowledge via the analysis of data. Neural networks are composed of linked nodes, which are arranged in layers. These networks have exceptional proficiency in addressing intricate issues, notably in domains like as pattern recognition, picture and audio recognition, natural language processing, and other related fields.

After obtaining a sub-model estimator in addition to the full model estimator, which contains all available variables , it becomes possible to assess the accuracy of the sub model using different criteria. The objective of this study is to identify estimators that are functions of both ${\widehat{\mathsf{\theta }}}^{RE}$ and ${\widehat{\mathsf{\theta }}}^{RER}$ with the purpose of mitigating the potential risks associated with either of these estimators over a wide range of parameter values. This objective can be accomplished by developing the pretest and shrinkage estimators in the subsequent section.

7. Numerical Illustrations

In this section, we will investigate the finite-sample properties of the proposed estimators using Monte Carlo simulation experiments and a real data example.

7.1. Simulation Experiments

The properties of the estimator under investigation and the standards being applied to evaluate the findings will determine how the simulation experiment should be prepared. Following Kibria and Lukman [3], the following equation was used to produce the regressor variables.

$$\begin{array}{ccc}\hfill {\mathbf{X}}_{ij}& =& \sqrt{1-{\rho }^2}{\mathbf{T}}_{ij}+\rho {\mathbf{T}}_{i,p+1},\hfill \\ & & i=1,2,\dots ,n,j=1,2,\dots ,p,\hfill \end{array}$$

(18)

where $T_{ij}$ are independent and identically standard normal random variables, and $\rho$ is the correlation between ${\mathbf{X}}_{ij}$ and ${\mathbf{X}}_{i^{\prime }j}$ for $i\ne i^{\prime }$. The response variable ${\mathbf{Y}}_i$ is then obtained via the following equation:

$$\begin{array}{ccc}\hfill {\mathbf{Y}}_i& =& {\theta }_0+{\theta }_1{\mathbf{X}}_{i1}+{\theta }_2{\mathbf{X}}_{i2}+\dots ++{\theta }_p{\mathbf{X}}_{ip}+{ϵ}_i,\hfill \\ & & i=1,2,\dots ,n,\hfill \end{array}$$

(19)

where ${ϵ}_i$ are $i.i.d$$N(0,{\tau }^2)$, the regression vector is partitioned as $\mathsf{\theta }=({\mathbf{1}}_{\mathbf{p}-\mathbf{q}}^{\mathbf{T}},{\mathbf{0}}_{\mathbf{q}-\mathbf{1}}^{\mathbf{T}},{\Delta }^{*})$, and let ${\Delta }^{*}$ varies over the set $\{0,0.1,0.2,\dots ,2\}$, which represents the degree of deviation from the null hypothesis in (3). In this simulation, we choose $k=0.50$, and set $d=0.10$, ${\tau }^2=1$. The correlation coefficient $\rho$ is chosen to vary over the set $\{0.5,0.75,0.95\}$, $p\in \{10,15,20,25,30\}$, and $q\in \{5,10,15,20,25\}$. We set $\alpha =0.05$ for testing the hypothesis in (3). It was seen that the performance of all estimators had a similar pattern when the values of $\rho ,p$, and q were varied. To save space, we have chosen $\rho \in \{0.75,0.95\}$, $p\in \{10,15\}$, $q\in \{5,10\}$, then we run the simulation for $B=1000$ iterations for $n\in \{50,100\}$.

For each estimator, we compute the mean squared error (MSE) as follows:

$$\begin{array}{ccc}\hfill MS{E}_j\left({\widehat{\mathsf{\theta }}}^{*}\right)& =& {\displaystyle \sum _{i=1}^p{({\theta }_i^{*}-{\theta }_i)}^2,j=1,2,\dots ,B,}\hfill \end{array}$$

(20)

where ${\widehat{\mathsf{\theta }}}^{*}$ is any of the proposed estimators in this study.

For the purpose of comparison, we use the relative efficiency of the mean squared error (RE) with respect to ${\widehat{\mathsf{\theta }}}^{RE}$ which is defined as follows:

$$\begin{array}{ccc}\hfill RE\left({\widehat{\mathsf{\theta }}}^{*},{\widehat{\mathsf{\theta }}}^{RE}\right)& =& {\displaystyle \frac{MSE\left({\widehat{\mathsf{\theta }}}^{RE}\right)}{MSE\left({\widehat{\mathsf{\theta }}}^{*}\right)}},\hfill \end{array}$$

(21)

A number grater than one of the RE indicates the superiority of ${\widehat{\mathsf{\theta }}}^{*}$ over ${\widehat{\mathsf{\theta }}}^{RE}$, and vice versa. Figure 1 shows the graphs for the cases we considered. The following conclusion can be obtained.

• The sub-model estimate ${\widehat{\mathsf{\theta }}}^{RRE}$ consistently beats the all other estimators when the null hypothesis in (3) is true or approximately true. However, its relative efficiency decreases and eventually approaches zero as $\Delta$ increases. Moreover, all estimators outperform the regular ${\widehat{\mathsf{\theta }}}^{RE}$estimator in terms of mean squared error across all values of $\Delta$.
• For all values of $\Delta$, the RE positive shrinkage estimator dominates all other estimators, except when the sub-model is true, in which case the RE sub-model and the pertest estimators outperformed it.
• The relative efficiencies exhibit a consistent pattern when the values of $p,q$, and $\rho$ are held constant for both sample sizes used in this simulation.

7.2. Data Examples

In this section, we will consider a case study discussed by [41] in two different chapters of their book about and also available for free within the "VisCollin"-R-package by [43]. The data originally studied by [42]. His goal was to detect the key soil factors that affect the aerial biomass in the marsh grass Spartina alterniflora in the Cape Fear Estuary of North Carolina at different three locations. At each location, three types of Spartina vegetation areas were sampled, namely devegetated “dead” areas, short" Spartina areas, and "tall" Spartina areas. Five samples of the soil substrate from different sites were collected within each location-vegetation type. These samples were then evaluated for 14 different physico-chemical parameters of the soil monthly for several months, resulting in a total of 45 samples. There are 45 observations in this data set covering the 17 variables, the location loc, area type type, hydrogen sulfide $H_{2}S$, salinity in percentage SAL, ester-hydrolase Eh7, soil acidity in water pH, buffer acidity at pH 6.6 BUF, concentration of following elements: phosphorus P, potassium K, calcium Ca, magnesium Mg, sodium Na, manganese Mn, zinc  Zn, copper Cu, ammonium NH4, and the arial biomass BIO as the repones variable.

One main objective is to employ the set of variables that can accurately predict the response variable. As a first investigation of multicollinearity, we construct the correlation plot among these variables which is given by below in Figure 2.

The plot shows that there are many significant relationships between these variables, and some of these relations are very strong and highly significant with $p-vallue<0.001$. This indicates a multicollinearity exists among theses variables, which can be easily detected by the variance inflation factor plot given below in Figure 3.

The VIF plot shows a serious multicollinearity problem, as some of the value are more than 5 in many cases. Hence, the RE estimator will reduce this problem and produce a better estimates of the parameters. Moreover, the prossed estimators will be an additional improvement to estimate and predict the target response variable BIO.

In the absence of prior knowledge, the limitation on the parameters is established either through the judgment of an expert or by the utilization of existing methodologies for variable selection such as the Akaike information criteria (AIC), forward selection(FW), backward elimination(BE), best subset selection(BS), Bayes information criterion (BIC), or by some penalization algorithms as Lasso, Adaptive lasso and others to produce a sub model. In this example, we at first employ the forward, backward and best subset selection methods to produces a sub model, then obtain the RE, pretest, and shrinkage estimators. Secondly, we apply the random forest, K-nearest neighbors, and the neural network as a machine learning algorithms to compare the prediction error with the seven proposed estimators. The sub models selected by the forward, backward, and best subset selection are summarized in Table 1 below.

In our analysis, we will examine two sub models: the forward sub model, which includes with the intercept the variables pH, Ca, Mg, and Cu. The second one is the backward/best subset that includes with the intercept the variables SAL, K, Zn, Eh7, Mg, Cu, and NH4. The two sub models will be designated as Sub.1 and Sub.2, respectively.

We fit a full model with all available variables and the selected sub-model. The whole model yields an estimated value of ${\widehat{\tau }}^2=s^2=333.9673$. Kibria and Lukman [3] presented several methods for estimating the biasing parameter $\left(k\right)$, but we choose the one with the lowest MSE, which also provided by [23]. The estimated value of $\left(k\right)$ is determined by writing the model in (1) in canonical form as follows:

$$\begin{array}{ccc}\hfill \mathbf{Y}& =& \mathbf{Z}\mathbf{\beta }+\mathbf{ϵ},\hfill \end{array}$$

where $\mathbf{Z}=\mathbf{X}\mathbf{T}$, $\mathbf{\beta }={\mathbf{T}}^T\mathsf{\theta }$, and ${\mathbf{Z}}^T\mathbf{Z}={\mathbf{T}}^T{\mathbf{X}}^T\mathbf{X}\mathbf{T}=\mathbf{D}=diag\left({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p\right)$ are the eigenvalues of $\left({\mathbf{X}}^T\mathbf{X}\right)$, and $\mathbf{T}$ is an orthogonal $(p\times p)$ matrix whose columns are the eigenvectors correspond to the eigenvalues in $\mathbf{D}$. In this case, $\widehat{\mathbf{\beta }}={\left({\mathbf{Z}}^T\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^T\mathbf{Y}={\mathbf{D}}^{-1}{\mathbf{Z}}^T\mathbf{Y}$, and the estimated value of k is given by:

$$\begin{array}{ccc}\hfill {\widehat{k}}_{min.}& =& \underset{i}{min}\left\{{\displaystyle \frac{s^2}{2{\widehat{\beta }}_i^2+(s^2/{\lambda }_i)}}\right\}\hfill \end{array}$$

Using the previous approach, we found ${\widehat{k}}_{min}=0.0000102$. Then the RE-type estimators were calculated. In order to evaluate the estimators’ performance, we implemented a bootstrap method see [24] and calculated the mean squared prediction error in the subsequent manner:

• Select with replacement a sample of size $n=97$ from the data set $K-$times, say $K=2000$.
• Partition each Sample in (1) into separate training and testing sets. The training and testing sets are divided at a ratio of 80% and 20%, respectively. Then, fit a full and sub models using the training data set, and obtain the values of all RE-type estimators.
• Evaluate the predicted response values using each estimator based on the testing data set as follows:

  $$\begin{array}{ccc}\hfill {\widehat{{\mathbf{Y}}_i}}_{,test}& =& {\mathbf{X}}_{i_{,test}}{\widehat{{\mathsf{\theta }}_i}}^{RE*},\hfill \end{array}$$

  where $i=1,2,\dots ,K$, ${\widehat{{\mathbf{Y}}_i}}_{,test}={\widehat{lps{a}_i}}_{,test}$, and ${\mathbf{X}}_{i_{,test}}$ is the matrix of other variables in the model, and ${\widehat{{\mathsf{\theta }}_i}}^{RE*}$ is any of the proposed RE estimators.
• Find the prediction error of each estimator for each sample as follows:

  $$\begin{array}{ccc}\hfill P{E}_i\left({\widehat{{\mathsf{\theta }}_i}}^{RE*}\right)& =& {\displaystyle \frac{1}{n_{i_{test}}}}{\mathbf{e}}_i^T{\mathbf{e}}_i,\hfill \end{array}$$

  where $\mathbf{e}={\widehat{{\mathbf{Y}}_i}}_{,test}-{\overline{{\mathbf{Y}}_i}}_{,train}-\left({\mathbf{X}}_{i_{test}}-{\overline{\mathbf{X}}}_{i,train}\right){\widehat{{\mathsf{\theta }}_i}}^{RE*}$
• Calculate the average prediction error of all estimators as follows:

  $$\begin{array}{ccc}\hfill AVPE\left({\widehat{\mathsf{\theta }}}^{RE*}\right)& =& {\displaystyle \frac{1}{K}\sum _{i=1}^{K}P{E}_i\left({\widehat{\mathsf{\theta }}}^{RE*}\right)}\hfill \end{array}$$
• Finally, calculate the relative efficiency of the prediction error with respect to ${\widehat{\mathsf{\theta }}}^{RE}$ as follows:

  $$\begin{array}{ccc}\hfill REPE\left({\widehat{\mathsf{\theta }}}^{RE*},{\widehat{\mathsf{\theta }}}^{RE}\right)& =& {\displaystyle \frac{AVPE\left({\widehat{\mathsf{\theta }}}^{RE}\right)}{AVPE\left({\widehat{\mathsf{\theta }}}^{R{E}^{*}}\right)}}\hfill \end{array}$$

We ran the program with several values of d and observed no noticeable effect when altering it, so we set it to $d=0.50$. The findings shown in Table 2 align with the outcomes of the simulations discussed in the preceding subsection.

The analysis shows that the sub model estimator ${\widehat{\mathsf{\theta }}}^{LSRE}$ , produced by Sub.1, outperforms all other estimators in terms of prediction error. Similarly, ${\widehat{\mathsf{\theta }}}^{RRE}$, in the case of Sub.2, exhibits the highest level of performance. This implies that in both sub models, the predictors that are eliminated from the full model are either irrelevant or nearly unrelated to the response variable. Furthermore, all RE-type estimators have superior performance in comparison to the usual RE estimator, as they all have a relative efficiency greater than one.

Next, we will employ various penalization and three machine learning algorithms to analyze the prostate cancer data. Our objective is to determine the prediction error. We first scale all variables including the repones (BIO), then apply the penalization or machine learning algorithms to avoid any differences of the variable units. The results of our investigation are summarized in the following table.

As shown in the Table 3 above, every machine learning algorithm outperformed the traditional RE estimator. However, the performance of the ridge, Lasso, and SCAD penalization methods was inferior to that of the RE estimator. This discrepancy may be attributed, in part, to the presence of muticolinearity among the predictor variables.

Upon careful examination of the numerical data, it becomes apparent that the prediction error relative efficiencies differ from one method to another. To gain a deeper comprehension of the range of values and potential outliers in our prediction, let us examine the associated prediction error using the following box plots, and based on 1000 replications.

The box plot in Figure 4 clearly illustrates the distribution of our data set, and with further analysis, it becomes apparent that there are clues of possible outliers. The extended whiskers and isolated data points beyond the usual range indicate variability and occurrences that diverge from the general pattern. Furthermore, using a shrinkage technique to the (RE) estimator reveals a significant effect on the suppression of outliers. Applying shrinkage not only improves the estimate process but also helps us identify and reduce the impact of outliers, resulting in a more flexible and dependable representation of the underlying data structure.

8. Conclusion

This article introduces the pretest and shrinkage of the new ridge estimators when there is multicollinearity among predictor variables in a multiple linear regression model. Since the pretest and shrinkage technique relies on prior information, we developed a linear hypothesis to assess the significance of certain regressor variables in the regression model. We then combine the result in our estimation process. Subsequently, we benefit from Kibria Lukman’s idea by implementing the novel estimating approach to mitigate the extent of multicollinearity.

In order to compare, we implemented several penalization approaches along with three machine learning algorithms. We assessed the comparative efficacy of each method’s prediction error in relation to the usual RE estimator.

Our findings demonstrated that the utilization of the shrinkage estimation technique resulted in a significant enhancement in the mean squared error when applied to simulated data under various configurations of the correlation coefficient among the predictor variables. Furthermore, these findings were supported by a real data illustration in which we employed the relative efficiency of the prediction errors as a benchmark for comparison.

The suggested estimation techniques can be expanded to encompass various different types of regression models, such as Poisson, logistic, beta regression models, among others. In addition, the analysis can also consider the case of high-dimensional data, evaluating the effectiveness of the suggested estimators in comparison to penalization methods and machine learning algorithms.

Appendix A.1

Proof of Theorem 1

 

1. Note that ${\widehat{\mathsf{\theta }}}^{RE}=(\mathbf{G}+k{\mathbf{I}}_p)(\mathbf{G}-k{\mathbf{I}}_p)\widehat{\mathsf{\theta }}={\mathbf{A}}_k\widehat{\mathsf{\theta }}$, where ${\mathbf{A}}_k=(\mathbf{G}+k{\mathbf{I}}_p)(\mathbf{G}-k{\mathbf{I}}_p)$, so

$$\begin{array}{ccc}\hfill B\left({\widehat{\mathsf{\theta }}}^{RE}\right)& =& E\{{\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\}\hfill \\ & =& E\{{\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\}\hfill \\ & =& ({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }.\hfill \end{array}$$

2. Note that

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\theta }}}^{RRE}& =& {\mathbf{A}}_k{\widehat{\mathsf{\theta }}}^R\hfill \\ & =& {\mathbf{A}}_k(\widehat{\mathsf{\theta }}-{\mathbf{G}}^{-1}{\mathbf{R}}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}\mathbf{R}\widehat{\mathsf{\theta }})\hfill \\ & =& {\mathbf{A}}_k({\mathbf{I}}_P-{\mathbf{G}}^{-1}{\mathbf{R}}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}\mathbf{R})\widehat{\mathsf{\theta }}\hfill \\ & =& {\mathbf{A}}_k({\mathbf{I}}_p-\mathbf{B})\widehat{\mathsf{\theta }},\mathbf{B}={\mathbf{G}}^{-1}{\mathbf{R}}^T{\left(\mathbf{R}{\mathbf{G}}^{-1}{\mathbf{R}}^T\right)}^{-1}\mathbf{R},\hfill \end{array}$$

Hence,

$$\begin{array}{ccc}\hfill B\left({\widehat{\mathsf{\theta }}}^{RRE}\right)& =& E\{{\mathbf{A}}_k({\mathbf{I}}_p-\mathbf{B})\widehat{\mathsf{\theta }}-\mathsf{\theta }\}\hfill \\ & =& \left({\mathbf{A}}_k({\mathbf{I}}_p-\mathbf{B})-{\mathbf{I}}_p\right)\mathsf{\theta }\hfill \\ & =& \left(({\mathbf{A}}_k-{\mathbf{I}}_p)-{\mathbf{A}}_k\mathbf{B}\right)\mathsf{\theta }\hfill \\ & =& B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }.\hfill \end{array}$$

For the proofs of (3)-(7), we can relay on the following general proof.

Note that, ${\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE}={\mathbf{A}}_k\widehat{\mathsf{\theta }}-{\mathbf{A}}_k({\mathbf{I}}_p-\mathbf{B})\widehat{\mathsf{\theta }}={\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}$, then by equation (6),

$$\begin{array}{ccc}\hfill B\left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}\right)& =& {\displaystyle E\left\{{\widehat{\mathsf{\theta }}}^{RE}-({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE})f\left({\varphi }_n\right)-\mathsf{\theta }\right\}}\hfill \\ & =& {\displaystyle E\{{\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\}-E\left\{({\widehat{\mathsf{\theta }}}^{RE}-{\widehat{\mathsf{\theta }}}^{RRE})f\left({\varphi }_n\right)\right\}}\hfill \\ & =& {\displaystyle B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}f\left({\varphi }_n\right)\right\}}\hfill \end{array}$$

Using the result of Theorem (1) Appendix B of [44] , we have

${\displaystyle E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}f\left({\varphi }_n\right)\right\}={\mathbf{A}}_k\mathbf{B}\mathsf{\theta }E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)}$. Therefore

$$\begin{array}{ccc}\hfill B\left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}\right)& =& B\left({\widehat{\mathsf{\theta }}}^{RE}\right)-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\hfill \end{array}$$

Using the suitable function ${\displaystyle f\left({\varphi }_n\right)}$, the bias functions of ${\widehat{\mathsf{\theta }}}^{REPT},{\widehat{\mathsf{\theta }}}^{RES}$, and ${\widehat{\mathsf{\theta }}}^{REPS}$ can be easily obtained.

Proof of Theorem 2

 

1.

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)& =& E\left\{{\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)}^T\mathbf{M}\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)\right\}\hfill \\ & =& tr\left\{\mathbf{M}E\left\{\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right){\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)}^T\right\}\right\}\hfill \end{array}$$

Let $E_1=E\left\{\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right){\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)}^T\right\}$, then

$$\begin{array}{ccc}\hfill E_1& =& E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right)}^T\right\}\hfill \\ & =& E\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-{\mathbf{A}}_k\mathsf{\theta }+{\mathbf{A}}_k\mathsf{\theta }-\mathsf{\theta }\right)\hfill \\ & & {\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-{\mathbf{A}}_k\mathsf{\theta }+{\mathbf{A}}_k\mathsf{\theta }-\mathsf{\theta }\right)}^T\}\hfill \\ & =& E\{\left[{\mathbf{A}}_k\left(\widehat{\mathsf{\theta }}-\mathsf{\theta }\right)+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }\right]\hfill \\ & & {\left[{\mathbf{A}}_k\left(\widehat{\mathsf{\theta }}-\mathsf{\theta }\right)+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }\right]}^T\}\hfill \\ & =& E\left\{{\mathbf{A}}_k(\widehat{\mathsf{\theta }}-\mathsf{\theta }){(\widehat{\mathsf{\theta }}-\mathsf{\theta })}^T{\mathbf{A}}_k^T\right\}\hfill \\ & +& 2E\left\{{\mathbf{A}}_k(\widehat{\mathsf{\theta }}-\mathsf{\theta }){\widehat{\mathsf{\theta }}}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right\}\hfill \\ & +& E\left\{({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right\}\hfill \\ & =& E_{11}+2E_{12}+E_{13},\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill E_{11}& =& E\left\{{\mathbf{A}}_k(\widehat{\mathsf{\theta }}-\mathsf{\theta }){(\widehat{\mathsf{\theta }}-\mathsf{\theta })}^T{\mathbf{A}}_k^T\right\}\hfill \\ & =& {\mathbf{A}}_{k}Var\left(\widehat{\mathsf{\theta }}\right){\mathbf{A}}_k^T\hfill \\ & =& {\tau }^2{\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{A}}_k^T,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E_{12}& =& 2E\left\{{\mathbf{A}}_k(\widehat{\mathsf{\theta }}-\mathsf{\theta }){\widehat{\mathsf{\theta }}}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right\}\hfill \\ & =& 2{\mathbf{A}}_{k}E\left\{(\widehat{\mathsf{\theta }}-\mathsf{\theta })\right\}{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p),\hfill \\ & & \mathrm{but}\widehat{\mathsf{\theta }}\mathrm{is}\mathrm{an}\mathrm{unbiased}\mathrm{estimator}\mathrm{for}\mathsf{\theta }\hfill \\ & =& {\mathbf{0}}_{\mathbf{p}},\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E_{13}& =& E\left\{({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right\}\hfill \\ & =& ({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p).\hfill \end{array}$$

Hence,

$$\begin{array}{ccc}\hfill E_1& =& {\tau }^2{\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{A}}_k^T+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p),\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)& =& tr\left\{\mathbf{M}\left[{\tau }^2{\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{A}}_k^T+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right]\right\}\hfill \\ & =& {\tau }^{2}tr\left({\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{A}}_k^T\mathbf{M}\right)+{\mathsf{\theta }}^T({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p)\mathsf{\theta }.\hfill \end{array}$$

2.

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RRE}\right)& =& E\left\{{\left({\widehat{\mathsf{\theta }}}^{RRE}-\mathsf{\theta }\right)}^T\mathbf{M}\left({\widehat{\mathsf{\theta }}}^{RRE}-\mathsf{\theta }\right)\right\}\hfill \\ & =& tr\left\{\mathbf{M}E\left\{\left({\widehat{\mathsf{\theta }}}^{RRE}-\mathsf{\theta }\right){\left({\widehat{\mathsf{\theta }}}^{RRE}-\mathsf{\theta }\right)}^T\right\}\right\}\hfill \end{array}$$

Let $E_2=E\left\{\left({\widehat{\mathsf{\theta }}}^{RRE}-\mathsf{\theta }\right){\left({\widehat{\mathsf{\theta }}}^{RRE}-\mathsf{\theta }\right)}^T\right\}$, then

$$\begin{array}{ccc}\hfill E_2& =& E\left\{\left({\mathbf{A}}_k({\mathbf{I}}_p-\mathbf{B})\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\left({\mathbf{A}}_k({\mathbf{I}}_p-\mathbf{B})\widehat{\mathsf{\theta }}-\mathsf{\theta }\right)}^T\right\}\hfill \\ & =& E\left\{\left(({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta })-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }\right){\left(({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta })-{\mathbf{A}}_k\mathbf{B}\mathsf{\theta }\right)}^T\right\}\hfill \\ & =& E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right)}^T\right\}\hfill \\ & -& E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\right\}\hfill \\ & -& E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta })}^T\right\}\hfill \\ & +& E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\right\}\hfill \\ & =& E_{21}-E_{22}-E_{23}+E_{24},\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill E_{21}& =& E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right)}^T\right\}\hfill \\ & =& E_1,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E_{22}& =& E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\right\}\hfill \\ & =& E\left\{{\mathbf{A}}_k\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T-\mathsf{\theta }{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\right\}\hfill \\ & =& {\mathbf{A}}_{k}E\left(\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T\right){\mathbf{B}}^T{\mathbf{A}}_k^T-\mathsf{\theta }{\mathsf{\theta }}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & =& {\mathbf{A}}_k\left({\tau }^2{\mathbf{G}}^{-1}+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right){\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & -& \mathsf{\theta }{\mathsf{\theta }}^T{\mathbf{B}}^T{\mathbf{A}}_k^T,\hfill \end{array}$$

similarly,

$$\begin{array}{ccc}\hfill E_{23}& =& E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta })}^T\right\}\hfill \\ & =& {\mathbf{A}}_k\mathbf{B}\left({\tau }^2{\mathbf{G}}^{-1}+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right){\mathbf{A}}_k^T\hfill \\ & -& {\mathbf{A}}_k\mathbf{B}\mathsf{\theta }{\mathsf{\theta }}^T,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E_{24}& =& E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\right\}\hfill \\ & =& {\mathbf{A}}_k\mathbf{B}E\left(\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T\right){\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & =& {\mathbf{A}}_k\mathbf{B}\left({\tau }^2{\mathbf{G}}^{-1}+({\mathbf{A}}_k-{\mathbf{I}}_p)\mathsf{\theta }{\mathsf{\theta }}^T({\mathbf{A}}_k^T-{\mathbf{I}}_p)\right){\mathbf{B}}^T{\mathbf{A}}_k^T.\hfill \end{array}$$

By little algebra, we can add the terms $E_{12},E_{22},E_{23},E_{24}$ to get the expression $E_2$, so

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RRE}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2{\tau }^{2}tr\left({\mathbf{A}}_k{\mathbf{G}}^{-1}{\mathbf{B}}^T{\mathbf{A}}_k^T\mathbf{M}\right)\hfill \\ & -& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p){\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }\hfill \\ & +& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}\right)\mathsf{\theta }\hfill \\ & +& {\tau }^{2}tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{G}}^{-1}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\hfill \\ & +& {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}({\mathbf{A}}_k-{\mathbf{I}}_p)\mathbf{M}({\mathbf{A}}_k^T-{\mathbf{I}}_p){\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }.\hfill \end{array}$$

Now, for the proofs of (3)-(7), we can do the general proof for any ${\widehat{\mathsf{\theta }}}^{RE-Shrinkage}$, then use the suitable function ${\displaystyle f\left({\varphi }_n\right)}$ to get the quadric risk expression as we did in the previous theorem’s proof.

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}\right)& =& E\{{\left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}-\mathsf{\theta }\right)}^T\hfill \\ & & \mathbf{M}\left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}-\mathsf{\theta }\right)\}\hfill \\ & =& tr\{\mathbf{M}\left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}-\mathsf{\theta }\right)\hfill \\ & & {\left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}-\mathsf{\theta }\right)}^T\}\hfill \\ & =& tr\{\mathbf{M}E\left({\displaystyle \left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)-{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}f\left({\varphi }_n\right)}\right)\hfill \\ & & {\left({\displaystyle \left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)-{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}f\left({\varphi }_n\right)}\right)}^T\}.\hfill \end{array}$$

Let ${\displaystyle E_3=E\left\{\left[\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)-{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}f\left({\varphi }_n\right)\right]{\left[\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)-{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}f\left({\varphi }_n\right)\right]}^T\right\}}$, then

$$\begin{array}{ccc}\hfill E_3& =& E\left[\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)({\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)}^T\right]\hfill \\ & -& {\displaystyle E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^{T}f\left({\varphi }_n\right)\right\}}\hfill \\ & -& {\displaystyle E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}\left({\widehat{\mathsf{\theta }}}^T{\mathbf{A}}_k^T-{\mathsf{\theta }}^T\right)f\left({\varphi }_n\right)\right\}}\hfill \\ & +& {\displaystyle E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T{f}^2\left({\varphi }_n\right)\right\}}\hfill \\ & =& E_{31}+E_{32}+E_{33}+E_{34},\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill E_{31}& =& E\left[\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)({\left({\widehat{\mathsf{\theta }}}^{RE}-\mathsf{\theta }\right)}^T\right]\hfill \\ \\ & =& E_1,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E_{32}& =& {\displaystyle E\left\{\left({\mathbf{A}}_k\widehat{\mathsf{\theta }}-\mathsf{\theta }\right){\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^{T}f\left({\varphi }_n\right)\right\}}\hfill \\ & =& {\displaystyle {\mathbf{A}}_{k}E\left\{\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^{T}f\left({\varphi }_n\right)\right\}{\mathbf{B}}^T{\mathbf{A}}_k^T-\mathsf{\theta }E\left\{{\widehat{\mathsf{\theta }}}^{T}f\left({\varphi }_n\right)\right\}{\mathbf{B}}^T{\mathbf{A}}_k^T.}\hfill \end{array}$$

Using the result of Theorems (1)-(3) Appendix B of [44] , we have

$$\begin{array}{ccc}\hfill E_{32}& =& {\mathbf{A}}_k\left[E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)+E\left({\displaystyle f\left(F_{q+4,r}\left({\Delta }^2\right)\right)}\right)\mathsf{\theta }{\mathsf{\theta }}^T\right]\hfill \\ & & {\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & -& \mathsf{\theta }\left({\mathsf{\theta }}^{T}E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\right){\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & =& {\mathbf{A}}_k{\mathbf{B}}^T{\mathbf{A}}_k^{T}E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & +& {\mathbf{A}}_k\mathsf{\theta }{\mathsf{\theta }}^T{\mathbf{B}}^T{\mathbf{A}}_k^{T}E\left({\displaystyle f\left(F_{q+4,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & -& \mathsf{\theta }{\mathsf{\theta }}^T{\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & & E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right),\hfill \end{array}$$

similarly,

$$\begin{array}{ccc}\hfill E_{33}& =& {\displaystyle E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T{f}^2\left({\varphi }_n\right)\right\}}\hfill \\ \\ & =& {\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^{T}E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & +& {\mathbf{A}}_k\mathbf{B}\mathsf{\theta }{\mathsf{\theta }}^T{\mathbf{A}}_k^{T}E\left({\displaystyle f\left(F_{q+4,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & -& {\mathbf{A}}_k\mathbf{B}\mathsf{\theta }{\mathsf{\theta }}^{T}E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E_{34}& =& {\displaystyle E\left\{{\mathbf{A}}_k\mathbf{B}\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{\mathbf{B}}^T{\mathbf{A}}_k^T{f}^2\left({\varphi }_n\right)\right\}}\hfill \\ & =& {\displaystyle {\mathbf{A}}_k\mathbf{B}E\left\{\widehat{\mathsf{\theta }}{\widehat{\mathsf{\theta }}}^T{f}^2\left({\varphi }_n\right)\right\}{\mathbf{B}}^T{\mathbf{A}}_k^T}\hfill \\ & =& {\mathbf{A}}_k\mathbf{B}E\left({\displaystyle f^2\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right){\mathbf{B}}^T{\mathbf{A}}_k^T\hfill \\ & +& {\mathbf{A}}_k\mathbf{B}\mathsf{\theta }{\mathsf{\theta }}^T{\mathbf{B}}^T{\mathbf{A}}_k^{T}E\left({\displaystyle f^2\left(F_{q+4,r}\left({\Delta }^2\right)\right)}\right).\hfill \end{array}$$

Using little algebra and combining the four expected values, we will get the value of $E_3$ . Now,

$$\begin{array}{ccc}\hfill \Gamma \left({\widehat{\mathsf{\theta }}}^{RE-Shrinkage}\right)& =& \Gamma \left({\widehat{\mathsf{\theta }}}^{RE}\right)\hfill \\ & -& 2tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & -& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}{\mathbf{A}}_k^T\mathbf{M}\right)\mathsf{\theta }E\left({\displaystyle f\left(F_{q+4,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & -& 2{\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}\right)\mathsf{\theta }E\left({\displaystyle f\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & +& tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k\right)E\left({\displaystyle f^2\left(F_{q+2,r}\left({\Delta }^2\right)\right)}\right)\hfill \\ & +& {\mathsf{\theta }}^{T}tr\left({\mathbf{A}}_k\mathbf{B}\mathbf{M}{\mathbf{B}}^T{\mathbf{A}}_k^T\right)\mathsf{\theta }\hfill \\ & & E\left({\displaystyle f^2\left(F_{q+4,r}\left({\Delta }^2\right)\right)}\right).\hfill \end{array}$$

Therefore, in order to obtain the quadric risk expression for the proofs of (3)-(7), we employ the appropriate function ${\displaystyle f\left({\varphi }_n\right)}$, just as we did in the proof of the previous Theorem (1).