Estimation of ill-conditioned models using penalized sums of squares of the residuals

Given the following multiple linear regression model for $n$ observations and $p$ independent variables (constant term included):

where $\mathbf{u}$ is the random disturbance (which is assumed to be spherical), the most common estimation method is Ordinary Least Squares (OLS), which consists of minimizing with respect to $\boldsymbol{\beta}$ the sum of squares of the residuals given by the expression:

However, there are several situations where the objective function to be minimized has been modified, as in the case of the method of Restricted Least Squares (RSM). In this case, the modification seeks to incorporate certain information into the model in order to obtain a more efficient estimation.

It is also common to modify it in order to achieve greater stability in the estimates of model (1). In this sense, according to O’Driscoll and Ramirez [34], it was Tikhonov [50] who first proposed the idea of regularization to estimate ill-conditioned problems by introducing additional information.

Other examples with alternative objective functions are the ridge estimator (Hoerl and Kennard [18, 19]) or LASSO regression (Tibshirani [48, 49]), widely used to mitigate the approximate multicollinearity existing in the linear model.

In the former, $SCR(\boldsymbol{\beta})$ is minimized subject to the constraint $||\boldsymbol{\beta}||_{2}^{2}=\sum\limits_{i=1}^{p}\beta_{i}^{2}=\boldsymbol{\beta}^{t}\boldsymbol{\beta}\leq r$ where $r$ is fixed, and implies minimizing with respect to $\boldsymbol{\beta}$ the Lagrangian:

The solution obtained corresponds to the following expression:

where $\mathbf{I}$ denoted the identity matrix of appropriate dimensions. In addition, its variance-covariance matrix is:

In the second, $SCR(\boldsymbol{\beta})$ is minimized subject to the constraint $||\boldsymbol{\beta}||_{1}^{2}=\sum\limits_{i=1}^{p}|\beta_{i}|\leq s$ where $s$ is fix which implies minimizing with respect to $\boldsymbol{\beta}$ the Lagrangian:

The difference between the two lies in the fact that due to the behavior of the $l_{1}$ norm, LASSO contracts and selects variables, while ridge only contracts. It can also be seen that both methods can be viewed as particular situations of the general case in which the penalty added in the objective function is considered to be given by the $l_{q}$ norm, i.e., it is of the form $||\boldsymbol{\beta}||_{q}^{2}=\sum\limits_{i=1}^{p}|\beta_{i}|^{q}$. Other shrinkage estimators can be found in Liu [22, 23], Gao and Liu [12] or Liu and Gao [24].

More recently, the methodology proposed by Zou and Hastie [55], known as elastic-net, has emerged as a particularly noteworthy approach. This technique allows the contraction of estimators and the automatic selection of variables by combining the penalty factors (in the sum of squares of the residuals) of the ridge and LASSO regressions. It is also particularly useful when the number of independent variables is much higher than the number of observations and to deal with the problem of multicollinearity (see, for example, Bingzhen et al. [4] or Zou and Zhang [56]).

Since its publication, this work has been received¹¹1 Figures as of May 8, 2024 in Scopus and Google Scholar. 11820 citations according to SCOPUS and 20777 according to Google Scholar (2260 in 2021, 2345 in 2022 and 2100 in 2023). These figures show that the development of techniques to deal with multicollinearity and variable selection is currently an interesting econometric research topic, probably due to the fact that are essential tools in high-dimensional data analysis (see Ding et al. [5]).

Within this line of research, Wang et al. [52] proposes a novel method called average least squares-centered penalized regression (ALPR). This penalized regression method shrinks the parameters toward the weighted average OLS estimator. This paper provides a theoretical development to obtain the mean square error of the proposed estimator, as well as a simulation study comparing this estimator with the estimators provided by OLS and ridge regression. Some issues to be considered with respect to this paper are the following:

• •

  If the OLS estimates have coefficients with different signs, it is possible that the weighted average is close to zero, so there is not much difference between the penalized estimator they propose and the ridge regression.

• •

  Shrinking the penalized estimators around the weighted average OLS estimator is just as arbitrary as shrinking them around zero (ridge regression).

• •

  ALPR shrinks all the parameters to the same single value. Again, this is a completely arbitrary question without any theoretical justification.

• •

  The intercept is not penalized by ALPR, then the approximate non-essential multicollinearity (see footnote in Example 1 about essential and non-essential multicollinearity) is ignored.

Recently, Lukman et al. [25] have integrated ALPR with the principal component dimension reduction method to improve its performance in terms of root mean square error.

Alternatively, the present work proposes to minimize the Lagrangian with respect to $\boldsymbol{\beta}$ the Langragian:

where $\boldsymbol{\alpha}$ is a vector that collects the estimates of the slope in all simple linear regressions of the form $\mathbf{y}=\alpha_{1}\cdot\mathbf{1}+\alpha_{i}\cdot\mathbf{X}_{i}+\mathbf{v}$ with $i=2,\dots,p$, where $\mathbf{1}$ denotes the constant term.

The choice of the parameter around which to shrink the penalized estimators, $\boldsymbol{\alpha}$, is intended to overcome the drawbacks of ALPR highlighted above (see section 2). In addition, the obtained estimator is analyzed in depth (sections 3 to 7), providing a consolidated methodology.

Specifically, the work is structured as follows: Section 2 justifies the usefulness of the proposal made to mitigate the approximate multicollinearity existing in the model (1) and achieve stability in the estimation of the coefficients of the independent variables, while in section 3 we obtain the expression of the estimator provided by minimizing the lagrangian given in (4) and its main characteristics: trace, norm, variance-covariance matrix, goodness-of-fit and mean square error. To highlight in this section the Monte Carlo simulation performed to analyze the behavior of the mean square error of the proposed penalized estimator. The implications of this methodology in multicollinearity mitigation are discussed in section 4. The above results are used in section 5 to propose procedures for determining the ideal value of the penalty parameter, $k$. An interpretation of the same is also provided. From these values, the section 6 shows how to perform inference from bootstrap methods. Section 7 analyzes the stability of the estimates obtained from the proposed penalized estimator and section 8 illustrates the results obtained by means of an example with real data. The paper ends with the section 9, where the main results obtained are highlighted. Finally, the code utilized in R [35] to generate the results presented in this study is accessible on GitHub, specifically at https://github.com/rnoremlas/Penalized-estimator.

Taking into account the following notions of partial and full effect given in Novales [31]:

• •

  What is the impact on $\mathbf{y}$ of a unit variation in $\mathbf{X}_{i}$ if the other explanatory variables in the model did not vary? Answer: partial effect (multiple linear regression).

• •

  What is the total impact on $\mathbf{y}$ of a unit variation in $\mathbf{X}_{i}$ if the rest of the explanatory variables vary as would be expected given the correlations observed between them throughout the sample? Answer: total effect (simple linear regression).

Both effects should coincide (see example 1) if the ceteris paribus, i.e., when one variable varies the rest remain constant, were really verified.

Developing the previous example from a theoretical point of view:

• •

  Estimates of simple linear regressions $\mathbf{y}=\alpha_{1}\cdot\mathbf{1}+\alpha_{i}\cdot\mathbf{X}_{i}+\mathbf{v}$ with $i=2,\dots,p$ are given by the expression:

  $$\widehat{\alpha}_{i}=\frac{cov(\mathbf{y},\mathbf{X}_{i})}{var(\mathbf{X}_{i})},\quad i=2,\dots,p.$$

• •

  Assuming orthogonality in the model (1), the OLS estimator responds to the expression:

  	$\displaystyle\widehat{\boldsymbol{\beta}}$	$\displaystyle=$	$\displaystyle\left(\mathbf{X}^{t}\mathbf{X}\right)^{-1}\mathbf{X}\mathbf{y}$
  		$\displaystyle=$	$\displaystyle\left(\begin{array}[]{cccc}n&0&\cdots&0\\ 0&\sum\limits_{j=1}^{n}X_{2j}^{2}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\sum\limits_{j=1}^{n}X_{kj}^{2}\\ \end{array}\right)^{-1}\left(\begin{array}[]{c}\sum\limits_{j=1}^{n}y_{j}\\ \sum\limits_{j=1}^{n}X_{2j}y_{j}\\ \vdots\\ \sum\limits_{j=1}^{n}X_{kj}y_{j}\\ \end{array}\right)^{-1}=\left(\begin{array}[]{c}\overline{\mathbf{y}}\\ \frac{cov(\mathbf{y},\mathbf{X}_{2})}{var(\mathbf{X}_{2})}\\ \vdots\\ \frac{cov(\mathbf{y},\mathbf{X}_{k})}{var(\mathbf{X}_{k})}\\ \end{array}\right)$

  since, for $i=2,\dots,p$, it is verified that $cov(\mathbf{y},\mathbf{X}_{i})=\frac{1}{n}\sum\limits_{j=1}^{n}X_{ij}y_{j}$ and $var(\mathbf{X}_{i})=\frac{1}{n}\sum\limits_{j=1}^{n}X_{ij}^{2}$ finding that $\sum\limits_{j=1}^{n}X_{2j}=\dots=\sum\limits_{j=1}^{n}X_{pj}=0$ by the orthogonality condition between $\mathbf{X}_{i}$ and the intercept.

For this reason, in the expression (4) it will be considered that:

that is, $\boldsymbol{\alpha}$ is a known fixed vector.

In short, in the restriction imposed on the sum of squares of the penalized residuals represented in the lagrangian of the expression (4) it will be considered that $\boldsymbol{\alpha}$ will take the value given in (6). In this way, the purpose is to make the estimates of the coefficients of the multiple linear regression model as close as possible to the estimates obtained in the corresponding simple linear regressions.

In this section we obtain the estimator that verifies the above conditions, as well as its norm, variance-covariance matrix, goodness-of-fit and mean square error. We also consider the augmented model whose OLS estimation leads to the penalized estimator obtained. In the latter case, the effect of this augmented model on the detection of the degree of multicollinearity existing once the proposed estimator has been applied is analyzed.

With this purpose, instead of working with the expression (4), expression (7) will be minimized where previously a parameter $h$ has been introduced as follows:

note that for $h=0$ the Hoerl and Kennard estimator would be obtained. [18, 19]. That is, the ridge estimator is a particular case of the one presented here.

While it makes sense that $h\in[0,\ 1]$, it will be considered only that $h\in\{0,1\}$ (for further details see subsection 3.1.1).

From expression (8) the penalized estimator can be obtained by OLS estimation of the augmented model $\mathbf{y}_{A}=\mathbf{X}_{A}\cdot\boldsymbol{\beta}+\mathbf{u}_{A}$, where:

However, the variance-covariance matrices would be different since:

Similarly, the goodness-of-fit would also change since the dependent variable is different, $\mathbf{y}\not=\mathbf{y}_{A}$.

It should be noted that the design matrix of the augmented model, $\mathbf{X}_{A}$, coincides with that of the augmented model of the ridge estimator (see Marquardt [26]), then all operations based on them will coincide. Consequently, papers presented by García et al. [13] and Salmerón et al. [39] to calculate, respectively, the extension of the Variance Inflation Factor (VIF) and the Condition Number (CN) to the ridge estimator are applicable to the penalized estimator. Thus, more specifically:

• •

  García et al. [13] states that to obtain an extension of the VIF, $VIF(i,k)$ with $i=2,\dots,p$, with the desirable properties of continuity (for $k=0$ the calculated VIF coincides with that obtained in OLS), decreasing monotonicity as $k$ increases (since it is assumed that when applying the ridge estimate the degree of approximate multicollinearity decreases) and that a value greater than 1 is always obtained (since this is the minimum value of the VIF); data must be standardized⁴⁴4 To standardize a set of data, center the data (subtract its mean) and divide it by its standard deviation multiplied by the square root of the number of observations. and subsequently the model is augmented. In other words, to calculate the VIF for the ridge estimator, the following matrix must be used:

  $$\mathbf{x}_{A}=\left(\begin{array}[]{c}\mathbf{x}\\ \sqrt{k}\cdot\mathbf{I}\end{array}\right)_{(n+p)\times p},$$

  where the matrix $\mathbf{x}$ is the matrix resulting from the standardization of the columns of $\mathbf{X}$.

• •

  Salmerón et al. [39] it is recommended that to obtain an extension of the CN, $CN(k)$, that verifies the above desirable properties (continuity at $k=0$, decreasing monotonicity at $k$ and being greater than 1 for any value of $k$) it has to be calculated from the expression:

  $$CN(k)=\sqrt{\frac{\xi_{max}+k}{\xi_{min}+k}},$$

  where$\xi_{min}$ and $\xi_{max}$ are, respectively, the minimum and maximum eigenvalues of $\mathbf{X}_{T}^{t}\mathbf{X}_{T}$ where $\mathbf{X}_{T}$ denotes the matrix $\mathbf{X}$ transformed. These authors recommend that the transformation should be typing⁵⁵5 To typify a set of data it is necessary to center it (subtract its mean) and divide it by its standard deviation. the data, since in such a case the results obtained are similar to those obtained from the VIF. However, this transformation ignores nonessential type multicollinearity (see, for example, Salmerón et al. [40, 41]), so this transformation is not going to be used in this paper. In accordance with the recommendations of Belsley et al. [1] or Belsley [3, 2] the eigenvalues of the matrix will be calculated. $\mathbf{X}_{T}^{t}\mathbf{X}_{T}$ where the transformation used is that of unit length⁶⁶6 To transform a set of data into unit length, divide it by the square root of the sum of each element squared. Note that if the data have zero mean (a property that is obtained by centering them), it is verified that this transformation coincides with the standardization. .

In order to obtain the value of the above expressions it is necessary to use a concrete value of the penalty parameter, $k$. For this reason, some possible criteria for determining a concrete value of the parameter, as well as its interpretation, are proposed below.

Although there are various attempts to address inference in the ridge estimator (see, e.g., Obenchain [32, 33], Halawa and El Bassiouni [17], Gökpınar and Ebegil [16], Sengupta and Sowell [45] or Vanhove [51]), in this paper we will focus on the use of bootstrap methods (see, for example, Efron [6, 7, 8, 9], Efron and Gong [10] or the a review of the above-mentioned works presented by Efron and Tibshirani [11]).

Thus, given a fixed value of $k$ obtained by one of the methods proposed in subsection 5.2 and a value of $\boldsymbol{\alpha}$, the following steps will be taken:

1. (i)

   Generate randomly and with replacement $m$ subsamples of equal size to the original one. The value of $m$ must be large.

2. (ii)

   For each sub-sample above, the statistic of interest $theta$ is calculated. Therefore, we have $m$ values for this statistic: $\theta_{1},\dots,\theta_{m}$.

3. (iii)

   Calculate the mean, $\overline{\theta}=\frac{1}{m}\sum\limits_{i=1}^{m}\theta_{i}$, and sample quasi-variance, $\sigma^{2}_{\theta}=\frac{1}{m-1}\sum\limits_{i=1}^{m}(\theta_{i}-\overline{\theta})^{2}$, for the $m$ values obtained in previous step. Both measures can be used to obtain the approximation of a confidence interval by using expression $\overline{\theta}\pm 1.96\cdot\sigma_{\theta}$.
   Similarly, $\overline{\theta}$ can be used as an approximation of the point estimator of $\theta$.

4. (iv)

   Another confidence region can be obtained alternatively by considering the lower and upper percentiles, 0.025 and 0.975, as the lower and upper extremes respectively, of the $m$ values calculated in the second step: $[P_{0.025}(\theta_{1},\dots,\theta_{m}),P_{0.975}(\theta_{1},\dots,\theta_{m})]$.

EWhat is interesting in this work are the cases in which $\theta$ is equal to $\widehat{\boldsymbol{\beta}}(k,1)$ or $GoF(k,1)$.