Robust function-on-function regression model with nonparametric random effects

1.   Introduction

As the development of science and technology, some data sets are recorded frequently with curves, surfaces and other types, which are usually called functional data that plays an important role in wide fields such as atmospheric science, engineering, medical research, see more details in Ramsay and Silverman^[1]. Functional regression models are useful tools in functional data analysis, where one of the most interesting and challenging cases is function-on-function regression, see Ramsay and Silverman^[1,2], Yao et al.^{[3, 4]}. In this paper, we consider the following functional model proposed by Wang et al.^[5], for $m=1,\cdots,M$,

where $y_m(t)$ is the functional response, ${\boldsymbol{{z}}}_{{m}}(t)$ is a $p$-vector of functional covariates, ${\boldsymbol{{\nu}}}$ is the corresponding parameters, ${\boldsymbol{{x}}}_{{m}}(s, t)$ is a $q$-dimensional of covariates depends on $s$ and $t$, and ${\boldsymbol{{\beta}}}(s,t)$ is a vector of the functional coefficients, $S_t$ is interval for $t$, $\varepsilon_m(t)$ is random error term for the $m$th curve. Model (1) is flexible, and includes some function-on-function models in Gervini^[6], Malfait and Ramsay^[7], Ramsay and Silverman^[2], as special cases. Note that $\tau_m$ is used to model the heterogeneity among the different subjects, which depends on ${\boldsymbol{{z}}}_{{m}}(t)$, ${\boldsymbol{{x}}}_{{m}}(\cdot, t)$. Wang et al.^[5] considered the above random effects model using Gaussian process priors. More on Gaussian process priors in functional model^[8,9].

However, when there exist outliers in the observations, it is not robust to use the model based on Gaussian process priors, see e.g. Wang et al.^[10]. Then in order to overcome the influence of outliers, various forms of student t-process have been developed to model a heavy-tailed process, e.g. Yu et al.^[11], Zhang and Yeung^[12]. Shah et al.^[13] pointed out that the t-distribution under addition is not closed to maintain the good properties of Gaussian models. Thus, Wang et al.^[10] developed an extended t-process regression, which has the following advantages: ① it can maintain the good properties of Gaussian process; ② it has flexible forms, and contains model in Shah et al.^[13] as a special case; ③ it is robust. More general discussions on t-process can see Refs. [10,14].

In this paper, we consider a functional nonparametric random effects model with extended t-process priors, and propose an estimation procedure. The proposed method has 3 merits. ① It applies the extended t-process prior to model the heterogeneity of individual effect in the function-on-function regression model such that the model has robustness; ② A basis expansion smoothing method and a penalized likelihood method are developed to estimate the parameter in the fixed effect and covariance function of random effects, which leads to estimation of the smoothing function and prediction of the random effect; ③ Information consistency of the parameter estimation is obtained.

The remainder of the paper is organized as follows. In Section 2, we present the nonparametric random effects model using extended t-process priors, and develop prediction distribution and estimation procedure. In Section 3, we conduct simulation studies and a real data example to evaluate the performance of the proposed method. The conclusions are given in Section 4. All the proofs are given in Appendix.

2.   Main results

2.1   Extended t-process

Extended t-process proposed by Wang et al.^[10] is briefly introduced as follows. Let $f(\cdot)$, a real-valued random function from ${\cal{X}}$ to $R$, satisfy that

where ${\rm{GP}}(\cdot,\cdot)$ and ${\rm{IG}}(\cdot,\cdot)$ stand for Gaussian process and inverse gamma distribution respectively. Then $f$ follows an extended t-process (ETP), and can be denoted by $f \sim {\rm{ETP}}(v, \omega, h, k)$. We call $h(\cdot): {\cal{X}} \rightarrow R$ mean function and $k(\cdot, \cdot): {\cal{X}} \times {\cal{X}} \rightarrow R$ covariance kernel function. From the definition of ETP, we show that for any points ${\boldsymbol{{X}}}=\left({\boldsymbol{{x}}}_{1}, \cdots, {\boldsymbol{{x}}}_{n}\right)^{\top}$, we have

meaning that ${\boldsymbol{{f}}}_{{n}}$ has an extended multivariate t-distribution (EMTD) with the following density function,

where ${\boldsymbol{{h}}}_{n}=\left(h\left({\boldsymbol{{x}}}_{1}\right), \cdots, h\left({\boldsymbol{{x}}}_{n}\right)\right)^{\top}, {\boldsymbol{{K}}}_{n}=\left(k_{i j}\right)_{n \times n}$ and $k_{i j}=k\left({\boldsymbol{{x}}}_{i}, {\boldsymbol{{x}}}_{j}\right)$.

2.2   Function-on-function regression model with random effects

In model (1), the random effect $\tau_m$ depicts individual effect. Considering robustness against outliers, an ETP process prior is applied to $\tau_m$. This paper assumes that $\tau_m$ and $\varepsilon_m$ have a joint extended t-process,

where $\delta_{\varepsilon}(t, s)=I(t=s)$ and $I(\cdot)$ is an indicator function.

Note that the random effect $\tau_m$ relies on ${\boldsymbol{z_m}}(t)$ and ${\boldsymbol{x_m}}(\cdot,t)$, then following Wang et al.^[5], the kernel function $k$ is an expression as

where ${\boldsymbol{{u}}}_{{m}}(t) = ({\boldsymbol{{z}}}_{{m}}^{{\top}}(t), {\boldsymbol{{x}}}_{{m}}^{{\top}}(\cdot, t))^\top$, ${\boldsymbol{{z}}}_{{m}}(t) = (z_{m1}(t),\cdots,z_{mp}(t))^\top$ and ${\boldsymbol{{x}}}_{{m}}(s,t) = (x_{m1}(s, t),\cdots,x_{mq}(s, t))^\top$. Let ${\boldsymbol{\theta}}=(\theta_{10}, \theta_{11}, \cdots , \theta_{1 Q}, \theta_{21}, \cdots , \theta_{2 Q})^{\top}$ represent a set of hyper-parameters with $Q=p+q$, and $\|g(\cdot)\|_{\Lambda}$ be a $\Lambda$ norm of function $g$. A choice of $\|\cdot\|_{\Lambda}$ is the $L_2$ norm of a function, that is, $\|g(\cdot)\|_{\Lambda}=\int g(s)^2{\rm{d}}s$ is a $\Lambda$ norm of function $g$.

Let observations $\{y_{m i}=y_{m}(t_{i})$, $i=1, \cdots, n, m=1, \cdots, M\}$, ${\boldsymbol{{u}}}_{{m}}(t_i)=({\boldsymbol{{z}}}_{{m}}^{{\top}}(t_i),{\boldsymbol{{x}}}_{{m}}^{{\top}}(\cdot,t_i))^{\top}$, error term $\varepsilon_{m i}=\varepsilon_{m}\left(t_{i}\right)$, where $\{t_{i}\}$ are observed times. Assume that true values of ${\boldsymbol{{\nu}}}$, ${\boldsymbol{{\beta}}}$, $\tau_m$ in model (1) are ${\boldsymbol{{\nu}}}_{{0}}$, ${\boldsymbol{{\beta}}}_{{0}}$, $\tau_{0 m}$ respectively. From model (1), we further consider the following (true) data model:

This paper aims to develop methods to estimate $\nu_0$, ${\boldsymbol{{\beta}}}_{{0}}$, and predict $\tau_{0m}$.

2.3   Prediction

Denote $c_m(t)={\boldsymbol{{z}}}_{{m}}^{{\top}}(t){\boldsymbol{{\nu}}}+\int_{S_t} {\boldsymbol{{x}}}_{{m}}^{{\top}}(s,t){\boldsymbol{{\beta}}}(s,t){\rm{d}}s$. From model (1), we have the following results,

It follows that for the observed data, we have the conditional distributions,

where ${\boldsymbol{{y}}}_{{m}}=(y_{m}(t_{1}),\cdots,y_{m}(t_{n}))^{\top}$ are observations for the $m$th subject at points $\{t_{1},\cdots,t_{n}\}$, similarly, ${\boldsymbol{{\tau}}}_{{m}}=(\tau_{m}({\boldsymbol{{u}}}_{{m}}(t_{1})), \cdots, \tau_{m}({\boldsymbol{{u}}}_{{m}}(t_{n})))^{\top}$, ${\boldsymbol{{c}}}_{{m}}=(c_{m}(t_{1}),\cdots,c_{m}(t_{n}))^{\top}$, ${\boldsymbol{{K}}}_{{m}}=(k_{{\boldsymbol{{\theta}}}}({\boldsymbol{{u}}}_{{m}}(t_{i}),{\boldsymbol{{u}}}_{{m}}(t_{j})))_{n \times n}$, ${\boldsymbol{I}}$ is the identity matrix.

Denoted by the data set ${\cal{D}}=\{(y_{m}(t_{j}), {\boldsymbol{{u}}}_{{m}}(t_{j})): j=1, \cdots,n, m=1, \cdots, M\}$. Since that

we obtain the posterior distribution of ${\boldsymbol{\tau}}_{m}$, that is

where $v^*=v + n/2$, $w^*=w + n/2$,

and

For prediction, at a new data point $t^*$, we have

where ${\boldsymbol{{k}}}_{{mt}}=(k({\boldsymbol{{u}}}_{{m}}(t), {\boldsymbol{{u}}}_{{m}}(t_{1})), \cdots, k({\boldsymbol{{u}}}_{{m}}(t), {\boldsymbol{{u}}}_{{m}}(t_{n})))^{\top}$. It indicates that

where

Therefore, we can use posterior mean

to predict $y_m(t^*)$, denoted by $\hat{y}_m(t^*)$. And using

estimate variance of $\hat{y}_m(t^*)$.

Moreover, we show $\tau_{m} \mid {\cal{D}} \sim {\rm{EMTD}} ( v^*,w^*, \tilde{\mu}_{m}, \tilde{\sigma}_{m}),$ where for data points $u$, $v$,

Similarly,

It follows that Eq. (4) is an estimation of the covariance function of $\hat{y}_{m}(\cdot)$.

2.4   Parameter estimation

Note that ${\boldsymbol{{\beta}}}(s, t)$ in model (1) is a smooth function and can be approximated based on basis functions $\{\phi_k(s),k=1,\cdots,K_s\}$, and $\{\psi_k(s),k=1,\cdots,K_t\}$,

where $\{b_{i k l}\}$ are coefficients, ${\boldsymbol{B}}_{{i}}=(b_{i k l})_{K_{s} \times K_{t}}$, ${\boldsymbol{{\phi}}}(s)=\left(\phi_{1}(s), \cdots, \phi_{K_{s}}(s)\right)^{\top}$,${\boldsymbol{{\psi}}}(t)=\left(\psi_{1}(t), \cdots, \psi_{K_{t}}(t)\right)^{\top}$. Let ${\boldsymbol{\phi}}_{{x m i}}(t)=\int_{S_{t}} {\boldsymbol{{\phi}}}(s) x_{m i}(s, t) {\rm{d}} s$,and

where “$\otimes$” represents the Kronecker product. Hence, $c_{m}(t)={\boldsymbol{\gamma}}_{{m}}(t)^{\top} {\boldsymbol{b}}$. Let ${\boldsymbol{\varGamma}}_{{m n}}=({\boldsymbol{\gamma}}_{{m}}(t_{1}), \cdots, {\boldsymbol{\gamma}}_{{m}}(t_{n})),$ then ${\boldsymbol{{c}}}_{{m}}={\boldsymbol{\varGamma}}_{{m n}}^{{\top}} {\boldsymbol{{b}}}$.

Next we estimate ${\boldsymbol{\theta}}$, ${\boldsymbol{b}}$ and ${\boldsymbol{\sigma}}^{{2}}$ via using a likelihood method. By Eq. (3), we obtain a likelihood function of ${\boldsymbol{{y}}}_{{m}}$,

where $H_m({\boldsymbol{{\theta}}}, {\boldsymbol{{b}}}, \sigma^{2})=({\boldsymbol{{y}}}_{{m}}-{\boldsymbol{\varGamma}}_{{m n}}^{{\top}} {\boldsymbol{{b}}})^{\top}({\boldsymbol{{K}}}_{{m}}+\sigma^{2} {\boldsymbol{{I}}})^{-1}({\boldsymbol{{y}}}_{{m}}-{\boldsymbol{\varGamma}}_{{m n}}^{{\top}} {\boldsymbol{{b}}})$. Then we have the following objective function based on the log-likelihood function,

Due to the smoothness of ${\boldsymbol{{\beta}}}(\cdot,\cdot)$, following from Ramsay and Silverman^[5], we consider the following penalty functions,

where

Therefore, we develop an objective function,

where $\lambda_s$ and $\lambda_t$ are tuning parameters. Take the derivative of $G({\boldsymbol{{\theta}}}, {\boldsymbol{{b}}}, \sigma^{2})$ with respect to ${\boldsymbol{b}}$, we can obtain the estimation equation

where ${\boldsymbol{\varLambda}}={\rm diag}({\bf{0}}_{{p\times p}},\lambda_s{\boldsymbol{J}}_{{{\boldsymbol{{\psi}}} {\boldsymbol{{\psi}}}}}\otimes{\boldsymbol{L}}_{{{\boldsymbol{{\phi}}}{\boldsymbol{{\phi}}}}}+\lambda_t{\boldsymbol{L}}_{{{\boldsymbol{{\psi}}} {\boldsymbol{{\psi}}}}}\otimes{\boldsymbol{J}}_{{{\boldsymbol{{\phi}}} {\boldsymbol{{\phi}}}}},\cdots,\lambda_s{\boldsymbol{J}}_{{{\boldsymbol{{\psi}}} {\boldsymbol{{\psi}}}}}\otimes{\boldsymbol{L}}_{{{\boldsymbol{{\phi}}}{\boldsymbol{{\phi}}}}}+ \lambda_t{\boldsymbol{L}}_{{{\boldsymbol{{\psi}}} {\boldsymbol{{\psi}}}}}\otimes{\boldsymbol{J}}_{{{\boldsymbol{{\phi}}} {\boldsymbol{{\phi}}}}})$ is a $(p+qK_sK_t)\times (p+qK_sK_t)$ matrix. Similarly, we can get estimation equations with respect to ${\boldsymbol{{\theta}}}$ and $\sigma^{2}$.

From these estimation equations, we construct an estimation procedure as follows.

Step 1 Given an initial estimate of ${\boldsymbol{\theta}}$;

Step 2 Given ${\boldsymbol{\theta}}$, we update the estimates of ${\boldsymbol{b}}$ and $\sigma^2$ via

Step 3 Given ${\boldsymbol{b}}$ and $\sigma^2$, we update the estimate of ${\boldsymbol{\theta}}$ via

Step 4 Repeat Step 2 and Step 3 until convergence.

Similar to Ref. [5], when the absolute value of relative difference of $l({\boldsymbol{{\theta}}}, {\boldsymbol{{b}}}, \sigma^{2})$ between two successive iterations is less than a given value, the procedure stops.

2.5   Information consistency

The common mean structure and its properties have been studied a lot in functional models, see Yao et al.^[4], Yuan and Cai^[15], Sun et al.^[16], and among others. Next we only consider the information consistency. Let ${\cal {X}}={\cal {X}}_1 \times {\cal{X}}_2$, where ${\cal {X}}_1$ and ${\cal{X}}_2$ are spaces covariates ${\boldsymbol{z}}_{m}(t)$ and ${\boldsymbol{x}}_{m}(\cdot,t)$ belonging to. Let $p_{\sigma _{0}}({\boldsymbol{{y}}}_{{m}}|\tau_{0m},{\boldsymbol{{u}}}_{{m}})$ be the density function to generate the data ${\boldsymbol{{y}}}_{{m}}$ given ${\boldsymbol{{u}}}_{{m}}$ and $\tau_{0m}$, where $\sigma_{0}$ is the true value of $\sigma$, $\tau_{0m}$ is the true value of $\tau_m$. Let $p_{{\boldsymbol{{\theta}}}}(\tau)$ be a measurement of the random process $\tau$ on space $\cal{F}=\{\tau(\cdot,\cdot): {\cal{X}} \rightarrow R\}$. Let

be the density function to generate the data ${\boldsymbol{{y}}}_{{m}}$ given ${\boldsymbol{{u}}}_m$ under model (1). Let $p_{\sigma_0,\hat{{\boldsymbol{\theta }}}}({\boldsymbol{{y}}}_{{m}}|{\boldsymbol{{u}}}_{{m}})$ be the estimated density function. Denote

as the Kullback-Leibler divergence between two densities $p_{1}$ and $p_{2}$. According to Ref. [6], we only need to show the Kullback-Leibler divergence between two density functions for ${\boldsymbol{{y}}}_{{m}}|{\boldsymbol{{u}}}_{{m}}$ from the true and the assumed models tends to zero when $n$ is large enough.

For information consistency of the parameter estimation, we need the following condition.

Condition (A): $\|\tau_{0m}\|_k$ is bounded and

where $\|\tau_{0m}\|_k$ is the reproducing kernel Hilbert space norm of $\tau_{0m}$ associated with $k(\cdot , \cdot ;{\boldsymbol{{\theta}}})$, ${\boldsymbol{{K}}}_{{m}}$ is covariance matrix of $\tau_{0m}$ over ${\boldsymbol{{u}}}_{{m}}$, ${\boldsymbol{{I}}}$ is the $n \times n$ identity matrix.

More details about Condition (A) can see Seeger et al.^[17] and Wang et al.^[5]. More on reproducing kernel Hilbert space can see Berlinet and Thomas^[18].

Proposition 2.1. Under the conditions in Lemma A.1 (Appendix) and condition (A), we have

where the expectation is taken over the distribution of ${\boldsymbol{{u}}}_{{m}}$.

3.   Numerical results

3.1   Simulations

Performance of the proposed method is investigated by numerical studies. Simulation data are generated by the following model,

where ${z}_{m}(\cdot) \sim {\rm{GP}}(h_1, k_1)$, $h_1 = h_1(t) = t$, for $t \in (0,1)$, $k_1 = k_1({z}_{m}(t_1), {z}_{m}(t_2)) \;=\; g(t_1,t_2) \;=\; 0.1\exp\{-5(t_1-t_2)^2\} + 0.1t_1t_2$, and ${x}_{m} (\cdot,\cdot) \sim GP(h_2, k_2)$, $h_2 = h_2(t) = t + {\rm{cos}}(s)(s)$, for $t,s \in (0,1)$, $k_2 = k_2({x}_{m}(s_1,t),{x}_{m}(s_2,t)) = g(s_1,s_2).$ Let ${\boldsymbol{{\nu}}} = 1.0, \theta_{10} = \theta_{12} =\theta_{21} = \theta_{22} = 0.1, \theta_{11} = 10,$ $\sigma^2 = 0.5$, and $t$ and $s$ take 20 points equally in (0,1). Consider four different combinations of $\tau_{m}$ and ${\boldsymbol{{\beta}}}(s,t)$,

S1: $\tau_{m} \sim {\rm{GP}}(0, {\rm Cov}(\tau_{m}({\boldsymbol{{u}}}_{{m}}(t_{1})), \tau_{m}({\boldsymbol{{u}}}_{{m}}(t_{2}))))$, and $h_2 = h_2(t) = t +{\rm{cos}}(s) (s)$, for $s,t \in (0,1)$;

S2: $\tau_{m} \sim {\rm{GP}}(0, {\rm Cov}(\tau_{m}({\boldsymbol{{u}}}_{{m}}(t_{1})), \tau_{m}({\boldsymbol{{u}}}_{{m}}(t_{2}))))$, and ${\boldsymbol{{\beta}}}(s,t) \;=\; \exp \{-(t^2 + s^2)\}/10$, for $s,t \in (0,1)$;

S3: $\tau_{m} =0$ and ${\boldsymbol{{\beta}}}(s,t) = (t^2 + \cos(s))/10$, for $s,t \in (0,1)$;

S4: $\tau_{m} =0$ and ${\boldsymbol{{\beta}}}(s,t) = \exp\{-(t^2 + s^2)\}/10$, for $s,t \in (0,1)$.

We take sample sizes $M$=10, 20, and 30. All simulations are repeated 500 times.

To show robustness of model (1) with random effect having ETPR, saying ETPR, we also compute model (1) with random effect having GPR, denoted by GPR. Two indices: prediction error (PE),

and average estimation bias (AB)

are applied to show performance of two methods: ETPR and GPR, where $\hat{f}(t)={\boldsymbol{{z}}}_{{m}}^{{\top}}(t)\hat{{\boldsymbol{{\nu}}}}+\int_{0}^{1} {\boldsymbol{{x}}}_{{m}}^{{\top}}(s,t)\hat{{\boldsymbol{{\beta}}}}(s,t){\rm{d}}s+ \hat{\tau}_{m}({\boldsymbol{{z}}}_{{m}}(t), {\boldsymbol{{x}}}_{{m}}(\cdot, t))$ is an estimator of the true regression function $f_0(t)= {\boldsymbol{{z}}}_{{m}}^{{\top}}(t){\boldsymbol{{\nu}}}_0+\int_{0}^{1} {\boldsymbol{{x}}}_{{m}}^{{\top}}(s,t){\boldsymbol{{\beta}}}_0(s,t){\rm{d}}s+\tau_{0m}({\boldsymbol{{z}}}_{{m}}(t), {\boldsymbol{{x}}}_{{m}}(\cdot, t))$. To show robustness of our method, one curve is randomly selected and added with an extra disturbance, $\delta t_3$, where $t_3$ stands for student $t$ distribution with degree of freedom 3. Table 1 presents the values of PE and AB from these two methods. We see that ETPR has smaller PE and AB than GPR, especially with $\delta = 1.0$ and small sample sizes. It shows that the proposed method ETPR has more robustness against outliers compared to GPR.

Table  1.  PE and AB of prediction from ETPR method and GPR method, where SDs are presented in parentheses.

Setup	$\delta$	Method	$M=10$			$M=20$			$M=30$
			${\rm{PE}}$	${\rm{AB}}$		${\rm{PE}}$	${\rm{AB}}$		${\rm{PE}}$	${\rm{AB}}$
S1	0.5	ETPR	0.420(0.200)	0.170(0.194)		0.375(0.073)	0.124(0.065)		0.360(0.040)	0.111(0.029)
		GPP	0.430(0.212)	0.188(0.207)		0.385(0.079)	0.135(0.072)		0.37(0.044)	0.12(0.035)
	1	ETPR	0.597(0.718)	0.3484(0.719)		0.443(0.26)	0.192(0.26)		0.395(0.119)	0.146(0.117)
		GPP	0.642(0.851)	0.3932(0.85)		0.466(0.296)	0.216(0.297)		0.414(0.137)	0.165(0.135)
S2	0.5	ETPR	0.392(0.183)	0.143(0.179)		0.354(0.071)	0.103(0.063)		0.341(0.038)	0.092(0.027)
		GPP	0.405(0.212)	0.156(0.21)		0.355(0.078)	0.106(0.072)		0.341(0.042)	0.092(0.034)
	1	ETPR	0.569(0.707)	0.321(0.712)		0.42(0.273)	0.172(0.273)		0.377(0.119)	0.129(0.117)
		GPP	0.608(0.852)	0.36(0.855)		0.437(0.295)	0.187(0.295)		0.386(0.138)	0.137(0.137)
S3	0.5	ETPR	0.345(0.229)	0.094(0.225)		0.308(0.088)	0.057(0.083)		0.293(0.041)	0.042(0.034)
		GPP	0.359(0.289)	0.109(0.283)		0.319(0.122)	0.069(0.119)		0.301(0.061)	0.051(0.056)
	1	ETPR	0.515(0.817)	0.263(0.814)		0.38(0.359)	0.129(0.359)		0.3275(0.191)	0.077(0.191)
		GPP	0.579(1.158)	0.328(1.15)		0.425(0.483)	0.175(0.482)		0.3568(0.24)	0.106(0.24)
S4	0.5	ETPR	0.311(0.22)	0.06(0.218)		0.276(0.091)	0.026(0.087)		0.263(0.04)	0.013(0.035)
		GPP	0.326(0.287)	0.074(0.283)		0.287(0.123)	0.036(0.121)		0.269(0.06)	0.019(0.057)
	1	ETPR	0.481(0.82)	0.229(0.817)		0.351(0.363)	0.1(0.363)		0.299(0.193)	0.049(0.194)
		GPP	0.547(1.154)	0.295(1.15)		0.395(0.485)	0.144(0.485)		0.326(0.241)	0.076(0.243)

 | Show Table

DownLoad: CSV

In addition, we also consider one constant disturbance for the abnormal curves with small sample sizes 10 and 20. Tables 2 and 3 present PE and AB of prediction from ETPR method and GPR method for one and two curves disturbed, respectively. We see that ETPR has better performance in prediction compared to GPR.

Table  2.  PE and AB of prediction from ETPR method and GPR method with one curve disturbed by constant 1.0, where SDs are presented in parentheses.

Setup	Method	$M=10$			$M=20$
		PE	AB		PE	AB
S1	ETPR	0.429(0.069)	0.18(0.042)		0.388(0.044)	0.138(0.025)
	GPP	0.457(0.073)	0.208(0.044)		0.406(0.046)	0.155(0.026)
S2	ETPR	0.401(0.066)	0.153(0.039)		0.365(0.042)	0.114(0.023)
	GPP	0.42(0.069)	0.173(0.04)		0.374(0.043)	0.123(0.023)
S3	ETPR	0.341(0.053)	0.09(0.027)		0.302(0.032)	0.051(0.01)
	GPP	0.368(0.065)	0.117(0.044)		0.312(0.039)	0.061(0.021)
S4	ETPR	0.303(0.048)	0.052(0.024)		0.268(0.029)	0.018(0.009)
	GPP	0.33(0.065)	0.08(0.045)		0.275(0.037)	0.025(0.022)

 | Show Table

DownLoad: CSV

Table  3.  PE and AB of prediction from ETPR method and GPR method with two curves disturbed by constant 1.0, where SDs are presented in parentheses.

Setup	Method	$M=10$			$M=20$
		PE	AB		PE	AB
S1	ETPR	0.512(0.082)	0.263(0.053)		0.425(0.05)	0.175(0.029)
	GPP	0.558(0.083)	0.308(0.055)		0.453(0.051)	0.202(0.03)
S2	ETPR	0.479(0.079)	0.23(0.049)		0.399(0.048)	0.148(0.027)
	GPP	0.517(0.082)	0.268(0.052)		0.421(0.05)	0.169(0.028)
S3	ETPR	0.423(0.071)	0.171(0.045)		0.329(0.037)	0.078(0.018)
	GPP	0.48(0.081)	0.227(0.055)		0.367(0.047)	0.115(0.032)
S4	ETPR	0.381(0.068)	0.129(0.044)		0.296(0.035)	0.045(0.017)
	GPP	0.441(0.078)	0.188(0.053)		0.333(0.048)	0.081(0.034)

 | Show Table

DownLoad: CSV

3.2   Real data example

The proposed method is applied to Canadian weather data, which is obtained from the R package fda. We aim to study fixed effect of temperature on precipitation by common temperature effect of stations in the same region, and random effect of temperature on precipitation by individual effect of each station. Generally, the 35 stations are divided into four regions: Arctic, Atlantic, Pacific and Continental. Obviously, there exists heterogeneity among the stations due to the spatial nature of the weather data. Then we propose the following model:

where ${y}_{ij}(t)$ represents precipitation and ${x}_{ij}(t)$ represents temperature, for time $t$, region $i$ and $j$th station. In this model, we have ${z}_{ij}(t) = 1$ and ${x}_{ij}(s,t) = {x}_{ij}(s)$ which effectively simplifies model fit.

Figs. 1 and 2 show random and fixed effects of the 4 regions: Arctic, Atlantic, Pacific and Continental from the proposed method. We see from the random effects that each station in the same region has different temperature effects on the precipitation. To compare performance of prediction from ETPR with GPR, 10-folds cross validation method is used to compute mean squares of prediction errors, 0.310 and 0.314, for ETPR and GPR, respectively. It shows that ETPR has a little better performance in prediction.

Figure  1.  Random and fixed effects of model using ETPR for Arctic and Atlantic.

DownLoad: Full-Size Img PowerPoint

Figure  2.  Random and fixed effects of model using ETPR for Continental and Pacific.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

A function-on-function random effects model with extended t-process prior in this paper is developed to analyze functional data which may include outliers. The proposed model is flexible, including various kinds of functional models, such as the function-on-function linear model^[2] and the historical functional regression model^[7] as special cases. The proposed extended t-process model is not only robust against outliers, but also inherits almost all the nice properties from Gaussian process regression, such as closed form of prediction and convenient computation procedure. The estimation procedure and computing algorithm are developed to estimate the parameters and predict the random effect in the regression model. The functional response considered in this paper has one dimension. In practical application, functional multi-response may consist of several correlated curves. It is interesting that the proposed method is extended to functional data with multi-response, which will be studied in our further work.

Appendix

Lemma A.1. Let $w=v-1$. Under model (1), assume that ${\boldsymbol{{y}}}_{{m}}$ are independently sampled, the covariance kernel function $k$ is bounded and continuous on the parameter ${\boldsymbol{{\theta}}}$, and $\hat{{\boldsymbol{\theta}}}$ converges to ${\boldsymbol{{\theta}}}$ when $n \rightarrow \infty$. Then, for a positive constant $c$ and any $\varepsilon>0$, when $n$ is large enough, we have

where $q_m^{2}=({\boldsymbol{{y}}}_m-{\boldsymbol{{c}}}_{{0m}}-{\boldsymbol{{\tau}}}_{{0m}})^{\top}({\boldsymbol{{y}}}_{{m}}-{\boldsymbol{{c}}}_{{0m}}-{\boldsymbol{{\tau}}}_{{0m}}) / \sigma_{0}^2$, ${\boldsymbol{{c}}}_{{0m}}$ is the true value of ${\boldsymbol{{c}}}_{{m}}$, $\|\tau_{0m}\|_k$ is the reproducing kernel Hilbert space norm of $\tau_{0m}$ associated with $k(\cdot , \cdot ;{\boldsymbol{{\theta}}})$, ${\boldsymbol{{K}}}_{{m}}$ is covariance matrix of $\tau_{0m}$ over ${\boldsymbol{{u}}}_{{m}}$, ${\boldsymbol{{I}}}$ is the $n \times n$ identity matrix.

Proof of Lemma A.1. Assume $r$ is a random variable following inverse gamma distribution ${\rm{IG}}(v,(v-1))$. Conditional on $r$, we have

where ${\rm{GP}}(h,k)$ stands for Gaussian process with mean function $h$ and covariance function $k$. Then conditional on $r_m$, the extended t-process regression model $y_m=c_m+\tau_m+\varepsilon_m$ becomes Gaussian process regression model

where $\tilde{\tau}_m=\tau_m|r_m \sim {\rm{G P}}(0, r_m k(\cdot, \cdot; {\boldsymbol{\theta}}))$, $\tilde{\varepsilon}_m=\varepsilon_m| r_m \sim {\rm{G P}}(0, r \sigma^2\delta_{\varepsilon})$, and $\tilde{\tau}_m$ and $\tilde{\varepsilon}_m$ are independent. Denoted the computation of conditional probability density for given $r_m$ by $\tilde{p}$. Let

where $\tilde{p}_{\boldsymbol{\theta}}$ is the induced measure from Gaussian process ${\rm{G P}}(0, r_m k(\cdot,\cdot ; \hat{\boldsymbol{\theta}}))$. Note that variable $r$ is independent of ${\boldsymbol{{u}}}_{{m}}$. We can show that

By similar procedures in Seeger et al.^[17] and Wang et al.^[10], for any given $r$, we have

then it follows that

Let $g^{*}(r)$ be the density function of ${\rm{IG}}(v+n / 2,v-1+ q_m^{2} / 2)$. It easily shows that

We have

which shows that Lemma A.1 holds.

Proof of Proposition 2.1. Obviously $q_m^{2}=({\boldsymbol{{y}}}_m-{\boldsymbol{{c}}}_{{0m}}-{\boldsymbol{{\tau}}}_{{0m}})^{\top} \cdot ({\boldsymbol{{y}}}_m-{\boldsymbol{{c}}}_{{0m}}-{\boldsymbol{{\tau}}}_{{0m}}) / \sigma_{0}^2=O(n)$. Under the conditions of Lemma A.1 and condition (A), by Lemma A.1, for a positive constant $c$ and any $\varepsilon>0$, when $n$ is large enough, we have

Thus, it completes the proof.

Acknowledgements

We thank the reviewers for their insightful comments and suggestions. This work was supported in part by the National Natural Science Foundation of China (11971457), Anhui Provincial Natural Science Foundation (1908085MA06) and the Fundamental Research Funds for the Central Universities (WK2040000035).

Conflict of interest

The authors declare that they have no conflict of interest.