Multivariate calibration with robust signal regression

Although much research has been done on robust estimation in the context of smoothing, as far as we know, little or no work has been done in implementing robust smoothing into the multivariate calibration problem, for example, using robust variants of PSR. An assortment of loss functions—in addition to squared error loss—have been applied to the penalized spline and regression splines. For instance, Huber (1979) and Cox (1983) proposed to use M-type robust smoothing with cubic regression splines. Other articles on robust smoothing and nonparametric regression include Hardle and Gasser (1984), Silverman (1985) and Hall and Jones (1990). Eilers and de Menezes (2005) applied quantile smoothing with L 1 loss on comparative genomic hybridization (CGH) data. Robust penalized regression spline was explored by Lee and Oh (2007), in which Huber loss was used for bivariate smoothing. In Section 6.3, we fully draw the connection between our work and that of Lee and Oh (2007). They also applied the robust penalized regression spline to an additive mixed modelling setting. Tharmaratnam et al. (2010) studied S-estimation for penalized regression splines. A recent paper by Kneib (2013) discussed the semiparametric regression models that go beyond the mean and squared error loss, such as quantile and expectile regressions.

4.1 Generalized Huber loss

The least-squares estimate works well in the presence of normal error or even in situations where the error distribution is reasonably continuous, symmetric and not too heavy tailed. However, with departures from the above, the least-squares criterion can lead to poor results, especially for heavy-tailed errors and outliers in the response. A ‘robust’ statistical procedure is often able to provide useful information even when some of the assumptions are not applicable. In regression analysis, a convenient approach is simply to remove the influential observations from the least-squares fit, whereas alternatively ‘robust regression’ employs a criterion that is more resistant (to unusual observations) than those found using least squares. A tried-and-true method of robust regression is ‘ M -estimation’, introduced by Huber (1964). The general M -estimator minimizes the objective function

i = 1 n ρ ( e i ) = i = 1 n ρ y i − β 0 − j = 1 p x j β j ,

(4.1)

where the function ρ has the following properties: (a) ρ ( e ) ≥ 0 , (b) ρ ( 0 ) = 0 , (c) ρ ( e ) = ρ ( − e ) and (d) ρ ( e i ) ≥ ρ ( e j ) for | e i | > | e j | . For example, for least-squares estimation uses ρ ( e ) = e 2 . (Huber 1981) further described a well-known robust M -estimator employing a loss function that is less affected by very large residual values. The Huber criterion can be written as

ρ H ( e ) = e 2 | e | < K 2 K | e | − e 2 = K 2 + 2 K ( | e | − K ) | e | ≥ K .

(4.2)

The parameter K defines the point of transition from quadratic to linear loss. Specifically, errors smaller than K get squared, while larger errors only increase the loss linearly.

Li and Yu (2009) generalized the Huber loss criterion described in (4.2) to a class of M -estimators, called ‘generalized Huber criterion’ as follows

ρ η ( e ) = e 2 | e | < K K 2 + 2 η K ( | e | − K ) | e | ≥ K ,

(4.3)

where 0 ≤ η ≤ 1 . The left panel of Figure 4 illustrates the family of generalized Huber loss with three different values of η at K = 2 . The black solid curve is the square loss function. The grey-shaded region represents the region for all the possible generalized Huber loss ρ η ( e ) with η ranging from zero to one. The right panel of the same figure shows the difference between squared error loss and ρ η ( e ) .

OPEN IN VIEWER Figure 4:

Illustration of generalized Huber loss with three different values of η at K = 2 (left) and the corresponding trailing convex functions (right). The three generalized Huber loss functions can be expressed as the differences of the squared error loss (solid curve, left) and the corresponding trailing functions (right)

The class of generalized Huber loss falls into the category of ‘redescending M -estimator’, which includes: Tukeys bi-weight, S -estimators (Maronna et al., 2006) and t -type scores (He et al., 2000). Clearly, the smaller the η , then the more robust the loss function is against the outliers in error. When η is one, the generalized Huber loss ρ η = 1 is the Huber loss ρ H . The truncated least-squares criterion corresponds to ρ η = 0 , which approximates the Tukey's bi-weight M -estimator. Although the generalized Huber criterion is not convex (in e ∈ ℝ ) for 0 ≤ η < 1 , it can be expressed as a difference of two convex functions (of e ) as follows:

ρ η ( e ) = e 2 − I ( | e | 〉 k ) [ 〉 e 2 + 2 η K ( K - | e | ) − k 2 ] ,

(4.4)

where I ( · ) is an indicator function, equal to one when condition holds, zero otherwise. Note that the leading convex function is the square loss function. In Figure 4, the right panel further shows the second term on the right side of (4.4), which is convex and K -insensitive (i.e., it is a constant within the range [ − K , K ] ).

The difference convex (d.c.) programming, developed by An and Tao (1997), addresses the problem of minimizing an objective function, which can be expressed as a difference of two convex functions, on the whole space. Consider minimizing an objective function h ( a ) which is a difference of two convex functions, that is, h ( a ) = h 1 ( a ) − h 2 ( a ) where both h 1 ( a ) and h 2 ( a ) are convex in a . The key idea of d.c. programming is to construct a sequence of subproblems, which are obtained by replacing the trailing convex function, for example, h 2 ( a ) , by its first order approximation function h 2 ( a ( o ) ) + ⟨ a − a ( o ) , ▽ h 2 ( a ( o ) ) ⟩ and solve them iteratively, where ▽ h 2 ( a ( o ) ) is the subgradient of h 2 ( a ) at a ( o ) with respect to a . Specifically, given the current solution for the subproblem a cur , the new subproblem solves

a new = arg min a h 1 ( a ) − h 2 ( a cur ) + ⟨ a − a cur , ▽ h 2 ( a cur ) ⟩ .

(4.5)

Note that after removing the constant terms in (4.4), minimizing the subproblem is equivalent to

a new = arg min a h 1 ( a ) − ⟨ a , ▽ h 2 ( a cur ) ⟩ .

(4.6)

In this article, we suggest to replace the least-square criterion by the generalized Huber criterion described in (4.3) within the PSR framework. Hence the robust PSR (rPSR) minimizes

Q ( α ) = i = 1 m ρ η ( y i − U i ′ α ) + λ α ′ D d ′ D d α ,

(4.7)

which can be represented as a difference of two convex functions as follows:

Q ( α ) = h 1 ( α ) − h 2 ( α ) , W H E R E

(4.8)

h 1 ( α ) = ∑ i = 1 m e i 2 + λ α ′ D ′ D α

(4.9)

h 2 ( α ) = ∑ i = 1 m I ( | e i | > K ) [ e i 2 + 2 η K ( K - | e i | ) − K 2 ] ,

(4.10)

and e i = y i − U i ′ α . The subgradient of h 2 ( α ) with respect to α is

▽ h 2 ( α ) = ∂ h 2 ∂ e · ∂ e ∂ α = − 2 i = 1 m I ( | e i | > K ) [ e i − η K Sign ( e i ) ] U i ,

(4.11)

where Sign ( a ) is the sign function, equal to 1 if a is positive, − 1 if a is negative and 0 otherwise. The vector U i ′ is the i th row of U = XB , with n elements { u ij } j = 1 n . Let V be a column vector of length m with elements { I ( | e i | > K ) [ e i − η K Sign ( e i ) ] } i = 1 m . It follows that the right side of (4.7) can be expressed as − 2 U ′ V . The inner product of α and subgradient ▽ h 2 ( α ) is then

⟨ α , ▽ h 2 ( α ) ⟩ = − 2 α ′ U ′ V .

(4.12)

Through d.c. programming, the minimization of the objective function (4.6) translates to the minimizing of a sequence of subproblems

α ̂ = arg min α ( Y − U α ) ′ ( Y − U α ) + λ α ′ D ′ D α + 2 α ′ U ′ V .

(4.13)

Setting the first order derivative of the right side of (4.9) to zero, we have the closed form of the solution

α ̂ = U ′ U + λ D ′ D − 1 U ′ ( Y − V ) = U ′ U + λ D ′ D − 1 U ′ Y A .

(4.14)

The right side of (4.10) further shows that the subproblem solution is itself a modified PSR solution, one with the adjusted responses Y A defined as

Y A = y 1 − I ( | e 1 | > K ) [ e 1 − η K Sign ( e 1 ) ] ⋮ y m − I ( | e m | > K ) [ e m − η K Sign ( e m ) ] m × 1 .

(4.15)

Note that only the observations with the residuals greater than K (in absolute value) will be ‘adjusted’. Further, if K is greater than all the residuals { e i }

OPEN IN VIEWER

Robust PSR Algorithm

1. Initializations: Choose the tuning parameter value λ and η . Construct B using a rich set of n B-spline basis functions of degree q on equally spaced knots and penalty order d . Default q = d = 3 . Calculate U = XB Calculate α ̂ = PSR ( U , Y , λ , d , n , q ) . 2. Cycle until convergence of α ̂ : Calculate residuals { e i } i = 1 m . Find the K based on residuals. Update the adjusted response vector Y A according to η and K . Update α ̂ = PSR ( U , Y A , λ , d , n , q ) . 3. Prediction: y ̂ new = x new ′ B α ̂

From the above algorithm, we have the following remarks: (a) For initial α ̂ , we use the PSR estimate (with the same value of λ ). Note that we deliberately fix the tuning parameter λ during the iterations to minimize the objective function (4.6); (b) Let α ̂ j pre and α ̂ j cur be the j th element of the α vector for the previous and current iteration. The algorithm terminates when max { | ( α ̂ j cur − α ̂ j pre ) / α ̂ j pre | } j = 1 n < ε , where ε is a pre-specified convergence tolerance. We set ε at 10 − 6 here; (c) The cut-off value K in the generalized Huber criterion is chosen based on the proportion, say γ , of the outliers among the residuals. In our algorithm, the 1.5 × ‘IQR rule’ (interquartile range) is used to identify the outlying residuals. Specifically, a residual is regarded as an outlier if it is either 1.5 × ‘IQR’ above the third quartile ( Q 3 ) or 1.5 × ‘IQR’ below the first quartile ( Q 1 ). The 1.5 × ‘IQR rule’ is commonly used for outlier detection (see, e.g., Rousseeuw and Hubert, 2011). The cut-off value K is set to be the γ -quantile of the distribution { | e i | } i = 1 m within each iteration. Other choices for the cut-off value K are available. For example, Lee and Oh (2007) and Tharmaratnam et al. (2010) chose the cut-off value for Huber loss based on the median absolute deviation (MAD) of the residuals. Although MAD is a popular robust estimator of scale and has the highest possible breakdown point (50%, twice as much as IQR), it achieves low efficiency at normal and other distributions (see, e.g., Randal, 2008). Furthermore, MAD assumes the symmetry on the distribution by taking equal importance to the positive and negative deviations from the median (Rousseeuw and Croux, 1993). Hence, MAD does not seem to be a good choice for asymmetric error distributions. On the other hand, IQR does not have this problem, since the quartiles need not be equally far away from the centre; (d) As with λ , the optimal value for η can be tuned through a grid search based on cross-validation performance. In this study, we only consider the value of η at 0, 0.5 and 1, unless specified otherwise; (e) The above algorithm usually converges within a few iterations. Note that updating α ̂ in Step 2 is not a refitting of PSR but conveniently only a matrix–vector multiplication. This follows since U ′ U + λ D ′ D − 1 U ′ in (4.10) does not change within Step 2. Hence, the above algorithm is computationally efficient.

5.1 Simulation studies

The purpose of this simulation study is twofold. First, we aim to show that the proposed rPSR is competitive with PSR for the normal errors. Second, when large errors exist (e.g., the error term has a heavy-tailed distribution), we aim to show that rPSR achieves better performance than PSR in terms of both prediction accuracy as well as model stability.

In this simulation study, we use the VisNIR spectra of the soil data described in Section 2 as the input variables. To generate responses, the PSR model with λ = 10 − 5 was fitted using the CEC variable as the response on the entire dataset. Note that 10 − 5 is a typical value of λ in PSR model on CEC for the soil dataset. The predicted values { f ̂ i } i = 1 675 from the PSR model are used as the ‘true’ values for this simulation study. Note the λ value is chosen based on the tenfold cross-validation on several random splits of the data. The data are then randomly split into a training set (506 observations or approximately 75% of the total sample size) and into a test set (169 observations). For the training samples, we created some ‘artificial responses’ y i * by adding random errors e i to the ‘true’ values f ̂ i (i.e., y i * = f ̂ i + e i ). Three types of error distributions are considered in this study: a normal distribution (i.e., e i ∼ N ( 0 , 2 . 39 2 ) , a mixed normal distribution and a slash distribution. Note that 2.39 is the standard deviation (SD) of the residuals from the PSR model, which is used to create the ‘true’ value f ̂ . The mixed normal errors are generated from 0.95 N ( 0 , 2 . 39 2 ) + 0.05 N ( 0 , 23 . 9 2 ) , that is, the error constituted with 95% from N ( 0 , 2 . 39 2 ) and 5% from N ( 0 , 23 . 9 2 ) . The slash distribution is defined as a standard normal variate divided by an independent standard uniform variate (i.e., N ( 0 , 1 ) / U ( 0 , 1 ) ). The slash distribution is well known for its heavy tail and extreme outliers.

We followed the ‘general recipe’ for PSR given in Marx and Eilers (1999) to choose the design parameters. Namely, the cubic B-splines (default value for q ) and the third-order difference penalty (default value for d ) were used on 100 equally spaced knots, providing sufficient (ample) flexibility. The optimal values for λ for rPSR and PSR were each found by minimizing the tenfold cross-validation error on the training set. For rPSR, we considered three levels on η : 0, 0.5 and 1. The simulation results are based on 50 random splits of the dataset. The prediction results are evaluated using ‘root mean square error’ (RMSE) and ‘mean absolute error’ (MAE) on the test set. The RMSE and MAE are defined as follows:

RMSE = 1 m test i = 1 m test ( f ̂ i − y ̂ i ) 2 ,

(5.1)

MAE = 1 m test i = 1 m test | f ̂ i − y ̂ i | ,

(5.2)

where m test is the number of observations on the test set and y ̂ i is the predicted response for the i th subject in the external test set, using the parameter estimates from the training set with the optimal λ . To compare the prediction performance, we use the ‘comparative’ test errors, defined by

c i , j = d ij min { d i , l } l = 1 , … , 4 , i = 1 , … , 50 , j = 1 , 2 , 3 , 4

where d ij is a performance measure (e.g., RMSE and MAE) over 50 replications for each of the four methods: PSR and rPSR with η = 1 , 0.5 and 0 . This quantity facilitates individual comparisons by using the test error of the best method for each dataset to calibrate the difficulty of the problem. Figure 5 shows the boxplots of the comparative RMSEs and MAE for PSR and rPSR based on 50 random replications. Table 1 shows the average of test RMSEs and MAEs based on 50 replications.

OPEN IN VIEWER Figure 5:

Boxplots of comparative RMSE for rPSR and PSR under three error distributions

OPEN IN VIEWER

Table 1:

Average of test RMSEs and MAEs in simulation study

RMSE	PSR	rPSR ( η = 1 )	rPSR ( η = 0.5 )	rPSR ( η = 0 )
Normal	0.696	0.694	0.699	0.709
Mixed	1.422	0.897	0.820	0.800
Slash	13.569	1.728	1.452	1.310
MAE	PSR	rPSR ( η = 1 )	rPSR ( η = 0.5 )	rPSR ( η = 0 )
Normal	0.440	0.437	0.442	0.446
Mixed	0.885	0.545	0.508	0.502
Slash	7.646	1.022	0.850	0.792

Based on Figure 5 and Table 1, we have the following remarks: (a) For the normal error case, all three rPSR models have very close prediction performance with their competitor PSR. In fact, rPSR with Huber loss (i.e., η = 1 ) performs slightly better than PSR; (b) For the mixed case, all three rPSR methods perform substantially better than PSR. In particular, the rPSR models with η = 0 and 1 perform better than rPSR with Huber loss; (c) For the slash distribution case, PSR performs much worse than all three rPSR models. In some extreme case, PSR performs about 100 times worse than rPSR. Among the three rPSR models, rPSR with truncated least-square loss ( η = 0 ) performs the best.

Model stability is an important issue in regression analysis. Generally, it is desired to have a stable estimate of regression coefficients and prediction on new data points. We consider the same three error distributions as above. For each error distribution, we randomly sample 95% of the data points (i.e., 641 out of 675), upon which the PSR and rPSR models (i.e., η = 0 , 0.5 and 1 ) with λ = 10 − 5 are fitted. To compare the model stability on regression coefficients, we use the L 2 distance standard deviation ( L 2 DSD ) criterion, defined by

L 2 DSD = StDev ( { | | β ̂ ( i ) − β ¯ | | 2 } i = 1 20 ) ,

where β ¯ is the average regression coefficient vector based on 20 random replications. Note that β is the regression coefficient vector on the wavelength, not on the basis functions, and | | β ̂ ( i ) − β ¯ | | is the Euclidean distance between the i th estimate of β and the average of the estimate over 20 trials. Table 2 shows the values of L 2 DSD for each method under different scenarios. We see that under the normal case, the variation of the regression coefficients on four methods is very close to each other. But PSR has much higher variation on regression coefficients in mixed and slash cases, where error distributions are heavily tailed.

OPEN IN VIEWER

Table 2:

Summary of L 2 DSD in four methods and three error distributions

Error distribution	PSR	rPSR ( η = 1 )	rPSR ( η = 0.5 )	rPSR ( η = 0 )
Normal	202	200	191	186
Mixed	1 247	402	359	377
Slash	2 495	442	416	512

To evaluate the model stability on prediction, we use each of the 20 models (for each method and error distribution) to predict all 675 observations in the dataset. The SD of the predicted values on each observation was then calculated. The average of these 675 SDs are used to evaluate the model's consistency on prediction under data perturbation (i.e., resample 95% of the data). Table 3 shows the average SD for each method under different scenarios. We see that under normal errors, all four methods have close variation on prediction. However, in mixed normal and slash cases, PSR has much higher variation on prediction than the rPSR method.

OPEN IN VIEWER

Table 3:

Summary of average SD in four methods and three error distributions

Error distribution	PSR	rPSR ( η = 1 )	rPSR ( η = 0.5 )	rPSR ( η = 0 )
Normal	0.126	0.126	0.135	0.148
Mixed	0.361	0.150	0.148	0.152
Slash	0.728	0.275	0.224	0.237

As with the simulation study, we randomly split the dataset into a training (sample size is 506) and a test set (the remaining 169 observations). PSR and rPSR with η = 1 , 0.5 and 0 are fitted on the training set and predict the test samples for all nine soil variables. Figure 6 shows the boxplots of comparative RMSE and MAE between rPSR (with η = 1 and 0 ) and PSR on nine soil variables. Table 4 shows the average test RMSE and MAE of PSR and rPSR on nine soil variables, and highlights the best result in each column. The results are based on 50 random splits of the dataset. From Figure 6 and Table 4, we see that rPSR achieves better overall performance than PSR. In general, rPSR with Huber loss performs the best in RMSE and rPSR with η = 0.5 performs the best in MAE.

OPEN IN VIEWER

Table 4:

Average of test RMSEs and MAEs on nine soil variables

RMSE	CEC	EC	Nitrogen	Carbon	SOM	Clay	Sand	Silt	pH
PSR	2.622	290.2	0.01848	0.1817	0.3383	3.240	5.425	4.473	0.3740
rPSR ( η = 1 )	2.590	286.7	0.01815	0.1794	0.3284	3.117	5.362	4.412	0.3729
rPSR ( η = 0.5 )	2.595	287.9	0.01806	0.1782	0.3269	3.138	5.397	4.413	0.3731
rPSR ( η = 0 )	2.624	290.1	0.01818	0.1795	0.3288	3.204	5.479	4.438	0.3782
MAE	CEC	EC	Nitrogen	Carbon	SOM	Clay	Sand	Silt	pH
PSR	1.795	211.0	0.01265	0.1313	0.2065	2.233	4.011	3.325	0.2842
rPSR ( η = 1 )	1.749	206.5	0.01240	0.1300	0.1946	2.105	3.940	3.291	0.2828
rPSR ( η = 0.5 )	1.737	205.5	0.01227	0.1287	0.1924	2.101	3.946	3.289	0.2821
rPSR ( η = 0 )	1.747	205.8	0.01232	0.1290	0.1933	2.124	3.977	3.302	0.2855

OPEN IN VIEWER Figure 6:

Boxplots of comparative RMSE and MAE for rPSR and PSR on nine responses

In order to see the identified outlying samples in the rPSR model and compare the coefficients between PSR and rPSR model, the PSR and rPSR model with η = 0.5 are fitted on all 675 spectra samples using carbon as the response variable. Since we already have a profile of optimal λ values using CV based on 50 replications in the previous study, the median is used for rPSR and PSR models to fit all the samples. Figure 7 shows the PCA plot of all 675 spectra samples on the leading two principal components, which explain 79.8% of the total variance. Different symbols are used to distinguish the three sample locations: Nebraska, California and Texas. The highlighted symbols are the samples identified as outliers in the rPSR model. The rPSR model identified 17 outliers, which is about 2.5% of the total sample size. We see that the majority of the outlying samples are from Nebraska and Texas. We did not see any obvious pattern for the outliers in the PCA plot. The middle panel shows the prediction plot for the rPSR model. We see that the outliers are evenly distributed above and below the 45 degree line. Interestingly, although there are two samples on the right end of the plot, only one of them is an outlier. The right panel shows the coefficient curves for the rPSR and PSR models. The coefficients of two models are close to each other. Overall, the coefficient curve of rPSR model is slightly smoother than PSR model.

OPEN IN VIEWER Figure 7:

PCA plot (left) of 675 samples from 3 locations with different symbols. The highlighted black points are the identified outliers from the rPSR model. The prediction plot (middle) with highlighted outliers (black). The coefficient plot (right) for rPSR (grey solid) and PSR (black dash) models.

6.1 Computation in rPSR

Computational efficiency is an important property of an algorithm. To examine the computational efficiency for the rPSR algorithm, we repeatedly applied rPSR ( η = 1 ) and PSR on the soil dataset ( m = 675 and p = 214 ) with CEC as the response variable 1 000 times. For each trial, the tuning parameter λ was randomly chosen (i.e., log 10 λ ∼ U ( − 8 , 5 ) ), and used for both rPSR and PSR. Both PSR and rPSR are in R. All the computing and timings were carried out on a computer node with two 10-core 2.8 GHz E5-2680v2 Xeon processors on Louisiana Optical Network Infrastructure (LONI) HPC systems.

Table 5 shows the total CPU timings of PSR and rPSR on 1 000 trials. For PSR, two computing modes are used: PSR i is the ‘inference’ mode of PSR using the signal.fit function from Marx and Eilers (1999); and PSR s is the ‘search’ mode of PSR. The signal.fit program, available at http://statweb.lsu., fits the PSR on multivariate calibration problems under the generalized linear model (GLM) framework. The signal.fit function outputs not only the fitted regression coefficients and predicted values on test samples but also provides a variety of useful model inference statistics and plots. However, we found the PSR s is more computationally efficient than PSR i . In addition, we see rPSR computes 1 000 trials within 7 seconds, which implies rPSR is computationally efficient and well suited for large datasets. The main difference between PSR i and PSR s algorithms are: (a) the latter fits several PSR models with different values of λ in a batch, while the former fits the PSR model on one λ each time. Hence, PSR s avoids unnecessarily recomputing several matrix multiplications in PSR fitting, for example, matrix products: U = XB , U ′ U , D ′ D ; (b) PSR s does not provide the model inference and plots. The timing results in Table 5 indicate that in model selection and cross-validation procedure, using search mode of PSR is highly recommended. After model selection, the final model should be fitted using the inference mode of PSR.

An important factor related to the computational efficiency for the rPSR algorithm is the number of iterations in Step 2. For the above timing study, we also examined the number of iterations for the rPSR algorithm. The five number summary for the number of iterations among 1 000 trials are 7 (minimum), 9 (the first quartile, Q 1 ), 13 (median), 16 ( Q 3 ) and 26 (maximum). We find that although the rPSR algorithm averages approximately 13 iterations, the CPU timing for rPSR is shorter by 5 times to that of PSR. Note that the rPSR comparisons are honest in that they also included the initial PSR solution of α ̂ for Step 2.

OPEN IN VIEWER

Table 5:

CPU timings (in seconds) on PSR and rPSR

	PSR i	PSR s	rPSR
CPU times	28.6	1.44	6.55

In generalized Huber criterion, there are two parameters: K and η . We have an automatic way to select the value of K adapted to the data, described in Section 4.3. For η , we can find the optimal value of η through a grid search. The reason for using generalized Huber criterion rather than the Huber criterion is that the former provides more flexibility on down-weighting the outliers in errors than the latter. In particular, for datasets that are subject to extreme outlying residual, which are commonly encountered in practice, rPSR can adaptively use more robust loss function ( η < 1 ) than the Huber criterion ( η = 1 ). This is partially shown in the simulation study on three types of error distributions. For illustration, we additionally investigated the η ‘effect’ based on a random split of the dataset as in the simulation study. We searched the optimal values of λ and η on a two-dimensional grid by minimizing the tenfold CV error on the training set. Figure 8 shows the standardized tenfold CV RMSE curve (dash) and test RMSE curve (solid) against η on three types of error distribution. From the figure, we see the rPSR using Huber loss ( η = 1 ) performs the best when errors are normally distributed. For mixed normal errors, the optimal value of η lies between 0 and 1. For slash distributed errors, the rPSR using the truncated least-square criterion ( η = 0 ) performs the best.

OPEN IN VIEWER

Figure 8:

Standardized cross-validation RMSE and test RMSE against different values of η on three types of error distribution

Robust penalized regression spline was explored by Lee and Oh (2007), in which Huber loss was used. In Lee and Oh (2007), they proposed an iterative fitting procedure based on the ‘pseudo-response’:

y ̃ i = y ̂ i + ψ ( e i ) 2 ,

where ψ ( · ) is the first order derivative of Huber loss ρ H ( · ) . We can show that the derivative of Huber loss is

ψ = 2 e i − 2 I ( | e i | > K ) [ e i − K Sign ( e i ) ] .

As a result, the pseudo-response y ̃ i is equivalent to

y ̃ i = y i − I ( | e i | > K ) [ e i − K Sign ( e i ) ] .

(6.1)

From Equation (6.1), we can see that the pseudo-response in Lee and Oh's approach is a special case in adjusted response Y A described in (4.11) with η = 1 . However, Lee and Oh's approach is theoretically supported by Cox's result (Cox, 1983; Theorem 3.1 on pp. 536–537), while our approach is motivated by the d.c. programming in An and Tao (1997). In fact, Cox's Theorem 3.1 requires ψ ( · ) to be continuously differentiable up to the second order derivative, (i.e., ψ ∈ C 2 ( − ∞ , + ∞ ) ). However, the proposed generalized Huber criterion ρ η ( · ) is not differentiable at K and − K . Compared with Lee and Oh's procedure, our proposed procedure is more general and motivated from a different perspective.