Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2915777

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Model and Data Structure
- 3 Model Estimation With Fixed Basis Functions
- 4 Estimation With Principal Components
- 5 Numerical Studies
- 6 Conclusions
- Supplementary Material
- References

Authors

Related links

J Am Stat Assoc. Author manuscript; available in PMC 2010 August 4.

Published in final edited form as:

J Am Stat Assoc. 2010 June 1; 105(490): 621–633.

doi: 10.1198/jasa.2010.tm09313PMCID: PMC2915777

NIHMSID: NIHMS184176

Yehua Li, Department of Statistics, University of Georgia, Athens, GA 30602, Email: ude.agu@ilauhey;.

See other articles in PMC that cite the published article.

We introduce a new class of functional generalized linear models, where the response is a scalar and some of the covariates are functional. We assume that the response depends on multiple covariates, a finite number of latent features in the functional predictor, and interaction between the two. To achieve parsimony, the interaction between the multiple covariates and the functional predictor is modeled semiparametrically with a single-index structure. We propose a two step estimation procedure based on local estimating equations, and investigate two situations: (a) when the basis functions are pre-determined, e.g., Fourier or wavelet basis functions and the functional features of interest are known; and (b) when the basis functions are data driven, such as with functional principal components. Asymptotic properties are developed. Notably, we show that when the functional features are data driven, the parameter estimates have an increased asymptotic variance, due to the estimation error of the basis functions. Our methods are illustrated with a simulation study and applied to an empirical data set, where a previously unknown interaction is detected. Technical proofs of our theoretical results are provided in the online supplemental materials.

In this paper, we are interested in functional data analysis for regression problems with a scalar response and where some of the covariates are functional. The most important model for such problems is the functional generalized linear model, where the response depends on a linear projection of the functional predictor. Some recent work on this problem includes Ramsay and Silverman (2005), Müller and Stadtmüller (2005), Yao, Müller and Wang (2005), Cardot and Sarda (2005), Cai and Hall (2006), Ferraty and Vieu (2006), Crainiceanu, Staicu and Di (2008) and Di, Crainiceanu, Caffo and Punjabi (2009), among many others.

However, most existing work does not readily accommodate the existence of an interaction between the functional predictor and other covariates. Our paper addresses this issue by introducing a novel class of models that accommodate such an interaction. A special feature of our work, which differs from the existing literature, is that we use a semiparametric single index structure to model the potential interaction between the functional predictor and other covariates. This single index structure not only gives flexibility in describing a complex interaction when there are two or more other covariates but it also allows us to propose a practically feasible two-stage procedure to address estimation and inference issues.

Our new *generalized functional linear models with semiparametric interactions* have the structure that the response depends on latent features of the functional data and their interaction with other possibly multivariate covariates. The “features” in functional data that we are considering are the projections of the functional data onto orthonormal basis functions. These basis functions can be either fixed, e.g., Fourier or wavelet basis functions, or data driven, e.g., principal components. The interaction is modeled semiparametrically, with the regression coefficient function for the functional predictor depending on a single-index function of the multivariate covariates.

After defining the modeling framework in Section 2, we split our study into two tracks. In the first track, in Section 3, we consider the situation that the functional predictor is fully observed and the basis functions are pre-determined, as might occur for Fourier basis functions or some versions of splines. Consequently, the score of the latent features can be evaluated. In this context we propose a backfitting estimation procedure based on local estimating equations. The procedure we propose consists of two stages, which are based on the philosophy of minimum average variance estimation (MAVE, Xia, Tong, Li and Zhu, 2002): in the first stage, we use multivariate kernel weights in the local estimating equations to get consistent initial values; in the second stage, we switch to univariate kernel methods to achieve more efficient estimators. We also derive the asymptotic properties of the methodology.

The second track of the paper, in Section 4, considers the case that the functional features/basis functions are data driven and need to be estimated. At present, the most common dimension reduction device in functional data analysis is functional principal component analysis (FPCA). Some recent work, including Yao, et al. (2005), Hall and Hosseini-Nasab (2006) and Hall, Müller and Wang (2006), have led to a deeper understanding of this method and its properties. We apply the FPCA method based on kernel smoothing as proposed by the authors mentioned above, and then plug the estimated PCA scores into the two stage estimation procedure of Section 3. With the assumption that the number of observations per curve goes to infinity at a sufficient rate, we show that the estimators of the proposed models are still root-n consistent and asymptotically normally distributed. An important finding is that, even when there is sufficient number of observations per curve, the estimation error in FPCA will increase the variance of the final estimator of the functional linear model. This fact is not sufficiently well appreciated in the literature.

We illustrate the numerical performance of our procedure in a simulation study given in Section 5, and by an empirical application to colon carcinogenesis data, where we exhibit a previously unknown interaction. Final remarks are provided in Section 6, and all proofs are sketched in the Web Appendix, which is available in the online supplemental materials.

The data consist of *i* = 1, …, *n* independent triples (*Y _{i}, X_{i}, Z_{i}*), where

$$\begin{array}{c}g\{{\mu}_{Y}(X,\mathit{Z})\}={\displaystyle {\int}_{\mathcal{T}}\U0001d504(t,{\mathit{Z}}_{1})X(t)\mathit{\text{dt}}+{\mathit{\beta}}^{\mathrm{T}}}\mathit{Z};\hfill \\ \text{var}(Y|X,\mathit{Z})={\sigma}_{Y}^{2}V\{{\mu}_{Y}(X,\mathit{Z})\},\hfill \end{array}$$

(1)

where *g*(·) and *V*(·) are known functions, (·,**Z**_{1}) and β are unknown and = [*a, b*] is a fixed interval. We allow the coefficient function (·, ·) to depend on a covariate vector **Z**_{1} which is a subset of * Z*.

Suppose ψ_{1}(*t*), ψ_{2}(*t*), , ψ_{p}(*t*) are *p* orthonormal functions on , and ξ* _{j}* = ∫ ψ

By assumption, *Y* depends on *X* (·) only through the features **ξ**. Let ψ(*t*) = {ψ_{1}(*t*), , ψ_{p}(*t*)}^{T}. To make the model parsimonious without imposing strong parametric structural assumptions, and to allow interactions between *X* and **Z**_{1}, we further assume that

$$\U0001d504(\xb7,{\mathit{Z}}_{1})={\mathit{\psi}}^{\mathrm{T}}(t)\mathit{\alpha}({\mathit{Z}}_{1};\mathit{\theta}),\text{and that}\mathit{\alpha}({\mathit{Z}}_{1};\mathit{\theta})={\mathit{\alpha}}_{1}+\mathit{S}({\mathit{Z}}_{1};\mathit{\theta}){\mathit{\alpha}}_{2},$$

(2)

where * S*(

For interpretation of such interaction, we focus on a given ψ* _{j}*(

Our model is innovative and flexible in the following ways. First, when * S*(

$${g}^{-1}\{\mathrm{E}(Y|\mathit{\xi},\mathit{Z})\}={\mathit{\alpha}}^{\mathrm{T}}({\mathit{Z}}_{1};\mathit{\theta})\mathit{\xi}+{\mathit{\beta}}^{\mathrm{T}}\mathit{Z}.$$

(3)

Therefore, our semiparametric interaction model is readily applicable to multivariate data, where **ξ** is an observed covariate vector.

In this section, we consider the case that the basis functions (ψ_{j}) are known, *X _{i}*(

We study this seemingly unrealistic scenario first for two reasons. First, this is a commonly used ideal situation typically considered first in the functional data literature in order to motivate new methodology. A more realistic situation will be studied in our Section 4. Second, as pointed out in the previous section, the semiparametric interaction model we proposed is not limited to functional data. The methods we propose below for this ideal situation are also applicable to the multivariate semiparametric model (3) when **ξ** is another observed multivariate covariate.

We estimate * S*(·) through a local polynomial smoothing approach. Let $\mathbf{\Theta}={({\mathit{\alpha}}_{1}^{\mathrm{T}},{\mathit{\alpha}}_{2}^{\mathrm{T}},{\mathit{\beta}}^{\mathrm{T}},{\mathit{\theta}}^{\mathrm{T}})}^{\mathrm{T}}$. Our strategy is to iterate between estimating

Let *Q*(*w, y*) be the quasi-likelihood function satisfying *Q*(*w, y*)/*w* = (*y* − *w*)/*V* (*w*) (McCullagh and Nelder, 1989). For a given value of **Θ**, we estimate * S*(υ) and

$${\sum}_{i=1}^{n}Q\left({g}^{-1}[{\mathit{\alpha}}_{1}^{\mathrm{T}}{\mathit{\xi}}_{i}+\{{a}_{0}+{a}_{1}({\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}-\upsilon )\}{\mathit{\alpha}}_{2}^{\mathrm{T}}{\mathit{\xi}}_{i}+{\mathit{\beta}}^{\mathrm{T}}{\mathit{Z}}_{i}],{Y}_{i}\right){w}_{i}(\upsilon ),$$

(4)

where *w _{i}*(υ) is the local weight for the

Model (1) is new, although in various special cases it shares the spirit of previous work on single index models such as Carroll, et al. (1997), Xia et al. (2002) and Xia and Härdle (2006). For computation, we adopt the philosophy of the latter two papers, describing a version of their MAVE method applicable to our problem. Consider υ = **θ**^{T} **Z**_{j1} in (4) for each *j* {1, , *n*} and let *a*_{0j} and *a*_{1j} be the estimate of * S*(

Summing over υ = **θ**^{T} **Z**_{j1} for all *j* in (4) while holding the values of (*a*_{0j}, *a*_{1j}) fixed, and denoting **Z**_{i1} − **Z**_{j1} by **Z**_{ij,1}, we can now update the rest of the parameters by minimizing

$$\mathcal{Q}={\displaystyle {\sum}_{j=1}^{n}{\displaystyle {\sum}_{i=1}^{n}Q\left[{g}^{-1}\{{\mathit{\alpha}}_{1}^{\mathrm{T}}{\mathit{\xi}}_{i}+({a}_{0j}+{a}_{1j}{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{\mathit{\text{ij}},1}){\mathit{\alpha}}_{2}^{\mathrm{T}}{\mathit{\xi}}_{i}+{\mathit{\beta}}^{\mathrm{T}}{\mathit{Z}}_{i}\},{Y}_{i}\right]}{w}_{\mathit{\text{ij}}},}$$

(5)

with respect to **Θ**, where *w*_{ij} are local weights specified in more detail below. As described above, for identifiability purposes we let **θ** and **α**_{2} satisfy ‖**θ**‖ = 1 and ‖**α**_{2}‖ = 1. The constraint E{* S*(

We have now restructured the estimation procedure into one that minimizes the local quasi-likelihood function (5) iteratively with respect to {*a*_{0j}, *a*_{1j}, *j* = 1, , *n*} and **Θ**. One remaining task is to specify the local weights *w*_{ij}. We use two sets of weights. In the initial stage, we let ${w}_{\mathit{\text{ij}}}^{[1]}={H}_{b}({\mathit{Z}}_{\mathit{\text{ij}},1})/{\displaystyle {\sum}_{\ell =1}^{n}{H}_{b}({\mathit{Z}}_{\ell j,1})}$, where *H*(·) is a *d*_{1}-dimensional kernel function and *H _{b}*(·) =

The estimation procedure then is as follows.

**Stage 1** (Initial estimator) Iteratively update * S*(·) and

**Stage 2** (Refined estimator) To improve efficiency, we replace the initial weights by ${w}_{\mathit{\text{ij}}}={w}_{\mathit{\text{ij}}}^{[2]}(\tilde{\mathit{\theta}})$, where is the consistent estimator of **θ** from Stage 1. The final estimators are denoted as and * Ŝ*(·).

In common with the MAVE procedure in Xia and Härdle (2006), the procedure described above does not require a root-n consistent pilot estimator of **Θ**, which is difficult to obtain in the functional data setting. We will show in Section 3.2 that the initial estimator using the multivariate kernel weights provides a consistent estimator, while the use of the univariate kernel weights with consistently estimated **θ** at Stage 2 yields a more efficient asymptotically normally distributed estimator.

In practice, we use a one-step Newton-Raphson update version of the iterative algorithm in Stages 1 and 2, which speeds up the computation considerably. As shown in our proofs for Theorem 1 and 2, the iterations in each stage generate a self-attraction process; that is, asymptotically, the distance between the current and previous estimates of **Θ** shrinks to zero along with the iteration. This is one of the standard ways to show convergence properties of an iterative algorithm. It is also used by Xia and Härdle (2006) in a setup simpler than ours. In practice, with *n* being finite, an algorithm with this asymptotic convergence property could still fail to converge. We have not encountered such numerical difficulties in our study though. We have provided more details about the algorithm in the web supplement.

Denote the true parameter as **Θ**_{0}, and the true interaction link function as **S**_{0}(·). Following the notation in Carroll et al. (1997), let *q*_{1}(*x, y*) = {*y* − *g*^{−1}(*x*)}ρ_{1}(*x*) and ${q}_{2}(x,y)=\{y-{g}^{-1}(x)\}{\rho}_{1}^{\prime}(x)-{\rho}_{2}(x)$, where ρ_{}(*x*) = {*dg*^{−1}(*x*)/*dx*}^{}/*V*{*g*^{−1}(*x*)}, = 1, 2. As mentioned previously, here we assume the (ξ_{ij}) are known. Estimation of ξ_{ij} and the effect of the estimation error will be discussed in Section 4. We make the following assumptions.

- (C1.1) The marginal density of
**Z**_{1}is positive and uniformly continuous, with a compact support^{d1}. Let $\mathcal{U}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}$ be the true single-index, with a density*f*_{}(*u*). Assume*f*is twice continuously differentiable._{} - (C1.2)
*g*^{(3)}(·) and*V*^{(2)}(·) exist and are continuous. - (C1.3)
**S**_{0}(·) is twice continuously differentiable. - (C1.4) For all
*x*, ρ_{2}(*x*) > 0 and*d*ρ_{2}(*x*)/*dx*is bounded.

- (C2.1)
*H*(·) is a symmetric*d*_{1}dimensional probability density function, ∫*H*(**x**)**xx**^{T}*d***x**=*I*.*b*→ 0, log(*n*)/(*nb*^{d1+2}) → 0. - (C2.2)
**V**_{1}(**z**;**Θ**,),**S****V**_{2}(**z**;**Θ**,) and**S****V**_{3}(**z**;**Θ**,) defined in Web Appendix B.1 are twice continuously differentiable in**S****z**, and Lipschitz continuous in**Θ**.

- (C3.1)
*K*(·) is a symmetric univariate probability density function, with support on [−1, 1], ${\sigma}_{K}^{2}={\displaystyle \int {u}^{2}K(u)\mathit{\text{du}},\phantom{\rule{thinmathspace}{0ex}}{\kappa}_{2}={\displaystyle \int {K}^{2}(u)\mathit{\text{du}}.\text{}h}}~{n}^{-\lambda},1/6\lambda 1/4$. - (C3.2)
_{1}(**z**;**θ**),_{2}(**z**;**θ**) and_{3}(**z**;**θ**) in Web Appendix B.2 are twice continuously differentiable in**z**, and Lipschitz continuous in**θ**.

Here are the two main results in this section.

THEOREM **1** *(Consistency of the initial estimator) Under conditions (C1)–(C3)*, ‖**Θ̃**−**Θ**_{0}‖ → 0 *with probability 1*.

THEOREM **2** *(Asymptotic normality of the refined estimator) Under conditions (C1)–(C3)*,

$$\begin{array}{c}\sqrt{n}(\widehat{\mathbf{\Theta}}-{\mathbf{\Theta}}_{0})\to \text{Normal}(\mathbf{0},{\mathcal{A}}^{-}\mathcal{B}{\mathcal{A}}^{-}),\hfill \\ \sqrt{nh}\{\widehat{\mathit{S}}(u)-{\mathit{S}}_{0}(u)-{h}^{2}{\mathit{S}}_{0}^{(2)}(u){\sigma}_{k}^{2}/2\}\to \mathit{\text{Normal}}\{0,{\sigma}_{\mathit{S}}^{2}(u)\},\hfill \\ \hfill \\ \hfill \\ \hfill \end{array}$$

*where and are defined in the* Appendix, ${\sigma}_{\mathit{S}}^{2}(u)={\kappa}_{2}\mathrm{E}\{({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i}){\in}_{i}^{2}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}=u\}/\{{\sigma}_{\upsilon}^{4}(u){f}_{\mathcal{U}}(u)\},{\sigma}_{\upsilon}^{2}(u)=\mathrm{E}\{{\rho}_{2}({\eta}_{i}){({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})}^{2}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}=u\}$.

**Remark:** The algorithm we use is similar to the MAVE method (Xia et al. 2002), which is a relatively new development in the literature on backfitting methods. Unlike many other back-fitting algorithms, the method we are using does not require any artificial undersmoothing for the nonparametric function * S*(·), see our condition (C3.1) on the rate of the bandwidth. The method also differs from the algorithm in Carroll et al. (1997) in that we do not need a root-

The most popular dimension reduction method in functional data analysis is principal component analysis, leading to estimated basis functions. Here we focus on the common scenario in which ψ_{j}(*t*) denotes the eigenfunction in functional principal component analysis. For the longitudinal process {*X*(*t*), *t* } with mean function μ_{X} (*t*), the covariance function is *R*(*s, t*) = cov{*X*(*s*), *X* (*t*)} with eigen-decomposition $R(s,t)={\displaystyle {\sum}_{j=1}^{\infty}{\omega}_{j}{\psi}_{j}(s){\psi}_{j}(t)}$, where the ω_{j} are the non-negative eigenvalues of *R*(·, ·), which, without loss of generality, satisfy ω_{1} > ω_{2} > > 0, and the ψ_{j}’s are the corresponding eigenfunctions. The Karhunen-Loève expansion of *X* (*t*) is

$$X(t)-{\mu}_{X}(t)={\displaystyle {\sum}_{j=1}^{\infty}{\xi}_{j}{\psi}_{j}(t),}$$

(6)

where ξ_{j} = ∫ ψ_{j}(*t*){*X*(*t*) − μ_{X}(*t*)}*dt* has mean zero, with cov(ξ_{j}, ξ_{k}) = *I* (*j* = *k*)ω_{j}.

Under this framework, our model (1) means that *Y* is dependent on the leading *p* principal components in *X*(·). The reason that we assume that only a finite number of PC’s are relevant is that, in functional data that is commonly seen in biological applications, estimation of high order PC’s is highly unstable and difficult to interpret, see the comments in Rice and Silverman (1991) and Hall, et al. (2006). Yao et al. (2005) proposed to use an AIC criterion to choose p. Even though, to the best of our knowledge, there is no theoretical support behind the use of AIC in the existing literature, we found it a very sensible criterion in our numerical investigations which we have reported in Section 5. We made a slight modification of their method in terms of counting the number of parameters, and at least in our simulations were able to select the correct number of PC in every case. More work on this general topic clearly needs to be done.

It often happens that the covariate process we observe contains additional random errors and instead we observe

$${W}_{\mathit{\text{ij}}}={X}_{i}({t}_{\mathit{\text{ij}}})+{U}_{\mathit{\text{ij}}},\text{}j=1,\cdots ,{m}_{i},$$

(7)

where *U _{ij}* are independent zero-mean errors, with $\text{var}\{{U}_{i}(t)\}={\sigma}_{u}^{2}$, and the

Many functions and parameters in the expressions given above can be estimated from the data. Thus, μ* _{X}*(·) and

$$({\widehat{a}}_{0},{\widehat{a}}_{1})={\text{argmin}}_{a0,a1}{\displaystyle {\sum}_{i=1}^{n}{\displaystyle {\sum}_{j=1}^{{m}_{i}}{\{{W}_{\mathit{\text{ij}}}-{a}_{0}-{a}_{1}({t}_{\mathit{\text{ij}}}-t)\}}^{2}K\{{t}_{\mathit{\text{ij}}}-t)/{h}_{\mu}\},}}$$

*K*(·) is a kernel function and *h*_{μ} is the bandwidth for estimating μ_{X}. Letting _{XX}(*s, t*) = E{*X*(*s*)*X*(*t*)}, define ϖ̂_{XX}(*s, t*) = *â*_{0}, where (*â*_{0}, *â*_{1}, *â*_{2}) minimizes

$${\sum}_{i=1}^{n}{\displaystyle {\sum}_{j=1}^{{m}_{i}}{\displaystyle {\sum}_{k\ne j}{\{{W}_{\mathit{\text{ij}}}{W}_{\mathit{\text{ik}}}-{a}_{0}-{a}_{1}({t}_{\mathit{\text{ij}}}-s)-{a}_{2}({t}_{\mathit{\text{ik}}}-t)\}}^{2}}K\left(\frac{{t}_{\mathit{\text{ij}}}-s}{{h}_{\varpi}}\right)K\left(\frac{{t}_{\mathit{\text{ik}}}-t}{{h}_{\varpi}}\right)}.$$

Then (*s, t*) = ϖ̂_{XX}(*s, t*) − _{X}(*s*)_{X}(*t*). The (ω_{k}) and {ψ_{k}(·)} can be estimated from an eigenvalue decomposition of (·, ·) by discretization of the smoothed covariance function, see Rice and Silverman (1991) and Capra and Müller (1997). By realizing that ${\sigma}_{w}^{2}(t)=\text{var}\{W(t)\}=R(t,t)+{\sigma}_{u}^{2}$, we let ${\widehat{\sigma}}_{w}^{2}(t)={\widehat{a}}_{0}-{\widehat{\mu}}_{X}^{2}(t)$, where (*â*_{0}, *â*_{1}) minimizes

$$\sum}_{i=1}^{n}{\displaystyle {\sum}_{j=1}^{{m}_{i}}{\{{W}_{\mathit{\text{ij}}}^{2}-{a}_{0}-{a}_{1}({t}_{\mathit{\text{ij}}}-t)\}}^{2}K\{{t}_{\mathit{\text{ij}}}-t)/{h}_{\sigma}\},$$

and ${\sigma}_{u}^{2}$ is estimated by ${\widehat{\sigma}}_{u}^{2}={(b-a)}^{-1}{\displaystyle {\int}_{a}^{b}\{{\widehat{\sigma}}_{w}^{2}(t)-\widehat{R}(t,t)\}}\mathit{\text{dt}}$.

There are two ways to predict the principal component scores, ξ_{ik}. The first is the PACE method proposed by Yao, et al. (2005). By assuming the covariate process *X _{i}*(

$${\widehat{\xi}}_{\mathit{\text{ik}}}={\widehat{\omega}}_{k}{\widehat{\mathit{\psi}}}_{\mathit{\text{ik}}}{\widehat{\mathrm{\Sigma}}}_{{W}_{i}}^{-1}({\mathit{W}}_{i}-{\widehat{\mathit{\mu}}}_{\mathit{\text{Xi}}}).$$

(8)

An alternative method is by numerical integration. Motivated by the definition that ξ_{ik} = ∫{*X _{i}*(

$${\widehat{\xi}}_{\mathit{\text{ik}}}={\displaystyle \sum _{j=2}^{{m}_{i}}\{{W}_{\mathit{\text{ij}}}-{\widehat{\mu}}_{x}({t}_{\mathit{\text{ij}}})\}{\widehat{\psi}}_{k}({t}_{\mathit{\text{ij}}})({t}_{\mathit{\text{ij}}}-{t}_{i,j-1}).}$$

(9)

In our numerical experience, when we have densely sampled functional data, i.e. *m _{i}* is large and the

In this section, we study the asymptotic properties of the estimators proposed in Section 3 when the principal component scores are replaced by their estimates. Without loss of generality, we assume that μ_{X} (*t*) = 0, and that there are the same number of discrete observations for each subject, i.e. *m _{i}* =

- (C4.1) (Hall, et al., 2006) The process
*X*(·) is independent of the measurement error*U*in (7), E(*U*) = 0 and $\text{var}(U)={\sigma}_{u}^{2}$. For some*c*> 0, ${\text{max}}_{j=0,1,2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{E}\left\{{\text{sup}}_{t\in \mathcal{T}}|{X}^{(j)}(t){|}^{c}\right\}+\mathrm{E}(|U{|}^{c})\infty $. - (C4.2) We have ω
_{1}> ω_{2}> > ω_{p}> ω_{p+1}≥ 0 and ϕ_{k}is twice continuously differentiable on , for*k*= 1, ,*p*. - (C4.3) The
*t*are equally spaced in , i.e._{ij}*t*=_{ij}*a*+ (*j*/*m*) × (*b − a*), for*j*= 1, ,*m. m*→ ∞. - (C4.4)
*h*_{}~*n*^{−λ}, 1/4 < λ_{}< 1/3.

In Assumption (C4.3) we assume the (*t _{ij}*) are fixed and equally spaced, thus allowing convenient mathematical derivations. The conclusion of Lemma 1 holds for random designs. The proofs of the following lemma and theorem are provided in the Web Appendix.

LEMMA **1** *Under Assumption (C4)*, $\Vert {\widehat{\psi}}_{k}-{\psi}_{k}\Vert ={O}_{p}\{{h}_{\varpi}^{2}+{({\mathit{\text{nh}}}_{\varpi})}^{-1/2}\}$, *for k* = 1, , *p. Let _{ik} be the estimated principal component score defined in* (9).

THEOREM **3** *Suppose* **Θ̃** *is the initial estimator of* **Θ** *defined in Section 3, using the estimated principal component scores. Assume that the bandwidth b for the multivariate kernel satisfies mb* → ∞. *Under conditions (C1), (C2) and (C4),* ‖**Θ̃** − **Θ**_{0}‖ → 0 *with probability 1*.

In this subsection we will discuss the scenario under which we can pretend that we know the entire trajectory of *X _{i}*(·). More precisely, we give conditions that ensure that the smoothing errors are asymptotically negligible. Consequently, we can smooth each curve and effectively eliminate effects from the measurement error in (7). Such a “smooth first, then perform estimation” procedure was considered and justified by Zhang and Chen (2007) in the functional linear model setting.

Let _{i} (*t*) be the estimated trajectory of *X _{i}*(

- (C5.1) (
*t*_{1},*t*_{2},*t*_{3},*t*_{4}) = E{*X*(*t*_{1})*X*(*t*_{2})*X*(*t*_{3})*X*(*t*_{4})} −*R*(*t*_{1},*t*_{2})*R*(*t*_{3},*t*_{4}) exists for all*t*._{j} - (C5.2) We have
*m*≥*Cn*, with κ > 5/4.^{κ} - (C5.3) The
*t*are independent identically distributed random design points with density function_{ij}*f*(·), where*f*is bounded away from 0 on and is continuously differentiable. - (C5.4)
*h*=_{X}*O*(*n*^{−κ/5}).

LEMMA **2** *Let* Ƶ (*s, t*) *be a zero-mean Gaussian random field on* {(*s, t*); *s, t* }, *with covariance function* *defined in Assumption (C5.1). Under assumptions (C4.1)–(C4.2) and (C5)*, $\Vert {\tilde{X}}_{i}-{X}_{i}\Vert ={o}_{p}({n}^{-1/2}),\sqrt{n}(\tilde{R}-R)\to \u01b5$, *and*

$$\begin{array}{c}{\widehat{\psi}}_{j}(t)-{\psi}_{j}(t)={n}^{-1/2}{\displaystyle {\sum}_{k\ne j}{({\omega}_{j}-{\omega}_{k})}^{-1}{\psi}_{k}(t)}{\displaystyle \int \u01b5{\psi}_{j}{\psi}_{k}+{o}_{p}({n}^{-1/2}),}\hfill \\ {\widehat{\xi}}_{\mathit{\text{ij}}}-{\xi}_{\mathit{\text{ij}}}={n}^{-1/2}{\displaystyle {\sum}_{k\ne j}{\xi}_{\mathit{\text{ik}}}{({\omega}_{j}-{\omega}_{k})}^{-1}}{\displaystyle \int \u01b5{\psi}_{j}{\psi}_{k}+{o}_{p}({n}^{-1/2}).}\hfill \end{array}$$

The limit distribution of in Lemma 2 was given in Theorem 4 in Zhang and Chen (2007), whereas the last two equations are direct results of (2.8) of Hall and Hosseini-Nasab (2006).

By plugging in the estimated principal component scores, we have the following asymptotic results for the refined estimator defined in Section 3.

THEOREM **4** *Let* *be the refined estimator with plugged in estimated PC scores. Under conditions (C1), (C3), (C4.1)–(C4.2) and (C5)*,

$$\sqrt{n}(\widehat{\mathbf{\Theta}}-{\mathbf{\Theta}}_{0})\to \mathit{\text{Normal}}\{\mathbf{0},{\mathcal{A}}^{-}(\mathcal{B}+{\mathcal{B}}_{1}){A}^{-}\},$$

*where* , and _{1} *are defined in the* Appendix. *In addition*, * Ŝ*(·)

By comparing the asymptotic distributions in Theorem 2 and Theorem 4, we can clearly see the additional variation, _{1}, which is the consequence of the estimation error in the estimated principal components.

**Remarks:**

- Taking the numerical integration method in (9) as an example, there are two sources of error in estimating FPCA scores: the first source is the prediction error for the PC scores given the true eigenfunctions; the second type of error is caused by plugging in the estimated eigenfunctions. Both sources of error tend to be ignored in the functional data literature.
- One important contribution of our work is that we bring attention to the second source of error mentioned above. Even if we have all the information on each curve, the eigenfunction can only be estimated in a root-
*n*rate (Hall and Hosseini-Nasab, 2006). This estimation error may seem to be small, but it affects PC score estimation in all subjects, and as illustrated in our theory this error will also affect the variance of the final estimator in regression. - Plugging in the predicted score is an idea similar to regression calibration (Carroll et al., 2006). With condition (C5.2), we basically assumed that the first type of error (prediction error for the PC scores given the true eigenfunctions) is negligible, which is a typical assumption in the functional data literature, see Müller and Stadtmüller (2005). This allows us to focus all attention on the second source of errors. Our model considers complex functional interactions and by its nature would need moderate to large sample sizes in order to reach sensible outcomes and reliable findings. When
*n*and/or*m*are small, whether one could or should consider such a complex model is an interesting yet hard problem, which calls for future research. - When condition (C5.2) is violated, the first type of error described above will prevail. Since these prediction errors are independent among subjects, they are similar to the classical error in the measurement error literature and tend to cause bias in the final estimator. Therefore, whether
*m*is sufficiently large can be empirically examined by checking the bias of the final estimator in a simulation study.

We performed a simulation study to illustrate the numerical performance of the proposed method. Our simulation setting resembles the empirical data described in Section 5.2. We generated the longitudinal process as a Gaussian process with mean function μ_{X}(*t*) = (*t* − 0.6)^{2} − 0.1 for *t* [0, 1]. The covariance function of the process had 2 principal components, ${\psi}_{1}(t)=1,{\psi}_{2}(t)=\sqrt{2}\phantom{\rule{thinmathspace}{0ex}}\text{sin}(2\pi t)$, and the eigenvalues were ω_{1} = 1 and ω_{2} = 0.6. There were *m* = 30 discrete observations on each curve, with random observation time points being uniformly distributed on the interval [0, 1]. The measurements are contaminated with zero-mean Gaussian error with variance ${\sigma}_{\epsilon}^{2}=0.1$. We generated a binary response *Y* from the model $\text{pr}(Y=1|\mathit{Z},\mathit{\xi})=H\{{\mathit{\alpha}}_{1}^{\mathrm{T}}\mathit{\xi}+\mathit{S}({\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{1}){\mathit{\alpha}}_{2}^{\mathrm{T}}\mathit{\xi}+{\mathit{\beta}}^{\mathrm{T}}\mathit{Z}\}$, where *H*(η) = 1/{1 + exp(η)} is the logistic distribution function, $\mathit{Z}={(1,{\mathit{Z}}_{1}^{\mathrm{T}},{Z}_{2})}^{\mathrm{T}},{\mathit{Z}}_{1}$ is a 2-dimensional random vector with a uniform distribution on [0, 1]^{2}, and *Z*_{2} is a binary variable with pr(*Z*_{2} = 1) = 0.5. We let ${\mathit{\alpha}}_{1}={(0.8,-1)}^{\mathrm{T}},{\mathit{\alpha}}_{2}={(1/\sqrt{2},1/\sqrt{2})}^{\mathrm{T}},\mathit{\theta}={(1/\sqrt{2},1/\sqrt{2})}^{\mathrm{T}}$, and **β** = (−1, 2, −2, 2). We let * S*(·) be a sine bump function similar to that used in Carroll et al. (1997),

We set *n* = 500 and repeated the simulation 200 times. For each simulated data set we performed FPCA and fit the model using the algorithm described in Section 3, where ξ_{ik} was predicted by the direct numerical integral predictor. We obtained almost identical outcomes using the PACE predictor (not reported). For the iteration procedure, we declared convergence if max_{k}|**Θ̃**_{k,curr} − **Θ̃**_{k,prev}| is less than a pre-determined value, 10^{−6}, where **Θ̃**_{k,prev} and **Θ̃**_{k,curr} denote the previous and current value of the *k*-th component of **Θ̃**. The number of principle component, *p*, was determined using AIC with the total number of parameters in the penalty equal to the number of PC scores.

The Monte Carlo (MC) means, standard deviations (SD) and biases of the refined estimator are presented in Table 1. All estimates behaved reasonably well. The comparisons of the MC biases and the MC SD’s indicate that the proposed estimator has a relatively smaller bias than the standard error. This seems to confirm the theoretical results that they are asymptotically unbiased. For all simulated data sets, the AIC criterion selected *p* = 2 which is the correct number of components. To confirm the efficiency gain of our refined estimator over the initial estimator **Θ̃**, we made a comparison of the mean squared error (MSE) of the two estimators. We found that MSE(**Θ̃**)/MSE() = 1.34. In other words, we have 34% efficiency gain by switching to refined kernel weights.

Estimation results for the simulation study. The reduced model is misspecified because it does not include the interaction, and the biases observed in the reduced model are thus to be expected.

We compare the proposed approach to two other alternatives. First, we fit a misspecified reduced model without interactions, where $P(Y=1|\mathit{Z},\mathit{\xi})=H({\mathit{\alpha}}_{1}^{\mathrm{T}}\mathit{\xi}+{\mathit{\beta}}^{\mathrm{T}}\mathit{Z})$. The results are reported in Table 1. An immediate observation is elevated biases in the estimates of **α**_{1} and **β** due to the misspecification of the model by not considering the interaction.

To assess the variation-inflation effect due to the PC scores being predicted, we also fit the model using the true principal component scores and compare the variances of the final estimators with those in Table 1. In Table 2, we present the relative variance inflation factor for using the true PC scores, which is defined as the ratio of the variance for each parameter when estimating the FPCA scores to the variance when the true FPCA scores are known. One clearly sees that the relative variance inflation factors are uniformly greater than one, with the highest inflation occurring for the coefficients of the second principal component, α_{12} and α_{22}. This phenomena can be explained by the fact that the second PC is much harder to estimate than the first.

The colon carcinogenesis data that we are using are similar to those analyzed in Morris et al. (2001, 2003), Li et al. (2007) and Baladandayuthapani et al. (2008). The biomarker of interest in this experiment is p27, which is a protein that inhibits the cell cycle. We used 12 rats randomly assigned to 2 diet groups (corn oil diet or fish oil diet) and 2 treatment groups (with/without butyrate supplement). Each rat was injected with a carcinogen, and then sacrificed 24 hours after the injection.

Beneath the colon tissue, there are pore structures called ‘colonic crypts’. A crypt typically contains 25 to 30 cells, lined up from the bottom to the top. The stem cells are at the bottom of the crypt, where daughter cells are generated. These daughter cells move towards the top as they mature. We sampled about 20 crypts from each of the 12 rats. The p27 expression level was measured for each cell within the sampled crypts. As previously noted in the literature (Morris et al. 2001, 2003), the p27 level measurements in the logarithm scale are natural functional data, since they are indexed by the relative cell location within the crypt. We have *m* = 25–30 observations (cells) on each function. In the literature, it has been noted that there is spatial correlation among the crypts within the same rat (Li et al., 2007, Baladandayuthapani et al., 2008). In this experiment, we sampled crypts sufficiently far apart so that the spatial correlations are negligible, and thus we can assume that the crypts are independent.

In this example, the response *Y* is the fraction of cells undergoing apoptosis (programmed cell death) in each crypt, to which we applied the usual arcsine-square root transformation. The functional data *X*(*t*) is the p27 level of a cell at relative location *t*. Write *s*_{0} = 1, *s*_{1} = indicator of a fish oil diet, *s*_{2} = indicator of butyrate supplementation, and *s*_{3} = mean proliferation level. Then ${\mathit{Z}}_{1}={({s}_{1},{s}_{2},{s}_{3})}^{\mathrm{T}},\mathit{Z}={({s}_{0},{\mathit{Z}}_{1}^{\mathrm{T}})}^{\mathrm{T}}$ and *Z*^{T}**β** = β_{0} + *s*_{1}β_{fish} + *s*_{2}β_{buty} + *s*_{3}β_{prolif}. Similarly, * S*(

We first performed FPCA on the p27 data. In Figure 1, we show _{X}(*t*), (*s, t*) and the first 3 estimated PCs. The first 3 eigenvalues are 0.871, 0.019 and 0.005 respectively. The estimated variance for the measurement error is ${\widehat{\sigma}}_{\epsilon}^{2}=0.103$. We again implemented the AIC criterion to choose the number of significant principal components in the p27 data, and *p* = 2 was picked. We thus only include the first 2 principal components of the p27 data in the regression model, so that **α**_{1} = (α_{11}, α_{12})^{T} and ${\mathit{\alpha}}_{2}=({\alpha}_{21},{\alpha}_{22}^{\mathrm{T}})$. The estimated parameters are presented in Table 3, and the estimated * S*(·) is shown in Figure 2.

Functional principal component analysis for the colon carcinogenesis p27 data. Top panel: local linear estimator of the mean curve μ_{X}(*t*); middle panel: local linear estimator of the covariance function; lower panel: the first 3 principal components. **...**

Estimated parameters in the colon carcinogenesis data. The standard errors were estimated by bootstrap, and the p-values are for the null hypothesis that the parameter = 0.0. Here the two principal component values have **α**_{1} = (α_{11}, α **...**

We estimated the standard error of our parameter estimates by a bootstrap procedure. Of course, after resampling the individual curves (crypts) we reran the functional principal component analysis in the bootstrap sample, so that the bootstrap standard error estimator automatically takes into account the estimation error in FPCA. We then refit the model to the new sample. The procedure was repeated 1, 000 times, and the standard deviations of the estimators in the bootstrap samples were used as our standard error estimators. We then tested each parameter in turn to see whether the true parameter is equal to 0.00. The estimated standard errors and the p-values are reported in Table 3.

To summarize, the parameters, β_{fish}, α_{22}, θ_{buty} and θ_{prolif}, have corresponding p-values less than 0.01. For the parametric main effect, the positive significant β_{fish} implies that the fish oil diet enhances the apoptotic index. The positive _{buty} and the negative _{prolif} suggest that the butyrate supplement and the increasing proliferation level are positively and negatively associated with the apoptotic index, respectively, even though neither is significant at the main effect level. The statistically significant results for α_{22}, θ_{buty} and θ_{prolif} imply that the second principal component of the p27 process is interacting with butyrate supplementation and proliferation.

In addition to testing if each individual parameter is equal to 0.00, it is important to have an overall test of whether interaction really exists in this data set. In our model, this is equivalent to testing the null hypothesis *H*_{0} : * S*(·) 0. We test this null hypothesis using the following parametric bootstrap procedure:

**Step 1.** Fit the model under the null hypothesis given as $\mathrm{E}({Y}_{i}|{\mathit{\xi}}_{i},{\mathit{Z}}_{i})={\mathit{\xi}}_{i}^{\mathrm{T}}{\mathit{\alpha}}_{1}+{\mathit{Z}}_{i}^{\mathrm{T}}\mathit{\beta}$. Denote the estimators under the reduced model as _{1,red} and _{red}, and call the estimated conditional variance of *Y* under the reduced model ${\widehat{\sigma}}_{Y,\text{red}}^{2}$. We assume the data are Gaussian, and calculate the log likelihood ratio of the full model with interaction to the reduced model above.

**Step 2.** Use the estimated reduced model as the true one to generate bootstrap samples. We resample with replacement *n* crypts from the original collection. For the *i ^{th}* resampled crypt, we retain the original covariate vector

To apply the estimation procedure to the bootstrap sample, we performed the principal component analysis to the regenerated observations $\{({t}_{\mathit{\text{ij}}},{W}_{\mathit{\text{ij}}}^{*});j=1,\cdots ,{m}_{i},i=1,\cdots ,n\}$, so that our test procedure automatically takes into account the estimation error in FPCA. After FPCA, we apply both the full model and the reduced model to the bootstrap sample and calculate the log likelihood ratio. The bootstrap procedure was repeated 1, 000 times.

**Step 3.** Compare the log likelihood ratio calculated from the data to the bootstrap quantities, and let the p-value be the upper percentile of the true log likelihood ratio among the bootstrap versions.

We applied this procedure to the colon carcinogenesis data, and the test yields a p-value of 0.044, which is indication that interaction does exist in these data. To see the effect of this interaction, we also present (*t*, **θ**^{T}**Z**_{1}) in Figure 3. We can see, as the single-index value **θ**^{T}**Z**_{1} changes, the coefficient function for the functional predictor changes quite dramatically. Other inference procedures have been proposed in different semiparametric models. Among them, Liang et al. (2009) proposed an empirical likelihood based inference for generalized partially linear models and Li and Liang (2008) proposed to use a semiparametric generalized likelihood ratio test to select significant variables in the nonparametric component in the generalized varying-coefficient partially linear models. How to adopt their testing ideas in our models warrants future research.

One interesting finding of our study is that when we apply the reduced model without interaction to the data, the only significant parameter is the main effect of diet. In other words, none of butyrate supplementation, proliferation and p27 show strong main effects. These factors only show their effects in the full model that includes an interaction.

To further confirm this conclusion, we also did another simulation. We bootstrapped crypts from the original data, but regenerated *Y* using our estimated full model with interactions. We then applied the reduced model to the simulated data sets. Again, none of the factors besides the fish oil diet showed significant main effects.

We now further look into the interactions between p27 and other covariates through the PC scores of p27 and * S*(·) values of the single index

Proportions of crypts receiving butyrate supplementation (left panel) and fish oil (right panel), depending on the levels of **S**(·). The left bar in each panel refers all crypts, the middle bar to the case that **S**(·) is high and the right **...**

Boxplots of proliferation depending on values of **S**(·). Here ”All”, ”High” and ”Low” refer to all observations, those with high values of **S**(·) and those with low values of **S**(·).

To understand the interacting effects of p27 and **θ**^{T}**Z**_{1} on the apoptotic levels through the p27 PC scores, we produce Figure 6. We first removed the estimated main effects by subtracting them from the apoptotic levels and refer to the adjusted values as apoptotic indices. We then dichotomized each of the two p27 PC scores according whether they belong to the top or bottom 50% of the scores; this practice allows us to produce 4 groups for the 2 PC-scores: PC1-Low, PC1-High, PC2-Low and PC2-High. An entry that belongs to PC1-Low tends to have a lower than average p27 and that in PC2-Low tends to have its p27 values lower in the top of the crypts. The equivalent implications apply to the PC1-High and PC2-High groups. We calculated the average apoptotic indices of these four PC groups against *S*_{low}, *S*_{mid} and *S*_{high} and plot them in Figure 6. The solid circle and square indicate the membership of PC1 and PC2, while the solid, dashed, dash-dotted and dotted lines indicate the membership of PC1-High, PC1-Low, PC2-High and PC2-Low, respectively. Figure 6 clearly shows that there are interactions between PC scores of p27 and * S*(

We have introduced a new class of functional generalized linear models that includes multivariate covariates and their interaction with the functional predictor. To achieve parsimony, the interaction is modeled semiparametrically with a single-index structure.

When the functional features related to the response variable are known basis functions, we proposed a two step procedure based on local estimation equations. We showed theoretically that the first step of the procedure produces a consistent estimator, while the second step estimator is more efficient and is asymptotically normal with a root *n* rate.

We further investigated the situation that the functional features are data-driven, i.e. the principal components of the process. We applied the functional principal component analysis of Yao et al. (2005) and Hall et al. (2006) to the discrete observations on the functional data, and plugged the estimated principal component scores into the two step procedure. We showed that when we have a large number of subjects, and enough number of observations on each curve, the initial estimator in our first step procedure is still consistent, while the second step estimator is still root-*n* consistent and asymptotically normal.

One important finding is that asymptotic variance of our second step estimator is larger than in the case that the functional features **ξ** are known. This extra variation in the estimator is due to the estimation error of the FPCA.

Our theoretical investigation is innovative since almost all the literature on functional data analysis that we are aware of tends to ignore the estimation error in FPCA. The only exception is perhaps the recent paper by Li, Li, Wang and Wang (2008), where estimation error in FPCA was taken into account, but their model did not include interactions between the functional predictor and other covariates. Our results verify that one should be aware of the extra variability when estimating FPCA scores.

Our methods were illustrated by a simulation study and applied to a colon carcinogenesis data set. In studying this FGLM model, we detected a previously unknown interaction of the response, p27, with the environmental covariates.

Throughout, η{**α, β, θ**, * S*(·)} = η{

Define ${\epsilon}_{i}={q}_{1}\left\{{\eta}_{i},{Y}_{i}\right\}$,

$$\begin{array}{c}{\mathcal{V}}_{1}(\mathbf{z};\mathit{\theta})=\mathrm{E}\{{\rho}_{2}({\eta}_{i}){({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})}^{2}|{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}={\mathit{\theta}}^{\mathrm{T}}\mathbf{z}\},\hfill \\ {\mathit{\pi}}_{i1}={\rho}_{2}({\eta}_{i})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\left\{\begin{array}{c}{\mathit{\xi}}_{i}\\ \mathit{S}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}){\mathit{\xi}}_{i}\\ {\mathit{Z}}_{i}\end{array}\right\},\text{}{\mathit{\pi}}_{i2}={\rho}_{2}({\eta}_{i}){({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})}^{2}{\mathit{Z}}_{i1},\hfill \\ {\widehat{\mathit{\pi}}}_{1}(\mathbf{z};\mathit{\theta})=\mathrm{E}({\mathit{\pi}}_{i1}|{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}={\mathit{\theta}}^{\mathrm{T}}\mathbf{z}),\text{}{\widehat{\mathit{\pi}}}_{2}(\mathbf{z};\mathit{\theta})=\mathrm{E}({\mathit{\pi}}_{i2}|{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}={\mathit{\theta}}^{\mathrm{T}}\mathbf{z}),\hfill \end{array}$$

and _{3}(**z**; **θ**) is a symmetric matrix whose (_{1}, _{2})^{th} block is denoted by ${\mathcal{V}}_{3}^{[{\ell}_{1},{\ell}_{2}]}\phantom{\rule{thinmathspace}{0ex}}\text{for}\phantom{\rule{thinmathspace}{0ex}}{\ell}_{1},{\ell}_{2}=1,2$, and

$$\begin{array}{c}{\mathcal{V}}_{3}^{[1,1]}(\mathbf{z};\mathit{\theta})=\mathrm{E}\text{}[{\rho}_{2}({\eta}_{i})\phantom{\rule{thinmathspace}{0ex}}{\left\{\mathit{S}\begin{array}{c}{\mathit{\xi}}_{i}\\ ({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1})\\ {\mathit{Z}}_{i}\end{array}{\mathit{\xi}}_{i}\right\}}^{\otimes 2}|{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}={\mathit{\theta}}^{\mathrm{T}}\mathbf{z}]\phantom{\rule{thinmathspace}{0ex}},\hfill \\ {\mathcal{V}}_{3}^{[1,2]}(\mathbf{z};\mathit{\theta})={\mathit{S}}^{(1)}({\mathit{\theta}}_{0}^{\mathrm{T}}\mathbf{z})\{\mathrm{E}({\mathit{\pi}}_{i1}{\mathit{Z}}_{i1}^{\mathrm{T}}|{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}={\mathit{\theta}}^{\mathrm{T}}\mathbf{z})-{\widehat{\mathit{\pi}}}_{1}(\mathbf{z};\mathit{\theta}){\mathit{z}}^{\mathrm{T}}\},\hfill \\ {\mathcal{V}}_{3}^{[2,1]}(\mathbf{z};\mathit{\theta})={\{{\mathcal{V}}_{3}^{[1,2]}(\mathbf{z};\mathit{\theta})\}}^{\mathrm{T}},\hfill \\ {\mathcal{V}}_{3}^{[2,2]}(\mathbf{z};\mathit{\theta})={\{{\mathit{S}}^{(1)}({\mathit{\theta}}_{0}^{\mathrm{T}}\mathbf{z})\}}^{2}[\mathrm{E}\{{\rho}_{2}({\eta}_{i}){({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})}^{2}{\mathit{Z}}_{i1}{\mathit{Z}}_{i1}^{\mathrm{T}}|{\mathit{\theta}}^{\mathrm{T}}{\mathit{Z}}_{i1}={\mathit{\theta}}^{\mathrm{T}}\mathbf{z}\}-{\widehat{\mathit{\pi}}}_{2}(\mathbf{z};\mathit{\theta}){\mathit{z}}^{\mathrm{T}}-\mathbf{z}{\widehat{\mathit{\pi}}}_{2}^{\mathrm{T}}(\mathbf{z};\mathit{\theta})+{\mathcal{V}}_{1}(\mathbf{z};\mathit{\theta}){\mathbf{\text{zz}}}^{\mathrm{T}}].\hfill \end{array}$$

Let

$${\mathcal{N}}_{n}={n}^{-1}{\displaystyle {\sum}_{j=1}^{n}{\epsilon}_{i}\left[\begin{array}{c}{\mathit{\xi}}_{i}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1};{\mathit{\theta}}_{0})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\mathrm{E}\{{\rho}_{2}(\eta )({\mathit{\alpha}}_{2,0}^{\mathrm{T}}\mathit{\xi})\mathit{\xi}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}\}\\ \mathit{S}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}){\mathit{\xi}}_{i}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1};{\mathit{\theta}}_{0})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\mathrm{E}\{{\rho}_{2}(\eta )({\mathit{\alpha}}_{2,0}^{\mathrm{T}}\mathit{\xi})\mathit{S}|({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1})\xi |{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}\}\\ {\mathit{Z}}_{i}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1};{\mathit{\theta}}_{0})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\mathrm{E}\{{\rho}_{2}(\eta )({\mathit{\alpha}}_{2,0}^{\mathrm{T}}\mathit{\xi})\mathit{Z}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}\}\\ ({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i}){\mathit{S}}^{(1)}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1})\{{\mathit{Z}}_{i1}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1;{\mathit{\theta}}_{0}}){\widehat{\mathit{\pi}}}_{2}({\mathit{Z}}_{i1};{\mathit{\theta}}_{0})\}\end{array}\right]\text{},}$$

and ${\mathcal{N}}_{1,n}={n}^{-1}{\displaystyle {\sum}_{i=1}^{n}{\epsilon}_{i}({\widehat{\mathit{\xi}}}_{i}-{\mathit{\xi}}_{i})}$ where

$$\begin{array}{cc}& {\epsilon}_{i}=\left[\begin{array}{c}{\mathit{\xi}}_{i}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1};{\mathit{\theta}}_{0})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\mathrm{E}\{{\rho}_{2}(\eta )({\mathit{\alpha}}_{2,0}^{\mathrm{T}}\mathit{\xi})\mathit{\xi}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}\}\\ \mathit{S}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}){\mathit{\xi}}_{i}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1};{\mathit{\theta}}_{0})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\mathrm{E}\{{\rho}_{2}(\eta )({\mathit{\alpha}}_{2,0}^{\mathrm{T}}\mathit{\xi})\mathit{S}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1})\mathit{\xi}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}\}\\ {\mathit{Z}}_{i}-{\mathcal{V}}_{1}^{-1}({Z}_{i1};{\mathit{\theta}}_{0})({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i})\mathrm{E}\{{\rho}_{2}(\eta )({\mathit{\alpha}}_{2,0}^{\mathrm{T}}\mathit{\xi})\mathit{Z}|{\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{1}={\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}\}\\ ({\mathit{\alpha}}_{2,0}^{\mathrm{T}}{\mathit{\xi}}_{i}){\mathit{S}}^{(1)}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1})\{{\mathit{Z}}_{i1}-{\mathcal{V}}_{1}^{-1}({\mathit{Z}}_{i1;{\mathit{\theta}}_{0}}){\widehat{\mathbf{\pi}}}_{2}({\mathit{Z}}_{i1},{\mathit{\theta}}_{0})\}\end{array}\right]\\ & \text{\hspace{1em}}\times {q}_{2}({\eta}_{i},{Y}_{i}){\{{\mathit{\alpha}}_{1,0}+{\mathit{S}}_{0}({\mathit{\theta}}_{0}^{\mathrm{T}}{\mathit{Z}}_{i1}){\mathit{\alpha}}_{2,0}\}}^{\mathrm{T}}.\hfill \end{array}$$

Finally, define , and *B*_{1} in Theorems 2 and 4 as:

$$\mathcal{A}=\mathrm{E}\{{\mathcal{V}}_{3}({\mathit{Z}}_{1};{\mathit{\theta}}_{0})\},\text{}\mathcal{B}=n\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}\text{cov}({\mathcal{N}}_{n}),\text{}{\mathcal{B}}_{1}=\underset{n\to \infty}{\text{lim}}n\phantom{\rule{thinmathspace}{0ex}}\times \phantom{\rule{thinmathspace}{0ex}}\text{cov}({\mathcal{N}}_{1,n}).$$

^{1}Li’s research was supported by the National Science Foundation (DMS-0806131). Wang’s research was supported by a grant from the National Cancer Institute (CA74552). Carroll’s research was supported by a grant from the National Cancer Institute (CA57030) and by Award Number KUS-CI-016-04, made by King Abdullah University of Science and Technology (KAUST).

**Supplemental Materials**

**Technical Details:** The web appendix describes the detailed estimation algorithm, and provides technical proofs for all of our theoretical results (pdf file).

Yehua Li, Department of Statistics, University of Georgia, Athens, GA 30602, Email: ude.agu@ilauhey..

Naisyin Wang, Department of Statistics, University of Michigan, Ann Arbor, MI 48109-1107, Email: ude.hcimu@aagnawn..

Raymond J. Carroll, Department of Statistics, Texas A&M University, TAMU 3143, College Station, TX 77843-3143, Email: ude.umat.tats@llorrac..

- Baladandayuthapani V, Mallick B, Hong M, Lupton J, Turner N, Carroll RJ. Bayesian hierarchical spatially correlated functional data analysis with application to colon carcinogenesis. Biometrics. 2008;64:64–73. [PMC free article] [PubMed]
- Cai T, Hall P. Prediction in functional linear regression. Annals of Statistics. 2006;34:2159–2179.
- Cardot H, Sarda P. Estimation in generalized linear models for functional data via penalized likelihood. Journal of Multivariate Analysis. 2005;92:24–41.
- Carroll RJ, Fan J, Gijbels I, Wand MP. Generalized partially linear single-index models. Journal of the American Statistical Association. 1997;92(438):477–489.
- Carroll RJ, Ruppert D, Stefanski LA, Crainiceanu CM. Measurement Error in Nonlinear Models: A Model Perspective. 2nd Edition. New York: Chapman & Hall; 2006.
- Capra WB, Müller HG. An accelerated-time model for response curves. Journal of the American Statistical Association. 1997;92:72–83.
- Crainiceanu CM, Staicu A, Di C. Generalized Multilevel Functional Regression. Johns Hopkins University, Dept. of Biostatistics Working Papers; 2008. Working Paper 173. http://www.bepress.com/jhubiostat/paper173.
- Di C, Crainiceanu CM, Caffo BS, Punjabi NM. Multilevel functional principal component analysis. Annals of Applied Statistics. 2009 to appear. [PMC free article] [PubMed]
- Ferraty F, Vieu P. Nonparametric Functional Data Analysis: Theory and Practice. New York: Springer-Verlag; 2006.
- Hall P, Hosseini-Nasab M. On properties of functional principal components analysis. Journal of the Royal Statistical Society, Series B. 2006;68:109–126.
- Hall P, Müuller H-G, Wang J-L. Properties of principal component methods for functional and longitudinal data analysis. Annals of Statistics. 2006;34:1493–1517.
- Li E, Li Y, Wang N-Y, Wang N. Functional latent feature models for data with longitudinal covariate processes. 2008 Preprint.
- Li KC. Sliced inverse regression for dimension reduction (with Discussion) Journal of the American Statistical Association. 1991;86 316342.
- Li R, Liang H. Variable selection in semiparametric regression modeling. Annals of Statistics. 2008;36:261–286. [PMC free article] [PubMed]
- Li Y, Wang N, Hong M, Turner N, Lupton J, Carroll RJ. Nonparametric estimation of correlation functions in spatial and longitudinal data, with application to colon carcinogenesis experiments. Annals of Statistics. 2007;35:1608–1643.
- Liang H, Qin YS, Zhang XY, Ruppert D. Empirical-likelihood-based inferences for generalized partially linear models. Scandinavian Journal of Statistics. 2009;36:433–443. [PMC free article] [PubMed]
- Mack YP, Silverman BW. Weak and Strong Uniform consistency of kernel regression estimates,
*Z. Wahrscheinlichkeitstheorie verw*. Gebiete. 1982;61:405–415. - McCullagh P, Nelder JA. Generalized Linear Models. 2nd Edition. New York: Chapman & Hall; 1989.
- Morris JS, Wang N, Lupton JR, Chapkin RS, Turner ND, Hong MY, Carroll RJ. Parametric and nonparametric methods for understanding the relationship between carcinogen-induced DNA adduct levels in distal and proximal regions of the colon. Journal of the American Statistical Association. 2001;96(455):816–826.
- Morris JS, Vannucci M, Brown PJ, Carroll RJ. Wavelet-based non-parametric modeling of hierarchical functions in colon carcinogenesis. Journal of the American Statistical Association. 2003;98(463):573–583.
- Müller H-G, Stadtmüller U. Generalized functional linear models. Annals of Statistics. 2005;33:774–805.
- Ramsay JO, Silverman BW. Functional Data Analysis. 2nd Edition. New York: Springer-Verlag; 2005.
- Rice J, Silverman B. Estimating the mean and covariance structure nonparametrically when the data are curves. Journal of the Royal Statistical Society, Series B. 1991;53:233–243.
- Xia Y, Tong H, Li WK, Zhu L. An adaptive estimation of dimension reduction space (with discussion) Journal of the Royal Statistical Society, Series B. 2002;64:363–410.
- Xia Y, Härdle W. Semi-parametric estimation of partially linear single-index models. Journal of Multivariate Analysis. 2006;97:1162–1184.
- Yao F, Müller HG, Wang JL. Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association. 2005;100:577–590.
- Zhang J-T, Chen J. Statistical inferences for functional data. Annals of Statistics. 2007;35:1052–1079.
- Zhou L, Huang JZ, Carroll RJ. Joint modeling of paired sparse functional data using principal components. Biometrika. 2008;95(3):601–619. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |