Robust functional ANOVA model with t-process

TL;DR

This paper introduces a robust functional ANOVA model using t-processes to improve estimation and prediction in the presence of outliers, with proven statistical properties and demonstrated effectiveness through simulations and real data.

Contribution

It develops a novel robust estimation method for functional ANOVA models utilizing t-processes, enhancing robustness and prediction accuracy.

Findings

Method performs well in simulations

Effective in real data applications

Ensures robustness and consistency

Abstract

Robust estimation approaches are of fundamental importance for statistical modelling. To reduce susceptibility to outliers, we propose a robust estimation procedure with t-process under functional ANOVA model. Besides common mean structure of the studied subjects, their personal characters are also informative, especially for prediction. We develop a prediction method to predict the individual effect. Statistical properties, such as robustness and information consistency, are studied. Numerical studies including simulation and real data examples show that the proposed method performs well.

Tables

Table 1. Table 4: Prediction errors of Canada precipitation data, market penetrations of the freezer and Hi-Fi stereo from TP, GP and TP( τ = 0 𝜏 0 \tau=0 ), where standard deviation is in parentheses.

Data	TP	GP	TP($\tau =0$)
Precipitation	0.0416(0.0029)	0.0488(0.0046)	0.0490(0.0045)
Freezer	0.1362(0.0060)	0.1480(0.0057)	0.1487(0.0055)
Hi-Fi stereo	0.0583(0.0149)	0.0641(0.0144)	0.0645(0.0144)

Equations

$$y_{ij}(t)=\mu (t)+{\alpha }_i(t)+{\tau }_{ij}({\mathbf{u}}_{ij}(t))+{\epsilon }_{ij}(t),j=1,\dots ,J_i,i=1,\dots ,I,$$

$$\sum _{i=1}^I{\alpha }_i(\cdot )=0.$$

$$f(\mathbf{z})=|2\pi (\nu -1){\sigma }^2{\mathbf{I}}_n{|}^{-\frac{1}{2}}\frac{\Gamma (n/2+\nu )}{\Gamma (\nu )}{\left(1+\frac{||\mathbf{z}-\mathbf{h}|{|}^2}{2(\nu -1){\sigma }^2}\right)}^{-(n/2+\nu )}.$$

$$y_{ij}(t)={\mathbf{z}}_{ij}^{\top }\mathbf{\beta }(t)+{\tau }_{ij}(t)+{\epsilon }_{ij}(t).$$

$$y_{ij}|{\tau }_{ij},{\sigma }^2,\nu$$

$${\tau }_{ij}$$

$$E(y_{ij}(t^{\ast })|{\mathcal{D}}_{\mathbf{n}})$$

$$=\int \{{\mathbf{z}}_{ij}^{\top }\mathbf{\beta }(t)+{\tau }_{ij}(t^{\ast })\}p({\tau }_{ij}(t^{\ast })|{\mathcal{D}}_{\mathbf{n}})d{\tau }_{ij}(t^{\ast }).$$

$$p({\tau }_{ij}(t^{\ast })|{\mathcal{D}}_{\mathbf{n}})$$

$$\approx \int p({\tau }_{ij}(t^{\ast })|{\mathbf{\tau }}_{ij})p_G({\mathbf{\tau }}_{ij}|{\mathcal{D}}_{\mathbf{n}})d{\mathbf{\tau }}_{ij},$$

$$\left(\begin{array}{c}{\mathbf{\tau }}_{ij}\\ {\tau }_{ij}(t^{\ast })\end{array}\right)\sim N\left(\begin{array}{c}0,\left(\begin{array}{cc}{\mathbf{K}}_{ij}& {\mathbf{k}}_{ij}(t^{\ast },\mathbf{t})\\ {\mathbf{k}}_{ij}(t^{\ast },\mathbf{t}{)}^{\top }& k_{ij}(t^{\ast },t^{\ast })\end{array}\right)\end{array}\right),$$

$$p_G({\mathbf{\tau }}_{ij}|{\mathcal{D}}_{\mathbf{n}})=N({\hat{\mathbf{\tau }}}_{ij},\mathbf{\Omega }),\mathbf{\Omega }\triangleq ({\mathbf{K}}_{ij}^{-1}+{\mathbf{D}}_m{)}^{-1},$$

$${\hat{y}}_{ij}(t^{\ast })=E(y_{ij}(t^{\ast })|{\mathcal{D}}_{\mathbf{n}})={\mathbf{z}}_{ij}^{\top }\mathbf{\beta }(t)+{\mathbf{a}}^{\top }{\hat{\mathbf{\tau }}}_{ij}.$$

$$\text{var}(y_{ij}(t^{\ast })|{\mathcal{D}}_{\mathbf{n}})$$

$$=$$

$$\text{var}(E(y_{ij}(t^{\ast })|{\tau }_{ij}(t^{\ast }))|{\mathcal{D}}_{\mathbf{n}}))={\mathbf{a}}^{\top }\mathbf{\Omega }\mathbf{a}+{\sigma }^{\ast 2},$$

$$E(\text{var}(y_{ij}(t^{\ast })|{\tau }_{ij}(t^{\ast }))|{\mathcal{D}}_{\mathbf{n}}))={\sigma }^2.$$

$$k_{ij}(t_1,t_2)=k(t_1,t_2;{\mathbf{\theta }}_{ij})$$

$$=$$

$${\mathbf{Y}}_{ij}|{\mathbf{\tau }}_{ij},\mathbf{B},{\sigma }^2,\nu$$

$${\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij}$$

$$\text{log }p({\mathcal{D}}_{\mathbf{n}}|\mathbf{\theta },\mathbf{B},\nu )$$

$$=\sum _{i,j}\text{log}\int p({\mathbf{Y}}_{ij}|{\mathbf{\tau }}_{ij},\mathbf{B},{\sigma }^2,\nu )p({\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij})d{\mathbf{\tau }}_{ij}.$$

$$p({\mathcal{D}}_{\mathbf{n}}|\mathbf{\theta },\mathbf{B},\nu )\dot{=}\prod _{i,j}{\frac{p({\mathbf{Y}}_{ij},{\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij},\mathbf{B},{\sigma }^2,\nu )}{p_G({\mathbf{\tau }}_{ij}|{\mathbf{Y}}_{ij},{\mathbf{\theta }}_{ij},\mathbf{B},{\sigma }^2,\nu )}|}_{{\mathbf{\tau }}_{ij}={\hat{\mathbf{\tau }}}_{ij}},$$

$$p({\mathbf{Y}}_{ij},{\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij},\mathbf{B},{\sigma }^2,\nu )=\text{exp}\{\text{log}(p({\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij}))+\sum _{k=1}^n\text{log}(p(y_{ij}(t_k)|{\tau }_{ijk},\mathbf{B},{\sigma }^2,\nu ))\}.$$

$$g_{ijk}({\tau }_{ijk})\approx g_{ijk}({\tau }_{ijk}^{(0)})+{\alpha }_{ijk}{\tau }_{ijk}-\frac{1}{2}{d}_{ijk}{\tau }_{ijk}^2$$

$${\alpha }_{ijk}=g_{ijk}^{{}^{\prime }}({\tau }_{ijk}^{(0)})-g_{ijk}^{{}^{\prime \prime }}({\tau }_{ijk}^{(0)}){\tau }_{ijk}^{(0)},d_{ijk}=-g_{ijk}^{{}^{\prime \prime }}({\tau }_{ijk}^{(0)}).$$

$$g_{ijk}^{{}^{\prime }}({\tau }_{ijk})=\frac{(1+2\nu )S_{ijk}}{2(\nu -1){\sigma }^2+S_{ijk}^2},g_{ijk}^{{}^{\prime \prime }}({\tau }_{ijk})=\frac{(1+2\nu )(S_{ijk}^2-2(\nu -1){\sigma }^2)}{(2(\nu -1){\sigma }^2+S_{ijk}^2{)}^2}.$$

$$p({\mathbf{Y}}_{ij},{\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij},\mathbf{B},{\sigma }^2,\nu )\propto \text{exp}\left\{-\frac{1}{2}{\mathbf{\tau }}_{ij}^{\top }{\mathbf{K}}_{ij}^{-1}{\mathbf{\tau }}_{ij}-\frac{{\mathbf{\tau }}_{ij}^{\top }{\mathbf{D}}_{ij}{\mathbf{\tau }}_{ij}}{2}+{\mathbf{\alpha }}_{ij}^{\top }{\mathbf{\tau }}_{ij}\right\},$$

$$l(\mathbf{B},\mathbf{\theta },\nu )=\sum _{m=1}^M\text{log}\int p({\mathbf{Y}}_{ij}|{\mathbf{\tau }}_{ij},\mathbf{B},{\sigma }^2,\nu )p({\mathbf{\tau }}_{ij}|{\mathbf{\theta }}_{ij})d{\mathbf{\tau }}_{ij}+{\text{PEN}}_{\lambda }(\mathbf{\beta }),$$

$${\text{PEN}}_{\lambda }(\mathbf{\beta })$$

$$=\lambda \sum _{i=1}^{I+1}\text{trace}[(({\mathbf{B}}^{\top }{)}_i{)}^{\top }({\mathbf{B}}^{\top }{)}_i{L}_{\Phi \Phi }],$$

$$p_{TP}({\mathbf{Y}}_{ij})={\int }_{\mathcal{F}}p(y_{ij}(t_1),...,y_{ij}(t_n)|{\tau }_{ij})d{p}_{\theta }({\tau }_{ij}),$$

$$p_0({\mathbf{Y}}_{ij})=p(y_{ij}(t_1),...,y_{ij}(t_n)|{\tau }_0),$$

$$\frac{1}{n}{E}_{{\mathbf{U}}_n}(D[p_0({\mathbf{Y}}_{ij}),p_{TP}({\mathbf{Y}}_{ij})])\to 0asn\to \infty ,$$

Full text

††journal: Computational Statistics & Data Analysis

1 Introduction

The analysis of variance (ANOVA) functional models are frequently applied to analyze data in many fields, such as econometrics, sociology and medical follow-up studies, such as application to the breast cancer data (Huang et al., 2000), the angles data of the hip over walking cycles (Rice and Wu, 2001), the vehicle road load data acquisition computer model (Dorin, 2010) and so on.

Functional ANOVA model is an adaptation of linear model or analysis of variance to analyze functional observations, which are also widely studied, such as Woods et al. (1996), Huang (1998), Brumback and Rice (1998) and therein reference. To the best of our knowledge, functional ANOVA models in literature rarely consider individual effect of the studied subject (each observed functional curve). However, for market penetration of new products (Ashish et al., 2009), data are collected from 86 countries in 7 regions. Market penetrations from different countries, such as the developing and developed countries, have dissimilar growth trend, because of discrepant economical levels. Hence, market penetration of each country has personal character besides the common mean structure. This paper not only considers the main structure of market penetration, but also focuses on individual effect of each country via a nonparametric random effect. Without loss of generality, we consider a one-way functional ANOVA concurrent regression model,

[TABLE]

where $t$ is the observed time, $y_{ij}(t)$ is a functional response for batch $j$ under level $i$, $\mu(t)$ is a grand function, $\alpha_{i}(t)$ measures the main effect of level $i$, $\tau_{ij}$ is a nonparametric random effect depicting personal character for the $j$-th subject at level $i$, $\bm{{u}}_{ij}(t)$ is a functional covariate with dimension $p$, and $\varepsilon_{ij}(t)$ is an error term. For identification of $\alpha_{i}(\cdot)$, it needs

[TABLE]

When random effect $\tau_{ij}=0$ and $\varepsilon_{ij}(t)$ has Gaussian distribution, model (1.1) becomes the common one way functional ANOVA model. For example, Brumback and Rice (1998) and Huang et al. (2000), respectively, used smoothing spline model and B-splines method to build this functional ANOVA model; Rice and Wu (2001) studied this model based on unequally sampled noisy curves.

When $\tau_{ij}$ follows Gaussian process (GP) with mean function 0 and kernel function $k_{ij}$, denoted by $GP(0,k_{ij})$ and $\varepsilon_{ij}(t)$ has Gaussian distribution, we denote model (1.1) by Gaussian process ANOVA model, saying GP model. Section 5.5 in Shi and Choi (2011) mentions the GP model. However, the GP is sensitive to outliers. For robustness, this paper assumes $\varepsilon_{ij}$ has heavy-tail process, such as t-process. For theoretical studies and application of heavy-tail processes, please see Yu et al. (2007), Zhang and Yeung (2010), Shah et al. (2014) and Wang et al. (2017) and so on. To illustrate the proposed method, this paper employs an extended t-process (ETP) (Wang et al., 2017) to depict the error term. Our method can be directly extended to other t-process.

This paper combines t-process error $\varepsilon_{ij}$ with Gaussian process random effect $\tau_{ij}$ in model (1.1), which leads to a robust functional ANOVA model. Independence of random effect and error term in this paper makes computation procedure more complicated compared to Wang et al. (2017) which uses a joint ETP process to build their model. We propose a computation method based on Gaussian approximation. Prediction of random effect is also an important for studying personal character. We develop a prediction method for functional ANOVA model. Statistical properties of the estimation method, such as robustness and information consistence, are investigated. Numerical studies including simulation and real examples show that the proposed method has robustness against outliers.

The remainder of this paper is organized as follows. Section 2 proposes a robust functional ANOVA model with t-process, and builds prediction and parameter estimation procedure. Asymptotic properties of estimation method are discussed in Section 3. Simulation studies and real examples are reported in Section 4, followed by the conclusion in Section 5. All the proofs are in Appendix.

2 Robust functional ANOVA model

Considering model (1.1), we assume that $\tau_{ij}(\bm{{u}}_{ij}(\cdot))\sim GP(0,k_{ij})$, where $k_{ij}$ is a kernel function. Hereafter, for convenience of notations, let $\tau_{ij}(\cdot)=\tau_{ij}(\bm{{u}}_{ij}(\cdot))$. For robustness to outlier, let $\varepsilon_{ij}$ be from the extended t-process $ETP(\nu,\nu-1,0,\sigma^{2}k_{\varepsilon})$ with $k_{\varepsilon}(u,v)=I(u=v)$ for $u,v\in R$ (Wang et al., 2017), where evaluated values of $\bm{{z}}=(\varepsilon_{ij}(t_{1}),...,\varepsilon_{ij}(t_{n}))^{\top}$ at any data points $t_{1},\ldots,t_{n}$, follows an extended multivariate t-distribution (EMTD) with the density function

[TABLE]

Denoted by $\bm{{\beta}}(t)=(\mu(t),\alpha_{1}(t),\ldots,\alpha_{I}(t))^{\top}$, $\bm{{z}}_{ij}=(1,0,\ldots,0,1,0,\ldots,0)^{\top}$ where the first and $(i+1)$-th elements are 1 and the others are 0. Model (1.1) is rewritten as

[TABLE]

From this model, it follows that

[TABLE]

2.1 Prediction

Let $\bm{{Y}}_{ij}=(y_{ij}(t_{1}),\ldots,y_{ij}(t_{n}))^{\top}$ be observed data for the $j-$th curve under level $i$, and $\bm{{\tau}}_{ij}=(\tau_{ij}(t_{1}),\ldots,\tau_{ij}(t_{n}))^{\top}$ at the observed times $\bm{{t}}=\{t_{1},\ldots,t_{n}\}$. Denoted by $\bm{{\mathcal{D}_{n}}}=\{\bm{{Y}}_{ij},j=1,\ldots,J_{i},i=1,...,I\}$. For a new point $t^{*}$, $E(y_{ij}(t^{*})|\bm{{\mathcal{D}_{n}}})$ is used to predict $y_{ij}$ at point $t^{*}$, denoted by $\hat{y}_{ij}(t^{*})$. From model (2.2), we show that

[TABLE]

Since the conditional distribution of $\tau_{ij}(t^{*})$ for given $\bm{{\mathcal{D}_{n}}}$ is involved in intractable integration, (2.3) is complicated to be computed. Gaussian approximation method is employed to calculate (2.3). It follows that

[TABLE]

where $p_{G}(\bm{{\tau}}_{ij}|\bm{{\mathcal{D}_{n}}})$ is the Gaussian approximation of $p(\bm{{\tau}}_{ij}|\bm{{\mathcal{D}_{n}}})$ (details seen in Section 2.2).

Since $\tau_{ij}\sim GP(0,k_{ij})$, we have

[TABLE]

where $\bm{{K}}_{ij}=(k_{ij}(t_{k},t_{l}))_{n\times n}$ is the covariance matrix of $\bm{{\tau}}_{ij}$, $\bm{{k}}_{ij}(t^{*},\bm{{t}})=(k_{ij}(t^{*},t_{1}),...,k_{ij}(t^{*},t_{n}))^{\top}$. It follows that $p(\tau_{ij}(t^{*})|\bm{{\tau}}_{ij})$ is normal distribution with mean $\bm{{a}}^{\top}\bm{{\tau}}_{ij}=\bm{{k}}_{ij}(t^{*},\bm{{t}})^{\top}\bm{{K}}_{ij}^{-1}\bm{{\tau}}_{ij}$ and variance $\sigma^{*2}=k_{ij}(t^{*},t^{*})-\bm{{k}}_{ij}(t^{*},\bm{{t}})^{\top}\bm{{K}}_{ij}^{-1}\bm{{k}}_{ij}(t^{*},\bm{{t}})$. From Section 2.2, we have

[TABLE]

where $\hat{\bm{{\tau}}}_{ij}$ and $\bm{{D}}_{m}$ are defined in the next subsection.

It easily shows that $p(\tau_{ij}(t^{*})|\bm{{\mathcal{D}_{n}}})$ has a normal density function $N(\bm{{a}}^{\top}\hat{\bm{{\tau}}}_{ij},\\ \bm{{a}}^{\top}\bm{{\Omega}}\bm{{a}}+\sigma^{*2}),$ which indicates that ${y}_{ij}(t^{*})$ can be estimated by

[TABLE]

We know that

[TABLE]

Easily, we show

[TABLE]

Hence, $\text{var}(y_{ij}(t^{*})|\bm{{\mathcal{D}_{n}}})=\bm{{a}}^{\top}\bm{{\Omega}}\bm{{a}}+\sigma^{*2}+\sigma^{2}$ is used to estimate the variance of $\hat{y}_{ij}(t^{*})$.

2.2 Parameter estimation

To apply the proposed method, it is necessary to estimate the unknown kernel function $k_{ij}(\cdot,\cdot)$, $\sigma^{2}$ and the function parameter $\bm{{\beta}}$. For kernel function, a function family such as a squared exponential kernel and Matérn class kernel can be applied (Shi and Choi, 2011). This paper takes a combination of a square exponential kernel and a non-stationary linear kernel,

[TABLE]

where $\bm{{\theta}}_{ij}=(\theta_{ij0},\theta_{ij1},...,\theta_{ijp},\eta_{ij1},...,\eta_{ijp})^{\top}$ is hyper-parameter. We use B-spline to approximate the functional parameter $\bm{{\beta}}(t)$, that is $\bm{{\beta}}(t)=\bm{{B}}^{\top}\bm{{\Phi}}(t)$, where $\bm{{\Phi}}(t)=(\phi_{1}(t),\ldots,\phi_{L}(t))^{\top},\{\phi_{k}(t),k=1,\ldots,L\}$ is a set of known basis functions and $\bm{{B}}$ is an unknown coefficient matrix with dimension $L\times(I+1)$.

Let $\bm{{\theta}}=(\bm{{\theta}}_{11}^{\top},\ldots,\bm{{\theta}}_{IJ_{I}}^{\top},\sigma^{2})^{\top},$ $\bm{{\Phi}}_{n}=(\bm{{\Phi}}(t_{1}),\ldots,\bm{{\Phi}}(t_{n}))_{L\times n}$, $\bm{{X}}_{ij}=(1,0,\ldots,$ $0,1,0,\ldots,0)^{\top}$ with the first and $(i+1)$-th elements are 1 and the others are 0. Then, we have

[TABLE]

where $\otimes$ represents the Kronecker product. Then a marginal log-likelihood is

[TABLE]

However, the marginal log-likelihood can not be computed directly, because of intractable integration. To avoid the complicated integration, some methods, such as Laplace approximation and the Markov Chain Monte Carlo(MCMC) method, are used to approximate the marginal log-likelihood. But the error rate of the Laplace approximation may be $O(1)$ as the dimension of $\bm{{\tau}}_{ij}$ increases with the sample size $n$ (Rue et al., 2009), and the MCMC method is too complicated to solve the problem conveniently, see Section 7.1.4 in Shi and Choi (2011). This paper uses a Gaussian approximation method to compute $p(\bm{{\mathcal{D}_{n}}}|\bm{{\theta}},\bm{{B}},\nu)$ as follows

[TABLE]

where $p_{G}(\bm{{\tau}}_{ij}|\bm{{Y}}_{ij},\bm{{\theta}}_{ij},\bm{{B}},\sigma^{2},\nu)$ is the Gaussian approximation to the full conditional density $p(\bm{{\tau}}_{ij}|\bm{{Y}}_{ij},\bm{{\theta}}_{ij},\bm{{B}},\sigma^{2},\nu)$ and $\hat{\bm{{\tau}}}_{ij}$ is the mode of the full conditional density for given $\bm{{\theta}}_{ij}$, $\bm{{B}},~{}\sigma^{2}$ and $\nu$. Let $\tau_{ijk}=\tau_{ij}(t_{k})$. We have

[TABLE]

From Taylor expansion, we approximate $g_{ijk}(\tau_{ijk})\triangleq\text{log}(p(y_{ij}(t_{k})|\tau_{ijk},\bm{{B}},\sigma^{2},\nu))$ to the second order

[TABLE]

where

[TABLE]

Let $S_{ijk}=y_{ij}(t_{k})-(\bm{{\Phi}}^{\top}(t_{i})\otimes\bm{{X}}_{ij}^{\top})\text{Vec}(\bm{{B}}^{\top})-\tau_{ijk}.$ From density function of EMTD, it gives

[TABLE]

Therefore,

[TABLE]

where $\bm{{D}}_{ij}=\text{diag}(d_{ij1},\ldots,d_{ijn})$ and $\bm{{\alpha}}_{ij}^{\top}=(\alpha_{ij1},\ldots,\alpha_{ijn})$.

Next, for some initial value $\bm{{\tau}}_{ij}^{(0)}$, we can use Fisher scoring algorithm to find the Gaussian approximation, and the $k$th iteration is given by

(i) Find the solution $\bm{{\tau}}_{ij}^{(k)}$ from $(\bm{{K}}_{ij}^{-1}+\bm{{D}}_{ij})\bm{{\tau}}_{ij}^{(k)}=\bm{{\alpha}}_{ij}$;

(ii) Update $\bm{{\alpha}}_{ij}$ and $\bm{{D}}_{ij}$ using $\bm{{\tau}}_{ij}^{(k)}$ and repeat (i).

After the process converges, say at $\hat{\bm{{\tau}}}_{ij}$, we can get the Gaussian approximation $p_{G}(\bm{{\tau}}_{ij}|\bm{{Y}}_{ij},\bm{{\theta}}_{ij},\bm{{B}},\sigma^{2},\nu)$ which is the density function of the normal distribution $MVN(\hat{\bm{{\tau}}}_{ij},(\bm{{K}}_{ij}^{-1}+\bm{{D}}_{ij})^{-1})$.

Considering smoothness of $\bm{{\beta}}(\cdot)$, we need to add a penalty $\text{PEN}_{\lambda}(\bm{{\beta}})$ into the log-likelihood function, denoted by

[TABLE]

where $\lambda$ is a tuning parameter. This paper takes the penalized function,

[TABLE]

where $L_{\Phi\Phi}=\int_{a}^{b}{[L\bm{{\Phi}}(t)][L\bm{{\Phi}}(t)^{\top}]dt}$ and $L$ is a derivative operator. The form of $L$ depends on the property of function $\bm{{\beta}}(t)$, such as for the second order continuously derivative function, $L\bm{{\Phi}}(t)=\omega_{0}\bm{{\Phi}}(t)+\omega_{1}D\bm{{\Phi}}(t)+D^{2}\bm{{\Phi}}(t)$, where the weights $\omega_{0},\omega_{1}$ may be either constant or function $\omega_{0}(t),\omega_{1}(t)$, and $D\bm{{\Phi}}(t),D^{2}\bm{{\Phi}}(t)$ mean the first and second derivative of $\bm{{\Phi}}(t)$.

From (2.6), (2.7) and (2.8),we obtain an estimation procedure for parameters $\bm{{\theta}}$, $\bm{{B}}$ and $\nu$. For initial values $\bm{{\theta}}=\bm{{\theta}}_{0}$, $\bm{{B}}=\bm{{B}}_{0}$ and $\nu=\nu_{0}$,

(I) For given $\bm{{\theta}}$, $\bm{{B}}$ and $\nu$, obtain$\hat{\bm{{\tau}}}=(\hat{\bm{{\tau}}}_{1n},\ldots,\hat{\bm{{\tau}}}_{Mn})$ by the Fisher scoring algorithm above-mentioned.

(II) For given $\bm{{\tau}}=\hat{\bm{{\tau}}}$, we update estimates of $\bm{{\theta}}$, $\bm{{B}}^{\top}$ and $\nu$ by maximizing the penalized log-likelihood (2.7) with (2.6) and (2.8).

(III) Repeat (I)-(II) until some convergence criteria hold.

3 Asymptotic properties

3.1 Robustness

Influence function is an important index to measure the robustness of estimation method to outliers. It is known that Gaussian distribution to fit data tends to be susceptive to outliers and has unbounded influence function, while t-distribution has heavy-tail and bound influence function. This paper shows that the proposed functional ANOVA models with t-process error (TP model) is more robust than those with Gaussian process error (GP model), see the next theorem (its proof seen in Appendix A).

Theorem 1. *If kernel function $k_{ij}$ is bounded for the $j$-th curve under level $i$, $j=1,...,J_{i},~{}i=1,...,I$, and continuously differentiable on $\bm{{\theta}}_{ij}$, then for given $\nu$, the estimators $\hat{\bm{{B}}}$ and $\hat{\bm{{\theta}}}$ from the TP model have bound influence functions, while that from the GP model does not. *

3.2 Information consistency

The consistency related to the estimation approach for model (2.2) involves two aspects: the common mean parameter function $\bm{{\beta}}(t)$ and the random effect curve $\tau_{ij}(t)$ (or the curve $y_{ij}(t)$). The common mean function is estimated by collected data from all curves and its consistency is already proved in many different functional linear models under suitable conditions, see Yao et al. (2005), Ramsay and Silverman (1997). Therefore, this paper focuses on the second issue, the information consistency of $\hat{y}_{ij}(\cdot)$ to $y_{ij}(\cdot)$. Information consistency of the estimation method is also an important issue. For Gaussian process regression model and the extended t-process regression model, this properties are obtained in Seeger et al. (2008), Wang and Shi (2014) and Wang et al. (2017). Here, we show this property for model (2.2).

Without loss of generality, suppose the mean function about the $j-$th curve under level $i$, denoted by $\bm{{\mu}}_{n}\triangleq(\bm{{X}}_{ij}^{\top}\otimes\bm{{\Phi}}_{n}^{\top})\text{Vec}(\bm{{B}})$, has already been known. Let covariate values $\bm{{U}}_{n}=\{\bm{{u}}_{ij}(t_{1}),\ldots,\bm{{u}}_{ij}(t_{n})\}$ where $\bm{{u}}_{ij}(t_{k})\in\mathcal{X}\subset R^{p}$ are independently drawn from a distribution $\mathcal{U}(\bm{{x}})$. Let $p_{\theta}(\tau)$ be a measure induced by the process $\tau(\cdot)$ on space $\mathcal{F}=\{\tau(\cdot):\mathcal{X}\rightarrow R\}.$

Denote

[TABLE]

[TABLE]

where $p_{TP}(\bm{{Y}}_{ij})$ is the Bayesian predictive distribution of $\bm{{Y}}_{ij}$ based on the t-process functional ANOVA model, and $p_{0}(\bm{{Y}}_{ij})$ is based on the true underlying function ${\tau}_{0}$ of $\tau_{ij}$. It is said that $p_{TP}(\bm{{Y}}_{ij})$ achieves information consistency if

[TABLE]

where $E_{\bm{{U}}_{n}}$ denotes the expectation under the distribution of $\bm{{U}}_{n}$ and $D[f,g]$ is the Kullback-Leibler divergence between functions $f$ and $g$, that is,

[TABLE]

Theorem 2. Under the appropriate conditions in Lemma 1 and condition (A) in Appendix B, we have for $j=1,\ldots,J_{i},~{}i=1,...,I$

[TABLE]

*where the expectation is taken over the distribution $\bm{{U}}_{n}$. *

The proof of Theorem 2 is in Appendix B.

4 Numerical studies

4.1 Simulation studies

Numerical simulations are constructed to evaluate the robustness of the t-process functional ANOVA model (2.2) (TP). The performance of TP is investigated by comparison of model (1.1) with error term having Gaussian process (GP). In addition, we also compute model (2.2) ignoring the random effect $\tau_{ij}$, denoted by TP($\tau=0$).

Simulation data are generated from three different models:

Model 1: The ANOVA model (1.1) with $I=3$ and $J=2$, where $\mu(t)=t^{2}$, $\alpha_{1}(t)=0.5\cos(t)$, $\alpha_{2}(t)=0.1\cos(t)$ and $\alpha_{3}(t)=-0.6\cos(t)$, $u_{ij}(t)=0.2t$ and $\theta=(0.1,10,0.1)$, $\varepsilon_{ij}\sim ETP(1.1,0.1,0,0.1)$.

• 2.

Model 2: $\varepsilon_{ij}\sim GP(0,0.1)$, the other setups are the same as Model 1.

• 3.

Model 3: $\tau_{ij}=0$, the other setups are the same as Model 1.

Model 1 is the true model for the TP method, while Model 2 is the true model for the GP method. Model 3 does not consider the random effect. Let observed time $t$ takes 61 points evenly spaced in [0,2]. Take $n$ sample points as training data with $n=11$ and $21$, the left are test data. Furthermore, to illustrate the performance of robustness, we randomly choose one point from the training data and disturb it by adding an extra error: a constant error 2 or random errors: $N(0,2)$ and $t(3)$, where $N(0,2)$ is normal distribution with mean 0 and variance 2, and $t(3)$ is student $t$ distribution with degree of freedom 3. All simulation are repeated 500 times.

Prediction errors (PE) and mean squared errors (MSE) of prediction from the TP, GP and TP($\tau=0$) methods are presented in Tables 1 - 3 for Model 1 - 3, respectively. From these tables, it follows that the TP method has the smallest PE and MSE, while the GP has comparable results with the TP($\tau=0$). It shows that the TP has more robustness against outliers than the GP, and better prediction than the TP($\tau=0$).

4.2 Real data examples

The robust functional ANOVA model is applied to two datasets, Canada weather data and market penetration of new products. Canadian weather data includes temperature and precipitation of weather at 35 Canadian weather observation stations, which is available in the R-package fda (Ramsay et al., 2011). These stations are divided into 4 regions: Arctic, Atlantic, Continental and Pacific. We study the week-average precipitation for 4 regions, where each region is a level and the week-average temperature is variable for random effect $\tau$. For the market penetration data (Ashish et al., 2009), we take the market penetration of Hi-Fi stereo and freezer as examples to illustrate our method. There are 86 countries in 7 continents from 1977 to 2015. Due to the missing data, we finally take data of 36 countries from 1983 to 2012 for the Hi-Fi stereo, and 35 countries from 1986 to 2015 for the freezer. With respect to relation with the Hi-Fi stereo and freezer, for random effect $\tau$, we take the penetration of personal computer as variable for the Hi-Fi stereo, and that of microwave oven for the freezer. Because 36 countries come from seven regions: Asia Pacific, Australasia, Eastern Europe, Latin America, Middle East and Africa, North America and Western Europe. So, the number of levels is 7.

Figures 1, 2 and 3 present the main effects and examples of the personal effect of some subjects for the weather data, the market penetration of Hi-Fi stereo and freezer, respectively. We can see that besides different main effects for the studied levels, each subject (weather station or country) also has different personal character.

For these two datasets, we respectively randomly choose 30 and 20 time points as training data and the rest as test data. Three estimation methods, TP, GP and TP($\tau=0$), are applied to fit the training data and then predict the test data. These procedures are repeated for 500 times. Table 4 shows prediction errors from TP, GP and TP($\tau=0$). We can see the TP has the smallest prediction errors, while GP and TP($\tau=0$) have the comparable ones.

5 Conclusion

The functional ANOVA model with Gaussian process offers a model for data with multi-dimensional covariances, and the specification of covariance kernels gives the ability to analyze a wide class of different type data to the model. However, the model with Gaussian process is not robust to outliers. Some heavy-tail processes can be used to overcome this shortcoming. This paper uses Gaussian process and t process to fit the personal random effect and error term, respectively, which results in a robust functional ANOVA model with t-process. Robustness property and information consistency of the proposed method are obtained. Numerical examples show that the robust functional ANOVA model has more robustness than that with Gaussian process.

Furthermore, it is a natural way to use t process to fit the random effect instead of Gaussian process, which leads to both of the personal random effect and error term having t processes. This extension makes estimation procedure more complicated because of more intractable integration.

Appendix A Robustness

**Proof of Theorem 1:

**From the log-likelihood function $l(\bm{{\theta}},\bm{{B}},\nu)$ defined in (2.7), we get the score functions of $\text{Vec}(\bm{{B}}^{\top})$,$~{}\bm{{\theta}}_{ij}$,$~{}\sigma^{2}$ as follows,

[TABLE]

[TABLE]

where for $q=1,\ldots,p,p=1,\ldots,L\times(I+1),k=1,\ldots,n$,

[TABLE]

[TABLE]

So, we can get that when $\bm{{Y}}_{ij}\rightarrow\infty$ for $j$-th object under level $i$, the score functions of $\text{Vec}(\bm{{B}}^{\top})$, $~{}\bm{{\theta}}_{ij}$ and $~{}\sigma^{2}$ are bounded, while that from the GPR model tend to $\infty$.

Let $T(F_{n})=T_{n}(\bm{{Y}}_{11},\ldots,\bm{{Y}}_{IJ_{I}})$ be an estimate of $\bm{{\beta}}$ (in this subsection let $\bm{{\beta}}=\bm{{\theta}},\bm{{B}}$), in which $F_{n}$ is the empirical distribution of $\{\bm{{Y}}_{11},\ldots,\bm{{Y}}_{IJ_{I}}\}$ and $T$ is a functional on some subset of all distributions. Let $\bm{{Y}}=(\bm{{Y}}_{11}^{\top},\ldots,$ $\bm{{Y}}_{IJ_{I}}^{\top})$, according to Hampel et al. (1986), influence function of $T$ at $F$ is defined as

[TABLE]

where $\delta_{\bm{{y}}}$ put mass 1 on point $\bm{{y}}$ and 0 on others.

Following Hampel et al. (1986) , we can get that estimator $\hat{\bm{{\beta}}}$ of $\bm{{\beta}}$ has the influence function

[TABLE]

The matrix $\frac{\partial^{2}l(\bm{{\theta}},\bm{{B}},\nu)}{\partial\bm{{\beta}}\partial\bm{{\beta}}^{\top}}$ is bounded according to $\bm{{Y}}_{ij},i=1,\ldots,I,~{}j=1,\ldots,J_{i}$, therefore the influence functions of $\hat{\bm{{\beta}}}$ are bounded under the propsoed functional ANOVA model. However, those under the GPR model are unbounded. $\sharp$

Appendix B Information consistency

Lemma 1**.**

Suppose $\bm{{Y}}_{ij}$ are independent samples from model (2.2) and $\tau_{0}\in\mathcal{F}$ has a Gaussian process with zero mean and bounded kernel function $k(\cdot,\cdot;\bm{{\theta}}_{ij})$. Let $k(\cdot,\cdot;\bm{{\theta}}_{ij})$ be continuous in $\bm{{\theta}}_{ij}$ and $\hat{\bm{{\theta}}}_{ijn}$ be one consistent estimator of $\bm{{\theta}}_{ij}$. Then we have

[TABLE]

where $\|\tau_{0}\|_{k}$ is the reproducing kernel Hilbert space(RKHS) norm of $\tau_{0}$ associated with $k(\cdot,\cdot;\bm{{\theta}}_{ij})$, $\bm{{K}}_{ij}$ is covariance matrix of $\tau_{0}$ over $\bm{{U}}_{n}$, $\bm{{I_{n}}}$ is the $n\times n$ identity matrix and $c,\delta$ are some positive constants.

Proof:  Let $\mathcal{H}$ is the reproducing kernel Hilbert space(RKHS) associated with $k(\cdot,\cdot;\bm{{\theta}}_{ij})$, and $\mathcal{H}_{n}$ is the span of $\{k(\cdot,\bm{{u}}_{ij}(t_{k});\bm{{\theta}}_{ij}\}$, that is, $\mathcal{H}_{n}=\{f(\cdot):f(\bm{{x}})=\sum\nolimits_{k=1}^{n}{\alpha_{k}k(\bm{{x}},\bm{{u}}_{ij}(t_{k});\bm{{\theta}}_{ij})},\text{for any}~{}\alpha_{k}\in R\}$ . First, we assume that $\tau_{0}\in\mathcal{H}_{n}$, then we can express $\tau_{0}(\cdot)$ as

[TABLE]

where $K(\cdot)=(k(\cdot,\bm{{u}}_{ij}(t_{1});\bm{{\theta}}_{ij}),\ldots,k(\cdot,\bm{{u}}_{ij}(t_{n});\bm{{\theta}}_{ij}))$ and $\bm{{\alpha}}=(\alpha_{1},\ldots,\alpha_{n})^{\top}$. With the property of RKHS, we can get $\|\tau_{0}\|_{k}^{2}=\bm{{\alpha}}^{\top}\bm{{K}}_{ij}\bm{{\alpha}}$, and $(\tau_{0}(\bm{{u}}_{ij}(t_{1})),$ $\ldots,\tau_{0}(\bm{{u}}_{ij}(t_{n})))^{\top}=\bm{{K}}_{ij}\bm{{\alpha}}$.

By the Fenchel-Legendre duality relationship, we can get

[TABLE]

where $p$ is the measure induced by $GP(0,k(\cdot,\cdot;\hat{\bm{{\theta}}}_{ij}))$, therefore its finite dimensional distribution at $\{y_{ij}(t_{1}),\ldots,y_{ij}(t_{n})\}$ is $N(0,\hat{\bm{{K}}}_{ij})$, where $\hat{\bm{{K}}}_{ij}$ is defined in the same way as $\bm{{K}}_{ij}$ but with $\bm{{\theta}}_{ij}$ being replaced by its estimator $\hat{\bm{{\theta}}}_{ijn}$, and $\overline{p}$ is the posterior distribution of ${\tau}_{ij}$ with a prior distribution $GP(0,k(\cdot,\cdot;\bm{{\theta}}_{ij}))$ and normal likelihood $\prod\nolimits_{k=1}^{n}{N(\hat{y}_{ij}(t_{k});{\tau}_{ij},\sigma^{2})}$, where $\hat{\bm{{y}}}_{ij}\triangleq(\hat{y}_{ij}(t_{1}),\ldots,\hat{y}_{ij}(t_{n}))^{\top}=(\bm{{K}}_{ij}+\sigma^{2}\bm{{I_{n}}})\bm{{\alpha}}$.

Therefore, by the Gaussian process regression, the posterior distribution of ${\tau}_{ij}$ is $\overline{p}({\tau}_{ij})=N(\bm{{K}}_{ij}\bm{{\alpha}},\bm{{K}}_{ij}\bm{{B}}^{-1})$, where $\bm{{B}}=\bm{{I_{n}}}+\sigma^{-2}\bm{{K}}_{ij}$. Then we have

[TABLE]

In addition, we show that

[TABLE]

and

[TABLE]

Hence, it follow from (B.1), (B.2), (B.3) and (B) that

[TABLE]

Since the covariance function is bounded and the covariance function is continous in $\bm{{\theta}}_{ij}$ and $\hat{\bm{{\theta}}}_{ijn}\rightarrow\bm{{\theta}}_{ij}$, we have $\hat{\bm{{K}}}_{ij}^{-1}\bm{{K}}_{ij}-\bm{{I_{n}}}\rightarrow 0$ as $n\rightarrow\infty$. Hence, there exist positive constants $c$ and $\varepsilon$ such that for $n$ large enough, $-\text{log}|\hat{\bm{{K}}}_{ij}^{-1}\bm{{K}}_{ij}|<c,$ $\bm{{\alpha}}^{\top}\bm{{K}}_{ij}(\hat{\bm{{K}}}_{ij}^{-1}\bm{{K}}_{ij}-\bm{{I_{n}}})\bm{{\alpha}}<c,$ and $\text{tr}(\hat{\bm{{K}}}_{ij}^{-1}\bm{{K}}_{ij}\bm{{B}}^{-1})<\text{tr}((\bm{{I_{n}}}+\varepsilon\bm{{K}}_{ij})\bm{{B}}^{-1}),$ which plugs in (B), there exist positive constant $\delta$, we have the inequality

[TABLE]

From the Representer Theorem, it can obtain

[TABLE]

for all $\tau_{0}(\cdot)\in\mathcal{H}$. The proof is over. $\sharp$

To prove Theorem 2, it need condition (A) $\|\tau_{0}\|_{k}$ is bounded and $E_{\bm{{U}}_{n}}(\text{log}|\bm{{I_{n}}}+\sigma^{-2}\bm{{K}}_{ij}|)=o(n)$.

**Proof of Theorem 2.

**From the definition of the information consistency, we get that

[TABLE]

From Lemma 1, we obtain that

[TABLE]

where $c$ and $\delta$ are two positive constants. Because $\|\tau_{0}\|_{k}$ is bounded and expected regret term $E_{\bm{{U}}_{n}}(\text{log}|\bm{{I_{n}}}+\sigma^{-2}\bm{{K}}_{ij}|)=o(n)$, Theorem 2 holds. $\sharp$

Reference