g.ridge: An R Package for Generalized Ridge Regression for Sparse and High-Dimensional Linear Models

Abstract

Ridge regression is one of the most popular shrinkage estimation methods for linear models. Ridge regression effectively estimates regression coefficients in the presence of high-dimensional regressors. Recently, a generalized ridge estimator was suggested that involved generalizing the uniform shrinkage of ridge regression to non-uniform shrinkage; this was shown to perform well in sparse and high-dimensional linear models. In this paper, we introduce our newly developed R package “g.ridge” (first version published on 7 December 2023) that implements both the ridge estimator and generalized ridge estimator. The package is equipped with generalized cross-validation for the automatic estimation of shrinkage parameters. The package also includes a convenient tool for generating a design matrix. By simulations, we test the performance of the R package under sparse and high-dimensional settings with normal and skew-normal error distributions. From the simulation results, we conclude that the generalized ridge estimator is superior to the benchmark ridge estimator based on the R package “glmnet”. Hence the generalized ridge estimator may be the most recommended estimator for sparse and high-dimensional models. We demonstrate the package using intracerebral hemorrhage data.

1. Introduction

In linear regression, the least squares estimator (LSE) may not be suitable to estimate regression coefficients if the number of regressors $p$ is higher than the sample size $n$ (i.e., $p>n$). Ridge regression [1,2] is an effective alternative for high-dimensional ($p>n$) data, and is widely employed in such data encountered in genetics [3,4], medical studies [5], electromagnetic signals [6], grain yields [7], and others.

Hoerl and Kennard [1] originally proposed the ridge estimator to reduce the problem in multicollinearity. Later, the ridge estimator was naturally adopted to high-dimensional ($p>n$) problems [8,9] as a way to avoid overfitting phenomena.

Ridge regression is a shrinkage estimator that shrinks all regression coefficients uniformly toward zero [1,9,10,11]. This approach is particularly suitable for modeling microarrays or single-nucleotide polymorphism (SNP) data, where many coefficients are close to zero (sparse models). For instance, Cule et al. [12] applied the ridge estimator to high-dimensional SNP data and performed significance testing for selecting an informative subset of SNPs. Similar applications of ridge regression to high-dimensional data have been performed by Whittaker et al. [13], Cule and De Lorio [14], and Yang and Emura [15].

Unlike ordinary ridge regression that shrinks all regression coefficients uniformly, the generalized ridge regression allows different degrees of shrinkage under multiple shrinkage parameters [1,16,17,18]. Loesgen [18] demonstrated that multiple shrinkage parameters in the generalized ridge estimator arise naturally from prior uncertainty (variance) for regression coefficients. Shen et al. [19] and Hofheinz and Frisch [20] adopted this idea and proposed a generalized ridge estimator for high-dimensional SNPs where multiple shrinkage parameters are estimated preliminarily. Readers are referred to a review paper of Yüzbaşı et al. [21] for more information about generalized ridge regression.

However, multiple shrinkage parameters are difficult to estimate for high-dimensional cases. To overcome this difficulty, Yang and Emura [15] suggested reducing the number of shrinkage parameters for the case of $p>n$ in their formulation of a generalized ridge estimator. The idea behind their approach was to assign two different weights (1 or 0.5) for shrinkage parameters via preliminary tests [22,23,24,25,26]. While this approach has sound statistical performance in sparse and high-dimensional models, no software packages were implemented for the generalized ridge estimator.

This paper introduces our R package “g.ridge” (https://cran.r-project.org/package=g.ridge, accessed on 5 February 2024) that implements both the ridge estimator of Hoerl and Kennard [1] and the generalized ridge estimator of Yang and Emura [15]. The package is equipped with generalized cross-validation for the automatic estimation of shrinkage parameters. The package also includes a convenient tool for generating a design matrix. By simulations, we test the performance of the R package with sparse and high-dimensional models, and normal and skew-normal distributions for error terms. We conclude that the generalized ridge estimator is superior to the benchmark ridge estimator based on the R package “glmnet”. Hence the generalized ridge estimator may be the most recommended estimator for sparse and high-dimensional models. We illustrate the package via the intracerebral hemorrhage dataset.

The remainder of this paper is structured as follows. Section 2 reviews the ridge and generalized ridge estimators. Section 3 introduces the proposed R package. Section 4 tests the package via simulations. Section 5 gives a real data example. Section 6 concludes the paper.

2. Ridge Regression and Generalized Ridge Regression

This section introduces the ridge regression method proposed by Hoerl and Kennard [1] and the generalized ridge regression proposed by Yang and Emura [15].

2.1. Linear Regression

Consider a linear regression model ($\mathit{y}=X\mathit{\beta }+\mathit{\epsilon }$) where

$$\mathit{y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right], \hspace{0.33em}X=\left[\begin{array}{c}{\mathit{x}}_1^{\mathrm{T}}\\ ⋮\\ {\mathit{x}}_n^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{ccc}{x}_{11}& \cdots & x_{1p}\\ ⋮& \ddots & ⋮\\ x_{n1}& \cdots & x_{np}\end{array}\right], \hspace{0.33em}\mathit{\beta }=\left[\begin{array}{c}{\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right], \hspace{0.33em}\mathit{\epsilon }=\left[\begin{array}{c}{\epsilon }_1\\ ⋮\\ {\epsilon }_n\end{array}\right];$$

$X$ is a design matrix, ${\mathit{x}}_i^{\mathrm{T}}=(x_{i1},\dots ,x_{ip})$ is the transpose of a vector ${\mathit{x}}_i$, $\mathit{\beta }$ is a vector of regression coefficients, and $\mathit{\epsilon }$ is a vector of errors with $E\left[\mathit{\epsilon }\right]=\hspace{0.33em}\mathbf{0}$. In some cases, we assume $\mathrm{C}\mathrm{o}\mathrm{v}\left[\mathit{\epsilon }\right]={\sigma }^2{I}_n$, where $I_n$ is the identity matrix of dimension $n$. The assumption of the covariance is necessary to obtain the standard error (SE) of regression coefficients and test their significance. Assume that $X$ is standardized such that ${\sum }_{i=1}^n{x}_{ij}=0$ and ${\sum }_{i=1}^n{x}_{ij}^2=c$ for $j=1,...,p$, where $c$ is $n$ or $n-1$. Also, we assume that $X$ does not include the intercept term (see Section 3.3 for details). These settings for $X$ are usually imposed in ridge regression [9].

If $X^{\mathrm{T}}X$ is invertible (non-singular), the LSE is defined as

$${\widehat{\mathit{\beta }}}^{\mathrm{L}\mathrm{S}\mathrm{E}}=(X^{\mathrm{T}}X{)}^{-1}{X}^{\mathrm{T}}\mathit{y}.$$

Clearly, the LSE is not available when $X^{\mathrm{T}}X$ is singular, especially when $p>n$.

2.2. Ridge Regression

The ridge estimator is an alternative to LSE, and can be computed even when $p>n$. Hoerl and Kennard [1] originally defined the ridge regression estimator as

$$\widehat{\mathit{\beta }}\left(\lambda \right)=(X^{\mathrm{T}}X+\lambda I_n{)}^{-1}{X}^{\mathrm{T}}\mathit{y},$$

where $\lambda >0$ is called the shrinkage parameter. A diagonal matrix $\lambda I_n$ added to $X^{\mathrm{T}}X$ makes all the components of $\widehat{\mathit{\beta }}\left(\lambda \right)$ shrink toward zero. Theorem 4.3 of Hoerl and Kennard [1] shows that there exists $\lambda >0$, such that the ridge estimator yields a strictly smaller mean squared error (MSE) compared to that of the LSE.

Golub et al. [8] suggested choosing $\lambda$ on the basis of the predicted residual error sum of squares, or Allen’s PRESS [17] statistics. Golub et al. [8] proposed the rotation-invariant version of the PRESS statistics, and called it the generalized cross-validation (GCV) criterion [8]. The GCV function defined by Golub et al. [8] is

$$V\left(\lambda \right)={\frac{1}{n}||\{I_n-A(\lambda \})\mathit{y}|{|}^2/\left[\frac{1}{n}\mathrm{Tr}\left\{I_n-A\left(\lambda \right)\right\}\right]}^2.$$

(1)

where $A\left(\lambda \right)=X(X^{\mathrm{T}}X+\lambda I_n)X^{\mathrm{T}}$ is a “hat matrix”. The GCV estimator of $\lambda$ is defined as

$$\widehat{\lambda }=\underset{\lambda \ge 0}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}V\left(\lambda \right).$$

Therefore, the resultant ridge estimator is

$$\widehat{\mathit{\beta }}\left(\widehat{\lambda }\right)=(X^{\mathrm{T}}X+\widehat{\lambda }{I}_n{)}^{-1}{X}^{\mathrm{T}}\mathit{y}.$$

The GCV theorem [8] justifies the above ridge estimator under the $p>n$ setup. Note that the GCV criterion is different from the cross-validation criterion that is adopted in the R package “glmnet”. There are many other criteria for choosing $\lambda$ [14,27,28,29,30,31], most of which are not applicable or justified for the $p>n$ setup. Therefore, we will adopt the GCV criterion for the following discussions.

2.3. Generalized Ridge Regression

Yang and Emura [15] suggest relaxing the uniform shrinkage to the non-uniform shrinkage to yield the generalized ridge regression. For this purpose, Yang and Emura [15] replaced the identity matrix $I_n$ with the diagonal matrix $\widehat{W}\left(\Delta \right)$ (defined later), and proposed

$$\widehat{\mathit{\beta }}(\lambda ,\Delta )=\{X^{\mathrm{T}}X+\lambda \widehat{W}(\Delta ){\}}^{-1}{X}^{\mathrm{T}}\mathit{y},$$

where $\lambda >0$ is a shrinkage parameter and $\Delta \ge 0$ is called the “threshold” parameter. The diagonal components of $\widehat{W}\left(\Delta \right)$ were suggested to be larger values (stronger shrinkage) for the components $\mathit{\beta }$ closer to zero, yielding $\widehat{W}(\Delta )=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{\widehat{w}}_1\right(\Delta ),\hspace{0.33em}...,{\widehat{w}}_p(\Delta \left)\right\}$, where

$${\widehat{w}}_j\left(\Delta \right)=\left\{\begin{array}{cc}1/2 & \mathrm{if}\hspace{1em}{z}_j\ge \Delta ,\\ 1& \mathrm{if}\hspace{1em}{z}_j<\Delta \end{array}\right.$$

where $z_j={\widehat{\beta }}_j^0/\mathrm{S}\mathrm{D}\left({\widehat{\mathit{\beta }}}^0\right)$, and $\mathrm{S}\mathrm{D}\left({\widehat{\mathit{\beta }}}^0\right)={\left\{{\sum }_{j=1}^p({\widehat{\beta }}_j^0-{\sum }_{j=1}^p{\widehat{\beta }}_j^0/p\hspace{0.33em}{)}^2/\left(p-1\right)\right\}}^{1/2}$ for $j=1,\hspace{0.33em}...,\hspace{0.33em}p$, and ${\widehat{\mathit{\beta }}}^0=({\widehat{\beta }}_1^0,...,{\widehat{\beta }}_p^0{)}^{\mathrm{T}}$, defined as ${\widehat{\beta }}_j^0=X_j^{\mathrm{T}}\mathit{y}/X_j^{\mathrm{T}}{X}_j$, where $X_j$ is the $j$-th column of $X$. Note that ${\widehat{\mathit{\beta }}}^0$ is called the “compound covariate estimator” [32].

The optimal value of $(\lambda ,\Delta )$ is estimated by the modified GCV function defined as

$$V(\lambda ,\Delta )={\frac{1}{n}||\{I_n-A(\lambda ,\Delta )\}\mathit{y}|{|}^2/\left[\frac{1}{n}\mathrm{Tr}\{I_n-A(\lambda ,\mathsf{\Delta }\left)\right\}\right]}^2,$$

(2)

where $A(\lambda ,\Delta )=X\{X^{\mathrm{T}}X+\lambda \widehat{W}(\Delta ){\}}^{-1}{X}^{\mathrm{T}}$. Then, estimators $(\widehat{\lambda },\widehat{\Delta })$ are defined as

$$(\widehat{\lambda },\widehat{\Delta })=\underset{\lambda \ge 0,\hspace{0.33em}\Delta \ge 0}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}V(\lambda ,\Delta ).$$

The computation of $(\widehat{\lambda },\widehat{\Delta })$ is not difficult. Given $\Delta$, the function $V(\lambda ,\Delta )$ is continuous in $\lambda$, and hence it is easily minimized using any optimization scheme, such as the R function “optim(.)”, to obtain $\widehat{\lambda }\left(\Delta \right)$. In a sparse ($\mathit{\beta }\approx \mathbf{0}$) and high-dimensional ($p>n$) model, the histogram of ${\widehat{\beta }}_j^0/\mathrm{S}\mathrm{D}\left({\widehat{\mathit{\beta }}}^0\right)$, $j=1,\hspace{0.33em}...,\hspace{0.33em}p$, is well approximated by $N\left(\mathrm{0,1}\right)$. This implies that $\left|{\widehat{\beta }}_j^0\right|/\mathrm{S}\mathrm{D}\left({\widehat{\mathit{\beta }}}^0\right)$ falls in the range $[0, 3]$ with nearly 99.73%, and hence a search range of $\Delta \in [0, 3]$ suffices. Since $V\left(\widehat{\lambda }\right(\Delta ),\Delta )$ is discontinuous in $\Delta$, a grid search was suggested on a grid $D=\{0, 3/100,..., 300/100\}$.

Finally, the generalized ridge estimator is defined as

$$\widehat{\mathit{\beta }}(\widehat{\lambda },\widehat{\Delta })=\{X^{\mathrm{T}}X+\widehat{\lambda }\widehat{W}(\widehat{\Delta }){\}}^{-1}{X}^{\mathrm{T}}\mathit{y}.$$

Also, the error variance can be estimated by

$${\widehat{\sigma }}^2={‖\mathit{y}-X\widehat{\mathit{\beta }}\left(\widehat{\lambda },\widehat{\Delta }\right)‖}^2/\nu ,$$

where $\nu =\mathrm{Tr}\{I_n-A(\widehat{\lambda },\widehat{\Delta }){\}}^2$ and $A(\widehat{\lambda },\widehat{\Delta })=X\{X^{\mathrm{T}}X+\widehat{\lambda }\widehat{W}(\widehat{\Delta }){\}}^{-1}{X}^{\mathrm{T}}$.

2.4. Significance Test

Ridge and generalized ridge estimators provide methods to test the following null hypothesis:

$$H_{0j}:\hspace{0.33em}{\beta }_j=0 \mathrm{vs}. H_{1j}:\hspace{0.33em}{\beta }_j\ne 0,$$

for $j=1,\hspace{0.33em}...,\hspace{0.33em}p$. One can perform significance tests, allowing one to access p-values of all the $p$ regressors [12,14,15]. Since the resultant tests are similar between the ridge and generalized ridge estimators, we will explain the tests based on the generalized ridge estimator below.

Let ${\widehat{\beta }}_j(\widehat{\lambda },\widehat{\Delta })$ be the $j$-th component of $\widehat{\mathit{\beta }}(\widehat{\lambda },\widehat{\Delta })$. As in Cule et al. [12] and Yang and Emura [15], we define a Z-value of $Z_j={\widehat{\beta }}_j(\widehat{\lambda },\widehat{\Delta })/\mathrm{S}\mathrm{E}\left\{\right(\widehat{\lambda },\widehat{\Delta }\left)\right\}$, where $\mathrm{S}\mathrm{E}\left\{{\widehat{\beta }}_j\right(\widehat{\lambda },\widehat{\Delta }\left)\right\}$ is the square root of the $j$-th diagonal of

$$\mathrm{C}\mathrm{o}\mathrm{v}\left\{\widehat{\mathit{\beta }}\right(\widehat{\lambda },\widehat{\Delta }\left)\right\}={\widehat{\sigma }}^2\{X^{\mathrm{T}}X+\widehat{\lambda }\widehat{W}(\widehat{\Delta })\hspace{0.33em}{\}}^{-1}{X}^{\mathrm{T}}X\{X^{\mathrm{T}}X+\widehat{\lambda }\widehat{W}\left(\widehat{\Delta }\right){\}}^{-1}.$$

We define the p-value as $P_j=2\{1-\mathsf{\Phi }\left(\left|Z_j\right|\right)\}$, where $\mathsf{\Phi }$(.) is the cumulative distribution function of $N\left(0, 1\right)$. One can then select a subset of regressors by specifying a significance level.

3. R Package: g.ridge

We implemented the methods of Section 2 in the R package “g.ridge”. This section explains two functions: “X.mat” for generating data and “g.ridge(.)“ for performing regression.

3.1. Generating Data

The design matrix $X$ and response vector $\mathit{y}$ are necessary to perform linear regression (Section 2.1). Therefore, following Section 5 of Yang and Emura ([15], p. 6093), we implemented a convenient R function that generates $X$ with three independent blocks:

$$\begin{gathered} X=\left[\begin{array}{c}{\mathit{x}}_1^{\mathrm{T}}\\ ⋮\\ {\mathit{x}}_n^{\mathrm{T}}\end{array}\right]=\left[\left.\begin{array}{ccc}{x}_{11}& \cdots & x_{1q}\\ ⋮& \ddots & ⋮\\ x_{n1}& \cdots & x_{nq}\end{array}\right| \begin{array}{ccc}{x}_{1,q+1}& \cdots & x_{1,q+r}\\ ⋮& \ddots & ⋮\\ x_{n,q+1}& \cdots & x_{n,q+r}\end{array} \left|\begin{array}{ccc}{x}_{1,q+r+1}& \cdots & x_{1p}\\ ⋮& \ddots & ⋮\\ x_{n,q+r+1}& \cdots & x_{np}\end{array}\right.\right], \\ = \mathbf{[} \mathrm{First} \mathrm{block}\mathbf{|} \mathrm{Second} \mathrm{block}\mathbf{|} \mathrm{Third} \mathrm{block}\mathbf{]}, \end{gathered}$$

(3)

where the first block consists of $q\ge 0$ correlated regressors (correlation = 0.5), the second block consists of $r\ge 0$ correlated regressors (correlation = 0.5), and the third block consists of $p-q-r>0$ independent regressors. That is,

$$\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}(x_{ij},\hspace{0.33em}{x}_{ik})=\left\{\begin{array}{cc} 0.5& \mathrm{if}\hspace{1em}j,\hspace{0.33em}k\in \{1,...,q\},\hspace{1em}\hspace{1em}\hspace{1em}\\ 0.5& \mathrm{if}\hspace{1em}j,\hspace{0.33em}k\in \{q+1,...,q+r\},\\ 0& \mathrm{otherwise}.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

The design matrix mimics the “gene pathway” for two blocks of correlated gene expressions often used in simulation studies [33,34,35,36]. The values $q=0$ and $r=0$ give a design matrix with all independent columns.

The marginal distributions of $p$ regressors follow $N\left(\mathrm{0,1}\right)$, which is achieved by

$${\mathit{x}}_i^{\mathrm{T}}=\left(\frac{1}{\sqrt{2}}\left(z_{i1}+u_i\right), ... ,\frac{1}{\sqrt{2}}\left(z_{iq}+u_i\right)\left| \frac{1}{\sqrt{2}}\left(z_{i,q+1}+v_i\right), ...,\frac{1}{\sqrt{2}}\left(z_{i,q+r}+v_i\right) \right| z_{i,q+r+1}, ...,z_{ip}\right),$$

where $z_{i1},\hspace{0.33em}...,\hspace{0.33em}{z}_{ip}$ and $u_i,\hspace{0.33em}{v}_i$ all independently follow $N\left(\mathrm{0,1}\right)$ for $i=1,...,n$.

Figure 1 shows an example for generating design matrices using “X.mat(.)”. One can simply input $n$, $p$, $q$, and $r$ into “X.mat(.)” by noting the constraint $q+r<p$.

3.2. Performing Regression

The ridge and generalized ridge estimators can be computed by the R function “g.ridge(.)” whose input requires a design matrix $X$ and a response vector $\mathit{y}$. Specifically, the command “g.ridge(X, Y, method = “HK”, kmax = 500)” can calculate $\widehat{\mathit{\beta }}\left(\widehat{\lambda }\right)$, where “HK” stands for “Hoerl and Kennard”, and “kmax = 500” means the range $0\le \lambda \le 500$ for estimating $\lambda$. The output displays $\widehat{\lambda }\in [0, 500]$ and $\widehat{\mathit{\beta }}\left(\widehat{\lambda }\right)$. Similarly, the command “g.ridge(X, Y, method = “YE”, kmax = 500)” can calculate $\widehat{\mathit{\beta }}(\widehat{\lambda },\widehat{\Delta })$, where “YE” stands for “Yang and Emura”. The output also displays the plot of $V\left(\lambda \right)$ against $\lambda \in [0, 500]$ and its minimizer $\widehat{\lambda }\in [0, 500]$ (the plot of $V(\lambda ,\widehat{\Delta })$ for the generalized ridge).

The R function “g.ridge(.)” can also compute the SE of $\widehat{\mathit{\beta }}$, Z-value and p-value for significance tests (Section 2.4), and the estimate of the error variance, ${\widehat{\sigma }}^2$ (Section 2.3). As in the typical practice of ridge regression, one should use the centered responses “Y-mean(Y)” rather than “Y” (this will be explained in Section 3.3). Also, if “X” were not generated by “X.mat”, the scaled design matrix “scale(X)” would be recommended rather than “X” itself.

Figure 2 shows the code and output for the ridge estimator. The output shows $\widehat{\lambda }=31.66314$, $\widehat{\mathit{\beta }}\left(\widehat{\lambda }\right)=\left(0.581485036,\dots 0.083260387\right)$, SE, Z-value, p-value, and ${\widehat{\sigma }}^2=1.778618$. The output also displays the GCV function $V\left(\lambda \right)$ against $\lambda \in [0, 200]$, showing its minimum at $\widehat{\lambda }=31.66314$. In the code, we changed “kmax” from the default value (500) to 200 for a better visualization of the plot. In many cases, users will need to try different values of “kmax” to reach a desirable plot for the GCV function. Appendix A explains another function, “GCV(.)”, that is used for internally computing the GCV function $V\left(\lambda \right)$.

The output of the generalized ridge estimator is similar to that of the ridge estimator. The only difference is an additional parameter estimate $\widehat{\Delta }\in [0, 3]$, shown at “delta$”.

3.3. Technical Remarks on Centering and Standardization

We assume that $X$ is standardized and does not have a column for an intercept term (Section 2.1). If $X$ is generated by “X.mat”, it is already standardized, and hence there is no need to perform standardization. However, in other cases, $X$ has to be standardized, e.g., by the R command “scale(X)”. By assumption, the design matrix $X$ does not include a column of ones (for an intercept) since one can always use the reduced model $\mathit{y}-({\sum }_{i=1}^n{y}_i/n)\mathbf{1}=X\mathit{\beta }+\mathit{\epsilon }$ derived from the intercept model $\mathit{y}={\beta }_0\mathbf{1}+X\mathit{\beta }+\mathit{\epsilon }$ (${\beta }_0$ is estimated by ${\sum }_{i=1}^n{y}_i/n$), where $\mathbf{1}$ is a vector of ones. This means that ${\widehat{\beta }}_0={\sum }_{i=1}^n{y}_i/n$ is applied to estimate ${\beta }_0$.

4. Simulations

We conducted simulation studies to test our R package “g.ridge”. The goals of the simulations are to show that the generalized ridge estimators in our package exhibit (a) superior performance over the ridge estimator in the R function “glmnet(.)” and (b) sound performance under normally and non-normally distributed errors (a skew-normal distribution will be considered). The Supplementary Materials include the R code to reproduce the results of the simulation studies.

4.1. Simulation Settings

We generated $X$ by “X.mat” with $n=100$, $p\in \{50,\hspace{0.33em}100,\hspace{0.33em}150,\hspace{0.33em}200\}$, and $q=r=10$ based on Equation (3). Given $X$, we set $\mathit{y}=X\mathit{\beta }+\mathit{\epsilon }$, with $\mathit{\beta }$ defined by a sparse model of

$$\mathit{\beta }=(\stackrel{q}{\overbrace{b/q, ...,b/q}},\hspace{0.33em}\stackrel{r}{\overbrace{d/r, ...,d/r}},\hspace{0.33em}\stackrel{p-q-r}{\overbrace{0,...,0}}{)}^{\mathrm{T}},$$

for four cases: (I) $b=d=5$; (II) $b=d=10$; (III) $b=5$, $d=-5$; (IV) $b=10$, $d=-10$. Errors $\mathit{\epsilon }$ were generated independently from the normal distribution or the skew-normal distribution [37]; both distributions had a mean of zero and standard deviation of one, and the skew-normal distribution had a slant parameter of ten (alpha = 10 in the R function “rsn(.)”). Figure 3 shows the remarkable difference between the two distributions. The skew-normal distribution was not previously examined in the simulation setting of Yang and Emura [15].

For a given dataset $(X,\mathit{y})$, we computed $\widehat{\mathit{\beta }}$, which can be (i) the ridge estimator by the R function “g.ridge(.)”, (ii) the generalized ridge estimator by the R function “g.ridge(.)”, or (iii) the ridge estimator by the R function “glmnet(.)”. Based on 500 replications (on random errors $\mathit{\epsilon }$), the performances of the three proposed estimators were compared by the total mean squared error (TMSE) criterion, defined as

$$\mathrm{T}\mathrm{M}\mathrm{S}\mathrm{E}\left(\widehat{\mathit{\beta }}\right)=\mathrm{E}\left[{‖\widehat{\mathit{\beta }}-\mathit{\beta }‖}^2\right]=\mathrm{E}\left[{\sum }_{j=1}^p{\left({\widehat{\beta }}_j-{\beta }_j\right)}^2\right],$$

where $\mathrm{E}[.]$ was performed by a Monte Carlo average over 500 replications.

4.2. Simulation Results

Table 1. The total mean squared error (TMSE) of the three estimators: (i) the ridge by “g.ridge(.)”, (ii) the generalized (g-) ridge by “g.ridge(.)”, and (iii) the ridge by “glmnet(.)”. The TMSE is computed by a Monte Carlo average over 500 replications.

Error Distribution	Regression Coefficients	$\mathit{p}$	(i) ridge	(ii) g-ridge	(iii) glmnet
Normal	(I) $b=d=5$	50	0.463	0.385	0.306
		100	0.950	0.682	2.182
		150	1.146	0.658	1.996
		200	1.520	0.920	2.199
	(II) $b=d=10$	50	0.855	0.681	0.545
		100	2.151	1.562	8.688
		150	3.008	1.482	7.904
		200	4.929	2.687	8.691
	(III) $b=5$ and $d=-5$	50	0.602	0.539	0.388
		100	0.990	0.628	2.025
		150	1.219	0.703	2.132
		200	1.589	0.953	2.226
	(IV) $b=10$ and $d=-10$	50	1.541	1.290	0.737
		100	2.398	1.580	8.046
		150	3.231	1.614	8.434
		200	4.651	2.770	8.804
Skew-normal	(I) $b=d=5$	50	0.440	0.361	0.294
		100	0.957	0.670	2.182
		150	1.162	0.678	2.000
		200	1.500	0.910	2.197
	(II) $b=d=10$	50	0.821	0.655	0.527
		100	2.285	1.705	8.691
		150	3.021	1.509	7.905
		200	4.883	2.673	8.686
	(III) $b=5$ and $d=-5$	50	0.576	0.519	0.376
		100	0.974	0.622	2.029
		150	1.233	0.721	2.137
		200	1.582	0.949	2.243
	(IV) $b=10$ and $d=-10$	50	1.504	1.273	0.720
		100	2.449	1.508	8.054
		150	3.224	1.616	8.453
		200	4.618	2.731	8.860

Note: We set the sample size $n=100$ for all cases.

compares the TMSE among the three estimators: (i) the ridge by “g.ridge(.)”, (ii) the generalized ridge by “g.ridge(.)”, and (iii) the ridge by “glmnet(.)”. We see that the generalized ridge estimator is superior to the ridge estimator since the former has smaller TMSE values for all cases. Also, the generalized ridge estimator is superior to the “glmnet(.)” estimator for cases of $p\in \left\{100, 150, 200\right\},$ while this is not the case for $p=50$. This conclusion holds consistently across the four parameter settings (I)–(IV) and two error distributions (normal and skew-normal). In conclusion, the generalized ridge estimator in the proposed R package seems to be the most recommended estimator for data with sparse and high-dimensional settings.

5. Data Analysis

This section analyzes a real dataset to illustrate the generalized ridge estimator in the proposed package. We retrospectively analyzed a dataset from patients with intracerebral hemorrhage who were hospitalized at Saiseikai Kumamoto Hospital, Kumamoto city, Japan. The data collection was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of Saiseikai Kumamoto Hospital on 29 March 2023. Saiseikai Kumamoto Hospital is a regional tertiary hospital that serves patients with stroke in southern Japan and provides acute care in a comprehensive stroke care unit.

Outcome variables ($\mathit{y}$) are the changes in blood volume (in mL) from the initial value to the follow-up value, which were measured by CT scans. Excluding patients with less than 10 mL blood volume at the initial CT scans and other reasonable inclusion/exclusion criteria, we arrive at $n=172$ patients. Regressor variables ($X$) consist of $p=35$ continuous measurements, including histological variables (e.g., age), vital signs (e.g., blood pressure, mmHG; respiratory rate, times/minute), and blood measurements (e.g., albumin, g/dL; Gamma-glutamyl transpeptidase (Gamma-GT), IU/L; lactate dehydrogenase, IU/L; prothrombin time, second; blood platelet count, 10⁻³/µL; C-reactive protein, mg/dL). The responses and regressors are centered and standardized before fitting the linear model, as explained in Section 2.1 and Section 3.3.

Table 2. Fitted results for estimated regression coefficients (only with p-value < 0.05) sorted by p-values applied to a dataset on patients with intracerebral hemorrhage.

	Ridge			Generalized Ridge
	${\widehat{\mathit{\beta }}}_{\mathit{j}}$	SE	p-Value	${\widehat{\mathit{\beta }}}_{\mathit{j}}$	SE	p-Value
Lactate dehydrogenase	0.122	0.047	0.008	0.145	0.055	0.008
Gamma-GT	0.116	0.048	0.016	0.143	0.056	0.010
Respiratory rate	−0.120	0.052	0.020	−0.140	0.059	0.018
Prothrombin time	0.077	0.036	0.031	0.083	0.040	0.038
Blood platelet count	−0.100	0.049	0.040	−0.114	0.056	0.044
C-reactive protein	None	None	>0.05	0.112	0.057	0.049

shows estimates for regression coefficients to compare their values between the ridge and generalized ridge estimators. Due to the limited space, we only show regression coefficients with p < 0.05 (5% significance). An interesting difference is found in the number of significant regressors: six regressors for the generalized ridge estimator and five regressors for the ridge estimator. That is, C-reactive protein was significant by the generalized ridge estimator and was not by the ridge estimator. A study conducted in a large population in China reported that a high C-reactive protein level was an independent risk factor for severe intracerebral hemorrhage [38]. The selection of C-reactive protein as a variable associated with increased hematoma volume may reflect the tendency of patients with severe intracerebral hemorrhage to have easier hematoma expansion. Hence, detecting the significance by the generalized ridge estimator provides an important advantage for this data example. Besides C-reactive protein, other estimated coefficients are convincing for the generalized ridge estimator. For instance, the most significant regressor is lactate dehydrogenase (also the most significant for the ridge estimator). This variable was previously reported as an important predictor of hematoma expansion [39].

Figure 4 displays scatter plots for centered responses $\mathit{y}-({\sum }_{i=1}^n{y}_i/n)\mathbf{1}$ against the predictors $X\widehat{\mathit{\beta }}$ based on the ridge estimator and generalized ridge estimator. While both predictors captured the variability in the changes in blood volume (response variables), this was not fully explained by the predictors since some residual errors remain. Figure 5 depicts the residuals $\mathit{y}-({\sum }_{i=1}^n{y}_i/n)\mathbf{1}-X\widehat{\mathit{\beta }}$ constructed by the generalized ridge estimator. We observe that the residuals are asymmetric around zero. The asymmetry in the errors does not yield any problem as the generalized ridge estimator was shown to work well under asymmetric errors (Section 4).

Therefore, we have demonstrated that the generalized ridge estimator in the proposed R package provides statistically and biologically sound conclusions in real data analysis.

6. Conclusions

This paper introduced the new R package “g.ridge” that can perform the ridge and generalized ridge regression methods. We showed that the generalized ridge estimator in the proposed package performs better than the widely used ridge estimator in the “glmnet” package. Furthermore, the ridge and generalized ridge estimators in the proposed package can deal with asymmetric error distributions. With an abundance of sparse and high-dimensional data [36,40,41,42,43,44] and asymmetrically distributed data [45,46,47,48], the proposed package may provide a reliable statistical tool for statisticians and data scientists.

We now compare our R package with some existing R packages for ridge regression. The “bigRR” package (https://cran.r-project.org/package=bigRR, accessed on 5 February 2024) implements the generalized ridge estimator developed by Shen et al. [19]. This package uses a different estimator for shrinkage parameters from the estimator in our package (Section 2.3). However, the bigRR package was already removed from the CRAN. The R package “genridge” is designed for the ridge estimator with a single shrinkage parameter, and hence it actually cannot calculate the generalized ridge estimator of multiple shrinkage parameters. Similarly, the R package “penalize” (https://cran.r-project.org/package=penalized, accessed on 5 February 2024) has the “optL2(.)” function to perform ridge regression with a single shrinkage, which does not allow multiple shrinkage parameters. In the same manner, an R package “lmridge” [49] was developed to compute the ridge estimator, allowing 25 different ways to estimate the single shrinkage parameter. In summary, these existing R packages are insufficient for using the generalized ridge estimator with multiple shrinkage parameters.

While our simulation study demonstrated the high estimation accuracy of the generalized ridge estimator, the computation time was not fast. For a single simulation run, computation took about 3 s ($p=50$), 6 s ($p=100$), 12 s ($p=150$), and 24 s ($p=200$). This means that the time grows steeply as the dimensions increase. In order to analyze very-high-dimensional regressors, such as SNPs, the package needs to be improved for computational efficiency.

The generalized ridge estimator considered in this article may be extended to logistic regression, Poisson regression, and Cox regression. Such extensions have not been explored yet. The extensions may require some efforts for formulating an estimation algorithm and efficient cross-validation criterion [50], performing extensive simulation studies, and developing a software package to fully justify the advantages and usefulness over the usual ridge estimator that is widespread via the “glmnet” package.