Penalty Strategies in Semiparametric Regression Models

1. Introduction

The Partially Linear Regression Model (PLRM), introduced by [16], serves as a versatile tool for utilizing the advantage of both parametric linear models and nonparametric regression models. It allows for modeling complex relationships by incorporating both linear and nonlinear components, thereby offering greater flexibility and interpretability across various scientific disciplines, including social sciences, biology, and economics. Formally, the PLRM is defined as:

$$y_i={\mathbf{x}}_i^{\text{'}}\mathsf{\beta }+f\left(t_i\right)+{\epsilon }_i, i=1,\dots ,n$$

(1.1)

where $y_i\text{'}s$ are the observations of the response variable, ${\mathbf{x}}_i=(x_{i1},\dots ,x_{ik})$ and $t_i\text{'}s$ are the values of the explanatory variable, $\mathsf{\beta }={({\beta }_1,\dots ,{\beta }_k)}^{\mathit{\text{'}}}$ is an unknown k-dimensional parameter vector to be estimated, $f\left(.\right)$ is an unknown smooth function, and ${\epsilon }_i\text{'}s$ are assumed to be uncorrelated random variables with mean zero and a common variance ${\sigma }^2$ and independent of the explanatory variables. Note also that the observations of $y_i$ depend linearly on the entries of ${\mathbf{x}}_i$ and nonlinearly on the values of a univariate variable $t_i$. The PLRM can be re-expressed in matrix and vector form

$$\mathbf{y}=\mathbf{X}\mathsf{\beta }+\mathbf{f}+\mathit{\epsilon }$$

(1.2)

where $\mathbf{y}={\left(y_1,\dots ,y_n\right)}^{\text{'}}$, $\mathbf{X}={[{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n]}^{\text{'}}$ is an ($n\times k$) design matrix with ${{\mathbf{x}}_i=(x_{i1},\dots ,x_{ik})}^{\text{'}}$ denoting the i.th$k-$dimensional row vector of the design matrix X, and $\mathbf{f}={\left(f\left(t_1\right),\dots ,f\left(t_n\right)\right)}^{\mathit{\text{'}}}$ and $\mathsf{\epsilon }=\left({\epsilon }_1,.{\epsilon }_2,\dots ,{\epsilon }_n\right)$ are random error vectors with $E\left(\mathit{\epsilon }\right)=0$ and $Var\left(\mathit{\epsilon }\right)={\mathsf{\sigma }}^2{\mathbf{I}}_n$. For more discussion on Model (1.1), see [20,35,39], among others. Since it was initially presented by [16], this model has been very popular and is frequently used in the social, biological, and economic sciences. The PLRM generalizes both the parametric linear regression model and nonparametric regression model. When the $\mathbf{f}=0$ model (1.1) reduces to the linear regression model$y_i={\mathbf{x}}_i\mathsf{\beta }+{\epsilon }_i$, when the $\mathsf{\beta }=0$ model (1.1) becomes a nonparametric regression model,

$$y_i=f\left(t_i\right)+{\epsilon }_i$$

(1.3)

with a univariate covariate. In addition, because PLRMs contain both linear and nonlinear components, they are more flexible than linear models are. In this work, we are interested in estimating the parameter vector $\mathsf{\beta }$ and the function $\mathbf{f}$.

Estimating the parameters and the smooth function in PLRM has been a focal point of research, with various methodologies proposed to enhance estimation accuracy and model interpretability. Early contributions by [22] and [39] utilized kernel smoothing techniques for estimating the nonparametric component, laying the groundwork for subsequent advancements.

The advent of regularization methods introduced by [41] and [54] revolutionized parameter estimation by addressing issues such as multicollinearity and high-dimensionality through penalty functions. In particular, [4,44] explored penalized least squares in semiparametric models, emphasizing the integration of penalty functions with kernel smoothing to enhance estimator performance. Further, [10,11,47] and [48].

Adaptive methods, such as the Adaptive Lasso (aLasso) introduced by [52], have been incorporated into PLRM frameworks to improve variable selection consistency. Moreover, SCAD proposed by [18] and MCP introduced by [51] penalties have been employed to mitigate the bias inherent in $L_1$ penalties, providing more accurate parameter estimates in complex models. Recent studies by [46] and [12] highlight the asymptotic properties of penalized estimators in semiparametric regression, emphasizing the significance of penalty strategies in PLRM. Challenges in balancing bias and variance persist, particularly in sparsity contexts with multicollinearity.

In regression problems, multicollinearity emerges when variables are highly correlated with each other and the dependent variable. This situation inflates confidence intervals, weakens parameter estimates, and leads to less reliable predictions. It also raises sampling variance, reducing accuracy for both inference and prediction. To counter these challenges in PLRMs, we integrate shrinkage estimation with penalty functions to bolster robustness and manage both multicollinearity and variable selection:

• We tackle variable selection and multicollinearity in Model (1.1), which features strongly correlated covariates, using six penalty functions: Ridge, Lasso, aLasso, SCAD, Elastic Net, and MCP.
• We then review and compare these methods, highlighting their strengths, limitations, and uses within PLRMs.
• To further enhance parameter estimation, we incorporate shrinkage estimation techniques—both standard and positive shrinkage estimators—into the penalty-based framework.
• Finally, kernel smoothing via penalized least squares is employed to estimate both the parametric (linear) and nonparametric (smooth) components of the model.

The remainder of this paper is structured as follows. Section 1.2 summarizes the background on penalty functions, partially linear models, and their integration. Section 2 describes the kernel smoothing method. Section 3 introduces the kernel-type ridge estimator. In Section 4, we discuss penalty functions and shrinkage estimators based on kernel smoothing. Section 5 presents the evaluation metrics for both parametric and nonparametric components. Section 6 outlines the estimators’ asymptotic properties. Section 7 details the Monte Carlo simulations. Section 8 compares the performance of our proposed estimators on real data. Finally, Section 9 provides concluding remarks and includes supplementary technical materials.

Background

In this section, we present a relevant literature review of partially linear regression models (PLRMs) and show why penalty functions and shrinkage estimators are essential for handling challenges like multicollinearity, high-dimensional data, and complex nonlinearities. We also explain how kernel smoothing complements these approaches by capturing intricate patterns in the data based on the previous studies. Altogether, these methods offer a strong foundation for building reliable and flexible statistical models.

Partially Linear Regression Models (PLRMs) are a cornerstone of semiparametric statistics, offering a powerful way to model relationships that aren’t fully captured by purely parametric or nonparametric approaches. They blend the interpretability of linear models with the flexibility of nonparametric methods, making them ideal for data with both linear and nonlinear components. Their ability to isolate certain covariate effects while accommodating nonlinearities has led to widespread use in fields such as economics, biostatistics, and environmental science [25]. Recent developments have expanded their scope to high-dimensional, longitudinal, and survival data ([38]; [37]).

Regularization methods, particularly those using penalty functions, have revolutionized regression analysis by tackling issues like multicollinearity, high dimensionality, and model selection. By constraining coefficient magnitudes, these methods shrink coefficients toward zero and improve model stability. Different penalty functions yield diverse estimators suited to various data characteristics and research goals.

Ridge regression, pioneered by [21], is well-established, yet its application in high-dimensional contexts and its ties to Bayesian methods continue to be studied [45]. [41] Lasso has also evolved substantially, with research on its consistency in variable selection, high-dimensional settings, and extensions to generalized linear models [14]. Ahmed and colleagues have contributed to Lasso’s development in PLRMs using Stein-type shrinkage estimators, especially under multicollinearity [48].

Building on [53] adaptive Lasso (aLasso), researchers have explored data-driven weights for improving variable selection and estimation. For SCAD, introduced by [18], recent work has focused on high-dimensional data, quantile regression, and robust estimation [42]. The researchers applied SCAD along with the other penalty functions and with shrinkage estimators in PLRMs, showing its effectiveness in handling outliers [10]. MCP, initially studied by [51], has inspired efficient algorithms and further theoretical exploration [13]. Elastic Net, proposed by [53], has also gained traction in areas like genomics, imaging, and finance, where it handles correlated predictors well [54]. its benefits are showed for balancing bias and variance in PLRMs [3].

Beyond these penalty-based approaches, shrinkage estimation techniques—particularly Stein-type estimators—offer another route. [5] investigate partially linear models by splitting the coefficient vector into main and nuisance effects and comparing shrinkage, pretest, and absolute penalty estimators. Their results show that a Stein-type shrinkage estimator uniformly outperforms the conventional least squares approach—especially when the nuisance parameter space is large—while the pretest estimator only offers benefits in a limited region. These methods shrink coefficients toward a chosen target, thereby reducing variance and mitigating noise ([49]; [23]; [28]). Further, Ahmed’s work demonstrates that these techniques can improve estimation accuracy and interpretability in partially linear regression models ([4]; [30]; [2]). Meanwhile, kernel smoothing remains central to nonparametric estimation for capturing nonlinear relationships, with advances such as adaptive bandwidth selection [24] and its integration with penalty and shrinkage methods [17]. [50] showed via simulation that applying Stein-type shrinkage to the parametric components of [26] semi-nonparametric regression model enhances accuracy. [9] proposed a two-stage approach using LASSO/Elastic-Net for variable selection followed by post-selection shrinkage estimation to improve predictions in high-dimensional sparse linear regression models, while [1] introduced shrinkage-type estimators for reliability and stress–strength parameters from progressively censored data that outperform MLEs in simulations and industrial applications. Additional insights on shrinkage estimators are provided in [6,7,29,31,32,33,36].

Combining penalty functions with kernel smoothing has become a potent strategy for complex data in PLRMs. This approach refines the parametric part through penalty-driven selection and shrinkage, while kernel smoothing captures nonlinearities. Recent studies highlight theoretical and practical benefits, including improved estimation accuracy and consistency in variable selection under multicollinearity or high dimensionality [8]. The advanced penalized kernel smoothing in PLRMs using shrinkage estimators are shown by [2], showing strong performance in simulations and real data. Their work illustrates the advantages of using penalties like aLasso and shrinkage estimation alongside kernel smoothing, especially with high-dimensional covariates and intricate nonlinearities.

2. Estimation Based on Kernel Smoothing

Let us first consider the nonparametric estimation of the unknown regression function $f\left(t\right)$ in Model (1.1). For simplicity, we assume that $\mathit{\beta }$ is known from equation (1.1). The relationship between $y_i- {\mathit{X}}_{\mathit{i}}\mathit{\beta }$ and $t_i$ in this instance may be represented as

$$\left(y_i-{\mathbf{X}}_{\mathbf{i}}\mathsf{\beta }\right)=f\left(t_i\right)+{\epsilon }_{i, } i=1, \dots , n$$

(2.1)

The nonparametric part of the semiparametric model is equivalent to equation (2.1). This results in the Nadaraya‒Watson estimator, commonly known as the kernel estimator, which was proposed by [27] and [43], as mentioned in [39]:

$$\widehat{f}\left(h,t\right)={\sum }_{i=1}^n{w}_i(h, t_i)(y_i-{\mathbf{X}}_{\mathbf{i}}\mathsf{\beta })={\mathbf{W}}_h(\mathbf{y}-\mathbf{X}\mathsf{\beta })$$

(2.2)

where ${\mathbf{W}}_h$ is a kernel smoothing matrix with $j-th$ entries${ w}_{hi}$ and where $h$ is a smoothing parameter (or bandwidth), as shown by:

$$w_i\left(h,t_i\right)=K\left({\displaystyle \frac{t-t_i}{h}}\right)/{\displaystyle \sum _{i=1}^{n}K\left(\frac{t-t_i}{h}\right)}=K\left(u_i\right)/{\displaystyle \sum _{i}K}\left(u_i\right)$$

(2.3)

To estimate $f\left(t\right)$, kernel smoothing (regression) uses appropriate weights $w_i(h,t)$, as shown in (2.1). The kernel function $K\left(u\right)$ with a smoothing parameter $h$, which defines the size of the neighborhood around $t$, directs the weights assigned to the observations $t_i$ ([40]). The kernel or weight function $K\left(u\right)$ in equation (2.3) has the properties of $\int K\left(u\right)du=1$ and $K\left(u\right)=K(-u)$. The kernel function is chosen to give the closest observations the most weight and the farthest from $t$ the least weight.

The following partial residuals may be obtained via the matrix and vector forms of model (1.2) in matrix form:

$$\mathsf{\epsilon }=\mathbf{y}-\mathbf{X}\mathsf{\beta }-\widehat{\mathbf{f}}=\left(\mathbf{I}-{\mathbf{W}}_h\right)\left(\mathbf{y}-\mathbf{X}\mathsf{\beta }\right)=\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }$$

(2.4)

where $\tilde{\mathbf{X}}=\left(\mathbf{I}-{\mathbf{W}}_h\right)\mathbf{X}$ and $\tilde{\mathbf{y}}=\left(\mathbf{I}-{\mathbf{W}}_h\right)\mathbf{y}.$ Hence, a transformed set of data based on kernel residuals is what we obtain as a result. The following weighted least squares (WLS) criterion is obtained by considering these partial residuals for the vector $\mathsf{\beta }$:

$$\mathbf{W}\mathbf{L}\mathbf{S}\left(\mathsf{\beta };h\right){=\left(\left(\mathbf{I}-{\mathbf{W}}_h\right)\left(\mathbf{y}-\mathbf{X}\mathsf{\beta }\right)\right)}^{\mathbf{T}}(\mathbf{I}-{\mathbf{W}}_h)(\mathbf{y}-\mathbf{X}\mathsf{\beta }{\left)\right)=(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta })}^{\mathbf{\text{'}}}(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta })$$

(2.5)

The result for the criterion $\mathbf{W}\mathbf{L}\mathbf{S}\left(\mathsf{\beta };h\right)$ provided in equation (2.4) is easily demonstrated to be similar to the solution of ordinary least squares (OLS)

$${{\widehat{\mathsf{\beta }}}_{KS}=\left({\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}\mathbf{X}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}\tilde{\mathbf{y}}$$

(2.6)

Additionally, revising the steps for $f\left(t\right)$ reduces equation (2.1) to the form:

$$\widehat{f}\left(h,t\right)=\sum _{i=1}^n{\displaystyle \frac{K\left(\frac{t-t_i}{h}\right)}{\sum _{i=1}^{n}K \left(\frac{t-t_i}{h}\right)}}(y_i-x_i\widehat{\beta })$$

(2.7)

Another way to express equation (2.7) in a matrix is as follows:

$${\widehat{\mathbf{f}}}_{KS}={\mathbf{W}}_h(\mathbf{y}-\mathit{X}{\widehat{\mathsf{\beta }}}_{KS}$$

(2.8)

Thus, our estimation for ${\mathbf{y}}_{KS}$ becomes:

$${\widehat{\mathbf{y}}}_{KS}= \mathit{X}{\widehat{\mathsf{\beta }}}_{KS}+{\widehat{\mathbf{f}}}_{KS}=\mathbf{X}({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\tilde{\mathbf{X}})}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}\tilde{\mathbf{y}}+{\mathbf{W}}_h\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{KS}\right)={\mathbf{H}}_{KS}\left(h\right)\mathbf{y}$$

for the hat (projection) matrix

$${\mathbf{H}}_{KS}\left(h\right)={\mathbf{W}}_h+\tilde{\mathbf{X}}{\left({\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}\mathbf{X}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}(\mathbf{I}-{\mathbf{W}}_h)$$

(2.9)

3. Kernel Type Ridge Estimator in Semiparametric Model

For the sake of this paper, we limit our attention to kernel smoothing estimators of the vector parameter $\mathsf{\beta }$ and the unidentified smooth function $f\left(.\right)$ in a semiparametric model. The corresponding estimators $\mathsf{\beta }$ and $\mathbf{f}$ are based on model (1.2) for a given bandwidth parameter $h$, and by multiplying $\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)$ on both sides of model (1.2), we obtain the following model:

$$\tilde{\mathbf{y}}=\tilde{\mathbf{X}}\mathsf{\beta }+\tilde{\mathit{\epsilon }}$$

(3.1)

where $\tilde{\mathbf{y}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathbf{y}$, $\tilde{\mathsf{\epsilon }}=\tilde{\mathbf{f}}+{\mathsf{\epsilon }}^{\mathrm{*}},$ $\tilde{\mathbf{f}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathbf{f}$ and ${\mathsf{\epsilon }}^{\mathrm{*}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathsf{\epsilon }$.

These considerations show that Model (3.1) is the optimal solution for obtaining the vector $\mathsf{\beta }$ estimator corresponding to the parametric part of the semiparametric model (1.2). For the ridge regression problem, this model provides the following penalized least squares (PLS) criterion:

$$PL{S}_{RK}=\underset{\mathsf{\beta }}{\mathit{arg\; min}}{\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)}^{\mathit{\text{'}}}\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)+\lambda {\mathsf{\beta }}^{\mathit{\text{'}}}\mathsf{\beta }$$

(3.2)

where $\lambda$ denotes a positive shrinkage parameter that regulates the severity of the penalty. The solution of the minimization problem (3.2) can be written as follows as ridge-type kernel smoothing estimator (see [46]) similar to [21] but modified based on partial residuals obtained using smoothing matrix ${\mathbf{W}}_h:$

$${{\widehat{\mathsf{\beta }}}_{RK}\left(\lambda \right)=\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}\tilde{\mathbf{X}}+\lambda {\mathbf{I}}_k\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}\tilde{\mathbf{y}}$$

(3.3a)

The ridge-type kernel smoothing estimate simplifies to an ordinary least squares estimation problem when $\lambda =0$ on the basis of the local residuals specified in equations (2.6) and (2.8). Additionally, to estimate the unknown function $\mathbf{f}$, we imitate equation (2.8) and define it.

$${\widehat{\mathbf{f}}}_{RK}={\mathbf{W}}_h\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{RK}\left(\lambda \right)\right)$$

(3.3b)

Thus, in the semiparametric model (1.2), the estimator (3.3b) is therefore defined as the ridge type kernel estimator of the unknown function $\mathbf{f}$

4. Penalty Functions and Shrinkage Estimators

Many penalty functions for linear and generalized regression models have been studied ([2,3]). Here, we consider the MCP, Lasso, SCAD, aLasso, and ElasticNet functions and shrinkage techniques additionally. Note that ElasticNet is a regularized regression technique that linearly combines the${ L}_1$ and $L_2$ penalties of the Lasso and ridge regression methods, respectively.

In this paper, we introduce the kernel smoothing estimators based on several penalties for the semiparametric regression model's components. For a given penalty function and tuning parameter $\lambda$, the general form of the penalty estimators is defined by the following penalized least squares $\left(PLS\right)$ criterion:

$$LS=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2\right\}+P_{\lambda }\left(\beta \right)={\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)}^{\mathit{\text{'}}}\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)+P_{\lambda }\left(\beta \right)$$

(4.1)

where $P_{\lambda }\left(\mathsf{\beta }\right)$ denotes the penalty function based on the shrinkage (or tuning) parameter$\lambda .$ This parameter needs to be selected with any selection criteria.

Note that [19] bridge estimator is the vector $\widehat{\mathsf{\beta }}$ that minimizes the Lasso and Ridge penalties in (4.1). However, SCAD, ElasticNet, aLasso, and MCP have different penalties, which are shortly introduced in the remainder of this study. Note that ${\left|{\beta }_j\right|}^q$ fulfills the $L_q$ norm of the regression coefficients${\beta }_j$ because of the penalty functions used.

Hence, for various values of the shrinkage parameter $\lambda$, different penalty estimators that correspond to the parametric and nonparametric parts of the semiparametric model can be defined.

4.1. Estimation Procedure for the Parametric Component

According to (4.1), it is possible to obtain ridge estimates for the parametric component for $q=2$ by minimizing the following penalized residual sum of squares.

$${\widehat{\mathsf{\beta }}}_{RK}=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2+\lambda \sum _{j=1}^k{\left({\beta }_j\right)}^2\right\}$$

(4.2)

where ${\tilde{\mathbf{x}}}_i^{\text{'}}$ is the $i.$th row of the matrix$\tilde{\mathbf{X}}$, and${\tilde{y}}_i$ is the $i^{th}$ observation of $\tilde{\mathbf{y}}$. Note that the regularization estimate in solution (4.2) matches that in (3.3a). Alternative methods to ridge are described below.

Lasso: [41] proposed Lasso, a regularization method that uses the${ L}_1$ penalty for estimation and variable selection. The modified kernel smoothing estimators based on the Lasso penalty can be defined as

$${\widehat{\mathsf{\beta }}}_{LK}=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2+\lambda \sum _{j=1}^k\left|{\beta }_j\right|\right\}$$

(4.3)

The absolute penalty term prevents Lasso from having an analytical solution, despite the apparent simplicity of equation (4.3).

aLasso: [52] proposed using adaptive weights on $L_1$ penalties on the regression coefficients to modify the Lasso penalty. The term "aLasso" refers to this weighted Lasso with oracle properties. The kernel smoothing estimator ${\widehat{\mathsf{\beta }}}_{aLK}$ using aLasso penalty is defined as follows:

$${\widehat{\mathsf{\beta }}}_{aLK}=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2+\lambda \sum _{j=1}^k{\widehat{w}}_j\left|{\beta }_j\right|\right\}$$

(4.4)

where $w$ is a weight function given by

$${\widehat{w}}_j={\displaystyle \frac{1}{{\left|{\widehat{\beta }}^{*}\right|}^q}},q>0$$

Notably, ${\widehat{\beta }}^{*}$ is an appropriate estimator of $\mathsf{\beta }$. In addition, an ordinary least squares (OLS) estimate can be used as a reference value. After the OLS estimate is obtained, $q>0$ is selected, and the weights are computed to obtain the adaptive-Lasso estimates in (4.4).

SCAD: The Lasso method's penalty term increases linearly with the regression coefficient's size, leading to highly biased estimates for coefficients with large values. To address this issue, [18] proposed a SCAD penalty that can be obtained by replacing $\left|{\beta }_j\right|$ in (4.3) with $P_{\alpha ,\lambda }\left|{\beta }_j\right|.$ A modified kernel smoothing estimator ${\widehat{\mathsf{\beta }}}_{SK}$ based on the SCAD penalty can be expressed as

$${\widehat{\mathsf{\beta }}}_{SK}=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2+\sum _{j=1}^k{P}_{\alpha ,\lambda }\left|{\beta }_j\right|\right\}$$

(4.5)

where $P_{\alpha ,\lambda }(.)$ is the SCAD penalty defined by

$$P_{\alpha ,\lambda }=\lambda \left\{\mathbf{I}\left(\left|\beta \right|\right)\le \lambda +{\displaystyle \frac{{\left(\alpha \lambda -\left|\beta \right|\right)}_{+}}{\left(\alpha -1\right)\lambda }}\mathbf{I}(\left|\beta \right|>\lambda )\right\}, for \beta \ge 0$$

(4.5a)

The penalty parameters in this case are $\lambda >0$ and $\alpha >2$, where ${\left(t\right)}_{+}=\mathrm{m}\mathrm{a}\mathrm{x}(t,0)$ and $\mathbf{I}(.)$ is the indicator function. When $\alpha =\infty$, $,\left(4.4a\right)$ is equivalent to the $L_1$ penalty.

ElasticNet:ElasticNet is a penalized least squares regression method proposed by [53] that is widely used for regularization and automatic variable selection to choose groups of correlated variables. It combines the $L_1$ and the $L_2$ penalty term, where the $L_1$ penalty term ensures sparsity, whereas the $L_2$ penalty term selects correlated variable groups. Accordingly, the modified kernel-smoothing estimator ${\widehat{\mathsf{\beta }}}_{ENK}$ using an ElasticNet penalty is the solution to the following minimization problem:

$${\widehat{\mathsf{\beta }}}_{ENK}=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2+{\lambda }_1\sum _{j=1}^k{\left|{\beta }_j\right|}^2+{\lambda }_2\sum _{j=1}^k\left|{\beta }_j\right|\right\}$$

(4.6)

In the expression above, the positive regularization parameters are ${\lambda }_1$ and ${\lambda }_2$. Note that (4.6) ensures that the estimates correspond to the parametric part of the semiparametric regression model (1.2).

MCP: [51] introduced MCP as a different method to obtain less biased estimates of nonzero regression coefficients in a sparse model. For given regularization parameters $\lambda >0$ and $\alpha >0$, the kernel smoothing estimator ${\widehat{\mathsf{\beta }}}_{MCK}$ based on the MCP penalty can be defined as

$${\widehat{\mathsf{\beta }}}_{MCK}=\underset{\mathsf{\beta }}{arg min}\left\{\sum _{i=1}^n{({\tilde{y}}_i-{\tilde{\mathbf{x}}}_i^{\text{'}}\mathsf{\beta })}^2+\sum _{j=1}^k{P}_{\alpha , \lambda }\left(\left|{\beta }_j\right|\right)\right\}$$

(4.7)

where $P_{\alpha , \lambda }(.)$ is the MCP penalty expressed as

$$P_{\alpha , \lambda }\left(\beta \right)={\int }_0^{\left|\beta \right|}{\left({\displaystyle \frac{\lambda -x}{\alpha }}\right)}_{+}dx=\left(\lambda \left|\beta \right|-{\displaystyle \frac{{\beta }^2}{2\alpha }}\right)\mathit{I}\left(0\le \left|\beta \right|<\lambda \alpha \right)+{\displaystyle \frac{{\lambda }^2\gamma }{2}}\mathit{I}\left(\left|\beta \right|\ge \lambda \alpha \right)$$

4.2. Shrinkage Estimators

Shrinkage estimation offers a robust method for parameter estimation, especially in sparse environments where traditional methods may face high variance or multicollinearity. This section introduces standard and positive shrinkage estimators to improve accuracy and stability. Shrinkage estimators balance full model estimators $\left({\widehat{\mathsf{\beta }}}^{FM}\right)$, which use extensive predictors, with sub model estimators $\left({\widehat{\mathsf{\beta }}}^{SM}\right)$, which focus on fewer, significant predictors using any of introduced penalty functions (see [2]). To better understand the shrinkage process, it is necessary to partition the design matrix X into subsets that correspond to significant and less significant predictors. Let X represent the n × k design matrix, where n is the number of observations, and k is the total number of predictors. We partition X into two submatrices: $\tilde{\mathbf{X}}= [\tilde{\mathbf{X}}₁, \tilde{\mathbf{X}}₂]$ where ${\tilde{\mathbf{X}}}_1$ corresponds to the subset of significant nonzero predictors ($p_1$ predictors) and ${\tilde{\mathbf{X}}}_2,$contains the remaining predictors ($p-p_1$) where $p=p_1+p_2$. This partitioning ensures that the shrinkage process selectively regularizes the coefficient associated with less significant predictors while preserving the interpretability of the model.

After optimizing ${\widehat{\mathsf{\beta }}}^{FM}$ and ${\widehat{\mathsf{\beta }}}^{SM}$ we can define the shrinkage estimator ${\widehat{\mathsf{\beta }}}_S$ and positive part of the shrinkage estimator to avoid the over-shrinking of the regression coefficients ${\widehat{\mathsf{\beta }}}_{PS}$ as follows:

$$\begin{gathered} {\widehat{\mathsf{\beta }}}_S={\widehat{\mathsf{\beta }}}^{FM}-{\displaystyle \frac{\left(p_2-2\right)}{\widehat{T}\left({\widehat{\mathsf{\beta }}}_{FM}-{\widehat{\mathsf{\beta }}}_{SM}\right)}} \\ {\widehat{\mathsf{\beta }}}_{PS}={\widehat{\mathsf{\beta }}}^{FM}-{\left[{\displaystyle \frac{\left(p_2-2\right)}{\widehat{T}}}\right]}_1\left({\widehat{\mathsf{\beta }}}_{FM}-{\widehat{\mathsf{\beta }}}_{SM}\right) \end{gathered}$$

(4.8)

where ${\left(\eta \right)}_1=\max \left[\eta ,1\right]$, $p_2$ represents the number of sparse parameters that are detected by the penalty functions to estimate ${\widehat{\mathsf{\beta }}}^{SM}$, $\widehat{\mathbf{T}}$ is the distance measure which can be defined as $\widehat{\mathbf{T}}={\sigma }^{-2}n{\widehat{\mathsf{\beta }}}_2^{\text{'}}\left({\tilde{\mathbf{X}}}_2^{\mathit{\text{'}}}\tilde{\mathbf{U}}{\tilde{\mathbf{X}}}_2\right){\widehat{\mathsf{\beta }}}_2$ with sparse subset of regression coefficients denoted by ${\widehat{\mathsf{\beta }}}_2$, model variance ${\sigma }^2$, ${\tilde{\mathbf{X}}}_2$ matrix of partial residuals of associated predictors with ${\widehat{\mathsf{\beta }}}_2$ and the projection matrix $\tilde{\mathbf{U}}={\mathbf{I}}_n-\left[{\tilde{\mathbf{X}}}_1{\left({\tilde{\mathbf{X}}}_1^{\mathit{\text{'}}}{\tilde{\mathbf{X}}}_1+{\mathsf{\lambda }}_{\mathrm{s}}{\mathbf{I}}_{p_1}\right)}^{-1}{\tilde{\mathbf{X}}}_1^{\mathit{\text{'}}}\right]$ where ${\tilde{\mathbf{X}}}_1$ is the partial residuals associated with the $p_1$ significant coefficients. Finally, ${\lambda }_s$ is again shrinkage parameter. Notice that some of the advantages of using shrinkage estimators can be mentioned as follows: introducing the bias through the shrinkage, these estimators reduce variance and leading the smaller Mean Squared Error (MSE) values which means better prediction accuracy as well. In addition, shrinkage techniques shrink the coefficients in a controlled manner which preserved the interpretability.

Building upon the Shrinkage Estimation framework introduced above, we integrate both standard and Positive Shrinkage Estimators into the penalty estimation strategies. These estimators provide an additional layer of regularization, enhancing the estimation accuracy for sparse models. The estimation procedure involves:

• Parametric Component Estimation: We utilize the respective penalty functions (Ridge, Lasso, aLasso, SCAD, ElasticNet, MCP) to estimate the parameter vector $\mathsf{\beta }$.
• Shrinkage Application: We apply the shrinkage techniques to obtain ordinary stein-type shrinkage ${\widehat{\mathsf{\beta }}}_S$ or positive shrinkage estimator ${\widehat{\mathsf{\beta }}}_{PS}$ to refine the parameter estimates.
• Nonparametric Component Estimation: We estimate the smooth function $f\left(t\right)$ using the refined parameter estimates via kernel smoothing as described in equation (4.8).

This comprehensive approach leverages the strengths of both penalty functions and shrinkage estimators, offering a robust estimation framework for PLRM. The incorporation of shrinkage estimators offers several advantages. By introducing a controlled bias, these estimators reduce the variance of parameter estimates, leading to more stable and reliable predictions. Additionally, they enhance model interpretability by focusing on significant predictors while mitigating the effects of multicollinearity. Moreover, the positive shrinkage estimator provides a practical safeguard against excessive regularization, striking a balance between bias and variance.

4.3. Estimation Procedure for the Nonparametric Component

For the parametric part of the semiparametric model in (1.2), modified kernel smoothing estimates based on various penalties are provided by equations (4.2– 4.7). Analogous to that of (3.3b), the estimation of the nonparametric part of the same model can be constructed via the vector of estimated parametric coefficients ${\widehat{\mathsf{\beta }}}_{LK}$ provided in (4.3). On the basis of the Lasso penalty of the unknown function, we obtain the transformed local estimates as follows:

$${\widehat{\mathbf{f}}}_{LK}={\mathbf{W}}_h\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{LK}\right)$$

(4.8)

as described in the previous section. Importantly, local estimates of the nonparametric component based on the adaptive-Lasso penalty are derived and are symbolically given as ${\widehat{\mathbf{f}}}_{aLK}$ in equation (4.8), where ${\widehat{\mathsf{\beta }}}_{aLK}$ defined in (4.4) is written in place of ${\widehat{\mathsf{\beta }}}_{LK}$. Similarly, for the nonparametric part of the semiparametric model (1.2), the modified local kernel estimators ${\widehat{\mathbf{f}}}_{SK}$, ${\widehat{\mathbf{f}}}_{ENK},$, and ${\widehat{\mathbf{f}}}_{MCK}$ based on the SCAD, ElasticNet, and MCP penalties, respectively, are obtained by replacing ${\widehat{\mathsf{\beta }}}_{LK}$ in (4.8) with ${\widehat{\mathsf{\beta }}}_{SK},$, ${\widehat{\mathsf{\beta }}}_{ENK}$, and ${\widehat{\mathsf{\beta }}}_{MCK}$.

To better understand the adaptation of the penalty functions to the estimation process of the partially linear model (1.1), a generic algorithm is introduced. As mentioned previously, to obtain a specific estimator on the basis of a determined penalty function, the corresponding penalty term is used in Algorithm 1. Therefore, the output estimators are denoted as $\left(\widehat{\mathsf{\beta }},\widehat{\mathbf{f}}\right)$ to generalize the algorithm. For example, if we obtain the Lasso estimator of the parametric components in model (1.1), first, ${\widehat{\mathsf{\beta }}}_{LK}$ is found to minimize (4.2), then ${\widehat{\mathbf{f}}}_{LK}$ is obtained by using the vector of estimated parametric coefficients ${\widehat{\mathsf{\beta }}}_{LK}$. A similar path is followed for the remaining estimators. The algorithm describing their findings is as follows:

Table 1. Computation of penalty estimators.

Input: Data matrix of parametric component $\mathbf{X}\in {\mathbb{R}}^{n\times k}$ , data vector $t_i\in {\mathbb{R}}^{n\times 1}$, and response vector $\mathbf{y}\in {\mathbb{R}}^{n\times 1}$. Output: Pair of estimates $\left(\widehat{\mathsf{\beta }},\widehat{\mathbf{f}}\right)$ based on a certain penalty function 1: Select an appropriate bandwidth $h$ using a predetermined criterion and compute the smoother matrix ${\mathbf{W}}_h$, as defined in (3.9). 2: Compute the partial residuals $\tilde{\mathbf{X}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathbf{X}$ and $\tilde{\mathbf{y}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathbf{y}$. 3: To minimize $PLS\left(\mathsf{\beta },\mathbf{f}\right),$determine the shrinkage parameter $\lambda$ by a predetermined criterion. 4: Partition the partial residuals of $\mathbf{X}$ in form $\tilde{\mathbf{X}}=({\tilde{\mathbf{X}}}_1,{\tilde{\mathbf{X}}}_2)$, as defined in Section 4.2. 5: Apply shrinkage estimators ${\widehat{\mathsf{\beta }}}^{PS}$, and ${\widehat{\mathsf{\beta }}}^S$ based on used penalty functions 6: Find the estimate of parametric component $\widehat{\mathsf{\beta }}$ associated with the $\tilde{\mathbf{X}}$ contains and $\tilde{\mathbf{y}}$. 7: Estimate the nonparametric smooth function $\mathbf{f}$ as follows: $\widehat{\mathbf{f}}={\mathbf{W}}_h\left(\mathbf{y}-\mathbf{X}\widehat{\mathsf{\beta }}\right)$ 8: Return $\left(\widehat{\mathsf{\beta }},\widehat{\mathbf{f}}\right)$.

Table 1. Computation of penalty estimators.

Input: Data matrix of parametric component $\mathbf{X}\in {\mathbb{R}}^{n\times k}$ , data vector $t_i\in {\mathbb{R}}^{n\times 1}$, and response vector $\mathbf{y}\in {\mathbb{R}}^{n\times 1}$. Output: Pair of estimates $\left(\widehat{\mathsf{\beta }},\widehat{\mathbf{f}}\right)$ based on a certain penalty function 1: Select an appropriate bandwidth $h$ using a predetermined criterion and compute the smoother matrix ${\mathbf{W}}_h$, as defined in (3.9). 2: Compute the partial residuals $\tilde{\mathbf{X}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathbf{X}$ and $\tilde{\mathbf{y}}=\left({\mathbf{I}}_n-{\mathbf{W}}_h\right)\mathbf{y}$. 3: To minimize $PLS\left(\mathsf{\beta },\mathbf{f}\right),$determine the shrinkage parameter $\lambda$ by a predetermined criterion. 4: Partition the partial residuals of $\mathbf{X}$ in form $\tilde{\mathbf{X}}=({\tilde{\mathbf{X}}}_1,{\tilde{\mathbf{X}}}_2)$, as defined in Section 4.2. 5: Apply shrinkage estimators ${\widehat{\mathsf{\beta }}}^{PS}$, and ${\widehat{\mathsf{\beta }}}^S$ based on used penalty functions 6: Find the estimate of parametric component $\widehat{\mathsf{\beta }}$ associated with the $\tilde{\mathbf{X}}$ contains and $\tilde{\mathbf{y}}$. 7: Estimate the nonparametric smooth function $\mathbf{f}$ as follows: $\widehat{\mathbf{f}}={\mathbf{W}}_h\left(\mathbf{y}-\mathbf{X}\widehat{\mathsf{\beta }}\right)$ 8: Return $\left(\widehat{\mathsf{\beta }},\widehat{\mathbf{f}}\right)$.

5. Measuring the Quality of Estimators

This section describes several performance metrics that can be used to assess the performance of the modified semiparametric kernel smoothing estimators based on penalty functions and shrinkage estimators defined in Section 4. Note that the estimators are indicated by the acronyms in parentheses. Each estimator's performance for the parametric component, nonparametric component, and overall estimated model is examined. Consequently, the performance metrics are explained in the following section.

5.1. Evaluation of the Parametric Component

The performance of a method is related to its ability to estimate the data. Note that the bias and variance of a penalty estimator ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ are measured simultaneously via the mean squared error $\left(MSE\right)$ matrix. Some metrics are defined in the following section.

MSE: This value is calculated by adding the variance and squared bias, which is provided by

$$\mathit{M}\mathit{S}\mathit{E}\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}},\mathsf{\beta }\right)=E\left[{\left(\mathsf{\beta }-{\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)}^{\text{'}}(\mathsf{\beta }-{\widehat{\mathsf{\beta }}}_{\mathrm{S}})\right]=Var\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)+{(E\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)-\mathsf{\beta })}^2$$

(5.1)

where ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ is the vector of the estimated regression coefficients determined by either of the methods. The key information regarding the quality standards for estimation is provided by equation (4.3). To evaluate the risk more accurately via the MSE, a scalar variant of the MSE (SMSE) is utilized, which yields an estimator's average quality score:

$$\mathit{S}\mathit{M}\mathit{S}\mathit{E}\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}},\mathsf{\beta }\right)= E\left[{\left(\mathsf{\beta }-{\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)}^{\text{'}}\left(\mathsf{\beta }-{\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)\right]=tr\left[\left(MSE\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}},\mathsf{\beta }\right)\right)\right]$$

(5.2)

Additionally, the relative efficiencies of the methods for estimating $\mathsf{\beta }$ can be easily obtained via (5.2). To do this, the relative efficiency (RE) is defined in Definition 5.1.

Definition 5.1: Let $S_1$ and $S_2$ denote the two distinct estimation methods for $\mathsf{\beta }$, and the estimators are shown as ${{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_1$ and ${{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_2,$, respectively. The formula for RE is as follows:

$$RE\left({{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_1,{{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_2\right)=SMSE\left({{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_1,\mathsf{\beta }\right)/SMSE({{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_2,\mathsf{\beta })$$

(5.3)

Here, $RE\left({{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_1,{{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_2\right)>1$ means that ${{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_2$ is more efficient than ${{\widehat{\mathsf{\beta }}}_{\mathrm{S}}}_1$.

Root mean squared error ($\mathbf{R}\mathbf{M}\mathbf{S}\mathbf{E}$): The RMSE is simply the square root of the MSE. The formula for calculating the estimated regression coefficients $\left(\widehat{\mathsf{\beta }}\right)$ is as follows:

$$MSE\left(\mathsf{\beta },{\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)=\sqrt{\mathrm{M}\mathrm{S}\mathrm{E}\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}},\mathsf{\beta }\right)}$$

(5.4)

where ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ is the estimate of $\mathsf{\beta }$ obtained via one of the methods described in this paper, as defined in (5.2). To obtain the RMSE score for each estimator, ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ is replaced with ${\widehat{\mathsf{\beta }}}_{RK},{\widehat{\mathsf{\beta }}}_{LK},{\widehat{\mathsf{\beta }}}_{aLK},{\widehat{\mathsf{\beta }}}_{SK},$ ${\widehat{\mathsf{\beta }}}_{ENK}$ and ${\widehat{\mathsf{\beta }}}_{MCK}$

5.2. Evaluation of the Nonparametric Component

To evaluate the nonparametric part of the model (1.1), $MSE$ can also be used. Suppose that $\widehat{\mathbf{f}}$ is any of the estimators of the nonparametric component in model (1.1). That is, let us assume that it is equivalent to one of the estimators ${\widehat{\mathbf{f}}}_{RK},{\widehat{\mathbf{f}}}_{LK},{\widehat{\mathbf{f}}}_{aLK},{\widehat{\mathbf{f}}}_{SK},$ ${\widehat{\mathbf{f}}}_{ENK}$ and ${\widehat{\mathbf{f}}}_{MCK}$ defined in Section 3. Therefore, the MSE is calculated as follows:

$$MSE\left(f,\widehat{f}\right)=n^{-1}\sum _{i=1}^n{\left(f\left(t_i\right)-\widehat{f}\left(t_i\right)\right)}^2=n^{-1}{‖(\mathbf{f}-\widehat{\mathbf{f}})‖}^2$$

(5.8)

To compare how well each of the six methods estimates the nonparametric component, the relative $MSE$ ($RMSE$) is used. The formula for $RMSE$ is computed as follows:

$$RMSE\left({\widehat{\mathbf{f}}}_i\right)=n_m^{-1}\left[{\sum }_{j=1}^{(n_m-1)}MSE\left({\widehat{\mathbf{f}}}_i\right)/MSE\left({\widehat{\mathbf{f}}}_j\right)\right] , i\ne j$$

(5.9)

6. Asymptotic Analysis

Suppose that ${\mathsf{\beta }}_1^{\mathbf{*}}\ne 0$ and${\mathsf{\beta }}_2^{\mathbf{*}}=0$ is the estimate of $\mathsf{\beta }=({\mathsf{\beta }}_1\ne 0,{\mathsf{\beta }}_2=0)$, where ${\mathsf{\beta }}_1^{\mathbf{*}}\ne 0$ is any of the proposed estimators, which are Ridge, Lasso, aLasso, SCAD, ElasticNet, and MCP. This section defines the asymptotic distributional risk (ADR) of our full model and selected model estimators. Our primary objective is the performance of these estimators when ${\mathsf{\beta }}_2=0$; for this reason, we consider a sequence $\left\{{\mathrm{K}}_{\mathrm{n}}\right\}$ given by

$${\mathrm{K}\mathrm{n}:\mathsf{\beta }}_2={\mathsf{\beta }}_{2\left(\mathbf{n}\right)}={\displaystyle \frac{\mathbf{w}}{\sqrt{\mathbf{n}}}},\mathrm{w}={({\mathrm{w}}_1,...,{\mathrm{w}}_{\mathit{p}2})}^{\mathit{\text{'}}}\in {\mathbb{R}}^{\mathit{p}2}$$

Now, using a positive definite matrix (p.d.m.) W, we define a quadratic loss function as follows:

$$\mathrm{L}\left({\mathsf{\beta }}_1^{\mathbf{*}}\right)=\mathbf{n}{({\mathsf{\beta }}_1^{\mathbf{*}}-{\mathsf{\beta }}_1)}^{\mathbf{\text{'}}}\mathbf{W}\left({\mathsf{\beta }}_1^{\mathbf{*}}-{\mathsf{\beta }}_1\right),$$

Now, under${\mathrm{K}}_{\mathrm{n}}$, we can define the asymptotic distribution function of ${\mathsf{\beta }}_1^{\mathbf{*}}$ as

$$F\left(\mathbf{x}\right)=\underset{n\to \infty }{\lim }P (\sqrt{n} ({\mathsf{\beta }}_1^{\mathbf{*}}-{\mathsf{\beta }}_1)\le \mathbf{x}|{\mathrm{K}}_{\mathrm{n}}),$$

where $F \left(\mathbf{x}\right)$ is nondegenerate. Hence, we can define the ADR of ${\mathsf{\beta }}_1^{\mathbf{*}}$ as follows.

$$\mathrm{A}\mathrm{D}\mathrm{R}\left({\mathsf{\beta }}_1^{\mathbf{*}}\right)=\mathrm{t}\mathrm{r}\left(\mathbf{W}{\int }_{{\mathbb{R}}^{\mathit{P}1}}\int {\mathbf{x}\mathbf{x}}^{\mathbf{\text{'}}}d F \left(\mathbf{x}\right)\right)=\mathrm{t}\mathrm{r}\left(\mathbf{W}\mathbf{V}\right)$$

where V is the dispersion matrix for the distribution $F\left(\mathbf{x}\right)$。

Assumption 6.1. 

We establish two regularity conditions as follows:

(i)

${\displaystyle \frac{1}{n}}\underset{0\le \mathrm{x}\le 1}{\max }{\tilde{\mathbf{x}}}_{\mathbf{i}}^{\mathbf{\text{'}}}{ \left({\tilde{\mathbf{X}}}_{\mathbf{i}}^{\mathbf{\text{'}}}\tilde{\mathbf{X}}\right)}^{-1}{\tilde{\mathbf{x}}}_{\mathbf{i}} \to 0\mathrm{a}\mathrm{s} n\to$ $\infty$ where ${\tilde{\mathbf{x}}}_{\mathbf{i}}^{\mathbf{\text{'}}}$ is the ith row of $\tilde{\mathbf{X}}$

(ii)

${\displaystyle \frac{1}{n}} \sum _{i=1}^n{\tilde{\mathbf{X}}}_{\mathbf{i}}^{\mathbf{\text{'}}}\tilde{\mathbf{X}}\to \widetilde{\mathbf{Q},}$where $\widetilde{\mathbf{Q},}$ is a finite positive-definite matrix.

By Lemma 6.1, which is defined below, with assumed regularity conditions and local alternatives, the ADRs of the estimators under $\left\{K_n\right\}$ are given as:

Theorem 6.1. 

ADR (${\widehat{\mathit{\beta }}}_1^{\mathit{F}\mathit{M}}$) =${\sigma }^2 tr\left(\mathit{W} {\widetilde{\mathit{Q},}}_{11.2}^{-1}\right)+{\mathit{\eta }}_{11.2}^{\text{'}} \mathit{W}{\mathit{\eta }}_{11.2},$ADR$({\widehat{\mathit{\beta }}}_1^{\mathit{S}\mathit{M}}$) =${\sigma }^2 tr\left(\mathit{W} {\widetilde{\mathit{Q},}}_{11}^{-1}\right)$+$\xi \text{'}\mathit{W}\xi$.

Lemma 6.1 enables us to derive the result of Theorem 6.1 in this study.

Lemma 6.1. 

$ifk/\sqrt{n}\to {\lambda }_0 \ge 0 and \widetilde{\mathit{Q},}is non-singular,then$

$$\sqrt{n} \left({\widehat{\mathit{\beta }}}_1^{\mathit{F}\mathit{M}}-\mathit{\beta }\right)\stackrel{d}{\to }\left(-{\lambda }_0{\tilde{\mathit{Q}}}^{-1}\mathit{\beta },{\sigma }^2 {\tilde{\mathit{Q}}}^{-1}\right).$$

where “$\stackrel{d}{\to }$” denotes the convergence in the distribution. See the proof in [2].

7. Simulation Studies

In this section, we conduct a thorough simulation study to assess the finite-sample performance of the aforementioned six semiparametric estimators proposed for a partially linear model. These estimators are compared to each other to evaluate their performance. R-software [34] was used for all calculations. A description of the simulation design and data generation is as follows.

Design of the Simulation: The simulations are carried out as follows.

i.

Samples of size n = 50, 100 and 200

ii.

Two numbers of parametric covariates k = 25 and 40

iii.

Two correlation levels ρ(rho) = 0.5 and 0.90

iv.

The number of replications is 1000.

In addition, we generated zero and nonzero coefficients of the model to enable our absolute penalty functions to select significant variables. As mentioned above, each possible simulation design was repeated 1000 times to detect significant findings. Finally, the evaluation metrics mentioned in Section 5 were used to assess the performance of the proposed estimators.

Data generation: We generated the explanatory variables ${\mathbf{x}}_i=x_{ij}=(x_{i1},\dots ,x_{ik})$ via the following equation:

$$x_{ij}={(1-{\rho }^2)}^{1/2}{x}_{ij}+\rho x_{i,k+1}, i=\mathrm{1,2},\dots ,n, j=\mathrm{1,2},\dots ,k$$

where $x_{ij}$ are independent standard normal pseudorandom numbers and where $\rho$ represents the level of collinearity among the explanatory variables. The two levels of collinearity assumed in this study are 0.50 and 0.90, as denoted above. These variables are standardized so that $\mathbf{X}\mathbf{\text{'}}\mathbf{X}$ and $\mathbf{X}\mathbf{\text{'}}\mathbf{y}$ are in correlation forms. The $n$ observations for the response variable $y_i$ are generated by

$$y_i={\mathbf{x}}_i^{\text{'}}\mathsf{\beta }+f\left(t_i\right)+{\epsilon }_i, 1\le i\le n$$

where $\mathsf{\beta }={\left({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4,{\beta }_5,{{\beta }_6,\dots ,\beta }_{25}\right)}^T={\left(-1,2,0.5,3,-2,0,0,\dots ,2\right)}^T$, the nonparametric function $f\left(t_i\right)={-t}_i\sin \left(t_i\right)$ with$t_i~U\left(\mathrm{0,1}\right)$, and ${\epsilon }_i~N\left(\mathrm{0,1}\right)$. It is understood from this model that when $k=25$ explanatory variables are considered, there are 6 nonzero ${\beta }_j$s to be estimated and $\left(k-6\right)=19$ sparse coefficients. If $k=40$ explanatory variables are considered, there will be 6 nonzero ${\beta }_j$s to be estimated and $\left(k-6\right)=34$ sparse coefficients.

In this study, to obtain the optimal value of $\lambda$, we consider modified cross-validation (CV), which is generalized cross-validation (GCV). The main idea of the $GCV$ criterion is to replace the factor $\left(1-{\left({\mathbf{W}}_h\right)}_{ii}\right)$ with the mean value, that is, $\left(1-n^{-1}tr\left({\mathbf{W}}_h\right)\right)$. In this case, the $GCV$ value is obtained as follows Craven and Wahba (1979):

$$GCV\left(h\right)={\displaystyle \frac{n{‖\left(\mathbf{I}-{\mathbf{W}}_h\right)\mathbf{y}‖}^2}{{\left[tr(\mathbf{I}-{\mathbf{W}}_h)\right]}^2}}={\displaystyle \frac{n{‖\left(\mathbf{I}-{\mathbf{W}}_h\right)\mathbf{y}‖}^2}{{\left[n-tr\left({\mathbf{W}}_h\right)\right]}^2}}$$

(7.1)

where $tr\left({\mathbf{W}}_h\right)$ represents the sum of the diagonal values of the matrix ${\mathbf{W}}_h$, that is, its trace. Finally, to obtain the optimum shrinkage parameter $\lambda$ in the objective function given in (4.1), we use a 5-fold cross-validation method.

In this simulation study, we analyzed 12 different simulation configurations and attempted to display all of them. In this sense, the simulation study results are shown separately for the parametric and nonparametric components.

7.1. Analysis of the Parametric Component

This subsection examines the estimation of the parametric component of a semiparametric model, with the results shown in Table 3 for all the simulations. In addition to the estimation, the performance of the estimators was evaluated via the performance metrics mentioned in Section 5. The results are presented in the following tables and figures.

Table 2, Table 3 and Table 4 detail the performance scores (MSE, RMSE, RE, SMSE) for the parametric component across different simulation setups.

In terms of general observations that can be concluded from the tables; as the sample size (n) increases, the MSE and RMSE values generally decrease for all estimators, indicating improved accuracy. Higher multicollinearity ($\rho$ = 0.9) tends to worsen the performance of all estimators compared to lower multicollinearity ($\rho$ = 0.5). The MCP and SCAD penalty functions, along with the positive shrinkage estimator (${\widehat{\mathsf{\beta }}}_{PS}$) and shrinkage estimator (${\widehat{\mathsf{\beta }}}_S$), often demonstrate superior performance in terms of lower MSE and RMSE, especially in scenarios with high multicollinearity or smaller sample sizes.

In detail, Ridge (${\widehat{\mathsf{\beta }}}_{RK}$), performs relatively poorly compared to other methods, particularly when multicollinearity is high. Lasso (${\widehat{\mathsf{\beta }}}_{LK}$) and aLasso (${\widehat{\mathsf{\beta }}}_{aLK}$), these methods show moderate performance. aLasso tends to be slightly better than Lasso. ElasticNet (${\widehat{\mathsf{\beta }}}_{ENK}$) performs similarly to Lasso and aLasso. MCP (${\widehat{\mathsf{\beta }}}_{MCK}$) and SCAD (${\widehat{\mathsf{\beta }}}_{SK}$), often exhibit the best performance among the penalty methods, with lower MSE and RMSE values, especially in challenging scenarios. Shrinkage Estimators (${\widehat{\mathsf{\beta }}}_{PS}$ and ${\widehat{\mathsf{\beta }}}_S$), both shrinkage estimators consistently perform well, often outperforming the other methods, especially when the sample size is small and multicollinearity is high.

Figure 1 illustrates how the Mean Squared Error (MSE) of the parametric component estimators changes with varying levels of multicollinearity ($\rho$). Notably, as multicollinearity increases, so does the MSE for all estimators, which aligns with the data presented in Table 2, Table 3 and Table 4. However, the ${\widehat{\mathsf{\beta }}}_{MCK}$, ${\widehat{\mathsf{\beta }}}_{SK}$, and the shrinkage estimators, particularly the ${\widehat{\mathsf{\beta }}}_{PS}$, consistently exhibit lower MSE values compared to ${\widehat{\mathsf{\beta }}}_{RK}$, ${\widehat{\mathsf{\beta }}}_{LK}$, ${\widehat{\mathsf{\beta }}}_{aLK}$, and ${\widehat{\mathsf{\beta }}}_{ENK}$. This trend is particularly pronounced when multicollinearity is high (ρ = 0.9). The superior performance of MCP, SCAD, and shrinkage estimators under these challenging conditions is further emphasized by the relatively flat lines they exhibit in Figure A3 in Appendix A.3, suggesting robustness to an increasing number of predictors.

7.2. Analysis of the Nonparametric Component

This subsection examines the estimation of the nonparametric component of a semiparametric model, with the results shown in Table 5 for all the simulations. Furthermore, the performance scores of the proposed methods were measured via the MSE and RMSE metrics. In addition, the estimated curves obtained from all methods are examined using all the distinct configurations in Figure 2 and Figure 3.

Table 5 presents the performance scores for the nonparametric component estimators in a similar format to Table 2, Table 3 and Table 4, which focused on the parametric component. This table would quantify the accuracy and efficiency of each estimator in estimating the nonparametric function $f\left(t\right)$ under various simulation scenarios, varying sample size ($n$), number of predictors ($p$), and multicollinearity levels ($\rho$).

Based on the trends observed in Figure 2 and Figure 3, Figure A1, and Figure A4, and the consistent findings in Table 2, Table 3 and Table 4, it is highly probable that Table 5 would further demonstrate the superior performance of ${\widehat{\mathbf{f}}}_{MCK}$, ${\widehat{\mathbf{f}}}_{SK}$, and the shrinkage estimators (${\widehat{\mathbf{f}}}_S$ and ${\widehat{\mathbf{f}}}_{PS}$), particularly in scenarios with high multicollinearity and smaller sample sizes. These estimators would likely exhibit lower MSE and RMSE values compared to others, indicating greater accuracy in estimating the nonparametric function, and potentially higher RE values, suggesting greater efficiency relative to a reference estimator like${\widehat{\mathbf{f}}}_{RK}$. The data in Table 5 would provide the numerical backbone to support the visual observations in Figure 2 and the trends shown in Figure 3, Figure A1, and Figure A4, solidifying the conclusion that${\widehat{\mathbf{f}}}_{MCK}$, ${\widehat{\mathbf{f}}}_{SK}$, ${\widehat{\mathbf{f}}}_S$ and ${\widehat{\mathbf{f}}}_{PS}$ offer significant advantages when estimating the nonparametric component of partially linear models.

Figure 2 provides a visual representation of how well each estimator captures the true nonparametric function, $f\left(t\right).$ When examining the fitted curves, especially under conditions of high multicollinearity ($\rho = 0.9$) and small sample size ($n=50$), it becomes apparent that MCP (${\widehat{\mathbf{f}}}_{MCK}$), SCAD ${\widehat{\mathbf{f}}}_{SK}$, and the shrinkage estimators (${\widehat{\mathbf{f}}}_{PS} \& {\widehat{\mathbf{f}}}_S$) produce curves that align more closely with the true function (represented by the black line) than those generated by Ridge (${\widehat{\mathbf{f}}}_{RK}$), Lasso (${\widehat{\mathbf{f}}}_{LK}$), aLasso (${\widehat{\mathbf{f}}}_{aLK}$), and ElasticNet (${\widehat{\mathbf{f}}}_{ENK}$). As the sample size increases, the differences become less visually striking, yet the advantage of ${\widehat{\mathbf{f}}}_{MCK}$, ${\widehat{\mathbf{f}}}_{SK}$, and shrinkage estimators persists, particularly when multicollinearity is present. This visual observation is consistent with the quantitative data shown in Figure 3, where these estimators consistently exhibit lower MSE values for the nonparametric component.

Figure 3 complements Figure 2 by quantifying the MSE of the nonparametric component against varying levels of multicollinearity. The trend observed mirrors that of Figure 1, MSE increases with multicollinearity for all estimators. However, ${\widehat{\mathbf{f}}}_{RK}$, consistently displays the highest MSE, while${\widehat{\mathbf{f}}}_{MCK}$,, ${\widehat{\mathbf{f}}}_{SK}$,, and the shrinkage estimators maintain lower MSE values across the board. This quantitative evidence corroborates the visual findings of Figure 2, demonstrating that ${\widehat{\mathbf{f}}}_{MCK}$,, ${\widehat{\mathbf{f}}}_{SK}$,, and shrinkage estimators (${\widehat{\mathbf{f}}}_S$ and ${\widehat{\mathbf{f}}}_{PS}$,,) provide more accurate estimates of the nonparametric function, especially when faced with correlated predictors. The relatively stable RMSE values of these better-performing estimators, even as the number of predictors increases, are shown in Figure A4 and further support their robustness in Appendix A.3.

8. Real Data

In this section, we implement our proposed six estimators, Ridge, Lasso, aLasso, SCAD, ElasticNet, and MCP, on real data. Our real data, Hitter’s dataset, can be obtained from the ISLR package in R. This dataset contains 322 rows and 20 variables. Three covariates, Division, League and New League, which are not scalers, were deleted. The CN value obtained from this dataset is 5830, which indicates the presence of multicollinearity in the dataset. The remaining variables used for the estimation of a partially linear model are $Salary$,$Atbat, Hits, Hmrun, Runs, RBI, Walks, CatBat, CHits, CHmrun$,$Cruns$, $CRBI$, $Cwalks$ and $PutsOuts$. In addition, the performance of the estimators is compared, and the results are shown in Table 6, and Figure 4 and Figure 5. Note that the logarithm of Salary is used as a response variable. From the visual inspection of the relationships between covariate and the response variable, the variable “$Years$” has a significant nonlinear relationship with the response $\mathrm{l}\mathrm{o}\mathrm{g}\left(sal\right)$. The remaining 15 predictor variables are added to the parametric component of the model. Hence, the partially linear model can be written as follows:

$$log({sal}_i )={\mathbf{X}}_i^T\mathsf{\beta }+f({years}_i )+{\epsilon }_i,1\le i\le 263$$

(8.1)

where

$${\mathbf{X}}_i^T=\left[\begin{array}{c}{AtBat}_i, {Hits}_i,{HmRun}_i, {Runs}_i,{RBI}_i,{Walks}_i,{CAtBat}_i{,CHits}_i,{CHmRun}_i{CRuns}_i,\\ {CRBI}_i,{CWalks}_i,{PutOuts}_i,{Assists}_{ i},{Errors}_i\end{array}\right]$$

Therefore, there is an $\left(263\times 15\right)$-dimensional covariate matrix for the parametric component of the model, and $\mathsf{\beta }={\left({\beta }_1,\dots ,{\beta }_{15}\right)}^T$ is the $\left(15\times 1\right)$-dimensional vector of the regression coefficients to be estimated.

Table 6 quantifies the performance of each estimator on the Hitters dataset using RMSE, MSE, and RE, all based on the prediction of log(Salary). The Shrinkage and Positive Shrinkage estimators stand out with the lowest RMSE and MSE values of 0.413 and 0.171, and 0.405 and 0.164, respectively. This indicates that these two methods provide the most accurate predictions for log(Salary) on this dataset. Their superior performance is further underscored by their RE values, 1.950 and 2.025 respectively, which demonstrate that they are almost twice as efficient as the Ridge estimator, which has an RE of 1.000 and performs the worst on this metric.

aLasso also performs well, with an RMSE of 0.443, MSE of 0.196, and RE of 1.697, outperforming Lasso, SCAD, ElasticNet, and MCP. Lasso, SCAD, ElasticNet, and MCP exhibit moderate performance, with RMSE values ranging from 0.474 to 0.512 and RE values between 1.270 and 1.484. These results are consistent with the simulation study, which suggested that shrinkage estimators, particularly positive shrinkage, and aLasso are effective in handling multicollinearity and achieving good prediction accuracy.

Figure 4 visually depicts the estimated nonparametric function f$\left(Years\right)$ for each estimator, revealing how they model the nonlinear relationship between a player's standardized years of experience and their $log\left(Salary\right)$. The curves generally capture an initial increase in log(Salary) as years of experience increase, followed by a plateau or a slight decrease for players with more years. Notably, the curves produced by Ridge, MCP, and the shrinkage estimators are relatively similar and smooth, suggesting a consistent estimation of the underlying nonlinear trend. Lasso, aLasso, and ElasticNet exhibit some deviations from this trend, particularly in the middle range of the "Years" variable, indicating potential sensitivity to noise or specific data points. SCAD produces a curve that is very different from the others. The black dots, representing the observed data points adjusted for the parametric component's effect, provide a visual reference for assessing the fit of each curve. The variation in the shapes of the curves highlights the influence of the chosen estimation method on the estimated nonparametric relationship, emphasizing the importance of considering different estimators and their potential impact on the interpretation of the results.

Figure 5 provides a visual comparison of the predicted log(Salary) values against the actual $log\left(Salary\right)$ values for each estimator, offering a direct assessment of their predictive accuracy. The points for Ridge, MCP, and the shrinkage estimators are more closely clustered around the diagonal line ($y=x$), indicating that their predictions are generally closer to the actual values. This observation aligns with the lower RMSE and MSE values reported for these estimators in Table 6. Conversely, Lasso, aLasso, and ElasticNet exhibit a wider spread of points around the diagonal, suggesting greater prediction errors and confirming their higher RMSE and MSE values in the table. The plot for SCAD has a different pattern that does not align with the diagonal. Overall, Figure 5 visually reinforces the quantitative findings in Table 6, demonstrating that Ridge, MCP, and particularly the shrinkage estimators achieve better predictive performance on the Hitters dataset, while Lasso, aLasso, and ElasticNet show relatively weaker performance. This visualization underscores the practical implications of estimator selection for achieving accurate predictions in partially linear models.

9. Conclusions

This paper focuses on the realm of semiparametric regression, specifically focusing on Partially Linear Regression Models (PLRMs) that incorporate both linear and nonlinear components. It considers a comprehensive review of six penalty estimation strategies: Ridge, Lasso, aLasso, SCAD, ElasticNet, and MCP. Recognizing the challenges posed by multicollinearity and the need for robust estimation in sparse models, they further introduce Stein-type shrinkage estimation techniques. The core of the study lies in evaluating the performance of these methods through both theoretical analysis and empirical investigation. A kernel smoothing technique, grounded in penalized least squares, is employed to estimate the semiparametric regression models.

The theoretical contributions include partly asymptotic analysis of the proposed estimators, providing insights into their long-term behavior and robustness. The empirical investigation encompasses a simulation study that examines the estimators' performance under various conditions, including different sample sizes, numbers of predictors, and levels of multicollinearity. Finally, the practical applicability of these methods is demonstrated through a real data example using the Hitters dataset, where the estimators are used to model baseball players' salaries based on their performance metrics. The paper concludes that aLasso and shrinkage estimators exhibit superior performance in terms of prediction accuracy and efficiency in the presence of multicollinearity, especially the positive shrinkage, performing even better than aLasso.

In detail, the following points can be emphasized in terms of theoretical inferences, simulation and real data studies:

• The paper establishes the asymptotic properties of the proposed estimators, including Ridge, Lasso, aLasso, SCAD, ElasticNet, MCP, and the Stein-type shrinkage estimators, providing a theoretical foundation for their use in PLRMs.
• The theoretical results highlight the advantages of aLasso and shrinkage estimation, particularly in scenarios with high multicollinearity and sparsity.
• The simulation study demonstrates that aLasso and the shrinkage estimators, especially the positive shrinkage estimator, consistently outperform other methods in terms of lower Mean Squared Error (MSE) and Root Mean Squared Error (RMSE) for both the parametric and nonparametric components of the PLRM.
• The superior performance of aLasso and shrinkage estimators is more pronounced when the sample size is small and multicollinearity is high, confirming their robustness in challenging conditions.
• MCP and SCAD also exhibit strong performance in the simulations, often outperforming Ridge, Lasso, and ElasticNet, particularly when multicollinearity is present.
• The simulation results reveal that the choice of estimator can significantly impact the estimation of the nonparametric function, with aLasso and shrinkage estimators generally producing smoother and more accurate curves.
• The analysis of the Hitters dataset confirms the practical advantages of aLasso and shrinkage estimation, particularly positive shrinkage, in a real-world scenario with multicollinearity, as indicated by the high condition number.
• The shrinkage and aLasso estimators achieve the lowest RMSE and MSE values when predicting log(Salary), demonstrating their superior predictive accuracy compared to Ridge, Lasso, SCAD, ElasticNet, and MCP.
• The fitted nonparametric curves for the "Years" variable reveal interesting differences in how each estimator captures the nonlinear relationship between experience and salary, with aLasso and shrinkage estimators providing a balance between flexibility and smoothness.
• The real data results align with the findings of the simulation study, further supporting the use of aLasso and shrinkage estimation, especially positive shrinkage, in PLRMs when multicollinearity is a concern. Also, SCAD has unexpected results which needs more investigation.

While this study provides valuable theoretical and empirical insights into the performance of various penalty and shrinkage estimators for partially linear models, several limitations should be acknowledged. The simulation study uses a GCV for the bandwidth which might not be the best one. The analysis is also limited to a single real-world dataset (Hitters), and the simulation study, while comprehensive, does not cover all possible scenarios, such as different error distributions.

The computational cost of some estimators, particularly SCAD and MCP, is not explicitly addressed. Additionally, the choice of tuning parameters for the penalized estimators is obtained based on CV criterion, whereas in practice, these need to be carefully selected, potentially impacting the results. The paper's theoretical results are based on regularity conditions that might not always hold in real-world applications. While shrinkage estimators are introduced and applied, a more in-depth investigation into their properties and performance would be beneficial. These limitations highlight the need for further research and careful consideration when applying these methods in practice.

Appendix A.1. Proof of Theorem 3.1

The most important step in obtaining the ridge penalty-based kernel smoothing estimator is to calculate the partial residuals and minimize (3.1). From that, as mentioned before, $\tilde{\mathbf{X}}=\left(\mathbf{I}-\mathbf{W}\right)\mathbf{X}$ and $\tilde{\mathbf{y}}=\left(\mathbf{I}-\mathbf{W}\right)\mathbf{y}$ are calculated, and the minimization of (3.1) is realized on the basis of response variable $\tilde{\mathbf{y}}$ as follows:

$$\begin{gathered} PL{S}_{RK}={\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)}^{\mathit{\text{'}}}\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)+\lambda {\mathsf{\beta }}^{\mathit{\text{'}}}\mathsf{\beta } \\ ={\tilde{\mathbf{y}}}^{\mathit{T}}\tilde{\mathbf{y}}-{\tilde{\mathbf{y}}}^{\mathit{T}}\tilde{\mathbf{X}}\mathsf{\beta }-{\mathsf{\beta }}^{\mathit{T}}{\tilde{\mathbf{X}}}^{\mathit{T}}\tilde{\mathbf{y}}+{\mathsf{\beta }}^{\mathit{T}}{\tilde{\mathbf{X}}}^{\mathit{T}}\tilde{\mathbf{X}}\mathsf{\beta }+\lambda {\mathsf{\beta }}^{\mathit{\text{'}}}\mathsf{\beta } \\ If, {\displaystyle \frac{\partial L}{\partial \beta }}=0 \\ -2{\tilde{\mathbf{y}}}^{\mathit{T}}\tilde{\mathbf{X}}+2{\tilde{\mathbf{X}}}^{\mathit{T}}\tilde{\mathbf{X}}\mathsf{\beta }+2\lambda \mathsf{\beta }=0 \\ {\tilde{\mathbf{X}}}^{\mathit{T}}\tilde{\mathbf{X}}\mathsf{\beta }+\lambda \mathsf{\beta }={\tilde{\mathbf{y}}}^{\mathit{T}}\tilde{\mathbf{X}} \\ {\widehat{\mathsf{\beta }}}_{RK}={\left({\tilde{\mathbf{X}}}^{\mathit{T}}\tilde{\mathbf{X}}\right)}^{-1}{\tilde{\mathbf{X}}}^T{\tilde{\mathbf{y}}}^{\mathit{T}} \end{gathered}$$

Appendix A.2. Proof of Lemma 1

${\mathrm{V}}_{\mathrm{n}}\left(\mathbf{u}\right)$ can be defined as follows:

$$\sum _{i=1}^n\left[{\left({\widetilde{ϵ}}_i-\mathbf{u}\text{'}{\overline{\mathbf{x}}}_{\mathbf{i}} /\sqrt{n}\right)}^2-{{\widetilde{ϵ}}_i}^2\right]+k\sum _{j=1}^p\left[{\left|{\beta }_j+u_j/\sqrt{n}\right|}^2-{\left|{\beta }_j\right|}^2\right],$$

where $\mathbf{u}= {\left(u_1, . . . , u_p\right)}^{\text{'}}$. According to Knight and Fu (2000), it can be shown that

$$\sum _{i=1}^n\left[{\left({\widetilde{ϵ}}_i-\mathbf{u}\text{'}{\overline{\mathbf{x}}}_{\mathbf{i}} /\sqrt{n}\right)}^2-{{\widetilde{ϵ}}_i}^2\right]\stackrel{d}{\to }-2{\mathbf{u}}^{\text{'}}\mathbf{D}+\mathbf{u}\text{'}\tilde{\mathbf{Q}}\mathbf{u}$$

where $\mathbf{D}~\mathcal{N}\left(0,{\sigma }^2 {\mathbf{I}}_P\right)$, with finite-dimensional convergence holding trivially. Hence,

$$k\sum _{j=1}^p\left[{\left|{\beta }_j+u_j/\sqrt{n}\right|}^2-{\left|{\beta }_j\right|}^2\right]\stackrel{d}{\to } {\lambda }_0\sum _{j=1}^p{u}_j\mathrm{s}\mathrm{g}\mathrm{n}\left({\beta }_j\right)\left|{\beta }_j\right|.$$

Hence, ${\mathbf{V}}_{\mathrm{n}}\left(\mathbf{u}\right) \stackrel{d}{\to } \mathbf{V}\left(\mathbf{u}\right)$. According to Geyer (1996), because ${\mathbf{V}}_n$ is convex and $\mathbf{V}$ has a unique minimum, it yields

$$\arg \min \left({\mathbf{V}}_n\right)=\sqrt{n} \left({\widehat{\mathsf{\beta }}}_1^{\mathrm{F}\mathrm{M}}-\mathit{\beta }\right)\stackrel{d}{\to }\arg \min \left(\mathbf{V}\right).$$

Hence,

$$\sqrt{n}\left({\widehat{\mathsf{\beta }}}_1^{\mathrm{F}\mathrm{M}}-\mathit{\beta }\right)\stackrel{d}{\to } {\tilde{\mathbf{Q}}}^{-1}\left(\mathit{D}-{\lambda }_0\mathit{\beta }\right)~\mathcal{N}\left({\lambda }_0{\tilde{\mathbf{Q}}}^{-1}\mathit{\beta },{\sigma }^2{\tilde{\mathbf{Q}}}^{-1}\right)$$

Appendix A.3. Additional Figures for Simulation and Real Data Studies

The boxplots in Figure A1 provide a detailed distributional view of the MSE for the nonparametric component. Here, we see that MCP, SCAD, and the shrinkage estimators not only have lower median MSE values but also smaller interquartile ranges, indicating greater consistency and less susceptibility to extreme errors. This is in contrast to methods like Ridge, Lasso, aLasso, and ElasticNet, which exhibit more outliers, especially under challenging conditions. These observations align with the findings presented in Figure 3 and further solidify the superior performance of MCP, SCAD, and shrinkage estimators in terms of both accuracy and precision when estimating the nonparametric part of the model.

Figure A1. Boxplots of MSE values obtained for nonparametric component of the model for all simulation configurations and all introduced estimators.

Figure A1. Boxplots of MSE values obtained for nonparametric component of the model for all simulation configurations and all introduced estimators.

Figure A2 shifts the focus to the relative efficiency (RE) of the estimators for the parametric component, comparing them to a reference estimator, likely Ridge. The plots reveal that MCP, SCAD, and the shrinkage estimators consistently demonstrate RE values greater than 1, often significantly so, particularly when the sample size is small and multicollinearity is high. This indicates that these methods are substantially more efficient than Ridge in estimating the parametric component under such conditions, a finding that resonates with the lower MSE values these estimators exhibit in Figure 1 and Table 2, Table 3 and Table 4.

Figure A2. Plots of RE values for estimated parametric components of the model including all simulation configurations.

Figure A2. Plots of RE values for estimated parametric components of the model including all simulation configurations.

Figure A3 examines the impact of the number of predictors on the RMSE of the parametric component. Notably, the lines for MCP, SCAD, and the shrinkage estimators remain relatively flat as the number of predictors increases. This suggests that these methods are less sensitive to model complexity compared to Ridge, Lasso, aLasso, and ElasticNet, which show a more pronounced upward trend in RMSE, especially when multicollinearity is present. This observation reinforces the robustness of MCP, SCAD, and shrinkage estimators, aligning with their overall superior performance highlighted in previous figures and tables.

Figure A3. RMSE values for estimated parametric component against the number of parameters (p) for all simulation configurations and introduced estimators.

Figure A3. RMSE values for estimated parametric component against the number of parameters (p) for all simulation configurations and introduced estimators.

Finally, Figure A4 investigates the relationship between RMSE and the number of predictors for the nonparametric component. Similar to Figure A3.3, the plots show that MCP, SCAD, and the shrinkage estimators exhibit more stable RMSE values as the number of predictors increases, particularly compared to Ridge, Lasso, aLasso, and ElasticNet. This demonstrates the robustness of these estimators when estimating the nonparametric part of the model, even as the complexity of the parametric component grows. The findings in this figure are consistent with the lower MSE values for the nonparametric component observed for MCP, SCAD, and shrinkage estimators in Figure 3 and visually confirmed in Figure 2.

Figure A4. RMSE values for estimated non-parametric component against the number of parameters (p) for all simulation configurations and introduced estimators.

Figure A4. RMSE values for estimated non-parametric component against the number of parameters (p) for all simulation configurations and introduced estimators.