A partial least squares approach for function-on-function interaction regression

Abstract

A partial least squares regression is proposed for estimating the function-on-function regression model where a functional response and multiple functional predictors consist of random curves with quadratic and interaction effects. The direct estimation of a function-on-function regression model is usually an ill-posed problem. To overcome this difficulty, in practice, the functional data that belong to the infinite-dimensional space are generally projected into a finite-dimensional space of basis functions. The function-on-function regression model is converted to a multivariate regression model of the basis expansion coefficients. In the estimation phase of the proposed method, the functional variables are approximated by a finite-dimensional basis function expansion method. We show that the partial least squares regression constructed via a functional response, multiple functional predictors, and quadratic/interaction terms of the functional predictors is equivalent to the partial least squares regression constructed using basis expansions of functional variables. From the partial least squares regression of the basis expansions of functional variables, we provide an explicit formula for the partial least squares estimate of the coefficient function of the function-on-function regression model. Because the true forms of the models are generally unspecified, we propose a forward procedure for model selection. The finite sample performance of the proposed method is examined using several Monte Carlo experiments and two empirical data analyses, and the results were found to compare favorably with an existing method.

Appendix

1.1 Identifiability of the proposed model

Model identifiability, which is an important theme in the estimation phase, for the function-on-function regression model with only main effects given in (1.1) was theoretically discussed by He et al. (2000) and Chiou et al. (2004). Later, Scheipl and Greven (2016) discussed the identifiability of such models with realistic applications. As is stated in the Proposition 3.1 of Scheipl and Greven (2016), the coefficient function $\beta _m(s,t)$ in (1.1) is identifiable if and only if the kernel space of the covariance operator of the functional predictors ($K^{{\mathcal {X}}}$) is empty except the zero function, i.e., $ke(K^{{\mathcal {X}}}) = \lbrace 0 \rbrace$. Luo and Qi (2019) extended the theoretical condition of model identifiability for function-on-function regression with only main effects to the model with interaction and quadratic effects. They defined the following $\lbrace p (p + 3)/2 \rbrace \times \lbrace p (p+3)/2 \rbrace$ dimensional covariance matrix of the predictor functions:

$$\begin{aligned} \pmb {\Sigma } = \begin{bmatrix} \pmb {A}(s,r) &{}\quad \pmb {B}(s,r,u) \\ \pmb {B}^\top (s,r,u) &{}\quad \pmb {C}(s, r, s^\prime , r^\prime ) \end{bmatrix}, \end{aligned}$$

(5.9)

where the $p \times p$ and $p \times \lbrace p (p+1)/2 \rbrace$ dimensional submatrices $\pmb {A}(s,r)$ and $\pmb {B}(s,r,u)$ equal to $\text {Cov} \left( {\mathcal {X}}(s) {\mathcal {X}}(r) \right)$ and $\text {Cov} \left( {\mathcal {X}}(s), {\mathcal {X}}(r) {\mathcal {X}}(u) \right)$, respectively. The $\lbrace p (p+1)/2 \rbrace \times \lbrace p (p+1)/2 \rbrace$ dimensional matrix $\pmb {C}(s, r, s^\prime , r^\prime )$ equals to $\text {Cov} \left( {\mathcal {X}}(s) {\mathcal {X}}(r), {\mathcal {X}}(s^\prime ) {\mathcal {X}}(r^\prime ) \right)$. Then, Luo and Qi (2019) proved that the coefficients in model (1.3) are identifiable if the kernel space of $\pmb {\Sigma }$ in (5.9) only has the zero function. Please see Proposition S.1. of Luo and Qi (2019) for the proof. Based on the above, our proposed method given in (2.4) is also identifiable since the proposed method is a reformulated version of the model (1.3).

1.2 Proof of Proposition (2.2)

Proof

The proof of Proposition (2.2) is an extended version of the proof of Proposition 2 of Aguilera et al. (2010). Denote by $\Lambda = \pmb {\Phi }^{1/2} \pmb {c}$ and $\Pi = \pmb {\Psi }^{1/2} \pmb {d}$ the column random vectors.

For any random variable $Z \in \mathcal {L}_2[0,1]$, the followings hold when $h = 1$:

$$\begin{aligned} W^{{\mathcal {Y}}} =\,&\int _0^1 \sum _{k=1}^{K_{{\mathcal {Y}}}} c_k \phi _k(t) \mathbb {E} \left[ \sum _{i=1}^{K_{{\mathcal {Y}}}} Z c_i \phi _i(t) \right] dt = \Lambda ^\top \mathbb {E} \left[ \Lambda Z \right] = W^{\Lambda } Z, \\ W^{\pmb {{\mathcal {X}}}} =\,&\int _0^1 \int _0^1 \sum _{j=1}^{K_{\pmb {{\mathcal {X}}}}} \sum _{l=1}^{K_{\pmb {{\mathcal {X}}}}} d_{jl} \psi _{jl}(s,r) \mathbb {E} \left[ \sum _{\iota = 1}^{K_{\pmb {{\mathcal {X}}}}} \sum _{\tau = 1}^{K_{\pmb {{\mathcal {X}}}}} Zd_{\iota \tau } \psi _{\iota \tau }(s,r) \right] \\&\quad ds dr = \Pi ^\top \mathbb {E} \left[ \Pi Z \right] = W^{\Pi } Z. \end{aligned}$$

The residuals obtained from the first iteration are given by:

$$\begin{aligned}&{\left\{ \begin{array}{ll} {\mathcal {Y}}_1(t) = {\mathcal {Y}}(t) - \zeta _1(t) \eta _1, &{} \zeta _1(t) \in \mathcal {L}_2[0,1] \\ \Lambda _1 = \Lambda - \widetilde{\zeta }_1 \eta _1, &{} \widetilde{\zeta _1} \in \mathbb {R}^{K_{{\mathcal {Y}}}} \end{array}\right. } \\&{\left\{ \begin{array}{ll} \pmb {{\mathcal {X}}}_1(s,r) = \pmb {{\mathcal {X}}}(s,r) - p_1(s,r) \eta _1, &{} p_1(s,r) \in \mathcal {L}_2[0,1] \\ \Pi _1 = \Pi - \widetilde{p}_1 \eta _1, &{} \widetilde{p}_1 \in \mathbb {R}^{K_{\pmb {{\mathcal {X}}}}^2} \end{array}\right. } \end{aligned}$$

Then, we have $\Lambda _1 = \pmb {\Phi }^{1/2} \pmb {c}_1$ and $\Pi _1 = \pmb {\Psi }^{1/2} \pmb {d}_1$, where $\pmb {c}_1$ and $\pmb {d}_1$ are the random vectors of the basis coefficients of ${\mathcal {Y}}_1(t)$ and $\pmb {{\mathcal {X}}}_1(s,r)$, respectively. For functional PLS, the residuals calculated from the first iteration are as follows:

$$\begin{aligned} {\mathcal {Y}}_1(t) =\,&\sum _{k=1}^{K_{{\mathcal {Y}}}} c_{1,k} \phi _k(t) = \pmb {c}_1^\top \pmb {\Phi }(t), \qquad \forall t \in \mathcal {L}_2[0,1],\\ \pmb {{\mathcal {X}}}_1(s,r) =\,&\sum _{j=1}^{K_{\pmb {{\mathcal {X}}}}} \sum _{l=1}^{K_{\pmb {{\mathcal {X}}}}} d_{1,jl} \psi _{jl}(s,r) = \pmb {d}_1^\top \pmb {\Psi }(s,r), \qquad \forall s,r \in \mathcal {L}_2[0,1]. \end{aligned}$$

On the other hand, the residuals from the PLS regression are obtained as follows:

$$\begin{aligned} {\mathcal {Y}}_1(t) =\,&\pmb {\Phi }^\top (t) \pmb {c} - \eta _1 \frac{\mathbb {E} \left[ \pmb {\Phi }^\top (t) \pmb {c} \eta _1 \right] }{\mathbb {E} \left[ \eta _1^2\right] },\\ \pmb {{\mathcal {X}}}_1(s,r) =\,&\pmb {\Psi }^\top (s,r) \pmb {d} - \eta _1 \frac{\mathbb {E} \left[ \pmb {\Psi }^\top (s,r) \pmb {d} \eta _1 \right] }{\mathbb {E} \left[ \eta _1^2 \right] }, \\ \Lambda _1 =\,&\pmb {\Phi }^{1/2} \pmb {c} - \eta _1 \frac{\mathbb {E} \left[ \pmb {\Phi }^{1/2} \pmb {c} \eta _1 \right] }{\mathbb {E} \left[ \eta _1^2 \right] },\\ \Pi _1 =\,&\pmb {\Psi }^{1/2} \pmb {d} - \eta _1 \frac{\mathbb {E} \left[ \pmb {\Psi }^{1/2} \pmb {d} \eta _1 \right] }{\mathbb {E} \left[ \eta _1^2 \right] } \end{aligned}$$

and thus we have $\pmb {c}_1 = \pmb {c} - \eta _1 \frac{\mathbb {E} \left( \pmb {c} \eta _1 \right) }{\mathbb {E} \left[ \eta _1^2 \right] }$ and $\pmb {d}_1 = \pmb {d} - \eta _1 \frac{\mathbb {E} \left[ \pmb {d} \eta _1 \right] }{\mathbb {E} \left[ \eta _1^2 \right] }$, that is, $W^{{\mathcal {Y}}_1} = W^{\Lambda _1}$ and $W^{\pmb {{\mathcal {X}}}_1} = W^{\Pi _1}$. Now, assuming $W^{{\mathcal {Y}}_{\ell }} = W^{\Lambda _{\ell }}$ and $W^{\pmb {{\mathcal {X}}}_{\ell }} = W^{\Pi _{\ell }}$ for each $\ell \le h$. Also, let us assume that $W^{{\mathcal {Y}}_{h+1}} = W^{\Lambda _{h+1}}$ and $W^{\pmb {{\mathcal {X}}}_{h+1}} = W^{\Pi _{h+1}}$. Then, at step $h+1$ the followings hold:

$$\begin{aligned}&{\left\{ \begin{array}{ll} {\mathcal {Y}}_{h+1}(t) = {\mathcal {Y}}(t) - \zeta _1(t) \eta _1 - \cdots - \zeta _h(t) \eta _h, &{} \zeta _i(t) \in \mathcal {L}_2[0,1] \\ \Lambda _{h+1} = \Lambda - \widetilde{\zeta }_1 \eta _1 - \cdots - \widetilde{\zeta }_h \eta _h, &{} \widetilde{\zeta _i} \in \mathbb {R}^{K_{{\mathcal {Y}}}} \end{array}\right. } \\&{\left\{ \begin{array}{ll} \pmb {{\mathcal {X}}}_{h+1}(s,r) = \pmb {{\mathcal {X}}}(s,r) - p_1(s,r) \eta _1 - \cdots - p_h(s,r) \eta _h, &{} p_i(s,r) \in \mathcal {L}_2[0,1] \\ \Pi _{h+1} = \Pi - \widetilde{p}_1 \eta _1 - \cdots - \widetilde{p}_h \eta _h, &{} \widetilde{p}_i \in \mathbb {R}^{K_{\pmb {{\mathcal {X}}}}^2} \end{array}\right. } \end{aligned}$$

As for $h=1$ and using the orthogonality of $\eta _i$ for $i = 1, \ldots , h$, it is shown that $\Lambda _{h+1} = \pmb {\Phi }^{1/2} \pmb {c}_{h+1}$ and $\Pi _{h+1} = \pmb {\Psi }^{1/2} \pmb {d}_{h+1}$, which concludes the proof. $\square$

1.3 Station names for the North Dakota weather data

See Table 4.

Table 4 Station names for the North Dakota weather data