Bayesian Analysis for Hidden Markov Factor Analysis Models

Abstract

The purpose of this chapter is to provide an introduction to Bayesian approach within a general framework and develop a Bayesian procedure for analyzing multivariate longitudinal data within the hidden Markov factor analysis framework.

1. Introduction

The Bayesian approach is now well recognized in the statistics literature as an attractive approach to analyzing a wide variety of models [1], and there is rich literature on this issue. Here, we are not going to present a full coverage on the general Bayesian theory, and readers may refer to excellent books, for example [2, 3], for more details for this general statistical method. This chapter provides an introduction to the Bayesian approach within a general framework and develops a specific Bayesian procedure for analyzing multivariate longitudinal data within the hidden Markov factor analysis framework. We begin with the basic ideas of the Bayesian approach and then describe the model under consideration in the second section. The following section considers Bayesian inferences including parameter estimation, model selection, and posterior density estimates. The final section demonstrates the practical value of proposed methodology to cocaine use data to get some Bayesian results. Some technical details are given in the Appendix.

Consider a data set $Y$ with the probability model $p\left(Y|\theta \right)$ where $\theta$ is a univariate or multivariate population parameters vector, which quantifies the uncertainty of data. In the statistical literature, $p\left(Y|\theta \right)$ is called likelihood or sampling distribution and often represented as $L\left(\theta \right)$. From the frequency statistics point of view, statistical inferences are carried out based on $L\left(\theta \right)$. In this case, $\theta$, though unknown, is treated as fixed. Unlike the frequency statistical inferences, the Bayesian approach for data analysis assumes that $\theta$ is random and has a distribution $\pi \left(\theta \right)$. This distribution, which represents the knowledge about $\theta$, is referred to as prior distribution or prior. When data are available, the information on $\theta$ is summarized within the posterior distribution or posterior, a conditional distribution $\theta$ given data, i.e.,

$$p\left(\theta |Y\right)=\frac{p\left(Y|\theta \right)\pi \left(\theta \right)}{p\left(Y\right)}\propto p\left(Y|\theta \right)p\left(\theta \right)$$E1

where $p\left(Y\right)=\int p\left(Y|\theta \right)\pi \left(\theta \right)\mathit{d\theta }$ is the marginal distribution of $Y$. The right-hand-side term in (1) omits the factor $p\left(Y\right)$ since given $Y$ it is a known constant. In Bayes literature, $p\left(Y|\theta \right)p\left(\theta \right)$ is also termed the unnormalized posterior. Analogous to the role of likelihood in frequency statistical inferences, posterior is the starting point of Bayesian inferences.

Selecting proper priors for parameters is fundamental to Bayesian analysis. Basically, there are two kinds of prior distributions, namely, the noninformative prior distributions and the informative prior distributions. Noninformative prior distributions associate with situations when the prior distributions have no population basis. They are used when we have little prior information on $\theta$ and desire that the prior distributions play a minimal role in the posterior distribution distribution. Informative prior distribution represents the distribution of possible parameter values, from which the parameter $\theta$ has been drawn. We may have prior knowledge about this distribution, either from closed related data or from the subjective knowledge of experts. A commonly used informative prior distribution in the general Bayesian approach to statistical problems is the conjugate prior distribution, a prior ensuring that the posterior distribution follows the same parametric form as the prior distribution [1, 3].

A potential difficulty underlying Bayesian inferences is the statistical computation when posterior distribution takes on the complicated form. This is particularly true in the situation where latent variables or other unobservable quantities are involved in the model, as discussed in this chapter. In such cases, statistical inferences usually recur to simulation-based methods. Among various sampling methods, Markov chains Monte Carlo methods (MCMC) provide powerful tools for simulating observations from posterior. The key to Markov chain simulation is to create a Markov sequence whose stationary distribution is a specified posterior $p\left(\theta |Y\right)$. Posterior inferences are carried out based on these simulated observations. There are many ways of constructing these Markov chains, but all of them, including the Gibbs sampler [4, 5], are special cases of the general framework of Metropolis et al. [6] and Hastings [7]. However, we do not intend to pursue this issue here, and details on simulation-based methods can be referenced to [2, 3, 8, 9].

In what follows, as an illustration, we will develop a Bayesian analysis procedure for multivariate data under longitudinal setting. Multivariate longitudinal or clustered data occur when multiple items are measured repeatedly over periods of time or across occasions. Under such setting, the primary interest is inference about the dependence of the multiple measurements and the temporal correlation resulting from the repeated measures on the same items. But more often, particular interest also focuses on exploring the potential heterogeneity of data and investigating its transition pattern over time. In these cases, hidden Markov latent variable model (HMLVM) [10–13] provides a feasible and unified framework to address these issues. HMLVM assumes that the overall model constitutes the observed process and the underlying hidden state process. The state process, as the convention in the classic HMM (see for example, [14–17]), is an univariate discrete process, which follows a first-order Markov chain, while the observed process, conditional on the state sequence, is an independent process with emission distribution specified via LVMs [18]. Hence, in this regard, HMLVM provides a unified way of describing the correlation of multiple items, temporal dependence, and heterogeneity among the data simultaneously. However, the current existing developments cited beforehand focus on the maximum likelihood analysis in which statistical inferences heavily depend on the asymptotic properties. As an illustration of Bayesian inferences on practical problems, in this chapter, we develop a Bayesian procedure to analyze cocaine use data within the hidden Markov factor analysis model framework. Compared to ML, a basic nice feature of a Bayesian approach is its flexibility to utilize useful prior information for achieving better results. Additionally, simulation-based Bayesian methods depend less on asymptotic theory and hence have the potential to produce more reliable results even with small samples.

2. Model description

2.1. Hidden Markov factor analysis model

Consider a set of multivariate longitudinal observations formed by $p$-dimensional observed vectors $y_{\mathit{it}}={\left(y_{\mathit{it}1},\dots ,y_{\mathit{itp}}\right)}^{⊺}$, which are recorded on $p$ items over periods of length $T$: $t=1,\cdots ,T$ across $N$ subjects: $i=1,\cdots ,N$. In the field of multivariate analysis, interest mainly focuses on exploring item dependence since measurements may be highly correlated arising from the multicollinearity problem. But more often, interest also concentrates on the heterogeneity resulting from the situation where the population of $y_{\mathit{it}}$ constitutes more than one component. This is particularly true in the situation where the data illustrate extreme behaviors such as multimodal and/or skewed characteristics. In these cases, a finite mixture factor analysis model (FMFAM) can provide a powerful tool to address these issues. Typically, FMFAM assumes that conditioning on an univariate discrete value state variable $z_{\mathit{it}}$ and an $m$-dimensional ($m<p$) continuous latent factor vector ${\mathit{\omega }}_{\mathit{it}}$, ${\mathbf{y}}_{\mathit{it}}$ are independent and distributed with a $p$-dimensional multivariate normal distribution, and meanwhile, given $z_{\mathit{it}}$, ${\mathit{\omega }}_{\mathit{it}}$ also follows an $m$-dimensional normal distribution, that is,

$$\left\{\begin{array}{l}\left({\mathbf{y}}_{\mathit{it}}|{\mathit{\omega }}_{\mathit{it}},z_{\mathit{it}}=r\right)\sim {\mathcal{N}}_p\left({\mathit{\mu }}_r+{\mathbf{\Lambda }}_r{\mathit{\omega }}_{\mathit{it}},{\mathbf{\Psi }}_{\mathit{\epsilon r}}\right)\\ \left({\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}}=r\right)\sim {\mathcal{N}}_m\left(\mathbf{0},{\mathbf{\Phi }}_r\right)\end{array}\right.$$E2

where ${\mathit{\mu }}_r={\left({\mu }_{r1},\dots ,{\mu }_{\mathit{rp}}\right)}^{⊺}$ is a $p$-dimensional intercept vector, which represents the baseline level of ${\mathbf{y}}_{\mathit{it}}$, ${\mathbf{\Lambda }}_r={\left({\mathbf{\Lambda }}_{r1}^{⊺},\dots ,{\mathbf{\Lambda }}_{\mathit{rp}}^{⊺}\right)}^{⊺}$ is a $p\times m$ factor loading matrix, ${\mathbf{\Psi }}_{\mathrm{ϵ}r}=\text{diag}\left\{{\Psi }_{\mathrm{ϵ}\mathit{kr}1},\cdots ,{\Psi }_{\mathrm{ϵ}\mathit{krp}}\right\}$ is a $p\times p$ diagonal matrix with the $j$th diagonal element ${\Psi }_{\mathrm{ϵ}\mathit{krj}}>0$, and ${\mathbf{\Phi }}_r$ is an $m\times m$ positive definite matrix.

Formulation given in (2) has two basic features: one is to characterize heterogeneity of population of ${\mathbf{y}}_{\mathit{it}}$ at the occasion level and the other is to establish the dependence among the multiple measurements. The heterogeneous population is specified via state-specific parameters contained in the model while the dependence between different measurements is identified via sharing the common factors in the manner of liner combinations. In particular, apart from explaining the idiosyncratic part of measurements, latent factors also characterize the association between any two measurements. As a matter of fact, one can show that the correlation coefficient between $y_{\mathit{itj}}$ and $y_{\mathit{itk}}$ at state $z_{\mathit{it}}$ is given by

$$\text{Corr}\left(y_{\mathit{itj}},y_{\mathit{itk}}|z_{\mathit{it}}=r\right)=\frac{\sum _{\mathrm{\ell }=1}^m\sum _{h=1}^m{\lambda }_{\mathit{rj}\mathrm{\ell }}{\lambda }_{\mathit{rkh}}{\mathrm{\Phi }}_{\mathit{rlh}}}{\sqrt{\sum _{\mathrm{\ell }=1}^m\sum _{h=1}^m{\lambda }_{\mathit{rj}\mathrm{\ell }}{\lambda }_{\mathit{rjh}}{\mathrm{\Phi }}_{r\mathrm{\ell }h}+{\Psi }_{\mathrm{ϵ}\mathit{rj}}}\sqrt{\sum _{\mathrm{\ell }=1}^m\sum _{h=1}^m{\lambda }_{\mathit{rk}\mathrm{\ell }}{\lambda }_{\mathit{rkh}}{\mathrm{\Phi }}_{r\mathrm{\ell }h}+{\Psi }_{\mathrm{ϵ}\mathit{rk}}}}$$E3

in which ${\lambda }_{\mathit{rjk}}$ is the $\left(j,k\right)$th element of ${\mathbf{\Lambda }}_r$ and ${\mathrm{\Phi }}_{r,\mathit{hk}}$ is the $\left(h,k\right)$th element in $\mathbf{\Phi }$, respectively. The strength of correlation is identified by the factor loadings and covariance of factors together. In the case when ${\mathit{\omega }}_{\mathit{it}}$ degenerates to zero (i.e., $\mathbf{\Phi }=\mathbf{0}$) or $\mathbf{\Lambda }=\mathbf{0}$, the association among items disappears and model (2) reduces to $p$-independent mean-variance models within cluster $r$. Hence, latent factors play a dominant role in characterizing association of multiple items. Note that, in actual applications, latent factors, though unobservable, often have their own physical interpretations. In psychology, for example, latent factors are often used to identify concepts such as treatment, temper, and anxiety, which are important within the framework of theoretical models. The measurements are just proxies for these unobserved concepts of interest. We will provide further interpretations in the real example.

The primary reason for collecting information on multiple occasions for each subject is that it allows investigation of change and/or temporal dependence over time within the subject. There exist various constructs for characterizing dynamic characteristics. A commonly used method is to construct proper dynamic structures for latent factors and establish dynamic factor models, see for example, [19–21]. An alternative choice we adopt here is specifying the joint distribution for state sequences. Following the common routine (see, for example, [22, 23]), we assume that each individual state sequence ${\mathbf{z}}_i=\left(z_{i1},\cdots ,z_{\mathit{iT}}\right)$ satisfies the following first-order hidden Markov model

$$p\left(z_i\right)=p\left(z_{i1}\right)\prod _{t=2}^{T}p\left(z_{\mathit{it}}|z_{i,t-1}\right)$$E4

where $p\left(z_{i1}\right)$ and $p\left(z_{\mathit{it}}|z_{i,t-1}\right)$ are, respectively, the initial distribution and transition probability given by

$$P\left(z_{i1}=r|\right)={\delta }_r,P\left(z_{\mathit{it}}=s|z_{i,t-1}=r\right)=Q_{\mathit{rs}}\left(r,s=1,\cdots ,S\right)$$E5

where $S$ is a positive integer, $\mathit{\delta }=\left({\delta }_1,\cdots ,{\delta }_S\right)$ is an $S\times 1$ vector satisfying ${\delta }_r\ge 0$ and ${\sum }_{r=1}^S{\delta }_r=1.0$, and $\mathbf{Q}=\left(Q_{\mathit{rs}}\right)$ is an $S\times S$ transition matrix with the $\left(r,s\right)$th entry being $Q_{\mathit{rs}}$, that is, $Q_{\mathit{rs}}\ge 0$ and ${\sum }_{s=1}^S{Q}_{\mathit{rs}}=1.0$ for $r=1,\cdots ,S$. Modeling state sequences into (5) allows us to explore the transition pattern of individuals across occasions exactly. For example, in the cocaine use data analysis, $z_{\mathit{it}}$ is often identified with the latent state of patient $i$ at time $t$, then $Q_{\mathit{rs}}$ specifies how individual $i$ being in state $r$ transfers to state $s$ on two successive occasions. Surely, we can relax the time-homogeneous assumption of transition probabilities by including relevant covariates to interpret the inhomogeneous transition behavior among observation data (see, for example, [12, 13, 16]) but at the expense of computational burden.

The current model defined in (2)–(5) provides a comprehensive framework for modeling the multivariate longitudinal data with the latent variables. It accommodates the dynamic behavior of observed sequences, heterogeneity of observed data at the occasion level, and dependence among the multiple items simultaneously. In particular, it makes sense to measure effects of latent factors on the manifest variables quantitatively.

Let $\mathbf{Y}$ be the collection of all observations, and $\mathbf{\Omega }$ be the set of corresponding factors. Denote $\mathbf{Z}=\left\{z_{\mathit{it}}:1\le i\le N,1\le t\le T\right\}$ be set of state variables. It follows from Eqs. (2), (4) and (5) that the joint sampling distribution of $\mathbf{Y},\mathbf{\Omega }$, and $\mathbf{Z}$ is given by

$$\begin{array}{ll}p\left(\mathbf{Y},\mathbf{\Omega },\mathbf{Z}|\mathit{\theta },\mathit{\delta },\mathit{Q}\right)& =\prod _{i=1}^{N}p\left({\mathbf{y}}_{i1},{\mathit{\omega }}_{\mathit{it}}|z_{i1},\mathit{\theta }\right)p\left(z_{i1}|\mathit{\delta }\right)\prod _{t=2}^{T}p\left({\mathbf{y}}_{\mathit{it}},{\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}},\mathit{\theta }\right)p\left(z_{\mathit{it}}|z_{i,t-1},\mathit{Q}\right)\\ & \propto \prod _{i=1}^N\prod _{t=1}^T\left(\frac{1}{{\left|{\mathbf{\Psi }}_{\mathrm{ϵ}{z}_{\mathit{it}}}\right|}^{1/2}}\text{exp}\left\{-\frac{1}{2}\mathrm{tr}{\mathbf{\Psi }}_{\mathrm{ϵ}{z}_{\mathit{it}}}{\left({\mathbf{y}}_{\mathit{it}}-{\mathit{\mu }}_{z_{\mathit{it}}}-{\mathbf{\Lambda }}_{z_{\mathit{it}}}^{⊺}{\mathit{\omega }}_{\mathit{it}}\right)}^{\otimes 2}\right\}\right.\hfill \\ & \times \left.\frac{1}{{\left|{\mathbf{\Phi }}_{z_{\mathit{it}}}\right|}^{1/2}}\text{exp}\left\{-\frac{1}{2}\mathrm{tr}{\mathbf{\Phi }}_{z_{\mathit{it}}}^{-1}{\mathit{\omega }}_{\mathit{it}}^{\otimes 2}\right\}\right)\times \prod _{i=1}^N\prod _{t=1}^T\left(\prod _{r=1}^S{\delta }_r^{I\left\{z_{i1}=r\right\}}\prod _{s=1}^S{Q}_{\mathit{rs}}^{I\left\{z_{i,t-1}=r,z_{\mathit{it}}=s\right\}}\right)\hfill \end{array}$$E6

where $\mathit{\theta }$ is formed by free parameters in ${\mathit{\mu }}_r,{\mathbf{\Lambda }}_r,{\mathbf{\Psi }}_r$, and ${\mathbf{\Phi }}_r$. Here, we write ${\mathbf{a}}^{\otimes 2}={\mathbf{aa}}^{⊺}$ and denote $I\left(A\right)$ the indicator function of a set $A$. The observed data likelihood is then achieved by taking integration of $p\left(\mathbf{Y},\mathbf{\Omega },\mathbf{Z}|\mathit{\theta },\mathit{\delta },\mathit{Q}\right)$ over $\mathbf{\Omega }$ and $\mathbf{Z}$, which involves high-dimensional integrations.

3. Posterior inferences

3.1. Prior specifications

Let $\mathit{\mu }=\left\{{\mathit{\mu }}_r\right\}$, $\mathbf{\Lambda }=\left\{{\mathbf{\Lambda }}_r\right\}$, ${\mathbf{\Psi }}_{\mathrm{ϵ}}=\left\{{\mathbf{\Psi }}_{\mathrm{ϵ}r}\right\}$, and $\mathbf{\Phi }=\left\{{\mathbf{\Phi }}_{\mathit{kr}}\right\}$. For the Bayesian analysis, we need to assign priors to the unknown parameters involved for completing model specification. Since $\mathit{\theta }$, $\mathit{\delta }$, and $\mathit{Q}$ are involved in different submodels, it is natural to assume that $\mathit{\theta },\mathit{\delta }$, and $\mathit{Q}$ are mutually independent and the components contained in $\mathit{\theta }$ are also mutually independent, that is,

$$p\left(\mathit{\theta },\mathit{\delta },\mathbf{Q}\right)=p\left(\mathit{\mu }\right)p\left(\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}\right)p\left(\mathbf{\Phi }\right)p\left(\mathit{\delta }\right)p\left(\mathbf{Q}\right).$$E7

For the convenience of conjugacy, we assume that the parameters are drawn from the following commonly used conjugate types prior distributions (see for example [24]).

$$\begin{array}{ll}p\left(\mathit{\mu }\right)\hfill & =\prod _{r=1}^{S}p\left({\mathit{\mu }}_r\right){\displaystyle \stackrel{D}{=}}\prod _{r=1}^S{\mathcal{N}}_p\left({\mathit{\mu }}_{0r},{\mathbf{\Sigma }}_{0r}\right),\mathbf{\Phi }\sim \prod _{r=1}^{S}p\left({\mathbf{\Phi }}_r\right){\displaystyle \stackrel{D}{=}}\prod _{r=1}^S{\mathcal{W}}_m^{-1}\left({\rho }_{0r},{\mathbf{R}}_{0r}^{-1}\right),\\ p\left(\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}\right)\hfill & =\prod _{r=1}^{S}p\left({\mathbf{\Lambda }}_r|{\mathbf{\Psi }}_{\mathrm{ϵ}r}\right)\times p\left({\mathbf{\Psi }}_{\mathrm{ϵ}r}\right) {\displaystyle \stackrel{D}{=}}\prod _{r=1}^S\prod _{j=1}^p{\mathcal{N}}_m\left({\mathbf{\Lambda }}_{0\mathit{rj}},{\mathrm{\psi }}_{\mathrm{ϵ}\mathit{rj}}{\mathbf{H}}_{\mathrm{ϵ}0\mathit{rj}}\right)\cdot \mathcal{G}{a}^{-1}\left({\alpha }_{\mathrm{ϵ}0\mathit{rj}},{\beta }_{\mathrm{ϵ}0\mathit{rj}}\right),\\ & \mathit{\delta }\mid {\delta }_0\sim \mathcal{D}{\mathit{ir}}_S\left({\gamma }_0,\dots ,{\gamma }_0\right),p\left(\mathbf{Q}\right)=\prod _{r=1}^{S}p\left({\mathbf{Q}}_r\right){\displaystyle \stackrel{D}{=}}\prod _{r=1}^S\mathcal{D}{\mathit{ir}}_S\left({\nu }_0,\cdots ,{\nu }_0\right)\hfill \end{array}$$E8

where ‘$\mathcal{G}{a}^{-1}\left(a,b\right)$’ denotes the inverse Gamma distributions with shape $a>0$ and scale $b>0$ and ‘${\mathcal{W}}_{m_2}^{-1}\left({\rho }_{0r},{\mathbf{R}}_{0r}^{-1}\right)$’ represents the $q$-dimensional inverse Wishart distribution with ${\rho }_{0r}$ degrees of freedom and $\left(m_2\times m_2\right)$ scale matrix ${\mathbf{R}}_{0r}$; ${\mathbf{Q}}_r$ is the $r$th row vector of $\mathbf{Q}$. The scalars ${\alpha }_{\mathrm{ϵ}0\mathit{rj}}$, ${\beta }_{\mathrm{ϵ}0\mathit{rj}}$, ${\rho }_{0r}$, ${\gamma }_0$, ${\nu }_0$, the vectors ${\mathit{\mu }}_{0r}$, ${\mathbf{\Lambda }}_{0\mathit{rj}}$, and the matrices ${\mathbf{R}}_{0r}$ and ${\mathbf{H}}_{\mathrm{ϵ}0\mathit{rj}}$ are assumed to be known. Thus, standard conjugate priors were specified for all parametric components in the model. The conjugate type prior distributions are sufficiently flexible in most applications, and for situations with a reasonable amount of data available, the hyperparameters scarcely affect the analysis. It should be noted that although Eq. (8) allows different hyperparameters for different latent states, in practice, we choose identical priors for all $s$. Details of hyperparameter choices are discussed later when we present the empirical results.

3.2. Gibbs sampling scheme and posterior analysis

Combining the sampling distribution for the observable ${\mathbf{y}}_{\mathit{it}}$’s and the prior distribution specified in (8) yields the joint posterior distribution of $\left\{\mathit{\theta },\mathit{\delta },\mathbf{Q}\right\}$ given by

$$p\left(\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)\propto p\left(\mathbf{Y}|\mathit{\theta },\mathit{\delta },\mathbf{Q}\right)p\left(\mathit{\theta }\right)p\left(\mathit{\delta }\right)p\left(\mathbf{Q}\right)$$E9

where we ignore the normalization constant $p\left(\mathbf{Y}\right)$. However, due to the latent factors and state variables present, the computation of $p\left(\mathbf{Y}|\mathit{\theta },\mathit{\delta },\mathbf{Q}\right)$ is intractable since it involves high-dimensional integrals. Consequently, no closed form can be available for the posterior $p\left(\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)$. This problem can be addressed via the data augmentation idea in Tanner and Wong [25]. Data augmentation technique treats the latent quantities $\left\{\mathbf{\Omega },\mathbf{Z}\right\}$ as the hypothetical missing data and augments them with the observed data to form complete data. The posterior analysis is now carried out based on the joint distribution $p\left(\mathbf{\Omega },\mathbf{Z},\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)$, which is proportional to $p\left(\mathbf{Y},\mathbf{\Omega },\mathbf{Z}|\mathit{\theta },\mathit{\delta },\mathbf{Q}\right)p\left(\mathit{\theta },\mathit{\delta },\mathbf{Q}\right)$, the product of likelihood of complete data and priors. Compared to the intractable observed data likelihood, the complete data likelihood has nice hierarchical structure based on conditional independent assumptions in (2) and (4) and hence is relatively easy to analyze. However, $p\left(\mathbf{\Omega },\mathbf{Z},\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)$ is still not in closed form and is thus difficult to deal with analytically. In this regard, simulation-based methods can be used to generate observations to carry out posterior analysis. In view of the multiple components involved, the usual independent sampling methods are not feasible. Note that, on the basis of complete data, the full conditional distributions of $\mathbf{\Omega },\mathbf{Z},\mathit{\theta },\mathit{\delta },$ and $Q$ have closed forms. This provides a solid foundation for Markov chain Monte Carlo methods. Markov chain Monte Carlo sampling does not draw observations from $p\left(\mathbf{\Omega },\mathbf{Z},\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)$ directly. On the contrary, it generates observations from the full conditionals of each component alternatively, thus forming the dependent sample, i.e., Markov chains. Specifically, as pointed out in the introduction, we use Gibbs sampler [4, 5] to draw observations from this target distribution. Obviously, the sampling scheme in the Gibbs sampler includes two types of moves: updating the components involved in the factor analysis model and updating the components related to the hidden Markov model. We propose using the following Gibbs sampler which iteratively simulates from the conditional distributions, where variables are removed from the conditioning set either by explicit integration or by conditional independence. The steps involved in the Gibbs sampler are

1. Step a: Generate $\mathbf{Z}$ from $p\left(\mathbf{Z}|\mathit{\theta },\mathit{\delta },\mathbf{Q},\mathbf{\Omega },\mathbf{Y}\right)$

2. Step b: Generate $\mathbf{\Omega }$ from $p\left(\mathbf{\Omega }|\mathbf{Z},\mathit{\theta },\mathbf{Y}\right)$

3. Step c: Generate $\left\{\mathit{\mu },\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}\right\}$ from $p\left(\mathit{\mu },\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}|\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)$

4. Step d: Generate $\mathbf{\Phi }$ from $p\left(\mathbf{\Phi }|\mathbf{Z},\mathbf{\Omega }\right)$

5. Step e: Generate $\mathit{\delta }$ from $p\left(\mathit{\delta }|\mathbf{Z}\right)$

6. Step f: Generate $\mathbf{Q}$ from $p\left(\mathbf{Q}|\mathbf{Z}\right)$

Under mild conditions and similar to [4] (see also, for example, [26]), one can show that for sufficiently large $b$, say $B_0$, the joint distribution of $\left\{{\mathbf{\Omega }}^{\left(b\right)},{\mathbf{Z}}^{\left(b\right)},{\mathit{\theta }}^{\left(b\right)},{\mathit{\delta }}^{\left(b\right)},{\mathbf{Q}}^{\left(b\right)}\right\}$ converges at an exponential rate to the desired posterior distribution $p\left(\mathbf{\Omega },\mathbf{Z},\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)$. Hence, $p\left(\mathbf{\Omega },\mathbf{Z},\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)$ can be approximated by the empirical distribution of $\left\{{\mathbf{\Omega }}^{\left(b\right)},{\mathbf{Z}}^{\left(b\right)},{\mathit{\theta }}^{\left(b\right)},{\mathit{\delta }}^{\left(b\right)},{\mathbf{Q}}^{\left(b\right)}\right\}:b=B_0+1,\cdots ,B_0+B\}$ where $B$ is chosen to give sufficient precision to the empirical distribution. The convergence of the Gibbs sampler can be monitored by the ‘estimated potential scale reduction (EPSR)’ values as suggested by Gelman and Rubin [27] or by plotting the traces of estimates against iterations under different starting values.

Simulated observations obtained from the posterior can be used for statistical inferences via straightforward analysis procedures. For brevity, let $\left\{{\mathit{\theta }}^{\left(b\right)},{\mathit{\delta }}^{\left(b\right)},{\mathbf{Q}}^{\left(b\right)},{\mathbf{\Omega }}^{\left(b\right)},{\mathbf{Z}}^{\left(b\right)}\right\}$ be the random observations generated by the Gibbs sampler from $p\left(\mathit{\theta },\mathit{\delta },\mathbf{Q},\mathbf{\Omega },\mathbf{Z}|\mathbf{Y}\right)$. The joint Bayesian estimate of $\mathit{\theta }$ and $\mathbf{\Omega }$ can be obtained easily via the corresponding sample means of the generated observations as follows:

$$\widehat{\mathit{\theta }}={\left(B-1\right)}^{-1}\sum _{b=1}^B{\mathit{\theta }}^{\left(b\right)},\widehat{\mathbf{\Omega }}={\left(B-1\right)}^{-1}\sum _{b=1}^B{\mathbf{\Omega }}^{\left(b\right)},\widehat{\mathbf{Z}}={\left(B-1\right)}^{-1}\sum _{b=1}^B{\mathbf{Z}}^{\left(b\right)}.$$E10

Clearly, these Bayesian estimates are consistent estimates of the corresponding posterior means, see [26]. The consistent estimates of covariance matrix of estimates can be obtained as follows:

$$\widehat{\mathrm{C}\mathrm{ov}(\mathbf{\theta }\mid }\mathbf{Y})={\left(B-1\right)}^{-1}\sum _{b=1}^B\left({\mathit{\theta }}^{\left(b\right)}-\widehat{\mathit{\theta }}\right){\left({\mathit{\theta }}^{\left(b\right)}-\widehat{\mathit{\theta }}\right)}^{⊺}$$E11

$$\mathrm{C}{\displaystyle \widehat{\mathrm{ov}(\mathbf{\Omega }\mid }}\mathbf{Y})={\left(B-1\right)}^{-1}\sum _{b=1}^B\left({\mathbf{\Omega }}^{\left(b\right)}-\widehat{\mathbf{\Omega }}\right){\left({\mathbf{\Omega }}^{\left(b\right)}-\widehat{\mathbf{\Omega }}\right)}^{⊺}$$E12

Hence, the standard error estimates can be obtained conveniently by the Gibbs sampler algorithm. Other statistical inferences about $\mathit{\theta }$ and $\mathbf{\Omega }$ such as deriving the confidence intervals and statistics for hypothesis testing can be achieved based on the simulative observations as well (see, for example, [28, 29]).

One important statistical inference beyond estimation is on testing of various hypotheses about the model. In the field of hidden Markov modeling, determining the proper number of states may be the first step towards data analysis. Too many states may overfit the observations, meaning that it can fit the training data accurately but may not be a good model for underlying data-generating process. On the other hand, too few states may not be flexible enough to approximate the underlying model. In the context of Bayesian model selection, Bayes factor (BF, [30]) is a popular choice for model comparison. BF is defined as the ratio of the marginal likelihoods of data under two competing models. However, the computation of BF is difficult since it often involves the high-dimensional integrations. It has also been shown that BF is sensitive to the choice of priors and will become infeasible when improper priors are used. A simple and more convenient alternative is the $L_{\nu }$-measures [31–34] which is based on the posterior predictive density. It has been shown [34] that this approach is conceptually and computationally simple and is useful in model checking for wide varieties of complicated situations. Moreover, the required computation is a by-product of the common Bayesian simulation procedures such as the Gibbs sampler or its related algorithms. Specifically, let ${\mathbf{Y}}^{\mathrm{rep}}$ denotes future values of $\mathbf{Y}$ in a replicate experiment, that is, ${\mathbf{Y}}^{\mathrm{rep}}$ has the same sampling density as that of $\mathbf{Y}$. The posterior predictive distribution $p\left({\mathbf{Y}}^{\mathrm{rep}}|\mathbf{Y}\right)$ is defined as

$$p\left({\mathbf{Y}}^{\mathrm{rep}}|\mathbf{Y}\right)=\int p\left({\mathbf{Y}}^{\mathrm{rep}}|\mathit{\theta },\mathit{\delta },\mathbf{Q}\right)p\left(\mathit{\theta },\mathit{\delta },\mathbf{Q}|\mathbf{Y}\right)d\mathit{\theta }d\mathit{\delta }d\mathbf{Q}$$E13

Naturally, if the posited model under consideration is the true model in the sense that from which the data are generated, then ${\mathbf{Y}}^{\mathrm{rep}}$ would behave like data $\mathbf{Y}$ and its squared biases and covariances should be small. With this notion in mind, Ibrahim, Chen, and Sinha [34] proposed an $L$ statistics to assess the fitness of posited models to the data by weighting the squared biases and covariance, which can be interpreted as a trade-off between them. Here, we extend it to the multivariate longitudinal setting. Let ${\mathbf{Y}}^{\mathrm{rep}}={\left({\mathbf{y}}_1^{\mathrm{rep}⊺},\cdots ,y_N^{\mathrm{rep}⊺}\right)}^{⊺}$ be a collection set of future responses in our proposal. For some $0\le \nu <1$, we consider the following multivariate version of $L_{\nu }$-measure:

$$L_{\nu }\left(\mathbf{Y}\right)=\sum _{i=1}^N\mathrm{tr}\left[\mathbb{C}\mathit{ov}\left({\mathbf{y}}_i^{\mathrm{rep}}|\mathbf{Y}\right)\right]+\nu \sum _{i=1}^N\mathrm{tr}\left[\left\{\mathbb{E}\left({\mathbf{y}}_i^{\mathrm{rep}}|\mathbf{Y}\right)-{\mathbf{y}}_i\right\}{\left\{\mathbb{E}\left({\mathbf{y}}_i^{\mathrm{rep}}|\mathbf{Y}\right)-{\mathbf{y}}_i\right\}}^{⊺}\right]$$E14

where the expectation is taken with respect to the posterior predictive distribution. Clearly, small values of the $L_{\nu }$-measure indicate that the model gives predictions close to the observed values, and the variability in the predictions is low as well. Hence, the model with the smallest $L_{\nu }$-measure is selected from a collection of competing models. It has been shown that $L_{\nu }$-measure with $\nu =0.5$ has nice theoretical properties [34]. Thus, this value of $\nu$ will be used in our empirical illustrations.

4. Cocaine use data analysis

In this section, a small portion of cocaine use data is analyzed to illustrate the practical value of the proposed methodology. The original data are collected from 321 cocaine use patients who were admitted in 1988–1989 to the West Los Angeles Veterans Affairs Medical Center. The whole data constitute 68 measurements of 17 items, which were recorded at four time points: at baseline, 1 year after the treatment, 2 years after the treatment, and 12 years after the treatment in 2002–2003. These measurements cover the information on the cocaine use, treatment received, psychological problems, social status, employments, and so on. As an illustration, three variables are selected to conduct data analysis: ‘$y_1:$ days of cocaine use per month at intake (CC)’, ‘$y_2:$ times per month in formal treatment (FT)’, and ‘$y_3:$ months in formal treatment (MFT)’, which, respectively, represent the severity of cocaine use and the levels of treatment received by a patient. Since these variables were measured in 0–120 points scale, to unify the scales, we take logarithms and standardize them. Among them, some measurements are missing. The missing proportion is about $8.4\%$. For brevity, we assume that the missing is missing at random [35]. A distinct characteristic underlying data are nonnormal and heavy tailed. Figure 1 gives the plots of histograms and the posterior predictive density estimates (see below) of logarithms of CC, FT, and MFT (with missing data removed) on four occasions. The histograms illustrate that the distributions of selected variables are deviated from normality in terms of multimodality and skewness. The skewness and kurtosis of CC on four occasions are $\left\{-1.631,5.031\right\}$, $\left\{-0.847,3.354\right\}$, $\left\{0.328,1.476\right\}$, and $\left\{-0.473,2.467\right\}$, respectively. Data set also demonstrates dynamic characteristics. The distribution of CC, for instance, is skewed to the left at baseline and moves to the right gradually on the following two occasions and becomes right-skewed eventually. This implies that a single factor analysis model may not be appropriate to fit the data at each time point.

Figure 1.

Plots of histograms and posterior predictive density estimates of ‘CC’, ‘FT’ and ‘MFT’ under FA model and hidden Markov CFA model with seven states in the cocaine use data analysis: the dashed lines denote CFA and the solid lines represent the hidden Markov FA.

In this analysis, one of the objectives is to explore the effects of latent factors on the observed variables and assess the dependence among latent factors. Based on the nature of the problem under consideration, it is natural to group the single variable ‘CC’ to reflect one latent factor ‘cocaine use’ ($\eta$) and to group ‘FT’ and ‘MFT’ to represent another latent factor ‘treatment’ ($\xi$). Let ${\mathbf{y}}_{\mathit{it}}={\left(y_{\mathit{it}1},y_{\mathit{it}2},y_{\mathit{it}3}\right)}^{⊺}$ and ${\mathit{\omega }}_{\mathit{it}}={\left({\eta }_{\mathit{it}},{\xi }_{\mathit{it}}\right)}^{⊺}$. To be convenient for interpretation and computation, ${\mathbf{\Phi }}_r$ and ${\mathbf{\Lambda }}_r$ are restricted to be invariant across states but leave the baseline level ${\mathit{\mu }}_r$ varying with $r$. Further, the following non-overlapped structure for factor loading matrix is considered

$${\mathbf{\Lambda }}^{⊺}=\left(\begin{array}{lll}{1}^{\ast }& 0^{\ast }& 0^{\ast }\\ 0^{\ast }& 1^{\ast }& {\mathrm{\Lambda }}_{32}\end{array}\right)$$E15

where parameters with an asterisk are treated as fixed for identification. Note that fixing ${\mathrm{\Lambda }}_{11}=1$ indicates that $\eta$ is identified with CC. This is similar to that in ${\mathrm{\Lambda }}_{22}$. Hence, in this case, ${\mathrm{\Phi }}_{12}$ in $\mathbf{\Phi }$ measures the magnitude of dependence of $\xi$ on $\eta$.

Data set is fitted to the proposed models with 10 different transition models: $S=1,\cdots ,10$. Although these state spaces are in nested forms, the corresponding models are not, since one cannot be reduced to another by constraining parameters in the interior of parameter space. This indicates that chi-square distribution may not be suitable for the classic likelihood ratio test statistic. We use $L$-measure to implement model selection. Obviously, if ${\mathbb{S}}_1$ is taken, then the proposed model reduces to common factor analysis model (CFA, [18]).

The following inputs are taken for the super-parameters involved in the prior distributions (8): for $r=1,\cdots ,S$, ${\mu }_{0\mathit{rj}}=\text{min}\left\{y_{\mathit{itj}}\right\}+r/S$, ${\mathbf{\Sigma }}_{0r}={\mathbf{S}}_{\mathit{yy}}/S$, where ${\mathbf{S}}_{\mathit{yy}}$ is the sample covariance matrix of data. The entries in ${\mathbf{\Lambda }}_0$ are set to be zeros, ${\rho }_0=10.0$, ${\mathbf{R}}_0^{-1}=7.0\times {\mathbf{I}}_2$, which leads to the mean of $\mathbf{\Phi }$ equal to ${\mathbf{I}}_2$, ${\mathbf{H}}_{\mathrm{ϵ}0}={\mathbf{I}}_3$, ${\alpha }_{\mathrm{ϵ}0j}=9.0$, ${\beta }_{\mathrm{ϵ}0j}=8.0$, ${\nu }_0={\gamma }_0=0.1$. Note that these values are the standard inputs in the latent variable analysis (see [24]). We also took other values for these inputs and found that the resulting estimates are scarcely affected.

We implement the proposed algorithm given in Section 3 to conduct Bayesian analysis. Let ${\mathbf{Y}}_{\mathit{obs}}$ be the collection of observed data and ${\mathbf{Y}}_{\mathit{mis}}$ be the set of missing data. Due to the missing data, we need to draw ${\mathbf{Y}}_{\mathit{mis}}$ from $p\left({\mathbf{Y}}_{\mathit{mis}}|\mathbf{\Omega },\mathbf{Z},\mathit{\theta },{\mathbf{Y}}_{\mathit{obs}}\right)$ in MCMC sampling. This can be implemented easily since conditioning on $\mathbf{\Omega }$, $\mathbf{Z}$, and $\mathit{\theta }$, $p\left({\mathbf{Y}}_{\mathit{mis}}|\mathbf{\Omega },\mathbf{Z},\mathit{\theta },{\mathbf{Y}}_{\mathit{obs}}\right)$, independent of ${\mathbf{Y}}_{\mathit{obs}}$, has the normal distribution. Hence, drawing ${\mathbf{Y}}_{\mathit{mis}}$ is rather straightforward and fast. To obtain some idea about the number of the Gibbs sampler iterations in getting convergence, we conducted a few test runs as a pilot study and found that in all these runs, the Gibbs sampler converged in about 1000–2000 iterations, where the EPSR values [27] are less than 1.2. So, for all cases under consideration, we collect 3000 random observations after initial 2000 iterations being removed for posterior analysis.

We calculate the values of $L_{0.5}$ under each fitting. For computation, we use simulation-based method by drawing predictive values ${\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}$ from $p\left({\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}|{\mathbf{Y}}_{\mathit{obs}}\right)$, where ${\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}$ is the hypothetical replication of ${\mathbf{Y}}_{\mathit{obs}}$. Note that $p\left({\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}|{\mathbf{Y}}_{\mathit{obs}}\right)=\int p\left({\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}|\mathbf{\Omega },\mathbf{Z},\mathit{\theta }\right)p\left(\mathbf{\Omega },\mathbf{Z},\mathit{\theta }|{\mathbf{Y}}_{\mathit{obs}}\right)d\mathbf{\Omega }d\mathbf{Z}d\mathit{\theta }$. Hence, drawing ${\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}$ is rather easy when $\mathbf{\Omega },\mathbf{Z}$, and $\mathit{\theta }$ are available. Given that we have $M$ simulations from the posterior of $\mathbf{\Omega },\mathbf{Z},\mathit{\theta }$ via MCMC sampling discussed before, we just draw one ${\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}$ from $p\left({\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}|\mathbf{\Omega },\mathbf{Z},\mathit{\theta }\right)$ for each $\mathbf{\Omega },\mathbf{Z}$, and $\mathit{\theta }$ and obtain $M$ simulations in the end for ${\mathbf{Y}}_{\mathit{obs}}^{\mathit{rep}}$. Based on these simulated observations, $L_{\nu }$ measures can be estimated consistently via sample means. We draw 3000 observations after convergence of MCMC algorithm for calculating $L_{0.5}$ and the results are reported in Table 1.

Model	$L_{0.5}$	Model	$L_{0.5}$
$S=1$	2322.447	$S=6$	590.448
$S=2$	2107.514	$S=7$	572.172
$S=3$	1030.264	$S=8$	597.843
$S=4$	941.230	$S=9$	932.763
$S=5$	839.726	$S=10$	1030.264

Table 1.

Summary of $L_{0.5}$ under competing models in the analysis of cocaine use data.

Examination of Table 1 indicates that the proposed model with six to eight latent states seems to give better fits to the data. Furthermore, we calculate the posterior predictive density estimates of $y_{\mathit{tj}}\left(t=1,\cdots ,4,j=1,\cdots ,3\right)$ under one state and seven states, respectively (see Figure 1). It can be seen clearly that our proposed method is successful in capturing the skewness and modes of data while factor analysis model fails. For the computation details, we choose 60–100 equally spaced grids in the interval $\left[\text{min}\left\{y_{\mathit{obs},\mathit{itj}}\right\}-1.0,\text{max}\left\{y_{\mathit{obs},\mathit{itj}}\right\}+1.0\right]$ and collect 3000 simulated observations from the Gibbs sampler at each point after removing initial 2000 iteration as burn-ins.

Table 2 presents the summary of Bayesian estimates of unknown parameters and their standard errors using the formula given in (11) with $S=1$ (denoted by FA) and $S=7$ (denoted by HMFA). For comparison, maximum likelihood estimates of unknown parameters with their standard deviations under HMFA are also presented in Table 2. The maximum likelihood analysis is conducted via MCECM algorithm [36] and the standard error estimates are calculated via Louis formula [37].

	FA		ML		HMFA
Para.	Est.	SD	Est.	SD	Est.	SD
${\mathrm{\Lambda }}_{32}$	0.001	0.014	0.196	0.029	0.752	0.045
${\mathrm{\Psi }}_{\epsilon 1}$	1.443	0.315	0.559	0.297	0.432	0.049
${\mathrm{\Psi }}_{\epsilon 2}$	0.439	0.056	0.204	0.039	0.339	0.034
${\mathrm{\Psi }}_{\epsilon 3}$	0.305	0.030	0.008	NAN	0.025	0.001
${\mathrm{\Phi }}_{11}$	0.770	0.315	0.510	0.132	0.346	0.049
${\mathrm{\Phi }}_{12}$	−0.018	0.018	−0.053	0.041	−0.182	0.052
${\mathrm{\Phi }}_{22}$	1.007	0.080	0.312	0.053	0.219	0.033

Table 2.

Summary statistics for Bayesian and ML estimates in the cocaine use data analysis.

Based on Table 2, we can find the following facts: First, three estimates of ${\mathrm{\Lambda }}_{32}$ give the positive effects of latent factor $\xi$ on the ‘MFT’. This is not surprising since $\xi$ is related to the treatment level of a patient received. But there are obvious differences in magnitudes among the three methods. For FA and HMFA, the former gives ${\widehat{\mathrm{\Lambda }}}_{32}=0.001$ associated with standard deviation 0.014, while the latter gives ${\widehat{\mathrm{\Lambda }}}_{32}=0.752$ with standard deviation 0.045. This reflects that the heterogeneity of data affects the estimates ${\widehat{\mathrm{\Lambda }}}_{32}$ seriously. Compared to the previous two methods, ML method produces that ${\widehat{\mathrm{\Lambda }}}_{32}=0.196$ with SD = 0.029, which are in between them. Second, the estimates of variance parameters ${\mathrm{\Psi }}_{\mathrm{ϵ}j}$ under $S=1$ are larger than those under $S=7$. This indicates that factor analysis model accommodates heavy tails of data at the expense of variance inflation. Further investigations on the estimates of ${\mathrm{\Phi }}_{\mathit{jj}}$ under FA and HMFA also reveal the same phenomenon as that of ${\mathrm{\Psi }}_{\mathrm{ϵ}j}$. However, we observe that the ML estimate of ${\mathrm{\Psi }}_{\mathrm{ϵ}3}$, the unique variance corresponding to the third item, is equal to 0.008 with SD = NAN, an illogical number, which is very close to an improper Heywood case. As pointed out by Lee [18], Heywood cases in the ML estimation can be avoided by imposing an inequality constraint on ${\mathrm{\Psi }}_{\mathrm{ϵ}3}$ with a penalty function. In the Bayesian approach, the conjugate prior distribution of ${\mathrm{\Psi }}_{\mathrm{ϵ}3}^{-1}$ specified ${\mathrm{\Psi }}_{\mathrm{ϵ}3}$ in a region of positive values and hence has a similar effect as adding a penalty function. Hence, no Heywood cases are found in the Bayesian solution because of the penalty function induced by the prior distribution on ${\mathrm{\Psi }}_{\mathrm{ϵ}3}^{-1}$. Third, three estimates give the negative correlation between $\eta$ and $\xi$, which is consistent with the fact that the improvement of treatment will decrease the intensity of cocaine use, thus leading to a decrease of cocaine use in days. ML estimates for ${\mathrm{\Phi }}_{\mathit{jk}}$ are very close to those under HMFA. However, the estimate of ${\mathrm{\Phi }}_{12}$ under $S=1$ is −0.018, which is quite different from $-0.182$ for $S=7$. Furthermore, the coefficients of correlation of $\xi$ and $\eta$ under $S=1$ and $S=7$ are −0.0204 and −0.6612, respectively. The former suggests that $\xi$ and $\eta$ are approximately independent while the latter implies stronger dependence between them.

Moreover, we computed the posterior probabilities $P\left(z_t=r|{\mathbf{Y}}_{\mathit{obs}}\right)$ for $r=1,\cdots ,7$ and $t=1,\cdots ,4$ under $S=7$ based on 10,000 simulated observations drawn from $p\left(\mathbf{Z}|{\mathbf{Y}}_{\mathit{obs}}\right)$ and found that the transition path corresponding to the maximum posterior probability is $7\to 1\to 1\to 1$. This implies that latent state of the patient being in is extremely serious at baseline and becomes moderate in the subsequent treatments. This also reflects a positive effect of intervention on the patient’s latent state. Note that unlike the common Viterbi algorithm in exploring the optimal transition path of states in ML analysis, calculating posterior probability $P\left(z_t=r|{\mathbf{Y}}_{\mathit{obs}}\right)$ within Bayesian framework is a by-product of the estimation procedure. This voids the complex computation of marginal likelihood of the observed data and hence is very fast.

5. Discussion

This chapter reviews Bayesian inferences within a general framework and proposes a Bayesian procedure for analyzing hidden Markov factor analysis model under multivariate longitudinal setting. Compared to ML method, the pragmatic advantage of Bayesian framework is its flexibility and generality for coping with very complex problems. When good prior information can be available, results obtained from Bayesian method are more reliable and accurate than that under ML. With increased access to computation advances in simulation-based approaches, in particular the MCMC methodology, Bayesian inferences provide enormous scope for realistic statistical modeling.

Although we concentrate our attention on applications of the hidden Markov factor analysis model, the methodology developed in this chapter can be extended to the case where the LVM is nonlinear. Another possible extension is to consider a dynamic LVM, wherein model parameters vary over time. These extensions will raise theoretical and computational challenges and certainly require further investigation.

Acknowledgments

The authors are grateful to the editor’s valuable suggestions and comments which have greatly improved the manuscript. Xia’s work was fully supported by grant from National Nature Science Foundation of China (No. 11471161) and Tang’s work was supported by National Science Fund for Distinguished Young Scholars of China (No. 11225103).

1. (a) $p\left(\mathbf{Z}|\mathit{\theta },\mathit{\delta },\mathbf{Q},\mathbf{\Omega },\mathbf{Y}\right)$

Let ${\mathit{\omega }}_i$ denote the sequence of latent factors across $T$ occasions for individual $i$. To draw state variables $\mathbf{Z}$ from $p\left(\mathbf{Z}|\mathit{\delta },\mathbf{Q},\mathit{\theta },\mathbf{\Omega },\mathbf{Y}\right)$, we first notice that

$$p\left(\mathbf{Z}|\mathit{\theta },\mathbf{\Omega },\mathbf{Y}\right)=\prod _{i=1}^{N}p\left({\mathbf{z}}_i|{\mathit{\omega }}_i,\mathit{\delta },\mathbf{Q},\mathit{\theta },{\mathbf{y}}_i\right)$$E16

Hence, drawing $\mathbf{Z}$ can be accomplished via single-component method by drawing ${\mathbf{z}}_i$ independently from $p\left({\mathbf{z}}_i|{\mathit{\omega }}_i,\mathit{\delta },\mathbf{Q},\mathit{\theta },{\mathbf{y}}_i\right)$. Furthermore, notice that the sequences $\left\{{\mathbf{y}}_i,{\mathit{\omega }}_i,{\mathbf{z}}_i\right\}$ are still the one-order Markov sequences. Hence, we can simulate ${\mathbf{z}}_i$ through a well-known forward filtering-backward sampling algorithm (see, for example, [38]). For notation clarity, we suppress $\mathit{\theta }$, $\mathit{\delta }$, and $\mathbf{Q}$ in the following derivations.

Forward filtering-backward sampling (FFBS) consists of first forward filtering (FF) and then backward sampling (BS). The forward filtering step recursively updates

$${\alpha }_{i,t\mid t}=p\left(z_{\mathit{it}}|{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right),t=1,\dots ,T.$$E17

Here ${\mathbf{y}}_{i,1:t}$ represents the set of observations of subject $i$ up to time $t$ and so are ${\mathit{\omega }}_{i,1:t}$ and ${\mathbf{z}}_{i,1:t}$. The backward sampling is to draw ${\mathbf{z}}_i$ from the joint distribution of the states given the data using

$$p\left({\mathbf{z}}_{i,1:T}|{\mathit{\omega }}_{i,1:T},{\mathbf{y}}_{i,1:T}\right)=p\left(z_{\mathit{iT}}|{\mathit{\omega }}_{i,1:T},{\mathbf{y}}_{i,1:T}\right) \dots p\left(z_{i1}|{\mathbf{z}}_{i,2:T},{\mathit{\omega }}_{i,1:T},{\mathbf{y}}_{i,1:T}\right).$$E18

That is, we first draw the last state given all the data and then work backwards in time drawing each state conditional on all the subsequent ones.

To implement forward filtering, let

$${\alpha }_{\mathit{it}}\left(r\right)=\mathbb{P}\left({\mathbf{y}}_{i,1:t},{\mathit{\omega }}_{i,1:t},z_{\mathit{it}}=r\right),t=1,\cdots ,T$$E19

Obviously, ${\alpha }_{i1}\left(r\right)={\delta }_{r}p\left({\mathbf{y}}_{i1},{\mathit{\omega }}_{i1}|z_{i1}=r\right)$. Moreover, it can be shown that

$${\alpha }_{\mathit{it}}\left(r\right)=\left(\sum _{s=1}^S{\alpha }_{\mathit{it}-1}\left(s\right)q_{\mathit{sr}}\right)p\left({\mathbf{y}}_{\mathit{it}},{\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}}=r\right),t=2,\cdots ,T$$E20

The outputs ${\left\{{\alpha }_{\mathit{it}}\right\}}_{t=1}^T$ from recursive Eq. (20) can be used to calculate the posterior probability

$${\alpha }_{i,t\mid t}\left(r\right)=\mathbb{P}\left(z_{\mathit{it}}=r|{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)=\frac{{\alpha }_{\mathit{it}}\left(r\right)}{\sum _{s=1}^S{\alpha }_{\mathit{it}}\left(s\right)}$$E21

which leads to the forward filtering (FF) iteration.

The backward sampling step depends on the observation that

$$\begin{array}{ll}& p\left(z_{\mathit{it}}|{\mathbf{z}}_{i,t+1:T},{\mathit{\omega }}_{i,1:T},{\mathbf{y}}_{i,1:T}\right)\propto p\left(z_{\mathit{it}},{\mathbf{z}}_{i,t+1:T},{\mathit{\omega }}_{i,1:T},{\mathbf{y}}_{i,1:T}\right)\\ =& p\left(z_{\mathit{it}},{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)p\left({\mathbf{z}}_{i,t+1:T},{\mathit{\omega }}_{i,t+1:T},{\mathbf{y}}_{i,t+1:T}|z_{\mathit{it}},{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)\\ =& p\left(z_{\mathit{it}},{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)p\left({\mathbf{z}}_{i,t+1:T}|z_{\mathit{it}},{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)p\left({\mathit{\omega }}_{i,t+1:T},{\mathbf{y}}_{i,t+1:T}|{\mathbf{z}}_{i,t+1:T},z_{\mathit{it}},{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)\\ =& p\left(z_{\mathit{it}},{\mathit{\omega }}_{i,1:t},{\mathbf{y}}_{i,1:t}\right)p\left({\mathbf{z}}_{i,t+1:T}|z_{\mathit{it}},\right)p\left({\mathit{\omega }}_{i,t+1:T},{\mathbf{y}}_{i,t+1:T}|{\mathbf{z}}_{i,t+1:T}\right)\end{array}$$E22

The last equation holds since given $z_{\mathit{it}}$, $\left\{{\mathbf{y}}_{i,t:T},{\mathit{\omega }}_{i,t:T},{\mathbf{z}}_{i,t+1:T}\right\}$ does not depend on the previous values due to the Markov Chain characteristics of $\left\{{\mathbf{y}}_{\mathit{it}},{\mathit{\omega }}_{\mathit{it}},{\mathbf{z}}_{\mathit{it}}\right\}$. This leads to

$$\mathbb{P}\left(z_{\mathit{it}}=r|{\mathbf{z}}_{i,t+1:T},{\mathbf{y}}_{i,1:T},{\mathit{\omega }}_{i,1:T},\right)=\frac{{\alpha }_{i,t\mid t}\left(r\right)q_{{\mathit{rz}}_{i,t+1}}}{\sum _{s=1}^S{\alpha }_{i,t\mid t}\left(s\right)q_{{\mathit{sz}}_{i,t+1}}}t=T-1,\cdots ,1.$$E23

Hence, FFBS algorithm for drawing ${\mathbf{z}}_i$ is implemented by

Algorithm:

1. running the recursion ${\alpha }_{\mathit{it}}$ and stored the conditional probabilities ${\alpha }_{i,t\mid t}$ for $t=1,\dots ,T$;

2. sampling $z_{\mathit{iT}}$ from the filtered conditional probability ${\alpha }_{i,T\mid T}$;

3. for $t=T-1,\cdots ,1$, sampling $z_{\mathit{it}}$ from the conditional probability

$$\mathbb{P}\left(z_{\mathit{it}}=r|{\mathit{\omega }}_{i,1:T},{\mathbf{y}}_{i,1:T},{\mathbf{z}}_{i,t+1:T}\right).$$E24

1. (b) $p\left(\mathbf{\Omega }|\mathbf{Z},\mathit{\theta },\mathbf{Y}\right)$

To draw $\mathbf{\Omega }$, we first note that

$$p\left(\mathbf{\Omega }|\mathbf{Z},\mathit{\theta },\mathbf{Y}\right)=\prod _{i=1}^N\prod _{t=1}^{T}p\left({\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}},\mathit{\theta },{\mathbf{y}}_{\mathit{it}}\right)$$E25

in which

$$\begin{array}{l}p\left({\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}},\mathit{\theta },{\mathbf{y}}_{\mathit{it}}\right)\propto \text{exp}\left\{-\frac{1}{2}{\left({\mathbf{y}}_{\mathit{it}}-{\mathit{\mu }}_r-{\mathbf{\Lambda }}_r{\mathit{\omega }}_{\mathit{it}}\right)}^{⊺}{\mathbf{\Psi }}_{\mathrm{ϵ}r}^{-1}\left({\mathbf{y}}_{\mathit{it}}-{\mathit{\mu }}_{\mathit{kr}}-{\mathbf{\Lambda }}_r{\mathit{\omega }}_{\mathit{it}}\right)\right.\\ \left.-\frac{1}{2}{\mathit{\omega }}_{\mathit{it}}^{⊺}{\mathbf{\Phi }}_r^{-1}{\mathit{\omega }}_{\mathit{it}}\right\}\end{array}$$E26

with $r=z_{\mathit{it}}$. Hence, similar to that in drawing $\mathbf{Z}$, updating $\mathbf{\Omega }$ can be achieved by drawing ${\mathit{\omega }}_{\mathit{it}}$ independently from $p\left({\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}},\mathit{\theta },{\mathbf{y}}_{\mathit{it}}\right)$ for $i=1,\cdots ,N$ and $t=1,\cdots ,T$. It can be shown that

$$p\left({\mathit{\omega }}_{\mathit{it}}|z_{\mathit{it}}=r,\mathit{\theta },{\mathbf{y}}_{\mathit{it}}\right){\displaystyle \stackrel{D}{=}}{\mathcal{N}}_m\left({{\displaystyle \widehat{\mathbf{m}}}}_{\mathit{it}},{\widehat{\mathbf{\Sigma }}}_r\right).$$E27

in which

$${\widehat{\mathbf{m}}}_{\mathit{it}}={\widehat{\mathbf{\Sigma }}}_r{\mathbf{\Lambda }}_r^{⊺}{\mathbf{\Psi }}_{\mathrm{ϵ}r}^{-1}\left({\mathbf{y}}_{\mathit{it}}-{\mathit{\mu }}_r\right),{\widehat{\mathbf{\Sigma }}}_r={\left({\mathbf{\Lambda }}_r^{⊺}{\mathbf{\Psi }}_{\mathrm{ϵ}r}^{-1}{\mathbf{\Lambda }}_r+{\mathbf{\Phi }}_r^{-1}\right)}^{-1}.$$E28

1. (c) $p\left(\mathit{\mu },\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}|\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)$

To draw $\left\{\mathit{\mu },\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}\right\}$, we can first draw $\mathit{\mu }$ from $p\left(\mathit{\mu }|\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}},\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)$ and then draw $\left\{\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}\right\}$ from $p\left(\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}|\mathit{\mu },\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)$. For this end, let ${\widehat{n}}^{\left(r\right)}=\#\left\{z_{\mathit{it}}=r\right\}$ be the size of cluster $z_{\mathit{it}}$, and let $w_{\mathit{itr}}=I\left\{z_{\mathit{it}}=r\right\}$. Denote

$$\begin{array}{l}{{\displaystyle \overline{\mathbf{Y}}}}^{\left(r\right)}=\sum _{i=1}^N\sum _{t=1}^T{w}_{\mathit{itr}}{\mathbf{y}}_{\mathit{it}}/{{\displaystyle \widehat{n}}}^{\left(r\right)},{{\displaystyle \overline{\mathbf{\Omega }}}}^{\left(r\right)}=\sum _{i=1}^N\sum _{t=1}^T{w}_{\mathit{itr}}{\mathit{\omega }}_{\mathit{it}}/{{\displaystyle \widehat{n}}}^{\left(r\right)},\hfill \\ {\mathbf{S}}_{\mathit{yy}}^{\left(r\right)}=\sum _{i=1}^N\sum _{t=1}^T{w}_{\mathit{itr}}{\mathbf{y}}_{\mathit{it}}{\mathbf{y}}_{\mathit{it}}^{⊺}/{{\displaystyle \widehat{n}}}^{\left(r\right)},{\mathbf{S}}_{\mathit{\omega y}}^{\left(r\right)}=\sum _{i=1}^N\sum _{t=1}^T{w}_{\mathit{itr}}{\mathit{\omega }}_{\mathit{it}}{\mathbf{y}}_{\mathit{it}}^{⊺}/{{\displaystyle \widehat{n}}}^{\left(r\right)},\hfill \end{array}$$E29

be the sample means and covariance matrices of $\mathbf{Y}$ and $\mathbf{\Omega }$ within the $r$th cluster, respectively. By some algebra calculations, it can be shown that

$$\begin{array}{l}p\left(\mathit{\mu }|\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}},\mathbf{K},\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)=\prod _{r=1}^{S}p\left({\mathit{\mu }}_r|{\mathbf{\Lambda }}_r,{\mathbf{\Psi }}_{\mathrm{ϵ}r},\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right),\text{and}\\ p\left(\mathbf{\Lambda },{\mathbf{\Psi }}_{\mathrm{ϵ}}|\mathit{\mu },\mathbf{K},\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)=\prod _{r=1}^{S}p\left({\mathbf{\Lambda }}_r,{\mathbf{\Psi }}_{\mathrm{ϵ}r}|{\mathit{\mu }}_r,\mathbf{\Omega },\mathbf{Z},\mathbf{Y}\right),\end{array}$$E30

where

$$\begin{array}{ll}p\left({\mathit{\mu }}_r|{\mathbf{\Lambda }}_r,{\mathbf{\Psi }}_{\mathrm{ϵ}r},\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)& ={\mathcal{N}}_p\left({\widehat{\mathbf{a}}}_{\mathit{\mu r}},{\widehat{\mathbf{\Sigma }}}_{\mathit{\mu r}}\right),\\ p\left({\mathbf{\Lambda }}_r,{\mathbf{\Psi }}_{\mathrm{ϵ}r}|{\mathit{\mu }}_r,\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)& =\prod _{j=1}^{p}p\left({\Psi }_{\mathrm{ϵ}\mathit{rj}}|{\mathit{\mu }}_r,\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)p\left({\mathbf{\Lambda }}_{\mathit{rj}}|{\Psi }_{\mathrm{ϵ}\mathit{rj}}|{\mathit{\mu }}_r,\mathbf{Z},\mathbf{\Omega },\mathbf{Y}\right)\\ & {\displaystyle \stackrel{D}{=}}\prod _{j=1}^p\mathcal{G}{a}^{-1}\left({\widehat{\alpha }}_{\mathrm{ϵ}\mathit{rj}},{\widehat{\beta }}_{\mathrm{ϵ}\mathit{rj}}\right){\mathcal{N}}_m\left({\widehat{\mathbf{\Lambda }}}_{\mathit{rj}},{\Psi }_{\mathrm{ϵ}\mathit{rj}}{\widehat{\mathbf{H}}}_{\mathit{rj}}\right),\end{array}$$E31

with

$$\begin{array}{ll}{{\displaystyle \widehat{\mathbf{a}}}}_{\mathit{\mu r}}\hfill & ={{\displaystyle \widehat{\mathbf{\Sigma }}}}_{\mathit{\mu r}}\left[{\mathbf{\Sigma }}_{0r}^{-1}{\mathit{\mu }}_{0r}+{{\displaystyle \widehat{n}}}^{\left(r\right)}{\mathbf{\Psi }}_{\mathrm{ϵ}r}^{-1}\left({{\displaystyle \overline{\mathbf{Y}}}}^{\left(r\right)}-{\mathbf{\Lambda }}_r{{\displaystyle \overline{\mathbf{\Omega }}}}^{\left(r\right)}\right)\right],{{\displaystyle \widehat{\mathbf{\Sigma }}}}_{\mathit{\mu r}}={\left({\mathbf{\Sigma }}_{0r}^{-1}+{{\displaystyle \widehat{n}}}^{\left(r\right)}{\mathbf{\Psi }}_{\mathrm{ϵ}r}^{-1}\right)}^{-1},\\ {{\displaystyle \widehat{\mathbf{\Lambda }}}}_{\mathit{rj}}\hfill & ={{\displaystyle \widehat{\mathbf{H}}}}_{\mathit{rj}}\left[{\mathbf{H}}_{\mathrm{ϵ}0\mathit{rj}}^{-1}{\mathbf{\Lambda }}_{0\mathit{rj}}+{{\displaystyle \widehat{n}}}^{\left(r\right)}\left({\mathbf{S}}_{\mathit{\omega y}\left(j\right)}^{\left(r\right)}-{\mu }_{\mathit{rj}}{{\displaystyle \overline{\mathbf{\Omega }}}}^{\left(r\right)}\right)\right],{{\displaystyle \widehat{\mathbf{H}}}}_{\mathit{rj}}^{-1}={\mathbf{H}}_{0\mathit{rj}}^{-1}+{{\displaystyle \widehat{n}}}^{\left(r\right)}{\mathbf{S}}_{\mathit{\omega \omega }}^{\left(r\right)},\\ {{\displaystyle \widehat{\alpha }}}_{\mathrm{ϵ}\mathit{rj}}\hfill & ={\alpha }_{\mathrm{ϵ}0\mathit{rj}}+{{\displaystyle \widehat{n}}}^{\left(r\right)}/2,\\ {{\displaystyle \widehat{\beta }}}_{\mathrm{ϵ}\mathit{rj}}\hfill & ={\beta }_{\mathrm{ϵ}0\mathit{rj}}+\left\{{\mathbf{\Lambda }}_{0\mathit{rj}}^{⊺}{H}_{0\mathit{rj}}^{-1}{\mathbf{\Lambda }}_{0\mathit{rj}}+{{\displaystyle \widehat{n}}}^{\left(r\right)}\left(S_{\mathit{yy}\left(j,j\right)}^{\left(r\right)}-2{\mu }_{\mathit{rj}}{{\displaystyle \overline{y}}}_j^{\left(r\right)}+{\mu }_{\mathit{rj}}^2\right)-{{\displaystyle \widehat{\mathbf{\Lambda }}}}_{\mathit{rj}}^{⊺}{{\displaystyle \widehat{\mathbf{H}}}}_{\mathit{rj}}^{-1}{{\displaystyle \widehat{\mathbf{\Lambda }}}}_{\mathit{rj}}\right\}/2,\end{array}$$E32

in which ${\overline{y}}_{\left(j\right)}^{\left(r\right)}$ is the $j$th element in ${\overline{\mathbf{Y}}}^{\left(r\right)}$, $S_{\mathit{yy}\left(j,j\right)}^{\left(r\right)}$ is the $j$th main diagonal element of ${\mathbf{S}}_{\mathit{yy}}^{\left(r\right)}$, and ${\mathbf{S}}_{\mathit{\omega y}\left(j\right)}^{\left(r\right)}$ is the $j$th column vector of ${\mathbf{S}}_{\mathit{\omega y}}^{\left(r\right)}$.

1. (d) $p\left(\mathbf{\Phi }|\mathbf{\Omega },\mathbf{Z}\right)$

From the prior distribution of ${\mathbf{\Phi }}_r^{-1}$ and the distribution of $\mathbf{\Omega }$, it can be shown that

$$p\left({\mathbf{\Phi }}_r|\mathbf{\Omega },\mathbf{Z}\right)\propto {\left|{\mathbf{\Phi }}_r\right|}^{\left({\widehat{n}}^{\left(r\right)}+{\rho }_{0r}+m+1\right)/2}\text{exp}\left\{-\frac{1}{2}\mathrm{tr}{\mathbf{\Phi }}^{-1}({\widehat{n}}^{\left(r\right)}{\mathbf{S}}_{\mathit{\omega \omega }}^{\left(r\right)}+{\mathbf{R}}_0^{-1}\right\}$$E33

where ${\widehat{n}}^{\left(r\right)}$ and ${\mathbf{S}}_{\mathit{\omega \omega }}$ are given in (c). Hence, $p\left({\mathbf{\Phi }}_r|\mathbf{\Omega },\mathbf{Z}\right)$ is the $m$-dimensional inverse Wishart distribution ${\mathcal{W}}_m^{-1}\left({\widehat{n}}^{\left(r\right)}+{\rho }_{0r},{\widehat{n}}^{\left(r\right)}{\mathbf{S}}_{\mathit{\omega \omega }}^{\left(r\right)}+{\mathbf{R}}_0^{-1}\right)$. It can be shown from exactly the same reasoning as before that drawing $\mathbf{\Phi }$ can be achieved by drawing ${\mathbf{\Phi }}_r$ from $p\left({\mathbf{\Phi }}_r|\mathbf{\Omega },\mathbf{Z}\right)$ independently.

1. (e) $p\left(\mathit{\delta }|\mathbf{Z}\right)$ and (f) $p\left(\mathbf{Q}|\mathbf{Z}\right)$

It can be verified directly that

$$\begin{array}{r}p\left(\mathit{\delta }|\mathbf{Z}\right)=p\left({\mathit{\delta }}_k|\mathbf{Z}\right)\text{and}\\ p\left(\mathit{\delta }|\mathbf{Z}\right){\displaystyle \stackrel{D}{=}}\mathcal{D}{\mathit{ir}}_S\left({\gamma }_0+{\widehat{n}}_{11},\dots ,{\gamma }_0+{\widehat{n}}_{1S}\right)\end{array}$$E34

in which ${\widehat{n}}_{1r}=\sum _{i=1}^{N}I\left\{z_{i1}=r\right\}$. Similarly, it can be shown that

$$\begin{array}{l}p\left(\mathbf{Q}|\mathbf{Z}\right)=\prod _{r=1}^{S}p\left({\mathbf{Q}}_r|\mathbf{Z}\right),\\ p\left({\mathbf{Q}}_r|\mathbf{Z}\right){\displaystyle \stackrel{D}{=}}\prod _{r=1}^S\mathcal{D}{\mathit{ir}}_S\left({\nu }_0+{\widehat{n}}_{r1},\dots ,{\nu }_0+{\widehat{n}}_{\mathit{rS}}\right).\end{array}$$E35

in which ${\widehat{n}}_{\mathit{rs}}=\sum _{i=1}^N\sum _{t=2}^{T}I\left\{z_{\mathit{it}-1}=r,z_{\mathit{it}}=s\right\}$.