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

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 (Hoerl and Kennard 1970; Montgomery et al. 2021) is an effective alternative for the high-dimensional ($p>n$) data, and is widely employed in such data encountered in genetics (Arashi et al. 2021; Vishwakarma et al. 2021), medical studies (Friedrich et al. 2023), electromagnetic signals (Gao et al. 2023), grain yields (Hernandez et al., 2015), and others.

Hoerl and Kennard (1970) originally proposed the ridge estimator to reduce the problem in multicollinearity. Later, the ridge estimator is naturally adopted to the high-dimensional ($p>n$) problem (Golub et al., 1979; Hastie et al., 2009) as a way to avoid overfitting phenomenon.

Ridge regression is a shrinkage estimator that shrinks all regression coefficients uniformly toward zero (Hoerl and Kennard 1970; Hastie et al., 2009). 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. (2011) applied the ridge estimator to the high-dimensional SNP data and performed significance testing for selecting an informative subset of SNPs. Similar applications of ridge regression to high-dimensional data include Whittaker et al. (2000), Cule and De Lorio (2013), and Yang and Emura (2017).

Unlike the ordinary ridge regression that shrinks all regression coefficients uniformly, the generalized ridge regression allows different degrees of shrinkage under multiple shrinkage parameters. Allen (1974) and Loesgen (1990) demonstrated that the multiple shrinkage parameters in the generalized ridge estimator arise naturally by utilizing prior information about regression coefficients; see also an extensive review paper of Yüzbaşı et al. (2020) for generalized ridge regression methods.

However, multiple shrinkage parameters are difficult to be estimated for high-dimensional cases. To overcome this difficulty, Yang and Emura (2017) 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 is to assign two different weights (1 or 0.5) for shrinkage parameters by preliminary tests (Saleh and Kibria 1993; 2019; Norouzirad and Arashi 2019; Shih et al., 2021;2023; Taketomi et al. 2023). While this approach is shown to be promising due to its sound statistical performance under sparse and high-dimensional models, none of the 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) that implements both the ridge estimator (Hoerl and Kennard 1970) and generalized ridge estimator (Yang and Emura 2017). The package equips with the generalized cross-validation for 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 the sparse and high-dimensional models, and the normal and skew-normal distributions for error terms. We conclude that the generalized ridge estimator is superior to the benchmark ridge estimator based on “glmnet”, and hence, it can be the most recommended estimator under sparse and high-dimensional models. We illustrate the package via the intracerebral hemorrhage data.

The remainder of the 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 (1970) and the generalized ridge regression proposed by Yang and Emura (2017).

2.1. Linear regression

Consider the linear regression model $\mathit{y}=X\mathit{\beta }+\mathit{\epsilon }$, where

$$\mathit{y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right],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],\mathit{\beta }=\left[\begin{array}{c}{\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right],\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 the 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]=0$. In some case, we assume $\mathrm{C}\mathrm{o}\mathrm{v}\left[\mathit{\epsilon }\right]={\sigma }^2{I}_n$, where $I_n$ is the identify matrix of dimension $n$. The assumption of the covariance is necessary to obtain the standard error (SE) of regression coefficients and testing 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 (Hastie et al. 2009).

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

Ridge estimator is an alternative to the LSE, which can be computed even when $p>n$. Hoerl and Kennard (1970) 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 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)$ shrunk toward zero. Theorem 4.3 of Hoerl and Kennard (1970) shows that there exists $\lambda >0$ such that the ridge estimator yields strictly smaller mean squared error (MSE) than that of the LSE.

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

$$V\left(\lambda \right)={\frac{1}{n}\left|\right|\{I_n-A(\lambda \left)\right\}\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 the “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 (Golub et al. 1979) guarantees the use of the above ridge estimator under the $p>n$ setup. Noted that the GCV criterion is different from the cross-validation criterion that is available in the widely used R function “cv.glmnet(.)” in the R package “glmnet”. There are many other available criteria for choosing $\lambda$ (Cule et al. 2013; Wong and Chiu 2015; Kibria and Banik 2016; Assaf et al. 2019; Michimae and Emura 2020), most of which are not applicable for the $p>n$ setup. Therefore, we will adopt the GCV criterion for the following discussions, which works well for both $p<n$ and $p>n$ setups.

2.3. Generalized ridge regression

Yang and Emura (2017) suggests relaxing the uniform shrinkage to the non-uniform shrinkage to yield the generalized ridge regression. For this purpose, Yang and Emura (2017) 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}\left(\Delta \right){\}}^{-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}\left(\Delta \right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{\widehat{w}}_1\right(\Delta ),...,{\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{)}^2/\left(p-1\right)\right\}}^{1/2}$ for $j=1,...,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 row of $X$. Note that ${\widehat{\mathit{\beta }}}^0$ is called “compound covariate estimator” (Chen and Emura 2017).

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

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

(2)

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

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

Computation of $(\hat{\lambda },\hat{\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 get $\hat{\lambda }\left(\Delta \right)$. Under the sparse model ($\mathit{\beta }\approx 0$), the histogram of ${\hat{\beta }}_j^0/\mathrm{S}\mathrm{D}\left({\hat{\mathit{\beta }}}^0\right)$, $j=1,...,p$, is well-approximated by $N\left(\mathrm{0,1}\right)$. This implies that $\left|{\hat{\beta }}_j^0\right|/\mathrm{S}\mathrm{D}\left({\hat{\mathit{\beta }}}^0\right)$ falls in the range $\left[\mathrm{0,3}\right]$ with nearly 99.73%, and hence, a search range $\Delta \in \left[\mathrm{0,3}\right]$ suffices. Since $V\left(\hat{\lambda }\right(\Delta ),\Delta )$ is discontinuous in $\Delta$, a grid search was suggested on the grid $D=\{0,3/100,...,300/100\}$.

Finally, the generalized ridge estimator is defined as

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

Also, the error variance can be estimated by

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

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

2.4. Significance test

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

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

for $j=1,...,p$. One can perform significance testing, allowing one to access P-values of all the $p$ regressors (Cule et al. 2011; Cule and De Lorio 2013; Yang and Emura 2017). Since the significance tests are similar between the ridge and generalized ridge estimators, we will explain the significance tests based on the generalized ridge estimator below.

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

$$\mathrm{C}\mathrm{o}\mathrm{v}\left\{\hat{\mathit{\beta }}\right(\hat{\lambda },\hat{\Delta }\left)\right\}={\hat{\sigma }}^2\{X^{\mathrm{T}}X+\hat{\lambda }\hat{W}(\hat{\Delta }\left){\}}^{-1}{X}^{\mathrm{T}}X\right\{X^{\mathrm{T}}X+\hat{\lambda }\hat{W}\left(\hat{\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 the main R function “g.ridge(.)“, and the other function “X.mat(.)”. Appendix A explains another function “GCV(.)” that may not be directly used, but is useful for internal computing.

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 (2017, p.6093), we implemented a convenient R function that generates $X$ having three independent blocks:

where the first block consists of $q\ge 0$ correlated regressors (correlation = 0.5), and 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},x_{ik})=\left\{\begin{array}{cc}0.5& \mathrm{if}\hspace{1em}j,k\in \{1,...,q\},\hspace{1em}\hspace{1em}\hspace{1em}\\ 0.5& \mathrm{if}\hspace{1em}j,k\in \{q+1,...,q+r\},\\ 0& \mathrm{otherwise}.\end{array}\right.$$

The design matrix mimics the “gene pathway” for two blocks of correlated gene expressions often used in simulation studies (Binder et al. 2009; Emura et al. 2012; 2016; 2023). 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},...,z_{ip}$ and $u_i,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 the ridge estimator $\hat{\mathit{\beta }}\left(\hat{\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 of the command includes $\hat{\lambda }\in [0,500]$ and $\hat{\mathit{\beta }}\left(\hat{\lambda }\right)$. Similarly, the command “g.ridge(X, Y, method = "YE", kmax = 500)” can calculate $\hat{\mathit{\beta }}(\hat{\lambda },\hat{\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 $\hat{\lambda }\in [0,500]$ (the plot of $V(\lambda ,\hat{\Delta })$ for the generalized ridge).

The R function “g.ridge(.)” can also compute the SE of $\hat{\mathit{\beta }}$, Z-value and P-value for significance tests (Section 2.4), and the estimate of the error variance ${\hat{\sigma }}^2$ (Section 2.3). As in the typical practice of ridge regression, we used the centered responses “Y-mean(Y)” rather than “Y” (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 $\hat{\lambda }=31.66314$, $\hat{\mathit{\beta }}\left(\hat{\lambda }\right)=\left(0.581485036,\dots 0.083260387\right)$, SE, Z-value, P-value, and ${\hat{\sigma }}^2=1.778618$. The output also displays the GCV function $V\left(\lambda \right)$ against $\lambda \in [0,200]$, showing its minimum at $\hat{\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.

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

3.3. Technical remarks on centering and standardization

We assume that $X$ is standardized and does not include the intercept term (Section 2.1). If $X$ is generated by “X.mat”, it is already standardized, and hence there is no need to do the 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 the intercept term (a column of ones) since one can always use the reduced model $\mathit{y}-({\sum }_{i=1}^n{y}_i/n)1=X\mathit{\beta }+\mathit{\epsilon }$ for the intercept model $\mathit{y}={\beta }_{0}1+X\mathit{\beta }+\mathit{\epsilon }$ (${\beta }_0$ is estimated by ${\sum }_{i=1}^n{y}_i/n$), where $1$ is a vector of ones. This means that the usual unbiased estimator is applied to the intercept.

4. Simulations

We conducted simulation studies to test our R package “g.ridge”. The goals of the simulations are to show: (a) the generalized ridge estimators in our package exhibits superior performance over the ridge estimator in the package “glmnet(.)”, and (b) the sound performance under the normally and non-normally distributed errors (the skew-normal distribution will be considered). 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,100,150,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 the sparse model:

$$\mathit{\beta }=(\stackrel{q}{\overbrace{b/q,...,b/q}},\stackrel{r}{\overbrace{d/r,...,d/r}},\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 (Azzalini 2014); both distributions have mean zero and standard deviation one, and the skew-normal distribution has the slant parameter ten (alpha=10 in the R function “rsn(.)”). Figure 3 shows the remarkable difference of the two distributions. The skew-normal distribution was not previously examined in the simulation setting of Yang and Emura (2017).

For a given dataset $(X,\mathit{y})$, we computed $\hat{\mathit{\beta }}$ that can be (i) the ridge estimator by “g.ridge(.)”, (ii) the generalized ridge estimator by “g.ridge(.)”, or (iii) the ridge estimator by “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(\hat{\mathit{\beta }}\right)=\mathrm{E}\left[{‖\hat{\mathit{\beta }}-\mathit{\beta }‖}^2\right]=\mathrm{E}\left[{\sum }_{j=1}^p{\left({\hat{\beta }}_j-{\beta }_j\right)}^2\right],$$

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

4.2. Simulation results

Table 1 compares the MSE 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 it 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 R proposed package is 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 study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of Saiseikai Kumamoto Hospital on March 29, 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.

The outcome variables ($\mathit{y}$) are the changes of the blood volume (in mL) from the initial value and to the follow-up value, which were measured by the CT scan. Excluding patients with less than 10 mL blood volume at the initial CT scan and other 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^-3/μL; C-reactive protein, mg/dL). The responses and regressors are centered and standardized before fitting a linear model as explained in Section 2.1 and Section 3.3.

Figure 3 displays scatter plots for the centered responses $\mathit{y}-({\sum }_{i=1}^n{y}_i/n)1$ against the predictors $X\hat{\mathit{\beta }}$ based on the ridge estimator and generalized ridge estimator. We observe that the predictors reasonably captured the variability of the changes in the blood volume (response variables). However, the figure also shows that the changes in the blood volumes were not fully explained by the predictors since some residual errors remain. Figure 4 depicts the residuals 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 proposed ridge estimators was verified for the asymmetric errors (Section 4).

Table 2 shows the fitted results for estimated regression coefficients sorted by P-values (only with P < 0.05). The fitted results are quite similar between the ridge and generalized ridge estimators. The most significant regressor is the lactate dehydrogenase. This variable was previously reported as an important predictor of hematoma expansion (Chu et al. 2019). An interesting difference is found in the number of significant regressors: 6 regressors for the generalized ridge estimator, and 5 regressors for the ridge estimator. That is, the C-reactive protein as a regressor selected solely by the generalized 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 the severe intracerebral hemorrhage (Wang et al. 2021). 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.

Finally, we compared the performance of the ridge and generalized ridge estimators by means of the residual mean squared error (RMSE) defined as

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={‖\mathit{y}-\left(\frac{{\sum }_{i=1}^n{y}_i}{n}\right)1-X\hat{\mathit{\beta }}‖}^2.$$

The RMSE were 2.288 (the ridge estimator) and 2.284 (the generalize ridge estimator). Therefore, the predictor constructed from the generalized ridge estimator had slightly better performance over the one from the ridge estimator.

Therefore, we have demonstrated that the ridge and generalize ridge estimators in the proposed R package provide statistically and biologically sound conclusions on the real data analysis.

6. Conclusion

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 the abundance of sparse and high-dimensional data (Kim et al. 2007; Zhang et al. 2021; Vishwakarma et al. 2021; Bhattacharjee 2022; Bhatnagar et al. 2023; Emura et al. 2023) and asymmetrically distributed data (Abe et al. 2012; Huynh et al. 2021; Yoshiba et al. 2023; Jimichi et al. 2023), the proposed package may provide a reliable statistical tool for statisticians and data scientists.

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. While the extensions might not be technically difficult, well-designed simulation studies and implementations in some software packages will be necessary to fully justify the advantage and usefulness over the usual ridge estimator that is widespread via the “glmnet” package.