The effect of the spatial domain in FANOVA models with ARH(1) error term

TL;DR

This paper investigates how the spatial domain shape influences FANOVA models with ARH(1) errors, introducing a new statistical test and analyzing boundary effects on dependence and estimation stability.

Contribution

It extends FANOVA models with ARH(1) errors to rectangular and circular supports under Dirichlet boundary conditions, and develops a new significance test for fixed effects.

Findings

Boundary conditions influence error dependence range.

Eigenvalue decay affects estimation stability.

Simulation and fMRI application demonstrate method effectiveness.

Abstract

Functional Analysis of Variance (FANOVA) from Hilbert-valued correlated data with spatial rectangular or circular supports is analyzed, when Dirichlet conditions are assumed on the boundary. Specifically, a Hilbert-valued fixed effect model with error term defined from an Autoregressive Hilbertian process of order one (ARH(1) process) is considered, extending the formulation given in Ruiz-Medina (2016). A new statistical test is also derived to contrast the significance of the functional fixed effect parameters. The Dirichlet conditions established at the boundary affect the dependence range of the correlated error term. While the rate of convergence to zero of the eigenvalues of the covariance kernels, characterizing the Gaussian functional error components, directly affects the stability of the generalized least-squares parameter estimation problem. A simulation study and a real-data application related to fMRI analysis are undertaken to illustrate the performance of the parameter estimator and statistical test derived.

Tables

Table 1. Table 1 : Scenarios for rectangular domain.

Cases	$a_1=a_2$	$b_1=b_2$	$h_x=h_y$	$p$	$TR$
(P1,a,C1)	$-2$	$3$	$0.05$	$4$	$4\times 4$
(P1,b,C2)	$-2$	$3$	$0.05$	$4$	$6\times 6$
(P1,c,C2)	$-2$	$3$	$0.05$	$4$	$8\times 8$
(P1,d,C1)	$-2$	$3$	$0.05$	$4$	$12\times 12$
(P2,a,C2)	$-2$	$3$	$0.05$	$9$	$4\times 4$
(P2,b,C1)	$-2$	$3$	$0.05$	$9$	$6\times 6$
(P2,c,C1)	$-2$	$3$	$0.05$	$9$	$8\times 8$
(P2,d,C2)	$-2$	$3$	$0.05$	$9$	$12\times 12$

Table 2. Table 4 : F 𝐹 F statistics ( 25 ) for rectangular domain.

Cases	(P1,a,C1)	(P1,b,C2)	(P1,c,C2)	(P1,d,C1)	(P2,a,$C_2$)	($P_2$,b,$C_1$)	(P2,c,C1)	(P2,d,C2)
$F$	$1.926$	$1.717$	$1.673$	$1.626$	$1.898$	$1.845$	$1.761$	$1.606$

Table 3. Table 5 : Rectangle . Percentage of successes for α = 0.05 , 𝛼 0.05 \alpha=0.05, at the left–hand side, and averaged p 𝑝 p –values at the right–hand side, for each one of the eight realizations considered of the Gaussian function h ∈ L 2 ​ ( [ − 2 , 3 ] × [ − 2 , 3 ] ) . ℎ superscript 𝐿 2 2 3 2 3 h\in L^{2}([-2,3]\times[-2,3]).

$D$	$\%$ Success	$p$
1	$100\%$	$0$
2	$100\%$	$0$
3	$99.75\%$	$1.998{\left(10\right)}^{-8}$
4	$100\%$	$0$
5	$99.8\%$	$7.541{\left(10\right)}^{-7}$
6	$100\%$	$0$
7	$100\%$	$0$
8	$100\%$	$6.441{\left(10\right)}^{-10}$

Table 4. Table 6 : Scenarios for disk domain.

Cases	$R$	$h_R$	$h_{\varphi }$	$TR$	$p$
(P1,a,C3)	$12$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$3$	$4$
(P1,b,C2)	$18$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$5$	$4$
(P1,c,C1)	$25$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$7$	$4$
(P1,d,C1)	$50$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$15$	$4$
(P1,e,C2)	$100$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$31$	$4$
(P1,f,C3)	$250$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$79$	$4$
(P2,a,C1)	$12$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$3$	$9$
(P2,b,C2)	$18$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$5$	$9$
(P2,c,C3)	$25$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$7$	$9$
(P2,d,C3)	$50$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$15$	$9$
(P2,e,C2)	$100$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$31$	$9$
(P2,f,C1)	$250$	$\frac{R}{145}$	$\frac{2\pi }{135}$	$79$	$9$

Table 5. Table 9 : F 𝐹 F statistics ( 25 ) over the disk domain.

Cases	(P1,a,C3)	(P1,b,C2)	(P1,c,C1)
$F$	$1.100\left({10}^2\right)$	$4.100\left({10}^3\right)$	$1.200\left({10}^5\right)$
Cases	(P1,d,C1)	(P1,e,C2)	(P1,f,C3)
$F$	$3.900\left({10}^6\right)$	$6.300\left({10}^6\right)$	$4.200\left({10}^6\right)$
Cases	(P2,a,C1)	(P2,b,C2)	(P2,c,C3)
$F$	$2.200\left({10}^3\right)$	$8.200\left({10}^3\right)$	$7.600\left({10}^7\right)$
Cases	(P2,d,C3)	(P2,e,C2)	(P2,f,C1)
$F$	$2.500\left({10}^7\right)$	$1.400\left({10}^7\right)$	$8.500\left({10}^7\right)$

Table 6. Table 10 : Disk . Percentage of successes for α = 0.05 , 𝛼 0.05 \alpha=0.05, at the left–hand side, and averaged p 𝑝 p –values at the right–hand side, for each one of the eight realizations of the Gaussian function h ∈ L 2 ​ ( D 25 ) . ℎ superscript 𝐿 2 subscript 𝐷 25 h\in L^{2}(D_{25}).

$D$	% Success	$p$
1	$99.95\%$	$1.672{\left(10\right)}^{-8}$
2	$99.5\%$	$9.746{\left(10\right)}^{-7}$
3	$100\%$	$0$
4	$99.9\%$	$8.546{\left(10\right)}^{-8}$
5	$97.45\%$	$7.400{\left(10\right)}^{-7}$
6	$100\%$	$0$
7	$100\%$	$8.775{\left(10\right)}^{-9}$
8	$100\%$	$0$

Table 7. Table 11 : Scenarios for circular sector domain.

Cases	$R$	$h_R$	$h_{\varphi }$	$TR$	$\theta$	$p$
(P1,a,C3)	$12$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$3$	$\frac{2}{3}$	$4$
(P1,b,C2)	$18$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$5$	$\frac{2}{3}$	$4$
(P1,c,C1)	$25$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$7$	$\frac{2}{3}$	$4$
(P1,d,C1)	$50$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$15$	$\frac{2}{3}$	$4$
(P1,e,C2)	$100$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$31$	$\frac{2}{3}$	$4$
(P1,f,C3)	$250$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$79$	$\frac{2}{3}$	$4$
(P2,a,C1)	$12$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$3$	$\frac{2}{3}$	$9$
(P2,b,C2)	$18$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$5$	$\frac{2}{3}$	$9$
(P2,c,C3)	$25$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$7$	$\frac{2}{3}$	$9$
(P2,d,C3)	$50$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$15$	$\frac{2}{3}$	$9$
(P2,e,C2)	$100$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$31$	$\frac{2}{3}$	$9$
(P2,f,C1)	$250$	$\frac{R}{145}$	$\frac{2\pi }{115}$	$79$	$\frac{2}{3}$	$9$

Table 8. Table 14 : F 𝐹 F statistics ( 25 ) for the circular sector.

Cases	(P1,a,C3)	(P1,b,C2)	(P1,c,C1)	(P1,d,C1)	(P1,e,C2)	(P1,f,C3)
$F$	$9.2\left({10}^2\right)$	$3.1\left({10}^3\right)$	$4.2\left({10}^6\right)$	$4.8\left({10}^8\right)$	$5.8\left({10}^6\right))$	$7.3\left({10}^8\right)$
Cases	(P2,a,C1)	(P2,b,C2)	(P2,c,C3)	(P2,d,C3)	(P2,e,C2)	(P2,f,C1)
$F$	$1.8\left({10}^3\right)$	$4.1\left({10}^3\right)$	$2.6\left({10}^7\right)$	$3.1\left({10}^9\right)$	$6.8\left({10}^6\right)$	$1.8\left({10}^9\right)$

Table 9. Table 15 : Circular Sector. Percentage of successes for α = 0.05 , 𝛼 0.05 \alpha=0.05, at the left–hand side, and averaged p 𝑝 p –values at the right–hand side, for each one of the eight realizations of the Gaussian function h ∈ L 2 ​ ( C ​ S ) . ℎ superscript 𝐿 2 𝐶 𝑆 h\in L^{2}(CS).

$D$	% Success	$p$
1	$97.5\%$	$6.504{\left(10\right)}^{-6}$
2	$100\%$	$0$
3	$100\%$	$3.600{\left(10\right)}^{-8}$
4	$100\%$	$0$
5	$98\%$	$2.006{\left(10\right)}^{-6}$
6	$99.5\%$	$9.807{\left(10\right)}^{-8}$
7	$100\%$	$0$
8	$99.5\%$	$4.111{\left(10\right)}^{-7}$

Table 10. Table 16 : E ​ F ​ M ​ S ​ E 𝒀 i f ​ M ​ R ​ I 𝐸 𝐹 𝑀 𝑆 subscript 𝐸 superscript subscript 𝒀 𝑖 𝑓 𝑀 𝑅 𝐼 EFMSE_{\boldsymbol{Y}_{i}^{fMRI}} and E ​ F ​ M ​ S ​ E 𝒀 i H 𝐸 𝐹 𝑀 𝑆 subscript 𝐸 superscript subscript 𝒀 𝑖 𝐻 EFMSE_{\boldsymbol{Y}_{i}^{H}} for case A.

Slices $S_i$	$EFMS{E}_{{𝒀}_i^{fMRI}}$	$EFMS{E}_{{𝒀}_i^H}$
1	$2.417{\left(10\right)}^{-}3$	$3.492{\left(10\right)}^{-}3$
2	$3.051{\left(10\right)}^{-}3$	$3.119{\left(10\right)}^{-}3$
3	$4.293{\left(10\right)}^{-}3$	$5.523{\left(10\right)}^{-}3$
4	$6.666{\left(10\right)}^{-}3$	$7.690{\left(10\right)}^{-}3$
5	$8.986{\left(10\right)}^{-}3$	$9.961{\left(10\right)}^{-}3$
6	$8.462{\left(10\right)}^{-}3$	$9.434{\left(10\right)}^{-}3$
7	$1.108{\left(10\right)}^{-}2$	$1.920{\left(10\right)}^{-}2$
8	$1.720{\left(10\right)}^{-}2$	$2.720{\left(10\right)}^{-}2$
9	$1.499{\left(10\right)}^{-}2$	$1.914{\left(10\right)}^{-}2$
10	$1.036{\left(10\right)}^{-}2$	$1.851{\left(10\right)}^{-}2$
11	$1.308{\left(10\right)}^{-}2$	$1.634{\left(10\right)}^{-}2$
12	$1.302{\left(10\right)}^{-}2$	$1.300{\left(10\right)}^{-}2$
13	$7.850{\left(10\right)}^{-}3$	$7.939{\left(10\right)}^{-}3$
14	$6.640{\left(10\right)}^{-}3$	$6.730{\left(10\right)}^{-}3$
15	$3.511{\left(10\right)}^{-}3$	$2.832{\left(10\right)}^{-}3$
16	$2.771{\left(10\right)}^{-}3$	$3.540{\left(10\right)}^{-}3$

Table 11. Table 17 : E ​ F ​ M ​ S ​ E 𝒀 i f ​ M ​ R ​ I 𝐸 𝐹 𝑀 𝑆 subscript 𝐸 superscript subscript 𝒀 𝑖 𝑓 𝑀 𝑅 𝐼 EFMSE_{\boldsymbol{Y}_{i}^{fMRI}} and E ​ F ​ M ​ S ​ E 𝒀 i H 𝐸 𝐹 𝑀 𝑆 subscript 𝐸 superscript subscript 𝒀 𝑖 𝐻 EFMSE_{\boldsymbol{Y}_{i}^{H}} for case B.

Slices $S_i$	$EFMS{E}_{{𝒀}_i^{fMRI}}$	$EFMS{E}_{{𝒀}_i^H}$
1	$2.417{\left(10\right)}^{-}3$	$2.592{\left(10\right)}^{-}3$
2	$3.051{\left(10\right)}^{-}3$	$3.119{\left(10\right)}^{-}3$
3	$4.293{\left(10\right)}^{-}3$	$4.733{\left(10\right)}^{-}3$
4	$6.666{\left(10\right)}^{-}3$	$7.671{\left(10\right)}^{-}3$
5	$8.986{\left(10\right)}^{-}3$	$9.065{\left(10\right)}^{-}3$
6	$8.462{\left(10\right)}^{-}3$	$8.435{\left(10\right)}^{-}3$
7	$1.108{\left(10\right)}^{-}2$	$1.120{\left(10\right)}^{-}2$
8	$1.720{\left(10\right)}^{-}2$	$1.919{\left(10\right)}^{-}2$
9	$1.499{\left(10\right)}^{-}2$	$1.524{\left(10\right)}^{-}2$
10	$1.036{\left(10\right)}^{-}2$	$1.040{\left(10\right)}^{-}2$
11	$1.308{\left(10\right)}^{-}2$	$1.481{\left(10\right)}^{-}2$
12	$1.302{\left(10\right)}^{-}2$	$1.299{\left(10\right)}^{-}2$
13	$7.849{\left(10\right)}^{-}3$	$7.929{\left(10\right)}^{-}3$
14	$6.640{\left(10\right)}^{-}3$	$6.719{\left(10\right)}^{-}3$
15	$3.511{\left(10\right)}^{-}3$	$2.829{\left(10\right)}^{-}3$
16	$2.771{\left(10\right)}^{-}3$	$3.540{\left(10\right)}^{-}3$

Table 12. Table 18 : Percentage of brain voxels per slice, where the real-valued fixed effect model with AR(1) error term, fitted by fmrilm.m MatLab function, is significative.

$S$	% voxels with rejection of $H_0$
$1$	$99.927\%$
$2$	$99.927\%$
$3$	$99.707\%$
$4$	$99.902\%$
$5$	$99.805\%$
$6$	$99.951\%$
$7$	$99.927\%$
$8$	$99.976\%$
$9$	$99.805\%$
$10$	$99.951\%$
$11$	$99.951\%$
$12$	$99.902\%$
$13$	$99.878\%$
$14$	$99.951\%$
$15$	$99.951\%$
$16$	$100\%$

Table 13. Table 19 : p 𝑝 p -values for T 𝒉 subscript 𝑇 𝒉 T_{\boldsymbol{h}} computed at the 16 16 16 slices, considering four random directions, for T ​ R = 16 𝑇 𝑅 16 TR=16 .

$S$	$D_1$	$D_2$	$D_3$	$D_4$
$1$	$0$	$0$	$0.082$	$0.023$
$2$	$0.590{\left(10\right)}^{-2}$	$0$	$0$	$0$
$3$	$0.018$	$0.066$	$0.049$	$0.030$
$4$	$0$	$0$	$0$	$0.170{\left(10\right)}^{-10}$
$5$	$0$	$0.026$	$0$	$0$
$6$	$0$	$0$	$0$	$0$
$7$	$0.710{\left(10\right)}^{-7}$	$0$	$0$	$0$
$8$	$0$	$0.006$	$0$	$0$
$9$	$0.049$	$0$	$0$	$0.023$
$10$	$0.390{\left(10\right)}^{-7}$	$0.031$	$0$	$0$
$11$	$0.004$	$0.006$	$0.660{\left(10\right)}^{-6}$	$0.052$
$12$	$0.046$	$0$	$0$	$0.034$
$13$	$0.340{\left(10\right)}^{-7}$	$0.028$	$0$	$0.440{\left(10\right)}^{-3}$
$14$	$0$	$0.180{\left(10\right)}^{-6}$	$0.021$	$0.050$
$15$	$0$	$0.140{\left(10\right)}^{-7}$	$0.044$	$0.052$
$16$	$0.110{\left(10\right)}^{-4}$	$0.230{\left(10\right)}^{-7}$	$0$	$0$

Table 14. Table 20 : p 𝑝 p -values for T 𝒉 subscript 𝑇 𝒉 T_{\boldsymbol{h}} computed at the 16 16 16 slices, considering four random directions, for T ​ R = 4 𝑇 𝑅 4 TR=4 .

$S$	$D_1$	$D_2$	$D_3$	$D_4$
$1$	$0$	$0.051$	$0.071$	$0.011$
$2$	$0.880{\left(10\right)}^{-4}$	$0$	$0$	$0$
$3$	$0.067$	$0.034$	$0$	$0.037$
$4$	$0$	$0.250{\left(10\right)}^{-4}$	$0.110{\left(10\right)}^{-4}$	$0.016$
$5$	$0.370{\left(10\right)}^{-6}$	$0$	$0.280{\left(10\right)}^{-6}$	$0$
$6$	$0.001$	$0$	$0$	$0.220{\left(10\right)}^{-4}$
$7$	$0.064$	$0.034$	$0.007$	$0.044$
$8$	$0.072$	$0.079$	$0.035$	$0$
$9$	$0.220{\left(10\right)}^{-5}$	$0.470{\left(10\right)}^{-4}$	$0.004$	$0.220{\left(10\right)}^{-9}$
$10$	$0$	$0.120{\left(10\right)}^{-3}$	$0.370{\left(10\right)}^{-4}$	$0.970{\left(10\right)}^{-7}$
$11$	$0.081$	$0.058$	$0$	$0$
$12$	$0.870{\left(10\right)}^{-4}$	$0$	$0$	$0.036$
$13$	$0.760{\left(10\right)}^{-3}$	$0$	$0$	$0.370{\left(10\right)}^{-3}$
$14$	$0.210{\left(10\right)}^{-6}$	$0$	$0$	$0.037$
$15$	$0$	$0.650{\left(10\right)}^{-4}$	$0.032$	$0$
$16$	$0.540{\left(10\right)}^{-6}$	$0$	$0$	$0.520{\left(10\right)}^{-3}$

Equations

$$Y_i(\mathbf{s})={\mathbf{x}}_i\mathbf{\beta }(\mathbf{s})+Z_i(\mathbf{s}){\sigma }_i(\mathbf{s}),\mathbf{s}\in S,$$

$$\mathbf{Y}\left(\cdot \right)=\mathbf{X}\mathbf{\beta }\left(\cdot \right)+\mathbf{\epsilon }\left(\cdot \right),$$

$$\mathbf{\beta }(\cdot )=[{\beta }_1(\cdot ),\dots ,{\beta }_p(\cdot ){]}^T\in H^p$$

$$\mathbf{Y}(\cdot )=[Y_1(\cdot ),\dots ,Y_n(\cdot ){]}^T$$

$$\mathbf{\epsilon }(\cdot )=[{\epsilon }_1(\cdot ),\dots ,{\epsilon }_n(\cdot ){]}^T$$

$${\epsilon }_m\left(\cdot \right)=\rho \left({\epsilon }_{m-1}\right)\left(\cdot \right)+{\nu }_m\left(\cdot \right),m\in \mathbb{Z},$$

$$\mathbf{\epsilon }(\cdot )=[{\epsilon }_1(\cdot ),\dots ,{\epsilon }_n(\cdot ){]}^T$$

$$\mathrm{E}\left\{{\epsilon }_i\otimes {\epsilon }_j\right\}=R_1,\text{if }j-i=1,\mathrm{E}\left\{{\epsilon }_i\otimes {\epsilon }_j\right\}=R_1^{\ast },\text{if }i-j=1,$$

$$\mathrm{E}\left\{{\epsilon }_i\otimes {\epsilon }_i\right\}=R_0,\text{if }i=j,$$

$${\mathbf{R}}_{\mathbf{\epsilon }\mathbf{\epsilon }}$$

$$\simeq$$

$$⟨\mathbf{f},\mathbf{g}{⟩}_{H^n}=\sum _{i=1}^n⟨f_i,g_i{⟩}_H,\text{leavevmode}\mathbf{f},\mathbf{g}\in H^n.$$

$$R_0=\sum _{k=1}^{\infty }{\lambda }_k{\varphi }_k\otimes {\varphi }_k,$$

$${\epsilon }_i=\sum _{k=1}^{\infty }{⟨{\epsilon }_i,{\varphi }_k⟩}_H{\varphi }_k=\sum _{k=1}^{\infty }\sqrt{{\lambda }_k}{\eta }_k(i){\varphi }_k,i=1,\dots ,n,$$

$$\sqrt{{\lambda }_k}{\eta }_k(i)={⟨{\epsilon }_i,{\varphi }_k⟩}_H,$$

$$\mathrm{E}\left\{{\eta }_k(i){\eta }_p(j)\right\}={\delta }_{k,p},k,p\in \mathbb{N},$$

$$R_1$$

$${\mathbf{\Phi }}^{\ast }\left(\mathbf{f}\right)$$

$$\mathbf{\Phi }\left(\left\{{\mathbf{f}}_k^T,k\ge 1\right\}\right)={\left(\sum _{k=1}^{\infty }{f}_{k1}{\varphi }_k,\dots ,\sum _{k=1}^{\infty }{f}_{kn}{\varphi }_k\right)}^T.$$

$$A_{i,j}=\sum _{k=1}^{\infty }{\gamma }_{kij}{\varphi }_k\otimes {\varphi }_k$$

$${\mathbf{\Phi }}^{\ast }\mathbf{A}\mathbf{\Phi }$$

$${\mathbf{\Phi }}^{\ast }{\mathbf{R}}_{\mathbf{\epsilon }\mathbf{\epsilon }}\mathbf{\Phi }$$

$${\mathbf{R}}_{\mathbf{\epsilon }\mathbf{\epsilon }}^{-1}\left(\mathbf{f},\mathbf{g}\right)$$

$$\Vert \mathbf{f}{\Vert }_{{\mathbf{R}}_{\mathbf{\epsilon }\mathbf{\epsilon }}^{-1}}^2$$

$$\boldsymbol{\Lambda}_{k}=\left[\begin{array}[]{cccccc}\lambda_{k}(R_{0})&\lambda_{k}(R_{1})&0&\dots&0&0\\ \lambda_{k}(R_{1}^{*})&\lambda_{k}(R_{0})&\lambda_{k}(R_{1})&\dots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\dots&\lambda_{k}(R_{1}^{*})&\lambda_{k}(R_{0})\\ \end{array}\right],$$

$${\Lambda }_{kij}$$

$${\Lambda }_{kii}$$

$${\gamma }_i\in (0,d/2),i=1,\dots ,n,$$

$$\widehat{R_0}$$

$${\mathbf{X}}^T\mathbf{X}={\mathbf{Id}}_p,{\mathbf{Id}}_p\in {\mathbb{R}}^{p\times p}.$$

$$\sum _{k=1}^{\infty }\mathrm{Tr}{\left({\mathbf{X}}^T{\mathbf{\Lambda }}_k^{-1}\mathbf{X}\right)}^{-1}<\infty .$$

$$\mathcal{H}\left(Z\right)=R_Z^{1/2}\left(\mathcal{H}\right)$$

Full text

Published in Statistics and Its Interface, 10, pp. 607-628. DOI: doi.org/10.4310/SII.2017.v10.n4.a7

1 Department of Statistics and O. R., University of Granada, Spain.

E-mail: javialvaliebana@ugr.es

Key words: ARH(1) error term; boundary value problems; Cram´er-Wold Theorem; functional analysis of variance; linear functional tests; fMRI data.

1 Introduction

In the last few decades, functional data analysis techniques have grown significantly given the new technologies available, in particular, in the field of medicine (see, for instance, Sorensen et al. [2013]). High–dimensional data, which are functional in nature, are generated, for example, from measurements in time, over spatial grids or images with many pixels (e.g., data on electrical activity of the heart, data on electrical activity along the scalp, data reconstructed from medical imaging, expression profiles in genetics and genomics, monitoring of continuity activity through accelerometers, etc). Effective experimental design and modern functional statistics have led to recent advances in medical imaging, improving, in particular, the study of human brain function (see, for example, Delzell et al. [2012]). Magnetic Resonance Imaging (MRI) data have been analysed with different aims. For example, we refer to the studies related with cortical thickness (see Lerch and Evans [2005]), where magnetic resonance imaging data are analysed to detect the spatial locations of the surface of the brain, where the cortical thickness is correlated with an independent variable, such as age or gender (see also Shaw et al. [2006]). Cortical thickness is usually previously smoothed along the surface of the brain (see Chung et al. [2005]). Thus, it can be considered as a functional random variable with spatial circular support. In general, the following linear model is considered, for cortical thickness $Y_{i}(s)$ on subjects $i=1,\dots,n,$:

[TABLE]

where $\boldsymbol{x}_{i}$ is a vector of known $p$ regressors, and for each $\mathbf{s}\in S,$ with $S$ denoting the surface of the brain, parameter $\boldsymbol{\beta}(\mathbf{s})$ is an unknown $p$–vector of regression coefficients. The errors $\left\{Z_{1},\dots,Z_{n}\right\}$ are independent zero-mean Gaussian random fields. In Taylor and Worsley [2007], this model is also considered to detect how the regressors are related to the data at spatial location $\mathbf{s},$ by testing contrasts in $\boldsymbol{\beta}(\mathbf{s}),$ for $\mathbf{s}\in S.$ The approach presented in this paper allows the formulation of model (1) in a functional (Hilbert–valued) framework, incorporating possible correlations between subjects, due to genetic characteristics, breed, geographic location, etc.

The statistical analysis of functional magnetic resonance image (fMRI) data has also generated an important activity in research about brain activity, where the functional statistical approach implemented in this paper could lead to important spatio–temporal analysis improvements. It is well–known that fMRI techniques have been developed to address the unobserved effect of scanner noise in studies of the auditory cortex. A penalized likelihood approach to magnetic resonance image reconstruction is presented in Bulaevskaya and Oehlert [2007]. A new approach which incorporates the spatial information from neighbouring voxels, as well as temporal correlation within each voxel, which makes use of regional kriging is derived in Christensen and Yetkin [2005]. Conditional autoregressive and Markov random field modelling involves some restrictions in the characterization of spatially contiguous effect regions, and, in general, in the representation of the spatial dependence between spatially connected voxels (see, for example, Banerjee et al. [2004], Besag [1986]). Multiscale adaptive regression models assume spatial independence to construct a weighted likelihood parameter estimate. At each scale, the weights determine the amount of information that observations in a neighborhood voxel provides on the parameter vector to be estimated at a given voxel, under the assumption of independence between the conditional distributions of the responses at the neighborhood voxels, for each scale. The weights are sequentially computed through different scales, for adaptively update of the parameter estimates and test statistics (see, for example, Li et al. [2011]).

In Zhu et al. [2012], a multivariate varying coefficient model is considered for neuroimaging data, under a mixed effect approach, to reflect dependence within–curve and between–curve, in the case where coefficients are one–parameter functions, although extension to higher dimension is straightforward. The approach presented in this paper adopts a functional framework to analyse multivariate varying coefficient models in higher dimensions (two–dimensional design points), under the framework of multivariate fixed effect models in Hilbert spaces. Namely, the response is a multivariate functional random variable reflecting dependence within-surface (between voxels), and between-surface (between different times), with Hilbert–valued multivariate Gaussian distribution. Hence, the varying coefficients are estimated from the application of an extended version of generalized least–squares estimation methodology, in the multivariate Hilbert–valued context (see Ruiz-Medina [2016]), while, in Zhu et al. [2012], local linear regression is applied to estimate the coefficient functions. The dependence structure of the functional response is estimated here from the moment–based parameter estimation of the ARH(1) error term (see Bosq [2000]). In Zhu et al. [2012], local linear regression technique is employed to estimate the random effects, reflecting dependence structure in the varying coefficient mixed effect model. An extended formulation of the varying coefficient model considered in Zhu et al. [2012] is given in Zhu et al. [2014], combining a univariate measurement mixed effect model, a jumping surface model, and a functional component analysis model. In the approach presented in this paper, we have combined a nonparametric surface model with a multivariate functional principal component approach in the ARH(1) framework. Thus, a continuous spatial variation of the fMRI response is assumed, incorporating temporal and spatial correlations (across voxels), with an important dimension reduction in the estimation of the varying coefficient functions.

The above–referred advances in medicine are supported by the extensive literature on linear models in function spaces developed in parallel in the last few decades. We particularly refer to the functional linear regression context (see, for example, Cai and Hall [2006], Cardot et al. [2003], Cardot and Sarda [2011], Chiou et al. [2004], Crambes et al. [2009], Cuevas et al. [2002], Ferraty et al. [2013], Kokoszka et al. [2008], among others). See also Bosq [2000, 2007], Ruiz-Medina [2011, 2012], in the functional time series context, and Ferraty and Vieu [2006, 2011] in the functional nonparametric regression framework. Functional Analysis of Variance (FANOVA) techniques for high–dimensional data with a functional background have played a crucial role, within the functional linear model literature as well. Related work has been steadily growing (see, for example, Angelini et al. [2003], Dette and Derbort [2001], Gu [2002], Huang [1998], Kaufman and Sain [2010], Kaziska [2011], Lin [2000], Ramsay and Silverman [2005], Spitzner et al. [2003], Stone et al. [1997], Wahba et al. [1995]). The paper Ruiz-Medina [2016] extends the results in Zoglat [2008] from the $L^{2}([0,1])$–valued context to the separable Hilbert–valued space framework, and from the case of independent homocedastic error components to the correlated heteroscedastic case. In the context of hypothesis testing from functional data, tests of significance based on wavelet thresholding are formulated in Fan [1996], exploiting the sparsity of the signal representation in the wavelet domain, for dimension reduction. A maximum likelihood ratio based test is suggested for functional variance components in mixed–effect FANOVA models in Guo [2002]. From classical ANOVA tests, an asymptotic approach is derived in Cuevas et al. [2004], for studying the equality of the functional means from $k$ independent samples of functional data. The testing problem for mixed–effect functional analysis of variance models is addressed in Abramovich and Angelini [2006], Abramovich et al. [2004], developing asymptotically optimal (minimax) testing procedures for the significance of functional global trend, and the functional fixed effects. The wavelet transform of the data is again used in the implementation of this approach (see also Antoniadis and Sapatinas [2007]). Recently, in the context of functional data defined by curves, considering the $L^{2}$–norm, an up–to–date overview of hypothesis testing methods for functional data analysis is provided in Zhang [2013], including functional ANOVA, functional linear models with functional responses, heteroscedastic ANOVA for functional data, and hypothesis tests for the equality of covariance functions, among other related topics.

In this paper, the model formulated in Ruiz-Medina [2016] is extended to the case where the error term is an ARH(1) process. Furthermore, an alternative test to contrast the significance of the functional fixed effect parameters is formulated, based on a sharp form of the Cramér–Wold’s Theorem derived in Cuesta-Albertos et al. [2007], for Gaussian measures on a separable Hilbert space. The simulation study undertaken illustrates the effect of the boundary conditions and the geometry of the domain on the spatial dependence range of the functional vector error term. Specifically, in that simulations, we consider the case where the Gaussian error components satisfy a stochastic partial differential equation, given in terms of a fractional power of the Dirichlet negative Laplacian operator. The autocovariance and cross–covariance operators of the functional error components are then defined in terms of the eigenvectors of the Dirichlet negative Laplacian operator. The eigenvectors of the Dirichlet negative Laplacian operator vanish continuously at the boundary, in the case of the regular domains studied (the rectangle, disk and circular sector), with decay velocity determined by the boundary conditions and the geometry of the domain. Thus, the boundary conditions and the geometry of the domain directly affect the dependence range of the error components, determined by the rate of convergence to zero of the Dirichlet negative Laplacian eigenvectors at the boundary. The influence of the truncation order is studied as well, since the rate of convergence to zero of the eigenvalues of the spatial covariance kernels, that define the matrix covariance operator of the error term, could affect the stability of the generalized least–squares estimation problem addressed here. Furthermore, in the fMRI data problem considered, the presented functional fixed effect model, with ARH(1) error term, is fitted. In that case, the temporal dependence range of the error term is controlled by the ARH(1) dynamics, while the spatial dependence range is controlled by the boundary conditions. Thus, the performance of the functional least–squares estimator and the functional significance test introduced in this paper is illustrated in both cases, the simulation study and the real–data example considered. A comparative study with the classical approach presented in Worsley et al. [2002] is also achieved for the fMRI data set analysed (freely available at http://www.math.mcgill.ca/keith/fmristat/).

The outline of this paper is as follows. The functional fixed effect model with ARH(1) error term is formulated in Appendix 2. The main results obtained on generalized least–squares estimation of the Hilbert–valued vector of fixed effect parameters, and the functional analysis of variance are also collected in this appendix. Linear hypothesis testing is derived in Appendix 3. The results obtained from the simulation study undertaken are displayed in Appendix 4. Functional statistical analysis of fMRI data is given in Appendix 5. Conclusions and open research lines are provided in Appendix 6. Finally, the Supplementary Material in Appendix 7 introduces the required preliminary elements on eigenvectors and eigenvalues of the Dirichlet negative Laplacian operator on the rectangle, disk and circular sector.

2 Multivariate Hilbert–valued fixed effect model with ARH(1) error term

This section provides the extended formulation of the multivariate Hilbert–valued fixed effect model studied in Ruiz-Medina [2016], to the case where the correlated functional components of the error term satisfy an ARH(1) state equation. In that formulation, compactly supported non–separable autocovariance and cross–covariance kernels are considered for the functional error components, extending the separable case studied in Ruiz-Medina [2016].

Denote by $H$ a real separable Hilbert space with the inner product $\left\langle\cdot,\cdot\right\rangle_{H},$ and the associated norm $\|\cdot\|_{H}.$ Let us first introduce the multivariate Hilbert–valued fixed effect model with ARH(1) error term

[TABLE]

where $\mathbf{X}$ is a real-valued $n\times p$ matrix, the fixed effect design matrix,

[TABLE]

represents the vector of fixed effect parameters,

[TABLE]

is the $H^{n}$-valued Gaussian response, with ${\rm E}\left\{\mathbf{Y}\right\}=\mathbf{X}\boldsymbol{\beta}$. The $H^{n}$-valued error term

[TABLE]

is assumed to be an ARH(1) process on the basic probability space $(\Omega,\mathcal{A},\mathcal{P})$; i.e., a stationary in time Hilbert–valued Gaussian process satisfying (see Bosq [2000])

[TABLE]

where ${\rm E}\left\{\varepsilon_{m}\right\}=0,$ for each $m\in\mathbb{Z},$ and $\rho$ denotes the autocorrelation operator of the error process $\varepsilon,$ which belongs to the space of bounded linear operators on $H.$ Here, $\nu=\left\{\nu_{m},\ m\in\mathbb{Z}\right\}$ is assumed to be a Gaussian strong white noise; i.e., $\nu$ is a Hilbert–valued zero–mean stationary process, with independent and identically distributed components in time, and with $\sigma^{2}={\rm E}\left\{\|\nu_{m}\|_{H}^{2}\right\}<\infty,$ for all $m\in\mathbb{Z}.$ Thus, in (2), the components of the vector error term $[\varepsilon_{1}(\cdot),\dots,\varepsilon_{n}(\cdot)]^{T}$ corresponding to observations at times $t_{1},\dots,t_{n},$ obey the functional state equation (3), under suitable conditions on the point spectrum of the autocorrelation operator $\rho$. Hence, the non–null functional entries of the matrix covariance operator $\mathbf{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}$ of

[TABLE]

are then constituted by the elements located at the three main diagonals. Specifically,

[TABLE]

and

[TABLE]

where $R_{1}$ and $R_{1}^{*}$ denote, respectively, the cross–covariance operator and its adjoint for the ARH(1) process $\varepsilon=\left\{\varepsilon_{i},\ i\in\mathbb{Z}\right\},$ and $R_{0}$ represents its autocovariance operator. Note that, in this appendix, it is assumed that $\rho$ is sufficiently regular. In particular, $\rho$ is such that $\|\rho^{2}\|_{\mathcal{L}(H)}\simeq 0$.

Equivalently, the matrix covariance operator $\mathbf{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}$ is given by

[TABLE]

where $0_{H}$ denotes the approximation by zero in the corresponding operator norm, given the conditions imposed on $\rho$.

In the space $\mathcal{H}=H^{n},$ we consider the inner product

[TABLE]

It is well–known that the autocovariance operator $R_{0}$ of an ARH(1) process is in the trace class (see [Bosq, 2000, pp. 27–36]). Therefore, it admits a diagonal spectral decomposition

[TABLE]

in terms of a complete orthogonal eigenvector system $\left\{\phi_{k},\ k\geq 1\right\},$ defining in $H$ a resolution of the identity $\displaystyle\sum_{k=1}^{\infty}\phi_{k}\otimes\phi_{k}$. Here, for each $k\geq 1,$ $\lambda_{k}=\lambda_{k}(R_{0})$ is the $k$–th eigenvalue of $R_{0},$ with $R_{0}\left(\phi_{k}\right)=\lambda_{k}(R_{0})\phi_{k}.$ The following series expansion then holds, in the mean–square sense:

[TABLE]

where $\eta_{k}(i)=\frac{1}{\sqrt{\lambda_{k}}}\langle\varepsilon_{i},\phi_{k}\rangle_{H},$ for $k\geq 1$ and $i\in\mathbb{N}.$

The following assumption is made:

Assumption A0. The standard Gaussian random variable sequences $\{\eta_{k}(i),\ k\geq 1,\leavevmode\nobreak\ i\in\mathbb{N}\},$ with, for each $k\geq 1,$

[TABLE]

for every $i\in\mathbb{N},$ satisfy the following orthogonality condition, for every $i,j\in\mathbb{N},$

[TABLE]

where $\delta$ denotes the Kronecker delta function, and

[TABLE]

Under Assumption A0, the computation of the generalized least–squares estimator of $[\beta_{1}(\cdot),\dots,\beta_{p}(\cdot)]^{T}$ is achieved by projection into the orthogonal basis of eigenvectors $\left\{\phi_{k},\ k\geq 1\right\}$ of the autocovariance operator $R_{0}$ of the ARH(1) process $\varepsilon=\left\{\varepsilon_{i},\ i\in\mathbb{Z}\right\}.$ Denote by $\boldsymbol{\Phi}^{*}$ the projection operator into the eigenvector system $\left\{\phi_{k},\ k\geq 1\right\},$ acting on a vector function $\mathbf{f}\in\mathcal{H}=H^{n}$ as follows:

[TABLE]

where $\boldsymbol{\Phi}\boldsymbol{\Phi}^{*}=\mathbf{Id}_{\mathcal{H}=H^{n}},$ with

[TABLE]

For $\mathbf{A}=\left\{A_{i,j}\right\}_{i=1,\ldots,n}^{j=1,\ldots,n}$ be a matrix operator such that, for each $i,j=1,\dots,n,$ its functional entries are given by

[TABLE]

with $\displaystyle\sum_{k=1}^{\infty}\gamma_{kij}^{2}<\infty.$ The following identities are straightforward:

[TABLE]

where, for each $k\geq 1,$ the entries of $\boldsymbol{\Gamma}_{k}$ are $\Gamma_{kij}=\gamma_{kij},$ for $i,j=1,\dots,n.$

Applying (14)–(15), we directly obtain

[TABLE]

where, for each $k\geq 1$, $\boldsymbol{\Lambda}_{k}=\boldsymbol{\Phi}^{*}_{k}\mathbf{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}\boldsymbol{\Phi}_{k}$ is given by

[TABLE]

with $\boldsymbol{\Lambda}_{k}^{-1}$ denoting its inverse matrix.

Remark 1

In Appendix 4, we restrict our attention to the functional error model studied in Ruiz-Medina [2016], considering the Hilbert–valued stochastic partial differential equation system framework. In that framework, matrices $\left\{\boldsymbol{\Lambda}_{k},\ k\geq 1\right\},$ are known, since they are defined from the eigenvalues of the differential operators involved in the equation system. Particularly, in that section, for each $k\geq 1,$ matrix $\boldsymbol{\Lambda}_{k}$ is considered to have entries $\Lambda_{kij}$ given by

[TABLE]

with

[TABLE]

and $(-\boldsymbol{\Delta}_{D_{l}})$ representing the Dirichlet negative Laplacian operator on domain $D_{l},$ for $l=1$ (the rectangle), $l=2$ (the disk) and $l=3$ (the circular sector). However, in practice, as shown in Appendix 5 in the analysis of fMRI data, matrices $\left\{\boldsymbol{\Lambda}_{k},\ k\geq 1\right\},$ are not known, and should be estimated from the data. Indeed, in that real–data example, we approximate the entries of $\left\{\boldsymbol{\Lambda}_{k},\ k\geq 1\right\},$ from the coefficients (eigenvalues and singular values), that define the diagonal spectral expansion of the empirical autocovariance $\widehat{R_{0}}$ and cross covariance $\widehat{R_{1}}$ operators, given by (see Bosq [2000])

[TABLE]

We also consider here the following semi–orthogonal condition for the non-square design matrix $\mathbf{X}:$

Assumption A1. The fixed effect design matrix $\mathbf{X}$ is a semi–orthogonal non–square matrix. That is,

[TABLE]

Remark 2

Assumption A1* implies (see Ruiz-Medina [2016])*

[TABLE]

The generalized least–squares estimation of $[\beta_{1}(\cdot),\ldots,\beta_{p}(\cdot)]^{T}$ is achieved by minimizing the loss quadratic function in the norm of the Reproducing Kernel Hilbert Space (RKHS norm). Note that, for an $\mathcal{H}$–valued zero–mean Gaussian random variable with autocovariance operator $R_{Z},$ the RKHS of $Z$ is defined by

[TABLE]

(see, for example, Prato and Zabczyk [2002]).

From equation (16) we get

[TABLE]

where, in the last identity, for each $k\geq 1,$ matrix $\widehat{\boldsymbol{\Lambda}}_{k}$ represents the empirical counterpart of $\boldsymbol{\Lambda}_{k},$ constructed from the eigenelements of $\widehat{R_{0}},$ $\widehat{R_{1}}$ and $\widehat{R_{1}}^{\ast},$ considered when $R_{0}$ and $R_{1}$ are unknown. Here,

[TABLE]

The minimum of equation (20) is attached if, for each $k\geq 1,$ the expectation

[TABLE]

is minimized, with, as before, $\boldsymbol{\Lambda}_{k}^{-1}$ defining the inverse of matrix $\boldsymbol{\Lambda}_{k}$ given in (17) (and approximated by $\widehat{\boldsymbol{\Lambda}}_{k},$ when $R_{0}$ and $R_{1}$ are unknown). That is,

[TABLE]

and given by

[TABLE]

in the case where $R_{0}$ and $R_{1}$ are unknown. Here, $\mathbf{Y}_{k}=\boldsymbol{\Phi}^{*}_{k}\left(\mathbf{Y}\right)$ is the vector of projections into $\phi_{k}$ of the components of $\mathbf{Y},$ for each $k\geq 1.$

In the remaining of this section, we restrict our attention to the case where $R_{0}$ and $R_{1}$ are known. In that case,

[TABLE]

The estimated response is then given by $\widehat{\mathbf{Y}}=\mathbf{X}\widehat{\boldsymbol{\beta}}$. Under Assumption A1,

[TABLE]

i.e., $\widehat{\boldsymbol{\beta}}\in H^{p}$ almost surely (see Ruiz-Medina [2016] for more details).

Remark 3

In the case where $R_{0}$ and $R_{1}$ are unknown, under the conditions assumed in [Bosq, 2000, Corollary 4.2, pp. 101–102], strong consistency of the empirical autocovariance operator $\widehat{R_{0}}$ holds. Moreover, under the conditions assumed in [Bosq, 2000, Theorem 4.8, pp. 116–117], the empirical cross–covariance operator $\widehat{R_{0}}$ is strongly–consistent. Therefore, the plug–in functional parameter estimator (21) satisfies (22), for $n$ sufficiently large.

The Functional Analysis of Variance of model in (2)–(3) can be achieved as described in Ruiz-Medina [2016]. Specifically, a linear transformation of the functional data should be considered, for the almost surely finiteness of the functional components of variance, in the following way:

[TABLE]

where $\mathbf{W}$ is such that

[TABLE]

and satisfies

[TABLE]

Here, for each $k\geq 1,$ $\boldsymbol{\Lambda}_{k}$ is defined in (17). The functional components of variance associated with the transformed model (23) are then given by

[TABLE]

where $\mathbf{M}_{k}=\mathbf{Id}_{n\times n}-\mathbf{X}\left(\mathbf{X}^{T}\boldsymbol{\Lambda}_{k}^{-1}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\boldsymbol{\Lambda}_{k}^{-1},$ for each $k\geq 1.$

The statistics

[TABLE]

provides information on the relative magnitude between the empirical variability explained by the functional transformed model and the residual variability (see Appendix 4).

3 Significance test from the Cramér–Wold’s Theorem

In Ruiz-Medina [2016], a linear functional statistical test is formulated, with explicit definition of the probability distribution of the derived functional statistics under the null hypothesis:

[TABLE]

against

[TABLE]

where $\mathbf{C}\in H^{m}$ and

[TABLE]

is a matrix operator such that its functional entries $\mathbf{K}=\left\{K_{ij}\right\}_{i=1,\dots,m}^{j=1,\dots,p},$ are given, for each $f,g\in H,$ by

[TABLE]

In particular,

[TABLE]

with

[TABLE]

At level $\alpha,$ there exists a test $\psi$ given by:

[TABLE]

where

[TABLE]

The constant $C(H_{0},\alpha)$ is such that

[TABLE]

where the probability distribution $\mathbf{F}$ on $\mathcal{H}=H^{n}$ has characteristic functional given in [Ruiz-Medina, 2016, Proposition 4, Eq. (66)].

Alternatively, as an application of [Cuesta-Albertos et al., 2007, Theorem 4.1], a multivariate version of the significance test formulated in Cuesta-Albertos and Febrero-Bande [2010] is considered here, for the fixed effect parameters (see, in particular, [Cuesta-Albertos and Febrero-Bande, 2010, Theorem 2.1]. Specifically, we consider

[TABLE]

for $\mathbf{h}=(h,\dots,h)^{T}_{p\times 1}$ defining a random vector in $H^{p},$ with $h$ generated from a zero–mean Gaussian measure $\mu$ in $H,$ with trace covariance operator $R_{\mu}$ (see, for example, Prato and Zabczyk [2002]). Here,

[TABLE]

$\mathbf{K}$ is given by

[TABLE]

and $\mathbf{C}$ is a null $\left(p-1\right)\times 1$ functional vector; i.e.,

[TABLE]

From equations (27)–(28), for any $\left(p\times 1\right)$–dimensional functional random vector $\mathbf{h}=(h,\dots,h)_{p\times 1}^{T}$ generated from a Gaussian measure $\mu$ on $H,$ $H_{0}^{\mathbf{h}}$ can then be equivalently expressed as

[TABLE]

The test statistic to contrast (29) is defined as

[TABLE]

where $\mathbf{K}$ and $\mathbf{C}$ are respectively given in equations (27)–(28), and

[TABLE]

with

[TABLE]

Here, $\boldsymbol{\Lambda}_{\mathbf{h}}$ is a $\left(n\times n\right)$–dimensional matrix with entries $\left\{\Lambda_{\mathbf{h}}(i,j)\right\}_{i=1,\dots,n}^{j=1,\ldots,n},$ given by

[TABLE]

where, as before, $\lambda_{k}(R_{ij})$ denotes the $k$–th coefficient in the diagonal expansion of the covariance operator $R_{ij}$ with respect to the basis $\{\phi_{k}\otimes\phi_{k},\ k\geq 1\}$; i.e., in the diagonal expansion

[TABLE]

Note that in the ARH(1) error term case described in Appendix 2, from equation (17),

[TABLE]

Assuming that the autocovariance and cross–covariance operator of the ARH(1) error terms are known, under the null hypothesis $H_{0}^{\mathbf{h}},$ the conditional distribution of $T_{h}$ in (30), given $Y=h,$ is a chi–square distribution with $p-1$ degrees of freedom. Here, $Y$ is a zero-mean $H$–valued random variable with Gaussian probability measure $\mu$ on $H,$ having trace covariance operator $R_{\mu}.$ Note that the last assertion directly follows from the fact that, in equation (31), the conditional distribution of $\widehat{\boldsymbol{\beta}}(\mathbf{h})$ given $Y=h$ is

[TABLE]

with $\mathbf{Q}_{\mathbf{h}}$ being introduced in equation (31); i.e., the conditional distribution of $\widehat{\boldsymbol{\beta}}(\mathbf{h})$, given $Y=h$, is a multivariate Gaussian distribution with mean vector $\boldsymbol{\beta}(\mathbf{h})$ and covariance matrix $\mathbf{Q}_{\mathbf{h}}.$

From [Cuesta-Albertos et al., 2007, Theorem 4.1] and [Cuesta-Albertos and Febrero-Bande, 2010, Theorem 2.1], if

[TABLE]

fails, then, for $\mu$-almost every function $h\in H,$ $H_{0}^{\mathbf{h}}$ in (26), or equivalently in (29), also fails. Thus, a statistical test at level $\alpha$ to test $H_{0}^{\mathbf{h}}$ is a statistical test at the same level $\alpha$ to test $H_{0}.$

4 Simulation study

In this section, we consider the real separable Hilbert space

[TABLE]

the closure, in the norm of the square integrable functions in $\mathbb{R}^{2},$ of the space of infinitely differentiable functions with compact support contained in $D_{l},$ for each $l=1,2,3.$ We restrict our attention to the family of error covariance operators given in (LABEL:A4:covopf). Thus, for each $i,j=1,\dots,n,$

[TABLE]

where $\delta_{i,j}^{*}=1-\delta_{i,j},$ and $\delta_{i,j}$ is the Kronecker delta function. As before, for each $i,j=1,\dots,n$ and $k\geq 1,$

[TABLE]

Note that the above error covariance operator models correspond to define, for $i=1,\dots,n,$ the functional Gaussian error component $\varepsilon_{i}$ as the solution, in the mean–square sense, of the stochastic partial differential equation

[TABLE]

with $\xi_{i}$ being spatial Gaussian white noise on $L^{2}(D_{l}),$ for $l=1,2,3.$

To approximate

[TABLE]

$\nu$ samples are generated for the computation of

[TABLE]

the empirical functional mean–square error ${\rm EFMSE}_{\boldsymbol{\beta}}$ associated with the functional estimates

[TABLE]

of $\boldsymbol{\beta},$ where $L$ is the number of nodes considered in the regular grid constructed over the domains $\left\{D_{l},\ l=1,2,3\right\}.$

Also, we will compute the following statistics:

[TABLE]

where

[TABLE]

and

[TABLE]

with $\|\cdot\|_{\infty}$ denoting the $L^{\infty}$–norm.

Let

[TABLE]

be the generated functional samples. The empirical approximation of

[TABLE]

with ${\rm FMSE}_{\mathbf{Y}}$ being the ${\rm FMSE}$ of $\mathbf{Y},$ can be computed as follows:

[TABLE]

Also, we will consider the statistics

[TABLE]

where

[TABLE]

for $i=1,\dots,n,$ $j=1,\dots,L,$ and $v=1,\dots,\nu.$

In the following numerical examples, the functional analysis of variance is implemented from a transformed functional data model, considering the matrix operator $\mathbf{W}$ such that, for each $k\geq 1,$ $\boldsymbol{\Phi}^{*}_{k}\mathbf{W}=\mathbf{W}_{k}$ compensates the divergence of the eigenvalues of $\boldsymbol{\Lambda}_{k}^{-1}.$ Thus, condition (24) is satisfied. Hence, for all $k\geq 1,\mathbf{W}_{k}$ can be defined as

[TABLE]

where $\boldsymbol{\Omega}\left(\mathbf{W}_{k}\right)=diag\left(\omega_{k11},\dots,\omega_{knn}\right)$ denoting a diagonal matrix, which elements are defined by

[TABLE]

under

[TABLE]

We have chosen $a_{k}=k^{2}.$ Here, for each $k\geq 1,$ $\boldsymbol{\Psi}_{k}$ denotes the projection operator into the system $\left\{\psi_{lk},l=1,\ldots,n\right\}$ of eigenvectors of matrix $\boldsymbol{\Lambda}_{k},$ and $\left\{\omega_{i}\left(\boldsymbol{\Lambda}_{k}\right),\ i=1,\ldots,n\right\}$ are the associated eigenvalues (see Ruiz-Medina [2016]).

In practice, the infinite series defining the generalized least–squares estimator, and the functional components of variance is truncated at $TR.$ Specifically, in the rectangle, we work with a two–dimensional truncation parameter $TR=TR_{1}\times TR_{2},$ and, for circular domains, we fix a one–dimensional parameter (the order $k$ of Bessel functions), thus, $TR_{1}=1,$ and move the second truncation parameter associated with the radius $R$ (see Appendices 7.2–7.3). We then have

[TABLE]

From the transformed model (23), the finite–dimensional approximations (37)–(41) of $\widetilde{SSE},$ $\widetilde{SST},$ and $\widetilde{SSR},$ respectively, are computed to obtain the values of the statistics (25), reflecting the relative magnitude between the empirical functional variability explained by the model and the residual variability.

In the computation of the test statistics $T_{\boldsymbol{h}},$ a truncation order is also considered in the calculation of the elements defining matrix $\boldsymbol{\Lambda}_{\boldsymbol{h}}.$

In all the subsequent sections, the truncation order $TR$ has been selected according to the following criteria:

• (i)

The percentage of explained functional variance. In all the subsequent numerical examples, the $TR$ values considered always ensure a percentage of explained functional variance larger or equal than $95\%.$

• (ii)

The rate of convergence to zero of the eigenvalues of the covariance operators, defining the functional entries of the matrix covariance operator of the $H^{n}$–valued error term. Specifically, in the simulation study undertaken, according to the asymptotic order (rate of convergence to zero) of such eigenvalues, we have selected the optimal $TR$ to remove divergence of the spectra of the corresponding inverse covariance operators.

• (iii)

The functional form of the eigenvectors, depending on the geometry of the domain and the Dirichlet conditions on the boundary. Small truncation orders or values of $TR$ are considered, when fast decay velocity to zero is displayed at the boundary, by the common eigenvectors of the autocovariance operators of the error components, since, in that case, the error dependence range is shorter.

Summarizing, lower truncation orders are required when a fast decay velocity to zero is displayed by the covariance kernel eigenvalues, since a sufficient percentage of explained variability is achieved with a few terms. Note that larger truncation orders can lead to a ill–posed nature of the functional parameter estimation problem, and associated response plug-in prediction. In the subsequent sections, applying criteria (i)–(iii), a smaller number of terms is required in circular domains than in rectangular domains.

Rectangular domain

The $H^{n}$–valued zero–mean Gaussian error term is generated from the matrix covariance operator $\boldsymbol{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}},$ whose functional entries $\left\{\boldsymbol{R}_{\varepsilon_{i}\varepsilon_{j}}\right\}_{i=1,\dots,n}^{j=1,\dots,n},$ are defined in equation (32), with for $i=1,\dots,n,$ $\lambda_{ki}=\lambda_{k}(R_{ii})$ being given in equations (LABEL:A4:covopf) and (LABEL:A4:19). Specifically, $\{\phi_{k},\ k\geq 1\}$ are the eigenvectors of the Dirichlet negative Laplacian operator on the rectangle, associated with the eigenvalues of such an operator (see equation (LABEL:A4:19) in the Supplementary Material in Appendix 7), arranged in decreasing order of their modulus magnitude.

Let us now define the scenarios studied for the rectangular domain

[TABLE]

where $\nu=20$ functional samples of size $n=200$ have been considered, for a given semi–orthogonal design matrix

[TABLE]

These scenarios are determined from the possible values of the vector variable $(P_{i},u,C_{i}),$ where $P_{i}$ refers to the number of components of $\boldsymbol{\beta},$ specifically, for $i=1,$ $p=4$ components, and for $i=2,$ $p=9$ components. Here, $u$ takes the values $a,b,c,d$ respectively corresponding to the truncation orders $TR=16$ ($u=a$), $TR=36$ ($u=b$), $TR=64$ ($u=c$) and $TR=144$ ($u=d$). In addition, $\left\{C_{i},\ i=1,2\right\}$ indicate the shape of $\boldsymbol{\beta}.$ Specifically, we have considered

• •

$\beta_{s}\left(x,y\right)=\sin\left(\frac{\pi sx_{b_{1}}}{l_{1}}\right)\sin\left(\frac{\pi sy_{b_{2}}}{l_{2}}\right)$ (C1)

• •

$\beta_{s}\left(x,y\right)=\cos\left(\frac{x_{b_{1}}+x_{a_{1}}}{l_{1}}\right)\cos\left(\frac{y_{b_{2}}+y_{a_{2}}}{l_{2}}\right)$  (C2),

where

[TABLE]

and $s=1,\dots,p$.

A summary of the generated and analysed scenarios are displayed in Table 1 below.

In Table 1, $h_{x}$ and $h_{y}$ refer to the discretization step size at each dimension. In the cases (P1,a,C1) and (P2,a,C2), a generation of a functional value (surface) of the response is respectively represented in Figures 1–2.

Figures 3–4 below show the respective functional estimates $\widehat{\mathbf{Y}}=\mathbf{X}\widehat{\boldsymbol{\beta}}$ of the responses displayed in Figures 1–2 above.

The statistics (33)–(34) are evaluated in all the cases displayed in Table 1 (see Tables 2–3 for the statistics $L_{\boldsymbol{\beta}}^{\infty}$ and $L_{\mathbf{Y}}^{\infty},$ respectively).

As expected, the results displayed in Table 2, corresponding to the empirical functional mean quadratic errors associated with the estimation of $\boldsymbol{\beta},$ are less than the ones obtained in Table 3 for the response, with order of magnitude $10^{-3}$ in all the scenarios generated. In Table 3, we can appreciate a better performance of the generalized least–squares estimator for the higher truncation orders. However, we have to note that, even for the smallest truncation order considered; i.e., for $TR=4\times 4=16,$ a good performance is observed with associated empirical functional mean quadratic errors having order of magnitude $10^{-2}$ in all the cases displayed in Table 1 (see the above truncation order criteria (i)–(iii)). It can also be observed that the number of components of parameter $\boldsymbol{\beta},$ and their functional shapes do not affect the accuracy of the least–squares generalized estimations of the functional values of the response. It can also be observed that the number of components of parameter $\boldsymbol{\beta},$ and their functional shapes do not affect the accuracy of the least–squares generalized estimations of the functional values of the response.

The statistics (25) is now computed, as an empirical approximation of the relative magnitude between the explained functional variability and the residual variability, after fitting the transformed Hilbert–valued fixed effect model (23). The results obtained are given in Table 4. It can be observed that, in all the cases studied, the explained functional variability exceeds the residual functional variability. The truncation order, the number of components of $\boldsymbol{\beta},$ and the functional shape of such components do not substantially affect the goodness of fit of the transformed Hilbert–valued fixed effect model in (23).

Let us now compute the statistics $T_{\boldsymbol{h}}$ in (30) to contrast the significance of parameter vector $\boldsymbol{\beta}$ in Case C${}_{1},$ when $p=4.$ To apply [Cuesta-Albertos et al., 2007, Theorem 4.1] and [Cuesta-Albertos and Febrero-Bande, 2010, Theorem 2.1], we have generated eight realizations of a Gaussian random function $\boldsymbol{h},$ from the trajectories of the Gaussian random field $\xi,$ solution, in the mean–square sense, of the following boundary value problem:

[TABLE]

where $\varsigma$ denotes a zero–mean Gaussian white noise on $L^{2}([-2,3]\times[-2,3])$; i.e., a zero–mean generalized Gaussian process satisfying

[TABLE]

Table 5 below reflects the percentage of successes, for $\alpha=0.05,$ and the averaged $p$–values over the $150$ samples of the response generated with parameter $\boldsymbol{\beta}$ of C1 type having $p=4$ components, and with size $n=150,$ for $TR=4\times 4.$

A high percentage of successes and very small $p$–values are observed in Table 5; i.e., a good performance of the test statistics is observed.

Disk domain

In the disk domain

[TABLE]

the zero–mean Gaussian $H^{n}$–valued error term is generated from the matrix covariance operator $\boldsymbol{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}},$ whose functional entries are defined in equation (32), considering the eigenvectors $\{\phi_{k},\ k\geq 1\}$ of the Dirichlet negative Laplacian operator on the disk (see equation (47) in the Supplementary Material in Appendix 7), arranged in decreasing order of the modulus magnitude of their associated eigenvalues. Specifically, for $i=1,\dots,n,$ $\lambda_{ki}=\lambda_{k}(R_{ii})$ in (32) is defined in equations (LABEL:A4:covopf) and (47). Again, $\nu=20$ functional samples of size $n=200$ of the response have been generated. The cases studied are summarized in terms of the vector (Pi, u, Cj), $i=1,2,$ $j=1,2,3,$ with variable $u=a,b,c,d,e,f.$ Namely, it is considered $u=a$ for $TR=3,$ $u=b$ for $TR=5,$ $u=c$ for $TR=7,$ $u=d$ for $TR=15,$ $u=e$ for $TR=31,$ and $u=f$ for $TR=79.$ Furthermore, Pi indicates the number of components of $\boldsymbol{\beta},$ with $p=4$ for $i=1,$ and $p=9$ for $i=2.$ Finally, the values of C${}_{j},$ $j=1,2,3,$ refer to the shape of the components of $\boldsymbol{\beta},$ defined from their projections, in terms of the following equations:

[TABLE]

Table 6 reflects a summary with all the cases analysed.

Figures 5–6 respectively reflect the generation of a functional value of the response in the cases (P1,c,C1) and (P1,f,C3).

The respective generalized least–squares functional estimates are displayed in Figures 7–8.

The empirical functional mean quadratic errors (see equations (33)–(34)) are displayed in Table 7, for the estimation of the functional parameter vector $\boldsymbol{\beta},$ and in Table 8 for the estimation of the response $\mathbf{Y}$. It can be observed, as in the rectangular domain, that the order of magnitude of the empirical functional mean quadratic errors associated with the estimation of $\boldsymbol{\beta}$ is of order $10^{-3},$ and for the estimation of the response is $10^{-2}$. However, the number of terms considered is less than in the case of the rectangle; i.e., a finite dimensional space with lower dimension than in the rectangle is required, according to criterion (iii) reflected in Appendix 4. It can also be appreciated that the number of components of $\boldsymbol{\beta}$ does not substantially affect the accuracy of the estimates.

The statistics (25) is now computed (see Table 9), as an empirical approximation of the relative magnitude between the explained functional variability and the residual variability, after fitting the transformed Hilbert-valued fixed effect model (23). It can be noticed that the values of $\frac{\widetilde{SSR}}{\widetilde{SST}}$ are very close to one in all the scenarios analysed. This fact induces large values of (25) (see Table 9), since

[TABLE]

It can be observed, one time more, from criterion (iii), reflected in Appendix 4, that the boundary conditions and the geometry of the domain allows in this case a more substantial dimension reduction than in the rectangular domain case, since with lower truncation orders a better model fitting is obtained.

The statistics $T_{\boldsymbol{h}}$ in (30) is computed to contrast the significance of the parameter vector $\boldsymbol{\beta}$ in case C${}_{1},$ with $p=4$ components. Again, eight realizations of Gaussian random functions $\boldsymbol{h}$ are considered, generated from a Gaussian random field $\xi,$ solution, in the mean–square sense, of the following boundary value problem on the disk:

[TABLE]

where $\varsigma$ denotes a zero–mean Gaussian white noise on $L^{2}(D_{25})$; i.e., a zero–mean generalized Gaussian process satisfying

[TABLE]

Table 10 reflects the percentage of successes, for $\alpha=0.05,$ and the averaged $p$–values over the $150$ samples, generated with size $n=150,$ of the functional response having parameter vector $\boldsymbol{\beta}$ of type C1 with $p=4$ components, for $TR=7.$

Table 10 again illustrates a good performance of the statistics $T_{\boldsymbol{h}}$ in (30). Indeed, we can appreciate a high percentage of successes, and very small $p$–values, very close to zero, that support the significance of the functional parameter vector, considered in the generation of the data set analysed.

Circular sector domain

In the circular sector

[TABLE]

of radius $R$ and angle $\pi\theta,$ the zero-mean Gaussian vector error term is generated from the matrix covariance operator $\boldsymbol{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}},$ whose functional entries are defined in equation (32). The eigenvectors $\{\phi_{k},\ k\geq 1\}$ of the Dirichlet negative Laplacian operator on the circular sector are considered (see equation (LABEL:A4:eqcsa) in the Supplementary Material in Appendix 7), arranged in decreasing order of the modulus magnitude of their associated eigenvalues. Specifically, here, $\boldsymbol{R}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}$ is defined in equation (32), with for $i=1,\dots,n,$ $\lambda_{ki}=\lambda_{k}(R_{ii})$ being given in equations (LABEL:A4:covopf) and (LABEL:A4:eqcsa).

As in the above examples, $\nu=20$ functional samples of size $n=200$ are generated. The cases studied are also summarized in terms of the values of the vector $(P_{i},u,C_{j})$, $i=1,2,\leavevmode\nobreak\ u=a,b,c,d,e,f,$ and $j=1,2,3,$ with the values of $u$ having the same meaning as in the disk domain. Again, values of $P_{i}$ provide the number $p$ of components of $\boldsymbol{\beta};$ i.e., $p=4$ if $i=1$, and $p=9$ if $i=2$. The values $C_{1},$ $C_{2}$ and $C_{3}$ respectively correspond to the following functions defining the components of $\boldsymbol{\beta},$ whose projections are given by:

[TABLE]

A summary of the cases analysed is given in Table 11.

Figures 9–10 display the generation of a functional value of the response in the cases (P2,e,C2) and (P1,f,C3), respectively.

The functional estimates obtained from the finite–dimensional approximation of the generalized least–squares estimator of $\boldsymbol{\beta}$ are now given in Figures 11–12, for the cases (P2,e,C2) and (P1,f,C3), respectively.

As in the previous sections, the empirical functional mean quadratic errors, associated with the estimation of $\boldsymbol{\beta}$ and $\mathbf{Y},$ are computed from equations (33)–(34). They are shown in Table 12, for $\boldsymbol{\beta},$ and in Table 13, for $\mathbf{Y}.$

These empirical functional mean quadratic errors are very stable through the different cases considered, and their order of magnitude is again $10^{-3}$ for the parameter $\boldsymbol{\beta},$ and $10^{-2}$ for the response. Here, the results displayed also correspond to the projection into lower finite–dimensional spaces than in the case of the rectangle, according to the functional form of the eigenvectors (see truncation order criterion (iii) in Appendix 4).

Statistics (25) is now computed. Its values are displayed in Table 14. Again, as in the disk, the proportion of explained functional variability is very close to one leading to large values of statistics (25), as it can be observed in Table 14 for all the cases analysed.

The statistics $T_{\boldsymbol{h}}$ in (30) is computed to contrast the significance of the parameter vector $\boldsymbol{\beta}$ in case C1 with $p=4$ functional components. Eight realizations of a Gaussian random function $\boldsymbol{h}$ are considered from a Gaussian random field $\xi,$ solution, in the mean-square sense, of the following boundary value problem on the circular sector

[TABLE]

where $\theta=2/3,$ $\varsigma$ denotes a zero–mean Gaussian white noise on the circular sector such that

[TABLE]

with $L^{2}(CS)$ denoting the space of square–integrable functions on the circular sector. Table 15 reflects the percentage of successes, for $\alpha=0.05,$ and the averaged $p$–values over the $150$ samples, with size $n=150,$ of the response, having C1–type functional parameter vector $\boldsymbol{\beta}$ with $p=4$ components, considering $TR=7.$

Table 15 again confirms the good performance of the test statistics $T_{\boldsymbol{h}},$ showing a high percentage of successes, and very small magnitudes for the averaged $p$–value (almost zero values), according to the significance of the parameter vector $\boldsymbol{\beta}$ considered in the generation of the analysed functional data set.

5 Functional statistical analysis of fMRI data

In this section, we compare the results obtained from the application of the MatLab function fmrilm.m (see Liao et al. [2012] and Worsley et al. [2002]) from fmristat.m function set (available at http://www.math.mcgill.ca/keith/fmristat), with those ones provided by the implementation of our proposed functional statistical methodology, based on the Hilbert-valued fixed effect models with ARH(1) error term above introduced. The fMRI data set analysed is also freely available in AFNI format at http://www.math.mcgill.ca/keith/fmristat/. (AFNI Matlab toolbox can be applied to read such a data set). In the next section, structural information about such fMRI data is provided (see BrikInfo.m Matlab function).

The first step in the statistical analysis of fMRI data is to modeling the data response to an external stimulus. Specifically, at each voxel, denote by $x(t)$ the (noise-free) fMRI response at time $t,$ and by $s(t)$ the external stimulus at that time. It is well–known that the corresponding fMRI response is not instantaneous, suffering a blurring and a delay of the peak response by about $6s$ (see, for example, Liao et al. [2012]). This fact is usually modelled by assuming that the fMRI response depends on the external stimulus by convolution with a hemodynamic response function $h(t)$ (which is usually assumed to be independent of the voxel), as follows:

[TABLE]

Several models have been proposed in the literature for the hemodynamic response function (hrf). For example, the gamma function (see Lange and Zeger [1997]), or the difference of two gamma functions, to model the slight intensity dip after the response has fallen back to zero (see Friston et al. [1998]).

The effects $\left(x_{i,1},\ldots,x_{i,p}\right)$ of $p$ different types of stimuli on data, in scan $i,$ is combined in terms of an additive model with different multiplicative coefficients $\left(\beta_{1},\ldots,\beta_{p}\right)$ that vary from voxel to voxel. The combined fMRI response is then modeled as the linear model (see Friston et al. [1995])

[TABLE]

for each voxel $v.$

An important drift over time can be observed in fMRI time series data in some voxels. Such a drift is usually linear, or a more general slow variation function. In the first case, i.e., for a linear function

[TABLE]

when the drift is not removed, it can be confounded with the fMRI response. Otherwise, it can be added to the estimate of the random noise $\varepsilon,$ which, in the simplest case is assumed to be an AR(1) process at each voxel. In that case, the linear model fitted to the observed fMRI data is usually given by

[TABLE]

for each one of the voxels $v,$ in the real–valued approach presented in Worsley et al. [2002]. In (44),

[TABLE]

where $\{\xi_{i}(v),\ i=1,\dots,n\}$ are $n$ random components of Gaussian white noise in time, for each voxel $v.$ This temporal correlation structure for the noise has sense, under the assumption that the scans are equally spaced in time, and that the error from the previous scan is combined with fresh noise to produce the error for the current scan. In the presented Hilbert–valued approach, a similar reasoning can be applied to arrive to the fixed effect model with ARH(1) error term, introduced in Appendix 2. This model allows the representation of fMRI data in a functional spatially continuous form. Specifically, for the scan $i,$ a continuous spatial variation is assumed underlying to the values of the noise across the voxels, reflected in the functional value of the ARH(1) process, representing the error term. In the same way, the $H$–valued components of the parameter vector $\boldsymbol{\beta}(\cdot)$ provide a continuous model to represent spatial variation over the voxels of the multiplicative coefficients $\beta_{1}(\cdot),\dots,\beta_{p}(\cdot),$ independently of time. Since the fMRI response is subsampled at the $n$ scan acquisition times $t_{1},\dots,t_{n},$ the fixed effect design matrix $\boldsymbol{X},$ constituted by the values of the fMRI response (43) at such times, under the $p$ different types of stimuli considered, has dimension $n\times p.$ Note that in (43) $x$ is assumed to be independent of the voxel, according to the definition of the hrf.

Description of the data set and the fixed effect design matrix

Brain scan measurements are represented on a set of $64\times 64\times 16$ voxels. Each one of such voxels represents a cube of $3.75\times 3.75\times 7\leavevmode\nobreak\ mm.$ At each one of the $16$ depth levels or slices $\left\{S_{i},\ i=1,\dots,16\right\}$, the brain is scanned in $68$ frames, $\left\{Fr_{h},\ h=1,\dots,68\right\}$. Equivalently, for $i=1,\dots,16,$ on the slice $S_{i},$ a $64\times 64$ rectangular grid is considered, where measurements at each one of the $68$ frames are collected.

We restrict our attention to the case $p=2,$ where two type of events are considered, respectively representing scans hot stimulus (with a height $h_{h}$) and scans warm stimulus (with a height $h_{w}$). The default parameters, chosen by Glover [1999], to generate the hrf as the difference of two gamma densities is the row vector $r=\left[5.4,5.2,10.8,7.35,0.35\right]$, where the first and third parameters represent the time to peak of the first and second gamma densities ($\Gamma_{1}$ and $\Gamma_{2}$), respectively; the second and fourth parameters represent the approximate full width at half maximum (FWHM) of the first and second gamma densities, respectively; and the fifth parameter (called also $DIP$) denotes the coefficient of the second gamma density, for more details, see Glover [1999], about modelling the hrf as the difference of two gamma density functions, in the following way:

[TABLE]

Considering $TR_{t}=5$ seconds as the temporal step between each frame $Fr_{h},\leavevmode\nobreak\ h=1,\dots,68$, in which all slices are scanned, frame times will be $Fr_{times}=\left(0,5,10,\dots,330,335\right)$ (see Figure 13). Remark that, for any of the $68$ scans, separated by $TR_{t}=5$ seconds, keeping in mind that the first $4$ frames are removed, $16$ slices $\left\{S_{i},\ i=1,\dots,16\right\}$, are interleaved every $0.3125$ seconds, approximately.

The events matrix $\boldsymbol{E}$, which will be convoluted with the hrf, is a matrix whose rows are the events, and whose columns are the identifier of the event type, the starting event time, the duration of the event, and the height of the response for the event, respectively. In our example, we have considered a block design of $4$ scans rest, $4$ scans hot stimulus, $4$ scans rest, $4$ scans warm stimulus, repeating $4$ times this block with $4$ last scans rest ($68$ scans total). As noted before, we remove the first 4 frames. The hot event is identified by $1$ and the warm event by $2,$ such that $h_{h}=0.5$ and $h_{w}=1$. Event times, for hot and warm stimulus, will be $\left[20,60,\dots,260,300\right]$, since there are $8$ frames between the beginning of events ($4$ frames for the previous event and $4$ frames rest). Then, our events matrix $\boldsymbol{E}$ considered is

[TABLE]

Convolution of matrix $\boldsymbol{E}$, in (45), with the hrf leads to the set of real–valued $64\times 2$ design matrices

[TABLE]

implemented by fmridesgin.m Matlab function (see Figure 14).

Hilbert–valued fixed effect model fitting to FMRI data. A comparative study

The estimation results obtained with the implementation of the classical and Hilbert–valued linear model methodology are now compared. Specifically, in the classical case, from the linear model approach presented in Worsley et al. [2002], we consider a fixed–effect model fitting, in the case where the error term is an AR(1) process, ignoring spatial correlation across the voxels. In particular, the MatLab function fmrilm.m is implemented to fit model (44) to a single run of fMRI data, allowing for spatially varying temporal correlated errors. The parameters of the spatial varying AR(1) models (from voxel to voxel) are estimated from the sample autocorrelation of the residuals, obtained after estimation of the fixed effect parameter by ordinary least–squares, ignoring temporal correlation of the errors, at each voxel. This procedure could be iterated. That is, the estimated autocorrelation coefficient can be used to pre–whitening the data at each voxel. Hence, the fixed effect parameter is estimated by ordinary least-squares, from such data. This iterative estimation procedure can be repeated several times. However, as pointed out in Worsley et al. [2002], such iterations do not lead to a substantial improvement in practice. A variance reduction technique is then applied in Worsley et al. [2002] to the estimated autocorrelation coefficient (reduced bias sample autocorrelation), consisting of spatial smoothing of the sample autocorrelations. This technique reduces variability, although slightly increases the bias.

In this subsection, we also implement the approach introduced in Appendix 2, from the fMRI data set described in Appendix 5.1. As commented before, our approach presents the advantage of providing a continuous spatial description of the variation of the fixed effect parameters, as well as of the parameters characterizing the temporal correlated error term, with autoregressive dynamics. Furthermore, the spatial correlations are also incorporated to our functional statistical analysis, computed from the spatial autocovariance and cross-covariance kernels, respectively defining the operators $R_{0}$ and $R_{1},$ characterizing the functional dependence structure of the ARH(1) error term.

Functional fixed effect model fitting is independently performed at each slice $S_{i},$ for $i=1,\dots,16.$ Specifically, for $i=1,\dots,16,$ as commented before, a real-valued $n\times p,$ with $p=2,$ fixed effect design matrix $\mathbf{X}_{i}$ is considered (see Appendix 5.1). The effects of the two different events studied are combined by the vector of functional fixed effect parameters

[TABLE]

Here, $H^{2}$ is the Hilbert space of $2$–dimensional vector functions, whose components are square–integrable over the spatial rectangular grid considered at each slice. Furthermore, for $i=1,\dots,16,$

[TABLE]

is the $H^{n}$–valued Gaussian fMRI data response, with $n$ representing the number of frames ($n=64,$ since the first $4$ frames are removed because they do not represent steady–state images). In the computation of the generalized least–squares estimate of $\boldsymbol{\beta},$ the empirical matrices $\left\{\widehat{\boldsymbol{\Lambda}}_{k},\ k=1,\dots,TR\right\}$ are computed from the empirical covariance operators (19), where $TR$ is selected according to the required conditions specified, in relation to the sample size $n,$ in Bosq [2000] (see, in particular, [Bosq, 2000, pp. 101–102 and pp. 116–117], and Remark 3).

In the subsequent developments, in the results obtained by applying the Hilbert–valued multivariate fixed effect approach, we will distinguish between cases A and B, respectively corresponding to the projection into two and five empirical eigenvectors. For each one of the $16$ slices, the temporal and spatial averaged empirical quadratic errors, associated with the estimates of the response, computed with the fmrilm.m MatLab function, and with the proposed multivariate Hilbert–valued mixed effect approach, respectively denoted as $EFMSE_{\boldsymbol{Y}_{i}^{fMRI}}$ and $EFMSE_{\boldsymbol{Y}_{i}^{H}},$ are displayed in Tables 16–17.

It can be observed, in Tables 16–17, that the performance of the two approaches is very similar. However, the advantage of the presented approach relies on the important dimension reduction it provides, since, as commented before, we have considered the truncations orders $TR=2$ (Case A) and $TR=5$ (Case B). Note that, for each slice, the parameter vector has dimension $2\times\times(64\times 64),$ in the model fitted by fmrilm.m Matlab function. While the presented approach fits the functional projected model, that, for the the cases A and B studied, is defined in terms of a parameter vector $\boldsymbol{\beta}$ with dimension $2\times 2$ and $2\times 5,$ respectively. Furthermore, the iterative estimation method implemented in fmrilm.m requires several steps, repeated at each one of the $64\times 64$ voxels in the $16$ slices (data pre–whitening, ordinary least-squares estimation of $\boldsymbol{\beta}$, and AR(1) correlation coefficient estimation iterations, jointly with the spatial smoothing of the temporal correlation - reduced bias - parameter estimates).

For the slices $1,$ $5,$ $10$ and $15,$ the temporal averaged (frames $5$–$68$) estimated values of the response, applying fmrilm.m MatLab function, and the fixed effect model with ARH(1) error term, in cases A and B, are respectively displayed in Figures 15–17. The corresponding empirical time-averaged quadratic errors are displayed in Figures 18–20, respectively.

Significance test

We are interested in contrast the significance of the spatial varying parameter vector $\boldsymbol{\beta}$ combining the effects of the two stimulus considered on the overall brain, in its real-valued, and $H^{2}$-valued version. The $F$ statistic in the MatLab function fmrilm.m (fMRI linear model), computed, as before, from a single run of fMRI data, leads to the results reflected in Table 18, on the percentage of brain voxels, where the real–valued fixed effect model with AR(1) term is significative, for each one of the $16$ slices considered.

As described in Appendix 3, for each slice, i.e., for $i=1,\dots,16,$ the value of the conditional chi–squared test statistics $T_{\boldsymbol{h}},$ in equation (30), is computed, from four realizations of a Gaussian random function $\boldsymbol{h},$ generated from a Gaussian random field $\xi,$ satisfying equation (LABEL:A4:eqdlrect) on the rectangle containing each brain slice. As before, the functional response sample size at each slice is $64,$ since the first four frames are discarded. It can be observed, in the numerical results displayed in Table 19, for $TR=16,$ and in Table 20, for $TR=4,$ that the null hypothesis is rejected, in most of the random directions in all the brain slices; i.e., the functional fixed effect model with ARH(1) error term is significative for $\alpha=0.05.$ Note that a very few $p$–values are slightly larger than $\alpha=0.05,$ with very small difference, that could be produced by the numerical errors accumulated, due to the presence of small values to be inverted. Thus, we can conclude the suitability of our approach, to combine the effects of the scans hot stimulus, and the scans warm stimulus, in a functional spatially continuous framework.

Comparing results in Tables 18–20, we can conclude that both methodologies, the one presented in Worsley et al. [2002], and the functional approach introduced here, lead to similar results regarding the significance of the models they propose, respectively based on spatial varying real–valued multiplicative coefficients with AR(1) error term, and Hilbert–valued coefficients with ARH(1) error term.

6 Conclusions

As shown in the simulation study, the boundary conditions affect the decay velocity at the boundary of the covariance kernels, defining the functional entries of the matrix covariance operator of the error term. Thus, the dependence range of the error components is directly affected by the boundary conditions. A better performance of the generalized least–squares estimator of the parameter vector $\boldsymbol{\beta}$ is observed, when a fast continuous decay is displayed by the error covariance kernels close to the boundary, as it is observed in the circular domains. Furthermore, in the simulation study undertaken, and in the real–data problem addressed, a good performance of the computed generalized least–squares estimator, and of the test statistics is observed for low truncation orders. Thus, an important dimension reduction is achieved with the presented approach. Summarizing, the proposed approach allows the incorporation of temporal and spatial correlations in the analysis, with an important dimension reduction.

The derivation of similar results under alternative boundary conditions like Neumann and Robin boundary conditions constitutes an open research problem (see, for example, Grebenkov and Nguyen [2013]). Another important research problem is to address the same analysis under a slow decay of the error covariance kernels at the boundary (see, for example, Frías et al. [2017], Jiang [2012, 2016], Tong [2011], beyond the Gaussian context).

7 Supplementary Material

The eigenvectors and eigenvalues of the Dirichlet negative Laplacian operator on the regular domains defined by the rectangle, disk and circular sector are described here (see, for example, Grebenkov and Nguyen [2013]). It is well–known that the negative Laplacian operator $\left(-\boldsymbol{\Delta}_{D}\right)$ on a regular bounded open domain $D\subset\mathbb{R}^{2},$ with Dirichlet boundary conditions, is given by

[TABLE]

where $\partial D$ is the boundary of $D.$ In the subsequent development, we will denote by $\left\{\phi_{k},\ k\geq 1\right\}$ and $\left\{\lambda_{k}(-\boldsymbol{\Delta}_{D}),\ k\geq 1\right\}$ the respective eigenvectors and eigenvalues of $\left(-\boldsymbol{\Delta}_{D}\right),$ that satisfy

[TABLE]

for $D$ being one of the following three domains:

[TABLE]

and

[TABLE]

Eigenelements of Dirichlet negative Laplacian operator on rectangles

Let us first consider domain

[TABLE]

The eigenvectors $\left\{\phi_{\mathbf{k}},\ \mathbf{k}\in\mathbb{N}^{2}_{*}\right\}$ and eigenvalues $\left\{\lambda_{\mathbf{k}}(-\boldsymbol{\Delta}_{D_{1}}),\ \mathbf{k}\in\mathbb{N}^{2}_{*}\right\}$ of $-\boldsymbol{\Delta}_{D_{1}}$ are given by (see Grebenkov and Nguyen [2013]):

[TABLE]

where $l_{i}=b_{i}-a_{i},$ for $i=1,2.$

Eigenelements of Dirichlet negative Laplacian operator on disks

In general, for the circular annulus

[TABLE]

its rotation symmetry allows us to define $-\boldsymbol{\Delta}_{\widetilde{D}_{2}}$ in polar coordinates as

[TABLE]

The application of variable separation method then leads to the following explicit formula of its eigenfunctions (see, for example, Grebenkov and Nguyen [2013])

[TABLE]

with

[TABLE]

where $\{J_{k}\left(z\right)\}$ and $\{Y_{k}\left(z\right)\}$ are the Bessel functions of order $k$ of first and second kind, respectively,

[TABLE]

are the corresponding eigenvalues, and the sets $\left\{\alpha_{k,h},\ k\geq 1,\leavevmode\nobreak\ h=1,\dots,M(k)\right\}$ and $\left\{c_{k,h},\ k\geq 1,\leavevmode\nobreak\ h=1,\dots,M(k)\right\}$ are defined from the boundary conditions at $r=R$ and $r=R_{0}.$

If we focus on domain $D_{2},$ the disk, i.e., $R_{0}=0$, the coefficients $\left\{c_{k,h},\ k\geq 1,\leavevmode\nobreak\ h=1,\dots,M(k)\right\}$ are set to $0.$ The eigenfunctions then adopt the following expression:

[TABLE]

with eigenvalues

[TABLE]

where $\left\{\alpha_{k,h},\ h=1,\dots,M(k)\right\}$ are the $M(k)$ positive roots of the Bessel function $J_{k}\left(z\right)$ of order $k.$ Note that we can also consider truncation at parameter $M(k)$ for $k\geq 1,$ since this parameter increases with the increasing of the radius $R.$

Eigenelements of Dirichlet negative Laplacian operator on circular sectors

Lastly, we consider domain $D_{3},$ the circular sector of radius $R$ and angle $0<\varphi<\pi\theta.$ The eigenvectors and eigenvalues are given by the following expression (see, for example, Grebenkov and Nguyen [2013]):

[TABLE]

with $M(k)$ and $\left\{\alpha_{k,h},\ k\geq 1,\leavevmode\nobreak\ h=1,\dots,M(k)\right\}$ being given as in the previous section.

Asymptotic behavior of eigenvalues

7.4.1 The rectangle

The functional data sets generated in Appendix 4 must have a covariance matrix operator with functional entries (operators) in the trace class. We then apply the results in Widom [1963] to study the asymptotic order of eigenvalues of the integral equation

[TABLE]

In our case, $V$ is the indicator function on the rectangle, i.e., on domain $D_{1},$ and $l_{\varepsilon_{i}}$ is the covariance kernel defining the square root

[TABLE]

of the autocovariance operator of the Hilbert-valued error component $\left\{\varepsilon_{i},\ i=1,\dots,n\right\},$ with

[TABLE]

Note that with the choice made of functions $V$ and $\left\{l_{\varepsilon_{i}},\ i=1,\dots,n\right\},$ the conditions assumed in Widom [1963] are satisfied. In particular, the following asymptotic holds:

[TABLE]

(see [Widom, 1963, p. 279, Eq. (2)]). Also, in general, the eigenvalues of the Dirichlet negative Laplacian operator on a regular bounded open domain $D$ satisfy

[TABLE]

7.4.2 Asymptotic behavior of zeros of Bessel functions.

As before, $J_{k}\left(z\right)$ denotes the Bessel function of the first kind of order $k.$ Let $\left\{j_{k,h},\ h=1,\dots,M(k)\right\}$ be its $M(k)$ roots. In Elbert [2001], Olver [1951, 1952], it is shown that, for a fixed $h$ and large $k$, the Olver’s expansion holds

[TABLE]

On the other hand, for fixed $k$ and large $h$, the McMahon’s expansion also is satisfied (see, for example, Watson [1966])

[TABLE]

These results will be applied in Appendix 4, in the definition of the eigenvalues of the covariance operators $\left\{\boldsymbol{R}_{\varepsilon_{i}\varepsilon_{i}},\ i=1,\dots,n\right\},$ on the disk and circular sector, to ensure their rapid decay to zero, characterizing the trace operator class.

Acknowledgments

This work has been supported in part by project MTM2015–71839–P (co-funded by Feder funds), of the DGI, MINECO, Spain.