compiles/166942874ebc3153602751f7ff246cfa/output.dvi


EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS

Vol. 6, No. 2, 2013, 222-238

ISSN 1307-5543 – www.ejpam.com

Smoothing Parameter Selection for Nonparametric Regression

Using Smoothing Spline

Dursun Aydin 1,∗, Memmedaga Memmedli2, Rabia Ece Omay2

1 Department of Statistics, Mugla Sitki Kocman University 48000 Mugla, Turkey
2 Department of Statistics, Anadolu University 26470 Eskisehir, Turkey

Abstract. In this paper, the smoothing parameter selection problem has been examined in respect to

a smoothing spline implementation in predicting nonparametric regression models. For this purpose,

a simulation study has been performed by using a program written in MATLAB. The simulation study

provides a comparison of the nine smoothing parameter selection methods. In this connection, 500

replications have been performed in simulation for sample sets with different sizes. Thus, the appro-

priate selection criteria are provided for a suitable smoothing parameter selection.

2010 Mathematics Subject Classifications: 62E17, 65C05, 65C20, 65C60

Key Words and Phrases: Nonparametric regression, Smoothing spline, Smoothing parameter, Selec-

tion criteria, Generalized cross validation.

1. Introduction

Smoothing spline method is one of the most popular methods used for the prediction of the

nonparametric regression models. The role of this method is to estimate the nonparametric

function that minimizes penalized least squares criterion. A roughness penalty term multiplied

by a positive smoothing parameter is added to the residual sum of squares in smoothing spline

regression. In the light of this approach, the estimation of the unknown function depends on

smoothing parameter λ. Therefore, the determination of an optimum smoothing parameter

in the interval (0,∞) was found to be an underlying complication. In the literature, different

selection methods are components of various studies for an appropriate smoothing parameter.

Indeed, to a considerable extent, Craven and Wahba [5], Hardle [8], Hardle, Hall and Marron

[9], Wahba [26], Hurvich, et al. [11], Eubank [6], Lee and Solo [17], Hastie and Tibshirani

[10], Schimek [22], Cantoni and Ronchetti [4], Ruppert, Wand and Carroll [21], Lee [15, 16],

and Kou [12] supplement on the selection of the smoothing parameter.

In this study, the empirical performances of the selection methods used in selection of the

smoothing parameter are compared. Selection methods used in our simulation study are an

∗Corresponding author.

Email addresses: duaydin@mu.edu.tr (D. Aydin)

http://www.ejpam.com 222 c© 2013 EJPAM All rights reserved.



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 223

improved version of Akaike information criterion (AICc), robustified cross-validation (RCV),

average predictive square error (PSE), parallel of Akaike’s information criterion (GFAIC), gen-

eralized cross-validation (GCV), cross-validation (CV), Mallows’ Cp criterion,risk estimation

using classical pilots (RCP) and local risk estimation (LRS). A simulation study was conducted

to find out which selection methods are the best in smoothing parameter selection. To throw

light on this issue, the samples differing in small and large sizes are secured by means of the

above mentioned simulation, and moreover, nine selection methods are evaluated.

This paper is mainly concerned with the selection of smoothing parameter (or penalty pa-

rameter) through Monte Carlo simulation study. Smoothing parameters play a crucial role in

this procedure. These parameters are said to control the trade off between fidelity to the data

and smoothness: too low values of smoothing parameter overfit the data, whereas too high

values oversmooth. Krivobokova and Kauermann [14] showed that using the REML to esti-

mate smoothing parameter outperforms other methods such as (generalized) CV or Akaike

criterion especially when the error correlation structure is misspecified. Krivobokova et al.

[13] formulated a hierarchical mixed model to estimate local smoothing parameter to achieve

adaptive penalized spline smoothing. Yanrong Cao et al. [27] discussed different methods of

choosing the important smoothing parameter and recommend GCV as the choice for penal-

ized spline smoothing parameter selection for both computational efficiency and accuracy of

the functional coefficient regression models. Aydin and Memmedli [3] recommended GCV

and REML as being good smoothing parameter selection criteria for small and medium sized

samples.

Nonparametric regression and its prediction are discussed in section 2. Section 3 reviews

nine different smoothing parameter selection methods. Section 4 compares these methods via

a simulation study, and finally, the conclusion and recommendations are presented in section

5.

2. Nonparametric Regression Model and its Prediction

Nonparametric regression model including a predictor (independent) variable and a re-

sponse variable is defined as

yi = f (x i) + εi , a < x1 < ...< xn < b (1)

where f ∈ C2[a, b] is an unknown smooth function, (yi)
n
i=1 are observation values of the

response variable y , (x i)
n
i=1 are observation values of the predictor variable x and (εi)

n
i=1 are

normal distributed random errors with zero mean and common variance σ2 (εi˜N(0,σ2)).

The basic aim of the nonparametric regression is to estimate unknown function f ∈
C2[a, b] (the class of all functions f with continuous first and second derivatives) in model

(1). Smoothing spline estimate of the f function appears as a solution to the following mini-

mization problem: Find f̂ ∈ C2[a, b] that minimizes the penalized residual sum of squares

S( f ) =

n
∑

i=1

�

yi − f (x i)
	2
+λ

b
∫

a

�

f ′′(x)
	2

d x (2)



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 224

for pre-specified value λ > 0. The first term in equation (2) denotes the residual sum of the

squares (RSS) and it penalizes the lack of fit. The second term which is weighted by λ denotes

the roughness penalty and it imposes a penalty on roughness. In other words, it penalizes the

curvature of the function f . The λ in (2) is recognized to be the smoothing parameter. As λ

varies from 0 to + , the solution varies from interpolation to a linear model. As λ→ +∞, the

roughness penalty dominates in (2) and the spline estimate is compelled to be a constant. As

λ→→ 0, the roughness penalty disappears in (2) and the spline estimation interpolates the

data. Thus, the smoothing parameter λ plays a key role in controlling the trade-off between

the goodness of fit represented by
n
∑

i=1

�

yi − f (x i)
	2

and smoothness of the estimate measured

by

b
∫

a

�

f ′′(x)
	2

d x .

The solution based on smoothing spline for minimum problem in the equation (2) is

known as a “natural cubic spline” with knots at x1, . . . , xn. From this point of view, a special

structured spline interpolation which depends on a chosen value λ develops into a suitable

approach of function f in model 1. Let f = ( f (x1), . . . , f (xn)) be the vector of values of

function f at the knot points x1, . . . , xn. The smoothing spline estimate f̂λ of this vector or

the fitted values for data y = (y1, . . . , yn)
T are projected by

f̂λ =













f̂λ(x1)

f̂λ(x2)
...

f̂λ(xn)













(n×1)

= Sλ













y1

y2
...

yn













(n×1)

or f̂λ = Sλy (3)

where f̂λ is a natural cubic spline with knots at x1, . . . , xn for a fixed λ > 0, and Sλ is

a well-known positive-definite (symmetrical) smoother matrix which depends on λ and the

knot points x1, . . . , xn, but not on y . Function f̂λ, the estimation of function f , is obtained by

cubic spline interpolation that rests on condition f̂ (x i) = ( f̂ )i , i = 1,2, . . . , n. To gain better

perspective on smoothing spline, Eubank [6], Green and Silverman [7] and Wahba [26] state

studied opinions.

3. Smoothing Parameter Selection Methods

Although smoothing spline estimator solves the problem of allowing fits with variable

slope, a new dilemma emerges. In fact, it generates the determination of the appropriate value

for the smoothing parameter λ for a given data set. The same value of λ is unlikely to work

equally well with every data set. As such, the estimation methods have been introduced for the

selection of smoothing parameter λ in equation (2). The positive value λ that minimizes any

smoothing parameter selection methods is selected as an appropriate smoothing parameter.



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 225

3.1. Selection Methods used in Simulation Study

Various smoothing parameter selection methods are featured in the literature. Most of

these suggested methods were implemented in our simulation study. Moreover, a selection

criterion from previous studies in the literature to provide an effective performance was also

introduced in this particular study. The selection criteria used in our simulation study are

classified as:

Average predictive squared error: In selection the smoothing parameter, it is essential

not to try to minimize the mean squared error at each x i , but instead, the focus should centre

on a global measures such as average predictive squared error

PSE(λ) =

¨

1+
t r(SλS

T
λ
)

n

«

σ 2+




(I − Sλ) f




2

n
(4)

Where t r(SλS
T
λ
) is trace of matrix SλS

T
λ

and



(I − Sλ) f



 is norm of matrix (I − Sλ) f . If σ2 is

not known, in practice an estimation for σ2 can be given by

σ̂2 =
RSS(λ∗)

�

n− t r
�

2Sλ∗ − Sλ∗S
T
λ∗
�	 =




(Sλ∗ − I)y




2

�

n− t r
�

2Sλ∗ − Sλ∗S
T
λ∗
�	 ,

where RSS(λ∗) is the residual sum of square from a smooth Sλ∗ y and λ∗ is a pilot λ selected

by any selection methods [see 10].

Cross-Validation: The basic idea of CV is to disregard one of the points
�

x i , yi

	n

i=1 sequen-

tially, to select the smoothing parameter λ that minimizes the residual sum of squares, and

to estimate the squared residual for a smooth function at x i based on the remaining (n− 1)

points. The CV score can be translated as

CV(λ) =
1

n

n
∑

i=1

n

yi − f̂
(−i)

λ
(x i)
o2 ≡ CV(λ) =

1

n

n
∑

i=1

¨

yi − f̂λ(x i)

1− (Sλ)ii

«2

, (5)

where f̂λ is the fit (spline smoother) for n pairs of measurements
�

x i , yi

	n

i=1 with smoothing

parameter λ, and f̂
(−i)

λ
is the fit calculated by leaving out the ith data point and (Sλ)ii is the

ith diagonal element of smoother matrix Sλ [see 26, 7].

Using the approximations (Sλ)ii ≈
¦

SλS
T
λ

©

ii
and 1/
�

1− (Sλ)ii
�2 ≈ 1+ 2(Sλ)ii , signifies

that E {CV (λ)} ≈ PSE(λ) + 2/n+
n
∑

i=1

(Sλ)ii b
2
i (λ) [see 10].

Generalized cross-validation: GCV is a modified form of the CV which is a conven-

tional method for choosing the smoothing parameter. The GCV score which is constructed by

analogy to CV score can be obtained from the ordinary residuals by dividing by the factors

1− (Sλ)ii. The underlying design of GCV is to replace the factors 1− (Sλ)ii in equation (5)

with the average score 1− n−1 t r(Sλ) Thus, by summing the squared corrected residual and

factor
¦

1− n−1 t r(Sλ)
©2

, by the analogy ordinary cross-validation, the GCV score function can



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 226

be procured as follow [see 5, 26]:

GCV (λ) =
1

n

n
∑

i=1

¦

yi − f̂λ(x1)
©2

�

1− n−1 t r(Sλ)
	2
=

n−1




�

I − Sλ
�

y




2

�

n−1 t r
�

I − Sλ
��2

(6)

Mallows’ CP criterion: In the literature, Cp criterion is referred to as an unbiased risk

estimate (UBR). This type of estimate was suggested by Mallows [18] in the regression case,

and applied to smoothing spline by Craven and Wahba [5]. If σ2 is recognized, an unbiased

estimate of the residual sum of squares is provided by Cp criterion:

Cp(λ) =
1

n

n



(Sλ − I)y




2
+ 2σ2 t r(Sλ)−σ2

o

=
1

n

n



 y − f̂λ




2
+ 2σ2 t r(Sλ)−σ2

o

(7)

Unless σ2 is known, in practice an estimation for σ2 can be given by

σ̂2 = σ̂2

λ̂
=

n
∑

i=1

�

yi − f̂λ̂(x i)
�2

t r
�

I − Sλ̂
� =




(Sλ̂− I)y




2

t r
�

I − Sλ̂
� (8)

where λ̂ is pre-chosen with any of the CV, GCV or AICC criteria (λ̂ is an estimate of λ) For

reference, see [15, 16, 26]. According to Hastie and Tibshirani [10] and Ruppert, et. al. [21],

GCV is approximately equal to Cp.

Improved Akaike information criterion: An improved version of a criterion based on the

classical Akaike information criterion (AIC), AICc criterion, is used for choosing the smoothing

parameter for nonparametric smoothers [11]. This improved criterion is defined as

AICC(λ) = log

∑¦

yi − f̂λ(x i)
©2

n
+ 1+

2
�

t r(Sλ) + 1
	

n− t r(Sλ)− 2

= log





�

Sλ − I
�

y




2

n
+ 1+

2
�

t r(Sλ) + 1
	

n− t r(Sλ)− 2
. (9)

This criterion is easy to apply for the selection of smoothing parameter, as can be seen from

the equation (9).

Robustified Cross-Validation: The smoothing parameter λ can be chosen by method of

robustified cross-validation (RCV). This serves to reduce risk of smoothing parameter mis-

specification in small sample sizes. Most of the selections methods have been proposed to

optimize λ are related to a distributional. RCV criterion that does not refer to a distributional

assumption is proposed by Robinson and Moyeed [20]. Robustified cross-validation score is

given by

RCV (λ) = n−1
1+ n−1+ t r
�

Sλ
�2

�

1+ n−1+ t r
�

Sλ
��2





�

I − Sλ
�

y




2
. (10)

The minimum of (10) was empirically obtained by successively searching interlocked grid

intervals of flexible sizes.



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 227

Parallel of Akaike’s information criterion: A result of Stone [24] states that in the

case of independent observation, the sum over i of the values of log-likelihood function in

y = (y1, . . . , yn) with parameters estimated by maximum likelihood based on
�

x i

	

is asymp-

totically equivalent to Akaike’s information criterion (AIC). This result is not directly applica-

ble to splines, instead of a simple, maximum-likelihood estimate and σ2 as additional noise

parameter has to be taken into account. However, there is parallel of AIC [22]. This parallel

criterion is given by

GFAIC = n−1




�

I − Sλ
�

y




2
+ exp(2n−1 t r(Sλ)) (11)

Minimization of GFAIC looks like maximization of the log-likelihood function in y .

Risk estimation using classical pilots: Risk function measures the distance between the

actual regression function ( f ) and its estimation ( f̂λ). Needless to say that, a good estimate

must contain minimum risk. A direct computation leads to the bias-variance decomposition

for R( f , f̂λ):

R( f , f̂λ) =
1

n
E



 f − f̂λ




2
=

1

n

n




�

Sλ − I
�

f




2
+σ2 t r(SλS

T
λ )
o

(12)

A clear-cut explanation shows that R( f , f̂λ) = E
¦

Cp(λ)
©

. Because the risk R( f , f̂λ) is an

unknown quantity, so-called risk is now estimated by computable quantity R( f̂λp
, f̂λ). The

obtained expression for R( f̂λp
, f̂λ) is

R( f̂λp
, f̂λ) =

1

n
E






 f̂λp
− f̂λ







2

=
1

n

�






�

Sλ − I
�

f̂λp







2

+ σ̂2
λp

t r(SλS
T
λ )

�

, (13)

where σ̂2
λp

and f̂λp
are the appropriate pilot estimates for σ2 and f , respectively. The pilot λp

selected by classical methods is used for computation of the pilot estimates.

Local risk estimation: The LRS method proposed by Lee [16], aims to select the f̂λ(x i)

that minimizes the local risk Rλ(x i) = E
¦

f (x i)− f̂λ(x i)
©2

for the each knot points x i A direct

computation leads to the bias-variance decomposition for Rλ(x i):

Rλ(x i) =
��

Sλ f
�

(x i)− f (x i)
	2
+σ2sλ(x i) (14)

In the above equation,
�

Sλ f
�

(x i) is the ith element of vector Sλ f and sλ(x i) is the ith diagonal

element of the square matrix SλS
T
λ

. An estimator for Rλ(x i) is firstly computed and the f̂λ(x i)

is selected in order to minimizes it. This process is repeated for all x i ’s and at the end of the

process a final mixed estimate for f is derived. The LRS method can be practically performed

with the following five steps [15]:

(i) For a set of pre-selected smoothing parameters λ1 < . . .< λm, calculate the correspond-

ing set of smoothing spline estimates: F =
¦

f̂λ1
, . . . , f̂λm

©

;

(ii) Select the pilot value λp from AICc criterion in (9) by using the elements in F ;

(iii) For λp, calculate the estimates f̂λp
, and σ̂2

λp
by using the (8);



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 228

(iv) Substitute the pilots f̂λp
and σ̂2

λp
into the expression

R̂λ(x i) =
n�

Sλ f̂λp

�

(x i)− f̂λp
(x i)
o2

+ σ̂2
λp

sλ(x i)

and obtain the estimates R̂λ(x i)

(v) For each x i , find the f̂λ(x i) from F which minimizes R̂λ(x i) and the final estimate accept

the appropriate values f̂λ(x i) for f (x i)

3.2. Other Selection Methods

We also explored most of the selection methods that are not used in our simulation. Har-

dle, Hall and Marron [9] recommended that AIC(λ), F PE(λ), and SH(λ) should not be used

because of their trivial minimum at the no-smoothing point. They are also preference T (λ)

for bandwidth selector, but the statement of Rice [19] that “these result are suggestive”, but

T (λ) hasn’t had a good performance in our simulation study. It is accepted that GLM(λ) es-

timates undersmooths relative to the GCV (λ) estimate (see, Wahba [25]). Moreover, Hastie

and Tibshirani [10], also point out that minimizing ASR(λ) over the smoothing parameter

leads to an interpolation estimate.

Furthermore, we have also performed a simulation study for all the selection methods.

In such cases, the selection methods produced invalidating results when compared with the

selection methods such as AICc, GFAIC, GCV, Cp and RCV used in our simulation study. For

this reason, in this paper, the simulation results of the AIC , F PE, SH, T , GLM and ASR criteria

are not displayed. For smoothing spline, these selection criteria can be given as:

(a) Akaike’s information criterion [2],

AIC(λ) =
1

n





�

I − Sλ
�

y




2
exp(1− n−1 t r(Sλ));

(b) finite prediction error [1],

F PE(λ) =
1

n





�

I − Sλ
�

y




2 1+ n−1 t r(Sλ)

1− n−1 t r(Sλ)
;

(c) a model selector of [23];

SH(λ) =
1

n





�

I − Sλ
�

y




2 �

1− n−12t r(Sλ)
�

(d) Rice’s T criterion [19];

T (λ) =
1

n





�

I − Sλ
�

y




2
/
�

1− n−12t r(Sλ)
�



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 229

(e) A generalization of the maximum likelihood (GLM),

GLM(λ) = y ′
�

I − Sλ
�

y/
�

det+
�

I − Sλ
�
�1/n−m

;

where det+
�

I − Sλ
�

is the product of n−m nonzero eigenvalues of
�

I − Sλ
�

.

(f) Average squared residual,

ASR(λ) =
1

n

n
∑

i=1

¦

yi − f̂λ(x i)
©

2

.

4. Simulation Study

This section reports the results of a Monte Carlo simulation study. This study was con-

ducted to evaluate the performances of the nine selection methods mentioned above. The

experimental setup in this paper is adopted from Professor Steve Marron. By using a program

coded in MATLAB, we generated the samples sized n = 25,50,100,150,200,350,400. The

number of replications was 500 for each of the samples. For each simulated data sets, the

mean squared-errors (MSE) was used for evaluate the quality of any curve estimate f̂ . To find

out if the difference between the MSE median values of any two selection methods is signifi-

cant or not, the paired Wilcoxon tests were assessed. In this way, methods which complement

the best smoothing parameter were determined by evaluating so-called selection methods.

4.1. Experimental Plan and Generation of the Data

The experimental setup applied at this stage was designed to study the effects of the

following three factors which vary an independent and effective approach: Noise level; Spatial

variation; Variance function.

The setup specification is listed Table 1. The simulation study was performed according to

MATLAB program, and the experimental setup was designed in the following framework:

• Totally three sets of numerical experiments were performed. For each set of experi-

ments, only one of the above three experimental factors (e.g., noise level, degree of

spatial variation and noise variance function) has been altered while the remaining two

have been kept unchanged.

• Within each set of experiments, the factor levels was modified four times (i.e., r =

1,2,3,4) to detect the effects of any experimental factor in Table 1.

• To see the performance of the small and large samples of the selection methods. For

each factor level r within each set of experiments, we generated 7 different samples

with sample sizes n= 25,50,100,150,200,350,400.

• The number of replications was 500 for each of the 84 numerical experiments.



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 230

• We computed the appropriate smoothing spline estimators f̂λ in equality (3) by selecting

the smoothing parameter λ which minimizes the selection methods.

• We used the MSE values to evaluate f̂λ computed according to each of the selection

criterion: MSE = 1

n

n
∑

i=1

¦

f (x i)− f̂λ(x i)
©

2

, ( f̂λ(x i) = ( f̂λ)i), (where f (x i) is value at

knots x i of the appropriate function f defined in Table 1)

• Paired Wilcoxon test was applied to test whether MSE values was considered as the

performance measure of any two methods are significant or not.

• The factor levels indicated as r was changed four times (i.e., r = 1,2,3,4) in order to

detect the effects of any factor from three factors in Table 1.

• By considering 3 factors, 4 factor levels and 7 samples, totally, 84 numerical experiments

were conducted.

Table 1: Averaged Wilcoxon test ranking values for the nine selection methods in small sample sizes

Factors Generic Form Particular Choices

Noise Level yir = f (x i) +σrεi Factors

Spatial variation yir = f (x i) +σrεi σ = 0.2, fr(x) =
p

x(1− x) sin
�

2π{1+2(9−4r)/5}
x+2(9−4r)/5

�

Variance function yir = f (x i) +
p

vr(x i)εi vr(x) = [0.15{1+ 0.4(2r − 7)(x − 0.5)}]2
r = 1, . . . , 4; x i =

i−0.5

n
;εi ∼ iid N(0,1); f (x) = 1.5θ

�

x−0.35

0.15

�

− θ
�

x−0.8

0.04

�

;θ (u) = 1p
2π

exp
�

−u2

2

�

n= 25,50,100,200,350,400 (it was taken seven sample size).

4.2. Experimental Evaluations

For each simulated data set used in the experiments, the MSE values were used in order

to evaluate the quality of any curve estimate f̂ . Paired Wilcoxon tests were applied to test

whether the difference between the median MSE values of any two methods is significant or

not. The significance level used was 5%. The selection methods were also ranked as follows: If

median MSE value of a method is significantly less than the remaining five, it will be assigned

a rank 1. If median MSE value of a method is significantly larger than one but less than the

remaining four, it will be assigned a rank 2, and similarly for ranks 3-9. Methods having non-

significantly different median values will share the same averaged rank, on the other hand,

method or methods having the smallest rank will be superior.

In this simulation study, because totally 84 different configurations are made, it is not

possible to display here all these configurations. Therefore, only 24 different configurations

are given in Figures A1-A6 (in the appendix) for different samples sized n. The head row plots

of the figures A1-A6 display the true regression function together with all typical simulated

data set. The bottom row plots display the boxplots of the logeMSE values for, from left to



D. Aydin, M. Memmedli, R. Omay / Eur. J. Pure Appl. Math, 6 (2013), 222-238 231

right, AICc, RCV, FSE, GFAIC, GCV, CV, Cp, RCP and LRS. The numbers below the boxplots

are the paired Wilcoxon test rankings. For 84 different simulation experiments, the averaged

ranking values of the selection methods according to Wilcoxon tests are tabulated in Tables A1

and A2. All results tables are shown in the appendix; (∗) indicates the selection methods with

the best rankings.

According to the results in Table A1, for small sized samples (for n = 25,50,100), GCV

has had the best empirical performance for all factors. Furthermore, GFAIC and RCV have

shared the better performance after GCV criterion. In accordance with the overall Wilcoxon

test rankings in Table A1, GCV, GFAIC and RCV have also displayed a good performance. As

shown Table A1, because of R
�

f , f̂λ
�

= E
¦

Cp(λ)
©

, Cp and RCP methods produced the same

results under all experimental factors. In this situation, for small sample sizes, for which

reason the effects of the replication of simulation, Cp is approximately equal to its E
¦

Cp

©

.

For small samples, it is observed that PSE has produced the worst performance. According

to Table A2, for large sized samples (for n = 150,200,350400), GFAIC criterion has had the

best empirical performance. Generally it is shown that AICc, RCV and GCV (Cp) criteria have

shared a better performance after than GFAIC. According to the overall Wilcoxon test rankings

in Table A2, GFAIC, AICc and RCV criteria can be ranked in terms of the performance. As

shown in Table A2, generally, Cp and GCV gave the same results. This can be interpreted to

follow the accepted view that GCV is asymptotically equal to GCV [see 10]. That is to say that,

for large sized samples, because of the effect of the replication of simulation, performance of

the GCV and Cp is approximately equal. PSE has also produced the worst performance for

large samples.

5. Conclusion and recommendations

The scores in Tables A3 and A4 in the appendix are obtained by taking the means of the

averaged Wilcoxon test ranking values tie with each of the selection methods in Table A1-A2

respectively.

As shown in Table A3, according to the means of the small samples for all factors, GCV,

GFAIC and RCV criteria have had the best empirical performance respectively. When it is

compared to the other criteria, PSE, Cp(RCP), AICc criteria have resulted in the worst perfor-

mance.

According to the means of the large samples for all factors in Table A4, it is observed that

GFAIC, AICc , GCV (Cp) have had the best empirical performance. However, PSE and CV have

produced the worst result

Finally, by considering the simulation results and evaluations given in the above, the fol-

lowing suggestions have to be taken into account:

• For both large and small samples, GFAIC and GCV are recommended as being the best

selection criteria;

• For especially large samples, the use GFAIC would seem to be more appropriate. As for

small samples, we propose the implementation of GCV criterion;



REFERENCES 232

• For large samples, the implementation of AICc criterion, in addition to GFAIC and GCV

(Cp) criteria would be more beneficial. For small samples, RCV criterion in addition to

GFAIC and GCV criteria should prove fruitful.

Naturally, the above recommended suggestions have to be considered with a fair amount of

caution as they are only an appraisal based on simulation results.

References

[1] H Akaike. Statistical predictor identification. Annals of the Institute of Statistical Mathe-

matics (AISM), 22:203–217, 1970.

[2] H Akaike. A new look at the statistical model identification. IEEE Transactions on Auto-

matic Control, 19(6):716–723, 1974.

[3] D Aydin and M Memmedli. Optimum smoothing parameter selection for penalized least

squares in form of linear mixed effect models. Optimization, iFirst:1–18, 2011.

[4] E Cantoni and E Ronchetti. Resistant selection of smoothing parameter for smoothing

splines. Statisics and Computing, 11:141–146, 2001.

[5] P Craven and G Wahba. Smoothing noisy data with spline functions. Numerische Math-

ematik, 31:377–403, 1979.

[6] R L Eubank. Nonparametric Regression and Smoothing Spline. Marcel Dekker Inc., New

York, 1999.

[7] P J Green and B W Silverman. Nonparametric Regression and Generalized Linear Model.

Chapman & Hall, New York, USA, 1994.

[8] W Hardle. Applied Nonparametric Regression. Cambridge University Press, Cambridge,

1991.

[9] W Hardle, P Hall, and J S Marron. How far are automatically chosen regression smooth-

ing parameters from their optimum? (with discussion). Journal of the American Statis-

tical Association, 83:86–89, 1988.

[10] T Hastie and R J Tibshirani. Generalized Additive Models. Chapman & Hall, London,

1999.

[11] C M Hurvich, J S Simonoff, and C L Tasi. Smoothing parameter selection in nonpara-

metric regression using an improved akaike information criterion. Journal of the Royal

Statistical Society. Series B (Statistical Methodology), 60:271–293, 1998.

[12] S C Kou. On the efficiency of selection criteria in spline regression. Probability Theory

and Related Fields, 127:153–176, 2003.



REFERENCES 233

[13] T Krivobokova, C M Crainiceanu, and G Kauermann. Fast adaptive penalized splines.

Journal of Computational and Graphical Statistics, 17:1–20, 2008.

[14] T Krivobokova and G Kauermann. A note on penalized spline smoothing with correlated

errors. Journal of the American Statistical Association, 102:1328–337, 2007.

[15] T C M Lee. Smoothing parameter selection for smoothing splines: A simulation study.

Computational Statistics & Data Analysis, 42:139–148, 2003.

[16] T C M Lee. Improved smoothing spline regression by combining estimates of different

smoothness. Statistics & Probability Letters, 67:133–140, 2004.

[17] T C M Lee and V Solo. Bandwidth selection for local linear regression: A simulation

study. Computational Statistics & Data Analysis, 14:515–532, 1999.

[18] C Mallows’. Some comments on cp. Technometrics, 15:661–375, 1973.

[19] J Rice. Bandwidth choice for nonparametric regression. Annals of Statistics, 12:1215–

1230, 1984.

[20] T Robinson and R Moyeed. Making robust the cross-validation choice of smoothing pa-

rameter in spline regression. Communications in Statistics - Theory and Methods, 18:523–

539, 1989.

[21] D Ruppert, M P Wand, and R J Carroll. Semiparametric Regression. Cambridge University

Press, Cambridge, 2003.

[22] G M Schimek. Smoothing and Regression: Approaches, Computation, and Application.

John Willey & Sons, Inc., USA, 2000.

[23] R Shibata. An optimal selection of regression variables. Biometrika, 68(1):45–45, 1981.

[24] M Stone. An asymptotic equivalence of choice of model by cross-validation and akaike’s

criterion. Journal of the Royal Statistical Society. Series B (Statistical Methodology),

39:44–47, 1977.

[25] G Wahba. A comparison of gcv and gml for choosing the smoothing parameter in the

generalized spline smoothing problem. Annals of Statistics, 13:1378–1402, 1985.

[26] G Wahba. Spline Model For Observational Data. Siam, Philadelphia, 1990.

[27] C Yanrong, Z W Tracy L Haiqun, and Y Yan. Penalized spline estimation for functional

coefficient regression models. Computational Statistics and Data Analysis, 54:891–905,

2010.

Simulation Results



REFERENCES 234

Figure A1: Simulation results correspond to the noise level factor for n= 25.

Figure A2: Simulation results correspond to the spatial variation factor for n= 50.

Figure A3: Simulation results correspond to the variance factor for n= 100.



REFERENCES 235

Figure A4: Simulation results correspond to the noise level factor for n= 200.

Figure A5: Simulation results correspond to the spatial variation factor for n= 350.

Figure A6: Simulation results correspond to the variance factor for n= 400.



REFERENCES 236

Table A1: Averaged Wilcoxon test ranking values for the nine selection methods in small sample sizes

Criteria Noise level Spatial variation Variance function Overall

For n= 25

AICc 4.75 5.375 5.25 5.125

RCV 4.875 3.250* 4.250* 4.125

PSE 8 8.5 7 7.83

GFAIC 4.625 3.250* 4.250* 4.042

GCV 3.750* 3.250* 4.250* 3.708*

CV 4.75 5.25 5.25 5.083

Cp 4.75 5.375 5.25 5.125

RCP 4.75 5.375 5.25 5.125

LRS 4.75 5.375 4.250* 4.792

For n= 50

AICc 6.125 6.25 5.625 6

RCV 3.625 2.625* 2.750* 3

PSE 7 7.625 9 7.875

GFAIC 3.625 2.625* 2.750* 3

GCV 2.625* 2.625* 2.750* 2.667*

CV 4.375 7.625 6.125 6.042

Cp 6.125 6.25 5.625 6

RCP 6.125 6.25 5.625 6

LRS 5.375 5.375 4.75 5.167

For n= 100

AICc 6.625 6.625 6.75 6.417

RCV 2.875 2.875 2.000* 2.583

PSE 8.125 8.5 9 8.542

GFAIC 2.875 2.000* 2.000* 2.292

GCV 2.000* 2.000* 2.000* 2.000*

CV 6.625 5.375 5.375 5.792

Cp 6.625 6.625 6.75 6.417

RCP 6.625 6.625 6.75 6.417

LRS 4.875 4.875 4.375 4.7



REFERENCES 237

Table A2: Averaged Wilcoxon test ranking values for the nine selection methods in large sample sizes

Criteria Noise level Spatial variation Variance function Overall

For n= 150

AICc 4.25 3.500* 3.250* 3.667

RCV 3.500* 4.75 3.250* 3.833

PSE 9 8.625 8.875 8.833

GFAIC 3.500* 3.500* 3.250* 3.417*

GCV 4.25 3.500* 4.875 4.208

CV 6.375 6.875 7 6.75

Cp 4.25 3.500* 3.250* 3.667

RCP 5.125 4.5 5.75 5.125

LRS 4.75 6.25 6 5.667

For n= 200

AICc 4.000* 3.750* 2.500* 3.417*

RCV 4.000* 3.875 2.500* 3.458

PSE 9 9 9 9

GFAIC 4.000* 3.750* 2.500* 3.417*

GCV 4.000* 3.750* 4.625 4.125

CV 6.25 6.75 6.5 6.5

Cp 4.000* 3.750* 4.625 4.125

RCP 5 4.75 6.5 5.417

LRS 4.75 5.625 5.25 5.208

For n= 350

AICc 3.875 3.375* 2.875* 3.375

RCV 6.5 6.625 7 6.708

PSE 9 8.75 9 8.917

GFAIC 2.625* 3.375* 2.875* 2.950*

GCV 3.875 3.375* 3.75 3.667

CV 7.75 7.25 5.875 6.958

Cp 3.875 3.375* 3.75 3.667

RCP 5 4.625 7 5.542

LRS 3 4.125 2.875* 3.333

For n= 400

AICc 4.5 3.375* 3.000* 3.625

RCV 6.375 6.75 7.25 6.792

PSE 9 9 9 9

GFAIC 3.500* 3.375* 3.000* 3.292*

GCV 4.5 3.375* 3.000* 3.625

CV 5.75 6 5.25 5.667

Cp 4.5 3.375* 3.000* 3.625

RCP 5.5 5 7.25 5.917

LRS 4.5 4.75 4.25 4.5



REFERENCES 238

Table A3: Means of the averaged Wilcoxon test ranking values for the nine selection methods

Criteria Noise level Spatial variation Variance function Overall

( Means for n= 25,50,100)

AICc 5.667 5.958 5.875 5.833

RCV 3.792 2.917 3.000* 3.236

PSE 7.708 8.208 8.333 8.083

GFAIC 3.708 2.625* 3.000* 3.111

GCV 2.792* 2.625* 3.000* 2.806*

CV 5.125 6.083 5.583 5.597

Cp 5.667 6.083 5.875 5.875

RCP 5.667 6.083 5.875 5.875

LRS 5 5.208 4.458 4.889

Table A4: Means of the averaged Wilcoxon test ranking values for the nine selection methods

Criteria Noise level Spatial variation Variance function Overall

( Means for n= 150,200,350,400)

AICc 4.156 3.406* 2.906* 3.489

RCV 5.094 5.5 5 5.198

PSE 9 8.313 8.969 8.761

GFAIC 3.281* 3.406* 2.906* 3.198*

GCV 4.156 3.406* 4.063 3.875

CV 6.531 6.719 6.156 6.469

Cp 4.156 3.406* 3.656 3.739

RCP 5.156 4.719 6.625 5.5

LRS 4.25 5.186 4.594 4.677