Local Quadratic Estimation of the Curvature in a
Functional Single Index Model1
Zi Ye2 , Giles Hooker3

arXiv:1803.09321v1 [math.ST] 25 Mar 2018

Abstract
The nonlinear effects of environmental variability on species abundance plays an important
role in the maintenance of ecological diversity. Nonetheless, many common models use parametric nonlinear terms pre-determining ecological conclusions. Motivated by this concern, we study
the estimate of the second derivative (curvature) of the link function g in a functional single inR

dex model: Y = g X (t) β 0 (t) dt + ǫ. Since the coefficient function β 0 and the link function g

are both unknown, the estimate is expressed as a nested optimization. For a fixed and unknown
β, the link function g and g ′′ are estimated by local quadratic approximation, then the coefficient function β 0 is estimated by minimizing the MSE of the model. In this paper, we derive the


h R

n
i2
R
P
rate of convergence of the estimation is n1
Xi β 0
E ĝ ′′
Xi β̂ − g ′′
= O h4n + nh1 4 ,
n

i=1

where hn is the bandwidth in the local quadratic approximation. In addition, we prove that the
R
argument of g, Xβ 0 , can be estimated root-n consistently. However, practical implementation
of the method requires solving a nonlinear optimization problem, and our results show that the

estimates of the link function and the coefficient function are quite sensitive to the choices of
starting values.

1

Introduction

1.1

Ecological Motivation

Within mathematical ecology, nonlinear responses to environmental variability play an important
role in maintaining the diversity of competing species. Species competing for the same resources
can nonetheless co-exist by exploiting differing environmental conditions; see Hutchinson (1961),
Chesson and Warner (1981) and Ellner (1987). For an individual species, environmental fluctuation
can accelerate growth rate (Drake, 2005; Koons et al., 2009) or sometimes decrease long-term pop1

Supported in part by NSF Grant DEB-1353039 and DMS-1407600
PhD Candidate, Department of Statistical Sciences, Cornell University, zy234@cornell.edu
3
Associate Professor, Department of Statistical Sciences, Cornell University, gjh27@cornell.edu

2

1

ulation growth rates (Lewontin and Cohen, 1969). The nonlinearity of these responses also plays
an important role in forecasting the effect of increased environmental variability under climate
change. The motivating data for this study come from long-term observations of communities of
prairie plants in which Artemesia Triparta – sage brush – is a dominant species and we wish to
understand its responses to climate given by temperature and precipitation.
Traditional statistical models for plant growth make parametric assumptions that imply specific
forms of nonlinearity, particularly in the presence of high-dimensional covariates. Instead, we use
a nonparametric growth model of an individual plant or animal:
G = g (E) + ǫ,
where G and E are the growth and environment of a plant, g is a link function to be estimated
and ǫ is the random error. To answer the ecological question, “Would the growth be higher if we
just gave the plant a constant environment at the average of E?”, we need to compare g [E (E)]
and E [g (E)].
If the link function g is convex, g [E (E)] ≤ E [g (E)] by Jensen’s inequality, and the plant grows
better in a varying environment. Otherwise, if the link function g is concave, a constant environment
is preferred. Assuming a smooth function g, convexity is equivalent to g′′ (s) > 0, for all s in the
domain of g. Therefore, in this paper we consider the problem of estimating the curvature of the
link function g.
To finalize this model, the environment E is described by the recent history of temperature and
rainfall recorded at up to daily resolution. Since plants may be impacted by climate events over
a long period of time (see Dahlgren and Ehrlén, 2011; Clark et al., 2011), we will consider the
past two years of data. Following Teller et al. (2016), these are thought of as functional covariates
leading to a representation of E as a functional linter term:
E=

Z

1

X (t) β 0 (t) dt,

0

where β 0 (t) is the coefficient function to be estimated, and X (t) is the covariate function we
observed, typically a measurement of climate history.

2

The growth model of a plant is now given by

Z 1
0
X (t) β (t) dt + ǫ.
Y =g
0

This is Functional Single Index model, introduced in Chen et al. (2011) and Ma (2016).
In functional data analysis, a functional linear model (FLM) is defined as
Z 1
X (t) β 0 (t) dt + ǫ,
Y =
0

which is often used in modeling the relationship between a functional covariate and a scalar response.
To assess curvature, we need a more flexible model than the FLM. A generalized functional linear
model (GFLM) is proposed in Müller and Stadtmüller (2005), James (2002) and Escabias et al.
(2007), which has the same form as the functional single index model but with a known link function g. The functional single index model could be considered as an extension to the GFLM, as it
is more flexible and could model a variety of real-world data.
Compared to the generalized functional linear model, estimation of the link function g based on
a unknown coefficient function β 0 is challenging. Even if β 0 is known, estimating the second
derivative of a nonparametric function directly is difficult. In this paper, we prove a theoretical
convergence rate for an estimate of g′′ in the functional single index model, even if there are some
bias in estimating the coefficient function β 0 .
The convergence rates that we derive are based on finding a global solution to a nonlinear optimization problem using a bandwidth that decreases at a known rate with n. However, this
requires overcoming several practical issues. First, to find an optimum for β 0 , we rely on nonlinear
optimization methods which require an initial value from which to search for a minimum. Our
experiments demonstrate that the performance of the estimate can depend critically on this choice
of initial condition and natural choices which provide good estimates of g do not necessarily work
well for g′′ . Further, the optimal choice of bandwidth can be quite different between estimation
targeting g and that targeting g′′ and we provide a heuristic post-cross-validation modification to
improve the estimate of bandwidth. We expect similar rates of convergence will hold for alternative
non-parametric estimators of g, penalized splines, for example, but that the specifics of smoothing
parameter selection and nonlinear optimization can be expected to be quite different. A detailed
analysis of the optimization problem is beyond the scope of this paper.
3

1.2

Previous Results

In this section, we will introduce previous theoretical and empirical results for Single Index model
and Functional Single Index model.
1.2.1

Single Index Model

There has been considerable research on the single index model, where the coefficient β 0 is finite
dimensional. The Single Index Model is defined as

Y = g Xβ 0 + ǫ,

where X is the covariate and β 0 is the coefficient vector. There are three methods to estimate the
link function g and the coefficient vector β 0 . The Projection Pursuit Regression (PPR) approach
introduced in Hardle et al. (1993) is a two-step estimation procedure:
1. Estimate the link function g by the kernel method


Xi β−Xj β
Y
K
j6=i j
h
 ,

ĝi (Xi β |β, h ) = P
Xi β−Xj β
K
j6=i
h
P

where h is the bandwidth.

2. Estimate the coefficient β 0 by minimizing the mean squared error
Ŝ (β, h) =

n
X
i=1

[Yi − ĝi (Xi β |β, h )]2 .

Hardle et al. (1993) proved that the coefficient vector β 0 can be estimated root-n consistently.
Ichimura (1993) showed the asymptotic normality of the estimator. The other two approaches
provide new methods to estimate the coefficient vector. The Average Derivative approach in
Hristache et al. (2001) showed that



 .

∂g (z)
∂g (z)
z =Xβ
.
= E
β =E
β = γβ.
∂z
∂z


∂g (Xβ0 )
If we could find a consistent estimator of the average derivative E
, we can get a consistent
∂X


∂g (Xβ)
E
∂X

estimator of the coefficient β 0 up to a scale. Normally, we require the coefficient vector to be norm
4

1. Stoker (1986) proposed two consistent estimators of the average derivative.
The sliced inverse regression method in Li (1991) considered the estimation of the coefficient vector
as a dimension-reduction problem. Any linear combination of the coefficient vector β 0 is assumed
to be an effective dimension-reduction (EDR) direction. They conduct a principle component
analysis on the inverse regression space E (X |Y ), and estimate the coefficient vector β 0 by the
largest component.
1.2.2

Functional Single Index Model

There are only a few papers in the functional single index model. In Chen et al. (2011), similar to
the projection pursuit regression in the single index model, the link function g and the coefficient
function β 0 are estimeated by a two-step procedure. The coefficient function β 0 is reduced to a
finite dimensional coefficient vector by a spline basis. Under some assumptions, Chen et al. (2011)
showed that

Z
2
n  Z

1X
0
ĝ
Xj β̂ − g
Xj β
= O n−c ,
n
j=1

for c > 0. In Ma (2016), two spline basis were used to represent the coefficient function and
the link function, respectively, and the MSE was minimized iteratively until convergence. Ma
(2016) constructed a asymptotic simultaneous confidence band for the coefficient function β 0 . Our
estimates follow Chen et al. (2011) but will examine the properties of ĝ′′ . By a clever decomposition
of squared error, Chen et al. (2011) were able to avoid the need to directly account for the estimate
R
√
of β 0 . Unlike that case, to examine g′′ we will need to obtain the n convergence rate for Xβ 0
directly, before we can examine our target.

2

Estimation Procedure

Suppose that we observe n environment histories and responses (X1 (t) , Y1 ) , · · · , (Xn (t) , Yn ), independent and identically distributed as (X (t) , Y ), where t ∈ [0, 1], with
Y =g

Z

1

0



X (t) β (t) dt + ǫ,
0

5

where Y is the scalar response variable, and X (t) is the covariate function. For the purpose of
simplification, we assume that the predicator X and the coefficient function β 0 are defined in the
domain [0, 1], and ǫ is a Gaussian random error.
To answer our ecological question, we are interested in estimating the second derivative (curvature)
of the link function g. The estimate of the coefficient function β 0 is denoted as β̂. Define a Hilbert
space B as the set of the coefficient functions β, where β̂ (t) , β 0 (t) ∈ B.
To estimate β 0 , g and g′′ , we use a local quadratic approximation. By Taylor’s expansion, at a
fixed point x, the link function g can be approximated by
g (x) ≈ g (x0 ) + g′ (x0 ) (x − x0 ) + g′′ (x0 )

(x − x0 )2
.
2

Fix u, where u is in the domain of the link function g, the curvature, denoted as ĝ′′ , is estimated
by minimizing the weighted sum of squares


â, b̂, ĉ =

R
2 2
1


Z
n

1
X 
0 Xi (t) β (t) dt − u

inf
Xi (t) β (t) dt − u − c
 Yi − a − b
 K

a,b,c
2
0

i=1

(1)
R1
0


!
Xi (t) β (t) dt − u 
,

hn


where K is a kernel function and hn is the bandwidth. The estimators are then (ĝ (u) , ĝ ′ (u) , ĝ ′′ (u)) =


â, b̂, ĉ .
The coefficient function β 0 (t) is unknown in the penalized weighted sum of square (1). We estimate
it by minimizing the MSE
β̂ = inf

β∈B

n 
X
i=1

Yi − ĝ

Z

Xi β

2

,

(2)

where

ĝ

Z

Xi β



R

R
Xi β− Xj β
Yj K
hn
j:j6=i
R
 .
R
= P
Xi β− Xj β
K
hn

P

(3)

j:j6=i

Since the kernel function K is only defined in [−1, 1], we constraint the domain of the estimate of
g or g ′′ to be in [−1, 1] by normalizing the coefficients of g or g ′′ under an orthonormal basis to be
1 after optimization procedure.
6

Denote a column vector Y = (Y1 , · · · , Yn )⊤ . Fix j ∈ {1, · · · , n}, and β ∈ B, the estimated



R
R
R
ĝ Xj β , ĝ′
Xj β and ĝ ′′
Xj β can be calculated as


⊤ 
 Z
Z
Z

−1 
⊤
K
(4)
K
X
X
= X⊤
Xj β , ĝ′
Xj β , ĝ′′
Xj β
ĝ
β,j β,j Y ,
β,j β,j β,j
where
ĝ ′
ĝ

′′

Z

Z

Xj β
Xj β







−1 
.
⊤
X
K
X⊤
K
X
β,j β,j Y = S1 (β; j) Y ,
β,j β,j β,j

=

2·


−1 

.
⊤
K
Y = S2 (β; j) Y ,
X
K
X
X⊤
β,j
β,j
β,j
β,j
β,j

=

3·

(5)
(6)

where Ak· denotes the kth row of a matrix A, and the (n × 3)-dimensional matrix Xβ,j is
R
R
 2 !
Z

Z
Xβ −
Xj β 1
Xβ,j = 1, Xβ −
Xj β 1,
,
2
with

R

R
⊤
. R
Xβ =
X1 β, · · · , Xn β , 1 is a n-dimensional column vector of ones, and the (n × n)-

dimensional matrix Kβ,j is



Kβ,j = diag K

R

R
R

R

X1 β − Xj β
Xn β − Xj β
,··· ,K
.
hn
hn

The estimation of the coefficient function β 0 and the link function g is therefore a nested procedure,
summarized in (2), (3) and (4). Following Ma (2016), the identifiability of the model is ensured by
R1
adding a constraint on the coefficient function, such that 0 β 2 (t) dt = 1.

3

Assumptions

In deriving a convergence rate for

1
n

n
P

i=1

h R

i2
R
E ĝ ′′
Xi β̂ − g ′′
Xi β 0 , we make the following as-

sumptions in the functional single index model.

1. The observations (Xi (t) , Yi ), where i = 1, · · · , n, are independent and identically distributed.
Each covariate function Xi (t) is a square-integrable random function defined in the interval
[0, 1]. The random error ǫ is independent from X, and has zero mean and variance σ 2 .
2. The dependent variable Y has the mth-order absolute moment, where m ≥ 2. This is an
assumption from Ichimura (1993). The finite moment m is used in establishing the main
convergence theorem.
7

3. The link function g and the curvature g′′ are bounded and satisfy the Lipschitz condition
such that
g (k) (u) − g (k) (v) ≤ D2 |u − v| ,
for all u and v, where D2 > 0 and k = 0, 2. The Lipschitz condition ensures that if β 0 can
R

R

X β̂ can be
be estimated root-n consistently, the distance between g′′
Xβ 0 and g ′′
controlled.

4. The kernel function K is nonnegative and symmetric with support [−1, 1], and

R1

−1 K

(s) ds =

1. Assume that K is three times continuously differentiable, with K (3) (s) ≤ D3 , for any
s ∈ [−1, 1] and D3 > 0. Since the kernel K satisfies a Lipschitz condition, the NadarayaWatson estimator of the link function g also has a Lipschitz condition.
5. For some orthonormal basis {φk (t) : k = 1, 2, · · · }, for each i = 1, · · · , n, there exists a sequence of random variables {cij }∞
j=1 , such that
Xi (t) =

∞
X

cij φj (t) ,

∞
X

bj φj (t) .

j=1

and
β (t) =

j=1

Assume that E (cij ) = 0.
In particular, we have
β̂ (t) =

∞
X

b̂j φj (t) ,

β 0 (t) =

∞
X

b0j φj (t) .

j=1

j=1

For any β ∈ B, we can write
Z

Xi β =

∞
X

cij bj .

j=1

We observe that an orthonormal basis approximation of the covariate function and coefficient
function transforms an integration to an infinite sum. In addition, define a sequence pn such

8

that pn → ∞ as n → ∞, we require
∞
X

j=pn +1



,
cij bj = O p−λ
n

(7)


R
where λ > 0, and pn = o h1n . Condition (7) ensures that the integration Xβ can be


approximated by a finite sum of coefficients under an orthonormal basis.

6. Assume that supβ∈B;x f (x |β ) < ∞, where f (x |β ) is the probability density of

4

R

Xβ.

Convergence Rates

By the definition of X⊤
β and Kβ , we can calculate


1⊤


 R
 ⊤
R

Xβ
−
X
β
1
j
K
=
X⊤

β
β
  R

R
2 ⊤

( Xβ−( Xj β )1)
2



Denote


 


=




R



 K






R

R

X1 β− Xj β
hn

0

..

0

0⊤ K

⊤
R
Xβ−( Xj β )1
K
hn
R
⊤
R
R
R
 
Xβ−( Xj β )1
Xβ −
Xj β 1 K
hn
⊤
R
 R
R
R
2
Xβ−( Xj β )1
( Xβ−( Xj β )1)
K
2
hn


Tjp (β)

=

n Z
X
i=1

Xi β −

Z

Xj β

p

0⊤

K

R



=

n
X
i=1

K

R

9

.
R

R

Xi β − Xj β
,
hn

R

Xi β − Xj β
.
hn

0

R
Xn β− Xj β
hn





.




and
Tj0 (β)

0







 

. R
Denote uj = Xj β, we have

n
X

R

Xi β − Xj β
K
=
hn
i=1


Z
z − uj
f (z |β ) dz + O (1)
=n K
hn
Z
=n K (m) f (uj + hn m |β ) hn dm + O (1)
Z

=nhn f (uj |β ) K (m) dm + O h2n + O (1)

Z
Xj β |β + O (1) ,
=nhn f

Tj0 (β)

R

and
Tjp (β)

=

n Z
X
i=1

=n

Z

Xi β −

Z

p

(z − uj ) K

Z

Xj β



p

z − uj
hn

K



R

R

Xi β − Xj β
hn

f (z |β ) dz + O (1)

(hn m)p K (m) f (uj + hn m |β ) hn dm + O (1)
Z

p+1
=nhn f (uj |β ) mp K (m) dm + O hp+2
+ O (1)
n

Z
X
β
|β
µp (K) + O (1) ,
=nhp+1
f
j
n
=n

. R1
where µp (K) = −1 mp K (m) dm. Since the kernel function K is symmetric, Tjp (β) = 0 if p is an

odd number.
We have

X⊤
β,j Kβ,j Xβ,j



Tj0 (β)

Tj1 (β)



1
2
≈
 Tj (β) Tj (β)
 T 2 (β)
T 3 (β)
j

j

2

2

Tj2 (β)
2
Tj3 (β)
2
Tj4 (β)
4

The determinant of the matrix X⊤
β,j Kβ,j Xβ,j is

X⊤
β,j Kβ,j Xβ,j =





 
 
=
 
 

0

0

Tj2 (β)

0

0

Tj4 (β)
4

Tj2 (β)
2


3
Tj0 (β) Tj2 (β) Tj4 (β) − Tj2 (β)
4

10

Tj2 (β)
2

Tj0 (β)

.





.



By definitions of S1 (β; j) and S2 (β; j) in (5) and (6), we can get
R
⊤
R

2   R
R
 
Xβ−( Xj β )1
Tj0 (β) Tj4 (β) − Tj2 (β)
Xβ −
Xj β 1 K
hn
,
S1 (β; j) =

3
Tj0 (β) Tj2 (β) Tj4 (β) − Tj2 (β)


S2 (β; j)

2Tj0 (β)
=

4.1



R

Xβ −

R

 2
Xj β 1 K

R

⊤
R
Xβ−( Xj β )1
hn

2Tj2 (β)

−

2
Tj0 (β) Tj4 (β) − Tj2 (β)



K

R

⊤
R
Xβ−( Xj β )1
hn

.

Convergence Rate of β

For any β ∈ B, define


2
Z
2
n 
n
∞
X
X
X
1
. 1
Yi − ĝ 
Jn (β) =
=
Yi − ĝ
Xi β
cij bj  ,
n
n
i=1

and

i=1

j=1



2
pn
n
X
X
1
.
Yi − ĝ 
Jn,pn (β) =
cij bj  .
n
i=1

j=1

Denote ci = (ci1 , · · · )⊤ and b = (b1 , · · · )⊤ , where i = 1, · · · , n. Define the subspaces C, B ⊂ R∞
such that ci ∈ C and b ∈ B.
Lemma 1. For a sequence of positive numbers Mn , suppose that nh2n ǫ20n Mn−2 → ∞, then
(
)

 Z
Z

n 
X
P sup
ĝni
xβ − E ĝni
≥ ǫ0n → 0,
xβ
C×B i=1

where
ĝni

Z

xβ



1
Yi I (Yi ∈ [−Mn , Mn ]) K
=
hn

as n → ∞.

11

R

R

xβ − Xi β
,
hn

Proof. We have
)

 Z
Z

n 
1X
P sup
ĝni
xβ − E ĝni
≥ ǫ0n
xβ
C×B n i=1
(
R

R
n 
X
xβ − Xi β
= P sup
Yi I (Yi ∈ [−Mn , Mn ]) K
hn
C×B i=1
R



R
xβ − Xi β
≥ nhn ǫ0n
− E Yi I (Yi ∈ [−Mn , Mn ]) K
hn
!
n
X
.
Ai ≥ nhn ǫ0n ,
= P sup
(

C×B i=1

where
Ai = Yi I (Yi ∈ [−Mn , Mn ]) K

R

R
R


R

xβ − Xi β
xβ − Xi β
− E Yi I (Yi ∈ [−Mn , Mn ]) K
.
hn
hn

By Assumption 4, the kernel function K is bounded. We will apply Bernstein’s inequality (see
Appendix A) to the above equation with
ηn = nhn ǫ0n ,
.
|Ai | ≤ 2 |Yi | K1 ≤ 2Mn C1 = cn ,
and
var (Ai ) ≤ Mn2 C2 ,
n

X
.
var (Ai ) ,
Vn = nMn2 C2 ≥
i=1

where K1 , C1 and C2 are constants. By Bernstein’s inequality, we can get
(
)

 Z
Z

n 
1X
P sup
ĝni
xβ − E ĝni
≥ ǫ0n
xβ
C×B n i=1
#
"
ηn2

≤ exp −
2 Vn + 31 cn ηn
"
#
(nhn ǫ0n )2

= exp −
2 nMn2 C2 + 32 Mn C1 nhn ǫ0n
"
#
nhn ǫ20n
= exp − 2M 2 C
.
4
n 2
+
M
C
ǫ
n
1
0n
hn
3
12

The assumption of Lemma A.5 in Ichimura (1993) is that the sequence {Mn }∞
n=1 should satisfy
ǫ0n hn Mnm−1 → ∞. Since hn → 0 and ǫ0n → 0, we need Mn → ∞. Therefore, in the denominator,
 2 
we have 34 Mn K1 ǫ0n = o 2MhnnK2 . If nh2n ǫ20n Mn−2 → ∞, then
− 2M 2 K
n

hn

nhn ǫ20n

2

+ 34 Mn K1 ǫ0n

→ −∞,

and
)

 Z
Z

n 
1X
P sup
ĝni
xβ − E ĝni
≥ ǫ0n → 0.
xβ
C×B n
(

i=1

Following similar arguments, we can derive Lemmas 2 and 3 below.
Lemma 2. For a sequence of positive numbers Mn , suppose that nh4n ǫ21n Mn−2 → ∞, then
(
)
Z

 Z

n 
X
′
′
ĝni
P sup
xβ − E ĝni
≥ ǫ1n → 0,
xβ
C×B i=1

where
′
ĝni

Z

xβ



R
R
′
xβ − Xi β
1
Yi I (Yi ∈ [−Mn , Mn ]) K
=
hn
hn
R
R

xβ − Xi β
. 1
′
= 2 Yi I (Yi ∈ [−Mn , Mn ]) H1 (x, Xi , β) K
,
hn
hn


as n → ∞.
Lemma 3. For a sequence of positive numbers Mn , suppose that nh6n ǫ22n Mn−2 → ∞, then
(
)

 Z

Z
n 
X
′′
′′
P sup
≥ ǫ2n → 0,
ĝni
xβ − E ĝni
xβ
C×B i=1

where
′′
ĝni

Z

xβ



R
R
′′
xβ − Xi β
1
Yi I (Yi ∈ [−Mn , Mn ]) K
=
hn
hn
R

R
xβ − Xi β
1
.
′′
= 3 Yi I (Yi ∈ [−Mn , Mn ]) H2 (x, Xi , β) K
,
hn
hn


as n → ∞.
13

Now, we show the root-n consistency of the estimator β̂.
Theorem 4. The estimator β̂ is consistent if Assumptions 1− 6 hold, and the bandwidth sequence
satisfies nh8n → 0 and nh6n → ∞.
Proof. Theorem 5.1 in Ichimura (1993) states the consistency of β̂, while β 0 is a coefficient vector.
Theorem 5.1 is based on Lemma 5.1, and Lemma 5.1 is based on Lemmas A.2 − A.10. The proof
of Lemmas A.2 − A.7 will be the same whenever β 0 is a vector or a function. Lemmas 1 − 3 above
are the functional version of Lemmas A.8 − A.10. We need to figure out the constraint on the
smoothing parameter hn . The constraints are
• Lemma A.2 − A.4
nh8n → 0.
• Lemma A.5 − A.7
ǫ0n hn Mnm−1 → ∞,

ǫ1n h2n Mnm−1 → ∞,

ǫ2n h3n Mnm−1 → ∞.

nǫ21n h4n Mn−2 → ∞,

nǫ22n h6n Mn−2 → ∞.

• Lemma A.8 − A.10
nǫ20n h2n Mn−2 → ∞,

We require nh2n ǫ20n Mn−2 → ∞. Since ǫ0n → 0 and Mn−2 → 0, we need to have nh2n → ∞. Following
the same argument, we need nh4n → ∞ and nh6n → ∞.
To prove the convergence rate of the functional single index model, we need to find a convergence
R
rate for β̂, b̂ or X β̂.
3+

Theorem 5. Suppose that

nh6n

→ ∞,

√

n

Z

nh8n

→ 0 and

Xi β̂ −

Z

Xi β

for i = 1, · · · , n.

14

3

nhn m−1
− log hn

0



→ ∞, then we have

= Op (1) ,

Proof. For each pn such that pn → ∞ as n → ∞, define



b̂pn = b̂1 , · · · , b̂pn , and, b0pn = b01 , · · · , b0pn .
Since the kernel function K satisfies the Lipschitz condition, and the estimated link function g is
R

R
P
Xi β− Xj β
Y
K
Z

j
hn
j:j6=i
R
 ,
R
ĝ
Xi β = P
Xi β− Xj β
K
hn
j:j6=i

∞
P

the estimated ĝ also satisfies the Lipschitz condition. Since

j=pn +1





pn
∞
X
X
cij bj  − ĝ 
cij bj  ≤ K2
ĝ 


j=1

j=1

as n → ∞, where K2 is a constant. Therefore, ĝ

pn
P

j=1

∞
X

j=pn +1

cij bj

!


cij bj = O p−λ
n , we have



,
cij bj = O p−λ
n

→ ĝ

Jn,pn (β), we can get Jn,pn (β) → Jn (β) as n → ∞, for any β ∈ B.

∞
P

(8)

!

cij bj . By the definition of

j=1

Hence, by Lemma 5.4 in Ichimura (1993), we have


√ 
.
n b̂pn − b0pn → D = N 0pn , V−1 ΣV−1 ,

where 0pn is a pn -dimensional mean vector, and V−1 ΣV−1 is a (pn × pn )-dimensional covariance
matrix. Suppose that X is the σ-algebra generated by (X1 , · · · , Xn ), the (k, m)th-term of the
matrix V is
"
#

∂ 2 Jn,pn β 0
Vkm =E
|X
∂b0k ∂b0m



 

2
 

pn
pn
pn
n
n
2 X

X
X
X
X
2
ĝ′ 
Yi − ĝ 
=E
cij b0j  cik cim |X
cij b0j  ĝ ′′ 
cij b0j  cik cim −
n

n
i=1
i=1
j=1
j=1
j=1





 

2
 

pn
pn
pn
n
n


2 X
2 X
X
X
X
ĝ′ 
Yi − ĝ 
cij b0j  cik cim |X .
cij b0j  ĝ′′ 
cij b0j  cik cim |X − E
=E


n
n
i=1

i=1

j=1

j=1

j=1

For any i = 1, · · · , n, we have

E (Yi |X ) = g

Z

Xi β 0



 Z
 

= E ĝ
Xi β 0 |X + O h2n ,
15

(9)

where the second equality is the bias property of the kernel density estimate. We can calculate

 


pn
n
1 X

X
Yi − ĝ 
E
cij b0j  |X
n

i=1
j=1


 

( n 
Z
 )
Z

pn
n
1 X

X
1X
ĝ
Yi − ĝ
Xi β 0 |X + E
Xi β 0 − ĝ 
cij b0j  |X ,
=E
n

n
i=1

i=1

j=1

where the first term converges to 0 by the equation (9), since hn → 0 as n → ∞, and the second
term converges to 0 by (8). Therefore, we have

 


pn
n
1 X

X
Yi − ĝ 
cij b0j  |X → 0,
E
n

i=1

j=1

as n → ∞. By the Slutsky’s Theorem, we can get


 



pn
pn
n
2 X

X
X
Yi − ĝ 
E
cij b0j  ĝ ′′ 
cij b0j  cik cim |X → 0.
n

i=1

j=1

j=1

Therefore,

Vkm



2
 
pn
n

2 X
X
ĝ′ 
cij b0j  cik cim |X .
=E

n
i=1

j=1

For any β ∈ B and any j ∈ {1, · · · , n}, we have
ĝ′

Z

Xj β



= |S1 (β; j) Y |


2  P
n R
  R Xi β−R Xj β 
R
0
4
2
Tj (β) Tj (β) − Tj (β)
Xi β − Xj β K
Yi
hn
i=1
=

3
Tj0 (β) Tj2 (β) Tj4 (β) − Tj2 (β)


2 
0
4
2
Tj (β) Tj (β) − Tj (β)
Tj1 (β) maxi∈{1,··· ,n} Yi
≤
 

3
Tj0 β̂ Tj2 (β) Tj4 (β) − Tj2 (β)
h
2 i 2
nhn f (u |β ) µ1 (K) maxi∈{1,··· ,n} Yi
nhn f (u |β ) nh5n f (u |β ) µ4 (K) − nh3n f (u |β ) µ2 (K)
≈
nhn f (u |β ) nh3n f (u |β ) µ2 (K) nh5n f (u |β ) µ4 (K) − (nh3n f (u |β ) µ2 (K))3
∼

1
,
hn

16

where u =

R

Xj β.

For any i = 1, · · · , n, denote c·· = cik , for any k ∈ {1, · · · , n}, and c·k = cik for a fix k. For any
k, m ∈ {1, · · · , n}, we have
Vkm ∼

E (cik cim |X )
cov (cik , cim )
cov (c·k , c·m )
∼
∼
,
2
2
hn
hn
h2n

since Xi are independent and identically distributed with E (cij ) = 0, for any i = 1, · · · , n and
j = 1, · · · .
The (k, m)th-term of the matrix Σ is
"
#

∂Jn,pn β 0
∂Jn,pn β 0
·σ
Σkm =E σ
∂b0m
∂b0k
 
 








pn
pn
pn
pn
n
n
X
X
X
X
X
X
2σ
2σ
Yi − ĝ 
Yi − ĝ 
cij b0j  ĝ′ 
cij b0j  ĝ ′ 
cij b0j  cik ·
cij b0j  cim 
=E 
n
n
i=1

∼

j=1

i=1

j=1

j=1

j=1

cov (c·k , c·m )
.
h2n

The diagonal term of the covariance matrix V−1 ΣV−1 is
V−1 ΣV−1



kk

∼

p2n h2n
.
var (c·· )

Define the truncated version of Xi (t) and β (t) as
pn
X

Xi,pn (t) =
βpn (t) =

cij φj (t) ,

j=1
pn
X

bj φj (t) ,

j=1

where i = 1, · · · , n.
We have for any i = 1, · · · , n, define ci,pn = (ci1 , · · · , cipn ),
Z

Z


√
√ 
√ 
0
Xi,pn β̂pn − Xi,pn βpn = n ci,pn b̂pn − ci,pn b0pn = ci,pn n b̂pn − b0pn .
n


R
√ R
n
Xi,pn β̂pn − Xi,pn βp0n converges to a normal distribution with the covariance
 
p2n h2n
p2n h2n
1
2 h2 . By Assumption 5, since p = o
matrix var (c·· ) · var(c
=
p
n
p
n
n
hn , we can get var (c·· ) · var(c·· ) =
·· )
Therefore,

op (1). When n → ∞, pn → ∞, we have
Z

Z
√
0
n
Xi β̂ − Xi β
= Op (1) ,

for i = 1, · · · , n.
17

4.2

Main Theorem
3

3+

nhn m−1
− log hn

Theorem 6. If nh6n → ∞, nh8n → 0 and

→ ∞, we have

 Z



Z
2
n
1
1X
′′
4
′′
0
E ĝ
= O hn +
.
Xi β̂ − g
Xi β
n
nh4n
i=1

Proof. Since Xi are independent and identically distributed, for any j ∈ {1, · · · , n}, we only need
h R

 i2
R
′′
′′
0
to find the convergence rate of E ĝ
Xj β̂ − g
Xj β
. We can decompose it into three
terms:



E ĝ

′′

Z


Xj β̂

≤ E ĝ

′′



Z

+E g

′′



−g

Xj β̂

Z

′′

Z



− ḡ

Xj β̂



Xj β
′′

Z

−g

′′

0

2
Xj β̂

Z

2

Xj β

0



+ E ḡ

2

′′

Z

Xj β̂



−g

′′

Z

.

where
ḡ ′′

Z

Xj β̂





.
= S2 β̂; j g,

and
.
g=

 Z

Z
⊤
0
0
X1 β , · · · , g
g
Xn β
.

By Lemma 7, 8 and 9 in the Appendix, we have


E ĝ


′′

Z

Z

Xj β̂



− ḡ



′′

Z

Xj β̂

2

2

Z

=O





1
nh4n



.

1
E ḡ
Xj β̂ − g
=O
+
Xj β̂
nh4n
 Z
 

Z
2
1
.
E g′′
=O
Xj β̂ − g′′
Xj β 0
n
′′

′′

h4n



.

Combining these three terms, we obtain
 Z



Z
2
n
1X
1
′′
4
′′
0
E ĝ
.
= O hn +
Xi β̂ − g
Xi β
n
nh4n
i=1

18

Xj β̂

2

5

Practical Implementation

5.1

Initialization

In (2), the coefficient function β 0 is estimated by minimizing the mean square error of g, which is
a nonlinear optimization problem. The coefficient function is approximated by a K-dimensional
Fourier basis, as
β 0 (t) = ψ ⊤ (t) c,
where c is a K-dimensional column vector, and t ∈ [0, 1]. In order to ensure identifiability, we
R
constraint the coefficient function β 0 to have β 0 = 0 and kβ 0 k2 = 1. Since the Fourier basis

is an orthonormal basis, a constraint on β 0 is equivalent to a constraint on the coefficient vector,

such that kck2 = 1. The first constraint can be enforced by dropping any constant terms from the
Fourier basis. To compensate for rescaling c, we also rescale the bandwidth h such that h = h kck
in the optimization step.
We use R function optim to minimize the MSE, and an initial value of the coefficient vector c is
needed, denoted as cinit . We use three different methods to select the initial value.
1. Assume each item in cinit is equal, such that cinit =



√1 , · · ·
K


, √1K . We don’t have any

previous knowledge about the start point, so for simplicity, choose a vector while all items
are equal to each other.

2. Assume g (s) = s, the coefficient vector cinit is estimated by minimizing the ordinary least
square, and normalize it to be kcinit k2 = 1. To obtain an initial for β 0 , we need to specify a
structure for the link function g. A linear structure of g is obviously the simplest, and could
be calculated easily.
3. Generate 1000 different standard normal distribution random initial vectors (with length c),
and select the best 10 initial vectors by penalized mean squared error, when the bandwidth
is equal to the mean of bandwidth sequence. For each bandwidth h, select the initial that
minimizes the penalized MSE. That means that in the cross-validation step, we select different

19

initial for different bandwidth h. The selected initial is determined by the bandwidth and
corresponding penalized MSE.

5.2

Cross-validation

We need to select the bandwidth hn in the kernel density estimation of g and g′′ . We examine two
cross-validation methods:
1. 10-fold cross-validation. We partition the dataset into 10 subsamples. Each time, we use
9 subsamples as the training set and the remaining subsample as the validation set. We
will observe that the 10-fold cross validation method produces similar results with the GCV
method introduced below.
2. Fix j ∈ {1, · · · , n}, define
Z
 

−1 
.
⊤
K
Y = S0 (β; j) Y ,
X
K
X
ĝ
Xj β = X⊤
β,j
β,j
β,j
β,j
β,j
1·

where S0 (β; j) is a n-dimensional row vector. Denote a (n × n)-dimensional smoother matrix

 

⊤
 R

R
⊤
.
Shn = S0 β̂; 1 , · · · , S0 β̂; n
. We can get Ŷ = Shn Y , where Ŷ = ĝ
X1 β̂ , · · · , ĝ
Xn β̂
.
The generalized cross-validation criterion is

2

1

. k(I − Shn ) Y k2
GCV (hn ) = n
2 ,
1
n tr (I − Shn )

where I is a n-dimensional identity matrix. Find a bandwidth hn minimizing the GCV (hn ).
Note that this does not account for estimating β 0 .

5.3

Simulation Study

In order to obtain simulated functional data, we defined a 25-dimensional Fourier basis ψ (t), where
t ∈ [0, 1]. The covariate function X (t) is defined in [0, 1], such that
X (t) =

25
X

ηi ψi (t) ,

i=1

where ηi =

i−1
24

· N (0, 1). The coefficient function is
β 0 (t) = a⊤ ψ (t) ,
20

where a = (0, 1, 1, 0.5, 0, · · · , 0)⊤ . We use three different link functions:
1. g (s) = e−s .
2. g (s) = −s2 .
3. g (s) = s.
In order to measure the performance of our estimators, we define the MSE of the estimated β 0 and
g(k) to be
RSE =

Z

0

2  12
β̂ (t) − β (t) dt ,

1

0

and
RASE(k) =
(k)

where Ŷi

= ĝ(k)

R
1
0

(

2 ) 21
Z 1
n 
1 X (k)
Ŷi − g(k)
Xi (t) β 0 (t) dt
,
n
0
i=1


Xi (t) β̂ (t) dt .

 1
 1
Theorem 6 shows that for consistency, the bandwidth must scale between O n− 6 and O n− 8 .
1

Assume that the optimal bandwidth for the curvature g′′ is n− 7 , by Chen et al. (2011), we can
1

expect that the optimal bandwidth for the link function g is n− 5 . We use either the GCV or the
5

10-fold cross-validation to select a bandwidth h, and rescale it to be h 7 . As we have discussed
before, we constrain the coefficient β 0

2

= 1; this was achieved by rescaling the bandwidth in our

objective function. After solving the nonlinear optimization problem, we rescale both c and the
bandwidth. The following table shows the simulation results of the GCV method. The number
of data points we use are n = 100 and n = 1000. In Table 2, we also show the RASE2 results
5

without rescaling the bandwidth h to be h kck and h 7 , if we start from random initials. We can
conclude that re-scaling does matter to the final results and it reduces the error of RASE2. The
10-fold cross-validation results in the Appendix C also confirm that by re-scaling the bandwidth,
we improve our estimate for g′′ .

21

Table 1: Simulation results by GCV using rescaled bandwidth h kck and constraint β̂.
g1
Initial
True

Linear

Equal

Random

g2

g3

n

RSE

RASE

RASE2

RSE

RASE

RASE2

RSE

RASE

RASE2

100

0.2887

0.0959

9.4156

0.2785

0.0945

7.2648

0.3442

0.0918

5.0951

1000

0.0829

0.0336

0.9003

0.0731

0.0334

0.5429

0.1095

0.0276

0.8620

100

1.1401

0.3401

1.3869

1.3569

0.4166

1.2014

0.7049

0.1393

0.2507

1000

1.9010

0.2340

0.5205

1.0719

0.3421

1.3075

0.5287

0.0891

0.0729

100

0.8389

0.2820

1.0933

0.8342

0.2566

1.4349

0.8329

0.1857

0.2589

1000

0.8502

0.3159

1.0133

0.8437

0.2838

1.4043

0.8508

0.1921

0.0905

100

1.1911

0.1461

0.5597

1.1617

0.1979

0.6186

1.2603

0.0942

0.2043

1000

1.2704

0.1320

0.4354

1.3021

0.1890

0.4413

1.2300

0.0796

0.0637

We observe that RSE and RASE, as expected, achieve best performance when we initialize our
optimizer at the true values. However, the more natural Linear initialization strategy does not
outperform initializing at Equal coefficients. For both RASE and RASE2, a more intensive search
over initializations pays off; in the case of RASE2 this even outperforms starting from true values.
We suspect that this is associated with differing optimal smoothness criteria. Table 2 compares
5

our results to when we do not use the h 7 re-scaling where applying this has a significant effect; by
starting from a position far from the optimum, a large bandwidth may have the effect of smoothing
the objective function, indirectly improving our re-estimate of g”.
Overall, the random initial produces a much better results compared to other initial strategies.
Starting from 1000 initial vectors, we have a chance to select the best 10 initial vectors, which
increases the probability of selecting a “good” starting points and decreases the chance of converging
to a local minimum for the non-linear optimization problem, although this comes at a significant
computational prices. In addition, we do observe an improvement of RASE and RASE2 when
n = 1000 compared to n = 100.

22

Table 2: Simulation results of random starting value with rescaled and original bandwidths selected
by GCV

g1

g2

g3

n

original

rescaled

original

rescaled

original

rescaled

100

1.2793

0.5597

2.1497

0.6186

0.2360

0.2043

1000

1.1297

0.4354

2.2081

0.4413

0.0703

0.0637

In order to provide a visual sense of the performance of our estimate, in Figures 1 and 2, we plot
the estimates of the link function g (s) = e−s and g (s) = s, respectively. We observe that the
estimate and the true curve of the link function g almost overlap with each other, but the second
derivative has significantly larger error relative to the truth.

6

Ecological Data

6.1

Model Formulation

Examining our ecological questions, the purpose of estimating the second derivative of the link
function in a functional single index model is to figure out whether the link function g is convex or
concave. Then, we can answer the question: in which environment, constant or varying, the plant
will grow better. We apply our nested estimation method to plant growth dataset. In this dataset,
there are several variables:
1. logarea.t1, logarea.t0 : the plant’s logarithm of area at time t0 and t1, where t0 is the
observation start time and t1 is the end time. A relatively large quantity indicates a high
growth rate of the plant at that time.
2. W: a measure of plant competition. Taken to be a scalar covariate.
3. p.00 − p.36: discrete aggregated temporal record of precipitation, denoted as p (s).
4. t.00 − t.36: discrete aggregated temporal record of temperature, denoted as t (s).
23

8 10

8

g’’

2

4

4

6

6

g

0

2
0
−1

0

⌠Xβ
⌡

1

2

−2

−1

−2

−1

0

1

2

0

1

2

⌠Xβ
⌡

8 10

−2

2

g’’

0

2

−1

4

6

1

^
⌠Xβ
⌡

0

−2
−2

−1

0

⌠Xβ0
⌡

1

2

⌠Xβ0
⌡

Figure 1: Example results using g (s) = e−s . The top-left and right panel are the plots of g and
g′′ over 1000 equally-spaced grid points, while the lower and upper bound are the minimum and
R

R1
1
maximum of 0 X (t) β̂ (t) dt. The bottom-right panel gives g′′ 0 X (t) β̂ (t) dt plotted against
R1
the true 0 X (t) β 0 (t) dt, providing a visual representation of the error we control in Theorem 6.
In all plots, the true curve is given by the black line, while the red line is the estimated curve. The
R1
R1
bottom-left panel is the plot of 0 X (t) β̂ (t) dt versus 0 X (t) β 0 (t) dt, while the red line is y = x.

24

1

g

0

g’’

−0.04

−1

0.00

0.04

2

−2
−2

−1

0

⌠Xβ
⌡

1

2

−2

−1

−2

−1

0

1

2

0

1

2

⌠Xβ
⌡

1

^
⌠Xβ
⌡

0

g’’

−0.04

−1

0.00

0.04

2

−2
−2

−1

0

⌠Xβ0
⌡

1

2

⌠Xβ0
⌡

Figure 2: Example results using g (s) = s. The top-left and right panel are the plots of g and
g′′ over 1000 equally-spaced grid points, while the lower and upper bound are the minimum and
R

R1
1
maximum of 0 X (t) β̂ (t) dt. The bottom-right panel gives g′′ 0 X (t) β̂ (t) dt plotted against
R1
the true 0 X (t) β 0 (t) dt, providing a visual representation of the error we control in Theorem 6.
In all plots, the true curve is given by the black line,, while the red line is the estimated curve. The
R1
R1
bottom-left panel is the plot of 0 X (t) β̂ (t) dt versus 0 X (t) β 0 (t) dt, while the red line is y = x.

25

The precipitation and temperature histories are modeled as two covariate functions. Assume that
the response variable is logarea.t1 − logarea.t0, a Functional Single Index model is:


logarea.t1 − logarea.t0 = g α · W +

Z

pβ1 +

Z



tβ2 ,

where the coefficient α, the functions g, β1 and β2 need to be estimated.

6.2

Results

We used three different starting values: linear, equal and random. For each starting point, we
selected the curve with minimum GCV value. Then, we selected the estimate with the minimum
GCV value among all starting points. The random starting point is selected. The plot of the
estimated g, g ′′ , and the coefficient functions β1 and β2 is in the Figure 3. Since the estimated g′′
is always negative, the link function g is concave. We could conclude that the species will grow
better with a constant environment.

7

Conclusion

To answer the ecological question, we need to figure out the convexity or concavity of the link
function g, or equivalently, find whether the second derivative is positive or negative. In this paper,
we used the local quadratic method to approximate the link function g, and estimated the curvature
of g and the coefficient function β 0 by a nested optimization procedure. Under some assumptions,
we showed that the coefficient function β 0 could be estimated root-n consistently. In addition, the

h R


n
i2
R
P
rate of convergence of the curvature g′′ is n1
E ĝ ′′
Xi β̂ − g ′′
Xi β 0
= O h4n + nh1 4 .
n

i=1

In the simulation study, we used three different link functions, convex, concave and neither convex

nor concave. While we derive convergence rates for the curvature of g, our simulation results
demonstrate the numerical challenges that accompany Functional Single Index models. We can
estimate g fairly well, but our estimates of g′′ are sensitive to the choice of initial condition, requiring
considerable care in optimization. We expect that these numerical challenges are specific to the
estimators employed, but they suggest that alternative means for the influence of environmental
variability on plant growth are warranted.
26

0.0

g −0.5

g’’

−1.0
−0.8

−0.4

0.0

−1.5
−2.0
−2.5
−3.0

0.4

−0.8

Environment

β1

10
5
0
−5
−10
0.2

0.0

0.4

Environment

β2

−0.2

−0.4

2
0
−2
−4

0.6

−0.4

I

0.0

0.4

I

Figure 3: The top-left panel is the estimate of g, while the top-right panel is the estimate of g ′′ .
The bottom panel are the plots of the estimated β1 and β2 .

27

A

Bernstein’s Inequality

Bernstein’s Inequality. Let Y1n , · · · , Ynn be independent random variables with means 0 and
2 for the variance of Y . Suppose V ≥ σ 2 + · · · + σ 2 .
bounded ranges, that is |Yin | ≤ cn . Write σin
in
n
nn
1n

Then for each ηn > 0,
#
ηn2
 .
P (|Y1n + · · · + Ynn | > ηn ) ≤ exp −
2 Vn + 31 cn ηn
"

B

Lemmas for Theorem 6
3+

Lemma 7. If

nh6n

→ ∞,

nh8n


E ĝ

→ 0 and

′′

Z

Xj β̂

3

nhn m−1
− log hn



− ḡ

′′

→ ∞, then we have

Z

Xj β̂

2

=O



1
nh4n



.

Proof. Denote ǫ = (ǫ1 , · · · , ǫn )⊤ . The term could be expressed as


E ĝ

′′

Z

Xj β̂



− ḡ

′′

Z

Xj β̂

2

h 
 i2
= E S2 β̂; j ǫ

R
R

 2
 P
R
 P
2  R
n
n R
R
Xi β̂− Xj β̂
Xi β̂− Xj β̂
2 β̂
ǫ
−
2T
ǫi
K
2Tj0 β̂
Xi β̂ − Xj β̂ K
i
j
hn
hn


i=1
i=1

= E
      2


2
4
0
Tj β̂ Tj β̂ − Tj β̂


R
R
 P
 2
 2
 P
2  R
R
n R
n
R
Xi β̂− Xj β̂
Xi β̂− Xj β̂
2 β̂
2Tj0 β̂
ǫ
ǫi
2T
Xi β̂ − Xj β̂ K
K
i
hn
hn



 j
i=1
i=1



≤ E
+ E
      2
      2 

 .
Tj0 β̂ Tj4 β̂ − Tj2 β̂
Tj0 β̂ Tj4 β̂ − Tj2 β̂


28

We bound the next two terms as:
  P
R
 2
2  R
n R
R
Xi β̂− Xj β̂
0
ǫi
T β̂
Xi β̂ − Xj β̂ K
hn

 j
i=1


E
      2

2
4
0
Tj β̂ Tj β̂ − Tj β̂


1

 1

= E 5 ·
 hn

T0
nhn f (u|β̂ ) j



1
nhn f (u|β̂ )

 
β̂ ·

1
n

 
Tj0 β̂ ·

n R
P

i=1

Xi β̂ −

1
nh5n f (u|β̂ )

R

Xj β̂

2

  
4
Tj β̂ −

2

K

R

R
 2
Xi β̂− Xj β̂
ǫi 
hn

1
nh3n f (u|β̂ )


  2 

Tj2 β̂

 
Z

 2

1
1
1
1
z−u
2



=
· 
·
·
f z β̂ dz · E ǫ2
(z − u) K
2
5
hn µ4 (K) − (µ2 (K)) f u β̂
n
hn
Z
2



1
.
2


(hn m) K (m) f u + hn m β̂ hn dm · E ǫ2
= M1 ·
2 u β̂
nh10
n f

2

Z

1
3
2
 hn f u β̂

∼ M1 ·
m K (m) dm · E ǫ2
2 u β̂
nh10
n f

∼

1
,
nh4n

and
R
 2
 P
R
n
Xi β̂− Xj β̂
ǫi
Tj2 β̂
K
hn


i=1

E
 0   4    2  2 
Tj β̂ Tj β̂ − Tj β̂





 1
= E
 h3 ·
n


1
T2
nh3n (u) j

1
nhn f (u|β̂ )

 
β̂ ·

 
Tj0 β̂ ·

1
n

n
P

K

i=1

1
nh5n f (u|β̂ )

R

R

Xi β̂− Xj β̂
ǫi
hn

2


 2 
  

Tj4 β̂ − nh31f (u) Tj2 β̂
n

2


 
Z
 2

z−u
µ2 (K)
1
1
1



·
K
f z β̂ dz · E ǫ2
=
· 
·
2
3
hn µ4 (K) − (µ2 (K)) f u β̂
n
hn
2
Z



1
.


= M2 ·
K (m) f u + hn m β̂ hn dm · E ǫ2
nh6n f 2 u β̂

2
Z


1

 hn f u β̂
∼ M2 ·
K (m) dm · E ǫ2
nh6n f 2 u β̂
∼

1
.
nh4n

where M1 , M2 are constants and u =

R

Xj β̂.
29

Therefore, we have


E ĝ

′′

Z

Xj β̂



− ḡ

′′

3+

Lemma 8. If nh6n → ∞, nh8n → 0 and


E ḡ

′′

Z

Xj β̂



Z

−g

=O



1
nh4n



.

3

nhn m−1
− log hn
′′

Xj β̂

2

Z

→ ∞, then we have
Xj β̂

2

=O



h4n

1
+
nh4n



.

Proof. Calculate
ḡ ′′

Write

Z




Xj β̂

Z


′′
= S2 β̂; j g − g
Xj β̂
Xj β̂

− g′′

Z

R
 P
2  R
 R
n R

R
Xi β̂− Xj β̂
2Tj0 β̂
Xi β̂ − Xj β̂ K
g Xi β 0
hn
i=1
=
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
R
 P
R
 R
n

Xi β̂− Xj β̂

Z
−2Tj2 β̂
K
g Xi β 0
hn
i=1
′′
−g
Xj β̂ .
−
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂

g

Z

Xi β

0



=g

Z

Xi β̂




 Z

Z
0
Xi β − g
Xi β̂
,
+ g

for i = 1, · · · , n, we decompose it into two terms A and B, where A is
R
 P
2  R
 R

n R
R
Xi β̂− Xj β̂
2Tj0 β̂
Xi β̂ − Xj β̂ K
g
X
β̂
i
hn
i=1

−

      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
R
 R

 P
R
n
Xi β̂− Xj β̂
g
X
β̂
2Tj2 β̂
K
i
hn
i=1

      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂

−g

′′

Z



Xj β̂ ,

and B is

R
R
  R

  R
R
  R

2  R
n

R
P
Xi β̂− Xj β̂
Xi β̂− Xj β̂
0
2
0
2Tj β̂
− 2Tj β̂ K
g Xi β − g
Xi β̂
Xi β̂ − Xj β̂ K
hn
hn
i=1

Tj0

      2
β̂ Tj4 β̂ − Tj2 β̂
30

.

.
Denote Kij = K

R

R

Xi β̂− Xj β̂
,
hn

the term A can be calculated as

 P
2
R

 P
R

n R
n
R

Z
2Tj0 β̂
Xi β̂ − Xj β̂ Kij g
Xi β̂ − 2Tj2 β̂
Kij g
Xi β̂
i=1
i=1
′′
−g
Xj β̂
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
 P
 P
R

2
R

n
n R
R

Z
2Tj2 β̂
2Tj0 β̂
Kij g
Xi β̂
Xi β̂ − Xj β̂ Kij g
Xi β̂
i=1
i=1
′′
−
Xj β̂ ,
=
      2
      2 − g
Tj0 β̂ Tj4 β̂ − Tj2 β̂
Tj0 β̂ Tj4 β̂ − Tj2 β̂

where
R
 P
2  R
 R

n R
R
Xi β̂− Xj β̂
2Tj0 β̂
Xi β̂ − Xj β̂ K
g
X
β̂
i
hn
i=1
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂

 P
2
R
2
R
R
4 

n R
1 ′′
0
′
2Tj β̂
Xi β̂ − u Kij g (u) + g (u)
Xi β̂ − u + O
Xi β̂ − u
Xi β̂ − u + 2 g (u)
i=1
=
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
   
 

 
2Tj0 β̂ Tj2 β̂ g (u) + Tj3 β̂ g ′ (u) + 12 Tj4 β̂ g′′ (u)

+ O h2n
=
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
    R

    R

2Tj0 β̂ Tj2 β̂ g
Xj β̂ + Tj0 β̂ Tj4 β̂ g′′
Xj β̂

=
+ O h2n ,
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
and

R
 P
 R

R
n
Xi β̂− Xj β̂
2Tj2 β̂
g
X
β̂
K
i
hn
i=1
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂

R
2
R
R

4 
 P
R
n
Xi β̂−u
1 ′′
′
2
Xi β̂ − u + O
g (u) + g (u)
Xi β̂ − u + 2 g (u)
Xi β̂ − u
2Tj β̂
K
hn
i=1
=
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
   
 
 

2Tj2 β̂ Tj0 β̂ g (u) + Tj1 β̂ g ′ (u) + 12 Tj2 β̂ g′′ (u)

=
+ O h2n
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
    R
   2 R

2Tj2 β̂ Tj0 β̂ g
Xj β̂ + Tj2 β̂
Xj β̂
g′′

+ O h2n .
=
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂

31

Therefore, the term A is
    R

    R

2Tj0 β̂ Tj2 β̂ g
Xj β̂ + Tj0 β̂ Tj4 β̂ g ′′
Xj β̂
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
    R
   2 R

Z

2Tj2 β̂ Tj0 β̂ g
Xj β̂ + Tj2 β̂
g′′
Xj β̂

′′
−g
−
Xj β̂ + O h2n
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂


Z
Z

′′
′′
=g
Xj β̂ − g
Xj β̂ + O h2n

= O h2n .

The term B is bounded by

R
R

  R
 
  R

2  R
n
R
P
Xi β̂− Xj β̂
Xi β̂− Xj β̂
1
2
0
√
−
2T
E
β̂
K
2Tj β̂
Xi β̂ − Xj β̂ K
1 n
j
hn
hn
|B| ≤ i=1
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂


   
    
2Tj0 β̂ Tj2 β̂ − 2Tj0 β̂ Tj2 β̂
E1 √1n
=
      2
Tj0 β̂ Tj4 β̂ − Tj2 β̂
1
∼√ 2 ,
nhn

where E1 is a constant.
Combining the terms A and B, we have


 Z

2
Z
1
4
′′
′′
.
= O hn +
E ḡ
Xj β̂ − g
Xj β̂
nh4n

3+

Lemma 9. If nh6n → ∞, nh8n → 0 and

3

nhn m−1
− log hn

→ ∞, then we have

 Z
 

Z
2
1
E g′′
=O
.
Xj β̂ − g′′
Xj β 0
n

 
R
R
Proof. By Theorem 5, we have
Xi β̂ − Xi β 0 = O √1n . Since g ′′ satisfies the Lipschitz
condition, we can get
 Z

Z
2
Z
 
2
Z
1
′′
′′
0
0
E g
,
Xj β̂ − g
Xj β
≤ K3
Xi β̂ − Xi β
=O
n

where K3 is a constant.
32

Table 4: Simulation results of random starting value with rescaled and original bandwidths selected
by 10-fold cross-validation

g1

C

g2

g3

n

original

rescaled

original

rescaled

original

rescaled

100

1.2710

0.5191

2.2005

0.7160

0.2768

0.2362

1000

1.0145

0.4209

2.2812

0.4418

0.0752

0.0744

10-fold CV Results

The simulation results using the 10-fold cross-validation is:
Table 3: Simulation results by 10-fold CV using rescaled bandwidth h kck and constraint β̂.
g1
Initial
True

Linear

Equal

Random

g2

g3

n

RSE

RASE

RASE2

RSE

RASE

RASE2

RSE

RASE

RASE2

100

0.1921

0.0868

1.2106

0.2052

0.0904

4.5719

0.1982

0.0690

0.7954

1000

0.0680

0.0290

0.4551

0.0620

0.0271

0.4016

0.0774

0.0243

0.2762

100

1.1383

0.3294

1.4303

1.3370

0.4245

1.3388

0.6995

0.1331

0.2692

1000

1.9015

0.2232

0.5060

1.0419

0.3185

1.3652

0.5297

0.0852

0.0786

100

0.8486

0.2789

1.0731

0.8385

0.2509

1.4728

0.8349

0.1836

0.3349

1000

0.8442

0.2977

1.0181

0.8500

0.2755

1.5396

0.8447

0.1827

0.1045

100

1.2998

0.1339

0.5191

1.2086

0.1979

0.7160

1.4296

0.0875

0.2362

1000

1.5607

0.1288

0.4209

1.3540

0.1687

0.4418

1.7622

0.0798

0.0744

For comparison purpose, we estimate the curvature with or without rescaling the bandwidth, where
the results are in Table 4.

33

Table 5: Simulation results of true starting value with rescaled and original bandwidths selected
by 10-fold cross-validation

g1

D

g2

g3

n

original

rescaled

original

rescaled

original

rescaled

100

5.9302

1.9266

33.811

7.9015

4.5630

0.9791

1000

2.4004

1.0660

2.2494

1.3416

1.2444

0.3558

True Starting Value

In Table 5 and 6, we compare the RASE2 results of true starting value with rescaled and original
bandwidths.
Table 6: Simulation results of true starting value with rescaled and original bandwidths selected
by GCV

g1

E

g2

g3

n

original

rescaled

original

rescaled

original

rescaled

100

51.397

7.3786

38.007

6.5712

50.816

5.9508

1000

9.0220

1.5416

8.9590

1.8170

7.8245

1.1493

CV Values

The CV values for both GCV and 10-fold cross-validation are in Table 7.

34

Table 7: CV values

g1
Initial
True

Linear

Equal

Random

g2

g3

n

GCV

10-fold

GCV

10-fold

GCV

10-fold

100

0.0398

3.8686

0.0399

3.3823

0.0377

1.4970

1000

0.0401

0.9468

0.0402

0.7882

0.0400

0.4781

100

0.1559

5.2259

0.2184

4.7373

0.0498

1.6148

1000

0.0793

2.0295

0.1272

2.2896

0.0439

0.7265

100

0.1401

4.7640

0.1046

3.9670

0.0795

1.8530

1000

0.1358

2.6925

0.1318

2.1995

0.0751

1.1081

100

0.0611

4.2920

0.0957

3.9025

0.0459

1.6237

1000

0.0591

2.4463

0.0850

2.2702

0.0468

1.0243

35

References
Chen, D., P. Hall, H.-G. Müller, et al. (2011). Single and multiple index functional regression
models with nonparametric link. The Annals of Statistics 39 (3), 1720–1747.
Chesson, P. L. and R. R. Warner (1981). Environmental variability promotes coexistence in lottery
competitive systems. The American Naturalist 117 (6), 923–943.
Clark, J. S., D. M. Bell, M. H. Hersh, and L. Nichols (2011). Climate change vulnerability of
forest biodiversity: climate and competition tracking of demographic rates. Global Change Biology 17 (5), 1834–1849.
Dahlgren, J. P. and J. Ehrlén (2011). Incorporating environmental change over succession in an
integral projection model of population dynamics of a forest herb. Oikos 120 (8), 1183–1190.
Drake, J. M. (2005). Population effects of increased climate variation. Proceedings of the Royal
Society of London B: Biological Sciences 272 (1574), 1823–1827.
Ellner, S. (1987). Alternate plant life history strategies and coexistence in randomly varying environments. In Theory and models in vegetation science, pp. 199–208. Springer.
Escabias, M., A. M. Aguilera, and M. J. Valderrama (2007). Functional pls logit regression model.
Computational Statistics & Data Analysis 51 (10), 4891–4902.
Hardle, W., P. Hall, H. Ichimura, et al. (1993). Optimal smoothing in single-index models. The
annals of Statistics 21 (1), 157–178.
Hristache, M., A. Juditsky, and V. Spokoiny (2001). Direct estimation of the index coefficient in a
single-index model. Annals of Statistics, 595–623.
Hutchinson, G. E. (1961). The paradox of the plankton. The American Naturalist 95 (882), 137–145.
Ichimura, H. (1993). Semiparametric least squares (sls) and weighted sls estimation of single-index
models. Journal of Econometrics 58 (1-2), 71–120.

36

James, G. M. (2002). Generalized linear models with functional predictors. Journal of the Royal
Statistical Society: Series B (Statistical Methodology) 64 (3), 411–432.
Koons, D. N., S. Pavard, A. Baudisch, J. E. Metcalf, et al. (2009). Is life-history buffering or lability
adaptive in stochastic environments? Oikos 118 (7), 972–980.
Lewontin, R. C. and D. Cohen (1969). On population growth in a randomly varying environment.
Proceedings of the National Academy of Sciences 62 (4), 1056–1060.
Li, K.-C. (1991). Sliced inverse regression for dimension reduction. Journal of the American
Statistical Association 86 (414), 316–327.
Ma, S. (2016). Estimation and inference in functional single-index models. Annals of the Institute
of Statistical Mathematics 68 (1), 181–208.
Müller, H.-G. and U. Stadtmüller (2005). Generalized functional linear models. Annals of Statistics,
774–805.
Stoker, T. M. (1986). Consistent estimation of scaled coefficients. Econometrica: Journal of the
Econometric Society, 1461–1481.
Teller, B. J., P. B. Adler, C. B. Edwards, G. Hooker, and S. P. Ellner (2016). Linking demography
with drivers: climate and competition. Methods in Ecology and Evolution 7 (2), 171–183.

37

