Penalized Spline Estimation for Nonparametric Multiplicative Regression Models

Wei Chen^1*, Mingzhen Wan², Jiahui Xu², Jing Zhong², Yuying Xia², Minmin Zhang²

¹School of Zhangjiagang, Jiangsu University of Science and Technology, Zhangjiagang, China.
²Suzhou Institute of Technology, Jiangsu University of Science and Technology, Zhangjiagang, China.

DOI: 10.4236/oalib.1111352                

Abstract

In this paper, we consider the estimation problem of the unknown link function in the nonparametric multiplicative regression model. Combining the penalized splines technique, the least product relative error estimation method is proposed, where an effective model degree of freedom is de-fined, then the smoothing parameter is chosen by some information criteria. Simulation studies show that these strategies work well. Some asymptotic properties are established. A real data set is analyzed to illustrate the usefulness of the proposed approach. Finally, some possible extensions are discussed.

Keywords

Multiplicative Regression, Nonparametric Estimation, Penalized Splines, Relative Error, Smoothing Parameter Selection

1. Introduction

In many applications, such as the studies of financial and biomedical data, the response variable usually is positive. For modelling the relationship between the positive response and a set of explanatory variables, a natural idea is that first take an appropriate transformation for the response, e.g., the logarithmic transformation, then some common regression models, such as the linear regression or quantile regression, can be employed based on the transformed data. As argued by [1] , the least square or least absolute deviation criteria are both based on absolute errors, which is not desirable in many practical applications. Rather, the relative errors are more of a concern. In the early literature, many authors have contributed fruitfully to this issue, see [2] [3] and [4] . Especially, since the work of [1] , more attention has focused on the multiplicative regression (MR) model, and various extensions have been investigated, for example, see [5] [6] [7] and references therein.

It is worth noting that all the existing studies on MR are in the framework of linear (parametric) or semi-parametric MR models. Once the parametric forms or structures are misspecified, the resulting estimation and inference may be biased and the resulting conclusions become unreliable. Rather, nonparametric modelling is conceptually appealing and more robust. To the best of our knowledge, there are seldom studies on this problem. To fill this gap, we will address this problem in detail in this paper.

When estimating the nonparametric function g(z) in some semi-parametric MR models, such as the partially linear MR model ( [6] [8] [9] ), single index MR model ( [10] [11] [12] ), varying coefficient MR model ( [13] ), and others ( [14] ), almost all the researchers use the local linear smoothing technique and approximate it in a neighborhood of z for obtaining its estimation, where a good choice of the bandwidth is quietly critical and its value is sensitive to the performance of the resulting estimation and inference. Besides, due to the fact that the value of the function at every observation of z is estimated separately, the optimal selection of bandwidth for all observations may be not the same and the numerical problem will become untractable then the sample size is large. As a result, the researchers had to make a compromise and assume that the bandwidths used for estimating g(z) are the same.

In the nonparametric regression literature, spline-based methods, such as regression splines, smoothing splines, and penalized splines, are also popular and applied extensively in many fields. Recently, [15] proposed the multiplicative additive models based on the least product relative error criterion (LPRE), where the B-spline basis functions are used to estimate the nonparametric functions. Approximating a smooth function by a spline function has many desirable benefits, as presented in [16] and [17] . It is well known that the number and location of internal knots should be well-addressed when using the B-spline approximation approach. Too many knots may cause overfitting, few knots then bring underfitting. Although some information criteria can be adopted to select an appropriate number of knots, the overall computation burden is rather heavy. Fortunately, the penalized spline (P-splines) can avoid this difficulty and have gained remarkable popularity, especially due to computational expediency and easy adaptability to more complex models, see [18] . However, this technique has not been studied in the context of MR, and this is the first work to employ P-splines estimation for MR, although this issue is of importance and meaningful. All these motivate us to conduct this study in a formal manner.

This paper is organized as follows. In Section 2, we first introduce the nonparametric multiplicative regression model. Combining the penalized splines and the least product relative error criterion, a new estimation method is proposed, and some remarks about the selection of smoothing parameters and knots, and the asymptotic properties of the proposed estimator are given. Some simulation studies are carried out to assess the performance of our method under finite sample situations in Section 3. To illustrate its usefulness, we also apply our method to one real data set in Section 4. Finally, some discussions in Section 5 conclude the paper.

2. Methodology

In this section, we mainly introduce the model and estimation approach. At the same time, some related issues will be answered in detail.

2.1. Model and Estimation

Consider the following nonparametric multiplicative regression model

$Y=\exp \left(g\left(Z\right)\right)\epsilon \mathrm{,}$ (1)

where Y is the response variable, $Z$ is a p-vector of covariates, $\epsilon$ is the positive random error and independent of $Z$ , $g(\cdot )$ is an unknown link function. Without loss of generality and for simplicity, we will assume that $p=1$ in the following, i.e., the covariate $Z$ is univariate. But the method discussed later can be easily be extended to the general case. Assume that an i.i.d. sample $\left(Y_i\mathrm{,}{Z}_i\right)\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n$ from model (1) are observed. For model identifiability, it is required that $E\left(\epsilon \right)=E\left({\epsilon }^{-1}\right)$ in [5] . On the other hand, after taking logarithmic transformation on bosh sides of Equation (1), it follows that

$\ln \left(Y\right)=g\left(Z\right)+\ln \left(\epsilon \right)=\left(g\left(Z\right)+c\right)+\left(\ln \left(\epsilon \right)-c\right)\triangleq \tilde{g}\left(Z\right)+\ln \left(\tilde{\epsilon }\right)\mathrm{,}$ (2)

holds for any real number c, which means that the former requirement is not enough to ensure the nonparametric component unique. This phenomenon has been found by [8] , where they imposed a strong condition on the model, namely, $E\left(\ln \left(\epsilon \right)\right)=0$ . In our opinion, this is not necessary. The unknown constant c can be set to be $E\left(\ln \left(\epsilon \right)\right)$ such that $E\left(\ln \left(\tilde{\epsilon }\right)\right)=0$ holds. Thus the function $g(\cdot )$ is identifiable. Meanwhile, suppose that one obtains a initial estimator of $g(\cdot )$ , denoted by ${\tilde{g}}_n(\cdot )$ . Combining the Equation (2), a more efficient estimator of $g(\cdot )$ , denoted by ${\hat{g}}_n(\cdot )$ , can be obtained by subtracting the mean of $E\left(\ln \left(\epsilon \right)\right)$ from the initial estimator ${\tilde{g}}_n(\cdot )$ , where

${\hat{g}}_n\left(z\right)={\tilde{g}}_n(\cdot )-\bar{c},\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\bar{c}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[\ln \left(Y_i\right)-{\tilde{g}}_n\left(Z_i\right)\right].$ (3)

In this process, it only needs that $E\left(\ln \left(\epsilon \right)\right)$ , the expectation of $\ln \left(\epsilon \right)$ is finite, which is allowed to be nonzero. Therefore, our condition is weaker than that stated in [8] .

Following the LPRE technique used in [5] , the estimator of $g(\cdot )$ in model (1) can be obtained by minimizing the objective function

$S_n\left(g\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left\{\left|\frac{Y_i-\exp \left(g\left(Z_i\right)\right)}{Y_i}\right|\times \left|\frac{Y_i-\exp \left(g\left(Z_i\right)\right)}{\exp \left(g\left(Z_i\right)\right)}\right|\right\}.$ (4)

When all $Z_i$ ’s are mutually different, the computation involved in Equation (4) can be accomplished individually for each $g\left(Z_i\right)$ . However, the resulting estimator is not continuous, which is unsatisfactory in some applications. To overcome these drawbacks, we approximate $g(\cdot )$ by B-spline functions. Exactly speaking, let $-\infty <a=t_0<t_1<\cdots <t_{k_n}<t_{k_n+1}=b<\infty$ be a sequence of knots, where $\left[a\mathrm{,}b\right]$ denotes the support of Z. As discussed in [19] , a set of B-spline basis function of degree d can be constructed and denoted by ${\phi }_j(\cdot ),j=1,\cdots ,K_n=k_n+1+d$ , where $d=2$ or 3 corresponds to the quadratic spline or the cubic spline. Then $g\left(z\right)$ can be approximated by a linear combination of basis functions, i.e., $g\left(z\right)\approx {\displaystyle {\sum }_{j=1}^{K_n}}{\alpha }_j{\phi }_j\left(z\right)\triangleq {\alpha }^{{\rm T}}B\left(z\right)$ , where $\alpha ={\left({\alpha }_1\mathrm{,}\cdots \mathrm{,}{\alpha }_{K_n}\right)}^{{\rm T}}$ is the spline coefficient vector and $B\left(z\right)={\left({\phi }_1\left(z\right)\mathrm{,}\cdots \mathrm{,}{\phi }_{K_n}\left(z\right)\right)}^{{\rm T}}$ . Further, the model (1) can be rewritten as

$Y=\exp \left({\alpha }^{{\rm T}}B\left(Z\right)\right)\epsilon \mathrm{.}$ (5)

Meanwhile, $\alpha$ can be estimated by $\stackrel{⌣}{\alpha }$ , the minimizer of the loss function

$\begin{array}{c}{L}_n\left(\alpha \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left\{\left|\frac{Y_i-\exp \left({\alpha }^{{\rm T}}{B}_i\right)}{Y_i}\right|\times \left|\frac{Y_i-\exp \left({\alpha }^{{\rm T}}{B}_i\right)}{\exp \left({\alpha }^{{\rm T}}{B}_i\right)}\right|\right\}\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[Y_i\exp \left(-{\alpha }^{{\rm T}}{B}_i\right)+Y_i^{-1}\exp \left({\alpha }^{{\rm T}}{B}_i\right)-2\right],\end{array}$

where $B_i=B\left(Z_i\right)$ . Then the B-spline estimator of $g(\cdot )$ can be obtained through Equation (3), where ${\tilde{g}}_n\left(z\right)$ is replaced by ${\stackrel{⌣}{\alpha }}^{{\rm T}}B\left(z\right)$ .

However, for reasons discussed in Section 1, here we estimate the nonparametric function in model (1) using penalized splines by adding a roughness penalty to the above minimization problem. In the literature, there are mainly two kinds of definition of roughness. One is based on the integrated squared q-th derivative of the spline function, the other is based on the q-th order difference of the spline coefficient vector, which is more popular for its simplicity and adopted in this paper. Exactly, let Δ denote the backward difference operator. Then we have $\Delta {\alpha }_j={\alpha }_j-{\alpha }_{j-1}$ and ${\Delta }^2{\alpha }_j={\alpha }_j-2{\alpha }_{j-1}+{\alpha }_{j-2}$ . In General, let ${\Delta }^q=\Delta {\Delta }^{q-1}$ . With these notations, the P-spline estimator of $g(\cdot )$ is defined as $\hat{g}$ , the corrected spline function in (3), where ${\tilde{g}}_n\left(z\right)$ is replaced by ${\hat{\alpha }}^{{\rm T}}B\left(z\right)$ , and $\hat{\alpha }$ is the minimizer of

$Q_n\left(\alpha \right)=L_n\left(\alpha \right)+\lambda {\displaystyle \underset{j=q+1}{\overset{K_n}{\sum }}}{\left|{\Delta }^q{\alpha }_j\right|}^{\gamma },$

where $\gamma >0$ and $\lambda >0$ is the smoothing parameter controlling the smoothness of the fitting vurve. If $\lambda =0$ , the resulting estimator becomes the B-spline estimator. In practice, common choices of q and $\gamma$ are $q=2$ and $\gamma =2$ for computational convenience. At this time, the above minimization problem can be rewritten as

$Q_n\left(\alpha \right)=L_n\left(\alpha \right)+\lambda {\alpha }^{{\rm T}}{D}^{{\rm T}}D\alpha \mathrm{,}$ (6)

where $D$ is the $\left(K_n-2\right)\times K_n$ matrix representation of the second order difference operator ${\Delta }^2$ .

2.2. Computation

Note that for given smoothing parameter and knots, the loss function (6) is differential and convex, minimizing it is not difficult, and some common algorithms can be applied. In the later simulation studies and real data analysis, we employ the function nmk in the R package dfoptim.

However, other alternatives are still helpful. As in the least square mean regression and penalized quantile regression (e.g. [20] ), an iterative reweighed least square (IRLS) algorithm can be developed. Meanwhile, the measure for the effective model degree of freedom can be introduced based on the trace of the hat matrix, as used in [20] . In the literature, the degree of freedom usually is adopted as the measure of model complexity to facilitate further model comparison, which has not been discussed for multiplicative regression using local smoothing technique. Thus, the smoothing parameter selection can be achieved naturally using some existing information criteria or cross-validation.

When the smoothing parameter and knots are prepared, $\hat{\alpha }$ can also be seen as the root of the first derivative of objective function (6), namely,

$-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[Y_i\exp \left(-{\hat{\alpha }}^{{\rm T}}{B}_i\right)-Y_i^{-1}\exp \left({\hat{\alpha }}^{{\rm T}}{B}_i\right)\right]B_i+2\lambda D^{{\rm T}}D\hat{\alpha }=\mathrm{0,}$ (7)

which can be rewritten as

$-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}2\left(\ln Y_i-{\hat{\alpha }}^{{\rm T}}{B}_i\right)W_i\left(\hat{\alpha }\right)B_i+2\lambda D^{{\rm T}}D\hat{\alpha }=\mathrm{0,}$ (8)

where

$W_i\left(\hat{\alpha }\right)=\frac{Y_i\exp \left(-{\hat{\alpha }}^{{\rm T}}{B}_i\right)-Y_i^{-1}\exp \left({\hat{\alpha }}^{{\rm T}}{B}_i\right)}{2\left(\ln Y_i-{\hat{\alpha }}^{{\rm T}}{B}_i\right)}\mathrm{.}$ (9)

Write $\tilde{Y}={\left(\ln Y_1\mathrm{,}\cdots \mathrm{,}\ln Y_n\right)}^{{\rm T}}$ , $B=\left(B_1\mathrm{,}\cdots \mathrm{,}{B}_n\right)$ and $W\left(\hat{\alpha }\right)=diag\left(W_1\left(\hat{\alpha }\right)\mathrm{,}\cdots \mathrm{,}{W}_n\left(\hat{\alpha }\right)\right)$ . The IRLS algorithm is implemented as follows.

Step 0. For given smoothing parameter $\lambda$ , provide an initial value of $\alpha$ , denoted as ${\hat{\alpha }}^0$ , which can be taken as the unpenalized B-spline estimator, or the ordinary least square estimator based on the logarithmic transformed data. Assume that we have the k-th iterative value ${\hat{\alpha }}^k$ .

Step 1. At the $k+1$ -th step, update the parameter by

${\hat{\alpha }}^{k+1}={\left(BW\left({\hat{\alpha }}^k\right)B^{{\rm T}}+2\lambda D^{{\rm T}}D\right)}^{-1}BW\left({\hat{\alpha }}^k\right)\tilde{Y}\mathrm{.}$

Repeat Step 1 until convergence. The value at the terminated step is defined as the final estimator of $\alpha$ , also denoted by $\hat{\alpha }$ . At the same time, as in the linear regression, the hat matrix here can be defined through

$H_{\lambda }=B^{{\rm T}}{\left(BW\left(\hat{\alpha }\right)B^{{\rm T}}+2\lambda D^{{\rm T}}D\right)}^{-1}BW\left(\hat{\alpha }\right)\mathrm{.}$

2.3. Choice of Smoothing Parameter and Knots

To implement the estimator, the smoothing parameter $\lambda$ and the number and location of knots had to be determined in advance. Compared with the knots, the smoothing parameter plays a more significant role in the penalized spline estimation. As done in [20] , we set the trace of the hat matrix mentioned above to be the model degrees of freedom, denoted by $df=tr\left(H_{\lambda }\right)$ , where $tr\left(A\right)$ denotes the trace of a square matrix A. In spline context, some criterions have been proposed to select the smoothing parameter $\lambda$ , such as the Bayesian Information Criterion (BIC), Generalized Approximation Cross-Validation criterion (GACV) and Generalized Cross-Validation criterion (GCV). Especially, in this paper, they are implemented by minimizing the functions:

$\text{BIC}\left(\lambda \right)=\ln \left(L_{1n}\left(\hat{\alpha }\right)\right)+\frac{df\cdot \ln n}{2n}\mathrm{,}$

$\text{GACV}\left(\lambda \right)=\frac{L_{1n}\left(\hat{\alpha }\right)}{1-df/n}\mathrm{,}$

$\text{GCV}\left(\lambda \right)=\frac{{\displaystyle {\sum }_{i=1}^n{W}_i\left(\hat{\alpha }\right)\left(\ln \left(Y_i\right)-{\hat{\alpha }}^{{\rm T}}{B}_i\right)/n}}{{\left(1-df/n\right)}^2}\mathrm{,}$

respectively.

[18] suggested some rules for the number and location of knots. In the same spirits, we follow the strategy in [20] , where quadratic B-splines with 20 equally-spaced quantiles knots are used unless otherwise specified. With the BIC, GACV, GCV criterions above, we use 51 equally-spaced log-scaled grid points on ${\log }_{10}\left(\lambda \right)\in \left[-\mathrm{5,5}\right]$ for the choice of $\lambda$ . All these settings work well in the numerical studies and a real example analysis.

2.4. Asymptotic Results

For any probability measure P, define $L_2$ -norm ${\Vert f\Vert }_2={\left({\displaystyle \int f^2\text{d}P}\right)}^{1/2}$ . Denote $p=d+1$ . To derive the asymptotic results, some regularity conditions are needed and listed as follows.

(A1) The true function in model (1) is $g_0$ . Suppose that $g_0\in \mathcal{G}=\left\{g\in C^r\left[a\mathrm{,}b\right]\mathrm{<mo>\; :\; </mo>}{\Vert g^{\left(j\right)}\Vert }_{\infty }\le M_0\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}r\mathrm{,}\left|g^{\left(r\right)}\left(z_1\right)-g^{\left(r\right)}\left(z_2\right)\right|\le M_1\left|z_1-z_2\right|\right\}$ , where $M_0$ and $M_1$ are some positive constants, ${\Vert \cdot \Vert }_{\infty }$ is the superior norm, $q\le r\le d+1$ .

(A2) Define $h_i=t_i-t_{i-1}$ . Assume that ${\max }_i\left|h_{i+1}-h_i\right|=o\left(1/k_n\right)$ . Moreover, the ratio of maximum and minimum spacings of knots is uniformly bounded.

(A3) $k_n=o\left(n\right)$ . (A4) The covariate Z has abounded support $\left[a\mathrm{,}b\right]$ with corresponding density $f_Z\left(\mathrm{.}\right)$ , which has second derivative and is bounded away from zero and infinity.

(A5) $\epsilon$ is independent of Z, and $E\left(\epsilon \right)=E\left({\epsilon }^{-1}\right)$ , $E\left(\ln \left(\epsilon \right)\right)$ is finite.

(A6) $E\left(\epsilon +{\epsilon }^{-1}\right)<+\infty$ .

Conditions (A1)-(A3) are common requirements in the penalized spline theory. (A4) is the regularization condition used in the study of MR. (A5) is an identification condition for the LPRE estimation, which is similar to the zero mean condition in the classical linear mean regression. (A6) is required for proof, which is also used in [15] .

Theorem 2.1 Under conditions (A1)-(A6) mentioned above, we have that

1) If $r=p$ , $k_n=O\left(n^{1/\left(2p+1\right)}\right)$ and $\lambda =O\left(n^{-\beta }\right)$ with $\beta >\left(1+p-q\right)/\left(2p+1\right)$ , then it yields that ${\Vert \hat{g}-g_0\Vert }_2^2=O_P\left(n^{-2p/\left(2p+1\right)}\right)$ .

2) If $r\mathrm{<mo>\; =\; </mo>}q$ , $\lambda =O\left(n^{-1/\left(2q+1\right)}\right)$ and $k_n=O\left(n^{\beta }\right)$ with $1>\beta >1/\left(2q+1\right)$ , then it yields that ${\Vert \hat{g}-g_0\Vert }_2^2=O_P\left(n^{-2q/\left(2q+1\right)}\right)$ .

Proof. Note that

$\begin{array}{l}{Q}_n\left(\alpha \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[Y_i\exp \left(-{\alpha }^{{\rm T}}{B}_i\right)+Y_i^{-1}\exp \left({\alpha }^{{\rm T}}{B}_i\right)-2\right]+\lambda {\alpha }^{{\rm T}}{D}^{{\rm T}}D\alpha \\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[{\epsilon }_i\exp \left(R_i+\left(g^{\ast }\left(Z_i\right)-{\alpha }^{{\rm T}}{B}_i\right)\right)+{\epsilon }_i^{-1}\exp \left(-R_i-\left(g^{\ast }\left(Z_i\right)-{\alpha }^{{\rm T}}{B}_i\right)\right)-2\right]\\ +\lambda {\alpha }^{{\rm T}}{D}^{{\rm T}}D\alpha \mathrm{,}\end{array}$

where ${\epsilon }_i=Y_i/\exp \left(g_0\left(Z_i\right)\right)$ , $R_i=g\left(Z_i\right)-g^{\ast }\left(Z_i\right)$ , and $g^{\ast }\left(Z\right)={\alpha }^{\ast }B\left(Z\right)$ denotes the best spline function approximation. Define $X=n^{-1}{\displaystyle {\sum }_{i=1}^n{B}_i{B}_i^{{\rm T}}}$ , $G_{\lambda }=X+\lambda D^{{\rm T}}D$ and $\theta =G_{\lambda }^{1/2}\left({\alpha }^{\ast }-\alpha \right)$ . Then we have that $g^{\ast }\left(Z_i\right)-{\alpha }^{{\rm T}}{B}_i=B_i^{{\rm T}}\left({\alpha }^{\ast }-\alpha \right)=B_i^{{\rm T}}{G}_{\lambda }^{-1/2}\theta$ . Let $a_n$ denote the convergence rate in the Theory. Since $\hat{\theta }\mathrm{<mo>\; =\; </mo>}{a}_n\left(\hat{\alpha }-{\alpha }^{\ast }\right)=\arg {\min }_{\theta }\left[Q_n\left({\alpha }^{\ast }+a_n^{-1}\theta \right)-Q_n\left({\alpha }^{\ast }\right)\right]$ , we only need to prove that for every $\nu >0$ , there exists a large constant $\delta >$ such that

$P\left(\underset{{\Vert \theta \Vert }_2=\delta }{\inf }{Q}_n\left({\alpha }^{\ast }+a_n^{-1}\theta \right)-Q_n\left({\alpha }^{\ast }\right)\right)\ge 1-\nu \mathrm{.}$

In the line of the proof of Theorem 1 in the Appendix of [17] , some similar techniques can be applied. Finally, by the Corollary 1 of [16] , the desirable results are derived. The proof is completed.

3. Numerical Studies

In this section, numerical studies were conducted to evaluate the finite sample performance of our proposed method under various situations. To fairly compare the un-penalized B-spline estimator with the penalized spline estimators obtained by the BIC, GCV, GACV criterions and for simplicity, we set the degree of spline basis to be $q=2$ , the number of internal knots $k_n=20$ , which are located on equally-spaced quantiles for all methods. All results below are based on 500 replicates, where the sample size $n=100,300,500$ , respectively. All simulations are implemented using the software R.

3.1. Model A and Results

Model A. We generated $\left(Y\mathrm{,}Z\right)$ from the following model

$Y=\exp \left(1+2\sin \left(\frac{\pi Z}{30}\right)\right)\epsilon \mathrm{,}$ (10)

where $Z~Unif\left(\mathrm{0,50}\right)$ and is independent with the random error $\epsilon$ . Further, three error distributions are considered, namely,

Case (i): $\epsilon =\exp \left(U\right)$ , $U~N\left(\mathrm{0,1}\right)$ ,

Case (ii): $\epsilon =\exp \left(U\right)$ , $U~Unif\left(-\mathrm{2,2}\right)$ ,

Case (iii): $\epsilon ~Unif\left(\mathrm{0.0001,4.635506}\right)$ .

Note that under case (i)-(ii), it holds that $E\left(\epsilon -{\epsilon }^{-1}\right)=0$ and $E\left(\log \left(\epsilon \right)\right)=0$ . But under case (iii), $E\left(\epsilon -{\epsilon }^{-1}\right)=0$ , $E\left(\log \left(\epsilon \right)\right)=0.5339771$ . We use the averaged Integrated Absolute Bias (IABIAS) and Mean Integrated Square Error (MISE), where for one estimator ${\hat{g}}_j\left(j=1,\cdots ,500\right)$ obtained from the j-th sample,

$\text{IABIAS}=\frac{1}{500}{\displaystyle \underset{j=1}{\overset{500}{\sum }}}\left[\frac{1}{ngrid}{\displaystyle \underset{k=1}{\overset{ngrid}{\sum }}}\left|{\hat{g}}_j\left(u_k\right)-g_0\left(u_k\right)\right|\right]\mathrm{,}$

$\text{MISE}=\frac{1}{500}{\displaystyle \underset{j=1}{\overset{500}{\sum }}}\left[\frac{1}{ngrid}{\displaystyle \underset{k=1}{\overset{ngrid}{\sum }}}{\left[{\hat{g}}_j\left(u_k\right)-g_0\left(u_k\right)\right]}^2\right]\mathrm{,}$

at the fixed grid points $\left\{u_k\right\}$ equally-spaced in [0, 50] and $ngrid=501$ . The values below them are the associated sample standard deviation. To study the different smoothing parameter selection methods at a finer scale, we also calculate the mean and standard deviation of effective model degree of freedom (DF).

Table 1 reported the results. It can be seen that all three penalized estimators

Table 1. Results of Model A under cases (i)-(iii) with different sample sizes and criterions (×10⁻²).

have smaller IABIAS and MISE than those of the unpenalized B-spline estimator. This indicates the benefits of the proposed method. As the sample size increases, the differences between them decrease and become comparable. The model degree using GACV are almost larger than those of using BIC and GCV in all settings, and those using BIC are the smallest, which is also observed in [20] . These imply that GACV tends to select more complex models, BIC tends to select simpler models, GCV makes a compromise and lies between them. Figure 1 displays the average estimated curves and their 95% point-wise confidence band under cases (i) with sample size $n=300$ for model A. We can see that all estimated curves are close to the true one, but the B-spline estimation has a drastic deviation near the boundaries. Figure 2 shows the boxplots of the estimates of $g\left(z\right)$ at $z=\mathrm{15,30,45}$ , respectively. We can see that most estimates are centered around the true values of $\left(\mathrm{3,1,}-1\right)$ with small deviations, although the BS estimator has a little larger standard deviation. The corresponding QQ-plots of them are presented in Figure 3, which show that the resulting estimators seem to converge asymptotically to the normal distribution. Estimated curves under other situations are similar, and are not displayed, but are available in the supplementary materials.

Figure 1. Estimate of $g\left(Z\right)$ for Model A under case (i) with sample size 300. The solid thick line (black) and solid thin line (red) correspond to the true and estimated curves, respectively. The dashed lines (blue) are the 95% point-wise confidence bands.

Figure 2. Box-plots of $g\left(Z\right)$ at $Z=\mathrm{15,30,45}$ , respectively, for Model A under case (i) with sample size 300.

3.2. Model B and Results

We generated $\left(Y\mathrm{,}Z\right)$ from the following model

$Y=\exp \left(\frac{{\left|Z-25\right|}^3}{5000}\right)\epsilon \mathrm{,}$ (11)

where $Z~Unif\left(\mathrm{0,50}\right)$ and is independent with the random error $\epsilon$ . Others are the same as defined in Model A. We reported the results in Table 2 and displayed the average estimated curves and their 95% point-wise confidence band under cases (i) with sample size $n=300$ for model B in Figure 4. Again, results are similar to those presented in Table 1 and Figures 1-3 are found. More related figures can be available in the supplementary materials.

4. Real Data Analysis

In this section, we will analyze the ethanol data to illustrate our proposed method. This data is available in the R package lattice, and aims to investigate the relationship between the emissions of nitrogen oxides, denoted by Y, and various settings of the engine compression ratio and equivalence ratio, denoted by Z, when ethanol fuel is burned in a single-cylinder engine. In this data set, 88 observations are collected. [13] analyzed it in the framework of varying coefficient multiplicative regression model based on the LPRE criterion and local linear approximation technique. His results demonstrate that various settings of the engine compression ratio have little effect in Y, as illustrated in page 280 therein. Therefore, here we only consider the covariate Z in the following model.

Figure 3. QQ-plots of $g\left(Z\right)$ at $Z=\mathrm{15,30,45}$ , respectively, for Model A under case (i) with sample size 300.

Table 2. Results of Model B under cases (i)-(iii) with different sample sizes and criterions (×10⁻²).

Figure 4. Estimate of $g\left(Z\right)$ for Model B under case (i) with sample size 300. The solid thick line (black) and solid thin line (red) correspond to the true and estimated curves, respectively. The dashed lines (blue) are the 95% point-wise confidence bands.

$Y=\exp \left(g\left(Z\right)\right)\epsilon \mathrm{.}$

The estimated curves of g with 8 and 20 internal knots quadratic P-splines are plotted in Figure 5, which have similar trends to the plot in [13] except a vertical shift. We can see that all the four curves are almost identical, which maybe caused by the data’s own well fitness to the above model. Scatted plot of Y and the corresponding estimated curves of $\exp \left(g\left(Z\right)\right)$ are presented in Figure 6, which indicate the proposed method fits the data well enough. The number of knots has some effects on the estimates of curves, but doesn’t affect the performance substantially, which is also supported by the mean and median of absolute prediction errors $\left|Y_i-{\hat{Y}}_i\right|$ and squared prediction errors ${\left(Y_i-{\hat{Y}}_i\right)}^2$ given in the Table 3, where ${\hat{Y}}_i=\exp \left(\hat{g}\left(Z_i\right)\right)$ .

Figure 5. Estimated curves of $g\left(Z\right)$ for the ethanol data. The solid (red) and dashed (blue) curves correspond to the estimated function curves with 8 and 20 internal knots, respectively.

Table 3. Prediction errors and degrees of freedom of the ethanol data.

Note: k_n: the nunber of interior knots; MAPE: mean absolute prediction errors; MDAPE: median absolute prediction errors; MSPE: mean squared prediction errors; MASPE: median squared prediction errors.

Figure 6. Scatted plot of Y and estimated curves of $\exp \left(g\left(Z\right)\right)$ for the ethanol data. The solid (red) and dashed (blue) curves correspond to the estimated function curves with 8 and 20 internal knots, respectively.

5. Conclusion and Discussions

In this study, we used the penalized spline method to estimate the nonparametric function in model (1) based on the LPRE loss function. Inspired by the iteratively reweighted least squares algorithm, we develop the effective model degree of freedom and propose three smoothing parameter selection methods. Some asymptotic results are established. Furthermore, numerical simulation studies and real data analysis found that the proposed approach works well in several settings. As indicated in Section 1, the approaches proposed in this paper may be adapted to the partially linear, even additive, single index or varying coefficients multiplicative regression models. Our future work will also consider extensions to fields with covariates measurement errors, censored data or longitudinal data analysis, which are meaningful for practitioners.

Data Availability Statement

The data set used in the real data analysis can be available publicly from the R package lattice.

Acknowledgements

This work was partly supported by a grant from Natural Science Foundation of Jiangsu Province under project ID BK20210889 and Start-up fund for doctoral research of Jiangsu University of Science and technology. The authors also thank the lecturer Feng-Ling Ren, School of Computer and Engineering, Xinjiang University of Finance & Economics, for helpful discussions.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Chen, K., Guo, S., Lin, Y. and Ying, Z. (2010) Least Absolute Relative Error Estimation. Journal of the American Statistical Association, 105, 1104-1112. https://doi.org/10.1198/jasa.2010.tm09307
[2]	Khoshgoftaar, T.M., Bhattacharyya, B.B. and Richardson, G.D. (1992) Predicting Software Errors, during Development, Using Nonlinear Regression Models: A Comparative Study. IEEE Transactions on Reliability, 41, 390-395. https://doi.org/10.1109/24.159804
[3]	Narula, S.C. and Wellington, J.F. (1977) Prediction, Linear Regression and the Minimum Sum of Relative Errors. Technometrics, 19, 185-190. https://doi.org/10.1080/00401706.1977.10489526
[4]	Park, H. and Stefanski, L.A. (1998) Relative-Error Prediction. Statistics and Probability Letters, 40, 227-236. https://doi.org/10.1016/S0167-7152(98)00088-1
[5]	Chen, K., Lin, Y., Wang, Z. and Ying, Z. (2016) Least Product Relative Error Estimation. Journal of Multivariate Analysis, 144, 91-98. https://doi.org/10.1016/j.jmva.2015.10.017
[6]	Zhang, Q. and Wang, Q. (2013) Local Least Absolute Relative Error Estimating Approach for Partially Linear Multiplicative Model. Statistica Sinica, 23, 1091-1116. https://doi.org/10.5705/ss.2012.133
[7]	Hirose, K. and Masuda, H. (2018) Robust Relative Error Estimation. Entropy, 20, Article 632. https://doi.org/10.3390/e20090632
[8]	Zhang, J., Feng, Z. and Peng, H. (2018) Estimation and Hypothesis Test for Partial Linear Multiplicative Models. Computational Statistics & Data Analysis, 128, 87-103. https://doi.org/10.1016/j.csda.2018.06.017
[9]	Chen, Y. and Liu, H. (2021) A New Relative Error Estimation for Partially Linear Multiplicative Model. Communications in Statistics-Simulation and Computation, 52, 4962-4980.
[10]	Liu, H. and Xia, X. (2018) Estimation and Empirical Likelihood for Single-Index Multiplicative Models. Journal of Statistical Planning and Inference, 193, 70-88. https://doi.org/10.1016/j.jspi.2017.08.003
[11]	Zhang, J., Zhu, J. and Feng, Z. (2019) Estimation and Hypothesis Test for Single-Index Multiplicative Models. Test, 28, 242-268. https://doi.org/10.1007/s11749-018-0586-2
[12]	Zhang, J., Cui, X. and Peng, H. (2020) Estimation and Hypothesis Test for Partial Linear Single-Index Multiplicative Models. Annals of the Institute of Statistical Mathematics, 72, 699-740. https://doi.org/10.1007/s10463-019-00706-6
[13]	Hu, D.H. (2019) Local Least Product Relative Error Estimation for Varying Coefficient Multiplicative Regression Model. Acta Mathematicae Applicatae Sinica, English Series, 35, 274-286. https://doi.org/10.1007/s10255-018-0794-2
[14]	Chen, Y., Liu, H. and Ma, J. (2022) Local Least Product Relative Error Estimation for Single-Index Varying-Coefficient Multiplicative Model with Positive Responses. Journal of Computational and Applied Mathematics, 415, Article ID: 114478. https://doi.org/10.1016/j.cam.2022.114478
[15]	Ming, H., Liu, H. and Yang, H. (2022) Least Product Relative Error Estimation for Identification in Multiplicative Additive Models. Journal of Computational and Applied Mathematics, 404, Article ID: 113886. https://doi.org/10.1016/j.cam.2021.113886
[16]	Kalogridis, I. and Van Aelst, S. (2021) M-Type Penalized Splines with Auxiliary Scale Estimation. Journal of Statistical Planning and Inference, 212, 97-113. https://doi.org/10.1016/j.jspi.2020.09.004
[17]	Kalogridis, I. and Van Aelst, S. (2021) Robust Penalized Spline Estimation with Difference Penalties. Econometrics and Statistics, 29, 169-188.
[18]	Ruppert, D., Wand, M.P. and Carroll, R.J. (2003) Semiparametric Regression. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511755453
[19]	Schumaker, L. (2007) Spline Functions: Basic Theory. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511618994
[20]	Wu, C. and Yu, Y. (2014) Partially Linear Modeling of Conditional Quantiles Using Penalized Splines. Computational Statistics & Data Analysis, 77, 170-187. https://doi.org/10.1016/j.csda.2014.02.020