Posterior propriety of an objective prior for generalized hierarchical normal linear models

ABSTRACT

Bayesian Hierarchical models has been widely used in modern statistical application. To deal with the data having complex structures, we propose a generalized hierarchical normal linear (GHNL) model which accommodates arbitrarily many levels, usual design matrices and ‘vanilla’ covariance matrices. Objective hyperpriors can be employed for the GHNL model to express ignorance or match frequentist properties, yet the common objective Bayesian approaches are infeasible or fraught with danger in hierarchical modelling. To tackle this issue, [Berger, J., Sun, D., & Song, C. (2020b). An objective prior for hyperparameters in normal hierarchical models. Journal of Multivariate Analysis, 178, 104606. https://doi.org/10.1016/j.jmva.2020.104606] proposed a particular objective prior and investigated its properties comprehensively. Posterior propriety is important for the choice of priors to guarantee the convergence of MCMC samplers. James Berger conjectured that the resulting posterior is proper for a hierarchical normal model with arbitrarily many levels, a rigorous proof of which was not given, however. In this paper, we complete this story and provide an user-friendly guidance. One main contribution of this paper is to propose a new technique for deriving an elaborate upper bound on the integrated likelihood, but also one unified approach to checking the posterior propriety for linear models. An efficient Gibbs sampling method is also introduced and outperforms other sampling approaches considerably.

KEYWORDS:

• Hierarchical linear model
• linear mixed-effect model
• objective Bayesian analysis
• posterior propriety
• Gibbs sampling

1. Introduction

Bayesian hierarchical models (or multilevel models) have been extensively used in the modern application, including education (Raudenbush & Bryk, 1986), psychology (Lindenberger & Pötter, 1998), clinical trials (Xia et al., 2011), economics (Shimotsu, 2010) and many other applied statistical fields. The fundamental idea of hierarchical modelling is to think of the lowest-level units (smallest and most numerous) as organized into a hierarchy of successively higher-level units. For example, students are in classes, classes are in schools, schools are in districts, and districts are in states. Accordingly, hierarchical models are naturally applicable to the survey, observational or experimental data involved with complicated nesting. However, the most commonly used and fully discussed hierarchical models are merely of two levels. Goldstein (2011) and Berger et al. (2020b) have ever defined 3-level hierarchical model and implemented statistical analysis on that. The hierarchical model with more levels was usually avoided by the researchers for the reason of analytical difficulty and intractable computation. To the best of authors' knowledge, a general hierarchical linear model with arbitrary levels seems to have never been defined or studied. In this paper, we will introduce the definition of a generalized hierarchical normal linear (GHNL) model and carry out an in-depth theoretical investigation of Bayesian inference for the GHNL models.

In order to implement fully Bayesian analysis, priors are supposed to be specified on the hyperparameters (parameters at higher levels of the hierarchical model). Improper (objective) priors are often used to express ignorance or to match frequentist properties (see the review article, Consonni et al., 2018). When using improper priors, an important issue whether the resulting posterior distributions are proper arises. As Hobert and Casella (1996) stated, without proper precaution, misuse of improper priors, sometimes unknowingly, will result in practical difficulties, such as the non-convergence of the Gibbs sampler. The enormous practical importance of posterior propriety motivates us to explore it in the framework of GHNL modelling afterwards. There is also a vast modern literature investigating the posterior propriety of improper priors applied to a large variety of models, such as, Sun et al. (2001), Speckman and Sun (2003), Berger et al. (2005) and Michalak and Morris (2016).

A great deal of efforts have been devoted to the development of objective hyperpriors in hierarchical modelling, such as, Daniels and Kass (1999), Everson and Morris (2000), Gelman (2006), Gustafson et al. (2006), Berger et al. (2005) and Berger et al. (2020b). Formal objective Bayesian approaches, like the Jeffreys-rule prior or reference prior are only feasible for the simple hierarchical settings. For instance, the exact Jeffreys-rule prior for covariance matrices at higher level depends on the parameters from the lower level of the model, leading to plenty of difficulties in formulation and computation. Therefore, a common way is to use less formal approaches, such as applying formal objective priors from non-hierarchical models to hierarchical modelling. Unfortunately, the non-hierarchical Jefferys-rule prior and reference prior typically yield improper posteriors in the hierarchical settings (cf. Berger et al., 2005). Those who can recognize this problem often use constant priors instead for higher level variance components. However, the constant prior is so diffuse that it requires twice as many observations as logically needed to achieve posterior propriety (cf. Berger et al., 2005 and Berger et al., 2020b). In other words, the extra observations required are wasted on correcting the over-diffuse tail of the constant prior. The most powerful tool known for detecting over-diffuse hyperpriors is by looking at the frequentist notion of admissibility of resulting estimators (see Berger et al., 2005 for discussions and references). Sensible choices of objective hyperpriors are on the boundary of admissibility, being as diffuse as possible without leading to inadmissible estimators.

Berger et al. (2005) studied the propriety and admissibility of a number of hyperprios, but no overall conclusion was reached as to a specific prior to recommend. The reasons are as follows: (a) the admissibility of the leading candidate prior was unable to get proved; (b) the proposed computation methods were only efficient for relatively low-dimensional covariance matrices and remained quite challenging for the candidate priors; (c) the hierarchical model discussed was merely of two levels, and the results are not adaptive to a general hierarchical model with many levels. To address this issue, Berger et al. (2020b) recommended a particular objective prior for use in all normal hierarchical models. Consider the following canonical form of 2-level hierarchical normal model. Suppose that, independently, ${\mathit{y}}_i\sim N_k({\mathit{\theta }}_i,{\mathit{I}}_k)$ and ${\mathit{\theta }}_i\sim N_k(\text{β},\mathit{V})$ for $i=1,\dots ,m$, where $N_k(\cdot ,\cdot )$ denotes the k-dimensional normal distribution, ${\mathit{y}}_i$ are $k\times 1$ observation vectors, ${\mathit{\theta }}_i$ are the $k\times 1$ unobserved mean vectors, $\text{β}$ is a $p\times 1$ ‘hypermean’ vector, and $\mathit{V}\in {\mathbb{R}}^{k\times k}$ is an unknown ‘hypercovariance’ matrix. Berger et al. (2020b) proposed a particular combination of independent priors on hyperparameters $\text{β}$ and $\mathit{V}$ as (1) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(k-1)/2}},\\ \pi \left(\mathit{V}\right)& \propto \frac{1}{|\mathit{V}{|}^{1-1/\left(2k\right)}\prod _{1\le s<t\le k}(v_s-v_t)},\end{array}$(1) where $v_1>v_2>\cdots >v_k>0$ are the ordered eigenvalues of $\mathit{V}$. The recommendation (1(1) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(k-1)/2}},\\ \pi \left(\mathit{V}\right)& \propto \frac{1}{|\mathit{V}{|}^{1-1/\left(2k\right)}\prod _{1\le s<t\le k}(v_s-v_t)},\end{array}$(1) ) for hyperpriors was justified by Berger et al. (2020b) from the aspects of admissibility, ease of computation and performance. Most importantly, prior (1(1) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(k-1)/2}},\\ \pi \left(\mathit{V}\right)& \propto \frac{1}{|\mathit{V}{|}^{1-1/\left(2k\right)}\prod _{1\le s<t\le k}(v_s-v_t)},\end{array}$(1) ) is adapted to being used at any level in hierarchical modelling, which is not true for other proposed objective priors as previously mentioned.

Since it is hazardous to skip demonstrating propriety at the risk of making inference from improper posterior, Berger et al. (2020b) has shown the posterior propriety of a 3-level hierarchical model using prior (1(1) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(k-1)/2}},\\ \pi \left(\mathit{V}\right)& \propto \frac{1}{|\mathit{V}{|}^{1-1/\left(2k\right)}\prod _{1\le s<t\le k}(v_s-v_t)},\end{array}$(1) ), while assuming square design matrices for a technical reason. Berger et al. (2020b) also conjectured that the posterior is proper with the recommended prior being utilized at all levels of a hierarchical normal model with arbitrarily many levels, a rigorous proof of which was unable to be provided, however. In this paper, we complete this story and prove the posterior propriety for the GHNL models in general situations. Besides, as pointed out in Michalak and Morris (2016), researchers have been finding it daunting and time-consuming to inspect posterior propriety when using improper priors, except in the simplest models. For this reason, we supply a user-friendly guidance for checking posterior propriety to practitioners in different practical situations.

In Section 2, we give an explicit definition of the GHNL model which accommodates arbitrarily many levels and usual design matrices. It is important to note that we are considering the ‘vanilla’ covariance matrix problem herein. We are not assuming any special structures or sparsity for hypercovariance matrices. The association between the GHNL model and a linear mixed-effect model is also discussed. In Section 3, we demonstrate that recommended prior yields a proper posterior in the framework of GHNL modelling. In addition, we exhibit a guidance for checking posterior propriety. An efficient MCMC algorithm for sampling from the posterior is introduced in Section 4. Section 5 provides some concluding remarks and further generalizations.

2. Generalized hierarchical normal linear model

In this section, we will introduce the definition of a GHNL model with $(r+1)$ levels, where $r\ge 1$. The association between the GHNL model and a linear mixed-effect model will be demonstrated as well, which brings an insight into the GHNL model. At last, the recommended prior on the hyperparameters of the GHNL model will be presented and discussed. Firstly, we introduce some notations to be used in the main body of this paper.

Notations Let $\left[k\right]=\{1,2,\dots ,k\}$ for a positive integer k; $1_{\{\cdot \}}$ stands for the indicator function; $N_k(\mathit{\mu },\mathbf{\Sigma })$ represents the k-dimensional normal distribution with mean $\mathit{\mu }$ and covariance $\mathbf{\Sigma }$; ${\mathcal{N}}_k(\mathit{\mu },\mathbf{\Sigma })$ denotes a k-dimensional normal random variable with mean $\mathit{\mu }$ and covariance $\mathbf{\Sigma }$; for a symmetric matrix $\mathit{A}$, $\mathit{A}>(<)0$ means that $\mathit{A}$ is a positive (negative) definite matrix, and $\mathit{A}\ge (\le )0$ denotes that $\mathit{A}$ is a non-negative (non-positive) definite matrix.

2.1. Model structure

Berger et al. (2020b) proposed a 3-level hierarchical model with the form: (2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) where ${\mathit{y}}_i$ are $k\times 1$ observation vectors, ${\mathit{\theta }}_i$ are the $k\times 1$ unobserved mean vectors, $\mathit{\eta }$ is a $p\times 1$ ‘hypermean’ vector, $\mathit{V}\in {\mathbb{R}}^{k\times k}$ and $\mathit{W}\in {\mathbb{R}}^{p\times p}$ are unknown ‘hypercovariance’ matrices, and ${\mathit{Z}}_i$ are $k\times sp$ known matrices. At last, all the normal random variables in model (2(2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) ) are mutually independent. Based on the 3-level hierarchical normal model, a more general hierarchical model with $(r+1)$ levels ($r\ge 1$) can be constructed as (3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) Firstly, all the normal random variables noted in the above model are mutually independent. Within model (3(3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) ), the output of level $(j+1)$ consists of $m_j$ units whose values are $k_j\times 1$ vectors for $j=0,1,\dots ,r$. By stacking the output units of level $(j+1)$ on top of one another, we can obtain the outcome vector of level $(j+1)$ as ${\mathit{\theta }}_j$ for $j\in \left[r\right]$ and $\mathit{y}={({\mathit{y}}_1^{\mathrm{\top }},\dots ,{\mathit{y}}_{m_0}^{\mathrm{\top }})}^{\mathrm{\top }}$ for level 1. Then ${\mathit{\theta }}_j$'s are $\left(m_j{k}_j\right)\times 1$ vectors and $\mathit{y}$ is a $\left(m_0{k}_0\right)\times 1$ vector. In fact, only the outcome of the lowest level can be observed, and the outcomes of higher levels are inaccessible and latent variables. Hence, the outcome variables of interest are always situated at the lowest level of the hierarchy. Different units in the same level share in common input effects (intercept can be included) which are exactly the outcome vectors from the upper level, except that the input effect of level $(r+1)$ is $\mathit{\eta }$ which is a $d\times 1$ vector of fixed effects. In addition, the units from the same level have the same variance component. The variance component within level $(j+1)$ is denoted by ${\mathit{V}}_j\in {\mathbb{R}}^{k_j\times k_j}$ for $j\in \left[r\right]$ and accounts for the magnitude of random variation within the corresponding level. The covariance matrices ${\mathit{V}}_j$ are unobserved for $j\in \left[r\right]$. The matrices ${\mathit{Z}}_{j{i}_j}$ are $k_j\times \left(m_{j+1}{k}_{j+1}\right)$ matrices and denote the matrices of observed covariates for unit $i_j$ in level j, where $j=0,1,\dots ,r$ and $i_j\in \left[m_j\right]$. It is natural to assume that there exist at least two units in each level and the dimensions of all units and $\mathit{\eta }$ are no less than 1, mathematically, $m_j\ge 2$ and $k_j\ge 1$ for $j=0,1,\dots ,r$, and $d\ge 1$. Table 1 summarizes several important notations that will mainly affect the results for the posterior propriety in Section 3.

Table 1. Summary of certain important notations within model (3(3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) ) and $j=0,1,\dots ,r$.

Notation	Meaning
r + 1	Total number of levels
$m_j$	Number of the output units in level $(j+1)$
$k_j$	Dimension of each output unit in level $(j+1)$
d	Dimension of the fixed effect $\mathit{\eta }$

The extensions from Berger et al. (2020b)'s model (2(2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) ) to model (3(3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) ) are two-fold, and model (3(3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) ) accommodates arbitrarily many levels and usual covariate matrices. Further define ${\mathit{Z}}_j={\{{\mathit{Z}}_{j1}^{\mathrm{\top }},\dots ,{\mathit{Z}}_{j{m}_j}^{\mathrm{\top }}\}}^{\mathrm{\top }}$ for $j=0,1,\dots ,r$. Then ${\mathit{Z}}_j$ are $\left(m_j{k}_j\right)\times \left(m_{j+1}{k}_{j+1}\right)$ matrices for $j=0,1,\dots ,(r-1)$, ${\mathit{Z}}_r$ is an $\left(m_r{k}_r\right)\times d$ matrix, and an alternative representation of the $(r+1)$-level hierarchical model (3(3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) ) is thereby given by (4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4)

Remark 2.1

If we assume that the covariance matrix for the units from level 1 in model (3(3) $\{\begin{array}{lll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \begin{array}{l}{\mathit{y}}_{i_0}={\mathit{Z}}_{0i_0}{\mathit{\theta }}_1\\ +{\mathcal{N}}_{k_0}(\mathbf{0},{\mathit{I}}_{k_0}),\end{array}& \begin{array}{l}{i}_0\in \left[m_0\right],\\ {\mathit{\theta }}_1^{\mathrm{\top }}=({\mathit{\theta }}_{11}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{1m_1}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& \begin{array}{l}{\mathit{\theta }}_{1i_1}={\mathit{Z}}_{1i_1}{\mathit{\theta }}_2\\ +{\mathcal{N}}_{k_1}(\mathbf{0},{\mathit{V}}_1),\end{array}& \begin{array}{l}{i}_1\in \left[m_1\right],{\mathit{\theta }}_2^{\mathrm{\top }}\\ =({\mathit{\theta }}_{21}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{2m_2}^{\mathrm{\top }});\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& \begin{array}{l}{\mathit{\theta }}_{r-1,i_{r-1}}\\ ={\mathit{Z}}_{r-1,i_{r-1}}{\mathit{\theta }}_r\\ +{\mathcal{N}}_{k_{r-1}}(\mathbf{0},{\mathit{V}}_{r-1}),\end{array}& \begin{array}{l}{i}_{r-1}\in \left[m_{r-1}\right],\\ {\mathit{\theta }}_r^{\mathrm{\top }}=({\mathit{\theta }}_{r1}^{\mathrm{\top }},\\ \dots ,{\mathit{\theta }}_{r{m}_r}^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& \begin{array}{l}{\mathit{\theta }}_{r{i}_r}={\mathit{Z}}_{r{i}_r}\mathit{\eta }\\ +{\mathcal{N}}_{k_r}(\mathbf{0},{\mathit{V}}_r),\end{array}& i_r\in \left[m_r\right].\end{array}$(3) ) is a positive definite matrix ${\mathbf{\Sigma }}_0$ instead of the identity matrix, when ${\mathbf{\Sigma }}_0$ is known, the two assumptions are actually equivalent by taking reparameterization that ${\mathit{y}}_{i_0}^{\ast }={\mathbf{\Sigma }}_0^{-\frac{1}{2}}{\mathit{y}}_{i_0}$ and ${\mathit{Z}}_{i_0}^{\ast }={\mathbf{\Sigma }}_0^{-\frac{1}{2}}{\mathit{Z}}_{i_0}$. Furthermore, for a technical reason, ${\mathbf{\Sigma }}_0$ must be assumed as known throughout this paper, and this reason will be explained in Section 5.

2.2. Connection with the linear mixed-effect model $($ LMM $)$

The two-level hierarchical normal models are often referred to as LMMs in many places. As for the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ), let $\mathbf{\Theta }=\{{\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_r\}$ denote the set of unobserved outcome vectors and $\mathcal{V}=\{{\mathit{V}}_1,\dots ,{\mathit{V}}_r\}$ represent the set of unknown covariance matrices. If we take ${\mathit{\theta }}_j$'s as intermediate variables, then marginalizing out over $\mathbf{\Theta }$ yields (5) $(\mathit{y}|\mathit{\eta },\mathcal{V})\sim N_{m_0{k}_0}({\mathit{X}}_r\mathit{\eta },\mathbf{\Delta }),$(5) where (6) $\begin{array}{rl}\mathbf{\Delta }& ={\mathit{I}}_{m_0{k}_0}+\sum _{t=1}^r{\mathit{X}}_{t-1}({\mathit{I}}_{m_t}\otimes {\mathit{V}}_t){\mathit{X}}_{t-1}^{\mathrm{\top }},\mathrm{a}\mathrm{n}\mathrm{d}\\ {\mathit{X}}_j& =\prod _{s=0}^j{\mathit{Z}}_s,j=0,1,\dots ,r.\end{array}$(6) $\mathbf{\Delta }$ is a $\left(m_0{k}_0\right)\times \left(m_0{k}_0\right)$ matrix and ${\mathit{X}}_j$ are $\left(m_0{k}_0\right)\times \left(m_{j+1}{k}_{j+1}\right)$ matrices for $j=0,1,\dots ,r$. Suppose that ${\mathit{Z}}_j$ are of full column ranks. Then by Sylvester's rank inequality ${\mathit{X}}_j$ are also of full column ranks, $j=0,1,\dots ,r$. In the rest of the paper, ${\mathit{Z}}_j$ are assumed to be of full column ranks for $j=0,1,\dots ,r$.

If we consider a particular LMM as (7) $\mathit{y}={\mathit{X}}_r\mathit{\eta }+{\mathit{X}}_0{\mathit{\theta }}_1^{\ast }+\cdots +{\mathit{X}}_{r-1}{\mathit{\theta }}_r^{\ast }+\mathit{ϵ},$(7) where $\mathit{\eta }$ is the fixed effect, ${\mathit{\theta }}_j$'s are random effects and independently distributed as $N_{m_j{k}_j}(\mathbf{0},{\mathit{I}}_{m_j}\otimes {\mathit{V}}_j)$ for $j\in \left[r\right]$, and $\mathit{ϵ}$ denotes the vector of random errors and is distributed as $N_{m_0{k}_0}(\mathbf{0},{\mathit{I}}_{m_0{k}_0})$. By integrating out the random effects, the marginal distribution of $\mathit{y}$ conditioning on $(\mathit{\eta },\mathcal{V})$ is identical to the distribution (5(5) $(\mathit{y}|\mathit{\eta },\mathcal{V})\sim N_{m_0{k}_0}({\mathit{X}}_r\mathit{\eta },\mathbf{\Delta }),$(5) ). In one word, the GHNL is equivalent to an LMM in the sense of the marginal distribution of observations after integrating out intermediate outcome vectors or random effects. The equivalence between the GHNL models and the LMMs can be illustrated by an example of mixed-effect ANOVA model as well.

Example 2.1

Mixed-effect ANOVA model

Suppose we can observe the scores of p courses for student (ijk) as ${\mathit{y}}_{ijk}$ for $i=1,\dots ,s_1$, $j=1,\dots ,s_2$ and $k=1,\dots ,s_3$. The observed data are within a hierarchy of three levels: student (ijk) is nested within class (ij), and class (ij) is nested within school i. Thus, we have total $s_1$ schools, each school has $s_2$ classes and each class has $s_3$ students. Consider a mixed-effect ANOVA model as (8) ${\mathit{y}}_{ijk}=\mathit{\eta }+{\mathit{\alpha }}_i+{\text{β}}_{ij}+{\mathit{ϵ}}_{ijk},$(8) where ${\mathit{y}}_{ijk}$, $\mathit{\eta },{\alpha }_i,{\text{β}}_{ij}$ and ${\mathit{ϵ}}_{ijk}$ are all $p\times 1$ vectors for $i=1,\dots ,s_1$, $j=1,\dots ,s_2$ and $k=1,\dots ,s_3$, $\mathit{\eta }$ denotes the overall mean and is fixed effect, ${\mathit{\alpha }}_i\sim N_p(\mathbf{0},{\mathit{V}}_{\alpha })$ is the effect of school, ${\text{β}}_{ij}$ is distributed as $N_p(\mathbf{0},{\mathit{V}}_{\text{β}})$ and represents the effect of class, the student-level independent random error is denoted by ${\mathit{ϵ}}_{ijk}$ and has distribution $N_p(\mathbf{0},{\mathbf{\Sigma }}_0)$, and ${\mathbf{\Sigma }}_0$ is a known matrix. At the same time, ${\mathit{\alpha }}_i$, ${\text{β}}_{ij}$ and ${\mathit{ϵ}}_{ijk}$ are independently distributed. Consequently, ${\mathit{V}}_{\alpha }$, ${\mathit{V}}_{\text{β}}$ and ${\mathbf{\Sigma }}_0$ are the variance components describing the school-level, class-level and student-level variations, respectively. Due to the hierarchical structure of the observations, we can naturally build a hierarchical model as (9) $\begin{array}{rl}{\mathit{y}}_{ijk}& \sim N_p({\text{β}}_{ij}^{\ast },{\mathbf{\Sigma }}_0),{\text{β}}_{ij}^{\ast }\sim N_p({\mathit{\alpha }}_i^{\ast },{\mathit{V}}_{\text{β}})\mathrm{a}\mathrm{n}\mathrm{d}\\ {\mathit{\alpha }}_i^{\ast }& \sim N_p(\mathit{\eta },{\mathit{V}}_{\mathit{\alpha }})\end{array}$(9) independently, for $i=1,\dots ,s_1$, $j=1,\dots ,s_2$ and $k=1,\dots ,s_3$. Denote that $\begin{array}{rl}{\mathit{Y}}_i& =\left(\begin{array}{ccc}{\mathit{y}}_{i11}& \cdots & {\mathit{y}}_{i{s}_{2}1}\\ ⋮& \ddots & ⋮\\ {\mathit{y}}_{i1s_3}& \cdots & {\mathit{y}}_{i{s}_2{s}_3}\end{array}\right),\\ {\mathit{E}}_i& =\left(\begin{array}{ccc}{\mathit{ϵ}}_{i11}& \cdots & {\mathit{ϵ}}_{i{s}_{2}1}\\ ⋮& \ddots & ⋮\\ {\mathit{ϵ}}_{i1s_3}& \cdots & {\mathit{ϵ}}_{i{s}_2{s}_3}\end{array}\right),\\ {\text{β}}_i& =\left(\begin{array}{c}{\text{β}}_{i1}\\ ⋮\\ {\text{β}}_{i{s}_2}\end{array}\right)\mathrm{a}\mathrm{n}\mathrm{d}\\ {\text{β}}_i^{\ast }& =\left(\begin{array}{c}{\text{β}}_{i1}^{\ast }\\ ⋮\\ {\text{β}}_{i{s}_2}^{\ast }\end{array}\right),\end{array}$where ${\mathit{Y}}_i$ and ${\mathit{E}}_i$ are both $\left(s_{3}p\right)\times s_2$ matrices, and ${\text{β}}_i$ and ${\text{β}}_i^{\ast }$ are both $\left(s_{2}p\right)\times 1$ vectors. Let ${\mathit{y}}_i=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{Y}}_i\right)$ and ${\mathit{ϵ}}_i=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{E}}_i\right)$, where $\mathrm{v}\mathrm{e}\mathrm{c}\left(\mathit{A}\right)$ denotes the column vector obtained by stacking the columns of the matrix $\mathit{A}$ on top of one another. Define that $\begin{array}{rl}\mathit{y}& =\left(\begin{array}{c}{\mathit{y}}_1\\ ⋮\\ {\mathit{y}}_{s_1}\end{array}\right),\mathit{ϵ}=\left(\begin{array}{c}{\mathit{ϵ}}_1\\ ⋮\\ {\mathit{ϵ}}_{s_1}\end{array}\right),\text{β}=\left(\begin{array}{c}{\text{β}}_1\\ ⋮\\ {\text{β}}_{s_1}\end{array}\right),\\ {\text{β}}^{\ast }& =\left(\begin{array}{c}{\text{β}}_1^{\ast }\\ ⋮\\ {\text{β}}_{s_1}^{\ast }\end{array}\right),\mathit{\alpha }=\left(\begin{array}{c}{\mathit{\alpha }}_1\\ ⋮\\ {\mathit{\alpha }}_{s_1}\end{array}\right),{\mathit{\alpha }}^{\ast }=\left(\begin{array}{c}{\mathit{\alpha }}_1^{\ast }\\ ⋮\\ {\mathit{\alpha }}_{s_1}^{\ast }\end{array}\right).\end{array}$Thus, $\mathit{y}$ and $\mathit{ϵ}$ are $\left(m_{0}p\right)\times 1$ vectors, and $\text{β}$ and ${\text{β}}^{\ast }$ are $\left(m_{1}p\right)\times 1$ vectors, and $\mathit{\alpha }$ and ${\mathit{\alpha }}^{\ast }$ are $\left(m_{2}p\right)\times 1$ vectors, where $m_0=s_1{s}_2{s}_3$, $m_1=s_1{s}_2$ and $m_2=s_1$. Then the hierarchical normal model (9(9) $\begin{array}{rl}{\mathit{y}}_{ijk}& \sim N_p({\text{β}}_{ij}^{\ast },{\mathbf{\Sigma }}_0),{\text{β}}_{ij}^{\ast }\sim N_p({\mathit{\alpha }}_i^{\ast },{\mathit{V}}_{\text{β}})\mathrm{a}\mathrm{n}\mathrm{d}\\ {\mathit{\alpha }}_i^{\ast }& \sim N_p(\mathit{\eta },{\mathit{V}}_{\mathit{\alpha }})\end{array}$(9) ) can be expressed as a GHNL model with the form (10) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& (\mathit{y}|{\text{β}}^{\ast },{\mathbf{\Sigma }}_0)& \sim & \begin{array}{l}{N}_{m_{0}p}({\mathit{Z}}_0{\text{β}}^{\ast },{\mathit{I}}_{m_0}\otimes {\mathbf{\Sigma }}_0);\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\text{β}}^{\ast }|{\mathit{\alpha }}^{\ast },{\mathit{V}}_{\text{β}})& \sim & \begin{array}{l}{N}_{m_{1}p}({\mathit{Z}}_1{\mathit{\alpha }}^{\ast },{\mathit{I}}_{m_1}\otimes {\mathit{V}}_{\text{β}});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:& ({\mathit{\alpha }}^{\ast }|\mathit{\eta },{\mathit{V}}_{\alpha })& \sim & \begin{array}{l}{N}_{m_{2}p}({\mathit{Z}}_2\mathit{\eta },{\mathit{I}}_{m_2}\otimes {\mathit{V}}_{\alpha }),\end{array}\end{array}$(10) where $\begin{array}{rl}{\mathit{Z}}_0& =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\underset{s_1{s}_2}{\underbrace{{\mathbf{1}}_{s_3}\otimes {\mathit{I}}_p,\dots ,{\mathbf{1}}_{s_3}\otimes {\mathit{I}}_p}}\right\},\\ {\mathit{Z}}_1& =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\underset{s_1}{\underbrace{{\mathbf{1}}_{s_2}\otimes {\mathit{I}}_p,\dots ,{\mathbf{1}}_{s_2}\otimes {\mathit{I}}_p}}\right\},{\mathit{Z}}_2={\mathbf{1}}_{s_1}\otimes {\mathit{I}}_p,\end{array}$where ${\mathbf{1}}_q$ denotes the $q\times 1$ vector with all elements being one, and ${\mathit{Z}}_0$, ${\mathit{Z}}_1$, ${\mathit{Z}}_2$ are $\left(m_{0}p\right)\times \left(m_{1}p\right)$, $\left(m_{1}p\right)\times \left(m_{2}p\right)$, $\left(m_{2}p\right)\times p$ matrices, respectively. Denote that $\begin{array}{rl}{\mathit{X}}_0& \triangleq {\mathit{Z}}_0,\\ {\mathit{X}}_1& \triangleq \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\underset{s_1}{\underbrace{{\mathbf{1}}_{s_2{s}_3}\otimes {\mathit{I}}_p,\dots ,{\mathbf{1}}_{s_2{s}_3}\otimes {\mathit{I}}_p}}\right\}={\mathit{Z}}_0{\mathit{Z}}_1,\\ {\mathit{X}}_2& \triangleq {\mathbf{1}}_{s_1{s}_2{s}_3}\otimes {\mathit{I}}_p={\mathit{Z}}_0{\mathit{Z}}_1{\mathit{Z}}_2,\end{array}$and ${\mathit{X}}_0$, ${\mathit{X}}_1$, ${\mathit{X}}_2$ are $\left(m_{0}p\right)\times \left(m_{1}p\right)$, $\left(m_{0}p\right)\times \left(m_{2}p\right)$, $\left(m_{0}p\right)\times p$ matrices, respectively. Thus, model (8(8) ${\mathit{y}}_{ijk}=\mathit{\eta }+{\mathit{\alpha }}_i+{\text{β}}_{ij}+{\mathit{ϵ}}_{ijk},$(8) ) can be summarized as (11) $\mathit{y}={\mathit{X}}_2\mathit{\eta }+{\mathit{X}}_1\mathit{\alpha }+{\mathit{X}}_0\text{β}+\mathit{ϵ},$(11) where $\mathit{\alpha }\sim N_{m_{2}p}(\mathbf{0},{\mathit{I}}_{m_2}\otimes {\mathit{V}}_{\alpha })$, $\text{β}\sim N_{m_{1}p}(\mathbf{0},{\mathit{I}}_{m_1}\otimes {\mathit{V}}_{\text{β}})$ and $\mathit{ϵ}\sim N_{m_{0}p}(\mathbf{0},{\mathit{I}}_{m_0}\otimes {\mathbf{\Sigma }}_0)$, independently. By integrating out $(\mathit{\alpha },\text{β})$ and $({\mathit{\alpha }}^{\ast },{\text{β}}^{\ast })$, the marginal distributions of $\mathit{y}$ for model (11(11) $\mathit{y}={\mathit{X}}_2\mathit{\eta }+{\mathit{X}}_1\mathit{\alpha }+{\mathit{X}}_0\text{β}+\mathit{ϵ},$(11) ) and model (10(10) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& (\mathit{y}|{\text{β}}^{\ast },{\mathbf{\Sigma }}_0)& \sim & \begin{array}{l}{N}_{m_{0}p}({\mathit{Z}}_0{\text{β}}^{\ast },{\mathit{I}}_{m_0}\otimes {\mathbf{\Sigma }}_0);\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\text{β}}^{\ast }|{\mathit{\alpha }}^{\ast },{\mathit{V}}_{\text{β}})& \sim & \begin{array}{l}{N}_{m_{1}p}({\mathit{Z}}_1{\mathit{\alpha }}^{\ast },{\mathit{I}}_{m_1}\otimes {\mathit{V}}_{\text{β}});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:& ({\mathit{\alpha }}^{\ast }|\mathit{\eta },{\mathit{V}}_{\alpha })& \sim & \begin{array}{l}{N}_{m_{2}p}({\mathit{Z}}_2\mathit{\eta },{\mathit{I}}_{m_2}\otimes {\mathit{V}}_{\alpha }),\end{array}\end{array}$(10) ) are identical and of the form: $\begin{array}{rl}& (\mathit{y}|\mathit{\eta },{\mathit{V}}_{\alpha },{\mathit{V}}_{\text{β}},{\mathbf{\Sigma }}_0)\sim N_{m_{0}p}({\mathit{X}}_2\mathit{\eta },\mathbf{\Omega }),\\ & \mathbf{\Omega }={\mathit{I}}_{m_0}\otimes {\mathbf{\Sigma }}_0+{\mathit{X}}_0({\mathit{I}}_{m_1}\otimes {\mathit{V}}_{\text{β}}){\mathit{X}}_0^{\mathrm{\top }}\\ & +{\mathit{X}}_1({\mathit{I}}_{m_2}\otimes {\mathit{V}}_{\alpha }){\mathit{X}}_1^{\mathrm{\top }}.\end{array}$

Example 2.1 provides a simple example to illustrate how the hierarchical model and the mixed-effect model can be constructed based on the nested data, and the equivalence between two models is also presented. In Appendix 1, we define a special LMM which is a special case of model (7(7) $\mathit{y}={\mathit{X}}_r\mathit{\eta }+{\mathit{X}}_0{\mathit{\theta }}_1^{\ast }+\cdots +{\mathit{X}}_{r-1}{\mathit{\theta }}_r^{\ast }+\mathit{ϵ},$(7) ) with $m_j=1$ for all $j\in \left[r\right]$, and the theoretical investigation of this special LMM is distinct from that of the GHNLM. This special LMM could be common in application, and we just want to provide some theoretical results for those who have interests to refer to. We will still focus on the GHNL model in the following sections.

2.3. Priors on the hyperparameters

In order to implement fully Bayesian analysis, we should specify hyperpriors on the parameters $(\mathit{\eta },\mathcal{V})$. It follows the recommendation from Berger et al. (2020b) that we can assume priors on $(\mathit{\eta },\mathcal{V})$ as: (12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) (13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) where ${\omega }_{j1}>{\omega }_{j2}>\cdots >{\omega }_{j{k}_j}>0$ are the decreasingly ordered eigenvalues of ${\mathit{V}}_j$, $j\in \left[r\right]$. Apart from prior (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ), common choices of prior on $\mathit{\eta }$ include the constant prior and conjugate prior. None of the three priors will result in improper priors or difficulties in computation. However, among the three priors, prior (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) is the most perfect for all k from the perspective of admissibility. Besides, it refers to Berger et al. (2005) that the prior (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) is a mixture-of-normal prior of the hierarchical structure as (14) $\left(\mathit{\eta }|\lambda \right)\sim N_d(\mathbf{0},\lambda {\mathit{I}}_d)\mathrm{a}\mathrm{n}\mathrm{d}\left[\lambda \right]\propto {\lambda }^{-1/2}\exp (-\frac{1}{2\lambda }),$(14) and those mixture-of-normal priors have shown great success in shrinkage estimation particularly (cf. Fourdrinier et al., 1998) and robust Bayesian estimation generally (cf. Berger, 1980). Therefore, prior (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) was actually recommended by Berger et al. (2005) for default use.

As for prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on the unknown covariance matrices ${\mathit{V}}_j$, $j\in \left[r\right]$, consider the transformation from ${\mathit{V}}_j$ to ${\mathbf{\Omega }}_j=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\omega }_{j1},\dots ,{\omega }_{j{k}_j})$ and the orthogonal matrix ${\mathbf{\Gamma }}_j$ of corresponding eigenvectors, the Jacobian is (15) $\left|\frac{\mathrm{\partial }{\mathit{V}}_j}{\mathrm{\partial }({\mathbf{\Omega }}_j,{\mathbf{\Gamma }}_j)}\right|=\prod _{s<t}({\omega }_{js}-{\omega }_{jt}).$(15) Consequently, the prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on ${\mathit{V}}_j$ becomes the prior density of $({\mathbf{\Omega }}_j,{\mathbf{\Gamma }}_j)$ as (16) $\pi ({\mathbf{\Omega }}_j,{\mathbf{\Gamma }}_j)\propto \frac{1}{|{\mathbf{\Omega }}_j{|}^{1-1/\left(2k_j\right)}}$(16) with respect to Lebesgue measure on $({\omega }_{j1},\dots ,{\omega }_{j{k}_j})$ and the invariant Haar measure over the space $\{\mathbf{\Gamma }:\mathbf{\Gamma }{\mathbf{\Gamma }}^{\mathrm{\top }}={\mathit{I}}_{k_j}\}$. Note that the prior on ${\mathbf{\Omega }}_j$ is improper and, independently, the prior on ${\mathbf{\Gamma }}_j$ is constant. Use of a uniform prior for ${\mathbf{\Gamma }}_j$ ranging over a compact space is natural and non-controversial and has no influence on the eigenvalues. The term $\prod _{s<t}({\omega }_{js}-{\omega }_{jt})$ is eliminated after changing variables for prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ). In contrast, the commonly used priors on the covariance matrix, such as inverse Wishart, Jeffreys-rule and constant priors, contain the term $\prod _{s<t}({\omega }_{js}-{\omega }_{jt})$ in the transformed space. This special term gives low mass to close eigenvalues and hence effectively force the eigenvalues apart. It is contrary to the common intuition which would suggest that one should choose a prior that pushes the eigenvalues closer together. As a result, prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) is essentially neutral as to expansion or shrinkage of the eigenvalues.

In the context of 2-level hierarchical normal model, Theorem 1 from Berger et al. (2020b) has demonstrated that the combination of priors (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) and (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on $(\mathit{\eta },\mathcal{V})$ is on the boundary of admissibility, being as diffuse as possible without yielding inadmissible estimators. Furthermore, it is shown that the generalization that allows covariates at all levels of the hierarchical model will not affect the result of admissibility (cf. Berger et al., 2020b). Nonetheless, the admissibility of the recommended prior for the $(r+1)$-level hierarchical model with $r\ge 2$ is not clear. Generally speaking, this is a very difficult question to answer, and we mainly justify the recommendation of hyperpriors from the angles of posterior propriety and computation afterwards in the framework of the GHNL model.

3. Posterior propriety

Berger et al. (2020b) has shown that the resulting posterior of the recommended prior is proper for the 3-level hierarchical model (2(2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) ), but under a narrow set of assumptions. They also conjectured the posterior propriety for a hierarchical model with any number of levels, a rigorous one of which was not given yet. In this section, we will comprehensively investigate the conditions for the posterior propriety of the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) using the recommended prior in more general situations. The dimension of $\mathit{\eta }$ affects the investigation of posterior propriety considerably, and two cases $d\ge 2$ or d = 1 will be discussed separately.

Based on (5(5) $(\mathit{y}|\mathit{\eta },\mathcal{V})\sim N_{m_0{k}_0}({\mathit{X}}_r\mathit{\eta },\mathbf{\Delta }),$(5) ) and (14(14) $\left(\mathit{\eta }|\lambda \right)\sim N_d(\mathbf{0},\lambda {\mathit{I}}_d)\mathrm{a}\mathrm{n}\mathrm{d}\left[\lambda \right]\propto {\lambda }^{-1/2}\exp (-\frac{1}{2\lambda }),$(14) ), by integrating out $\mathit{\eta }$, we can obtain the marginal distribution of $\mathit{y}$ conditioning on $(\mathcal{V},\lambda )$ as (17) $(\mathit{y}|\mathcal{V},\lambda )\sim N_{m_0{k}_0}(\mathbf{0},\mathbf{\Delta }+\lambda {\mathit{X}}_r{\mathit{X}}_r^{\mathrm{\top }}).$(17) The posterior propriety of the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) employing priors (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) and (14(14) $\left(\mathit{\eta }|\lambda \right)\sim N_d(\mathbf{0},\lambda {\mathit{I}}_d)\mathrm{a}\mathrm{n}\mathrm{d}\left[\lambda \right]\propto {\lambda }^{-1/2}\exp (-\frac{1}{2\lambda }),$(14) ) is defined as (18) $m\left(\mathit{y}\right)=\int f(\mathit{y}|\mathcal{V},\lambda )\pi \left(\lambda \right)\prod _{s=1}^r\pi \left({\mathit{V}}_s\right)\mathrm{d}\lambda \prod _{t=1}^r\mathrm{d}{\mathit{V}}_t<\mathrm{\infty },$(18) where $m\left(\mathit{y}\right)$ denotes the marginal density of the observation vector. Next, we display some definitions and additional notations which are frequently used in this section.

More Notations Let $\text{card}\left(A\right)$ denote the cardinality of the set A. For $0<I_1,I_2\le \mathrm{\infty }$, denote that $I_1\simeq I_2$ if there exist constants $0<C_1\le C_2$ such that $C_1{I}_2\le I_1\le C_2{I}_2$; For a symmetric matrix $\mathit{A}\in {\mathbb{R}}^{n\times n}$, ${\lambda }_i\left(\mathit{A}\right)$ represents the i-th largest eigenvalue of $\mathit{A}$, namely, ${\lambda }_1\left(\mathit{A}\right)\ge \cdots \ge {\lambda }_n\left(\mathit{A}\right)$; Let ${\lambda }_{max}\left(\mathit{A}\right)$ and ${\lambda }_{min}\left(\mathit{A}\right)$ denote the maximum and minimum eigenvalues of an arbitrary symmetric matrix $\mathit{A}$.

Definition

For convenience, let ${\omega }_{r+1,1}\triangleq \lambda$, $m_{r+1}\triangleq d$ and $k_{r+1}\triangleq 1$. Define a function of an arbitrary non-empty set E as $\mathcal{F}\left(E\right)\triangleq \{D|D\subset E,D\ne \mathrm{\varnothing }\}$, such that $\mathcal{F}\left(E\right)$ denotes the set of all non-empty subsets of E. For any $R\in \mathcal{F}\left(\right[r+1\left]\right)$, let $H\left(R\right)\triangleq \left\{\right(j,l)|j\in R,l\in [k_j\left]\right\}$. Further define the composition of $\mathcal{F}$ and H as $\mathcal{S}\left(R\right)\triangleq (\mathcal{F}\circ H)\left(R\right)=\{D|D\subset H(R),D\ne \mathrm{\varnothing }\}$for any $R\in \mathcal{F}\left(\right[r+1\left]\right)$. Define that (19) $c_{jl,s}=1_{\left\{m_j\right(l-1)<s\le m_{j}l\}},$(19) for $j\in [r+1]$, $l\in \left[k_j\right]$ and $s\in \left[m_0{k}_0\right]$.

3.1. Two key lemmas

Before we formally investigate the posterior propriety, we first introduce two important lemmas which dominate in the process of proving the main theorems in this paper.

Lemma 3.1

Assume that ${\mathit{A}}_j$ are $p_j\times p_j$ positive definite matrices, $j\in \left[r\right]$. Let ${\mathit{X}}_j$ be $n\times p_j$ matrices of full column ranks, $j\in \left[r\right]$. Define that (20) $\mathit{H}=\sum _{j=1}^r{\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}.$(20) Then

1. ${\lambda }_{max}\left(\mathit{H}\right)\le C_1\sum _{j=1}^r{\lambda }_{max}\left({\mathit{A}}_j\right)$, where $C_1=\underset{j\in \left[r\right]}{max}{\lambda }_{max}\left({\mathit{X}}_j^{\mathrm{\top }}{\mathit{X}}_j\right)$.

2. Also, (21) $\left|\mathit{H}\right|\ge {\left(\frac{C_2}{r}\right)}^n\left|\sum _{j=1}^r{\mathit{D}}_j\right|,$(21) where $C_2=\underset{j\in \left[r\right]}{min}{\lambda }_{min}\left({\mathit{X}}_j^{\mathrm{\top }}{\mathit{X}}_j\right)>0$. For any $j\in \left[r\right]$, ${\mathit{D}}_j=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_{j1},\dots ,a_{jn})$, where $a_{jk}={\lambda }_k\left({\mathit{A}}_j\right)$ for $k\in \left[p_j\right]$ and $a_{jk}=0$ for $p_j<k\le n$.

Lemma 3.1 mainly demonstrates two inequalities with respect to the summation of a series of quadratic forms which have matrical inputs. Since ${\mathit{X}}_i$'s have full column rank and ${\mathit{A}}_j$'s are positive definite, then $n\ge p_j$ and $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)=p_j$. It is worthwhile to note that the non-zero diagonal elements of ${\mathit{D}}_j$ are the decreasingly ordered eigenvalues of ${\mathit{A}}_j$ in the lower bound of $\left|\mathit{H}\right|$, and this relation will deeply influence the last result when we derive the sufficient condition for posterior propriety afterwards. Besides, we can never find a constant $C^{\ast }>0$ such that ${\left|\mathit{H}\right|}_{+}\le C^{\ast }{\left|\sum _{j=1}^r{\mathit{D}}_j\right|}_{+}$, where ${\left|\mathit{M}\right|}_{+}$ denotes the product of all non-zero eigenvalues of $\mathit{M}$. The proof of Lemma 3.1 can be found in Appendix A.2.

Lemma 3.2

Assume a positive integer k. Suppose $\mathcal{M}$ is a subset of $\mathcal{F}\left(\right[k\left]\right)$ with cardinality n and let $C_1,\dots ,C_n$ denote the sequence of all the elements within $\mathcal{M}$. Define an integral $I={\int }_{\mathrm{\Omega }}\frac{1}{[\prod _{j=1}^k{\lambda }_j^{a_j}][\prod _{i=1}^n(1+\sum _{r\in C_i}{\lambda }_r{)}^{b_i}]}\mathrm{d}\mathit{\lambda },$where $a_{j}j=1,2,\dots ,k,$ and $b_i,i=1,2,\dots ,n,$ are real constants, and $\mathit{\lambda }=({\lambda }_1,\dots ,{\lambda }_k{)}^{\mathrm{\top }}\in \mathrm{\Omega }\equiv [0,\mathrm{\infty }{]}^k$. Then the integral I is finite if and only if (iff) the following two conditions are both satisfied.

• $a_j<1$, $j\in \left[k\right]$;

• inequalities (22) $\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right)$(22) for all $D\in \mathcal{F}\left(\right[k\left]\right)$ hold, where ${\mathcal{G}}_D=\{i|D\bigcap C_i\ne \mathrm{\varnothing },i\in \left[n\right]\}$.

Here, Lemma 3.2 (b) may not be straightforward enough, so we will use the following example to elaborate on how Lemma 3.2 can be employed.

Example 3.1

Consider integral $I_0={\int }_{\mathrm{\Omega }}\frac{1}{\begin{array}{c}[\prod _{j=1}^3{\lambda }_j^{a_j}](1+{\lambda }_1{)}^{b_1}(1+{\lambda }_1+{\lambda }_2{)}^{b_2}(1+{\lambda }_2+{\lambda }_3{)}^{b_3}\end{array}}\mathrm{d}{\lambda }_1\mathrm{d}{\lambda }_2\mathrm{d}{\lambda }_3,$where $a_1,a_2,a_3,b_1,b_2,b_3$ are all real constants. Similar to Lemma 3.2, we can define $\mathcal{M}=\{\left\{1\right\},\{1,2\},\{2,3\}\}$. Then $I_0<\mathrm{\infty }$ iff all the following inequalities hold: (a) $a_j<1,j\in \left[3\right]$; (b) $\begin{array}{rl}& \mathrm{f}\mathrm{o}\mathrm{r}D=\left\{1\right\},a_1+b_1+b_2>1;\\ & \mathrm{f}\mathrm{o}\mathrm{r}D=\left\{2\right\},a_2+b_2+b_3>1;\\ & \mathrm{f}\mathrm{o}\mathrm{r}D=\left\{3\right\},a_3+b_3>1;\\ & \mathrm{f}\mathrm{o}\mathrm{r}D=\{1,2\},a_1+a_2+b_1+b_2+b_3>2;\\ & \mathrm{f}\mathrm{o}\mathrm{r}D=\{1,3\},a_1+a_3+b_1+b_2+b_3>2;\\ & \mathrm{f}\mathrm{o}\mathrm{r}D=\{2,3\},a_2+a_3+b_2+b_3>2;\\ & \mathrm{f}\mathrm{o}\mathrm{r}D=\{1,2,3\},a_1+a_2+a_3+b_1+b_2+b_3>3.\end{array}$

Note that, no matter how $\mathcal{M}$ is defined, we always need to check all the inequalities corresponding to all $D\in \mathcal{F}\left(\right[k\left]\right)$. Even though some inequalities could be trivial after being written down, for the sake of assurance, we would better take all non-empty subsets of $\left[k\right]$ into account in the early stage. Lemma 3.2 plays a crucial role in obtaining the follow-up theorems, and detailed proof of Lemma 3.2 sees Appendix A.2.

3.2. Conditions for the posterior to be proper when $d\ge 2$

In this subsection, the case $d\ge 2$ is mainly considered.

Theorem 3.1

Consider the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with priors (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) and (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on $\mathit{\eta }$ and ${\mathit{V}}_i\in \mathcal{V}$, respectively. When $d\ge 2$, a sufficient condition for the posterior propriety is given by (23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) for any $D\in \mathcal{S}\left(\right[r+1\left]\right)$.

Proof.

It follows (17(17) $(\mathit{y}|\mathcal{V},\lambda )\sim N_{m_0{k}_0}(\mathbf{0},\mathbf{\Delta }+\lambda {\mathit{X}}_r{\mathit{X}}_r^{\mathrm{\top }}).$(17) ) that the integrated likelihood of $(\mathcal{V},\lambda )$ after marginalizing out $({\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_r,\eta )$ is given by (24) $\begin{array}{rl}L(\mathcal{V},\lambda ;\mathit{y})& \propto \frac{1}{|\mathbf{\Delta }+\lambda {\mathit{X}}_r{\mathit{X}}_r^{\mathrm{\top }}{|}^{1/2}}\\ & \times \exp \{-\frac{1}{2}{\mathit{y}}^{\mathrm{\top }}{(\mathbf{\Delta }+\lambda {\mathit{X}}_r{\mathit{X}}_r^{\mathrm{\top }})}^{-1}\mathit{y}\}.\end{array}$(24) By dropping the exponent term involving $\mathit{y}$ (since it is less than one), we have $L(\mathcal{V},\lambda ;\mathit{y})<\frac{1}{|\mathbf{\Delta }+\lambda {\mathit{X}}_r{\mathit{X}}_r^{\mathrm{\top }}{|}^{1/2}}.$By applying Lemma 3.1 (b), we can further bound the integrated likelihood as (25) $L(\mathcal{V},\lambda ;\mathit{y})\le \frac{C_1}{{\left|{\mathit{M}}_1\right|}^{1/2}}$(25) where $C_1$ is a positive constant that only depends on ${\mathit{X}}_j$'s and (26) $\begin{array}{rl}{\mathit{M}}_1& ={\mathit{I}}_{m_0{k}_0}+\sum _{j=1}^{r+1}{\mathit{D}}_j\mathrm{a}\mathrm{n}\mathrm{d}\\ {\mathit{D}}_j& ={\left(\begin{array}{cc}{\mathbf{\Omega }}_j\otimes {\mathit{I}}_{m_j}& \\ & {\mathit{O}}_{q_j}\end{array}\right)}_{\left(m_0{k}_0\right)\times \left(m_0{k}_0\right)},\end{array}$(26) where ${\mathbf{\Omega }}_j$ are the diagonal matrices of the decreasingly ordered eigenvalues of ${\mathit{V}}_j$ for $j\in \left[r\right]$, ${\mathbf{\Omega }}_{r+1}=\left({\omega }_{r+1,1}\right)$, ${\mathit{O}}_{q_j}$ are $q_j\times q_j$ zero matrices and $q_j=m_0{k}_0-m_j{k}_j$ for $j\in [r+1]$. Since ${\mathbf{\Omega }}_{r+1}=\left({\omega }_{r+1,1}\right)$, then ${\mathbf{\Omega }}_{r+1}$ is a degenerate matrix as scalar λ and the prior on λ becomes (27) $\pi \left({\mathbf{\Omega }}_{r+1}\right)\propto \frac{1}{|{\mathbf{\Omega }}_{r+1}{|}^{1-1/\left(2k_{r+1}\right)}}\exp \{-\frac{1}{2}\mathrm{t}\mathrm{r}\left({\mathbf{\Omega }}_{r+1}^{-1}\right)\}.$(27) Combining (16(16) $\pi ({\mathbf{\Omega }}_j,{\mathbf{\Gamma }}_j)\propto \frac{1}{|{\mathbf{\Omega }}_j{|}^{1-1/\left(2k_j\right)}}$(16) ), (18(18) $m\left(\mathit{y}\right)=\int f(\mathit{y}|\mathcal{V},\lambda )\pi \left(\lambda \right)\prod _{s=1}^r\pi \left({\mathit{V}}_s\right)\mathrm{d}\lambda \prod _{t=1}^r\mathrm{d}{\mathit{V}}_t<\mathrm{\infty },$(18) ), (27(27) $\pi \left({\mathbf{\Omega }}_{r+1}\right)\propto \frac{1}{|{\mathbf{\Omega }}_{r+1}{|}^{1-1/\left(2k_{r+1}\right)}}\exp \{-\frac{1}{2}\mathrm{t}\mathrm{r}\left({\mathbf{\Omega }}_{r+1}^{-1}\right)\}.$(27) ) and (25(25) $L(\mathcal{V},\lambda ;\mathit{y})\le \frac{C_1}{{\left|{\mathit{M}}_1\right|}^{1/2}}$(25) ), we have $\begin{array}{rl}m\left(\mathit{y}\right)& \le \int \frac{\begin{array}{c}{C}_1\exp \{-\mathrm{t}\mathrm{r}\left({\mathbf{\Omega }}_{r+1}^{-1}\right)/2\}\prod _{j=1}^{r+1}{1}_{\{{\omega }_{j1}>\cdots >{\omega }_{j{k}_j}>0\}}\end{array}}{{\left|{\mathit{M}}_1\right|}^{1/2}\prod _{j=1}^{r+1}{\left|{\mathbf{\Omega }}_j\right|}^{1-\frac{1}{2k_j}}}\\ & \times \left[\prod _{j=1}^r\mathrm{d}{\mathbf{\Omega }}_j\mathrm{d}{\mathbf{\Gamma }}_j\right]\mathrm{d}{\mathbf{\Omega }}_{r+1}\\ & <\int \frac{C_1}{{\left|{\mathit{M}}_1\right|}^{1/2}\prod _{j=1}^{r+1}{\left|{\mathbf{\Omega }}_j\right|}^{1-\frac{1}{2k_j}}}\prod _{j=1}^{r+1}\mathrm{d}{\mathbf{\Omega }}_j\triangleq I_0.\end{array}$The definition of $c_{jl,s}$ in (19(19) $c_{jl,s}=1_{\left\{m_j\right(l-1)<s\le m_{j}l\}},$(19) ) yields $\left|{\mathit{M}}_1\right|=\prod _{s=1}^{m_0{k}_0}(1+\sum _{j=1}^{r+1}\sum _{l=1}^{k_j}{c}_{jl,s}{\omega }_{jl}).$Therefore, $I_0\simeq \int \frac{\prod _{j=1}^{r+1}\mathrm{d}{\mathbf{\Omega }}_j}{\begin{array}{c}\prod _{j=1}^{r+1}{\left|{\mathbf{\Omega }}_j\right|}^{1-\frac{1}{2k_j}}\prod _{s=1}^{m_0{k}_0}{(1+\sum _{j=1}^{r+1}\sum _{l=1}^{k_j}{c}_{jl,s}{\omega }_{jl})}^{\frac{1}{2}}\end{array}},$which is finite iff (28) $\sum _{(j,l)\in D}(1-\frac{1}{2k_j})+\frac{1}{2}\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right)$(28) for any $D\in \mathcal{S}\left(\right[r+1\left]\right)$ by employing Lemma 3.2. It's obvious that the inequality (28(28) $\sum _{(j,l)\in D}(1-\frac{1}{2k_j})+\frac{1}{2}\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right)$(28) ) is equivalent to that in (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ).

Applying Lemma 3.1 yields the upper bound on the integrated likelihood function of hyperparameters as $C_0|{\mathit{M}}_1{|}^{-\frac{1}{2}}$, where ${\mathit{M}}_1$ is an $\left(m_0{k}_0\right)\times \left(m_0{k}_0\right)$ matrix and defined in (26(26) $\begin{array}{rl}{\mathit{M}}_1& ={\mathit{I}}_{m_0{k}_0}+\sum _{j=1}^{r+1}{\mathit{D}}_j\mathrm{a}\mathrm{n}\mathrm{d}\\ {\mathit{D}}_j& ={\left(\begin{array}{cc}{\mathbf{\Omega }}_j\otimes {\mathit{I}}_{m_j}& \\ & {\mathit{O}}_{q_j}\end{array}\right)}_{\left(m_0{k}_0\right)\times \left(m_0{k}_0\right)},\end{array}$(26) ). Therefore, the special notation $c_{jl,s}$ can be understood as the indicator of whether eigenvalue ${\omega }_{jl}$ appears in the s-th diagonal element of ${\mathit{M}}_1$ for $j\in [r+1]$, $l\in \left[k_j\right]$ and $s\in \left[m_0{k}_0\right]$. At the same time, the left-hand side of inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) actually stands for the cardinality of set $\left\{s|\mathrm{\exists }\right(j,l)\in D,\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}{c}_{jl,s}>0\}$ for any $D\in \mathcal{S}\left(\right[r+1\left]\right)$.

The cardinality of $\mathcal{S}\left(\right[r+1\left]\right)$ is $2^{\sum _{j\in [r+1]}{k}_j}-1$, which means that the total number of the inequalities to be checked in (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) is exponential with r and the dimensions $k_j$. It has to be admitted that this will impose considerably heavy computational burden in common practice by applying Theorem 3.1 directly. Nevertheless, the researchers have no need to be anxious about the heavy computational burden, because most of the inequalities in Theorem 3.1 are trivial. To conclude this point, we have the following corollary.

Corollary 3.1

Recursively define that $\begin{array}{rl}{R}_1& =\{j|m_j\le r+1,j\in [r+1\left]\right\};\\ R_2& =\{j|m_j\le \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(R_1),j\in R_1\};\\ & ⋮\\ R_p& =\{j|m_j\le \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(R_{p-1}),j\in R_{p-1}\},\end{array}$where p is the smallest positive integer i such that $\{j|m_j>\text{card}(R_i),j\in R_i\}=\mathrm{\varnothing }$. We call the levels within $R_p$ as kernel levels and denote $R_p$ by $R_{\mathrm{k}\mathrm{e}\mathrm{r}}$. The inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) holds for any $D\in \mathcal{S}\left(\right[r+1\left]\right)$ iff the inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) holds for any $D\in \mathcal{S}\left(R_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$. Consequently, if $R_{\mathrm{k}\mathrm{e}\mathrm{r}}=\mathrm{\varnothing }$, then the posterior is always proper.

Proof.

Let $R_1^c=[r+1]-R_1$, thus, $m_j>r+1$ for $j\in R_1^c$ by the definition of $R_1$. For any $D\in \mathcal{S}\left(\right[r+1\left]\right)$, if there exists $j^{\ast }\in R_1^c$ such that $(j^{\ast },l)\in D$ for any $l\in \left[k_j^{\ast }\right]$, then $\begin{array}{rl}\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}& \ge \underset{(j,l)\in D}{max}{m}_j>r+1\\ & =\underset{D\in \mathcal{S}\left(\right[r+1\left]\right)}{max}\sum _{(j,l)\in D}\frac{1}{k_j},\end{array}$which is a trivial one. As a result, inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) holds for any $D\in \mathcal{S}\left(\right[r+1\left]\right)$ iff inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) holds for any $D\in \mathcal{S}\left(R_1\right)$. Since $\underset{D\in \mathcal{S}\left(R_i\right)}{max}\sum _{(j,l)\in D}\frac{1}{k_j}=\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(R_i\right),$it can be recursively shown that inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) holds for any $D\in \mathcal{S}\left(\right[r+1\left]\right)$ iff inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) holds for any $D\in \mathcal{S}\left(R_i\right)$ and $i\in \left[p\right]$, where p is the smallest positive integer i such that $R_i-R_{i+1}=\mathrm{\varnothing }$.

By using the technique of extracting kernel levels, we dramatically narrow down the checking region for posterior propriety. Since we should only check the inequalities for the levels within $R_{\mathrm{k}\mathrm{e}\mathrm{r}}$, substantially reducing the number of inequalities to be checked from $2^{\sum _{j\in [r+1]}{k}_j}-1$ to $2^{\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}{k}_j}-1$. Moreover, Corollary 3.1 also indicates two interesting conclusions depicted as follows.

1. First, it reveals the mechanism how three roles, number of levels, numbers of units in levels and dimensions of levels, affect the posterior propriety simultaneously. Roughly speaking, in the context of GHNL, the more levels with fewer units having lower dimensions, the less likely the posterior is to be proper. For example, if $m_{r-2}=m_{r-1}=m_r=2$ and $k_{r-2}=k_{r-1}=k_r=1$, for $D=\left\{\right(j,l)|j=r-2,r-1,r,l=1\}$, we have (29) $2=\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}<\sum _{(j,l)\in D}\frac{1}{k_j}=3.$(29) Therefore, the posterior propriety can hardly be guaranteed by Theorem 3.1. Conversely, if the units in each level are adequate enough such that the set of kernel levels is empty, i.e., $R_{\mathrm{k}\mathrm{e}\mathrm{r}}=\mathrm{\varnothing }$, then the posterior is always proper. As a consequence, more attention should be focused on the levels with small number of units, namely, the kernel levels.

2. Second, the recommended prior for use at any level in hierarchical modelling is further justified from the aspect of posterior propriety and ease of implementation. For instance, if we switch the prior on ${\mathit{V}}_j$, $j\in \left[r\right]$ from (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) to (30) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-a}\prod _{1\le s<t\le k_j}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(30) where $0<a\le 1$ (a has to be larger than zero by Lemma (3.2)). Then the condition in Theorem 3.1 becomes (31) $\begin{array}{rl}& \sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right),\\ & D\in \mathcal{S}\left(\right[r+1\left]\right).\end{array}$(31) On one hand, when a is greater than but close to zero (denoted by $a⪆0$), all the inequalities in (31(31) $\begin{array}{rl}& \sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right),\\ & D\in \mathcal{S}\left(\right[r+1\left]\right).\end{array}$(31) ) always hold. Thus, the posterior will be always proper. However, it is impractical to decide how small is for a and find such a fixed value that fits all levels. On the other hand, when $a\lesssim 1$, which denotes that a is less than and close to 1, then $2a\ast \text{card}\left(D\right)>\sum _{(j,l)\in D}\frac{1}{k_j}$, concluding that inequality (31(31) $\begin{array}{rl}& \sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right),\\ & D\in \mathcal{S}\left(\right[r+1\left]\right).\end{array}$(31) ) is harder to get reached than inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ), especially for large dimensions $k_j$'s. Therefore, the posterior using prior (30(30) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-a}\prod _{1\le s<t\le k_j}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(30) ) is rather unlikely to be proper than using prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ). Similar to Corollary 3.1, we can recursively define that $\begin{array}{rl}{R}_1^{\ast }& =\{j|m_j\le 2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(E_0\right),j\in [r+1\left]\right\};\\ R_2^{\ast }& =\{j|m_j\le 2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(E_1^{\ast }),j\in R_1^{\ast }\};\\ & ⋮\\ R_{p^{\ast }}^{\ast }& =\{j|m_j\le 2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(E_{p^{\ast }-1}^{\ast }),j\in R_{p^{\ast }-1}^{\ast }\},\end{array}$where $E_0=H\left(\right[r+1\left]\right)$, $E_i^{\ast }=H\left(R_i^{\ast }\right)$ for $i\in \left[p^{\ast }\right]$ and $p^{\ast }$ is the smallest positive integer l such that $\{j|m_j>2a\ast \text{card}(E_l^{\ast }),j\in R_l^{\ast }\}=\mathrm{\varnothing }$. Then the posterior using prior (30(30) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-a}\prod _{1\le s<t\le k_j}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(30) ) instead is proper if (31(31) $\begin{array}{rl}& \sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>2a\times \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right),\\ & D\in \mathcal{S}\left(\right[r+1\left]\right).\end{array}$(31) ) holds for any $D\in \mathcal{S}\left(R_{p^{\ast }}^{\ast }\right)$. When $a\lesssim 1$, $\text{card}\left(R_{p^{\ast }}^{\ast }\right)$ will be remarkably larger than $\text{card}\left(R_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$ for large values of $k_j$'s, imposing a dramatically heavier burden of checking inequalities than that for prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ). Above all, one sensible choice for a is that let a be inversely proportional to $k_j$ for level j, which takes both practical and theoretical considerations into account.

The upper bound on $\sum _{(j,l)\in D}\frac{1}{k_j}$ as $\text{card}\left(R_i\right)$ for $D\in \mathcal{S}\left(R_i\right)$, leads to an effective way to extract kernel levels as presented in Corollary 3.1, but this upper bound is still too rough. Next, an elaborate upper bound on $\sum _{(j,l)\in D}\frac{1}{k_j}$ is demonstrated and a sufficient condition of clean form for posterior propriety is derived then.

Theorem 3.2

Consider the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with priors (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) and (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on $\mathit{\eta }$ and ${\mathit{V}}_i\in \mathcal{V}$, respectively. Denote that $m^{\ast }=\underset{j\in [r+1]}{min}{m}_j=\underset{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}{min}{m}_j$. When $d\ge 2$, the posterior is always proper if (32) $\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}<m^{\ast }.$(32)

Proof.

For $D\in \mathcal{S}\left(\right[r+1\left]\right)$, define that (33) $\begin{array}{rl}L\left(D\right)& =\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}},\mathrm{a}\mathrm{n}\mathrm{d}\\ R\left(D\right)& =\left\{j|\right(j,l)\in D\}.\end{array}$(33) It follows from Corollary 3.1 that we only need to prove $\sum _{(j,l)\in D}\frac{1}{k_j}<L\left(D\right),\mathrm{\forall }D\in \mathcal{S}\left(R_{\mathrm{k}\mathrm{e}\mathrm{r}}\right).$Also, for any D belonging to $\mathcal{S}\left(R_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$, we have (34) $\begin{array}{rl}\sum _{(j,l)\in D}\frac{1}{k_j}& =\sum _{j\in R\left(D\right)}\frac{1}{k_j}\sum _{l=1}^{k_j}{1}_{\left\{\right(j,l)\in D\}}\\ & \le \left(\sum _{j\in R\left(D\right)}\frac{1}{k_j}\right)\underset{j\in R\left(D\right)}{max}\sum _{l=1}^{k_j}{1}_{\left\{\right(j,l)\in D\}}.\end{array}$(34) Distinct eigenvalues from the same level ($\le r$-th) never occur in the same row of matrix ${\mathit{M}}_1$. Mathematically, $c_{j{l}_1,s}{c}_{j{l}_2,s}=0$, for $1\le l_1<l_2\le k_j$, $j\in \left[r\right]$, $s\in \left[m_0{k}_0\right]$. Thus, for any $j\in \left[r\right]$, (35) $\sum _{l=1}^{k_j}{1}_{\left\{\right(j,l)\in D\}}\le \frac{1}{m_j}L\left(D\right),$(35) which is also true for j = r + 1 since $k_{r+1}=1$. It is obvious that $\underset{j\in [r+1]}{min}{m}_j=\underset{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}{min}{m}_j$ by the definition of $R_{\mathrm{k}\mathrm{e}\mathrm{r}}$, and that is denoted by $m^{\ast }$. Combining (35(35) $\sum _{l=1}^{k_j}{1}_{\left\{\right(j,l)\in D\}}\le \frac{1}{m_j}L\left(D\right),$(35) ) with (34(34) $\begin{array}{rl}\sum _{(j,l)\in D}\frac{1}{k_j}& =\sum _{j\in R\left(D\right)}\frac{1}{k_j}\sum _{l=1}^{k_j}{1}_{\left\{\right(j,l)\in D\}}\\ & \le \left(\sum _{j\in R\left(D\right)}\frac{1}{k_j}\right)\underset{j\in R\left(D\right)}{max}\sum _{l=1}^{k_j}{1}_{\left\{\right(j,l)\in D\}}.\end{array}$(34) ) yields (36) $\begin{array}{rl}\sum _{(j,l)\in D}\frac{1}{k_j}& \le \left(\sum _{j\in R\left(D\right)}\frac{1}{k_j}\right)\frac{L\left(D\right)}{\underset{j\in R\left(D\right)}{min}{m}_j}\\ & \le \frac{1}{m^{\ast }}\left(\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}\right)L\left(D\right)<L\left(D\right).\end{array}$(36)

According to the proof procedure of Theorem 3.2, it can be deduced that $\sum _{j\in [r+1]}\frac{1}{k_j}<m^{\ast }$ is also sufficient for posterior propriety. Obviously, the condition in Theorem 3.2 is easier to be satisfied. Theorem 3.2 reveals that for fixed $m^{\ast }$, the posterior is more likely to be proper for higher dimensions of the units in the kernel levels. Theorem 3.2 also provides the researchers with a powerful tool to check the posterior propriety quickly.

Remark 3.1

Consider the model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with r = 1, namely, a two-level hierarchical model. When $d\ge 2$, we have $m^{\ast }=min\{d,m_1\}\ge 2$. Then the posterior using the recommended prior is always proper for $k_1\ge 2$ referring to Theorem 3.2.

Example 3.2

Continue with Example 2.1

Consider the GHNL modelling of the mixed-effect ANOVA as (10(10) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& (\mathit{y}|{\text{β}}^{\ast },{\mathbf{\Sigma }}_0)& \sim & \begin{array}{l}{N}_{m_{0}p}({\mathit{Z}}_0{\text{β}}^{\ast },{\mathit{I}}_{m_0}\otimes {\mathbf{\Sigma }}_0);\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\text{β}}^{\ast }|{\mathit{\alpha }}^{\ast },{\mathit{V}}_{\text{β}})& \sim & \begin{array}{l}{N}_{m_{1}p}({\mathit{Z}}_1{\mathit{\alpha }}^{\ast },{\mathit{I}}_{m_1}\otimes {\mathit{V}}_{\text{β}});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:& ({\mathit{\alpha }}^{\ast }|\mathit{\eta },{\mathit{V}}_{\alpha })& \sim & \begin{array}{l}{N}_{m_{2}p}({\mathit{Z}}_2\mathit{\eta },{\mathit{I}}_{m_2}\otimes {\mathit{V}}_{\alpha }),\end{array}\end{array}$(10) ), which is a 3-level hierarchical model with r = 2, $m_0=s_1{s}_2{s}_3$, $m_1=s_1{s}_2$, $m_2=s_1$, $m_3=p$, $k_0=k_1=k_2=p$ and $k_3=1$. It is natural to assume that we have at least two schools, each school has at least two classes and each class has at least two students, namely, $s_1\ge 2$, $s_2\ge 2$ and $s_3\ge 2$. If (a) p>2 and $s_1\ge 2$ or (b) $p\ge 2$ and $s_1>2$ holds, it can be readily derived that the set of kernel levels is empty. Thus, the posterior is always proper according to Corollary 3.1. When $s_1=p=2$, the set of kernel levels is $R_2=\{2,3\}$, since $\frac{1}{k_2}+\frac{1}{k_3}<2$, then the posterior is proper by applying Theorem 3.2. In conclusion, the posterior using the recommended prior is always proper when $p\ge 2$.

In Berger et al. (2020b)'s work, for a technical reason, they assumed k = sp for 3-level hierarchical normal model (2(2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) ) such that the design matrices for units within level 2 and 3 are square matrices. They eventually reached a conclusion that the posterior employing the recommended prior in model (2(2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) ) is always proper for k = sp and $p\ge 2$ (p is the dimension of hypermean in model (2(2) $\{\begin{array}{ll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:{\mathit{y}}_i={\mathit{\theta }}_i+{\mathcal{N}}_k(\mathbf{0},{\mathit{I}}_k),& i\in \left[m\right];\\ \begin{array}{l}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:{\mathit{\theta }}_i={\mathit{Z}}_i\text{β}+{\mathcal{N}}_k(\mathbf{0},\mathit{V}),\end{array}& \begin{array}{l}{\text{β}}^{\mathrm{\top }}=({\text{β}}_1^{\mathrm{\top }},\\ \dots ,{\text{β}}_s^{\mathrm{\top }});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:{\text{β}}_j=\mathit{\eta }+{\mathcal{N}}_p(\mathbf{0},\mathit{W}),& j\in \left[s\right],\end{array}$(2) )). However, the assumption that the design matrices for the units at high levels are square matrices appears to be unnatural and hard to be interpreted in practice. Nevertheless, we still generalize Berger et al. (2020b)'s result to the GHNL model in the following so as to draw a consistent conclusion with theirs.

Corollary 3.2

Consider the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with priors (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) and (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on $\mathit{\eta }$ and ${\mathit{V}}_i\in \mathcal{V}$, respectively. Assume that $d\ge 2$ and $k_j=m_{j+1}{k}_{j+1}$, $j\in \left[r\right]$. Then the posterior is always proper.

Proof.

Since $m_j\ge 2$, $j\in [r+1]$, then $m^{\ast }\ge 2$. It remains to show that $\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}<2$ by Theorem 3.2. By utilizing the condition that $k_j=m_{j+1}{k}_{j+1}$ for $j\in \left[r\right]$ and $k_{r+1}=1$, we have $k_j=\prod _{s=j+1}^{r+1}{m}_s\ge 2^{r+1-j}$ for $j\in \left[r\right]$. Thus, $\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}\le \sum _{j\in [r+1]}\frac{1}{k_j}\le \sum _{j\in [r+1]}\frac{1}{2^{r+1-j}}=2-\frac{1}{2^r}<2,$which completes the proof.

Above all, a general procedure for checking the posterior propriety of the GHNL models (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) employing the recommended prior for $d\ge 2$ can be summarized as follows.

Guidance for checking the posterior propriety when $d\ge 2$:

1. If the design matrices for each unit in each level are square matrices, then the posterior is proper, otherwise, turn to (b).

2. Derive the set of kernel levels, $R_{\mathrm{k}\mathrm{e}\mathrm{r}}$. If $R_{\mathrm{k}\mathrm{e}\mathrm{r}}=\mathrm{\varnothing }$ or inequality (32(32) $\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}<m^{\ast }.$(32) ) holds, the the posterior is proper. If neither, turn to (c).

3. Check inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ) for all D belonging to $\mathcal{S}\left(R_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$. If that always holds, the the posterior is proper. If not, the posterior propriety can hardly be guaranteed.

3.3. Conditions for the posterior to be proper when d = 1

It is quite common that the dimension of the fixed effect $\mathit{\eta }$ is only one in practice. However, when $d=1$, note that $\begin{array}{rl}1& =m_{r+1}=\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}\\ & =\sum _{(j,l)\in D}\frac{1}{k_j}=\frac{1}{k_{r+1}}=1\end{array}$for $D=\left\{\right(r+1,1\left)\right\}$, resulting in the failure of the sufficient condition in Theorem 3.1. Therefore, in this subsection, we mainly reinvestigate the conditions for the the posterior to be proper for d = 1.

Theorem 3.3

Consider the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with constant prior on η and prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on ${\mathit{V}}_j\in \mathcal{V}$. When d = 1, the posterior is proper if (37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) holds for all $D\in \mathcal{S}\left(\right[r\left]\right)$.

Proof.

When d = 1, vector $\mathit{\eta }$ will degenerate into a scalar η. Hence the prior (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ) on $\mathit{\eta }$ becomes a constant prior on η. Based on (5(5) $(\mathit{y}|\mathit{\eta },\mathcal{V})\sim N_{m_0{k}_0}({\mathit{X}}_r\mathit{\eta },\mathbf{\Delta }),$(5) ), by integrating out over η, we can get the upper bound on the integrated likelihood of $\mathcal{V}$ after dropping the exponential term (less than one) $L\left(\mathcal{V}\right)<\frac{1}{{\left|\mathbf{\Delta }\right|}^{\frac{1}{2}}{\left|{\mathit{X}}_r^{\mathrm{\top }}{\mathbf{\Delta }}^{-1}{\mathit{X}}_r\right|}^{\frac{1}{2}}}.$Since $\left|{\mathit{X}}_r^{\mathrm{\top }}{\mathbf{\Delta }}^{-1}{\mathit{X}}_r\right|\ge {\mathit{X}}_r^{\mathrm{\top }}{\mathit{X}}_r{\lambda }_{min}\left({\mathbf{\Delta }}^{-1}\right)={\mathit{X}}_r^{\mathrm{\top }}{\mathit{X}}_r{\lambda }_{max}{\left(\mathbf{\Delta }\right)}^{-1}$, using Lemma 3.1 (a) and (b), we have $\begin{array}{rl}L\left(\mathcal{V}\right)& <C_0{\left|\mathbf{\Delta }\right|}^{-\frac{1}{2}}{(1+\sum _{j=1}^r{\omega }_{j1})}^{\frac{1}{2}}\\ & \le C_1{(1+\sum _{j=1}^r{\omega }_{j1})}^{\frac{1}{2}}{\left|{\mathit{M}}_2\right|}^{-\frac{1}{2}},\end{array}$where $C_0$ and $C_1$ are constants that are independent of the ${\mathbf{\Omega }}_j$ for $j\in \left[r\right]$, ${\mathit{M}}_2={\mathit{I}}_{m_0{k}_0}+\sum _{j=1}^r{\mathit{D}}_j$ and ${\mathit{D}}_j$'s are defined in (26(26) $\begin{array}{rl}{\mathit{M}}_1& ={\mathit{I}}_{m_0{k}_0}+\sum _{j=1}^{r+1}{\mathit{D}}_j\mathrm{a}\mathrm{n}\mathrm{d}\\ {\mathit{D}}_j& ={\left(\begin{array}{cc}{\mathbf{\Omega }}_j\otimes {\mathit{I}}_{m_j}& \\ & {\mathit{O}}_{q_j}\end{array}\right)}_{\left(m_0{k}_0\right)\times \left(m_0{k}_0\right)},\end{array}$(26) ). Similar to the proof of Theorem 3.1, we can derive the upper bound on $m\left(\mathit{y}\right)$ as $m\left(\mathit{y}\right)\le \int \frac{C_1{(1+\sum _{j=1}^r{\omega }_{j1})}^{\frac{1}{2}}}{{\left|{\mathit{M}}_2\right|}^{1/2}\prod _{j=1}^r{\left|{\mathbf{\Omega }}_j\right|}^{1-\frac{1}{2k_j}}}\prod _{j=1}^r\mathrm{d}{\mathbf{\Omega }}_j\triangleq I_1.$It follows from the definition of $c_{jl,s}$ in (19(19) $c_{jl,s}=1_{\left\{m_j\right(l-1)<s\le m_{j}l\}},$(19) ) that $\left|{\mathit{M}}_2\right|=\prod _{s=1}^{m_0{k}_0}(1+\sum _{j=1}^r\sum _{l=1}^{k_j}{c}_{jl,s}{\omega }_{jl}).$Thus, $I_0\simeq \int \frac{{(1+\sum _{j=1}^r{\omega }_{j1})}^{\frac{1}{2}}\prod _{j=1}^r\mathrm{d}{\mathbf{\Omega }}_j}{\begin{array}{c}\prod _{j=1}^r{\left|{\mathbf{\Omega }}_j\right|}^{1-\frac{1}{2k_j}}\prod _{s=1}^{m_0{k}_0}{(1+\sum _{j=1}^r\sum _{l=1}^{k_j}{c}_{jl,s}{\omega }_{jl})}^{\frac{1}{2}}\end{array}},$which is finite iff (38) $\begin{array}{rl}& \sum _{(j,l)\in D}(1-\frac{1}{2k_j})+\frac{1}{2}\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}\\ & -\frac{1}{2}{1}_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right)\end{array}$(38) for any $D\in \mathcal{S}\left(\right[r\left]\right)$ by employing Lemma 3.2. It's obvious that inequality (38(38) $\begin{array}{rl}& \sum _{(j,l)\in D}(1-\frac{1}{2k_j})+\frac{1}{2}\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}\\ & -\frac{1}{2}{1}_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right)\end{array}$(38) ) is equivalent to that in (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ).

Resembling the interpretation of Theorem 3.1, the left-hand side of inequality (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ) actually denotes the cardinality of set $\left\{s|\mathrm{\exists }\right(j,l)\in D,\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}{c}_{jl,s}>0\}$ for any $D\in \mathcal{S}\left(\right[r\left]\right)$. To reduce the burden of checking inequalities, we have the following corollary.

Corollary 3.3

Recursively define that $\begin{array}{rl}{\tilde{R}}_1& =\{j|m_j\le r+1,j\in [r\left]\right\};\\ {\tilde{R}}_2& =\{j|m_j\le \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}({\tilde{R}}_1)+1,j\in {\tilde{R}}_1\};\\ & ⋮\\ {\tilde{R}}_q& =\{j|m_j\le \mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}({\tilde{R}}_{p-1})+1,j\in {\tilde{R}}_{p-1}\},\end{array}$where q is the smallest positive integer i such that $\{j|m_j>\text{card}({\tilde{R}}_i)+1,j\in {\tilde{R}}_i\}=\mathrm{\varnothing }$. We call the levels within ${\tilde{R}}_q$ as kernel levels and denote ${\tilde{R}}_q$ by ${\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}$. Inequality (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ) holds for any $D\in \mathcal{S}\left(\right[r\left]\right)$ iff inequality (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ) holds for any $D\in \mathcal{S}\left({\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$. Consequently, if ${\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}=\mathrm{\varnothing }$, then the resulting posterior is always proper.

In the process of extracting kernel levels, the thresholds of $m_j$ to split up levels are increased by one in Corollary 3.3, when compared with that in Corollary 3.1, and this is because the upper bound on the right-hand side of inequality (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ) is increased by one than that of inequality (23(23) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{k_j},$(23) ). Except for this point, the proof of Corollary 3.1 is same as that of Corollary 3.3. For good measure, a simple tool to check the posterior propriety is shown as follows, which is a counterpart of Theorem 3.2 for d = 1.

Theorem 3.4

Consider the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with constant prior on η and prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on ${\mathit{V}}_j\in \mathcal{V}$. When $d=1$, the posterior is always proper if (39) $\sum _{j\in {\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}<m^{\ast }-1,$(39) where $m^{\ast }=\underset{j\in \left[r\right]}{min}{m}_j$ and ${\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}$ is the derived set of kernel levels.

Proof.

For any $D\in \mathcal{S}\left(\right[r\left]\right)$, define $L\left(D\right)$ and $R\left(D\right)$ in the same way as that in (33(33) $\begin{array}{rl}L\left(D\right)& =\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}},\mathrm{a}\mathrm{n}\mathrm{d}\\ R\left(D\right)& =\left\{j|\right(j,l)\in D\}.\end{array}$(33) ). According to Corollary (3.3), it suffices to show that $L\left(D\right)>\sum _{(j,l)\in D}\frac{1}{k_j}+1$for any $D\in \mathcal{S}\left({\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$. Similar to (36(36) $\begin{array}{rl}\sum _{(j,l)\in D}\frac{1}{k_j}& \le \left(\sum _{j\in R\left(D\right)}\frac{1}{k_j}\right)\frac{L\left(D\right)}{\underset{j\in R\left(D\right)}{min}{m}_j}\\ & \le \frac{1}{m^{\ast }}\left(\sum _{j\in R_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}\right)L\left(D\right)<L\left(D\right).\end{array}$(36) ), we have $\begin{array}{rl}\sum _{(j,l)\in D}\frac{1}{k_j}& \le \frac{1}{m^{\ast }}\left(\sum _{j\in R\left(D\right)}\frac{1}{k_j}\right)L\left(D\right)\\ & \le \frac{1}{m^{\ast }}\left(\sum _{j\in {\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}\right)L\left(D\right)<L\left(D\right)-\frac{L\left(D\right)}{m^{\ast }},\end{array}$for any $D\in \mathcal{S}\left({\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$. Since $L\left(D\right)\ge m^{\ast }$ always holds, then the proof is completed.

Remark 3.2

In model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ), when r = 1 and d = 1, suppose $m_1\ge 2$ and $k_1\ge 2$. The posterior using the recommended prior is always proper by applying Theorem 3.4.

Remark 3.3

If all $k_j$'s are equal to one, the sufficient condition in Theorem 3.3 can be simplified as $\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(D\right)<\underset{j\in D}{max}{m}_j-1,$where $D\in \mathcal{F}\left(\right[r\left]\right)$. By employing the technique of extracting kernel levels, Theorem 3.4 is equivalent to Theorem 3.3, rather than a mere sufficient condition.

Example 3.3

Continue with Example 2.1

Consider model (10(10) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& (\mathit{y}|{\text{β}}^{\ast },{\mathbf{\Sigma }}_0)& \sim & \begin{array}{l}{N}_{m_{0}p}({\mathit{Z}}_0{\text{β}}^{\ast },{\mathit{I}}_{m_0}\otimes {\mathbf{\Sigma }}_0);\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\text{β}}^{\ast }|{\mathit{\alpha }}^{\ast },{\mathit{V}}_{\text{β}})& \sim & \begin{array}{l}{N}_{m_{1}p}({\mathit{Z}}_1{\mathit{\alpha }}^{\ast },{\mathit{I}}_{m_1}\otimes {\mathit{V}}_{\text{β}});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}3:& ({\mathit{\alpha }}^{\ast }|\mathit{\eta },{\mathit{V}}_{\alpha })& \sim & \begin{array}{l}{N}_{m_{2}p}({\mathit{Z}}_2\mathit{\eta },{\mathit{I}}_{m_2}\otimes {\mathit{V}}_{\alpha }),\end{array}\end{array}$(10) ) with assuming that $s_1\ge 2$, $s_2\ge 2$ and $s_3\ge 2$. If p = 1, we have $k_0=k_1=k_2=k_3=m_3=1$. When $s_1>2$, we can easily derive the set of kernel levels as empty set. Thus, the posterior is always proper by Corollary 3.3. If $s_1=2$, inequality (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ) fails, and the posterior propriety can hardly be guaranteed. Consequently, when p = 1, the posterior using the recommended prior is always proper for $s_1>2$.

Next, we generalize Berger et al. (2020b)'s result to the GHNL models for d = 1, assuming that the design matrices ${\mathit{Z}}_{j{i}_j}$'s for units are square matrices.

Corollary 3.4

When $d=1$, consider the same model and prior as Theorem 3.3. Suppose $m_j\ge 3$ and $k_j=m_{j+1}{k}_{j+1}$, $j\in \left[r\right]$. Then the posterior is always proper.

Proof.

It follows from Theorem 3.4 that we should only present $\sum _{j\in \left[r\right]}\frac{1}{k_j}<m^{\ast }-1.$According to the conditions that $k_j=m_{j+1}{k}_{j+1}$, $j\in \left[r\right]$ and $m_{r+1}=k_{r+1}=1$, we have $k_j=\prod _{s=j+1}^r{m}_s$, $j\in \left[r\right]$. Thus, $\sum _{j\in \left[r\right]}\frac{1}{k_j}\le \sum _{j\in \left[r\right]}\frac{1}{3^{r-j}}=\frac{3}{2}(1-\frac{1}{3^r})<\frac{3}{2}<2\le m^{\ast }-1.$

Summing up the theoretical results above, a general procedure for checking the posterior propriety of the GHNL models (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) employing the recommended prior for d = 1 can be depicted as follows.

Guidance for checking the posterior propriety when d = 1:

1. If the design matrices for each unit in each level are square matrices and $m_j\ge 3$, $j\in \left[r\right]$, then the posterior is proper, otherwise, turn to (b).

2. Derive the set of kernel levels ${\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}$. If ${\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}=\mathrm{\varnothing }$ or inequality (39(39) $\sum _{j\in {\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}}\frac{1}{k_j}<m^{\ast }-1,$(39) ) holds, then the posterior is proper. If neither, turn to (c).

3. Check inequality (37(37) $\sum _{s=1}^{m_0{k}_0}{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{k_j}$(37) ) for all D belonging to $\mathcal{S}\left({\tilde{R}}_{\mathrm{k}\mathrm{e}\mathrm{r}}\right)$. If that always holds, then the posterior is proper. If not, the posterior propriety can hardly be guaranteed.

4. Computation

In this section, we consider the MCMC sampling from the posterior arising from the model in Section 2. For the GHNL model (4(4) $\{\begin{array}{llll}\mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}1:& \left(\mathit{y}|{\mathit{\theta }}_1\right)& \sim & \begin{array}{l}{N}_{m_0{k}_0}({\mathit{Z}}_0{\mathit{\theta }}_1,{\mathit{I}}_{m_0{k}_0});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}2:& ({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1)& \sim & \begin{array}{l}{N}_{m_1{k}_1}({\mathit{Z}}_1{\mathit{\theta }}_2,\\ {\mathit{I}}_{m_1}\otimes {\mathit{V}}_1);\end{array}\\ ⋮& ⋮& ⋮\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r:& ({\mathit{\theta }}_{r-1}|{\mathit{\theta }}_r,{\mathit{V}}_{r-1})& \sim & \begin{array}{l}{N}_{m_{r-1}{k}_{r-1}}({\mathit{Z}}_{r-1}{\mathit{\theta }}_r,{\mathit{I}}_{m_{r-1}}\\ \otimes {\mathit{V}}_{r-1});\end{array}\\ \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}r+1:& ({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)& \sim & \begin{array}{l}{N}_{m_r{k}_r}({\mathit{Z}}_r\mathit{\eta },\\ {\mathit{I}}_{m_r}\otimes {\mathit{V}}_r).\end{array}\end{array}$(4) ) with prior (14(14) $\left(\mathit{\eta }|\lambda \right)\sim N_d(\mathbf{0},\lambda {\mathit{I}}_d)\mathrm{a}\mathrm{n}\mathrm{d}\left[\lambda \right]\propto {\lambda }^{-1/2}\exp (-\frac{1}{2\lambda }),$(14) ) and (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on $\mathit{\eta }$ and $\mathcal{V}$, respectively, the joint posterior of $(\mathbf{\Theta },\mathcal{V},\mathit{\eta },\lambda )$ can be written as (40) $\begin{array}{rl}& \pi (\mathbf{\Theta },\mathcal{V},\mathit{\eta },\lambda |\mathit{y})\propto f\left(\mathit{y}|{\mathit{\theta }}_1\right)\prod _{j=1}^{r-1}f({\mathit{\theta }}_j|{\mathit{\theta }}_{j+1},{\mathit{V}}_j)f({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)\\ & \times \prod _{s=1}^r\pi \left({\mathit{V}}_s\right)\pi \left(\mathit{\eta }|\lambda \right)\pi \left(\lambda \right).\end{array}$(40) Sampling $(\mathbf{\Theta },\mathcal{V},\mathit{\eta },\lambda )$ from the posterior density (40(40) $\begin{array}{rl}& \pi (\mathbf{\Theta },\mathcal{V},\mathit{\eta },\lambda |\mathit{y})\propto f\left(\mathit{y}|{\mathit{\theta }}_1\right)\prod _{j=1}^{r-1}f({\mathit{\theta }}_j|{\mathit{\theta }}_{j+1},{\mathit{V}}_j)f({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)\\ & \times \prod _{s=1}^r\pi \left({\mathit{V}}_s\right)\pi \left(\mathit{\eta }|\lambda \right)\pi \left(\lambda \right).\end{array}$(40) ) can be handled by Gibbs sampling method. The main difficulty of the computation is to sample the covariance matrices ${\mathit{V}}_j$'s efficiently.

4.1. Gibbs sampling for input effects

The full conditionals of the input effects $(\mathbf{\Theta },\mathit{\eta })$ can be derived from the joint posterior (40(40) $\begin{array}{rl}& \pi (\mathbf{\Theta },\mathcal{V},\mathit{\eta },\lambda |\mathit{y})\propto f\left(\mathit{y}|{\mathit{\theta }}_1\right)\prod _{j=1}^{r-1}f({\mathit{\theta }}_j|{\mathit{\theta }}_{j+1},{\mathit{V}}_j)f({\mathit{\theta }}_r|\mathit{\eta },{\mathit{V}}_r)\\ & \times \prod _{s=1}^r\pi \left({\mathit{V}}_s\right)\pi \left(\mathit{\eta }|\lambda \right)\pi \left(\lambda \right).\end{array}$(40) ) and are illustrated as follows.

1. Conditioning on ${\mathit{\theta }}_2$ and ${\mathit{V}}_1$, the posterior distribution of ${\mathit{\theta }}_1$ is (41) $({\mathit{\theta }}_1|{\mathit{\theta }}_2,{\mathit{V}}_1;\mathit{y})\sim N_{m_1{k}_1}({\tilde{\mathit{\theta }}}_1,{\tilde{\mathit{V}}}_1),$(41) where $\begin{array}{rl}{\tilde{\mathit{V}}}_1& ={({\mathit{Z}}_0^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}{\mathit{Z}}_0+{\mathit{I}}_{m_1}\otimes {\mathit{V}}_1^{-1})}^{-1},\\ {\tilde{\mathit{\theta }}}_1& ={\tilde{\mathit{V}}}_1[{\mathit{Z}}_0^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathit{y}+({\mathit{I}}_{m_1}\otimes {\mathit{V}}_1^{-1}){\mathit{Z}}_1{\mathit{\theta }}_2].\end{array}$

2. The full conditional posteriors of ${\mathit{\theta }}_j$, $j=2,\dots ,r$ have the forms: (42) $({\mathit{\theta }}_j|{\mathit{\theta }}_{j-1},{\mathit{\theta }}_{j+1},{\mathit{V}}_{j-1},{\mathit{V}}_j)\sim N_{m_j{k}_j}({\tilde{\mathit{\theta }}}_j,{\tilde{\mathit{V}}}_j),$(42) where $\begin{array}{rl}{\tilde{\mathit{V}}}_j& =[{\mathit{Z}}_{j-1}^{\mathrm{\top }}({\mathit{I}}_{m_{j-1}}\otimes {\mathit{V}}_{j-1}^{-1}){\mathit{Z}}_{j-1}\\ & +{{\mathit{I}}_{m_j}\otimes {\mathit{V}}_j^{-1}]}^{-1},\\ {\tilde{\mathit{\theta }}}_j& ={\tilde{\mathit{V}}}_j[{\mathit{Z}}_{j-1}^{\mathrm{\top }}({\mathit{I}}_{m_{j-1}}\otimes {\mathit{V}}_{j-1}^{-1}){\mathit{\theta }}_{j-1}\\ & +({\mathit{I}}_{m_j}\otimes {\mathit{V}}_j^{-1}){\mathit{Z}}_j{\mathit{\theta }}_{j+1}],\end{array}$ and ${\mathit{\theta }}_{r+1}\triangleq \mathit{\eta }$.

3. By using (14(14) $\left(\mathit{\eta }|\lambda \right)\sim N_d(\mathbf{0},\lambda {\mathit{I}}_d)\mathrm{a}\mathrm{n}\mathrm{d}\left[\lambda \right]\propto {\lambda }^{-1/2}\exp (-\frac{1}{2\lambda }),$(14) ), the full conditional of $\mathit{\eta }$ can be derived as (43) $(\mathit{\eta }|{\mathit{\theta }}_r,\lambda ,{\mathit{V}}_r)\sim N_d(\tilde{\mathit{\eta }},{\tilde{\mathit{V}}}_{\mathit{\eta }}),$(43) where $\begin{array}{rl}{\tilde{\mathit{V}}}_{\mathit{\eta }}& ={[{\mathit{Z}}_r^{\mathrm{\top }}({\mathit{I}}_{m_r}\otimes {\mathit{V}}_r^{-1}){\mathit{Z}}_r+{\lambda }^{-1}{\mathit{I}}_d]}^{-1},\\ \tilde{\mathit{\eta }}& ={\tilde{\mathit{V}}}_{\mathit{\eta }}{\mathit{Z}}_r^{\mathrm{\top }}({\mathit{I}}_{m_r}\otimes {\mathit{V}}_r^{-1}){\mathit{\theta }}_r.\end{array}$

Input effects ${\mathit{\theta }}_j\in \mathbf{\Theta }$ and $\mathit{\eta }$ can be readily sampled from their conditionals during the Gibbs sampling procedure, as their full conditional posterior distributions are all standard distributions.

4.2. Gibbs sampling for variance components

The variance components which include ${\mathit{V}}_j$'s and λ can be updated from their full conditionals, and these conditionals have densities as follows.

1. Given $\mathit{\eta }$, the conditional posterior density of λ is (44) $\pi \left(\lambda |\mathit{\eta }\right)\propto {\lambda }^{-\frac{d+1}{2}}\exp \{-\frac{1+\Vert \mathit{\eta }{\Vert }^2}{2\lambda }\},$(44) which is actually an inverse gamma distribution as $IG(\frac{d-1}{2},\frac{1+\Vert \mathit{\eta }{\Vert }^2}{2})$.

2. For $j\in \left[r\right]$, define that $t_j=\frac{m_j}{2}+1-\frac{1}{2k_j}$. The marginal posterior density of ${\mathit{V}}_j$ given $({\mathit{\theta }}_j,{\mathit{\theta }}_{j+1})$ is (45) $\begin{array}{rl}\pi ({\mathit{V}}_j|{\mathit{\theta }}_j,{\mathit{\theta }}_{j+1})& \propto \frac{1}{{\left|{\mathit{V}}_j\right|}^{t_j}\prod _{1\le s<t\le k_j}({\omega }_{js}-{\omega }_{jt})}\\ & \times \mathrm{e}\mathrm{t}\mathrm{r}\{-\frac{1}{2}{\mathit{V}}_j^{-1}{\mathit{H}}_j\},\end{array}$(45) where $\mathrm{e}\mathrm{t}\mathrm{r}\left(\mathit{A}\right)$ denotes $\exp (\mathrm{t}\mathrm{r}(\mathit{A}\left)\right)$ for a square matrix $\mathit{A}$, and $\begin{array}{rl}{\mathit{H}}_j& \triangleq {\mathit{H}}_j({\mathit{\theta }}_j,{\mathit{\theta }}_{j+1})\\ & =\sum _{i=1}^{m_j}({\mathit{\theta }}_{ji}-{\mathit{Z}}_{ji}{\mathit{\theta }}_{j+1})({\mathit{\theta }}_{ji}-{\mathit{Z}}_{ji}{\mathit{\theta }}_{j+1}{)}^{\mathrm{\top }}.\end{array}$

The updating of λ can be simply carried out by sampling from an inverse gamma distribution. The full conditional posteriors of ${\mathit{V}}_j$, (45(45) $\begin{array}{rl}\pi ({\mathit{V}}_j|{\mathit{\theta }}_j,{\mathit{\theta }}_{j+1})& \propto \frac{1}{{\left|{\mathit{V}}_j\right|}^{t_j}\prod _{1\le s<t\le k_j}({\omega }_{js}-{\omega }_{jt})}\\ & \times \mathrm{e}\mathrm{t}\mathrm{r}\{-\frac{1}{2}{\mathit{V}}_j^{-1}{\mathit{H}}_j\},\end{array}$(45) ) are actually distributed as a recently proposed class of prior distributions by Berger et al. (2020a) for the covariance matrix, which is called the Shrinkage Inverse Wishart (SIW) distributions. The new class $\mathrm{S}\mathrm{I}\mathrm{W}(a,\mathit{H})$ for a $k\times k$ covariance matrix $\mathit{W}$ has the density as (46) ${\pi }^{\mathrm{S}\mathrm{I}\mathrm{W}}(\mathit{W}|a,\mathit{H})\propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}{\mathit{W}}^{-1}\mathit{H})}{{\left|\mathit{W}\right|}^a\prod _{i<j}({\nu }_i-{\nu }_j)},$(46) where ${\nu }_1>{\nu }_2>\cdots >{\nu }_k>0$ are the ordered eigenvalues of $\mathit{W}$, a is a real constant and $\mathit{H}$ is a $k\times k$ non-negative definite matrix. Thus, ${\mathit{V}}_j$ are distributed as $\mathrm{S}\mathrm{I}\mathrm{W}(t_j,{\mathit{H}}_j)$, $j\in \left[r\right]$. To sample the covariance matrices from the full conditional posteriors, the previously suggested methods include the Metropolis-Hastings algorithm (cf. Berger et al., 2005) and Hit-and-run method (cf. Yang & Berger, 1994). The two methods both generate full candidate matrices by utilizing full-parameter proposal distributions, resulting in that they only work for moderate dimensions of the covariance matrices. To tackle this issue, Berger et al. (2020a) proposed a powerful Gibbs method for efficiently sampling the covariance matrices from their conditional densities and this new method works for higher dimensions k. The audience can refer to Berger et al. (2020a) or Appendix 3 for details of this Gibbs sampling method. According to the simulation results of Berger et al. (2020a), the new Gibbs method outperforms the Metropolis-Hastings and Hit-and-run methods for moderate dimensions and work for k up to 100, while the other two algorithms break down in much lower dimensions.

In the framework of 2-level HNLM, Berger et al. (2020b) compared the numerical performance, from the mean square error (MSE) perspective, of a dozen of objective hyperpriors, which are the product of three objective hyperpriors for the hypermean and four objective hyperpriors for the hypercovariance matrix. Priors on the hypermean include constant prior, conjugate prior and the recommended prior (12(12) $\begin{array}{rl}\pi \left(\mathit{\eta }\right)& \propto \frac{1}{(1+\Vert \mathit{\eta }{\Vert }^2{)}^{(d-1)/2}},\mathit{\eta }\in {\mathbb{R}}^d,\end{array}$(12) ). Priors on the hypercovariance matrix include constant prior, hierarchical Jefferys prior, hierarchical reference prior and the recommended prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ). Their simulation results have shown that the recommended combination of hyperpriors dominates all the others in terms of Bayes risk, and the constant prior on the hypercovariance performs the worst. However, neither of the two remaining choices for hypercovariance is computationally easy. Considering the 4-level HNLM, Song et al. (2020) performed numerical experiment to compare the recommended prior with constant prior for the hypercovariance matrices, and the other two priors were canceled due to intractable computation. Also, Song et al. (2020)'s result presented the domination of the recommended hyperpriors over other priors. In conclusion, both Berger et al. (2020b) and Song et al. (2020) have provided strong numerical evidence of the superiority of the recommended hyperpriors for use in GHNLM, since 2-level and 4-level HNLMs both are specific GHNLM.

5. Discussions

We have proposed a generalized hierarchical normal linear model applicable to the nested data with complex structures. The GHNL model proves to be equivalent to a LMM model, while the GHNL model is more natural for researchers to model nested data from scratch, especially when incorporating covariates at high levels. Like generalizations to the simple normal linear model, the GHNL can be generalized to the hierarchical model with generalized liner model in the first level, and thus discrete observations can be handled. Besides, the first level (or even higher levels) of the GHNL model can be also extended to the setting of semiparametric regression models, such as, single index model and partially linear model. The technique of modelling and investigation in this paper can be applied to the linear part of the models mentioned above. The statistical analysis could be complicated and such explorations are beyond the scope of this paper, however.

Berger et al. (2020b) put an end to the endless search for the appropriate hyperpriors in hierarchical modelling and investigated properties comprehensively to justify the recommendation. Nonetheless, when it came to the propriety of the resulting posterior, they only suspected that that is true for use at any level for a general hierarchical normal model, the conditions for which were not given. To complete the story, we have studied the conditions for the posterior to be proper in more general situations than Berger et al. (2020b), when employing the recommended prior for the GHNL model. Theorems 3.1 and 3.3 demonstrate the main result, and Corollaries 3.1 and 3.3 reduce the computational burdens by defining kernel sets for $d\ge 2$ and d = 1, respectively. In addition, Theorems 3.2 and 3.4 provide powerful tools of simple forms for checking propriety of posterior for $d\ge 2$ and d = 1, separately. The user-friendly guidance for checking posterior propriety is eventually supplied. Note that our results only present sufficient conditions, and necessary conditions have never been discussed. The reason is because the derivation of the lower bound on the integrated likelihood of hyperparameters is intractable. Moreover, it is not worthwhile to investigate necessary conditions, as the derived upper bounds are tight enough such that the corresponding sufficient conditions are very modest, according to the remarks and examples in Section 3. At last, an efficient and powerful Gibbs sampling method for sampling from the posterior is introduced, overcoming the bottleneck of computation that the previously proposed sampling method only works for low dimensions or moderate dimensions inefficiently. The numerical evidence supporting the superiority of the recommended prior for hierarchical models was presented in Berger et al. (2020b) and Song et al. (2020).

Though we have made much progress in the hierarchical linear modelling, a major obstacle to applying our results is that the variance component for the first level is supposed to be known, which can hardly be satisfactory in practice. If we assume an unknown covariance matrix ${\mathbf{\Sigma }}_0$ for the first level and specify prior (13(13) $\begin{array}{rl}\pi \left({\mathit{V}}_j\right)& \propto \frac{1}{|{\mathit{V}}_j{|}^{1-1/\left(2k_j\right)}\prod _{s<t}({\omega }_{js}-{\omega }_{jt})},\\ & {\mathit{V}}_j>0,j\in \left[r\right],\end{array}$(13) ) on it, the exponential term within the likelihood can not be dropped simply any longer when deriving the upper bound, otherwise, the resulting integral will be always infinite. The upper bound on the exponential term with respect to the eigenvalues of covariance matrices is very tricky to be obtained, and the condition for the integrability of the resulting integral remains to be further studied. Thus, the GHNL modelling with unknown ${\mathbf{\Sigma }}_0$ can be taken as a sequential study of this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The research was supported by the National Natural Science Foundation of China [grant number 11671146].

Appendices

Appendix 1. A special LMM

Consider a special LMM of the form: (A1) $\mathit{y}=\mathit{X}\text{β}+{\mathit{Z}}_1{\mathit{u}}_1+\cdots +{\mathit{Z}}_r{\mathit{u}}_r+\mathit{ϵ},$(A1) where $\mathit{y}$ denotes the observations and is an $n\times 1$ vector, $\text{β}$ is the vector of fixed effects and is an $p\times 1$ vector. For $i\in \left[r\right]$, ${\mathit{u}}_i$'s are $q_i\times 1$ vectors and represent the vectors of random effects, and ${\mathit{u}}_i$'s are assumed to be independently distributed as $N_{q_i}(\mathbf{0},{\mathit{W}}_i)$, where ${\mathit{W}}_i$'s are $q_i\times q_i$ positive definite matrices and unknown. $\mathit{X}$ is an $n\times p$ matrix, ${\mathit{Z}}_i$'s are $n\times q_i$ matrices, and $\mathit{X}$ and ${\mathit{Z}}_i$'s are known design matrices. $\mathit{ϵ}$ is the vector of random errors and distributed as $N_n(\mathbf{0},\mathbf{\Sigma })$, $\mathbf{\Sigma }$ is an $n\times n$ positive definite matrix and given.

It follows from Berger et al. (2020b) that we can assume independent priors on $(\text{β},{\mathit{W}}_1,\dots ,{\mathit{W}}_r)$ as (A2) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(p-1)/2}},\text{β}\in {\mathbb{R}}^p,\\ \pi \left({\mathit{W}}_j\right)& \propto \frac{1}{|{\mathit{W}}_j{|}^{1-1/\left(2q_j\right)}\prod _{1\le s<t\le q_j}({\nu }_{js}-{\nu }_{jt})},\\ & {\mathit{W}}_j>0,j\in \left[r\right],\end{array}$(A2) where ${\nu }_{j1}>{\nu }_{j2}>\cdots >{\nu }_{j{q}_j}>0$ are the ordered eigenvalues of ${\mathit{W}}_j$, $j\in \left[r\right]$. The prior on $\text{β}$ has a hierarchical structure of the form $\left(\text{β}|\tau \right)\sim N_d(\mathbf{0},\tau {\mathit{I}}_d)\mathrm{a}\mathrm{n}\mathrm{d}\left[\tau \right]\propto {\tau }^{-1/2}\exp (-\frac{1}{2\tau }).$The posterior propriety results for the special LMM (A1(A1) $\mathit{y}=\mathit{X}\text{β}+{\mathit{Z}}_1{\mathit{u}}_1+\cdots +{\mathit{Z}}_r{\mathit{u}}_r+\mathit{ϵ},$(A1) ) is displayed as follows. Firstly, let $\tau \triangleq {\nu }_{01}$ and $q_0=1$. Denote the index set of the variance scale or the eigenvalues of the covariance matrices by $F=\left\{\right(j,l)|j=0,1,\dots ,r,l\in [q_j\left]\right\}$, and $\mathcal{T}=\{D|D\subseteq F,D\ne \mathrm{\varnothing }\}$ represents the set of the non-empty subsets of F. Define that $c_{jl,s}=1_{\{l=s\}}$ for $j\in \left[r\right]$, $l\in \left[q_j\right]$ and $s\in \left[n\right]$.

Theorem A.1

Consider linear mixed effect model (A1(A1) $\mathit{y}=\mathit{X}\text{β}+{\mathit{Z}}_1{\mathit{u}}_1+\cdots +{\mathit{Z}}_r{\mathit{u}}_r+\mathit{ϵ},$(A1) ) with prior (A2(A2) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(p-1)/2}},\text{β}\in {\mathbb{R}}^p,\\ \pi \left({\mathit{W}}_j\right)& \propto \frac{1}{|{\mathit{W}}_j{|}^{1-1/\left(2q_j\right)}\prod _{1\le s<t\le q_j}({\nu }_{js}-{\nu }_{jt})},\\ & {\mathit{W}}_j>0,j\in \left[r\right],\end{array}$(A2) ) on $(\text{β},{\mathit{W}}_1,\dots ,{\mathit{W}}_r)$. Assume p>1, and then the posterior is proper if (A3) $\sum _{s=1}^n{1}_{\{1_{\left\{\right(0,1)\in D\}}{1}_{\{s\le p\}}+\sum _{j\ne 0,(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}$(A3) holds for any $D\in \mathcal{T}$.

The proof of Theorem A.1 is similar to that of Theorem 3.1 and is omitted here.

Fact A.1

When p>1, (A3(A3) $\sum _{s=1}^n{1}_{\{1_{\left\{\right(0,1)\in D\}}{1}_{\{s\le p\}}+\sum _{j\ne 0,(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}$(A3) ) holds for any $D\in \mathcal{T}$ iff (A4) $\sum _{j\in \left[r\right]}\frac{1}{q_j}<1\mathrm{a}\mathrm{n}\mathrm{d}p>1+\sum _{j\in \left[r\right]}\frac{min(p,q_j)}{q_j}.$(A4)

Proof.

For any $D\in \mathcal{T}$ and $(0,1)\notin D$, (A3(A3) $\sum _{s=1}^n{1}_{\{1_{\left\{\right(0,1)\in D\}}{1}_{\{s\le p\}}+\sum _{j\ne 0,(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}$(A3) ) is equivalent to (A5) $\sum _{s=1}^n{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}.$(A5) It can be deduced that inequality (A5(A5) $\sum _{s=1}^n{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}.$(A5) ) holds for any $D\in \mathcal{T}$ and $(0,1)\notin D$ iff (A6) $L>\sum _{j\in \left[r\right]}\frac{min(L,q_j)}{q_j},\mathrm{f}\mathrm{o}\mathrm{r}L\in \left[n\right],$(A6) which is equivalent to $\sum _{j\in \left[r\right]}\frac{1}{q_j}<1$ since $q_j\ge 1$ for $j\in \left[r\right]$.

Inequality (A3(A3) $\sum _{s=1}^n{1}_{\{1_{\left\{\right(0,1)\in D\}}{1}_{\{s\le p\}}+\sum _{j\ne 0,(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}$(A3) ) holds for any $D\in \mathcal{T}$ with $(0,1)\in D$ iff (A7) $L>1+\sum _{j\in \left[r\right]}\frac{min(L,q_j)}{q_j},\mathrm{f}\mathrm{o}\mathrm{r}L=p,\dots ,n.$(A7) Under the condition that $\sum _{j\in \left[r\right]}\frac{1}{q_j}<1$, (A7(A7) $L>1+\sum _{j\in \left[r\right]}\frac{min(L,q_j)}{q_j},\mathrm{f}\mathrm{o}\mathrm{r}L=p,\dots ,n.$(A7) ) is equivalent to $p>1+\sum _{j\in \left[r\right]}\frac{min(p,q_j)}{q_j}.$

Corollary A.1

Consider model (A1(A1) $\mathit{y}=\mathit{X}\text{β}+{\mathit{Z}}_1{\mathit{u}}_1+\cdots +{\mathit{Z}}_r{\mathit{u}}_r+\mathit{ϵ},$(A1) ) with prior (A2(A2) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(p-1)/2}},\text{β}\in {\mathbb{R}}^p,\\ \pi \left({\mathit{W}}_j\right)& \propto \frac{1}{|{\mathit{W}}_j{|}^{1-1/\left(2q_j\right)}\prod _{1\le s<t\le q_j}({\nu }_{js}-{\nu }_{jt})},\\ & {\mathit{W}}_j>0,j\in \left[r\right],\end{array}$(A2) ) on parameters. The posterior is proper if one of the following condition holds,

1. p>1 + r and $\sum _{j\in \left[r\right]}\frac{1}{q_j}<1$;

2. p>1 and $\sum _{j\in \left[r\right]}\frac{1}{q_j}<1-\frac{1}{p}$.

Proof.

Since $\sum _{j\in \left[r\right]}\frac{min(p,q_j)}{q_j}\le r\mathrm{a}\mathrm{n}\mathrm{d}\sum _{j\in \left[r\right]}\frac{min(p,q_j)}{q_j}\le p\sum _{j\in \left[r\right]}\frac{1}{q_j},$(a) and (b) follow from Fact A.1 directly.

Remark A.1

Consider model (A1(A1) $\mathit{y}=\mathit{X}\text{β}+{\mathit{Z}}_1{\mathit{u}}_1+\cdots +{\mathit{Z}}_r{\mathit{u}}_r+\mathit{ϵ},$(A1) ) with r = 1. The posterior using prior (A2(A2) $\begin{array}{rl}\pi \left(\text{β}\right)& \propto \frac{1}{(1+\Vert \text{β}{\Vert }^2{)}^{(p-1)/2}},\text{β}\in {\mathbb{R}}^p,\\ \pi \left({\mathit{W}}_j\right)& \propto \frac{1}{|{\mathit{W}}_j{|}^{1-1/\left(2q_j\right)}\prod _{1\le s<t\le q_j}({\nu }_{js}-{\nu }_{jt})},\\ & {\mathit{W}}_j>0,j\in \left[r\right],\end{array}$(A2) ) is prior if either (a) $p\ge 2$, $q_1\ge 3$ or (b) $p\ge 3$, $q_1\ge 2$ holds. If p = 1, the posterior propriety can hardly be satisfied, the reason of which is two-fold. First, inequality (A3(A3) $\sum _{s=1}^n{1}_{\{1_{\left\{\right(0,1)\in D\}}{1}_{\{s\le p\}}+\sum _{j\ne 0,(j,l)\in D}{c}_{jl,s}>0\}}>\sum _{(j,l)\in D}\frac{1}{q_j}$(A3) ) fails for $D=\left\{\right(0,1\left)\right\}$. Second, if we follow the thread of deriving the condition in Theorem 3.3, a sufficient condition can be derived as (A8) $\sum _{s=1}^n{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{q_j},$(A8) for any $D\in \mathcal{T}$ with $(0,1)\notin D$. However, inequality (A8(A8) $\sum _{s=1}^n{1}_{\{\sum _{(j,l)\in D}{c}_{jl,s}>0\}}>1_{\{\mathrm{\exists }j\in [r],(j,1)\in D\}}+\sum _{(j,l)\in D}\frac{1}{q_j},$(A8) ) does not hold for $D=\left\{\right(j,1\left)\right\}$, $j\in \left[r\right]$.

Appendix 2

Proof of lemmas in Section 3.1

Lemma A.1

min–max theorem, cf. Horn & Johnson, (2012): For an $n\times n$ symmetric matrix $\mathit{A}$ and a non-zero $n\times 1$ vector $\mathit{x}$, the Rayleigh quotient for $\mathit{A}$ and $\mathit{x}$ can be defined as $R(\mathit{A},\mathit{x})=\frac{〈\mathit{A}\mathit{x},\mathit{x}〉}{〈\mathit{x},\mathit{x}〉},$where $〈\cdot ,\cdot 〉$ denotes the Euclidean inner product. Then (A9) $\begin{array}{rl}{\lambda }_k\left(\mathit{A}\right)& =\underset{U}{max}\{\underset{\mathit{x}}{min}\left\{R\right(\mathit{A},\mathit{x})|\mathit{x}\in U,\\ & \Vert \mathit{x}\Vert \ne 0\}|\dim \left(U\right)=k\},k\in \left[n\right]\end{array}$(A9) where U denotes the linear subspace of ${\mathbb{R}}^n$. Especially, (A10) $\begin{array}{rl}{\lambda }_{max}\left(\mathit{A}\right)& =\underset{\mathit{x}}{max}\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}\mathrm{a}\mathrm{n}\mathrm{d}\\ {\lambda }_{min}& =\underset{\mathit{x}}{min}\left(\mathit{A}\right)\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}.\end{array}$(A10)

Lemma A.2

For $n\times n$ symmetric matrices ${\mathit{A}}_j$, $j\in \left[r\right]$, we have

• ${\lambda }_{max}\left(\sum _{j=1}^r{\mathit{A}}_j\right)\le \sum _{j=1}^r{\lambda }_{max}\left({\mathit{A}}_j\right)$;

• supposing ${\mathit{A}}_j\ge 0$, $j\in \left[r\right]$, then ${\lambda }_k\left(\sum _{j=1}^r{\mathit{A}}_j\right)\ge \frac{1}{r}\sum _{j=1}^r{\lambda }_k\left({\mathit{A}}_j\right),k\in \left[n\right].$

Proof.

For (a), given an $n\times 1$ non-zero vector $\mathit{x}$, by (A10(A10) $\begin{array}{rl}{\lambda }_{max}\left(\mathit{A}\right)& =\underset{\mathit{x}}{max}\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}\mathrm{a}\mathrm{n}\mathrm{d}\\ {\lambda }_{min}& =\underset{\mathit{x}}{min}\left(\mathit{A}\right)\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}.\end{array}$(A10) ), (A11) $\begin{array}{rl}R(\sum _{j=1}^r{\mathit{A}}_j,\mathit{x})& =\sum _{j=1}^{r}R({\mathit{A}}_j,\mathit{x})\le \sum _{j=1}^r\underset{{\mathit{x}}_j\ne \mathbf{0}}{max}R({\mathit{A}}_j,{\mathit{x}}_j)\\ & =\sum _{j=1}^r{\lambda }_{max}\left({\mathit{A}}_j\right).\end{array}$(A11) The proof for (a) is completed by using (A10(A10) $\begin{array}{rl}{\lambda }_{max}\left(\mathit{A}\right)& =\underset{\mathit{x}}{max}\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}\mathrm{a}\mathrm{n}\mathrm{d}\\ {\lambda }_{min}& =\underset{\mathit{x}}{min}\left(\mathit{A}\right)\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}.\end{array}$(A10) ) again.

For (b), it suffices to prove that for any $j\in \left[r\right]$ (A12) ${\lambda }_k\left(\sum _{j=1}^r{\mathit{A}}_j\right)\ge {\lambda }_k\left({\mathit{A}}_j\right),i\in \left[n\right].$(A12) Since ${\mathit{A}}_l\ge 0$, $l\in \left[r\right]$, for any $j\in \left[r\right]$, $\mathit{x}\in {\mathbb{R}}^n$ and $\mathit{x}\ne \mathbf{0}$, $R(\sum _{l=1}^r{\mathit{A}}_l,\mathit{x})=\sum _{l=1}^{r}R({\mathit{A}}_l,\mathit{x})\ge R({\mathit{A}}_j,\mathit{x}).$Minimize the both sides of the inequality above over $\{\mathit{x}|\mathit{x}\in U,\mathit{x}\ne \mathbf{0}\}$ first. Then take the maximum over $\{U|U\subseteq {\mathbb{R}}^n,\dim (U)=k\}$. (A12(A12) ${\lambda }_k\left(\sum _{j=1}^r{\mathit{A}}_j\right)\ge {\lambda }_k\left({\mathit{A}}_j\right),i\in \left[n\right].$(A12) ) can be easily obtained by using Lemma A.1.

A.1. Proof of Lemma 3.1

For (a), by Lemma A.2 (a), it suffices to prove that ${\lambda }_{max}\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)\le C_1{\lambda }_{max}\left({\mathit{A}}_j\right)$, $j\in \left[r\right]$. For any $j\in \left[r\right]$, applying (A10(A10) $\begin{array}{rl}{\lambda }_{max}\left(\mathit{A}\right)& =\underset{\mathit{x}}{max}\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}\mathrm{a}\mathrm{n}\mathrm{d}\\ {\lambda }_{min}& =\underset{\mathit{x}}{min}\left(\mathit{A}\right)\left\{R\right(\mathit{A},\mathit{x})|\Vert \mathit{x}\Vert \ne 0\}.\end{array}$(A10) ), yields $\begin{array}{rl}0<{\lambda }_{max}\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)& =R({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},{\mathit{x}}^{\ast })\\ & =\underset{\mathit{x}\ne \mathbf{0}}{max}R({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},\mathit{x}).\end{array}$It is obvious that ${\mathit{X}}_j^{\mathrm{\top }}{\mathit{x}}^{\ast }\ne \mathbf{0}$, otherwise, it will result in $R({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},{\mathit{x}}^{\ast })=0$, which contradicts. In addition, since $〈{\mathit{X}}_j^{\mathrm{\top }}{\mathit{x}}^{\ast },{\mathit{X}}_j^{\mathrm{\top }}{\mathit{x}}^{\ast }〉\le {\lambda }_{max}\left({\mathit{X}}_j^{\mathrm{\top }}{\mathit{X}}_j\right)〈{\mathit{x}}^{\ast },{\mathit{x}}^{\ast }〉\le C_1〈{\mathit{x}}^{\ast },{\mathit{x}}^{\ast }〉$, we have $\begin{array}{rl}R({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},{\mathit{x}}^{\ast })& \le C_{1}R({\mathit{A}}_j,{\mathit{X}}_j^{\mathrm{\top }}{\mathit{x}}^{\ast })\\ & \le C_{1}max\{R({\mathit{A}}_j,\mathit{z})|\mathit{z}\in {\mathbb{R}}^{p_j},\mathit{z}\ne \mathbf{0}\}\\ & =C_1{\lambda }_{max}\left({\mathit{A}}_j\right).\end{array}$Therefore, we have proved part (a).

For (b), we only need to prove that ${\lambda }_k\left(\mathit{H}\right)\ge \frac{C_2}{r}\sum _{j=1}^r{a}_{jk},k\in \left[n\right].$It follows from Lemma A.2 (b) that for any $k\in \left[n\right]$ ${\lambda }_k\left(\mathit{H}\right)\ge \frac{1}{r}\sum _{j=1}^r{\lambda }_k\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right).$Since $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{A}}_j\right)$, then ${\lambda }_k\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)=0$, $p_j<k\le n$. Thus, it remains to show that (A13) ${\lambda }_k\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)\ge C_2{\lambda }_k\left({\mathit{A}}_j\right),j\in \left[r\right],k\in \left[p_j\right].$(A13) Firstly, for any $j\in \left[r\right]$, we introduce a linear transformation ${\mathcal{L}}_j:{\mathbb{R}}^n\to {\mathbb{R}}^{p_j}$ defined as ${\mathcal{L}}_j\left(\mathit{\nu }\right)={\mathit{X}}_j^{\mathrm{\top }}\mathit{\nu }$, $\mathit{\nu }\in {\mathbb{R}}^n$. The kernel space of ${\mathcal{L}}_j$ is denoted by $\mathrm{K}\mathrm{e}\mathrm{r}\left({\mathcal{L}}_j\right)=\{\mathit{\nu }\in {\mathbb{R}}^n:{\mathcal{L}}_j(\mathit{\nu })=\mathbf{0}\}$. Since ${\mathit{X}}_j$ is of full column rank $p_j\le n$, then the dimension of the complementary space of $\mathrm{K}\mathrm{e}\mathrm{r}\left({\mathcal{L}}_j\right)$ is $p_j$, i.e., $\dim (\mathrm{K}\mathrm{e}\mathrm{r}({\mathcal{L}}_j{)}^{\perp })=p_j$. Thus, the mapping ${\mathcal{L}}_j:\mathrm{K}\mathrm{e}\mathrm{r}({\mathcal{L}}_j{)}^{\perp }\mapsto {\mathbb{R}}^{p_j}$ is a one-to-one mapping. For any $U\subseteq {\mathbb{R}}^{p_j}$ with $\dim \left(U\right)=k$, $k\in \left[p_j\right]$, define that ${\mathcal{L}}_j^{\ast }\left(U\right)=\{\mathit{\nu }\in \mathrm{K}\mathrm{e}\mathrm{r}({\mathcal{L}}_j{)}^{\perp }:{\mathcal{L}}_j\left(\mathit{\nu }\right)=\mathit{x},\mathit{x}\in U\}.$It is obvious that ${\mathcal{L}}_j^{\ast }\left(U\right)\subseteq \mathrm{K}\mathrm{e}\mathrm{r}({\mathcal{L}}_j{)}^{\perp }$ and $\dim \left({\mathcal{L}}_j^{\ast }\right(U\left)\right)=k$. For any $U\subseteq {\mathbb{R}}^{p_j}$ with $\dim \left(U\right)=k$ and any $\mathit{x}\in U$, there exists one and only one $\mathit{\nu }\in {\mathcal{L}}_j^{\ast }\left(U\right)$ such that ${\mathcal{L}}_j\left(\nu \right)=\mathit{x}$. Since $〈\mathit{x},\mathit{x}〉={\mathit{\nu }}^{\mathrm{\top }}\left({\mathit{X}}_j^{\mathrm{\top }}{\mathit{X}}_j\right)\mathit{\nu }\ge C_2〈\mathit{\nu },\mathit{\nu }〉$, we have (A14) $R({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},\mathit{\nu })\ge C_{2}R({\mathit{A}}_j,\mathit{x}).$(A14) It refers to Lemma A.1 that (A15) $\begin{array}{rl}{\lambda }_k\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)& =\underset{V}{max}\{\underset{\mathit{v}}{min}\left\{R\right({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},\mathit{\nu })|\mathit{\nu }\in V,\Vert \mathit{\nu }\Vert \ne 0\}\\ & |V\subseteq \mathrm{K}\mathrm{e}\mathrm{r}({\mathcal{L}}_j{)}^{\perp },\dim (V)=k\},\end{array}$(A15) for $k\in \left[p_j\right]$. Minimize the both sides of inequality (A14(A14) $R({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},\mathit{\nu })\ge C_{2}R({\mathit{A}}_j,\mathit{x}).$(A14) ) over $\{\mathit{x}|\mathit{x}\in U,\mathit{x}\ne \mathbf{0}\}$ first. Then take the maximum over $\{U|U\subseteq {\mathbb{R}}^{p_j},\dim (U)=k\}$, (A13(A13) ${\lambda }_k\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)\ge C_2{\lambda }_k\left({\mathit{A}}_j\right),j\in \left[r\right],k\in \left[p_j\right].$(A13) ) can be easily obtained by using Lemma A.1 and (A15(A15) $\begin{array}{rl}{\lambda }_k\left({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }}\right)& =\underset{V}{max}\{\underset{\mathit{v}}{min}\left\{R\right({\mathit{X}}_j{\mathit{A}}_j{\mathit{X}}_j^{\mathrm{\top }},\mathit{\nu })|\mathit{\nu }\in V,\Vert \mathit{\nu }\Vert \ne 0\}\\ & |V\subseteq \mathrm{K}\mathrm{e}\mathrm{r}({\mathcal{L}}_j{)}^{\perp },\dim (V)=k\},\end{array}$(A15) ).

A.2. Proof of Lemma 3.2

The Domain of the integral can be divided into $\begin{array}{rl}& {\mathrm{\Omega }}_0=\{\mathit{\lambda }|0\le {\lambda }_j\le 1,j\in [k\left]\right\},\\ & {\mathrm{\Omega }}_D=\{\mathit{\lambda }|{\lambda }_j>1,j\in D\mathrm{a}\mathrm{n}\mathrm{d}0\le {\lambda }_i\le 1,i\in [k\left]/D\right\},\\ & D\in \mathcal{F}\left(\right[k\left]\right)\end{array}$i.e., $\mathrm{\Omega }=\left(\bigcup _{D\in \mathcal{F}\left(\right[k\left]\right)}{\mathrm{\Omega }}_D\right)\bigcup {\mathrm{\Omega }}_0$. Thus, the integral I is finite iff the integrals over ${\mathrm{\Omega }}_0$ and ${\mathrm{\Omega }}_D$ for each $D\in \mathcal{F}\left(\right[k\left]\right)$ are finite.

Denote the integrand as $F\left(\mathit{\lambda }\right)$. Then ${\int }_{{\mathrm{\Omega }}_0}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\simeq {\int }_{{\mathrm{\Omega }}_0}\frac{1}{\prod _{j=1}^k{\lambda }_j^{a_j}}\mathrm{d}\mathit{\lambda },$which is finite iff condition (a) is satisfied.

To verify condition (b), we only need to justify the following statement. Also, we assume condition (a) is always satisfied hereafter.

Fact A.2

For all $D\in \mathcal{F}\left(\right[k\left]\right)$ with $\text{card}\left(D\right)=L$, $1\le L\le k$, the integrals ${\int }_{{\mathrm{\Omega }}_D}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }$ are finite iff inequalities (A16) $\sum _{j\in E}{a}_j+\sum _{i\in {\mathcal{G}}_E}{b}_i>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(E\right)$(A16) hold for all $E\in \mathcal{F}\left(\right[k\left]\right)$ with $\text{card}\left(E\right)\le L$ and ${\mathcal{G}}_E=\{i|E\bigcap C_i\ne \mathrm{\varnothing },i\in \left[n\right]\}$. Under the condition above, (A17) $\begin{array}{rl}& {\int }_{{\mathrm{\Theta }}_D\left(t\right)}{G}_D\left(\mathit{\lambda }\right)\left(\prod _{r\in D}\mathrm{d}{\lambda }_r\right)\\ & \simeq \exp \{-\log t(\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i-\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(D\left)\right)\}\end{array}$(A17) always holds, where $\begin{array}{rl}{\mathrm{\Theta }}_D\left(t\right)& =\{\mathit{\lambda }|{\lambda }_j\ge t,j\in D\mathrm{a}\mathrm{n}\mathrm{d}0\le {\lambda }_i\le 1,i\in [k\left]/D\right\},\\ G_D\left(\mathit{\lambda }\right)& =\frac{1}{[\prod _{j\in D}{\lambda }_j^{a_j}][\prod _{i\in {\mathcal{G}}_D}(1+\sum _{r\in C_i}{\lambda }_r{)}^{b_i}]}.\end{array}$

The reason why formula (A17(A17) $\begin{array}{rl}& {\int }_{{\mathrm{\Theta }}_D\left(t\right)}{G}_D\left(\mathit{\lambda }\right)\left(\prod _{r\in D}\mathrm{d}{\lambda }_r\right)\\ & \simeq \exp \{-\log t(\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i-\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(D\left)\right)\}\end{array}$(A17) ) is required is that it plays an important role in verifying condition (A16(A16) $\sum _{j\in E}{a}_j+\sum _{i\in {\mathcal{G}}_E}{b}_i>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(E\right)$(A16) ).

Proof.

We prove the result by the technique of mathematical induction. First, we assume that the statement in Fact A.2 is true for L = l, $1\le l\le (k-1)$. With this assumption, we must show that the statement is true for its successor, $L=(l+1)$. Write an arbitrary set $D\in \mathcal{F}\left(\right[k\left]\right)$ with cardinality $(l+1)$ as $\{j_1,\dots ,j_{l+1}\}$, where $1\le j_1<\cdots <j_{l+1}\le k$. Denote that $D_{-j_i}=D/\left[\right\{j_i\left\}\right]$, $i=1,\dots ,(l+1)$.

Step 1: We first prove that ${\int }_{{\mathrm{\Omega }}_D}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }$ is finite iff inequalities (A16(A16) $\sum _{j\in E}{a}_j+\sum _{i\in {\mathcal{G}}_E}{b}_i>\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}\left(E\right)$(A16) ) hold for $L=(l+1)$.

Region ${\mathrm{\Omega }}_D$ can be divided into $\begin{array}{rl}{\mathrm{\Sigma }}_1& =\{\mathit{\lambda }|{\lambda }_j\ge {\lambda }_{j_1}>1,j\in D_{-j_1}\mathrm{a}\mathrm{n}\mathrm{d}\\ & 0\le {\lambda }_i\le 1,i\in \left[k\right]/D\}\\ & ⋮\\ {\mathrm{\Sigma }}_{l+1}& =\{\mathit{\lambda }|{\lambda }_j\ge {\lambda }_{j_{l+1}}>1,j\in D_{-j_{l+1}}\\ & \mathrm{a}\mathrm{n}\mathrm{d}0\le {\lambda }_i\le 1,i\in \left[k\right]/D\}.\end{array}$Therefore, integral ${\int }_{{\mathrm{\Omega }}_D}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }$ is finite iff ${\int }_{{\mathrm{\Sigma }}_i}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }<\mathrm{\infty }$ for any $i=1,\dots ,(l+1)$. For $i=1,\dots ,(l+1)$, we have ${\int }_{{\mathrm{\Sigma }}_i}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\simeq {\int }_{{\mathrm{\Sigma }}_i}\frac{1}{\prod _{s\notin D}{\lambda }_s^{a_s}}{\left(\frac{1}{{\lambda }_{j_i}}\right)}^{a_{j_i}+\sum _{r\in {\mathcal{H}}_i}{b}_r}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda },$where ${\mathcal{H}}_i=\{r|C_r\bigcap D_{j_i}=\mathrm{\varnothing }\mathrm{a}\mathrm{n}\mathrm{d}{j}_i\in C_r,r\in [n\left]\right\}$, and it's easy to see that ${\mathcal{G}}_D={\mathcal{G}}_{D_{-j_i}}\bigcup {\mathcal{H}}_i$ and ${\mathcal{G}}_{D_{-j_i}}\bigcap {\mathcal{H}}_i=\mathrm{\varnothing }$. Since (A18) $\begin{array}{rl}& {\int }_{{\mathrm{\Sigma }}_i}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right)\\ & \simeq {\int }_{{\mathrm{\Theta }}_{D_{-j_i}}\left({\lambda }_{j_i}\right)}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right),\end{array}$(A18) (A19) $\begin{array}{rl}& {\int }_{{\mathrm{\Omega }}_{D_{-j_i}}}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\simeq {\int }_{{\mathrm{\Omega }}_{D_{-j_i}}}\frac{1}{\prod _{s\notin D_{-j_i}}{\lambda }_s^{a_s}}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\end{array}$(A19) and the RHS (Right-Hand Side) of (A19(A19) $\begin{array}{rl}& {\int }_{{\mathrm{\Omega }}_{D_{-j_i}}}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\simeq {\int }_{{\mathrm{\Omega }}_{D_{-j_i}}}\frac{1}{\prod _{s\notin D_{-j_i}}{\lambda }_s^{a_s}}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\end{array}$(A19) ) and (A18(A18) $\begin{array}{rl}& {\int }_{{\mathrm{\Sigma }}_i}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right)\\ & \simeq {\int }_{{\mathrm{\Theta }}_{D_{-j_i}}\left({\lambda }_{j_i}\right)}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right),\end{array}$(A18) ) are finite simultaneously under condition (a). Hence, The LHS (Left-Hand Side) of (A18(A18) $\begin{array}{rl}& {\int }_{{\mathrm{\Sigma }}_i}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right)\\ & \simeq {\int }_{{\mathrm{\Theta }}_{D_{-j_i}}\left({\lambda }_{j_i}\right)}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right),\end{array}$(A18) ) is finite iff ${\int }_{{\mathrm{\Omega }}_{D_{-j_i}}}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }$ is finite.

Furthermore, by assumption and (A17(A17) $\begin{array}{rl}& {\int }_{{\mathrm{\Theta }}_D\left(t\right)}{G}_D\left(\mathit{\lambda }\right)\left(\prod _{r\in D}\mathrm{d}{\lambda }_r\right)\\ & \simeq \exp \{-\log t(\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i-\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(D\left)\right)\}\end{array}$(A17) ), we have $\begin{array}{rl}& {\int }_{{\mathrm{\Sigma }}_i}{G}_{D_{-j_i}}\left(\mathit{\lambda }\right)\left(\prod _{r\in D_{-j_i}}\mathrm{d}{\lambda }_r\right)\\ & \simeq \exp \{-\log {\lambda }_{j_i}(\sum _{j\in D_{-j_i}}{a}_j+\sum _{r\in {\mathcal{G}}_{D_{-j_i}}}{b}_r-l)\}.\end{array}$Thus, under condition (a) and assumption, we have ${\int }_{{\mathrm{\Sigma }}_i}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\simeq {\int }_1^{\mathrm{\infty }}{\left(\frac{1}{{\lambda }_{j_i}}\right)}^{\sum _{j\in D}{a}_j+\sum _{r\in {\mathcal{G}}_D}{b}_r-l}\mathrm{d}{\lambda }_{j_i},$the RHS of which is finite iff $\sum _{j\in D}{a}_j+\sum _{r\in {\mathcal{G}}_D}{b}_r>1+l=\text{card}\left(D\right)$.

In conclusion, ${\int }_{{\mathrm{\Omega }}_D}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }$ is finite iff $\sum _{j\in D}{a}_j+\sum _{r\in {\mathcal{G}}_D}{b}_r>\text{card}\left(D\right)$ and ${\int }_{{\mathrm{\Omega }}_{D_{-j_i}}}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }$ is finite for any $i\in [l+1]$. Since D is arbitrary and $\text{card}\left(D_{-j_i}\right)=l$, we have accomplished the goal of Step 1.

Step 2: Next, we prove that formula (A17(A17) $\begin{array}{rl}& {\int }_{{\mathrm{\Theta }}_D\left(t\right)}{G}_D\left(\mathit{\lambda }\right)\left(\prod _{r\in D}\mathrm{d}{\lambda }_r\right)\\ & \simeq \exp \{-\log t(\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i-\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(D\left)\right)\}\end{array}$(A17) ) holds for D with cardinality $(l+1)$.

Region ${\mathrm{\Theta }}_D\left(t\right)$ can be divided into $\begin{array}{rl}{\mathrm{\Theta }}_D^{\left(1\right)}\left(t\right)& =\{\mathit{\lambda }|{\lambda }_j\ge {\lambda }_{j_1}\ge t,j\in D_{-j_1}\mathrm{a}\mathrm{n}\mathrm{d}\\ & 0\le {\lambda }_i\le 1,i\in \left[k\right]/D\}\\ & ⋮\\ {\mathrm{\Theta }}_D^{(l+1)}\left(t\right)& =\{\mathit{\lambda }|{\lambda }_j\ge {\lambda }_{j_{l+1}}\ge t,j\in D_{-j_{l+1}}\mathrm{a}\mathrm{n}\mathrm{d}\\ & 0\le {\lambda }_i\le 1,i\in \left[k\right]/D\}.\end{array}$Similar to the proof of Step 1, we can prove that for $i=1,\dots ,(l+1)$, $\begin{array}{rl}& {\int }_{{\mathrm{\Theta }}_D^{\left(i\right)}\left(t\right)}{G}_D\left(\mathit{\lambda }\right)\left(\prod _{r\in D}\mathrm{d}{\lambda }_r\right)\\ & \simeq {\int }_t^{\mathrm{\infty }}{\left(\frac{1}{{\lambda }_{j_i}}\right)}^{\sum _{j\in D}{a}_j+\sum _{r\in {\mathcal{G}}_D}{b}_r-l}\mathrm{d}{\lambda }_{j_i}\\ & =\exp \{-\log t(\sum _{j\in D}{a}_j+\sum _{r\in {\mathcal{G}}_D}{b}_r-\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(D\left)\right)\}.\end{array}$Therefore, we get Step 2 proved.

Step 3: We need to present that the statement is true for L = 1 to complete the proof, on the basis of mathematical induction.

Denote that $D=\left\{r\right\}$, $r=1,\dots ,k$. Then ${\int }_{{\mathrm{\Omega }}_D}F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\simeq {\int }_1^{\mathrm{\infty }}{\left(\frac{1}{{\lambda }_r}\right)}^{a_r+\sum _{i\in {\mathcal{G}}_D}{b}_i}\mathrm{d}{\lambda }_r,$which is finite iff $\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i>1=\text{card}\left(D\right)$, under which, $\begin{array}{rl}& {\int }_{{\mathrm{\Theta }}_D\left(t\right)}{G}_D\left(\mathit{\lambda }\right)\left(\prod _{i\in D}\mathrm{d}{\lambda }_i\right)\\ & \simeq {\int }_t^{\mathrm{\infty }}{\left(\frac{1}{{\lambda }_r}\right)}^{a_r+\sum _{i\in {\mathcal{G}}_D}{b}_i}\mathrm{d}{\lambda }_r\\ & =\exp \{-\log t(\sum _{j\in D}{a}_j+\sum _{i\in {\mathcal{G}}_D}{b}_i-\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{d}(D\left)\right)\},\end{array}$which accomplishes the proof of Fact A.2.

Appendix 3

Gibbs sampling from the SIW distributions

As for the SIW distribution (46(46) ${\pi }^{\mathrm{S}\mathrm{I}\mathrm{W}}(\mathit{W}|a,\mathit{H})\propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}{\mathit{W}}^{-1}\mathit{H})}{{\left|\mathit{W}\right|}^a\prod _{i<j}({\nu }_i-{\nu }_j)},$(46) ), we first consider the change of variables from $\mathit{W}$ to $\mathbf{\Xi }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }_1,\dots ,{\nu }_k)$ and the orthogonal matrix $\mathit{O}$ of corresponding eigenvectors. The Jacobian is (A20) $\left|\frac{\mathrm{\partial }\mathit{W}}{\mathrm{\partial }(\mathbf{\Xi },\mathit{O})}\right|=\prod _{i<j}({\nu }_i-{\nu }_j).$(A20) According to (A20(A20) $\left|\frac{\mathrm{\partial }\mathit{W}}{\mathrm{\partial }(\mathbf{\Xi },\mathit{O})}\right|=\prod _{i<j}({\nu }_i-{\nu }_j).$(A20) ) and Lemma 4 in Berger et al. (2020a), (46(46) ${\pi }^{\mathrm{S}\mathrm{I}\mathrm{W}}(\mathit{W}|a,\mathit{H})\propto \frac{\mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}{\mathit{W}}^{-1}\mathit{H})}{{\left|\mathit{W}\right|}^a\prod _{i<j}({\nu }_i-{\nu }_j)},$(46) ) can be transformed to (A21) $\pi (\mathbf{\Xi },\mathit{O})\propto \frac{1}{{\left|\mathbf{\Xi }\right|}^a}\mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}{\mathbf{\Xi }}^{-1}{\mathit{O}}^{\prime }\mathit{H}\mathit{O}).$(A21) Gibbs sampling of $\mathbf{\Xi }$: We first sample $\mathbf{\Xi }$ given $(\mathit{O},\mathit{H})$ from $\begin{array}{rl}\pi (\mathbf{\Xi }|\mathit{O},\mathit{H})& \propto \frac{1}{\prod _{i=1}^k{\nu }_i^a}\mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}{\mathbf{\Xi }}^{-1}{\mathit{O}}^{\prime }\mathit{H}\mathit{O})\\ & =\prod _{i=1}^k\frac{1}{{\nu }_i^a}\mathrm{e}\mathrm{t}\mathrm{r}(-\frac{c_i}{{\nu }_i}),\end{array}$where $c_i$ is the i-th diagonal element of ${\mathit{O}}^{\prime }\mathit{H}\mathit{O}/2$, $i\in \left[k\right]$. Therefore, we can sample ${\nu }_i$ independently from $IG(a-1,c_i)$.

Gibbs sampling of $\mathit{O}$: Given $(\mathbf{\Xi },\mathit{H})$, the marginal density of $\mathit{O}$ has the form: $\pi (\mathit{O}|\mathbf{\Xi },\mathit{H})\propto \mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}\mathit{H}\mathit{O}{\mathbf{\Xi }}^{-1}{\mathit{O}}^{\prime }).$Let $\mathit{H}=\mathit{L}\mathit{U}{\mathit{L}}^{\mathrm{\top }}$, where $\mathit{L}{\mathit{L}}^{\mathrm{\top }}={\mathit{I}}_k$ and $\mathit{U}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(u_1,\dots ,u_k)$ is the diagonal matrix of corresponding eigenvalues with $u_1\ge \cdots \ge u_k$. Define $\mathit{G}={\mathit{L}}^{\mathrm{\top }}\mathit{O}$. Since the invariant right Haar measure is invariant to the orthonormal transformation, the conditional density of $\mathit{G}$ is (A22) $\pi (\mathit{G}|\mathbf{\Xi },\mathit{H})\propto \mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}\mathit{U}\mathit{G}{\mathbf{\Xi }}^{-1}{\mathit{G}}^{\prime }).$(A22) The updating of $\mathit{G}$ from (A22(A22) $\pi (\mathit{G}|\mathbf{\Xi },\mathit{H})\propto \mathrm{e}\mathrm{t}\mathrm{r}(-\frac{1}{2}\mathit{U}\mathit{G}{\mathbf{\Xi }}^{-1}{\mathit{G}}^{\prime }).$(A22) ) can be implemented by applying a Gibbs update to two randomly selected columns (cf. Hoff, 2009b) or rows (cf. Berger et al., 2020a). The two ways are essentially equivalent when $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{H}\right)=k$, but Berger et al. (2020a)'s method is considerably faster if $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{H}\right)<k$. Without any loss, assume that the two randomly selected rows are the first and second rows. The updated value of $\mathit{G}$ can be written as ${\mathit{G}}^{\mathrm{n}\mathrm{e}\mathrm{w}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{\Phi },{\mathit{I}}_{k-2})\left(\begin{array}{c}{\mathit{G}}_{12}^{\mathrm{o}\mathrm{l}\mathrm{d}}\\ {\mathit{G}}_{-12}^{\mathrm{o}\mathrm{l}\mathrm{d}}\end{array}\right)$, where ${\mathit{G}}_{12}^{\mathrm{o}\mathrm{l}\mathrm{d}}$ denotes the first two rows of the old value of $\mathit{G}$ which is ${\mathit{G}}^{\mathrm{o}\mathrm{l}\mathrm{d}}$, ${\mathit{G}}_{-12}^{\mathrm{o}\mathrm{l}\mathrm{d}}$ is the remaining k−2 rows of ${\mathit{G}}^{\mathrm{o}\mathrm{l}\mathrm{d}}$ and $\mathbf{\Phi }={\mathit{D}}_{ϵ}{\mathbf{\Phi }}_0=\left(\begin{array}{cc}{ϵ}_1& 0\\ 0& {ϵ}_2\end{array}\right)\left(\begin{array}{cc}\cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{array}\right),$with $\varphi \in (-\frac{\text{π}}{2},\frac{\text{π}}{2}]$ and ${ϵ}_i=\pm 1$ for i = 1, 2. Let ${\mathit{U}}_1=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(u_1,u_2)$. The full conditional density of ϕ has the form: $\pi (\varphi |{\mathit{G}}^{\mathrm{o}\mathrm{l}\mathrm{d}},\mathbf{\Xi },\mathit{H})\propto \mathrm{e}\mathrm{t}\mathrm{r}\{-\frac{1}{2}{\mathit{U}}_1{\mathbf{\Phi }}_0{\mathit{G}}_{12}^{\mathrm{o}\mathrm{l}\mathrm{d}}{\mathbf{\Xi }}^{-1}{\left({\mathit{G}}_{12}^{\mathrm{o}\mathrm{l}\mathrm{d}}\right)}^{\mathrm{\top }}{\mathbf{\Phi }}_0^{\mathrm{\top }}\}.$Write $\begin{array}{rl}{\mathit{G}}_{12}^{\mathrm{o}\mathrm{l}\mathrm{d}}{\mathbf{\Xi }}^{-1}{\left({\mathit{G}}_{12}^{\mathrm{o}\mathrm{l}\mathrm{d}}\right)}^{\mathrm{\top }}& =\left(\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{cc}{s}_1& 0\\ 0& s_2\end{array}\right)\\ & \times \left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right)\end{array}$where $\theta \in (-\frac{\text{π}}{2},\frac{\text{π}}{2}]$ and $s_1>s_2$. Then the conditional density of ϕ can be rewritten as $\pi (\varphi |{\mathit{G}}^{\mathrm{o}\mathrm{l}\mathrm{d}},\mathbf{\Xi },\mathit{H})\propto \exp \{-c_0{\cos }^2(\varphi +\theta \left)\right\},$where $c_0=\frac{1}{2}(s_1-s_2)(u_1-u_2)\ge 0$. Define $\alpha ={\cos }^2(\varphi +\theta )$. Then the full conditional density of α has the form: $\begin{array}{rl}\pi (\alpha |{\mathit{G}}^{\mathrm{o}\mathrm{l}\mathrm{d}},\mathbf{\Xi },\mathit{H})& \propto \exp \{-c_0\alpha \}{\alpha }^{-\frac{1}{2}}(1-\alpha {)}^{-\frac{1}{2}},\\ & \alpha \in [0,1].\end{array}$Simulating $\alpha \in [0,1]$ can proceed with a rejection sampler by setting the proposal distribution as $\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(\frac{1}{2},\frac{1}{2})$.