Mean Equality Tests for High-Dimensional and Higher-Order Data with k-Self Similar Compound Symmetry Covariance Structure

Abstract

An extension of the $D^2$ test statistic to test the equality of mean for high-dimensional and k-th order array-variate data using k-self similar compound symmetry (k-SSCS) covariance structure is derived. The k-th order data appear in many scientific fields including agriculture, medical, environmental and engineering applications. We discuss the property of this k-SSCS covariance structure, namely, the property of Jordan algebra. We formally show that our $D^2$ test statistic for k-th order data is an extension or the generalization of the $D^2$ test statistic for second-order data and for third-order data, respectively. We also derive the $D^2$ test statistic for third-order data and illustrate its application using a medical dataset from a clinical trial study of the eye disease glaucoma. The new test statistic is very efficient for high-dimensional data where the estimation of unstructured variance-covariance matrix is not feasible due to small sample size.

1. Introduction

We study the hypotheses testing problems of equality of means for high-dimensional and higher-order (multi-dimensional arrays) data. Standard multivariate techniques such as Hotelling’s ${\mathcal{T}}^2$ distribution with one big unstructured variance-covariance matrix (assuming a large sample size) do not work for these higher-order data, as Hotelling’s ${\mathcal{T}}^2$ distribution cannot incorporate any higher-order information into the test statistic and thus draws wrong conclusions [1]. Higher-order data are formed by representing the additional associations that are inherent in repetition across several dimensions. To obtain a better understanding of higher-order data, we first share a simple and interesting example of higher-order data:

• Traditional multivariate (vector-variate) data are the first-order data. For example, a clinical trial study of glaucoma, where several factors $\left(m_1\right)$ such as intraocular pressure (IOP) and central corneal thickness (CCT) are effective in the diagnosis of glaucoma. This example is an illustration of the $(m_1\times 1)$ vector-variate first-order data.
• When the first-order data are measured at various locations/sites or time points, the data become two-dimensional matrix-variate data, and we name it as second-order data. These data are also recognized as multivariate repeated measures data, or doubly multivariate data; e.g., multivariate spatial data or multivariate temporal data. In the above example of the clinical trial study, an ophthalmologist or optometrist diagnoses glaucoma by measuring IOP and CCT $\left(m_1\right)$ in both the eyes $\left(m_2\right)$. So, we see how the vector-variate first-order dataset discussed in the previous paragraph becomes a $(m_1\times m_2)$ matrix-variate second-order dataset by measuring $m_1$ variables repeatedly over another dimension.
• When the second-order data are measured at various sites, or over various time points, the data become three-dimensional array-variate data, and we name it as third-order data. In addition, these are recognized as triply multivariate data, e.g., multivariate spatio-temporal data or multivariate spatio-spatio data. In the previous example, if the IOP and CCT $(m_1=2)$ are measured in both eyes $(m_2=2)$ as well as over, say, three time points $(m_3=3)$, the dataset would become third-order data.
• When the third-order data are measured at various directions, the data become four-dimensional $(m_1\times m_2\times m_3\times m_4)$ array-variate fourth-order data, e.g., multivariate directo-spatio-temporal data or multivariate directo-spatio-spatio data.
• When the fourth-order data are measured at various depths, the data become five-dimensional $(m_1\times m_2\times m_3\times m_4\times m_5)$ array-variate fifth-order data, and so on, e.g., multivariate deptho-directo-spatio-temporal data.

In the above glaucoma data example, the dataset has $m_1$ variables, and these $m_1$ variables are repeatedly measured over various dimensions adding higher-order information to the dataset. Now, the question is, what is the higher-order information? Higher-order information is embedded in the higher-order covariance structures that are formed by signifying the additional associations that are inherent in repetition of the variables across several dimensions. The other question is how can we measure and capture the higher-order information? For this, one needs to understand how to read these structured higher-order data and how to use the appropriate variance-covariance structure to incorporate the higher-order information that are integral to the higher-order data.

Higher-order data have been studied by many authors in the last 20 years using various variance-covariance structures to reduce the number of unknown parameters, which is very important for high-dimensional data. Second-order data are studied using matrix-variate normal distribution [2,3]. Second-order data can also be analyzed vectorially using a two-separable (Kronecker product) variance-covariance structure [4,5], or a block compound symmetry (BCS), also named a block exchangeable (BE) or a 2-SSCS covariance structure [6]. Two-separable covariance structure for second-order data has two covariance matrices, one for each order of the data; in other words, one covariance matrix for within-subject information and the other covariance matrix for between-subject information. Combining the covariance structures of within-subject information and between-subject information results in a second-order model for second-order data. Ignoring this information often influences the test statistic, and if not properly taken care of this information, test statistic will end up yielding wrong conclusions [1]. To obtain a picture of the third-order data, see [7]. Manceur and Dutilleul [7] used tensor normal distribution with double separable covariance structure. 2-SSCS and 3-SSCS covariance structures are useful tools for the analyses of the second- and third-order datasets, respectively. Manceur and Dutilleul [7] also studied fourth-order data with four-separable covariance structure. In the same way, k-th order data can be analyzed vectorially with some structured variance-covariance matrix to integrate the higher-order information into the model, e.g., k-separable covariance structure [8,9] for the k-th order data. However, k-separable covariance structure may not be appropriate for all datasets; thus, we investigate some other structure, namely, k-SSCS covariance structure (defined in Section 3) for the k-th order data in this article. See [10].

The high-dimensionality of a dataset needs to exploit the structural properties of the data to reduce the number of estimated degrees of freedom for more accurate conclusion for the k-th order data, and k-SSCS covariance structure is one of them. For example, for the third-order glaucoma data, the number of unknown parameters in the $(12\times 12)$-dimensional unstructured variance-covariance matrix is 78, whereas the number of unknown parameters for 3-SSCS covariance structure is just 9 and thus may help in providing the correct information about the true association of the structured third-order data. The data quickly become high-dimensional with the increase in the order of the data, and thus, the variance-covariance matrix becomes singular for small samples, and thus testing of mean is not possible. This necessitates the development of new statistical methods with a suitable structured variance-covariance matrix, which can integrate the existing correlations information of the higher-order data into the test statistic and can take care of the high-dimensionality of the data as well.

Rao [11,12] introduced 2-SSCS covariance structure while classifying genetically different groups. Olkin and Press [13] examined a circular stationary model. The problem of estimation in balanced multilevel models with a block circular symmetric covariance structure was studied by Liang et al. [14]. Olkin [15] studied the hypothesis testing problem of the equality of the mean vectors of multiple populations of second-order data using a 2-SSCS covariance structure, which is reminiscent of a model of Wilks [16]. Arnold [17] studied normal testing problems that mean is patterned when the variance-covariance matrix has a 2-SSCS structure. Arnold [17] also studied multivariate analysis of variance problem when the variance-covariance matrix has a 2-SSCS structure. Arnold [18] later developed linear models with a 2-SSCS structure as the error matrix for one matrix-variate observation. Roy et al. [19] and Žežula et al. [1] studied the hypotheses testing problems on the mean for the second-order data using a 2-SSCS covariance structure. There are few studies on third-order data using the 3-SSCS covariance structure. See Leiva and Roy [20] for classification problems and Roy and Fonseca [21] for linear models with a 3-SSCS covariance structure on the error vectors. Recently, Žežula et al. [22] studied the mean value test for third-order data using a 3-SSCS covariance structure.

A majority of the above-mentioned authors only studied the second-order matrix-variate data and used a 2-SSCS covariance structure where the exchangeability (invariance) property in one factor was present. However, we obtain datasets these days with more than one factor, and the assumption of exchangeability on the levels of factors is appropriate for these datasets. A k-SSCS structured matrix results from the exchangeability property of the $k-1$ factors of a dataset. Employing a 2-SSCS covariance structure would be wrong for the datasets with more than one factor. One may construct a second-order data from a k-th order data by summing the observations, however, it would result in a loss of detailed information of particular characteristics that may be of interest. One may also consider matricization of the k-th order data to a second-order data and then using the 2-SSCS covariance structure, but then, once again, all the correlation information will be wiped out. So, the development of new statistical methods are in demand to handle the k-th order data using k-SSCS variance-covariance matrix.

The aim of this paper is to derive a test statistic for mean for high-dimensional k-th order data using k-SSCS covariance matrix by generalizing $D^2$ test statistics developed in Žežula et al. [1]. In doing so, we exploit the distributions of the eigenblocks of the k-SSCS covariance matrix. We obtain the test statistic to test the mean for one sample case, paired samples case and two independent samples case. We show through Remark 2 that our generalized $D^2$ test statistic for the k-th order data generalizes the test statistic for the second-order data, and we derive the test statistic for the third-order data, which is largely motivated by the work of Žežula et al. [22] and in Remark 2 as well.

This article is organized as follows. In Section 2, we set up some preliminaries about some matrix notations and definitions related to block matrices. Section 3 has the definition of k-SSCS covariance matrix. Section 3 has properties of the k-SSCS covariance matrix, such as Jordan algebra. Section 4 discusses the estimation of the eigenblocks and their distributions. The test for the mean for one population is proposed in Section 5. Tests for the equality of means for two populations are proposed in Section 6, and an example of a dataset exemplifying our proposed method is presented in Section 7. Finally, Section 8 concludes with some discussion and the scope for the future research.

2. Preliminaries

Let $m_g,$ for $g=1,\dots ,k,$ be natural numbers greater than $1,$ and $p_{i,j}$ be given by:

$$p_{i,j}=\left\{\begin{array}{ccc}{\displaystyle \prod _{g=i}^j}{m}_g& \mathrm{if}& 0\le j-i<k\\ 1& \mathrm{if}& j-i=-1\\ 0& \mathrm{if}& j-i=-2\end{array}\right.,$$

(1)

with $i=1,\dots ,k.$ We denote by $F_g$ the set $F_g=\left\{1,\dots ,m_g\right\}$, for $g=1,\dots ,k.$

Definition 1.

We say that a matrix $\underset{p_{1,k}\times p_{1,k}}{{\mathit{A}}_k}$ is a k-th order block matrix according to the factorization $p_{1,k}={\displaystyle \prod _{g=1}^k}{m}_g$ to point out that it can be expressed as k different “natural" partitioned matrix forms, that is:

$$\underset{p_{1,k}\times p_{1,k}}{{\mathit{A}}_k}={\left(\underset{p_{1,j}\times p_{1,j}}{{\mathit{A}}_{f_k,\dots ,f_{k+1-\left(k-j\right)};f_k^{*},\dots ,f_{k+1-\left(k-j\right)}^{*}}}\right)}_{f_k,f_k^{*}\in F_k;\dots ;f_{j+1},f_{j+1}^{*}\in F_{j+1}}:j=0,\dots ,k-1.$$

Note that for the case $j=0$, the matrix $\mathit{A}$ is a $(k\times k)-$dimensional matrix with $1\times 1$ blocks. Clearly, both $m_1\ge 2$ and $m_2\ge 2$ for second-order data, and $m_1\ge 2$, $m_2\ge 2$ and $m_3\ge 2$ for third-order data, and so on. Next, we define matrix operators that will be useful tools in working with these k-th order block matrices, where $k\ge 2$. Let ${\mathcal{M}}_{p_{1,g}}$ denote the set of $p_{1,g}\times p_{1,g}$-matrices.

Definition 2.

Let ${\mathrm{BS}}_{p_{1,g}}$ and ${\mathrm{BT}}_{p_{1,g}}$ denote the $p_{1,g}-Sum$ and $p_{1,g}-Trace$ block operators from ${\mathcal{M}}_{p_{1,h}}$ to ${\mathcal{M}}_{p_{1,g}}$ for $1\le g\le h\le k$, respecitvely, where ${\mathcal{M}}_{p_{1,h}}$ will always be evident from the context. These block operators applied to a matrix

$$\begin{array}{c}\hfill \underset{p_{1,g}{p}_{g+1,h}\times p_{1,g}{p}_{g+1,h}}{{\mathit{G}}_k}={\left(\underset{p_{1,g}\times p_{1,g}}{{\mathit{G}}_{\mathit{f},{\mathit{f}}^{*}}}\right)}_{\mathit{f},{\mathit{f}}^{*}\in F_{h,g+1}=F_h\times \cdots \times F_{g+1}=\stackrel{h-g}{\underset{j=1}{\times }}{F}_{h+1-j}}\end{array}$$

give the following $p_{1,g}\times p_{1,g}$-matrices:

$$\begin{array}{c}\hfill {\mathrm{BS}}_{p_{1,g}}\left({\mathit{G}}_k\right)=\sum _{\mathit{f}\in F_{h,g+1}}\sum _{{\mathit{f}}^{*}\in F_{h,g+1}}{\mathit{G}}_{\mathit{f},{\mathit{f}}^{*}}and{\mathrm{BT}}_{p_{1,g}}\left({\mathit{G}}_k\right)=\sum _{\mathit{f}\in F_{h,g+1}}{\mathit{G}}_{\mathit{f},\mathit{f}}.\end{array}$$

The subindex $p_{1,g}$ in these block matrix operators represents $p_{1,g}\times p_{1,g}-$dimensional blocks in a partitioned square matrix $\mathit{G}$, and thus their use results in $p_{1,g}\times p_{1,g}-$dimensional matrices. Many useful properties of these block operators, which we will use later in this article, are examined in Leiva and Roy [10]. For any natural number $a>1,$ we use the following additional notations:

$$\begin{array}{ccc}\hfill \underset{a\times a}{{\mathit{Q}}_a}& =& {\mathit{I}}_a-{\mathit{P}}_a,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \mathrm{and}{\mathit{P}}_a& =& \frac{1}{a}{\mathit{J}}_a,\hfill \end{array}$$

(3)

where ${\mathit{J}}_a={\mathbf{1}}_a{\mathbf{1}}_a^{\prime }$, with ${\mathbf{1}}_a$ be the $a\times 1$ vector of ones, and ${\mathit{I}}_a=\left[{\mathit{e}}_{a,1},\cdots ,{\mathit{e}}_{a,a}\right]$ being the $a\times a-$ identity matrix with ${\mathit{e}}_{a,i}$ the ith column vector of ${\mathit{I}}_a.$ Observe that ${\mathit{P}}_a$ and ${\mathit{Q}}_a$ are idempotent matrices and mutually orthogonal to each other, that is:

$${\left({\mathit{P}}_a\right)}^2={\mathit{P}}_a,{\left({\mathit{Q}}_a\right)}^2={\mathit{Q}}_a\mathrm{and}{\mathit{P}}_a{\mathit{Q}}_a=\mathbf{0}.$$

For a fixed natural number $k\ge 2,$ let ${\mathit{R}}_{k,j}$ be the $p_{1,k}\times p_{1,k}$-matrix

$${\mathit{R}}_{k,j+1}={\mathit{R}}_{k,j+1}^{*}\otimes {\mathit{I}}_{m_1},$$

(4)

where, the symbol ⊗ represents the Kronecker product operator and for each $j=1,\dots ,k-1,$

$$\begin{array}{ccc}\hfill \underset{p_{2,k}\times p_{2,k}}{{\mathit{R}}_{k,j+1}^{*}}& =& \left(\underset{h=1}{\overset{k-\left(j+1\right)}{⨂}}{\mathit{I}}_{m_{k+1-h}}\right)\otimes {\mathit{Q}}_{m_{j+1}}\otimes \left(\underset{h=k-\left(j-1\right)}{\overset{k-1}{⨂}}{\mathit{P}}_{m_{k+1-h}}\right)\hfill \\ & =& {\mathit{I}}_{p_{j+2,k}}\otimes {\mathit{Q}}_{m_{j+1}}\otimes {\mathit{P}}_{m_{j,}{m}_2},\hfill \end{array}$$

(5)

with

$${\mathit{P}}_{m_i,m_{i^{*}}}=\left\{\begin{array}{ccc}{\displaystyle \underset{h=k-\left(i-1\right)}{\overset{k-\left(i^{*}-1\right)}{⨂}}}{\mathit{P}}_{m_{k+1-h}}={\mathit{P}}_{m_i}\otimes {\mathit{P}}_{m_{i-1}}\otimes \cdots \otimes {\mathit{P}}_{m_{i^{*}}}& \mathrm{if}& i\ge i^{*}\\ 1& \mathrm{if}& i<i^{*}\end{array}\right.$$

and ${\displaystyle \underset{h=1}{\overset{0}{⨂}}}{\mathit{I}}_{m_{k+1-h}}=1={\displaystyle \underset{h=k}{\overset{k-1}{⨂}}}{\mathit{P}}_{m_{k+1-h}}$. Also, let ${\mathit{R}}_{k,k+1}$ be the $p_{1,k}\times p_{1,k}$-matrix such that

$${\mathit{R}}_{k,k+1}={\mathit{R}}_{k,k+1}^{*}\otimes {\mathit{I}}_{m_1},$$

(6)

where

$$\underset{p_{2,k}\times p_{2,k}}{{\mathit{R}}_{k,k+1}^{*}}=\underset{h=1}{\overset{k-1}{⨂}}{\mathit{P}}_{m_{k+1-h}}={\mathit{P}}_{m_k,m_2}.$$

(7)

3. Properties of the Self Similar Compound Symmetry Covariance Matrix

Let ${\mathit{x}}_{r;\mathit{f}}$ be an $m_1-$variate vector of measurements on the rth replicate (individual) at the $\mathit{f}=\left(f_k,\dots ,f_2\right)\in F=F_{k,2}=F_k\times \cdots \times F_2=\stackrel{2}{\underset{g=k}{\times }}{F}_g$ factor combination. Let ${\mathit{x}}_r$ be the $p_{1,k}={\prod }_{j=1}^k{m}_j-$variate vector of all measurements corresponding to the rth sample unit of the population, that is, ${\mathit{x}}_r={\left(x_{r;1,\dots ,1},\dots ,x_{r,m_k,\dots ,m_1}\right)}^{\prime }$. Thus, the unstructured covariance matrix ${\mathsf{\Gamma }}_k$ has $q=p_{1,k}\left(p_{1,k}+1\right)/2$ unknown parameters that can be large for random values of either of the $m_j$’s. Consequently, if the data are high-dimensional, k-SSCS covariance matrix (defined bellow in Definition 3) with number of unknown parameters $k{m}_1\left(m_1+1\right)/2$ is a good choice if the exchangeable feature is present in the data.

Definition 3.

We say that ${\mathit{x}}_r$ has a k-SSCS covariance matrix if ${{\Gamma }}_k=\mathrm{cov}\left[{\mathit{x}}_r\right]$ is of the form:

$${\mathsf{\Gamma }}_k=\left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{k,j}-{\mathit{U}}_{k,j+1}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{k,k},$$

(8)

where ${\mathit{U}}_{k,j},$ for $j=1,\dots ,k,$ are $m_1\times m_1$-matrices called SSCS-component matrices, with the assumption that ${\mathit{J}}_{p_{2,1}}$ is equal to the real number 1.

An alternative expression to (8) is:

$${\mathsf{\Gamma }}_k=\left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes {\mathit{T}}_{k,j}\right]+\underset{h=1}{\overset{k-1}{⨂}}{\mathit{J}}_{m_{k+1-h}}{\mathit{T}}_{k,k},$$

(9)

with

$$\begin{array}{ccc}\hfill {\mathit{T}}_{k,k}& =& {\mathit{U}}_{k,k},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \mathrm{and}{\mathit{T}}_{k,j}& =& {\mathit{U}}_{k,j}-{\mathit{U}}_{k,j+1}\mathrm{for}j=1,\dots ,k-1\hfill \\ \hfill \mathrm{or\hspace{0.17em}equivalently}{\mathit{U}}_{k,k}& =& {\mathit{T}}_{k,k},\mathrm{and}{\mathit{U}}_{k,i}=\sum _{h=i}^k{\mathit{T}}_{k,h}\mathrm{for}i=1,\dots ,k-1.\hfill \end{array}$$

(11)

The covariance matrix ${\mathsf{\Gamma }}_k$ given in (8) is called k-self similar compound symmetry covariance matrix because if we consider the $p_{1,k}$-dimensional vector $\mathit{x}={\left(x_{1,\dots ,1},\dots ,x_{m_k,\dots ,m_1}\right)}^{\prime }$ with a k-SSCS covariance matrix ${\mathsf{\Gamma }}_k$ and for each fixed $gϵ\left\{2,\dots ,k-1\right\}$, we also consider the partition of $\mathit{x}$ in $p_{1,g}$-subvectors. Then, its corresponding covariance matrix ${\mathsf{\Gamma }}_k$ is partitioned in $p_{1,g}\times p_{1,g}$-submatrices, which is $\left(k+1-g\right)$-SSCS matrix ${\mathsf{\Gamma }}_{k+1-g}^{*}$ (see Leiva and Roy [10]) as follows:

$$\begin{array}{ccc}\hfill \underset{p_{1,k}\times p_{1,k}}{{\mathsf{\Gamma }}_k}& =& {\mathsf{\Gamma }}_{k+1-g}^{*}\hfill \\ & =& \left[\sum _{j=1}^{k+1-g-1}{\mathit{I}}_{p_{j+1,k+1-g}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{k+1-g,j}^{*}-{\mathit{U}}_{k+1-g,j+1}^{*}\right)\right]\hfill \\ & & +{\mathit{J}}_{p_{2,k+1-g}}\otimes {\mathit{U}}_{k+1-g,k+1-g}^{*},\hfill \end{array}$$

where ${\mathit{U}}_{k+1-g,1}^{*}$ is the g-SSCS matrix given by:

$$\begin{array}{ccc}\hfill {\mathit{U}}_{k+1-g,1}^{*}& =& \left\{\sum _{f=1}^{g-1}{\mathit{I}}_{p_{j+1,g}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{k,f}-{\mathit{U}}_{k,f+1}\right)\right\}+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{k,g},\hfill \\ \hfill \mathrm{and}{\mathit{U}}_{k+1-g,j}^{*}& =& \left(\underset{h=1}{\overset{g-1}{⨂}}{\mathit{J}}_{m_{g+1-h}}\right)\otimes {\mathit{U}}_{k,g+j-1}\mathrm{for}jϵ\left\{2,\dots ,k+1-g\right\}.\hfill \end{array}$$

The existence of ${\mathsf{\Gamma }}_k^{-1}$ can be proved using the principle of Mathematical Induction and to derive its expression as well. For the expression of ${\mathsf{\Gamma }}_k^{-1}$, we need matrices ${\mathsf{\Delta }}_{k,j}$ for $j=1,\dots ,k$, which are defined as follows:

$${\mathsf{\Delta }}_{k,j}=\sum _{i=1}^j{p}_{2,i}\left({\mathit{U}}_{k,i}-{\mathit{U}}_{k,i+1}\right),$$

(12)

where ${\mathit{U}}_{k,k+1}=\mathbf{0}$ and $p_{2,1}=1.$ Note that

$${\mathsf{\Delta }}_{k,j}=\left\{\begin{array}{ccc}{\mathit{U}}_{k,1}-{\mathit{U}}_{k,2}& if& j=1\\ {\mathsf{\Delta }}_{k,j-1}+p_{2,j}\left({\mathit{U}}_{k,j}-{\mathit{U}}_{k,j+1}\right)& if& j=2,\dots ,k\end{array}\right..$$

(13)

It can be proved that if matrices ${\mathsf{\Delta }}_{k,j},$$j=1,\dots ,k$ are non-singular, then ${\mathsf{\Gamma }}_k^{-1}$ exists and is given by:

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k^{-1}& =& \left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \frac{1}{p_{2,j}}\left({\mathsf{\Delta }}_{k,j}^{-1}-{\mathsf{\Delta }}_{k,j-1}^{-1}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes \frac{1}{p_{2,k}}\left({\mathsf{\Delta }}_{k,k}^{-1}-{\mathsf{\Delta }}_{k,k-1}^{-1}\right),\hfill \end{array}$$

(see Leiva and Roy [10]), where the symbol ${\mathsf{\Delta }}_{k,0}^{-1}$ indicates the $m_1\times m_1$ zero matrix $({\mathsf{\Delta }}_{k,0}^{-1}={\mathbf{0}}_{m_1\times m_1})$. It is worthwhile to note that the structure of ${\mathsf{\Gamma }}_k^{-1}$ is the same as the structure of ${\mathsf{\Gamma }}_k,$ that is, it has the k-SSCS structure given in (9) with (10) and (11) and

$${\mathsf{\Gamma }}_k^{-1}=\left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes {\mathit{T}}_{k,j}\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{T}}_{k,k},$$

where, in this formula ${\mathsf{\Gamma }}_k^{-1}$, ${\mathit{T}}_{k,j}$ is as follows

$${\mathit{T}}_{k,j}=\frac{1}{p_{2,j}}\left({\mathsf{\Delta }}_{k,j}^{-1}-{\mathsf{\Delta }}_{k,j-1}^{-1}\right),\mathrm{for}\mathrm{each}j=1,\dots ,k.$$

Using a similar inductive arguments, it can be proved that:

$$\left|{\mathsf{\Gamma }}_k\right|=\prod _{j=1}^k{\left|{\mathsf{\Delta }}_{k,j}\right|}^{p_{j+1,k}-p_{j+2,k}},$$

where the matrices ${\mathsf{\Delta }}_{k,j}$ are given by (12), and it is assumed that $p_{k+1,k}=1$ and $p_{k+2,k}=0$. The matrices ${\mathsf{\Delta }}_{k,j},j=1,\dots ,k$ are the k eigenblocks of the k-SSCS covariance structure. See Lemma 4 of Leiva and Roy [10] for proof. The matrix ${\mathsf{\Gamma }}_k$ can be written as the following sum of k orthogonal parts:

$${\mathsf{\Gamma }}_k=\sum _{j=1}^k{\mathit{R}}_{k,j+1}^{*}\otimes {\mathsf{\Delta }}_{k,j}$$

(14)

and if ${\mathsf{\Gamma }}_k^{-1}$ exists, then it can be written as:

$${\mathsf{\Gamma }}_k^{-1}=\sum _{j=1}^k{\mathit{R}}_{k,j+1}^{*}\otimes {\mathsf{\Delta }}_{k,j}^{-1},$$

(15)

where ${\mathit{R}}_{k,j+1}^{*}$ is given in (5), for each $j=1,\dots ,k-1,$ and, for $j=k,$${\mathit{R}}_{k,k+1}^{*}$ is given in (7).

The conventional Hotelling’s $T^2$ statistic to test the mean is based on the unbiased estimate of the unstructured variance-covariance matrix, which follows a Wishart distribution. Nevertheless, the unbiased estimate of the k-SSCS covariance matrix does not follow a Wishart distribution, and thus the test statistic to test the equality of mean does not follow Hotelling’s $T^2$ statistic. We thus make a canonical transformation of the data to block diagonalize the k-SSCS covariance matrix, and show that a scalar multiple of the estimates of the diagonal blocks (eigenblocks) follow independent Wishart distributions and use this property in our advantage to obtain test statistics to test the mean for the k-th order data $(k\ge 2)$. We see from Leiva and Roy [10] that the k-SSCS matrix ${\mathsf{\Gamma }}_k$ given by (8) can be transformed into an $m_1\times m_1$-block diagonal matrix (blocks in the diagonal are $m_1\times m_1$-matrices) by pre- and post-multiplying ${\mathsf{\Gamma }}_k$ by appropriate orthogonal matrices.

For $1\le j<g<i\le k,$ let ${\mathit{I}}_{m_i,m_j}$ denote the $m_i{m}_{i-1}\cdots m_j$ identity matrix, that is:

$${\mathit{I}}_{m_i,m_j}={\mathit{I}}_{m_i{m}_{i-1}\cdots m_j}={\mathit{I}}_{m_i}\otimes {\mathit{I}}_{m_{i-1}}\otimes \cdots \otimes {\mathit{I}}_{m_j},$$

and let

$$\underset{p_{g,i}\times p_{g,i}}{{\mathit{H}}_{m_i,m_g}}={\mathit{H}}_{m_i}\otimes {\mathit{H}}_{m_{i-1}}\otimes \cdots \otimes {\mathit{H}}_{m_g},$$

where ${\mathit{H}}_{m_h}$ is an $m_h\times m_h$ Helmert matrix for each $h=2,\dots ,k,$ i.e., each ${\mathit{H}}_{m_h}$ is an orthogonal matrix whose first column is proportional to ${\mathbf{1}}_{m_h}.$ Then:

$${\mathit{L}}_h^{\prime }={\mathit{H}}_{m_h,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}$$

is an orthogonal matrix (note that ${\mathit{L}}_h$ are not function of either of the ${\mathit{U}}_i$’s), and in particular

$${\mathit{L}}_k^{\prime }={\mathit{H}}_{m_h,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}=\left({\mathit{I}}_{m_k}\otimes {\mathit{L}}_{k-1}^{\prime }\right)\left({\mathit{H}}_{m_k}^{\prime }\otimes {\mathit{I}}_{m_{k-1},m_1}\right).$$

(16)

Lemma 4 of Leiva and Roy [10] states and proves the block diagonalization result of the k-SSCS matrix ${\mathsf{\Gamma }}_k$ by using the orthogonal matrix ${\mathit{L}}_k^{\prime }$ as defined in (16), that is:

$$\begin{array}{ccc}\hfill {\mathit{L}}_k^{\prime }{\mathsf{\Gamma }}_k{\mathit{L}}_k& =& \mathrm{diag}\left\{{\mathit{D}}_{\mathit{f}};\mathit{f}={\left(f_k,f_{k-1},\dots ,f_2\right)}^{\prime }\in F\right\}\hfill \\ & =& \mathrm{diag}\left\{{\mathit{D}}_{f_k,f_{k-1},\dots ,f_2};{\left(f_k,f_{k-1},\dots ,f_2\right)}^{\prime }\in F_k\times F_{k-1}\times \dots \times F_2\right\},\hfill \end{array}$$

(17)

where, for each $j=1,\dots ,k,$ the $m_1\times m_1$-diagonal matrices ${\mathit{D}}_{\mathit{f}}={\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}$ are given by:

$$\begin{array}{ccc}{\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}={\mathsf{\Delta }}_{k,j}& if& f_2=1,\dots ,f_j=1,f_{j+1}\ne 1,\end{array}$$

where $f_{k+1}\ne 1$ is not taken into consideration, that is:

$$\begin{array}{ccc}\hfill {\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}& =& \left\{\begin{array}{ccc}{\mathsf{\Delta }}_{k,k}& \mathrm{if}& f_2=1,\dots ,f_{k-1}=1,f_k=1\\ {\mathsf{\Delta }}_{k,k-1}& \mathrm{if}& f_2=1,\dots ,f_{k-1}=1,f_k\ne 1\\ {\mathsf{\Delta }}_{k,k-2}& \mathrm{if}& f_2=1,\dots ,f_{k-2}=1,f_{k-1}\ne 1\\ ⋮& ⋮& ⋮\\ {\mathsf{\Delta }}_{k,2}& \mathrm{if}& f_2=1,f_3\ne 1\\ {\mathsf{\Delta }}_{k,1}& \mathrm{if}& f_2\ne 1\end{array}\right..\hfill \end{array}$$

(18)

Thus, ${\mathsf{\Delta }}_{k,j}forj=1,\dots ,k$ are the k eigenblocks of the k-SSCS covariance matrix ${\mathsf{\Gamma }}_k$. We will obtain the estimators of the eigenblocks ${\mathsf{\Delta }}_{k,j},j=1,\dots ,k$ in Section 4. In the following section, we briefly discuss that the k-SSCS covariance structure is of the Jordan algebra type.

k-SSCS Covariance Structure Is of the Jordan Algebra Type

The k-SSCS covariance structure is of the Jordan algebra type (Jordan et al. [23]). Let ${\mathcal{G}}_{p_{1,k}}$ be the set of all k-SSCS $p_{1,k}\times p_{1,k}$ matrices. It is clear that under the usual matrix addition and scalar multiplication, ${\mathcal{G}}_{p_{1,k}}$ is a subspace of the linear vectorial space ${\mathcal{S}}_{p_{1,k}}$ of the $p_{1,k}\times p_{1,k}$ symmetric matrices. For any natural number $k\ge 2$, it is easy to prove the following proposition:

Proposition 1.

If ${\mathsf{\Gamma }}_k$ is a k-SSCS matrix given by (9) (or equivlently by (8)), then${\mathsf{\Gamma }}_k^2={\mathsf{\Gamma }}_k{\mathsf{\Gamma }}_k$ is also a k-SSCS matrix given by:

$${\mathsf{\Gamma }}_k^2={\mathsf{\Gamma }}_k{\mathsf{\Gamma }}_k={\mathsf{\Gamma }}_k^{*}=\left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes {\mathit{T}}_{k,j}^{*}\right]+\underset{h=1}{\overset{k-1}{⨂}}{\mathit{J}}_{m_{k+1-h}}{\mathit{T}}_{k,k}^{*},$$

where

$$\begin{array}{ccc}\hfill {\mathit{T}}_{k,k}^{*}& =& \left[\sum _{j=1}^{k-1}{p}_{2,j}\left({\mathit{T}}_{k,j}{\mathit{T}}_{k,k}+{\mathit{T}}_{k,k}{\mathit{T}}_{k,j}\right)\right]+p_{2,k}{\mathit{T}}_{k,k}^2\hfill \\ & =:& {\mathit{U}}_{k,k}^{*},\hfill \end{array}$$

and for $h=k-1,k-2,\dots ,2,$ ${\mathit{T}}_{k,h}^{*}=:{\mathit{U}}_{k,h}^{*}-{\mathit{U}}_{k,h+1}^{*},$ with

$${\mathit{U}}_{k,h}^{*}=\left[\sum _{j=1}^{h-1}{p}_{2,j}\left({\mathit{T}}_{k,j}{\mathit{T}}_{k,h}+{\mathit{T}}_{k,h}{\mathit{T}}_{k,j}\right)\right]+p_{2,k}{\mathit{U}}_{k,h}^2+{\mathit{U}}_{k,h+1}^{*}$$

and where ${\mathit{T}}_{k,1}^{*}=:{\mathit{U}}_{k,1}^{*}-{\mathit{U}}_{k,2}^{*}={\left({\mathit{U}}_{k,1}-{\mathit{U}}_{k,2}\right)}^2={\mathit{T}}_{k,1}^2,$ that is, with ${\mathit{U}}_{k,1}^{*}=$${\mathit{T}}_{k,1}^2+{\mathit{U}}_{k,2}^{*}.$

Therefore, we conclude that ${\mathcal{G}}_{p_{1,k}}$ is a Jordan Algebra. See Lemma 4.1 on Page 10 in Malley [24], which states that ${\mathcal{G}}_{p_{1,k}}$ is a Jordan Algebra if and only if ${\mathsf{\Gamma }}_k^2={\mathsf{\Gamma }}_k{\mathsf{\Gamma }}_k\in {\mathcal{G}}_{p_{1,k}}$ for all ${\mathsf{\Gamma }}_k\in {\mathcal{G}}_{p_{1,k}}$. See Roy et al. [25] and Kozioł et al. [26] for proofs that 2-SSCS and 3-SSCS covariance structures are of Jordan algebra types.

4. Estimators of the Eigenblocks

Let ${\mathit{x}}_r:r=1,\dots ,n$ be $p_{1,k}\times 1$ random vectors partitioned into $m_1\times 1$ subvectors as follows:

$$\begin{array}{ccc}\hfill {\mathit{x}}_r& =& {\left({\mathit{x}}_{r;\mathit{f}}^{\prime }:\mathit{f}=\left(f_k,f_{k-1},\dots ,f_2\right)\in F=F_{k,2}=\stackrel{2}{\underset{j=k}{\times }}{F}_j\right)}^{\prime }\hfill \\ & =& {\left({\mathit{x}}_{r;f_k,f_{k-1},\dots ,f_2}^{\prime }:f_j\in F_j=\left\{1,\dots ,m_j\right\},\mathrm{for}j=k,k-1,\dots ,2\right)}^{\prime }.\hfill \end{array}$$

The vectors ${\mathit{x}}_r:r=1,\dots ,n$ are a random sample from a population with distribution ${\mathcal{N}}_{p_{1,k}}\left(\mathbf{\mu };{\mathsf{\Gamma }}_k\right),$ where ${\mathsf{\Gamma }}_k$ is a positive definite k-SSCS structured covariance matrix as given in (8) in Definition 3. Let $\mathit{X}$ be the $n\times p_{1,k}$-sample data matrix as follows:

$$\underset{n\times p_{1,k}}{\mathit{X}}=\left(\begin{array}{c}{\mathit{x}}_1^{\prime }\\ ⋮\\ {\mathit{x}}_n^{\prime }\end{array}\right)=\left(\underset{n\times m_1}{{\mathit{X}}_{·;1,\dots ,1}},\dots ,\underset{n\times m_1}{{\mathit{X}}_{·;f_k,\dots ,f_2}},\dots ,\underset{n\times m_1}{{\mathit{X}}_{·;m_k,\dots ,m_2}}\right)$$

with

$$\underset{n\times m_1}{{\mathit{X}}_{·,\mathit{f}}}=\underset{n\times m_1}{{\mathit{X}}_{·,f_k,f_{k-1},\dots ,f_2}=}\left(\begin{array}{c}{\mathit{x}}_{1,\mathit{f}}^{\prime }\\ ⋮\\ {\mathit{x}}_{n,\mathit{f}}^{\prime }\end{array}\right)={\left(\underset{1\times m_1}{{\mathit{x}}_{r,\mathit{f}}^{\prime }}\right)}_{r=1}^n.$$

In this section, we prove that certain unbiased estimators (to be defined) of the matrix parameters ${\mathit{U}}_{k,j}:j=1,\dots ,k-1$ can be written as functions of the usual sample variance-covariance matrix $\mathit{S}$ as follows:

$$\begin{array}{ccc}\hfill \mathit{S}& =& \frac{1}{n-1}{\mathit{X}}^{\prime }{\mathit{Q}}_n\mathit{X}\hfill \\ & =& \frac{1}{n-1}\left(\begin{array}{c}\underset{m_1\times n}{{\mathit{X}}_{·;1,\dots ,1}^{\prime }}\\ ⋮\\ \underset{m_1\times n}{{\mathit{X}}_{·;m_k,\dots ,m_2}^{\prime }}\end{array}\right){\mathit{Q}}_n\left(\underset{n\times m_1}{{\mathit{X}}_{·;1,\dots ,1}},\dots ,\underset{n\times m_1}{{\mathit{X}}_{·;m_k,\dots ,m_2}}\right)={\left({\mathit{S}}_{\mathit{f},{\mathit{f}}^{*}}\right)}_{\mathit{f},{\mathit{f}}^{*}\in \mathit{F}},\hfill \end{array}$$

where ${\mathit{Q}}_n$ is given in (2) with (3). Now the sample mean $\overline{\mathit{x}}$ can be expressed as:

$$\begin{array}{ccc}\hfill \underset{p_{1,k}\times 1}{\overline{\mathit{x}}}& =& \frac{1}{n}\underset{p_{1,k}\times n}{{\mathit{X}}^{\prime }}{\mathbf{1}}_n=\frac{1}{n}\left(\begin{array}{c}{\mathit{X}}_{·;1,\dots ,1}^{\prime }\\ ⋮\\ {\mathit{X}}_{·;m_k,\dots ,m_2}^{\prime }\end{array}\right){\mathbf{1}}_n=\left(\begin{array}{c}\frac{1}{n}{\mathit{X}}_{·;1,\dots ,1}^{\prime }{\mathbf{1}}_n\\ ⋮\\ \frac{1}{n}{\mathit{X}}_{·;m_k,\dots ,m_2}^{\prime }{\mathbf{1}}_n\end{array}\right)\hfill \\ & =& {\left(\frac{1}{n}\underset{m_1\times n}{{\mathit{X}}_{·;\mathit{f}}^{\prime }}{\mathbf{1}}_n\right)}_{\mathit{f}\in \mathit{F}}={\left(\frac{1}{n}\sum _{r=1}^n{\mathit{x}}_{r;\mathit{f}}\right)}_{\mathit{f}\in \mathit{F}}={\underset{p_{1,k}\times 1}{\left(\underset{m_1\times 1}{{\overline{\mathit{x}}}_{\mathit{f}}}\right)}}_{\mathit{f}\in \mathit{F}}.\hfill \end{array}$$

Thus, ${\mathit{S}}_{\mathit{f},{\mathit{f}}^{*}}$ in $\mathit{S}$ can be expressed as:

$$\begin{array}{ccc}\hfill {\mathit{S}}_{\mathit{f},{\mathit{f}}^{*}}& =& \frac{1}{n-1}{\mathit{X}}_{·;\mathit{f}}^{\prime }{\mathit{Q}}_n{\mathit{X}}_{·;{\mathit{f}}^{*}}^{\prime }=\left\{\begin{array}{ccc}\frac{1}{n-1}{\displaystyle \sum _{r=1}^n}\left(\underset{m_1\times 1}{{\mathit{x}}_{r;\mathit{f}}-{\overline{\mathit{x}}}_{\mathit{f}}}\right)\left(\underset{1\times m_1}{{\mathit{x}}_{r;\mathit{f}}^{\prime }-{\overline{\mathit{x}}}_{\mathit{f}}^{\prime }}\right)& \mathrm{if}& \mathit{f}={\mathit{f}}^{*}\\ \frac{1}{n-1}{\displaystyle \sum _{r=1}^n}\left(\underset{m_1\times 1}{{\mathit{x}}_{r;\mathit{f}}-{\overline{\mathit{x}}}_{\mathit{f}}}\right)\left(\underset{1\times m_1}{{\mathit{x}}_{r;{\mathit{f}}^{*}}^{\prime }-{\overline{\mathit{x}}}_{{\mathit{f}}^{*}}^{\prime }}\right)& \mathrm{if}& \mathit{f}\ne {\mathit{f}}^{*}\end{array}\right..\hfill \end{array}$$

Since $\mathit{S}$ is an unbiased estimator of ${\mathsf{\Gamma }}_k$, we have:

$$\begin{array}{ccc}\hfill \mathrm{E}\left[{\mathit{S}}_{\mathit{f};{\mathit{f}}^{*}}\right]& =& \mathrm{E}\left[{\mathit{S}}_{f_k,f_{k-1},\dots ,f_2;f_k^{*},f_{k-1}^{*},\dots ,f_2^{*}}\right]\hfill \\ & =& \left\{\begin{array}{ccccc}{\mathit{U}}_{k,1}& if& f_2=f_2^{*},\dots ,f_k=f_k^{*}& & \\ {\mathit{U}}_{k,j}& if& f_j\ne f_j^{*},f_{j+1}=f_{j+1}^{*},\dots ,f_k=f_k^{*}& \mathrm{for}& j=2,\dots ,k\end{array}\right..\hfill \end{array}$$

Therefore, to find a better unbiased estimator of ${\mathit{U}}_{k,j}$, we average all the above random matrices that are unbiased estimator of the same ${\mathit{U}}_{k,j}$. The unbiased estimators ${\widehat{\mathit{U}}}_{k,j}$ of ${\mathit{U}}_{k,j}$ for each $j=1,\dots ,k$ are derived in Lemma 5 in Leiva and Roy [10] with $q_{k,j}$ defined in Lemma 3 in Leiva and Roy [10] as:

$$q_{k,j}=\left\{\begin{array}{ccc}{p}_{2,k}& if& j=1\\ \begin{array}{c}{p}_{2,k}\left(m_j-1\right)p_{2,j-1}\end{array}& if& j=2,\dots ,k\end{array}\right..$$

(19)

Unbiased estimators of the eigenblocks ${\mathsf{\Delta }}_{k,j}$ can be obtained from (13). Then, using (14), the unbiased estimators of ${\mathsf{\Gamma }}_k$ can be obtained as the following othogonal sums:

$${\widehat{\mathsf{\Gamma }}}_k=\sum _{j=1}^k{\mathit{R}}_{k,j+1}^{*}\otimes {\widehat{\mathsf{\Delta }}}_{k,j},$$

and if ${\widehat{\mathsf{\Gamma }}}_k^{-1}$ exists, it can be obtained from (15) as follows:

$${\widehat{\mathsf{\Gamma }}}_k^{-1}=\sum _{j=1}^k{\mathit{R}}_{k,j+1}^{*}\otimes {\widehat{\mathsf{\Delta }}}_{k,j}^{-1},$$

(20)

where ${\mathit{R}}_{k,j+1}^{*}$ is given in (5), for each $j=1,\dots ,k-1,$ and, for $j=k,$ ${\mathit{R}}_{k,k+1}^{*}$ is given in (7).

The computation of the unbiased estimates of the component matrices ${\mathit{U}}_{k,j}$ for each $j=1,\dots ,k$ is easy, as all of them have explicit solutions. At this point, we want to mention that for k-separable covariance structure the estimates of the component matrices are not easy, as the MLEs have implicit equations, and therefore are not tractable analytically. Now, from Theorem 1 of Leiva and Roy [10], we see that a multiple of the unbiased estimators of the eigenblocks ${\widehat{\mathsf{\Delta }}}_{k,j}$ for each $j=1,\dots ,k$, have Wishart distributions as follows:

$$\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right){\widehat{\mathsf{\Delta }}}_{k,j}=\left(n-1\right){\mathrm{BT}}_{p_{1,1}}\left({\mathit{R}}_{k,j+1}\mathit{S}{\mathit{R}}_{k,j+1}\right):j=1,\dots ,k-1,$$

$$\mathrm{and}\left(n-1\right){\widehat{\mathsf{\Delta }}}_{k,k}=\left(n-1\right){\mathrm{BT}}_{p_{1,1}}\left({\mathit{R}}_{k,k+1}\mathit{S}{\mathit{R}}_{k,k+1}\right),$$

where ${\mathit{R}}_{k,j+1}={\mathit{R}}_{k,j+1}^{*}\otimes {\mathit{I}}_{m_1}$ given by (4) with (5) and (6) with (7), are independent and

$$\begin{array}{ccc}\hfill \left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right){\widehat{\mathsf{\Delta }}}_{k,j}& \sim & {\mathcal{W}}_{m_1}\left(\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right);{\mathsf{\Delta }}_{k,j}\right)\mathrm{for}1,\dots ,k-1\hfill \\ \hfill \mathrm{and}\left(n-1\right){\widehat{\mathsf{\Delta }}}_{k,k}& \sim & {\mathcal{W}}_{m_1}\left(\left(n-1\right);{\mathsf{\Delta }}_{k,k}\right).\hfill \end{array}$$

From Corollary 1 of Leiva and Roy [10], the 2-SSCS covariance matrix for second-order data or multivariate repeated measures data has two eigenblocks, ${\mathsf{\Delta }}_{2,2}$ and ${\mathsf{\Delta }}_{2,1}$ with multiplicity $m_2-1$, and their distributions are as follows:

$$\begin{array}{ccc}\hfill \left(n-1\right)\left(m_2-1\right){\widehat{\mathsf{\Delta }}}_{2,1}& \sim & {\mathcal{W}}_{m_1}\left((n-1)(m_2-1),{\mathsf{\Delta }}_{2,1}\right),\hfill \\ \hfill \mathrm{and}\left(n-1\right){\widehat{\mathsf{\Delta }}}_{2,2}& \sim & {\mathcal{W}}_{m_1}\left((n-1),{\mathsf{\Delta }}_{2,2}\right).\hfill \end{array}$$

The 3-SSCS covariance matrix for third-order data has three eigenblocks, ${\mathsf{\Delta }}_{3,3}$, ${\mathsf{\Delta }}_{3,2}$ with multiplicity $m_3-1$ and ${\mathsf{\Delta }}_{3,1}$ with multiplicity $m_3(m_2-1)$, and their distributions are as follows:

$$\begin{array}{ccc}\hfill \left(n-1\right)m_3\left(m_2-1\right){\widehat{\mathsf{\Delta }}}_{3,1}& \sim & {\mathcal{W}}_{m_1}\left((n-1)m_3(m_2-1),{\mathsf{\Delta }}_{3,1}\right),\hfill \\ \hfill and\left(n-1\right)\left(m_3-1\right){\widehat{\mathsf{\Delta }}}_{3,2}& \sim & {\mathcal{W}}_{m_1}\left((n-1)(m_3-1),{\mathsf{\Delta }}_{3,2}\right),\hfill \\ \hfill and\left(n-1\right){\widehat{\mathsf{\Delta }}}_{3,3}& \sim & {\mathcal{W}}_{m_1}\left((n-1),{\mathsf{\Delta }}_{3,3}\right).\hfill \end{array}$$

5. Test for the Mean

5.1. One Sample Test

Using the notation and assumptions in Section 4, let $\underset{n\times p_{1,k}}{\mathit{X}}={\left({\mathit{x}}_1,\dots ,{\mathit{x}}_n\right)}^{\prime }$ be a $n\times p_{1,k}$-dimensional data matrix formed from the random samples ${\mathit{x}}_1,\dots ,{\mathit{x}}_n$ from ${\mathcal{N}}_{p_{1,k}}\left(\mathbf{\mu };{\mathsf{\Gamma }}_k\right)$. Let $\overline{\mathit{x}}$ be the sample mean, then $\overline{\mathit{x}}\sim {\mathcal{N}}_{p_{1,k}}\left(\mathbf{\mu },\frac{1}{n}{\mathsf{\Gamma }}_k\right)$. We are interested in testing the following hypothesis:

$${\mathcal{H}}_0:\mathbf{\mu }={\mathbf{\mu }}_{0}vs{\mathcal{H}}_1:\mathbf{\mu }\ne {\mathbf{\mu }}_0,$$

(21)

for known ${\mathbf{\mu }}_0$. For testing hypothesis (21), we use the test statistic $D^2$ defined as:

$$D^2=n{\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)}^{\prime }{\widehat{\mathsf{\Gamma }}}_k^{-1}\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right).$$

(22)

5.1.1. Distribution of Test Statistic $D^2$ under $H_0$

Now, let ${\mathit{L}}_k^{\prime }={\mathit{H}}_{m_k,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}$ be the matrix as given in (16). We use here the following canonical transformation:

$$\mathit{z}={\mathit{L}}_k^{\prime }\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)=\left({\mathit{H}}_{m_k,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right)\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right).$$

Therefore, $\mathit{z}={\left(\mathit{z}{{}_{1,\dots ,1}}^{\prime },\cdots ,\mathit{z}{{}_{m_k,\dots ,m_2}}^{\prime }\right)}^{\prime }={\mathit{L}}_k^{\prime }\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)\sim {\mathcal{N}}_{p_{1,k}}\left(\mathbf{0},{\mathsf{\Omega }}_k\right)$. Therefore, according to (17) with (18), we have:

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_k& =& \frac{1}{n}{\mathit{D}}_k={\mathit{L}}_k^{\prime }\left(\frac{1}{n}{\mathsf{\Gamma }}_k\right){\mathit{L}}_k\hfill \\ & =& \frac{1}{n}\left({\mathit{H}}_{m_k,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right){\mathsf{\Gamma }}_k\left({\mathit{H}}_{m_k,m_2}\otimes {\mathit{I}}_{m_1}\right)\hfill \\ & =& \frac{1}{n}\mathrm{diag}\left\{{\mathit{D}}_{f_k,f_{k-1},\dots ,f_2};{\left(f_k,f_{k-1},\dots ,f_2\right)}^{\prime }\in F_k\times F_{k-1}\times \dots \times F_2\right\},\hfill \end{array}$$

where, for each $j=1,\dots ,k,$ the diagonal $m_1\times m_1$- matrices ${\mathit{D}}_{\mathit{f}}={\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}$ are given by:

$$\begin{array}{ccc}{\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}={\mathsf{\Delta }}_{k,j}& if& f_2=1,\dots ,f_j=1,f_{j+1}\ne 1,\end{array}$$

where $f_{k+1}\ne 1$ is not taken into consideration, and the $m_1\times 1$ component vectors ${\mathit{z}}_{f_k,\dots ,f_2},$ with $\mathit{f}=\left(f_k,\dots ,f_2\right)\in F,$ are independent. The distribution of ${\mathit{z}}_{f_k,\dots ,f_2}$, under ${\mathcal{H}}_0$ is given by:

$$\begin{array}{c}\hfill {\mathit{z}}_{f_k,\dots ,f_2}\sim {\mathcal{N}}_{m_1}\left(\mathbf{0},\frac{1}{n}{\mathsf{\Delta }}_{k,j}\right)if\end{array}$$

$${\mathit{f}}_{2,j}=\left(f_2,\dots ,f_j\right),f_{j+1}\in F_{j+1}-\left\{1\right\},\mathrm{and}{\mathit{f}}_{j+2,k}=\left(f_{j+2},\dots ,f_k\right)\in F_{j+2}\times \cdots \times F_k.$$

Since

$$\begin{array}{ccc}\hfill {\mathit{H}}_{m_i}{\mathit{I}}_{m_i}{\mathit{H}}_{m_i}^{\prime }& =& {\mathit{I}}_{m_i}=\left[{\mathit{e}}_{m_i,1},{\mathit{e}}_{m_i,2},\dots ,{\mathit{e}}_{m_i,m_i}\right]\hfill \\ \hfill {\mathit{H}}_{m_i}{\mathit{P}}_{m_i}{\mathit{H}}_{m_i}^{\prime }& =& {\mathit{e}}_{m_i,1}{\mathit{e}}_{m_i,1}^{\prime }=\mathrm{diag}\left(1,\underset{m_i-1}{\underbrace{0,\dots ,0}}\right)\hfill \\ \hfill \mathrm{and}{\mathit{H}}_{m_i}{\mathit{Q}}_{m_i}{\mathit{H}}_{m_i}^{\prime }& =& \sum _{j=2}^{m_i}{\mathit{e}}_{m_i,j}{\mathit{e}}_{m_i,j}^{\prime }=\mathrm{diag}\left(0,\underset{m_i-1}{\underbrace{1,\dots ,1}}\right),\hfill \end{array}$$

for $j=1,\dots ,k-1,$ we have:

$$\begin{array}{ccc}\hfill {\mathit{H}}_{m_k,m_2}{\mathit{R}}_{k,j+1}^{*}{\mathit{H}}_{m_k,m_2}^{\prime }& =& \left({\mathit{H}}_{m_k,m_{j+2}}\otimes {\mathit{H}}_{m_{j+1}}\otimes {\mathit{H}}_{m_j,m_2}\right)\left({\mathit{I}}_{p_{j+2,k}}\otimes {\mathit{Q}}_{m_{j+1}}\otimes {\mathit{P}}_{m_{j,}{m}_2}\right)\hfill \\ & & \left({\mathit{H}}_{m_k,m_{j+2}}^{\prime }\otimes {\mathit{H}}_{m_{j+1}}^{\prime }\otimes {\mathit{H}}_{m_j,m_2}^{\prime }\right)\hfill \\ & =& \underset{i=1}{\overset{k-j-1}{⨂}}\left[{\mathit{e}}_{m_{k+1-i},1},\dots ,{\mathit{e}}_{m_{k+1-i},m_{k+1-i}}\right]\otimes \left(\sum _{i=2}^{m_{j+1}}{\mathit{e}}_{m_{j+1},i}{\mathit{e}}_{m_{j+1},i}^{\prime }\right)\otimes \left({\mathit{e}}_{p_{j,2},1}{\mathit{e}}_{p_{j,2},1}^{\prime }\right)\hfill \\ & =& \mathrm{diag}\left(\mathbf{⋮}\underset{p_{j+2,k}}{\underbrace{\underset{m_{j+1}}{\underbrace{\underset{p_{2,j}}{\underbrace{0,\dots ,0}⋮}1,\underset{p_{2,j}-1}{\underbrace{0,\dots ,0}⋮}\cdots ⋮1,\underset{p_{2,j}-1}{\underbrace{0,\dots ,0}}}}\mathbf{⋮}\cdots \mathbf{⋮}\underset{m_{j+1}}{\underbrace{\underset{p_{2,j}}{\underbrace{0,\dots ,0}⋮}1,\underset{p_{2,j}-1}{\underbrace{0,\dots ,0}⋮}\cdots ⋮1,\underset{p_{2,j}-1}{\underbrace{0,\dots ,0}}}}}}\mathbf{⋮}\right),\hfill \end{array}$$

and for $j=k,$ we have:

$$\begin{array}{ccc}\hfill {\mathit{H}}_{m_k,m_2}{\mathit{R}}_{k,k+1}^{*}{\mathit{H}}_{m_k,m_2}^{\prime }& =& {\mathit{H}}_{m_k,m_2}{\mathit{P}}_{m_{k,}{m}_2}{\mathit{H}}_{m_k,m_2}^{\prime }={\mathit{e}}_{p_{2,k},1}{\mathit{e}}_{p_{2,k},1}^{\prime }=\mathrm{diag}\left(1,\underset{p_{2,k}-1}{\underbrace{0,\dots ,0}}\right).\hfill \end{array}$$

Therefore, using (20), the statistic $D^2$ in (22) can be written as:

$$\begin{array}{ccc}\hfill D^2& =& n{\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)}^{\prime }{\widehat{\mathsf{\Gamma }}}_k^{-1}\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)\hfill \\ & =& n{\mathit{z}}^{\prime }\left({\mathit{H}}_{m_k,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right)\left\{\sum _{j=1}^k{\mathit{R}}_{k,j+1}^{*}\otimes {\widehat{\mathsf{\Delta }}}_{k,j}^{-1}\right\}\left({\mathit{H}}_{m_k,m_2}\otimes {\mathit{I}}_{m_1}\right)\mathit{z}\hfill \\ & =& n{\mathit{z}}^{\prime }\sum _{j=1}^{k-1}\left\{\mathrm{diag}\left(\stackrel{p_{j+2,k}}{\overbrace{\mathbf{⋮}\underset{m_{j+1}}{\underbrace{⋮{\mathbf{0}}_{p_{2,j}}^{\prime }⋮1,{\mathbf{0}}_{p_{2,j-1}}^{\prime }⋮\cdots ⋮1,{\mathbf{0}}_{p_{2,j-1}}^{\prime }⋮}}\mathbf{⋮}\cdots \mathbf{⋮}\underset{m_{j+1}}{\underbrace{⋮{\mathbf{0}}_{p_{2,j}}^{\prime }⋮1,{\mathbf{0}}_{p_{2,j-1}}^{\prime }⋮\cdots ⋮1,{\mathbf{0}}_{p_{2,j-1}}^{\prime }⋮}}\mathbf{⋮}}}\right)\right.\hfill \\ & & \left.\underset{}{\stackrel{}{\stackrel{}{\underset{}{}}}}\otimes {\widehat{\mathsf{\Delta }}}_{k,j}^{-1}\right\}\mathit{z}+n{\mathit{z}}^{\prime }\left\{\mathrm{diag}\left(1,\underset{p_{2,k}-1}{\underbrace{0,\dots ,0}}\right)\otimes {\widehat{\mathsf{\Delta }}}_{k,k}^{-1}\right\}\mathit{z},\hfill \end{array}$$

that is,

$$\begin{array}{c}\hfill D^2=:\sum _{j=1}^k{T}_{0j}^2,\end{array}$$

(23)

where for $j=1,\dots ,k-1$

$$\begin{array}{c}\hfill T_{0j}^2=tr\left(n\underset{k-jSums}{\underbrace{\sum _{f_k=1}^{m_k}\dots \sum _{f_{j+2}=1}^{m_{j+2}}\sum _{f_{j+1}=2}^{m_{j+1}}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\widehat{\mathsf{\Delta }}}_{k,j}^{-1}\right)\end{array}$$

(24)

and for $j=k$, we assume

$$\begin{array}{ccc}\hfill T_{0k}^2& =& tr\left[n\underset{k-jSums}{\underbrace{\sum _{f_k=1}^{m_k}\dots \sum _{f_{j+2}=1}^{m_{j+2}}\sum _{f_{j+1}=2}^{m_{j+1}}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\widehat{\mathsf{\Delta }}}_{k,k}^{-1}\right]\hfill \\ & =& ntr\left(\mathit{z}{}_{\underset{k-1ones}{\underbrace{1,\dots ,1}}}\mathit{z}{}_{\underset{k-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\widehat{\mathsf{\Delta }}}_{k,k}^{-1}\right).\hfill \end{array}$$

(25)

Note that the subsets of vectors involved in $T_{01}^2,\dots ,T_{0k}^2,$ respectively, form a partition of the set of independent vectors $\left\{\mathit{z}{}_{f_k,\dots ,f_2}:f_k,\dots ,f_2\in F\right\}$. Therefore, $T_{01}^2,\dots ,T_{0k}^2$ are mutually independent. Moreover, since for $j=1,\dots ,k,$

$${\mathit{H}}_j=n\underset{k-jSums}{\underbrace{\sum _{f_k=1}^{m_k}\dots \sum _{f_{j+2}=1}^{m_{j+2}}\sum _{f_{j+1}=2}^{m_{j+1}}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}^{\prime }\sim {\mathcal{W}}_{m_1}\left(k_{\left(j\right)};{\mathsf{\Delta }}_{k,j}\right),$$

where

$$k_{\left(j\right)}=p_{j+2,k}\left(m_{j+1}-1\right),$$

and

$${\mathit{E}}_j=\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right){\widehat{\mathsf{\Delta }}}_{k,j}\sim {\mathcal{W}}_{m_1}\left(d_j;{\mathsf{\Delta }}_{k,j}\right),$$

with

$$d_j=\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right).$$

Therefore, $T_{0j}^2$ given by (24) reduces to

$$\begin{array}{ccc}\hfill T_{0j}^2& =& tr\left[{\mathit{H}}_j{\left(\frac{{\mathit{E}}_j}{\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)}\right)}^{-1}\right]\hfill \\ & =& \left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)tr\left({\mathit{H}}_j{\mathit{E}}_j^{-1}\right)=d_{j}tr\left({\mathit{H}}_j{\mathit{E}}_j^{-1}\right),\hfill \end{array}$$

and has a Lawley–Hotelling trace (LH-trace) distribution denoted by ${\mathcal{T}}_0^2\left(m_1;p_{j+2,k}\left(m_{j+1}-1\right),d_j\right)$ if $d_j=\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)\ge m_1.$ Note that, using (25), the case $j=k$ reduces to ${\mathit{H}}_k=n\mathit{z}{}_{1,\dots ,1}\mathit{z}{{}_{1,\dots ,1}}^{\prime }\sim {\mathcal{W}}_{m_1}\left(1,{\mathsf{\Delta }}_{k,k}\right)$ and ${\mathit{E}}_k=\left(n-1\right){\widehat{\mathsf{\Delta }}}_{k,k}\sim {\mathcal{W}}_{m_1}\left(n-1;{\mathsf{\Delta }}_{k,k}\right)$. Then, $T_{0k}^2$ has the LH-trace distribution ${\mathcal{T}}_0^2\left(m_1;1,\left(n-1\right)\right)$ if $\left(n-1\right)\ge m_1.$

Thus, the distribution of $D^2$ given by (23) is the convolution of k-independent LH-trace distributions:

$$\underset{j=1}{\overset{k}{⨁}}{\mathcal{T}}_0^2\left(m_1;p_{j+2,k}\left(m_{j+1}-1\right),d_j=\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)\ge m_1\right).$$

The critical values of this distribution can be obtained using simulations. However, LH-trace distribution is usually aproximated by $\mathcal{F}$ distribution, and we use here the second approximation suggested in McKeon [27]. For jth case, i.e., for $j=1,\dots ,k-1,$ let us use the notations $m_{\left(j\right)}=d_j-m_1-1=\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)-m_1-1$, $k_{\left(j\right)}=p_{j+2,k}\left(m_{j+1}-1\right)$ and

$$\begin{array}{ccc}\hfill B_j& =& \frac{\left[m_{\left(j\right)}+k_{\left(j\right)}\right]}{\left[m_{\left(j\right)}-2\right]}·\frac{\left[m_{\left(j\right)}+m_1\right]}{\left[m_{\left(j\right)}+1\right]}\hfill \\ & =& \frac{\left[n{p}_{j+2,k}\left(m_{j+1}-1\right)-m_1-1\right]\left[\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)-1\right]}{\left[\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)-m_1-3\right]\left[\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)-m_1\right]}.\hfill \end{array}$$

Then, the distribution

$${\mathcal{T}}_0^2\left(m_1;p_{j+2,k}(m_{j+1}-1),\left(n-1\right)p_{j+2,k}(m_{j+1}-1)\right)$$

of $T_{0j}^2=d_{j}tr\left({\mathit{H}}_j{\mathit{E}}_j^{-1}\right)$ can be approximated by $g_j\mathcal{F}\left(K_j,D_j\right)$ where $K_j=k_{\left(j\right)}{m}_1=p_{j+2,k}$$\left(m_{j+1}-1\right)$$m_1,$$D_j=4+\frac{K_j+2}{B_j-1}=4+\frac{p_{j+2,k}\left(m_{j+1}-1\right)m_1+2}{B_j-1}$ and $g_j=\frac{d_j{K}_j}{m_{\left(j\right)}}\frac{D_j-2}{D_j}=\frac{\left(n-1\right)p_{j+2,k}^2{\left(m_{j+1}-1\right)}^2{m}_1}{\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right)-m_1-1}\frac{D_j-2}{D_j}.$

Finally, for $j=k,$ the distribution

$${\mathcal{T}}_0^2\left(m_1;1,(n-1)\right)$$

is the usual Hotelling ${\mathcal{T}}_{1,n-1}^2,$ that is, distributed as an exact distribution as follows:

$$\frac{\left(n-1\right)m_1}{n-m_1}\mathcal{F}\left(m_1,n-m_1\right).$$

This means that the distribution of $D^2$ can be approximated by the convolution of the above k distributions ($\left(k-1\right)$ approximated $\mathcal{F}$ distribution and one exact $\mathcal{F}$ distribution), where its critical values are obtained by the method suggested by Dyer [28].

Remark 1.

The statistic $T_{0j}^2,j=1,\dots ,k-1$ has LH-trace distribution ${\mathcal{T}}_0^2\left(m_1;p_{j+2,k}(m_{j+1}-1),d_j\right)$ if $d_j=(n-1)p_{j+2,k}(m_{j+1}-1)\ge m_1$. We also note that $T_{0k}^2$ has the LH-trace distribution ${\mathcal{T}}_0^2\left(m_1;1,\left(n-1\right)\right)$ if $\left(n-1\right)\ge m_1$. Now, for k-th order data, all $m_j\ge 2$ for $j=1,\dots ,k$ and $k\ge 2$. See Definition 1. Now, $k>j$. Therefore, $k-j\ge 1$ and then $k-j-2\ge -1$. Thus, $p_{j+2,k}\ge 1$. Since $m_j\ge 2$ for all k-th order data, $m_{j+1}\ge 2$, we have $(n-1)p_{j+2,k}(m_{j+1}-1)\ge m_1$ when $\left(n-1\right)\ge m_1$. Therefore, the only constraint needed on sample size in order to have all $T_{0j}^2,j=1,\dots ,k$ to follow LH-trace distribution is $\left(n-1\right)\ge m_1$, i.e., $n\ge m_1+1$, regardless of any $m_j,j=2,\dots ,k$. In essence, the minimum sample size needed to compute the $D^2$ test statistic is $m_1+1$, although the minimum sample size needed to compute the Hotelling’s ${\mathcal{T}}^2$ test statistic is $p_{1,k}+1$, where $p_{1,k}$ is the full dimension of the observations. For this reason, one cannot compute Hotelling’s ${\mathcal{T}}^2$ test statistic for a small sample data set where $n\le p_{1,k}$, which is doable for the $D^2$ test statistic.

We will now discuss some special cases of the $D^2$ statistic in the following remark.

Remark 2.

For second-order data or multivariate repeated measures data, $D^2=T_{01}^2+T_{02}^2$. Now, $T_{02}^2$ is distributed as LH-trace distribution ${\mathcal{T}}_0^2(m_1;1;n-1)$, and $T_{01}^2$ is distributed as LH-trace distribution as follows

$$\begin{array}{ccc}& & {\mathcal{T}}_0^2\left(m_1;p_{1+2,k}(m_{1+1}-1);(n-1)p_{1+2,k}(m_{1+1}-1)\right)\hfill \\ \hfill or& & {\mathcal{T}}_0^2\left(m_1;p_{3,2}(m_2-1);(n-1)p_{3,2}(m_2-1)\right)ask=2\hfill \\ \hfill or& & {\mathcal{T}}_0^2\left(m_1;(m_2-1);(n-1)(m_2-1)\right)as{p}_{3,2}=1by\left(1\right).\hfill \end{array}$$

Thus, $T_{01}^2$ is distributed as LH-trace distribution ${\mathcal{T}}_0^2\left(m_1;(m_2-1);(n-1)(m_2-1)\right)$.

So, we see that this test exactly matches the test obtained by Žežula et al. [1] for multivariate repeated measures data (second-order data) with 2-SSCS or BCS covariance structure. Therefore, we can say that our mean test statistic in this article is an extension or generalization of Žežula et al.’s [1] mean test statistic for k-th order data with k-SSCS covariance structure.

We will now derive the mean test statistic for third-order data with 3-SSCS covariance structure. For third-order data, $D^2=T_{01}^2+T_{02}^2+T_{03}^2$. Now, $T_{03}^2$ is distributed as LH-trace distribution ${\mathcal{T}}_0^2(m_1;1;n-1)$, and $T_{02}^2$ is distributed as LH-trace distribution as follows:

$$\begin{array}{ccc}& & {\mathcal{T}}_0^2\left(m_1;p_{2+2,k}(m_{2+1}-1);(n-1)p_{2+2,k}(m_{2+1}-1)\right)\hfill \\ \hfill or& & {\mathcal{T}}_0^2\left(m_1;p_{4,3}(m_3-1);(n-1)p_{4,3}(m_3-1)\right)ask=3\hfill \\ \hfill or& & {\mathcal{T}}_0^2\left(m_1;(m_3-1);(n-1)(m_3-1)\right)as{p}_{4,3}=1by\left(1\right),\hfill \end{array}$$

that can be approximated by $g_2\mathcal{F}\left(K_2,D_2\right),$ where $K_2=\left(m_3-1\right)m_1,$ $D_2=4+\frac{K_2+2}{B_2-1}$ and $g_2=\frac{d_2{K}_2}{m_{\left(2\right)}}\frac{D_2-2}{D_2}=\frac{\left(n-1\right){\left(m_3-1\right)}^2{m}_1}{\left(n-1\right)\left(m_3-1\right)-m_1-1}\frac{D_2-2}{D_2},$ with $d_2=\left(n-1\right)\left(m_3-1\right),$ $m_{\left(2\right)}=d_2-m_1-1$ and $B_2=\frac{\left[m_{\left(2\right)}+\left(m_3-1\right)m_1\right]\left[m_{\left(2\right)}+m_1\right]}{\left[m_{\left(2\right)}-2\right]\left[m_{\left(2\right)}+1\right]},$ and $T_{01}^2$ is distributed as LH-trace distribution as follows

$$\begin{array}{ccc}& & {\mathcal{T}}_0^2\left(m_1;p_{1+2,k}(m_{1+1}-1);(n-1)p_{1+2,k}(m_{1+1}-1)\right)\hfill \\ \hfill or& & {\mathcal{T}}_0^2\left(m_1;p_{3,3}(m_2-1);(n-1)p_{3,3}(m_2-1)\right)ask=3\hfill \\ \hfill or& & {\mathcal{T}}_0^2\left(m_1;m_3(m_2-1);(n-1)m_3(m_2-1)\right)as{p}_{3,3}=m_{3}by\left(1\right),\hfill \end{array}$$

that can be approximated by$g_1\mathcal{F}\left(K_1,D_1\right)$where$K_1=m_3\left(m_2-1\right)m_1,$$D_1=4+\frac{K_1+2}{B_1-1}$ and $g_1=\frac{d_1{K}_1}{m_{\left(1\right)}}\frac{D_1-2}{D_1}=\frac{\left(n-1\right)m_3^2{\left(m_2-1\right)}^2{m}_1}{\left(n-1\right)m_3\left(m_2-1\right)-m_1-1}\frac{D_1-2}{D_1},$ with $d_1=\left(n-1\right)m_3\left(m_2-1\right),$ $m_{\left(1\right)}=d_1-m_1-1$ and $B_1=\frac{\left[m_{\left(1\right)}+m_3\left(m_2-1\right)m_1\right]\left[m_{\left(1\right)}+m_1\right]}{\left[m_{\left(1\right)}-2\right]\left[m_{\left(1\right)}+1\right]}.$

So, one can easily derive the test statistic for j-th order data for $j=2,\dots ,k$ from our generalized $D^2$ statistic. The distribution of $D^2$ under $H_0$ for second-order data with 2-SSCS covariance structure is discussed in detail in Žežula et al. [1]. We will discuss the distribution of $D^2$ under $H_0$ for third-order data in detail in the following section.

5.1.2. Distribution of Statistic $D^2$ under $H_0$ for Third-order Data with 3-SSCS Covariance Structure

This section is adopted from the work of Žežula et al. [22]. However, we use a much simpler, straightforward approach so that the practitioners or the analysts can appreciate and apply the method easily. Let ${\mathit{L}}_3={\mathit{H}}_{m_3}\otimes {\mathit{H}}_{m_2}\otimes {\mathit{I}}_{m_1}={\mathit{H}}_{m_3,m_2}\otimes {\mathit{I}}_{m_1}$ be a matrix such that for each $j=2,3$${\mathit{H}}_{m_j}$ is an $(m_j\times m_j)$ Helmert matrix, that is, an $(m_j\times m_j)$ orthogonal matrix with the first column proportional to the $m_j\times 1$ vector of 1’s. We use here the following canonical transformation:

$$\mathit{z}={\mathit{L}}_3^{\prime }\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)=\left({\mathit{H}}_{m_3,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right)\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right).$$

Therefore, $\mathit{z}={\left(\mathit{z}{{}_{1,1}}^{\prime },\cdots ,\mathit{z}{{}_{m_3,m_2}}^{\prime }\right)}^{\prime }={\mathit{L}}_3^{\prime }\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)\sim {\mathcal{N}}_{p_{1,k}}\left(\mathbf{0},{\mathsf{\Omega }}_3\right),$ where

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_3& =& \frac{1}{n}{\mathit{D}}_3={\mathit{L}}_3^{\prime }\left(\frac{1}{n}{\mathsf{\Gamma }}_3\right){\mathit{L}}_3=\frac{1}{n}\left({\mathit{H}}_{m_3,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right){\mathsf{\Gamma }}_3\left({\mathit{H}}_{m_3,m_2}\otimes {\mathit{I}}_{m_1}\right)\hfill \\ & =& \frac{1}{n}\mathrm{diag}\left\{D_{f_3{f}_2};{\left(f_3,f_2\right)}^{\prime }\in \mathit{F}=F_3\times F_2\right\},\hfill \end{array}$$

where, for each $j=1,2,3,$ the diagonal $m_1\times m_1$-matrices ${\mathit{D}}_{\mathit{f}}={\mathit{D}}_{f_3,f_2}$ are given by

$$\begin{array}{ccc}{\mathit{D}}_{f_3,f_2}={\mathsf{\Delta }}_{3,j}& \mathrm{if}& \left\{\begin{array}{cc}{f}_2\ne 1& j=1\\ f_2=1,f_3\ne 1& j=2\\ f_2=1,f_3=1& j=3\end{array}\right.,\end{array}$$

and the $m_1\times 1$ component vectors ${\mathit{z}}_{f_3,f_2},$ with $\mathit{f}=\left(f_3,f_2\right)\in \mathit{F},$ are independent, with distributions (under the null hypothesis) ${\mathit{z}}_{f_3,f_2}\sim {\mathcal{N}}_{m_1}\left(\mathbf{0},\frac{1}{n}{\mathsf{\Delta }}_{3,j}\right)$ if

$$\begin{array}{c}\hfill {\mathit{f}}_{2,j}=\left(f_2,\dots ,f_j\right)={\mathbf{1}}_{j-1},f_{j+1}\in F_{j+1}-\left\{1\right\},\mathrm{and}{\mathit{f}}_{j+2,3}=\left(f_{j+2},\dots ,f_k\right)\in F_{j+2}\times \cdots \times F_3.\end{array}$$

Therefore, particularizing $D^2$, given by (23), for $k=3$ we have

$$\begin{array}{ccc}\hfill D^2& =& \sum _{j=1}^k\left[tr\left\{n\underset{k-jSums}{\underbrace{\sum _{f_k=1}^{m_k}\dots \sum _{f_{j+2}=1}^{m_{j+2}}\sum _{f_{j+1}=2}^{m_{j+1}}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}\mathit{z}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\widehat{\mathsf{\Delta }}}_{k,j}^{-1}\right\}\right]\hfill \\ & =& tr\left[n\underset{3-1Sums}{\underbrace{\sum _{f_k=1}^{m_3}\sum _{f_2=2}^{m_2}}}\mathit{z}{}_{f_3,f_2}\mathit{z}{}_{f_3,f_2}^{\prime }{\widehat{\mathsf{\Delta }}}_{3,1}^{-1}\right]+tr\left[n\underset{3-2Sums}{\underbrace{\sum _{f_k=2}^{m_3}}}\mathit{z}{}_{f_3,1}\mathit{z}{}_{f_3,1}^{\prime }{\widehat{\mathsf{\Delta }}}_{3,2}^{-1}\right]+tr\left[n\mathit{z}{}_{1,1}\mathit{z}{}_{1,1}^{\prime }{\widehat{\mathsf{\Delta }}}_{3,3}^{-1}\right]\hfill \\ & =:& T_{01}^2+T_{02}^2+T_{03}^2.\hfill \end{array}$$

Since the subsets of vectors involved in $T_{01}^2,T_{02}^2,T_{03}^2,$ respectively, form a partition of the set of independent vectors $\left\{\mathit{z}{}_{f_3,f_2}:f_3,f_2\in F_3\times F_2=F\right\},$$T_{01}^2,T_{02}^2,T_{03}^2$ are mutually independent. Moreover, since, for $j=1,$

$$\begin{array}{ccc}\hfill {\mathit{H}}_1& =& n\underset{3-1Sums}{\underbrace{\sum _{f_k=1}^{m_3}\sum _{f_2=2}^{m_2}}}\mathit{z}{}_{f_3,f_2}\mathit{z}{}_{f_3,f_2}^{\prime }\sim {\mathcal{W}}_{m_1}\left(m_3\left(m_2-1\right);{\mathsf{\Delta }}_{k,j}\right)={\mathcal{W}}_{m_1}\left(k_1;{\mathsf{\Delta }}_{k,j}\right),\mathrm{and}\hfill \\ \hfill {\mathit{E}}_1& =& \left(n-1\right)m_3\left(m_2-1\right){\widehat{\mathsf{\Delta }}}_{3,1}\sim {\mathcal{W}}_{m_1}\left(\left(n-1\right)m_3\left(m_2-1\right);{\mathsf{\Delta }}_{3,1}\right)={\mathcal{W}}_{m_1}\left(d_1;{\mathsf{\Delta }}_{3,1}\right),\mathrm{and}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill T_{01}^2& =& \mathrm{tr}\left[n\sum _{f_k=1}^{m_3}\sum _{f_2=2}^{m_2}\mathit{z}{}_{f_3,f_2}\mathit{z}{}_{f_3,f_2}^{\prime }{\widehat{\mathsf{\Delta }}}_{3,1}^{-1}\right]=\mathrm{tr}\left[{\mathit{H}}_1{\left(\frac{{\mathit{E}}_1}{\left(n-1\right)m_3\left(m_2-1\right)}\right)}^{-1}\right]\hfill \\ & =& \left(n-1\right)m_3\left(m_2-1\right)\mathrm{tr}\left({\mathit{H}}_1{\mathit{E}}_1^{-1}\right)=d_1\mathrm{tr}\left({\mathit{H}}_1{\mathit{E}}_1^{-1}\right)\hfill \end{array}$$

has a LH-trace distribution denoted by ${\mathcal{T}}_0^2\left(m_1;m_3\left(m_2-1\right),d_1\right)$ if $d_1=(n-1)m_3(m_2-1)>m_1.$ Similarly, for $j=2,$

$$\begin{array}{ccc}\hfill {\mathit{H}}_2& =& n\underset{3-2Sums}{\underbrace{\sum _{f_k=2}^{m_3}}}\mathit{z}{}_{f_3,1}\mathit{z}{}_{f_3,1}^{\prime }\sim {\mathcal{W}}_{m_1}\left(\left(m_3-1\right);{\mathsf{\Delta }}_{3,2}\right)\mathrm{and}\hfill \\ \hfill {\mathit{E}}_2& =& \left(n-1\right)\left(m_3-1\right){\widehat{\mathsf{\Delta }}}_{3,2}\sim {\mathcal{W}}_{m_1}\left(\left(n-1\right)\left(m_3-1\right);{\mathsf{\Delta }}_{3,2}\right),\mathrm{and}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill T_{02}^2& =& tr\left[n\sum _{f_3=2}^{m_3}\mathit{z}{}_{f_3,1}\mathit{z}{}_{f_3,1}^{\prime }{\widehat{\mathsf{\Delta }}}_{3,2}^{-1}\right]=\mathrm{tr}\left[{\mathit{H}}_2{\left(\frac{{\mathit{E}}_2}{\left(n-1\right)\left(m_3-1\right)}\right)}^{-1}\right]\hfill \\ & =& \left(n-1\right)\left(m_3-1\right)\mathrm{tr}\left({\mathit{H}}_2{E}_2^{-1}\right)=d_2\mathrm{tr}\left({\mathit{H}}_2{\mathit{E}}_2^{-1}\right),\hfill \end{array}$$

has a LH-trace distribution denoted by ${\mathcal{T}}_0^2\left(m_1;\left(m_3-1\right),d_2\right)$ if $d_2=\left(n-1\right)\left(m_3-1\right)>m_1.$ Note that the case $j=3$ reduces to ${\mathit{H}}_3=n\mathit{z}{}_{1,1}\mathit{z}{{}_{1,1}}^{\prime }\sim {\mathcal{W}}_{m_1}\left(1,{\mathsf{\Delta }}_{3,3}\right)$ and ${\mathit{E}}_3=\left(n-1\right){\widehat{\mathsf{\Delta }}}_{3,3}\sim \mathcal{W}\left(n-1;{\mathsf{\Delta }}_{3,3}\right),$ then $T_{03}^2$ has the LH-trace distribution ${\mathcal{T}}_0^2\left(m_1;1,d_3\right)$ if $d_3=\left(n-1\right)>m_1.$

Therefore, the distribution of $D^2={\displaystyle \sum _{j=1}^3}{T}_{0j}^2$ is the following convolution of three independent LH-trace distributions:

$$\begin{array}{ccc}& & \underset{j=1}{\overset{3}{⨁}}{\mathcal{T}}_0^2\left(m_1;p_{j+2,3}\left(m_{j+1}-1\right),d_j\right)=\hfill \\ & & {\mathcal{T}}_0^2\left(m_1;m_3\left(m_2-1\right),d_1\right)\oplus {\mathcal{T}}_0^2\left(m_1;\left(m_3-1\right),d_2\right)\oplus {\mathcal{T}}_0^2\left(m_1;1,d_3\right),\hfill \end{array}$$

(26)

if for $j=1,2,3$, $\left(n-1\right)p_{j+2,3}\left(m_{j+1}-1\right)\ge m_1.$ The critical values of this distribution can be obtained using simulations. LH-trace distribution is usually approximated by $\mathcal{F}$ distribution as mentioned before, however, we use here the second approximation suggested in McKeon [27].

For $j=1,$ denoting by $m_{\left(1\right)}=d_1-m_1-1=\left(n-1\right)m_3\left(m_2-1\right)-m_1-1,$ by $k_1=m_3\left(m_2-1\right)$ and by

$$\begin{array}{ccc}\hfill B_1& =& \frac{\left[m_{\left(1\right)}+k_1\right]}{\left[m_{\left(1\right)}-2\right]}·\frac{\left[m_{\left(1\right)}+m_1\right]}{\left[m_{\left(1\right)}+1\right]}\hfill \\ & =& \frac{\left[n{m}_3\left(m_2-1\right)-m_1-1\right]\left[\left(n-1\right)m_3\left(m_2-1\right)-1\right]}{\left[\left(n-1\right)m_3\left(m_2-1\right)-m_1-3\right]\left[\left(n-1\right)m_3\left(m_2-1\right)-m_1\right]},\hfill \end{array}$$

the distribution

$${\mathcal{T}}_0^2\left(m_1;m_3\left(m_2-1\right),\left(n-1\right)m_3\left(m_2-1\right)\right)$$

of $T_{01}^2=d_{1}tr\left({\mathit{H}}_1{\mathit{E}}_1^{-1}\right)$ can be approximated by $g_1\mathcal{F}\left(K_1,D_1\right),$ where $K_1=m_3\left(m_2-1\right)m_1,$ $D_1=4+\frac{m_3\left(m_2-1\right)m_1+2}{B_1-1}$ and $g_1=d_1\frac{K_1}{m_{\left(1\right)}}\frac{D_1-2}{D_1}=\frac{\left(n-1\right)m_3^2{\left(m_2-1\right)}^2{m}_1}{\left(n-1\right)m_3\left(m_2-1\right)-m_1-1}\frac{D_1-2}{D_1}.$

For $j=2,$ denoting by $m_{\left(2\right)}=d_2-m_1-1=\left(n-1\right)\left(m_3-1\right)-m_1-1,$ by $k_{\left(2\right)}=m_3-1$ and by

$$\begin{array}{ccc}\hfill B_2& =& \frac{\left[m_{\left(2\right)}+k_{\left(2\right)}\right]}{\left[m_{\left(2\right)}-2\right]}·\frac{\left[m_{\left(2\right)}+m_1\right]}{\left[m_{\left(2\right)}+1\right]}\hfill \\ & =& \frac{\left[n\left(m_3-1\right)-m_1-1\right]\left[\left(n-1\right)\left(m_3-1\right)-1\right]}{\left[\left(n-1\right)\left(m_3-1\right)-m_1-3\right]\left[\left(n-1\right)\left(m_3-1\right)-m_1\right]},\hfill \end{array}$$

the distribution

$${\mathcal{T}}_0^2\left(m_1;\left(m_3-1\right),d_2\right)$$

of $T_{02}^2=d_{2}tr\left({\mathit{H}}_2{\mathit{E}}_2^{-1}\right)$ can be approximated by $g_2\mathcal{F}\left(K_2,D_2\right),$ where $K_2=\left(m_3-1\right)m_1,$ $D_2=4+\frac{\left(m_3-1\right)m_1+2}{B_2-1}$ and $g_2=\left(n-1\right)\left(m_3-1\right)\frac{K_2}{m_{\left(2\right)}}\frac{D_2-2}{D_2}=\frac{\left(n-1\right){\left(m_3-1\right)}^2{m}_1}{\left(n-1\right)\left(m_3-1\right)-m_1-1}\frac{D_2-2}{D_2}.$

Finally, for our last case corresponding to $j=3,$ the distribution ${\mathcal{T}}_0^2\left(m_1;1,\left(n-1\right)\right)$ is the usual Hotelling ${\mathcal{T}}_{m_1,n-1}^2,$ that is, an exact distribution as follows:

$$\frac{\left(n-1\right)m_1}{n-m_1}\mathcal{F}\left(m_1,n-m_1\right).$$

This means that the distribution of $D^2$ can be approximated by the convolution of the above three $\mathcal{F}$ distributions (two approximated $\mathcal{F}$ distribution and one exact $\mathcal{F}$ distribution), where its critical values are obtained by the method suggested by Dyer [28].

Now, we need to perform the convolution of three distribution functions. Since convolution is associative, for three distribution functions ${\mathcal{F}}_1,{\mathcal{F}}_2,$ and ${\mathcal{F}}_3$, the associative law of random variables implies that $\left({\mathcal{F}}_1\otimes {\mathcal{F}}_2\right)\otimes {\mathcal{F}}_3={\mathcal{F}}_1\otimes \left({\mathcal{F}}_2\otimes {\mathcal{F}}_3\right)$, so we can dispense with the parentheses and can write ${\mathcal{F}}_1\otimes {\mathcal{F}}_2\otimes {\mathcal{F}}_3$. In the following section, we present the unbiased estimates of the eigenblocks for a 3-SSCS covariance matrix.

5.2. The Expressions of the $\mathsf{\Delta }$’s Estimators for the Case $k=3$

1.

From Lemma 5 in Leiva and Roy [10], the unbiased estimators ${\widehat{\mathit{U}}}_{3,j}$ of ${\mathit{U}}_{3,j}$ for each $j=1,2,3,$ are written as follows:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{U}}}_{3,1}& =& \frac{1}{p_{2,3}}{\mathrm{BT}}_{p_{1,1}}\left(\mathit{S}\right)=\frac{1}{p_{2,3}}\sum _{\mathit{f}\in \mathit{F}}{\mathit{S}}_{\mathit{f};\mathit{f}}\hfill \\ & =& \frac{1}{\left(n-1\right)q_{3,1}}\stackrel{2sums}{\overbrace{\sum _{f_3\in F_3}\sum _{f_2\in F_2}}}\sum _{r=1}^n\left(\underset{m_1\times 1}{{\mathit{x}}_{r;f_3,f_2}-{\overline{\mathit{x}}}_{f_3,f_2}}\right)\left(\underset{m_1\times 1}{{\mathit{x}}_{r;f_3,f_2}^{\prime }-{\overline{\mathit{x}}}_{f_3,f_2}^{\prime }}\right),\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{U}}}_{3,2}& =& \frac{{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,2}}\left(\mathit{S}\right)\right]-{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,1}}\left(\mathit{S}\right)\right]}{q_{3,2}}=\frac{{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,2}}\left(\mathit{S}\right)\right]-{\mathrm{BT}}_{p_{1,1}}\left(\mathit{S}\right)}{q_{3,2}}\hfill \\ \hfill & =& \frac{1}{(n-1)q_{3,2}}\stackrel{1sum}{\overbrace{\sum _{f_3\in F_3}}}\stackrel{1specialsumpair}{\overbrace{\left(\sum _{f_2\in F_2}\sum _{f_2\ne f_2^{*}\in F_2}\right)}}\sum _{r=1}^n\left(\underset{m_1\times 1}{{\mathit{x}}_{r;f_3,f_2}-{\overline{\mathit{x}}}_{f_3,f_2}}\right)\left(\underset{1\times m_1}{{\mathit{x}}_{r;f_3^{*},f_2^{*}}^{\prime }-{\overline{\mathit{x}}}_{f_3^{*},f_2^{*}}^{\prime }}\right),\hfill \end{array}$$

(28)

and

$$\begin{array}{ccc}\hfill {\widehat{\mathit{U}}}_{3,3}& =& \frac{{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,3}}\left(\mathit{S}\right)\right]-{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,2}}\left(\mathit{S}\right)\right]}{q_{3,3}}\hfill \\ \hfill & =& \frac{1}{(n-1)q_{3,3}}\stackrel{1specialsumpair}{\overbrace{\left(\sum _{f_3\in F_3}\sum _{f_3\ne f_3^{*}\in F_3}\right)}}\stackrel{1sumpair}{\overbrace{\left(\sum _{f_2\in F_2}\sum _{f_2^{*}\in F_2}\right)}}\sum _{r=1}^n\left(\underset{m_1\times 1}{{\mathit{x}}_{r;f_3,f_2}-{\overline{\mathit{x}}}_{f_3,f_2}}\right)\left(\underset{1\times m_1}{{\mathit{x}}_{r;f_3^{*},f_2^{*}}^{\prime }-{\overline{\mathit{x}}}_{f_3^{*},f_2^{*}}^{\prime }}\right),\hfill \\ \hfill \end{array}$$

(29)

where $q_{3,j},j=1,2,3$ are given in (19). Therefore, an unbiased estimator of ${\mathsf{\Gamma }}_3$ is given by:

$$\begin{array}{c}\hfill {\widehat{\mathsf{\Gamma }}}_3=\left[\sum _{j=1}^2{\mathit{I}}_{p_{j+1,3}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\widehat{\mathit{U}}}_{3,j}-{\widehat{\mathit{U}}}_{3,j+1}\right)\right]+{\mathit{J}}_{p_{2,3}}\otimes {\widehat{\mathit{U}}}_{3,3}.\end{array}$$

Since k-SSCS matrix ${\mathsf{\Gamma }}_k$ in (9) is of Jordan algebra type, following Kozioł et al. [26] one can show that the above estimate ${\widehat{\mathsf{\Gamma }}}_3$ is the best unbiased, consistent and complete estimator for ${\mathsf{\Gamma }}_3$.

2.

For each $j=1,2,3$ an unbiased estimator of ${\mathsf{\Delta }}_{3,j},$${\widehat{\mathsf{\Delta }}}_{3,j}$ is given by:

$${\widehat{\mathsf{\Delta }}}_{3,j}=\sum _{i=1}^j{p}_{2,i}\left({\widehat{\mathit{U}}}_{3,i}-{\widehat{\mathit{U}}}_{3,i+1}\right),$$

where ${\widehat{\mathit{U}}}_{3,4}=\mathbf{0}$ and $p_{2,1}=1,$ or equivalently:

$${\widehat{\mathsf{\Delta }}}_{3,j}=\left\{\begin{array}{ccc}{\widehat{\mathit{U}}}_{3,1}-{\widehat{\mathit{U}}}_{3,2}& if& j=1\\ {\widehat{\mathsf{\Delta }}}_{3,j-1}+p_{2,j}\left({\widehat{\mathit{U}}}_{3,j}-{\widehat{\mathit{U}}}_{3,j+1}\right)& if& j=2,3\end{array}\right.,$$

The above unbiased estimators admit the following expressions as functions of $\mathit{S}$:

$${\widehat{\mathsf{\Delta }}}_{3,j}=\frac{m_{j+1}·{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,j}}\left(\mathit{S}\right)\right]-{\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,j+1}}\left(\mathit{S}\right)\right]}{p_{2,3}\left(m_{j+1}-1\right)},$$

for $j=1,2$ and where ${\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,1}}\left(\mathit{S}\right)\right]={\mathrm{BT}}_{p_{1,1}}\left(\mathit{S}\right)$ and ${\mathrm{BS}}_{p_{1,1}}\left[{\mathrm{BT}}_{p_{1,j+1}}\left(\mathit{S}\right)\right]={\mathrm{BS}}_{p_{1,1}}\left(\mathit{S}\right),$ and an unbiased estimator of ${\mathsf{\Delta }}_{3,3}$, ${\widehat{\mathsf{\Delta }}}_{3,3}$ is given by:

$${\widehat{\mathsf{\Delta }}}_{3,3}=\frac{1}{p_{2,3}}{\mathrm{BS}}_{p_{1,1}}\left(\mathit{S}\right).$$

6. Test for the Equality of Two Means

6.1. Paired Observation Model

In this section, we consider that in each one of the n individuals, a $p_{1,k}$-variate vector is measured at two different times (e.g., before and after a treatment). These measurements are k-th order (array-variate) measurements from each individual. To be more precise, for each $\mathit{f}={\left(f_k,\dots ,f_2\right)}^{\prime }\in F,$ let ${\mathit{v}}_{r;f_k,\dots ,f_2}$ and ${\mathit{w}}_{r;f_k,\dots ,f_2}$ be the paired $m_1$-dimensional vectors measured at the $\left(f_k,\dots ,f_2\right)$ site of the r individuals, for $r=1,\dots ,n.$ Let ${\mathit{u}}_r$ be the partitioned $2m$-variate vectors ${\mathit{u}}_r^{\prime }=\left({\mathit{v}}_r^{\prime },{\mathit{w}}_r^{\prime }\right)=\left(\left({\mathit{v}}_{r,1}^{\prime },\dots ,{\mathit{v}}_{r,m_k}^{\prime }\right),\left({\mathit{w}}_{r,1}^{\prime },\dots ,{\mathit{w}}_{r,m_k}^{\prime }\right)\right),$ where ${\mathit{v}}_{r;f_k}={\left({\mathit{v}}_{r;f_k,1}^{\prime },\dots ,{\mathit{v}}_{r;f_{k-1},m_{k-1}}^{\prime }\right)}^{\prime }$ and ${\mathit{w}}_{r;f_k}={\left({\mathit{w}}_{r;f_k,1}^{\prime },\dots ,{\mathit{w}}_{r;f_k,m_{k-1}}^{\prime }\right)}^{\prime }$ where ${\mathit{v}}_{r;f_k,f_{k-1}}={\left({\mathit{v}}_{r;f_k,f_{k-1},1}^{\prime },\dots ,{\mathit{v}}_{r;f_k,f_{k-1},m_{k-2}}^{\prime }\right)}^{\prime }$ and ${\mathit{w}}_{r;f_k,f_{k-1}}={\left({\mathit{w}}_{r;f_k,f_{k-1},1}^{\prime },\dots ,{\mathit{w}}_{r;f_k,f_{k-1},m_{k-2}}^{\prime }\right)}^{\prime },$ with $\mathit{v}{}_{r;f_k,\dots ,f_3}$$={\left({\mathit{v}}_{r;f_k,\dots ,f_3,1}^{\prime },\dots ,{\mathit{v}}_{r;f_k,\dots ,f_3,m_2}^{\prime }\right)}^{\prime }$ and ${\mathit{w}}_{r;f_k,\dots ,f_3}={\left({\mathit{w}}_{r;f_k,\dots ,f_3,1}^{\prime },\dots ,{\mathit{w}}_{r;f_k,\dots ,f_3,m_2}^{\prime }\right)}^{\prime }$, respectively, where ${\mathit{v}}_{r;f_k,\dots ,f_3,f_2}^{\prime }=\left(v_{r;f_k,\dots ,f_3,f_2,1},\dots ,v_{r;f_k,\dots ,f_3,f_2,m_1}\right)$ and ${\mathit{w}}_{r;f_k,\dots ,f_3,f_2}^{\prime }$$=\left(w_{r;f_k,\dots ,f_3,f_2,1},\dots ,w_{r;f_k,\dots ,f_3,f_2,m_1}\right)$ are the $m_1$ paired measurements taken from the rth individual, for $r=1,\dots ,n.$ We assume that ${\mathit{u}}_1,\dots ,{\mathit{u}}_n\stackrel{i.i.d.}{\sim }{\mathcal{N}}_{2p_{1,k}}\left({\mathbf{\mu }}_{\mathit{u}},{\mathsf{\Gamma }}_{\mathit{u}}\right),$ where i.i.d. stands for independent and identically distributed, and ${\mathbf{\mu }}_{\mathit{u}}={\left({\mathbf{\mu }}_{\mathit{v}}^{\prime },{\mathbf{\mu }}_{\mathit{w}}^{\prime }\right)}^{\prime }$ and ${\mathsf{\Gamma }}_{\mathit{u}}$ is the partitioned $2p_{1,k}\times 2p_{1,k}$-matrix

$${\mathsf{\Gamma }}_{\mathit{u}}=\left(\begin{array}{cc}\underset{p_{1,k}\times p_{1,k}}{{\mathsf{\Gamma }}_{\mathit{v}}}& \underset{p_{1,k}\times p_{1,k}}{{\mathsf{\Gamma }}_{\mathit{vw}}}\\ \underset{p_{1,k}\times p_{1,k}}{{\mathsf{\Gamma }}_{\mathit{wv}}}& \underset{p_{1,k}\times p_{1,k}}{{\mathsf{\Gamma }}_{\mathit{w}}}\end{array}\right),$$

where

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_{\mathit{v}}& =& \left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{\mathit{v},j}-{\mathit{U}}_{\mathit{v},j+1}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{\mathit{v},k}\mathrm{and}\hfill \\ \hfill {\mathsf{\Gamma }}_{\mathit{w}}& =& \left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{\mathit{w},j}-{\mathit{U}}_{\mathit{w},j+1}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{\mathit{w},k},\mathrm{and}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_{\mathit{vw}}& =& {\mathsf{\Gamma }}_{\mathit{wv}}^{\prime }=\left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{\mathit{vw},j}-{\mathit{U}}_{\mathit{vw},j+1}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{\mathit{vw},k}.\hfill \end{array}$$

The matrices ${\mathsf{\Gamma }}_{\mathit{vw}}$ and ${\mathsf{\Gamma }}_{\mathit{wv}}$ are accountable for the linear dependence among the considered $m_1$ paired measurements. Particular cases of ${\mathsf{\Gamma }}_{\mathit{vw}}$ could be of interest, e.g., ${\mathit{U}}_{\mathit{vw},1}=\cdots ={\mathit{U}}_{\mathit{vw},k},$ that is, ${\mathsf{\Gamma }}_{\mathit{vw}}={\mathsf{\Gamma }}_{\mathit{wv}}^{\prime }={\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{\mathit{vw},k}$(see, for example, when $k=2,$ Definition 2 on Page 388 in Leiva [6]). Under this set up, we are interested in testing the following hypothesis:

$${\mathcal{H}}_0:{\mathbf{\mu }}_{\mathit{v}}={\mathbf{\mu }}_{\mathit{w}}vs.{\mathcal{H}}_1:{\mathbf{\mu }}_{\mathit{v}}\ne {\mathbf{\mu }}_{\mathit{w}}.$$

If we define $\mathit{x}=\mathit{v}-\mathit{w}$ the above hypothesis is equivalent to

$${\mathcal{H}}_0:{\mathbf{\mu }}_{\mathit{x}}=\mathbf{0}vs{\mathcal{H}}_1:{\mathbf{\mu }}_{\mathit{x}}\ne \mathbf{0},$$

as ${\mathbf{\mu }}_{\mathit{x}}=\mathrm{E}(\mathit{v}-\mathit{w})={\mathbf{\mu }}_{\mathit{v}}-{\mathbf{\mu }}_{\mathit{w}}.$ Moreover, ${\mathit{x}}_r=$${\mathit{v}}_r-{\mathit{w}}_r,$$r=1$$,\dots ,n,$ are i.i.d. $\mathcal{N}\left({\mathbf{\mu }}_{\mathit{x}};\mathsf{\Gamma }\right)$ where ${\mathbf{\mu }}_{\mathit{x}}={\mathbf{\mu }}_{\mathit{v}}-{\mathbf{\mu }}_{\mathit{w}}$ and

$$\begin{array}{ccc}\hfill \mathsf{\Gamma }& =& \mathrm{cov}\left(\mathit{x}\right)=\mathrm{cov}(\mathit{v}-\mathit{w})={\mathsf{\Gamma }}_{\mathit{v}}+{\mathsf{\Gamma }}_{\mathit{w}}-{\mathsf{\Gamma }}_{\mathit{vw}}-{\mathsf{\Gamma }}_{\mathit{wv}}\hfill \\ & =& \left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{\mathit{x},j}-{\mathit{U}}_{\mathit{x},j+1}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{\mathit{x},k},\hfill \end{array}$$

$$\mathrm{where}{\mathit{U}}_{\mathit{x},j}={\mathit{U}}_{\mathit{v},j}+{\mathit{U}}_{\mathit{w},j}-{\mathit{U}}_{\mathit{vw},j}-{\mathit{U}}_{\mathit{vw},j}^{\prime }\mathrm{for}j=1,\dots ,k.$$

Assuming $\mathsf{\Gamma }$ is a positive definite matrix and that $n>p_{1,k},$ one may consider the likelihood ratio test for the above hypothesis testing problem for k- level multivariate data assuming the mean vectors ${\mathbf{\mu }}_{\mathit{v}}$ and ${\mathbf{\mu }}_{\mathit{w}}$ are unstructured. Note that this test problem reduces to the one sample mean case of the previous section where ${\mathbf{\mu }}_0=\mathbf{0}.$ Therefore, all the results obtained in the previous section are valid for this case. Following the same logic as in Remark 1, the needed sample size for the test is $n\ge m_1+1$, regardless of any $m_j,j=2,\dots ,k$.

6.2. Independent Observation Model

In this section, we consider the case where we have two independent samples: one random sample of size $n_1$ of $p_{1,k}$ vectors ${\mathit{v}}_{r_1}:r_1=1,\dots ,n_1\stackrel{i.i.d.}{\sim }{\mathcal{N}}_{p_{1,k}}\left({\mathbf{\mu }}_{\mathit{v}},{\mathsf{\Gamma }}_{\mathit{v}}\right),$ with ${\mathit{v}}_{r_1}^{\prime }=\left({\mathit{v}}_{r_1,1}^{\prime },\dots ,{\mathit{v}}_{r_1,m_k}^{\prime }\right),$ where ${\mathit{v}}_{r_1;f_k}^{\prime }=\left({\mathit{v}}_{r_1;f_k,1}^{\prime },\dots ,{\mathit{v}}_{r_1;f_k,m_{k-1}}^{\prime }\right)$, ${\mathit{v}}_{r_1;f_k,f_{k-1}}^{\prime }=\left({\mathit{v}}_{r_1;f_k,f_{k-1},1}^{\prime },\dots ,{\mathit{v}}_{r_1;f_k,f_{k-1},m_{k-2}}^{\prime }\right),$ …, and where ${\mathit{v}}_{r_1;f_k,\dots ,f_3}^{\prime }=\left({\mathit{v}}_{r_1;f_k,\dots ,f_3,1}^{\prime },\dots ,{\mathit{v}}_{r_1;f_k,\dots ,f_3,m_2}^{\prime }\right)$ with ${\mathit{v}}_{r_1;f_k,\dots ,f_3,f_2}^{\prime }=\left(v_{r_1;f_k,\dots ,f_3,f_2,1},\dots ,v_{r_1;f_k,\dots ,f_3,f_2,m_1}\right)$ being the $m_1$ measurements taken from the $r_1^{th}$ individual, for $r_1=1,\dots ,n_1,$ and another random sample of size $n_2$ of $p_{1,k}$ vectors ${\mathit{w}}_{r_2}:r_2=1,\dots ,n_2\stackrel{i.i.d.}{\sim }{\mathcal{N}}_{p_{1,k}}\left({\mathbf{\mu }}_{\mathit{w}},{\mathsf{\Gamma }}_{\mathit{w}}\right)$, with ${\mathit{w}}_{r_2}^{\prime }=\left({\mathit{w}}_{r_2,1}^{\prime },\dots ,{\mathit{w}}_{r_2,m_k}^{\prime }\right),$ where ${\mathit{w}}_{r_2;f_k}^{\prime }=\left({\mathit{w}}_{r_2;f_k,1}^{\prime },\dots ,{\mathit{w}}_{r_2;f_k,m_{k-1}}^{\prime }\right)$, where ${\mathit{w}}_{r_2;f_k,f_{k-1}}^{\prime }=\left({\mathit{w}}_{r_2;f_k,f_{k-1},1}^{\prime },\dots ,{\mathit{w}}_{r_2;f_k,f_{k-1},m_{k-2}}^{\prime }\right)$, …, and where ${\mathit{w}}_{r_2;f_k,\dots ,f_3}^{\prime }=\left({\mathit{w}}_{r_2;f_k,\dots ,f_3,1}^{\prime },\dots ,{\mathit{w}}_{r_2;f_k,\dots ,f_3,m_1}^{\prime }\right)$ with ${\mathit{w}}_{r_2;f_k,\dots ,f_3,f_2}^{\prime }=\left(w_{r_2;f_k,\dots ,f_3,f_2,1},\dots ,w_{r_2;f_k,\dots ,f_3,f_2,m_1}\right)$ being the $m_1$ measurements taken from the $r_2^{th}$ individual, for $r_2=1,\dots ,n_2,$ for $\mathit{f}={\left(f_k,\dots ,f_2\right)}^{\prime }\in F$.

Our interest is in testing the following hypothesis:

$${\mathcal{H}}_0:{\mathbf{\mu }}_{\mathit{v}}={\mathbf{\mu }}_{\mathit{w}}vs.{\mathcal{H}}_1:{\mathbf{\mu }}_{\mathit{v}}\ne {\mathbf{\mu }}_{\mathit{w}},$$

(30)

under the assumption that ${\mathsf{\Gamma }}_{\mathit{v}}={\mathsf{\Gamma }}_{\mathit{w}}=\mathsf{\Gamma }={\mathsf{\Gamma }}_k$ is an unknown k-SSCS covariance matrix of the form (8). Let $\mathit{V}=\left({\mathit{v}}_1,\dots ,{\mathit{v}}_{n_1}\right)$ and $\mathit{W}=\left({\mathit{w}}_1,\dots ,{\mathit{w}}_{n_2}\right)$ denote the corresponding two sample matrix data. We know that the sample means $\overline{\mathit{v}}=\frac{1}{n_1}{\displaystyle \sum _{r_1=1}^{n_1}}{\mathit{v}}_{r_1}$ and $\overline{\mathit{w}}=\frac{1}{n_2}{\displaystyle \sum _{r_2=1}^{n_2}}{\mathit{w}}_{r_2}$ are independent of the covariance matrix estimators ${\mathit{S}}_1=\frac{1}{n_1-1}\mathit{V}{\mathit{Q}}_{n_1}{\mathit{V}}^{\prime }$ and ${\mathit{S}}_2=\frac{1}{n_2-1}\mathit{W}{\mathit{Q}}_{n_2}{\mathit{W}}^{\prime }$ respectively. Therefore, they are also independent of $S^{\mathrm{pl}},$ the pooled unbiased estimator $\mathsf{\Gamma }$ (convex linear combination of unbiased estimators of $\mathsf{\Gamma }),$ which is given by:

$$\begin{array}{ccc}\hfill {\mathit{S}}^{\mathrm{pl}}& =& \frac{n_1-1}{n_1+n_2-2}{\mathit{S}}_1+\frac{n_2-1}{n_1+n_2-2}{\mathit{S}}_2\hfill \\ & =& {\left({\mathit{S}}_{f_k,\dots ,f_2;f_k^{*},\dots ,f_2^{*}}^{\mathrm{pl}}\right)}_{f_2,f_2^{*}=1;\dots ;f_k,f_k^{*}=1}^{m_2;\dots ;m_k}={\left({\mathit{S}}_{\mathit{f};{\mathit{f}}^{*}}^{\mathrm{pl}}\right)}_{\mathit{f};{\mathit{f}}^{*}\in F},\hfill \end{array}$$

where

$${\mathit{S}}_{\mathit{f};{\mathit{f}}^{*}}^{\mathrm{pl}}=\frac{{\displaystyle \sum _{r_1=1}^{n_1}}\left({\mathit{v}}_{r_1;\mathit{f}}-{\overline{\mathit{v}}}_{\mathit{f}}\right)\left({\mathit{v}}_{r_1;{\mathit{f}}^{*}}^{\prime }-{\overline{\mathit{v}}}_{{\mathit{f}}^{*}}^{\prime }\right)+{\displaystyle \sum _{r_2=1}^{n_2}}\left({\mathit{w}}_{r_2;\mathit{f}}-{\overline{\mathit{w}}}_{\mathit{f}}\right)\left({\mathit{w}}_{r_2;{\mathit{f}}^{*}}^{\prime }-{\overline{\mathit{w}}}_{{\mathit{f}}^{*}}^{\prime }\right)}{\left(n_1+n_2-2\right)}.$$

Now

$$\overline{\mathit{v}}-\overline{\mathit{w}}\sim {\mathcal{N}}_{p_{1,k}}\left({\mathbf{\mu }}_{\mathit{v}}-{\mathbf{\mu }}_{\mathit{w}},\frac{1}{n^{\mathrm{pl}}}{\mathsf{\Gamma }}_k\right),\mathrm{where}{n}^{\mathrm{pl}}=\frac{n_1{n}_2}{n_1+n_2}.$$

We know that under $H_0$:

$$\begin{array}{c}\hfill \overline{\mathit{v}}-\overline{\mathit{w}}\sim {\mathcal{N}}_{p_{1,k}}\left(\mathbf{0},\frac{1}{n^{\mathrm{pl}}}{\mathsf{\Gamma }}_k\right),\mathrm{and}{\mathit{S}}^{\mathrm{pl}}\sim {\mathcal{W}}_{p_{1,k}}\left(n_1+n_2-2,\frac{1}{n_1+n_2-2}{\mathsf{\Gamma }}_k\right).\end{array}$$

Due to the exchangeable form of ${\mathsf{\Gamma }}_k,$ it is clear that we again have:

$$\begin{array}{ccc}\hfill \mathrm{E}\left[{\mathit{S}}_{\mathit{f};{\mathit{f}}^{*}}^{\mathrm{pl}}\right]& =& \mathrm{E}\left[{\mathit{S}}_{f_k,\dots ,f_2;f_k^{*},\dots ,f_2^{*}}^{\mathrm{pl}}\right]\hfill \\ & =& \left\{\begin{array}{ccccc}{\mathit{U}}_{k,1}& if& f_2=f_2^{*},\dots ,f_k=f_k^{*}& & \\ {\mathit{U}}_{k,j}& if& f_j\ne f_j^{*},f_{j+1}=f_{j+1}^{*},\dots ,f_k=f_k^{*}& for& j=2,\dots ,k\end{array}\right..\hfill \end{array}$$

(31)

Note that each of the following expressions is the arithmetic mean of all submatrices of ${\mathit{S}}^{\mathrm{pl}}$, which, according to (31), have the same expectation. It is easy to prove that for each $j=1,\dots ,k,$ an unbiased estimator of ${\mathit{U}}_{k,j},$${\widehat{\mathit{U}}}_{k,j}^{\mathrm{pl}},$ for $j=1$ is given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{U}}}_{k,1}^{\mathrm{pl}}& =& \frac{1}{p_{2,k}}{\mathrm{BT}}_{p_{1,1}}\left({\mathit{S}}^{\mathrm{pl}}\right)=\frac{1}{p_{2,k}}\sum _{\mathit{f}\in \mathit{F}}{\mathit{S}}_{\mathit{f};\mathit{f}}^{\mathrm{pl}}\hfill \\ & =& \frac{n_1-1}{n_1+n_2-2}\frac{1}{p_{2,k}}\sum _{\mathit{f}\in \mathit{F}}{\mathit{S}}_{1:\mathit{f};\mathit{f}}+\frac{n_2-1}{n_1+n_2-2}\frac{1}{p_{2,k}}\sum _{\mathit{f}\in \mathit{F}}{\mathit{S}}_{2:\mathit{f};\mathit{f}},\hfill \end{array}$$

and for $j=2,\dots ,k,$${\widehat{\mathit{U}}}_{k,j}^{\mathrm{pl}}$ after some algebraic simplification is given by:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{U}}}_{k,j}^{\mathrm{pl}}& =& \frac{{\mathrm{BS}}_{p_{1,1}}\left[\left({\mathrm{BT}}_{p_{1,j}}\left({\mathit{S}}^{\mathrm{pl}}\right)\right)-{\mathrm{BT}}_{p_{1,j-1}}\left({\mathit{S}}^{\mathrm{pl}}\right)\right]}{q_{k,j}}\hfill \\ & =& \frac{n_1-1}{n_1+n_2-2}\frac{1}{q_{k,j}}\sum _{f_k\in F_k}\cdots \sum _{f_{j+1}\in F_{j+1}}\left(\sum _{f_j\in F_j}\sum _{f_j\ne f_j^{*}\in F_j}\right)\left(\sum _{f_{j-1}\in F_{j-1}}\sum _{f_{j-1}^{*}\in F_{j-1}}\right)\cdots \left(\sum _{f_2\in F_k}\sum _{f_2^{*}\in F_k}\right){\mathit{S}}_{1:\mathit{f};{\mathit{f}}^{*}}\hfill \\ & +& \frac{n_2-1}{n_1+n_2-2}\frac{1}{q_{k,j}}\sum _{f_k\in F_k}\cdots \sum _{f_{j+1}\in F_{j+1}}\left(\sum _{f_j\in F_j}\sum _{f_j\ne f_j^{*}\in F_j}\right)\left(\sum _{f_{j-1}\in F_{j-1}}\sum _{f_{j-1}^{*}\in F_{j-1}}\right)\cdots \left(\sum _{f_2\in F_k}\sum _{f_2^{*}\in F_k}\right){\mathit{S}}_{2:\mathit{f};{\mathit{f}}^{*}},\hfill \end{array}$$

where $q_{k,j}$ is given in (19). Therefore, we can use the following unbiased estimator of variance and covariance matrices

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\Gamma }}}^{\mathrm{pl}}& =& {\widehat{\mathsf{\Gamma }}}_k^{\mathrm{pl}}=\left[\sum _{j=1}^{k-1}{\mathit{I}}_{p_{j+1,k}}\otimes {\mathit{J}}_{p_{2,j}}\otimes \left({\mathit{U}}_{k,j}^{\mathrm{pl}}-{\mathit{U}}_{k,j+1}^{\mathrm{pl}}\right)\right]+{\mathit{J}}_{p_{2,k}}\otimes {\mathit{U}}_{k,k}^{\mathrm{pl}},\hfill \\ \hfill \mathrm{and}{\widehat{\mathsf{\Omega }}}_k^{\mathrm{pl}}& =& \frac{1}{n^{\mathrm{pl}}}{\widehat{\mathit{D}}}_k^{\mathrm{pl}}={\mathit{L}}_k^{\prime }\left(\frac{1}{n^{\mathrm{pl}}}{\widehat{\mathsf{\Gamma }}}_k^{\mathrm{pl}}\right){\mathit{L}}_k=\frac{1}{n^{\mathrm{pl}}}\left({\mathit{H}}_{m_k,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right){\widehat{\mathsf{\Gamma }}}_k^{\mathrm{pl}}\left({\mathit{H}}_{m_k,m_2}\otimes {\mathit{I}}_{m_1}\right)\hfill \\ & =& \frac{1}{n^{\mathrm{pl}}}\mathrm{diag}\left\{{\widehat{\mathit{D}}}_{f_k,f_{k-1},\dots ,f_2}^{\mathrm{pl}};{\left(f_k,f_{k-1},\dots ,f_2\right)}^{\prime }\in F_k\times F_{k-1}\times \dots \times F_2\right\},\hfill \end{array}$$

where, for each $j=1,\dots ,k,$ the diagonal $m_1\times m_1$-matrices ${\mathit{D}}_{\mathit{f}}^{\mathrm{pl}}={\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}^{\mathrm{pl}}$ are given by:

$$\begin{array}{ccc}{\widehat{\mathit{D}}}_{f_k,f_{k-1},\dots ,f_2}^{\mathrm{pl}}={\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}& if& f_2=1,\dots ,f_j=1,f_{j+1}\ne 1,\end{array}$$

where $f_{k+1}\ne 1$ is not taken into consideration and where

$${\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}=\sum _{i=1}^j{p}_{2,i}\left({\widehat{\mathit{U}}}_{k,i}^{\mathrm{pl}}-{\widehat{\mathit{U}}}_{k,i+1}^{\mathrm{pl}}\right),$$

with ${\widehat{\mathit{U}}}_{k,k+1}^{\mathrm{pl}}=\mathbf{0}$, or equvalently:

$${\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}=\left\{\begin{array}{ccc}{\widehat{\mathit{U}}}_{k,1}^{\mathrm{pl}}-{\widehat{\mathit{U}}}_{k,2}^{\mathrm{pl}}& if& j=1\\ {\widehat{\mathsf{\Delta }}}_{k,j-1}^{\mathrm{pl}}+p_{2,j}\left({\widehat{\mathit{U}}}_{k,j}^{\mathrm{pl}}-{\widehat{\mathit{U}}}_{k,j+1}^{\mathrm{pl}}\right)& if& j=2,\dots ,k\end{array}\right..$$

The usual likelihood ratio test of (30) is to reject $H_0$ if:

$$T^2=\frac{n_1{n}_2}{n_1+n_2}{\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)}^{\prime }{\left({\mathit{S}}^{\mathrm{pl}}\right)}^{-1}\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)>c^2.$$

Since $\left(n_1+n_2-2\right){\mathit{S}}^{\mathrm{pl}}=\left(n_1-1\right){\mathit{S}}_1+\left(n_2-1\right){\mathit{S}}_2\sim {\mathcal{W}}_{m_1}\left(\left(n_1+n_2-2\right);{\mathsf{\Gamma }}_k\right)$,

$$\begin{array}{ccc}\hfill T^2& =& \underset{\sim Norma{l}_{p_{1,k}}}{\underbrace{{\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}^{-\frac{1}{2}}{\left[\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)-\left({\mathbf{\mu }}_{\mathit{v}}-{\mathbf{\mu }}_{\mathit{w}}\right)\right]}^{\prime }}}{\left(\underset{\sim Wishart/d.f.}{\underbrace{\frac{\left(n_1+n_2-2\right){\mathit{S}}^{\mathrm{pl}}}{\left(n_1+n_2-2\right)}}}\right)}^{-1}\underset{\sim Norma{l}_{p_{1,k}}}{\underbrace{{\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}^{-\frac{1}{2}}\left[\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)-\left({\mathbf{\mu }}_{\mathit{v}}-{\mathbf{\mu }}_{\mathit{w}}\right)\right]}}.\hfill \end{array}$$

Nevertheless, we cannot use the above result, as in our case, ${\widehat{\mathsf{\Gamma }}}_k^{\mathrm{pl}}$ is an estimator of ${\mathsf{\Gamma }}_k$. However, by Theorem 1 of Leiva and Roy [10], we know that the random vectors $\left(n-1\right)p_{j+2,k}\left(m_{j+1}-1\right){\widehat{\mathsf{\Delta }}}_{k,j}=\left(n-1\right){\mathrm{BT}}_{p_{1,1}}\left({\mathit{R}}_{k,j+1}\mathit{S}{\mathit{R}}_{k,j+1}\right):j=1,\dots ,k-1,$ and $\left(n-1\right){\widehat{\mathsf{\Delta }}}_{k,k}=\left(n-1\right){\mathrm{BT}}_{p_{1,1}}\left({\mathit{R}}_{k,j+1}\mathit{S}{\mathit{R}}_{k,j+1}\right),$ where ${\mathit{R}}_{k,j+1}={\mathit{R}}_{k,j+1}^{*}\otimes {\mathit{I}}_{m_1}$ are given by (4) with (5) and (6) with (7) are independent and

$$\begin{array}{ccc}\hfill \left(n_1+n_2-2\right)p_{j+2,k}\left(m_{j+1}-1\right){\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}& \sim & {\mathcal{W}}_{m_1}\left({\mathsf{\Delta }}_{k,j};\left(n_1+n_2-2\right)p_{j+2,k}\left(m_{j+1}-1\right)\right)\hfill \\ & & \mathrm{for}j=1,\dots ,k-1,\hfill \\ \hfill \mathrm{and}\left(n_1+n_2-2\right){\widehat{\mathsf{\Delta }}}_{k,k}^{\mathrm{pl}}& \sim & {\mathcal{W}}_{m_1}\left({\mathsf{\Delta }}_{k,k};\left(n_1+n_2-2\right)\right),\mathrm{where}{p}_{k+1,k}=1.\hfill \end{array}$$

Since the estimators ${\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}:j=1,\dots ,k$ are functions of ${\mathit{S}}^{\mathrm{pl}},$ they are independent of $\overline{\mathit{v}}-\overline{\mathit{w}}.$ Therefore, using a similar procedure as in the one sample case where we used the transformation $\mathit{z}={\mathit{L}}_k\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)=\left({\mathit{H}}_{m_k}\otimes \cdots \otimes {\mathit{H}}_{m_2}\otimes {\mathit{I}}_{m_1}\right)\left(\overline{\mathit{x}}-{\mathbf{\mu }}_0\right)$, we now use the following transformation:

$$\begin{array}{ccc}\hfill \mathit{d}& =& {\left(\mathit{d}{{}_{1,\dots ,1}}^{\prime },\cdots ,\mathit{d}{{}_{m_k,\dots ,m_2}}^{\prime }\right)}^{\prime }={\mathit{L}}_k^{\prime }\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)\left({\mathit{H}}_{m_k}\otimes \cdots \otimes {\mathit{H}}_{m_2}\otimes {\mathit{I}}_{m_1}\right)\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right).\hfill \end{array}$$

According to the previous result, $\mathit{d}={\mathit{L}}_k^{\prime }\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)\sim {\mathcal{N}}_{p_{1,k}}\left(\mathbf{0},{\mathsf{\Omega }}_k\right),$ where

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_k& =& \frac{1}{n}{\mathit{D}}_k={\mathit{L}}_k^{\prime }\left(\frac{1}{n}{\mathsf{\Gamma }}_k\right){\mathit{L}}_k=\frac{1}{n}\left({\mathit{H}}_{m_k,m_2}^{\prime }\otimes {\mathit{I}}_{m_1}\right){\mathsf{\Gamma }}_k\left({\mathit{H}}_{m_k,m_2}\otimes {\mathit{I}}_{m_1}\right)\hfill \\ & =& \frac{1}{n}\mathrm{diag}\left\{{\mathit{D}}_{f_k,f_{k-1},\dots ,f_2};{\left(f_k,f_{k-1},\dots ,f_2\right)}^{\prime }\in F_k\times F_{k-1}\times \dots \times F_2\right\},\hfill \end{array}$$

where, for each $j=1,\dots ,k,$ the diagonal $m_1\times m_1$- matrices ${\mathit{D}}_{\mathit{f}}={\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}$ are given by:

$$\begin{array}{ccc}{\mathit{D}}_{f_k,f_{k-1},\dots ,f_2}={\mathsf{\Delta }}_{k,j}& \mathrm{if}& f_2=1,\dots ,f_j=1,f_{j+1}\ne 1.\end{array}$$

Using a similar result as the one used in the one sample case, we obtain the statistic ${\left({\mathit{D}}^{\mathrm{pl}}\right)}^2$ as follows:

$$\begin{array}{ccc}\hfill {\left({\mathit{D}}^{\mathrm{pl}}\right)}^2& =& \frac{n_1{n}_2}{n_1+n_2}{\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)}^{\prime }{\left({\widehat{\mathsf{\Gamma }}}^{\mathrm{pl}}\right)}^{-1}\left(\overline{\mathit{v}}-\overline{\mathit{w}}\right)=n^{\mathrm{pl}}{\overline{\mathit{d}}}^{\prime }{\widehat{\mathsf{\Gamma }}}_k^{-1}\overline{\mathit{d}}\hfill \\ & =& n^{\mathrm{pl}}\sum _{j=1}^{k-1}\left\{\underset{k-jSums}{\underbrace{\sum _{f_k=1}^{m_k}\dots \sum _{f_{j+2}=1}^{m_{j+2}}\sum _{f_{j+1}=2}^{m_{j+1}}}}\mathit{d}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\left({\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}\right)}^{-1}\mathit{d}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}\right\}\hfill \\ & & +n^{\mathrm{pl}}\mathit{d}{}_{f_k,\dots ,f_{j+1},\underset{k-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\left({\widehat{\mathsf{\Delta }}}_{k,k}^{\mathrm{pl}}\right)}^{-1}\mathit{d}{}_{\underset{k-1ones}{\underbrace{1,\dots ,1}}}\hfill \\ & =& \sum _{j=1}^k\left[tr\left\{\frac{n_1{n}_2}{n_1+n_2}\underset{k-jSums}{\underbrace{\sum _{f_k=1}^{m_k}\dots \sum _{f_{j+2}=1}^{m_{j+2}}\sum _{f_{j+1}=2}^{m_{j+1}}}}\mathit{d}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}\mathit{d}{}_{f_k,\dots ,f_{j+1},\underset{j-1ones}{\underbrace{1,\dots ,1}}}^{\prime }{\left({\widehat{\mathsf{\Delta }}}_{k,j}^{\mathrm{pl}}\right)}^{-1}\right\}\right]\hfill \\ & =:& \sum _{j=1}^k{T}_{0j}^2.\hfill \end{array}$$

Then, the distribution of ${\left({\mathit{D}}^{\mathrm{pl}}\right)}^2={\displaystyle \sum _{j=1}^k}{T}_{0j}^2$ is the convolution of k independent LH-trace distributions as follows:

$$\underset{j=1}{\overset{k}{⨁}}{\mathcal{T}}_0^2\left(m_1;p_{j+2,k}\left(m_{j+1}-1\right),d_j=\left(n_1+n_2-2\right)p_{j+2,k}\left(m_{j+1}-1\right)\ge m_1\right).$$

The only condition needed on the sample size in order to have the above convolution of k independent LH-trace distribution is $n_1+n_2\ge m_1+2$. However, LH-trace distributions are usually approximated by $\mathcal{F}$ distribution (We use here the second approximation suggested in McKeon [27]). This means that the distribution of ${\left(D^{\mathrm{pl}}\right)}^2$ can be aproximated by the convolution of k$\mathcal{F}$ distributions, where its critical values are obtained by the method sugested by Dyer [28].

7. An Example

We apply our proposed extended $D^2$ test statistic to a third-order ($k=3$) medical dataset as described in the Introduction, where the interest is in testing the equality of mean of a population of glaucoma patients to a target mean of another population of glaucoma patients [29]. Several studies showed that the central corneal thickness (CCT) plays a major role in the diagnosis of glaucoma. Intraocular pressure (IOP) is positively correlated with CCT and may therefore affect diagnosis. Therefore, CCT should be measured along with IOP in all patients for verification of glaucoma. CCT and IOP vary from individual to individual, from right eye to left eye, and from time to time. We have a sample of 30 glaucoma patients. Measurements on IOP and CCT were taken from both the eyes (sites) and were observed over three time points at an interval of three months. Clearly, then, this dataset is a third-order dataset with $m_1=2,m_2=2$, and $m_3=3$. This dataset was studied by Leiva and Roy [20] by assuming a 3-SSCS covariance structure. Here, we also assume that this dataset has a 3-SSCS covariance structure. The $(2\times 1)$-dimensional sample partitioned mean vector in our sample of 30 glaucoma patients is presented in

Table 1. The (2 × 1) dimensional sample partitioned mean vector in our sample of 30 glaucoma patients.

t	s	(IOP, CCT)
1	1	${(24.333,527.367)}^{\prime }$
1	2	${(23.567,534.633)}^{\prime }$
2	1	${(20.233,525.333)}^{\prime }$
2	2	${(19.567,532.500)}^{\prime }$
3	1	${(19.233,527.133)}^{\prime }$
3	2	${(18.933,534.867)}^{\prime }$

.

Additionally, using the Formulas (27)–(29) presented in Section 5.2, the unbiased estimates ${\widehat{\mathit{U}}}_{3,1}$, ${\widehat{\mathit{U}}}_{3,2}$, and ${\widehat{\mathit{U}}}_{3,3}$ are:

$${\widehat{\mathit{U}}}_{3,1}=\left[\begin{array}{cc}\hfill 12.230& \hfill 12.061\\ \hfill 12.061& \hfill 426.155\end{array}\right],{\widehat{\mathit{U}}}_{3,2}=\left[\begin{array}{cc}\hfill 5.826& \hfill 6.939\\ \hfill 6.939& \hfill 164.156\end{array}\right],\mathrm{and}{\widehat{\mathit{U}}}_{3,3}=\left[\begin{array}{cc}\hfill 3.528& \hfill 9.268\\ \hfill 9.268& \hfill 288.684\end{array}\right],$$

respectively. Using the above estimates, the unbiased estimate of ${\mathsf{\Gamma }}_3$ is:

The $2\times 2$ block diagonals ${\widehat{\mathit{U}}}_{3,1}$ represent the estimate of the variance-covariance matrix of the two response variables IOP and CCT at any given eye and at any given time point, whereas the $2\times 2$ block off diagonals ${\widehat{\mathit{U}}}_{3,2}$ represent the estimate of the covariance matrix of the two response variables IOP and CCT between the two eyes and at any given time point. The $2\times 2$ block off diagonals ${\widehat{\mathit{U}}}_{3,3}$ represent the covariance matrix of the two response variables IOP and CCT between any two time points.

Iester et al. [29] reported the mean and standard deviation (SD) of the IOP and CCT measurements for both the eyes from 794 Italian Caucasian glaucoma patients (see

Table 2. IOP and CCT measurements from 794 Italian Caucasian glaucoma patients.

	Mean	SD
Right IOP	16.16	3.55
Right CCT	545.68	35.82
Left IOP	16.28	3.31
Left CCT	546.89	36.09

). We deem these means as the means of the IOP and CCT at the first time point and then randomly generate four samples within three standard errors (SD of mean) from these reported means of IOP and CCT to represent the means of IOP and CCT for the left and right eyes in the third and sixth months, respectively. These randomly generated means of IOP and CCT for the left and right eyes at three time points in vector form are reported in

Table 3. The (2 × 1) dimensional targeted partitioned mean vector ${\mathbf{\mu }}_0$ in the Italian Caucasian glaucoma patients.

t	s	(IOP, CCT)
1	1	(16.16, 545.68)
1	2	(16.28, 546.89)
2	1	(15.97, 546.18)
2	2	(16.25, 550.30)
3	1	(16.20, 546.90)
3	2	(16.07, 549.64)

, and we will take this mean vector as the targeted mean ${\mathbf{\mu }}_0$ in (21). The sample mean vector in Table 1 appears to be very different from the targeted population mean vector ${\mathbf{\mu }}_0$ in Table 3.

The aim of our study is to see whether our sample of 30 glaucoma patients has the same mean vector as the Italian Caucasian glaucoma patients. Our main intention of the analysis of our glaucoma dataset is to illustrate the use of our new hypotheses testing procedures rather than giving any insight into the dataset itself.

The calculated $D^2$ statistic (26), which is a convolution of three independent L–H distributions, $T_0^2(2;3,87)$, $T_0^2(2;2,58)\mathrm{and}{T}_0^2(2;1,29)$, respectively, which in turn is approximated by two approximated F distributions and one exact F distribution, is 317.2971, and the corresponding p value is 0. So, we reject the null hypothesis that the population mean of our dataset is equal to the Italian Caucasian glaucoma patients, and this conclusion was expected from the data.

8. Conclusions and Discussion

We study the tests of hypotheses of equality of means for one population as well as for two populations for high-dimensional and higher-order data with k-SSCS covariance structure. Such a structure is natural and a credible assumption in many research studies. MLEs and the unbiased estimates of the matrix parameters of the k-SSCS covariance structure have closed-form solutions. On the other hand, the MLEs and the unbiased estimates of the matrix parameters of the $k-$separable covariance structure are not tractable and are computationally intensive. So, k-SSCS covariance structure is a desirable covariance structure for k-th order data. Aghaian et al. [30] examined differences in CCT of 801 subjects, establishing the fact that the CCT of Japanese participants was significantly lower than that of Caucasians, Chinese, Filipinos, and Hispanics, and greater than that of African Americans. African American individuals have thinner corneas compared to white individuals [31]. So, CCT and IOP in glaucoma patients vary with race, and our result confirms this fact. Our proposed new hypotheses testing procedures are perfect for high-dimensional array-variate data, which are ubiquitous in this century. In discriminant analysis [32], the first step is to test the equality of means for the two populations. Therefore, our new method developed in this article will have important applications in the analysis of modern multivariate datasets with higher-order structure. Our new method can be extended to non-normal datasets. In addition, it can be extended in testing the equality of means for more than two populations and simultaneous hypotheses testing in models with k-SSCS covariance structure.