Bayesian quantile nonhomogeneous hidden Markov models

Abstract

Hidden Markov models are useful in simultaneously analyzing a longitudinal observation process and its dynamic transition. Existing hidden Markov models focus on mean regression for the longitudinal response. However, the tails of the response distribution are as important as the center in many substantive studies. We propose a quantile hidden Markov model to provide a systematic method to examine the entire conditional distribution of the response given the hidden state and potential covariates. Instead of considering homogeneous hidden Markov models, which assume that the probabilities of between-state transitions are independent of subject- and time-specific characteristics, we allow the transition probabilities to depend on exogenous covariates, thereby yielding nonhomogeneous Markov chains and making the proposed model more flexible than its homogeneous counterpart. We develop a Bayesian approach coupled with efficient Markov chain Monte Carlo methods for statistical inference. Simulations are conducted to assess the empirical performance of the proposed method. The proposed methodology is applied to a cocaine use study to provide new insights into the prevention of cocaine use.

Keywords

Asymmetric Laplace distribution Bayesian analysis hidden Markov model nonhomogeneous transition quantile regression

1 Introduction

Heterogeneous longitudinal data are frequently encountered in medical, behavioral, socioeconomic, environmental, and psychological sciences. Hidden Markov models (HMMs), which consist of a transition model to describe the dynamic transition of hidden states and a conditional regression model to examine state-specific covariate effects on the response of interest, provide a useful model framework to accommodate heterogeneous and longitudinal features.^1–3 HMMs and their variants have elicited extensive attention in various disciplines due to their ability to simultaneously reveal the longitudinal association structure and dynamic heterogeneity of the observed process.^4–12 An R package for HMMs was also developed by Visser and Speekenbrink.¹³

Despite the rapid development and wide application of HMMs, existing literature has focused on mean regression-based HMMs,^1–15 in which the conditional model is a mean regression for examining the effects of potential covariates on the average of the response variable. One problem of mean regression-based models is that they are sensitive to the violation of distributional assumptions, such as nonnormality of the response variable and presence of outliers. Moreover, mean regression-based models only assess covariate effects on the central tendency of the response distribution. However, in many circumstances, covariate effects on the tails of the response distribution are of even more interest than those on the center. A highly comprehensive analysis should perform quantile regression (QR) in the context of HMMs so that the dynamic heterogeneity in the effects of potential covariates on the entire distribution of the response variable can be investigated.

QR^16,17 is a valuable alternative to its mean counterpart. It reveals the effects of covariates on the entire conditional distributions of the response variable instead of the effects on its average only. Moreover, QR makes minimal assumptions on the error distribution and thus enables the accommodation of nonnormal data that are frequently encountered in substantive research. QR has received considerable attention in statistical applications and has recently been introduced into HMMs due to its appealing features. For example, Farcomeni¹⁸ used QR to analyze longitudinal data with latent Markov subject-specific parameters and adopted the expectation–maximization algorithm to perform estimation. Marino and coworkers^19–21 considered informative missingness and latent dropout-based transitions in the context of quantile HMMs and proposed maximum likelihood approaches for statistical inference. However, the preceding analyses were conducted within a frequentist framework and based on a restrictive assumption that the between-state transitions are homogeneous, which implies that the probability of transitioning from one state to another is independent of subject-specific characteristics and observation times. This homogeneous assumption is unrealistic because the dynamic transition between hidden states is actually influenced by certain covariates in many circumstances. To our knowledge, no existing study has investigated quantile HMMs with nonhomogeneous between-state transitions.

We aim to fill this gap and consider a novel quantile HMM that allows nonhomogeneous transitions of hidden states. A continuation-logit model is proposed to examine how subject- and/or time-specific covariates influence the transition probabilities. We develop a Bayesian approach coupled with Markov chain Monte Carlo (MCMC) sampling techniques to conduct statistical inference because of this approach’s potential in managing highly complex data/model structures and capability to incorporate additional model inputs that stem from the prior distributions of model parameters. The majority of extant studies on Bayesian QR concentrated on cross-sectional data analyses. Yu and Moyeed²² proposed the asymmetric Laplace distribution (ALD) approach to facilitate the Bayesian inference of QR. Kottas and Gelfand²³ conducted Bayesian QR based on Dirichlet process priors. Dunson et al.²⁴ developed a substitution likelihood-based approach for QR with latent variables. Readers can also refer to other studies in relevant fields^25–32 and the R package bayesQR.³³ The only exception is the study of Koutsourelis,³⁴ who considered homogeneous quantile HMMs in the Bayesian framework. However, these available methods cannot jointly accommodate all of the abovementioned features. This study is the first to provide a Bayesian method to analyze nonhomogeneous quantile HMMs. We utilize the ALD approach proposed by Yu and Moyeed²² to facilitate Bayesian inference for the proposed quantile HMMs not only because the approach provides a natural and effective means to model Bayesian QR but also because it enables highly efficient MCMC algorithms for posterior sampling.³⁵ Moreover, Sriram et al.³⁶ provided theoretical guarantees for this approach by establishing the posterior consistency of Bayesian QR on the basis of the misspecified asymmetric Laplace density in the case of QR on independent responses.

The proposed model is motivated by the longitudinal analysis of cocaine use as described in the Longitudinal Analysis of Cocaine Use section, in which the response of our primary interest, $Y_{i, t}$ , was measured as days of cocaine use per month. The main goal of this study is to investigate the potential influential factors of $Y_{i, t}$ , such as treatment received and patients’ mental health. Figure 1 presents the histograms of $Y_{i, t}$ at different time points and indicates that the distribution of $Y_{i, t}$ is highly skewed. Moreover, the magnitudes of covariate effects may vary across different states, leading to state- and quantile-specific covariate effects. Our proposed quantile HMM accommodates the nonnormality of $Y_{i, t}$ and provides a more complete picture of the associations between cocaine use and their influential factors than mean regression does. The proposed model also accommodates the heterogeneity in the transition probabilities, which cannot be identified by existing homogeneous HMMs.

Figure 1.

Histograms of the days of cocaine use per month (y_it) at different time points.

The rest of this article is organized as follows. The next section defines the proposed nonhomogeneous quantile HMM. This is followed by the Bayesian estimation section that describes the Bayesian approach and the implementation of the MCMC algorithm. The model identifiability issue and the associated solution are also discussed. Then, simulation studies to evaluate the empirical performance of the proposed methodology are presented in the Simulation study section. The penultimate section provides an application of the proposed method to a cocaine use study, and the final section concludes the paper.

2 Model settings

Let ${D_{i, t} = (Y_{i, t}, X_{i, t})}_{t = 1}^{T}$ be the observation process for the ith subject, $i = 1, \dots, N$ , where $X_{i, t}$ is the covariate vector in the model, and $Y_{i, t}$ is the response variable. ${Z_{i, t}}_{t = 1}^{T}$ is the associated hidden state process with hidden state space ${1, \dots, S}$ . In a given hidden state s, the observations $D_{i, t}$ are assumed to be independent for all i and t. In real applications, the number of hidden states S is usually unknown. We assume that S is invariant to quantile levels so that we can decide on S before the quantile analysis. A commonly used method for determining S is the utilization of a modified deviance information criterion.^15,37 In the subsequent context, we assume that S is predetermined.

We describe the conditional model given the hidden state. For a given quantile level $0 < τ < 1$ , similar to Marino et al.,²¹ we assume that the conditional density of $Y_{i, t}$ given ${X_{i, t}, Z_{i, t} = s, τ}$ is

f (Y_{i, t} | X_{i, t}, Z_{i, t} = s, τ) = \frac{τ (1 - τ)}{σ_{s}} \exp {- ρ_{τ} (Y_{i, t} - β_{τ, s}^{'} X_{i, t}) / σ_{s}}

(1)

where

β_{τ, s}^{'} X_{i, t}

is the location of the distribution, with

X_{i, t} = {(1, X_{i, t, 1}, \dots, X_{i, t, q})}^{'}

and

β_{τ, s} = {(β_{τ, s, 1}, β_{τ, s, 2}, \dots, β_{τ, s, q + 1})}^{T}

being the covariate and coefficient vectors, σ_s is the scale parameter, and

ρ_{τ} (t) = t {τ - I (t < 0)}

is the quantile check loss function. Let

ε_{i, t} (τ) | {Z_{i, t} = s} = Y_{i, t} - β_{τ, s}^{'} X_{i, t}

(2)

If S = 1, then models (1) and (2) are the same as the model of Yu and Moyeed,²² which assumes an ALD for noise. In the following, we omit τ from $ε_{i, t} (τ)$ for notational simplicity.

Next, we present the model setting for the hidden state process. Following existing literature,⁷,^14,15 we assume that ${Z_{i, t}}_{t = 1}^{T}$ is Markovian of order 1 and satisfies

P (Z_{i, t} = s | Z_{i, 1}, \dots, Z_{i, t - 1} = u) = P (Z_{i, t} = s | Z_{i, t - 1} = u) ≐ p_{i, t; u, s}, t \geq 2

(3)

where

p_{i, t; u, s}

is the transition probability from state

Z_{i, t - 1} = u

at occasion t – 1 to state

Z_{i, t} = s

at occasion t for the ith individual. The initial distribution of

Z_{i, 1}

is assumed to be a multinomial distribution with probabilities

(p_{1}, p_{2}, \dots, p_{S})

such that

p_{s} \geq 0

and

\sum_{s = 1}^{S} p_{s} = 1

. The distribution of

{Z_{i, t}}_{t = 1}^{T}

is then fully determined by the transition probabilities and the distribution of the initial state.

To model the transition probabilities, we assume that the hidden states ${1, \dots, S}$ are ordered. Instead of modeling p_{i,t;u,s} directly, we consider a continuation-ratio logistic model for $η_{i, t; u, s} = P (Z_{i, t} = s | Z_{i, t} \geq s, Z_{i, t - 1} = u)$ with $s \leq S - 1$ . Specifically, we assume that

\begin{array}{l} logit (η_{i, t; u, s}) = \log \frac{P (Z_{i, t} = s | Z_{i, t - 1} = u)}{P (Z_{i, t} > s | Z_{i, t - 1} = u)} = \log \frac{p_{i, t; u, s}}{p_{i, t; u, s + 1} + \dots + p_{i, t; u, S}} = ζ_{u, s} + α^{'} b_{i, t} \\ i = 1, \dots, N, t = 2, \dots, T, u = 1, \dots, S, s = 1, \dots, S - 1 \end{array}

(4)

where

{ζ_{u, s} : u = 1, \dots, S, s = 1, \dots, S - 1}

are state-specific intercepts,

b_{i, t}

is an h-dimensional covariate vector that may influence the transition probabilities, and

α

is an h-dimensional coefficient vector that can be interpreted as a conditional log odds ratio in a logistic regression. Such parameterization aims to facilitate the interpretation of transition to a state rather than to a better one. Similar to that in the proportional odds model,

α

in (4) is assumed to be independent of u and s to maintain the order of the hidden states and avoid a tedious transition model, in which every transition elicits a set of parameters for all possible states of origination and destination. This outcome greatly reduces the complexity and enhances the interpretability of the transition model. On the basis of (4), we obtain the following:

\begin{array}{l} p_{i, t; u, s} = \frac{\exp (ζ_{u, s} + α^{'} b_{i, t})}{\prod_{k = 1}^{s} {1 + \exp (ζ_{u, k} + α^{'} b_{i, t})}}, s = 1, \dots, S - 1, \\ p_{i, t; u, S} = \frac{1}{\prod_{k = 1}^{S - 1} {1 + \exp (ζ_{u, k} + α^{'} b_{i, t})}} \end{array}

(5)

Notably, although we assume an ordering of the underlying states, we do not restrict transitions between states. The transition process is allowed to be reversible with transition probability $p_{i, t; u, s} = P (Z_{i, t} = s | Z_{i, t - 1} = u)$ , where s and u can have arbitrary values in ${1, 2, \dots, S}$ . Thus, the transition from one state to another is unnecessarily successive or proceeds only in a single direction. For any given subject i and time t, we can use (5) to calculate a transition matrix $P (i, t) = {[p_{i, t; u, s}]}_{S \times S}$ , which allows for forward and reverse transitions.

3 Bayesian estimation

3.1 Identification

Labeling the component membership is an important issue when studying mixture-type models, including HMMs. If label switching occurs, the posterior samples may jump between labeling subspaces in a balanced manner when MCMC iterations are performed, leading to meaningless estimation results. Several methods have been proposed to address the label switching problem.^38–40 Among them, imposing an ordering restriction on component means that $μ_{1} < μ_{2} < \dots < μ_{n}$ is an efficient method for preventing label switching and ensuring model identifiability.^41–43 Following this common practice, we use $β_{τ, 1, 1} < β_{τ, 2, 1} < \dots < β_{τ, S, 1}$ to obtain an identified model, where $β_{τ, s, 1}$ is the intercept of model (1) conditional on $Z_{i, t} = s$ .

3.2 Full likelihood with scale mixture normal representation

Let $D = {D_{i, t} : i = 1, \dots, N, t = 1, \dots, T}, Z = {Z_{i, t} : i = 1, \dots, N, t = 1, \dots, T},$ $P = {p_{s} : s = 1, \dots, S}, ζ = {ζ_{u, s} : u = 1, \dots, S, s = 1, \dots, S - 1},$ $β = {β_{τ, s} : s = 1, \dots, S}, σ = {σ_{s} : s = 1, \dots, S}, θ_{M} = {β, σ}$ be the vector of unknown state-specific parameters involved in the conditional likelihood. Let $θ_{H} = {P, ζ, α}$ be the vector of unknown parameters involved in the transition model, and $θ = (θ_{M}, θ_{H})$ .

On the basis of the model defined by (1) and (2), we have

f (ε_{i, t} | Z_{i, t} = s) = \frac{τ (1 - τ)}{σ_{s}} \exp {- ρ_{τ} (ε_{i, t}) / σ_{s}}

To simplify MCMC sampling, we represent the ALD of $ε_{i, t}$ by using a scale mixture normal distribution^30,35 as follows:

[ε_{i, t} | Z_{i, t} = s] = ξ_{1} v_{i, t} + ξ_{2} \sqrt{σ_{s} v_{i, t}} e_{i, t}

where

ξ_{1} = \frac{1 - 2 τ}{τ (1 - τ)}, ξ_{2}^{2} = \frac{2}{τ (1 - τ)}, v_{i, t} \sim Exp (σ_{s})

, and

e_{i, t} \sim N (0, 1)

. Thus, the conditional models (1) and (2) can be presented as

[Y_{i, t} | X_{i, t}, Z_{i, t} = s, τ] = β_{τ, s}^{'} X_{i, t} + ξ_{1} v_{i, t} + ξ_{2} \sqrt{σ_{s} v_{i, t}} e_{i, t}

(6)

Accordingly, the conditional likelihood of $D_{i, t}$ given $Z_{i, t} = s$ can be written as

\begin{array}{l} p (D_{i, t} | θ_{M}, Z_{i, t} = s) = p (D_{i, t} | θ_{M}, Z_{i, t} = s, v_{i, t}) p (v_{i, t} | θ_{M}, Z_{i, t} = s) \\ = \frac{1}{\sqrt{2 π ξ_{2}^{2} σ_{s} v_{i, t}}} \exp {- \frac{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}}} \cdot \frac{1}{σ_{s}} \exp (- v_{i, t} / σ_{s}) \\ = \frac{1}{\sqrt{2 π ξ_{2}^{2} σ_{s}^{3 / 2} v_{i, t}}} \exp {- \frac{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2} + 2 ξ_{2}^{2} v_{i, t}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}}} \end{array}

(7)

Let $V = {v_{i, t} : i = 1, \dots, N, t = 1, \dots, T}$ . Then, the complete dataset is ${D, V, Z}$ , where $D$ is the observed data, whereas $V$ and $Z$ are not directly observable. Recall that $p_{s} = P (Z_{i, 1} = s)$ and $p_{i, t; u, s} = P (Z_{i, t} = s | Z_{i, t - 1} = u)$ . By combining (5) and (7), we obtain the complete-data likelihood function as follows:

\begin{array}{l} L_{c} (D, V, Z | θ) = \prod_{s = 1}^{S} \prod_{i = 1}^{N} \prod_{t = 1}^{T} [\frac{1}{\sqrt{2 π ξ_{2}^{2} σ_{s}^{3 / 2} v_{i, t}}} \times \exp {- \frac{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2} + 2 ξ_{2}^{2} v_{i, t}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}}}]^{Z_{i, t} (s)} \\ \times \prod_{s = 1}^{S} p_{s}^{n_{s, 1}} \times \prod_{u = 1}^{S} \prod_{s = 1}^{S} \prod_{i = 1}^{N} \prod_{t = 2}^{T} p_{i, t; u, s}^{Z_{i, t - 1} (u) Z_{i, t} (s)} \end{array}

where

Z_{i, t} (s) = I (Z_{i, t} = s)

and

n_{s, 1} = \sum_{i = 1}^{N} Z_{i, 1} (s)

3.3 Prior specification and posterior sampling

To conduct a Bayesian analysis, we must initially specify prior distributions for the unknown parameters. Following existing literature,^14,15,35 we consider commonly used prior distributions for the unknown parameters as follows:

\begin{array}{l} (p_{1}, p_{2}, \dots, p_{S}) \overset{D}{=} Dirichlet (a, \dots, a), \\ ζ_{u, s} \overset{D}{=} N (ζ_{u, s}^{0}, δ_{ζ}^{2}), \\ α \overset{D}{=} N (α^{0}, Σ_{α}), \\ β_{τ, s} \overset{D}{=} N (β_{τ, s}^{0}, B_{τ, s}^{0}), \\ σ_{s} \overset{D}{=} I G (n_{s}^{0} / 2, b_{s}^{0} / 2) \end{array}

(8)

where

Dirichlet (\cdot, \dots, \cdot)

denotes the Dirichlet distribution, and

I G (\cdot, \cdot)

stands for the inverse gamma distribution.

The main task of posterior inference is to sample from $p (θ | D)$ , which is intractable due to the existence of unobserved $V$ and $Z$ . Therefore, instead of directly working on $p (θ | D)$ , we deal with $p (Z, V, θ | D) = p (Z, P, ζ, α, β, σ | D)$ . However, this joint posterior distribution is still complicated. Thus, we use the Gibbs sampler to update each unknown by sampling from its full conditional distribution through the following steps.

Step 1: Update hidden states $Z = {Z_{i, t}, i = 1, \dots, N, t = 1, \dots, T}$ : We use the recursive forward-filtering backward-sampling (FFBS) algorithm proposed by Baum et al.⁴⁴ to sample $Z_{i, t}$ from $p (Z_{i, t} | \cdot)$ . The FFBS algorithm is implemented as follows:

\begin{array}{l} p (Z_{i, t} | \cdot) \propto p (Y_{i}, Z_{i, t} | θ) = p (Y_{i, 1 : t}, Z_{i, t} | θ) p (Y_{i, t + 1 : T} | Z_{i, t}, θ) \\ ≐ q_{i, t} (Y_{i}, Z_{i, t} | θ) {\bar{q}}_{i, t} (Y_{i} | Z_{i, t}, θ) \end{array}

where

Y_{i, 1 : t}

denotes the observations up to time t for subject i. Moreover, we initialize

q_{i, 1} (Y_{i}, Z_{i, 1} | θ) = p (Y_{i, 1}, Z_{i, 1} | θ) = p (Y_{i, 1} | Z_{i, 1}, θ) p (Z_{i, 1} | θ)

and calculate

q_{i, t} (Y_{i}, Z_{i, t} | θ)

for

t = 2, \dots, T

in a recursive manner as follows:

\begin{array}{l} q_{i, t} (Y_{i}, Z_{i, t} | θ) = p (Y_{i, 1 : t}, Z_{i, t} | θ) = \sum_{u = 1}^{S} p (Y_{i, 1 : t}, Z_{i, t}, Z_{i, t - 1} = u | θ) \\ = \sum_{u = 1}^{S} {p (Y_{i, 1 : t - 1}, Z_{i, t - 1} = u | θ) p (Z_{i, t} | Z_{i, t - 1} = u, θ) p (Y_{i, t} | Z_{i, t}, θ)} \\ = \sum_{u = 1}^{S} {q_{i, t - 1} (Y_{i}, Z_{i, t - 1} = u | θ) p (Z_{i, t} | Z_{i, t - 1} = u, θ) p (Y_{i, t} | Z_{i, t}, θ)} \end{array}

where

p (Z_{i, t} | Z_{i, t - 1} = u, θ)

is the transition probability, and

p (Y_{i, t} | Z_{i, t}, θ)

is the likelihood in hidden state

Z_{i, t}

Similarly, we initialize ${\bar{q}}_{i, T} (Y_{i} | Z_{i, t}, θ) = 1$ and calculate ${\bar{q}}_{i, t} (Y_{i} | Z_{i, t}, θ)$ for $t = T - 1, \dots, 1$ as follows:

\begin{array}{l} {\bar{q}}_{i, t} (Y_{i} | Z_{i, t}, θ) = p (Y_{i, t + 1 : T} | Z_{i, t}, θ) = \sum_{u = 1}^{S} p (Y_{i, t + 1 : T}, Z_{i, t + 1} = u | Z_{i, t}, θ) \\ = \sum_{u = 1}^{S} {p (Y_{i, t + 2 : T} | Z_{i, t + 1} = u, θ) p (Z_{i, t + 1} = u | Z_{i, t}, θ) p (Y_{i, t + 1} | Z_{i, t + 1} = u, θ)} \\ = \sum_{u = 1}^{S} {{\bar{q}}_{i, t + 1} (Y_{i} | Z_{i, t + 1} = u, θ) p (Z_{i, t + 1} = u | Z_{i, t}, θ) p (Y_{i, t + 1} | Z_{i, t + 1} = u, θ)} \end{array}

where

p (Z_{i, t + 1} = u | Z_{i, t}, θ)

is the transition probability, and

p (Y_{i, t + 1} | Z_{i, t + 1} = u, θ)

is the likelihood in hidden state

Z_{i, t + 1} = u

Step 2: Update $V = {v_{i, t}, i = 1, \dots, N, t = 1, \dots, T}$ : Given $Z_{i, t} = s, Y_{i, t}$ , and $θ$ , the full conditional distribution of $v_{i, t}$ is

\begin{array}{l} p (V | \cdot) \propto v_{i, t}^{- 1 / 2} \exp {- \frac{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2} + 2 ξ_{2}^{2} v_{i, t}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}}}, \\ \propto v_{i, t}^{- 1 / 2} \exp [- {(\frac{ξ_{1}^{2}}{ξ_{2}^{2} σ_{s}} + \frac{2}{σ_{s}}) v_{i, t} + \frac{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t})}^{2}}{ξ_{2}^{2} σ_{s}} v_{i, t}^{- 1}} / 2] \end{array}

or equivalently,

[v_{i, t} | \cdot] \sim GIG ({\tilde{γ}}_{i, t}, {\tilde{ω}}_{i, t}, \frac{1}{2})

where

GIG (\cdot, \cdot)

denotes the generalized inverse Gaussian distribution,

{\tilde{γ}}_{i, t} = 2 / σ_{s} + ξ_{1}^{2} / (ξ_{2}^{2} σ_{s})

, and

{\tilde{ω}}_{i, t} = {(Y_{i, t} - β_{τ, s}^{T} X_{i, t})}^{2} / (ξ_{2}^{2} σ_{s})

Step 3: Update initial probability $P = {p_{1}, \dots, p_{S}}$ : Initial probability P appears in the likelihood as $\prod_{s = 1}^{S} p_{s}^{n_{s, 1}}$ . By incorporating the prior distribution $Dirichlet (a, \dots, a)$ , we obtain the full conditional distribution of P as

P = (p_{1}, p_{2}, \dots, p_{S}) \sim Dirichlet (a + n_{1, 1}, \dots, a + n_{S, 1})

Step 4: Update transition parameters $ζ = {ζ_{u, s}, u = 1, \dots, S, s = 1, \dots, S - 1}$ : In the likelihood function of the transition model, $ζ_{u, s}$ appears in $p_{i, t; u, v}$ for $v \geq s$ . By incorporating the prior distribution $N (ζ_{u, s}^{0}, δ_{ζ}^{2})$ , we derive the full conditional distribution of $ζ_{u, s}$ as

p (ζ_{u, s} | \cdot) \propto \exp {\sum_{v = s}^{S} \sum_{i = 1}^{N} \sum_{t = 2}^{T} Z_{i, t - 1} (u) Z_{i, t} (v) \log (p_{i, t; u, v}) - \frac{{(ζ_{u, s} - ζ_{u, s}^{0})}^{2}}{2 δ_{ζ}^{2}}}

Step 5: Update transition parameter $α$ : Parameter $α$ appears in $p_{i, t; u, s}$ [see (4) or (5)]. By incorporating the prior distribution $N (α^{0}, Σ_{α})$ , we obtain the full conditional distribution of $α$ as

p (α | \cdot) \propto \exp {\sum_{u = 1}^{S} \sum_{s = 1}^{S} \sum_{i = 1}^{N} \sum_{t = 2}^{T} Z_{i, t - 1} (u) Z_{i, t} (s) \log (p_{i, t; u, s}) - \frac{1}{2} {(α - α^{0})}^{T} Σ_{α}^{- 1} (α - α^{0})}

Step 6: Update QR coefficient $β = {β_{τ, s}, s = 1, \dots, S}$ : Parameter $β_{τ, s}$ appears in the likelihood as

\exp {- \sum_{i = 1}^{N} \sum_{t = 1}^{T} \frac{Z_{i, t} (s) {(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}}}

By incorporating the prior distribution $N (β_{τ, s}^{0}, B_{τ, s}^{0})$ , we derive the full conditional distribution of $β_{τ, s}$ as

\exp {- \sum_{i = 1}^{N} \sum_{t = 1}^{T} \frac{Z_{i, t} (s) {(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}} - \frac{1}{2} {(β_{τ, s} - β_{τ, s}^{0})}^{T} B_{τ, s}^{0}^{- 1} (β_{τ, s} - β_{τ, s}^{0})}

That is,

[β_{τ, s} | \cdot] \sim N ({\tilde{β}}_{τ, s}, {\tilde{B}}_{τ, s})

where

\begin{array}{l} {\tilde{B}}_{τ, s}^{- 1} = \sum_{i = 1}^{N} \sum_{t = 1}^{T} Z_{i, t} (s) X_{i, t} X_{i, t}^{'} / (ξ_{2}^{2} σ_{s} v_{i, t}) + B_{τ, s}^{0}^{- 1}, \\ {\tilde{β}}_{τ, s} = {\tilde{B}}_{τ, s} {\sum_{i = 1}^{N} \sum_{t = 1}^{T} Z_{i, t} (s) X_{i, t} (Y_{i, t} - ξ_{1} v_{i, t}) / (ξ_{2}^{2} σ_{s} v_{i, t}) + B_{τ, s}^{0}^{- 1} β_{τ, s}^{0}} \end{array}

Step 7: Update scale parameter $σ = {σ_{1}, \dots, σ_{S}}$ : Parameter σ_s appears in the likelihood as

\begin{array}{l} Π_{i = 1}^{N} Π_{t = 1}^{T} {[σ_{s}^{- 3 / 2} \exp {- \frac{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2} + 2 ξ_{2}^{2} v_{i, t}^{2}}{2 ξ_{2}^{2} σ_{s} v_{i, t}}}]}^{Z_{i, t} (s)} \\ = σ_{s}^{- 3 n_{s} / 2} \exp [- σ_{s}^{- 1} \sum_{i = 1}^{N} \sum_{t = 1}^{T} Z_{i, t} {^{{(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2}} / (2 ξ_{2}^{2} v_{i, t}) + v_{i, t}}] \end{array}

By incorporating the prior distribution $I G (n_{s}^{0} / 2, b_{s}^{0} / 2)$ , we acquire the full conditional distribution of σ_s as

[σ_{s} | \cdot] \sim I G ({\tilde{n}}_{s} / 2, {\tilde{b}}_{s} / 2)

where

{\tilde{n}}_{s} = n_{s}^{0} + 3 n_{s}

and

{\tilde{b}}_{s} = b_{s}^{0} + 2 \sum_{i = 1}^{N} \sum_{t = 1}^{T} Z_{i, t} (s) v_{i, t} + \sum_{i = 1}^{N} \sum_{t = 1}^{T} Z_{i, t} (s) {(Y_{i, t} - β_{τ, s}^{'} X_{i, t} - ξ_{1} v_{i, t})}^{2} / (ξ_{2}^{2} v_{i, t})

4. Simulation study

4.1 Simulation 1

In this section, we assess the finite sample performance of the proposed method through numerical studies. We consider an HMM with S = 2 as follows:

\begin{array}{l} [Y_{i, t} | Z_{i, t} = s] = β_{s}^{'} X_{i, t} + ε_{i, t}, s = 1, 2 \\ logit (η_{i, t; u, 1}) = \log (p_{i, t; u, 1} / p_{i, t; u, 2}) = ζ_{u, 1} + α^{'} b_{i, t}, u = 1, 2 \end{array}

(9)

where

b_{i, t} = {(b_{i, t, 1}, b_{i, t, 2})}^{'}

, in which

b_{i, t, 1}

and

b_{i, t, 2}

are independently generated from Bernoulli (0.5);

X_{i, t} = {(1, X_{i, t, 1}, X_{i, t, 2})}^{'}

, in which

X_{i, t, 1}

and

X_{i, t, 2}

are independently generated from N(0, 1);

β_{1} = {(β_{1, 1}, β_{1, 2}, β_{1, 3})}^{'} = {(- 2.0, - 2.0, - 2.0)}^{'}

and

β_{2} = {(β_{2, 1}, β_{2, 2}, β_{2, 3})}^{'} = {(2.0, 2.0, 2.0)}^{'}

;

ζ_{1, 1} = ζ_{2, 1} = - 0.1, α = (α_{1}, α_{2})' = (- 0.5, 0.5)'

; and the parameters for the initial distribution of the hidden states are

(p_{1}, p_{2}) = (0.4, 0.6)

The hyperparameters of the prior distributions in (8) are assigned as follows (Prior 1): a = 1; $ζ_{u, s}^{0} = 0, δ_{ζ}^{2} = 1.0$ ; $α^{0} = {(0, 0)}^{'}, Σ_{α} = I$ ; $β_{τ, s}^{0} = {(0, 0, 0)}^{'}, B_{τ, s}^{0} = 100 I$ ; $n_{s}^{0} = 3, b_{s}^{0} = 0.1$ ; and s = 1, 2. We consider three different error distributions, namely, N(0, 1), t(3), and $Γ (1, 2)$ , where $Γ (1, 2)$ denotes the gamma distribution with shape parameter 1 and scale parameter 2 (centralized by subtracting its mean value) and four combinations of (N, T) as follows: (200, 4), (400, 4), (400, 8), and (800, 8). A total of 1000 replicated datasets are generated for each scenario.

We use estimated potential scale reduction (EPSR)⁴⁵ to check the convergence of the MCMC algorithm. The MCMC algorithm converges within 2000 iterations. Thus, we discard 2000 burn-ins and collect the subsequent 10,000 posterior samples for posterior inference. Tables 1 to 3 present the biases (Bias) and root mean squared errors (RMSE) between the Bayesian estimates and true population values of the unknown parameters at τ = 0.5, 0.25, and 0.75, respectively. Bias and RMSE are close to zero in all the considered scenarios, thus revealing the good performance of the proposed method regardless of the distribution of $ε_{i, t}$ and the sizes of N and T. Figure 2 presents boxplots of the summations of the absolute biases of all unknown parameters under various combinations of (N, T). As expected, the performance of the estimates improves when N and/or T increases.

Table 1.

Bayesian estimates of the parameters under $τ = 0.50$ .

Error	(N, T)	(200,4)		(400,4)		(400,8)		(800,8)
Error	Par	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
N(0, 1)	p ₁	0.008	0.041	0.002	0.029	0.001	0.027	0.001	0.018
	p ₂	–0.008	0.041	–0.002	0.029	–0.001	0.027	–0.001	0.018
	$β_{1, 1}$	0.012	0.060	0.007	0.040	0.003	0.028	0.001	0.021
	$β_{1, 2}$	0.006	0.056	0.007	0.041	0.001	0.027	0.002	0.020
	$β_{1, 3}$	0.003	0.057	0.002	0.040	0.000	0.028	0.001	0.019
	$β_{2, 1}$	0.002	0.056	0.002	0.037	0.000	0.026	0.000	0.018
	$β_{2, 2}$	–0.005	0.051	–0.002	0.036	–0.001	0.026	0.001	0.019
	$β_{2, 3}$	0.000	0.048	–0.002	0.037	0.001	0.028	0.000	0.019
	$ζ_{1, 1}$	0.000	0.173	0.003	0.126	–0.001	0.089	0.006	0.062
	$ζ_{2, 1}$	0.024	0.135	0.028	0.109	0.018	0.074	0.024	0.055
	α ₁	0.014	0.174	0.004	0.133	0.005	0.088	–0.002	0.065
	α ₂	–0.016	0.185	–0.011	0.141	–0.002	0.093	–0.004	0.067
t(3)	p ₁	0.009	0.039	0.011	0.031	0.008	0.028	0.008	0.021
	p ₂	–0.009	0.039	–0.011	0.031	–0.008	0.028	–0.008	0.021
	$β_{1, 1}$	–0.057	0.141	–0.059	0.096	–0.056	0.077	–0.056	0.068
	$β_{1, 2}$	–0.019	0.106	–0.021	0.078	–0.030	0.060	–0.026	0.042
	$β_{1, 3}$	–0.025	0.118	–0.021	0.073	–0.027	0.057	–0.028	0.044
	$β_{2, 1}$	0.054	0.103	0.061	0.088	0.059	0.074	0.060	0.068
	$β_{2, 2}$	0.029	0.091	0.032	0.066	0.033	0.055	0.031	0.044
	$β_{2, 3}$	0.024	0.090	0.031	0.068	0.032	0.053	0.033	0.045
	$ζ_{1, 1}$	0.014	0.169	–0.015	0.126	0.011	0.083	0.010	0.060
	$ζ_{2, 1}$	0.030	0.140	0.025	0.113	0.026	0.072	0.027	0.054
	α ₁	0.046	0.187	0.043	0.132	0.038	0.090	0.035	0.073
	α ₂	–0.059	0.186	–0.030	0.148	–0.039	0.096	–0.037	0.074
$Γ (1, 2)$	p ₁	–0.004	0.040	–0.001	0.029	–0.004	0.028	–0.002	0.020
	p ₂	0.004	0.040	0.001	0.029	0.004	0.028	0.002	0.020
	$β_{1, 1}$	0.008	0.029	0.008	0.021	0.006	0.016	–0.005	0.010
	$β_{1, 2}$	0.012	0.031	0.011	0.024	0.009	0.017	0.002	0.009
	$β_{1, 3}$	0.019	0.035	0.012	0.023	0.010	0.018	0.013	0.016
	$β_{2, 1}$	–0.017	0.032	–0.017	0.025	–0.017	0.022	–0.024	0.025
	$β_{2, 2}$	–0.009	0.029	–0.007	0.019	–0.004	0.013	–0.011	0.013
	$β_{2, 3}$	–0.007	0.027	–0.005	0.018	–0.004	0.014	–0.001	0.007
	$ζ_{1, 1}$	–0.015	0.172	–0.003	0.129	–0.006	0.082	–0.001	0.065
	$ζ_{2, 1}$	0.018	0.138	0.010	0.101	0.011	0.068	0.019	0.054
	α ₁	–0.019	0.176	–0.008	0.133	–0.016	0.087	–0.020	0.068
	α ₂	0.010	0.186	0.007	0.137	0.023	0.089	0.019	0.068

RMSE: root mean squared error.

Table 2.

Bayesian estimates of the parameters under $τ = 0.25$ .

Error	(N,T)	(200,4)		(400,4)		(400,8)		(800,8)
Error	Par	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
N(0, 1)	p ₁	0.004	0.041	0.002	0.029	0.002	0.028	0.003	0.020
	p ₂	–0.004	0.041	–0.002	0.029	–0.002	0.028	–0.003	0.020
	$β_{1, 1}$	–0.029	0.063	–0.013	0.042	–0.006	0.029	–0.001	0.020
	$β_{1, 2}$	0.008	0.057	0.005	0.040	0.002	0.027	0.001	0.019
	$β_{1, 3}$	0.008	0.058	0.007	0.041	0.002	0.028	0.001	0.020
	$β_{2, 1}$	–0.023	0.056	–0.011	0.039	–0.005	0.028	0.001	0.017
	$β_{2, 2}$	–0.002	0.050	–0.002	0.036	–0.002	0.026	–0.001	0.016
	$β_{2, 3}$	–0.005	0.052	–0.001	0.035	0.003	0.027	0.001	0.019
	$ζ_{1, 1}$	–0.004	0.171	0.012	0.126	–0.003	0.080	0.000	0.054
	$ζ_{2, 1}$	0.023	0.141	0.023	0.104	0.019	0.063	0.017	0.046
	α ₁	0.015	0.178	0.003	0.129	0.011	0.087	0.005	0.065
	α ₂	–0.008	0.183	–0.017	0.129	–0.008	0.089	0.002	0.069
t(3)	p ₁	0.008	0.038	0.009	0.029	0.009	0.029	0.009	0.021
	p ₂	–0.008	0.038	–0.009	0.029	–0.009	0.029	–0.009	0.021
	$β_{1, 1}$	–0.091	0.136	–0.072	0.106	–0.066	0.084	–0.059	0.070
	$β_{1, 2}$	–0.026	0.101	–0.021	0.075	–0.024	0.055	–0.026	0.044
	$β_{1, 3}$	–0.017	0.101	–0.025	0.075	–0.028	0.056	–0.030	0.047
	$β_{2, 1}$	0.031	0.093	0.048	0.078	0.052	0.069	0.058	0.067
	$β_{2, 2}$	0.029	0.092	0.030	0.067	0.031	0.054	0.033	0.046
	$β_{2, 3}$	0.035	0.088	0.033	0.071	0.033	0.053	0.033	0.045
	$ζ_{1, 1}$	0.001	0.176	0.005	0.131	0.006	0.085	0.014	0.062
	$ζ_{2, 1}$	0.024	0.135	0.025	0.107	0.024	0.070	0.029	0.059
	α ₁	0.041	0.181	0.030	0.135	0.043	0.095	0.036	0.073
	α ₂	–0.057	0.184	–0.038	0.137	–0.040	0.099	–0.042	0.077
$Γ (1, 2)$	p ₁	–0.003	0.038	–0.004	0.029	–0.001	0.029	0.000	0.019
	p ₂	0.003	0.038	0.004	0.029	0.001	0.029	0.000	0.019
	$β_{1, 1}$	–0.021	0.035	–0.005	0.021	0.000	0.015	–0.007	0.010
	$β_{1, 2}$	0.012	0.031	0.009	0.022	0.010	0.017	0.002	0.008
	$β_{1, 3}$	0.013	0.032	0.010	0.022	0.008	0.017	0.013	0.016
	$β_{2, 1}$	–0.041	0.048	–0.026	0.032	–0.022	0.026	–0.025	0.026
	$β_{2, 2}$	–0.004	0.025	–0.004	0.019	–0.005	0.014	–0.010	0.012
	$β_{2, 3}$	–0.005	0.026	–0.004	0.019	–0.004	0.012	–0.001	0.007
	$ζ_{1, 1}$	–0.014	0.172	0.002	0.122	–0.005	0.084	–0.002	0.062
	$ζ_{2, 1}$	0.010	0.141	0.014	0.099	0.009	0.070	0.012	0.045
	α ₁	–0.018	0.178	–0.021	0.130	–0.024	0.086	–0.021	0.064
	α ₂	0.009	0.190	0.022	0.133	0.010	0.087	0.017	0.058

RMSE: root mean squared error.

Table 3.

Bayesian estimates of the parameters under $τ = 0.75$ .

Error	(N,T)	(200,4)		(400,4)		(400,8)		(800,8)
Error	Par	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
N(0, 1)	p ₁	0.008	0.041	0.002	0.028	0.004	0.029	0.002	0.020
	p ₂	–0.008	0.041	–0.002	0.028	–0.004	0.029	–0.002	0.020
	$β_{1, 1}$	0.033	0.068	0.013	0.042	0.008	0.029	0.002	0.020
	$β_{1, 2}$	0.000	0.058	0.003	0.040	–0.001	0.026	0.000	0.018
	$β_{1, 3}$	0.004	0.056	0.001	0.040	0.001	0.027	0.001	0.020
	$β_{2, 1}$	0.022	0.058	0.013	0.042	0.005	0.027	0.002	0.019
	$β_{2, 2}$	–0.004	0.052	–0.002	0.038	0.000	0.026	0.002	0.019
	$β_{2, 3}$	–0.004	0.053	–0.003	0.038	0.000	0.027	–0.002	0.019
	$ζ_{1, 1}$	0.018	0.175	–0.001	0.126	0.003	0.085	0.009	0.062
	$ζ_{2, 1}$	0.028	0.139	0.024	0.104	0.024	0.072	–0.025	0.058
	α ₁	0.006	0.184	0.007	0.135	0.005	0.087	0.001	0.066
	α ₂	–0.019	0.187	–0.005	0.133	–0.004	0.092	–0.003	0.063
t(3)	p ₁	0.010	0.040	0.008	0.029	0.009	0.030	0.010	0.021
	p ₂	–0.010	0.040	–0.008	0.029	–0.009	0.030	–0.010	0.021
	$β_{1, 1}$	–0.021	0.107	–0.041	0.089	–0.053	0.077	–0.055	0.067
	$β_{1, 2}$	–0.023	0.103	–0.027	0.076	–0.029	0.062	–0.028	0.045
	$β_{1, 3}$	–0.024	0.102	–0.027	0.081	–0.030	0.057	–0.029	0.046
	$β_{2, 1}$	0.087	0.124	0.071	0.095	0.060	0.075	0.060	0.068
	$β_{2, 2}$	0.028	0.087	0.028	0.065	0.030	0.053	0.035	0.047
	$β_{2, 3}$	0.022	0.087	0.030	0.067	0.030	0.053	0.031	0.044
	$ζ_{1, 1}$	0.008	0.172	0.015	0.127	0.012	0.075	0.013	0.063
	$ζ_{2, 1}$	0.029	0.142	0.026	0.105	0.022	0.068	0.028	0.058
	α ₁	0.051	0.184	0.038	0.133	0.040	0.090	0.036	0.073
	α ₂	–0.062	0.193	–0.040	0.135	–0.035	0.090	–0.037	0.072
$Γ (1, 2)$	p ₁	–0.001	0.039	–0.004	0.030	–0.002	0.027	0.000	0.023
	p ₂	0.001	0.039	0.004	0.030	0.002	0.027	0.000	0.023
	$β_{1, 1}$	0.036	0.046	0.016	0.027	0.010	0.017	–0.005	0.009
	$β_{1, 2}$	0.009	0.030	0.008	0.025	0.008	0.016	0.002	0.009
	$β_{1, 3}$	0.008	0.030	0.006	0.021	0.008	0.016	0.015	0.016
	$β_{2, 1}$	0.006	0.027	–0.008	0.020	–0.012	0.018	–0.025	0.026
	$β_{2, 2}$	–0.006	0.026	–0.003	0.020	–0.004	0.013	–0.010	0.012
	$β_{2, 3}$	–0.006	0.026	–0.005	0.019	–0.004	0.013	–0.001	0.007
	$ζ_{1, 1}$	0.001	0.167	–0.010	0.114	–0.005	0.083	0.011	0.061
	$ζ_{2, 1}$	0.013	0.139	0.013	0.010	0.018	0.071	0.019	0.054
	α ₁	–0.017	0.178	–0.023	0.133	–0.020	0.092	–0.027	0.070
	α ₂	0.019	0.185	0.015	0.133	0.017	0.090	0.010	0.071

RMSE: root mean squared error.

Figure 2.

Boxplots of the summations of the absolute biases of all unknown parameters under various combinations of (N, T) in the simulation study.

Moreover, we conduct a sensitivity analysis to determine if the Bayesian results are sensitive to the prior specification. We disturb the hyperparameters in Prior 1 as follows (Prior 2): a = 2; $ζ_{u, s}^{0} = 2, δ_{ζ}^{2} = 4.0$ ; $α^{0} = {(0, 0)}^{'}, Σ_{α} = I$ ; $β_{τ, s}^{0} = {(2, 2, 2)}^{'}, B_{τ, s}^{0} = 10 I$ ; $n_{s}^{0} = 1, b_{s}^{0} = 0.05$ ; and s = 1, 2. Then, we reanalyze the 1000 replicated datasets that are simulated in the scenarios of $(N, T) = (200, 4), τ = 0.5$ , and $ε_{i, t} \sim N (0, 1)$ , t(3), or $Γ (1, 2)$ . The estimation results and those obtained under Prior 1 are presented in Table A1 of the Online Appendix. The estimation performance under Priors 1 and 2 is similar. Thus, the Bayesian estimates are insensitive to the prior inputs under consideration.

As suggested by an anonymous referee, we compare the results obtained by the proposed method with those obtained by a standard HMM assuming the conditional distribution of $Y_{i, t}$ . Table A2 of the Online Appendix shows that the two methods perform similarly, except for the intercepts in the conditional model at $τ = 0.25$ and 0.75. This finding coincides with the referee’s point that the simulation focuses on a location model where the error distributions are allowed to be nonnormal and the effect of covariates on the mean is the same across any particular quantile; therefore, a naive model assuming a normal distribution for the errors is expected to perform as well as the location model.

To demonstrate the capability of the quantile HMM in revealing quantile-specific covariate effects, we consider another scenario where the effect of covariates varies with τ as follows: at $τ = 0.25, β_{1} = {(- 1.9, 0.8, 1.5)}^{'}, β_{2} = {(1.7, - 1.8, - 1.3)}^{'}$ ; at $τ = 0.5, β_{1} = {(- 2.2, - 1.4, - 1.1)}^{'}, β_{2} = {(1.6, 1.3, 1.7)}^{'}$ ; and at $τ = 0.75, β_{1} = {(- 1.8, 0.9, 1.5)}^{'}, β_{2} = (2.4, - 1.7,$ $- 0.8)^{'}$ . The parameters in the transition model are set the same as before. Table A3 of the Online Appendix presents the estimation results under $(N, T) = (200, 4), ε_{i, t} \sim Γ (1, 2)$ , and Prior 1. The results demonstrate that the proposed model satisfactorily estimates quantile-specific covariate effects that cannot be discovered via mean regression.

4.2 Simulation 2

In this section, we expand the proposed model by relaxing the assumption of ordered hidden states. The conditional model is defined by (1) and (2), and the transition model is a multinomial logit model as follows:

\begin{array}{l} logit (η_{i, t; u, s}) = \log \frac{P (Z_{i, t} = s | Z_{i, t - 1} = u)}{P (Z_{i, t} = S | Z_{i, t - 1} = u)} = \log \frac{p_{i, t; u, s}}{p_{i, t; u, S}} = ζ_{u, s} + α_{u, s}^{'} b_{i, t}, \\ i = 1, \dots, N, t = 2, \dots, T, u = 1, \dots, S, s = 1, \dots, S - 1 \end{array}

(10)

where the hidden states are unordered, and

ζ_{u, S}

and

α_{u, S}

are fixed to 0 for identification purpose. On the basis of model (10), the transition probabilities can be obtained as follows:

\begin{array}{l} p_{i, t; u, s} = \frac{\exp (ζ_{u, s} + α_{u, s}^{'} b_{i, t})}{1 + \sum_{k = 1}^{S - 1} \exp (ζ_{u, k} + α_{u, s}^{'} b_{i, t})}, s = 1, \dots, S - 1, \\ p_{i, t; u, S} = \frac{1}{1 + \sum_{k = 1}^{S - 1} \exp (ζ_{u, k} + α_{u, s}^{'} b_{i, t})} \end{array}

(11)

With this transition model, the posterior sampling of the parameters in (11) is slightly modified. For example, the full conditional distribution of $ζ_{u, s}$ and $α_{u, s}$ in Steps 4 and 5 become the following:

p (ζ_{u, s} | \cdot) \propto \exp {\sum_{v = 1}^{S} \sum_{i = 1}^{N} \sum_{t = 2}^{T} Z_{i, t - 1} (u) Z_{i, t} (v) \log (p_{i, t; u, v}) - \frac{{(ζ_{u, s} - ζ_{u, s}^{0})}^{2}}{2 δ_{ζ}^{2}}}

and

p (α_{u, s} | \cdot) \propto \exp {\sum_{v = 1}^{S} \sum_{i = 1}^{N} \sum_{t = 2}^{T} Z_{i, t - 1} (u) Z_{i, t} (v) \log (p_{i, t; u, v}) - \frac{1}{2} {(α_{u, s} - α_{u, s}^{0})}^{T} Σ_{α_{u, s}}^{- 1} (α_{u, s} - α_{u, s}^{0})}

where

p_{i, t; u, s}

is defined by (11).

We consider S = 3 and other settings as follows: $b_{i, t} = {(b_{i, t, 1}, b_{i, t, 2})}^{'}$ and $X_{i, t} = (1, X_{i, t, 1},$ $X_{i, t, 2})^{'}$ are generated in the same manner as before. The parameters in the conditional and transitional models are set as follows: $β_{1} = {(β_{1, 1}, β_{1, 2}, β_{1, 3})}^{'} = {(- 2.0, - 2.0, - 2.0)}^{'}, β_{2} = {(β_{2, 1}, β_{2, 2}, β_{2, 3})}^{'} = {(0.5, 0.5, 0.5)}^{'}, β_{3} = {(β_{3, 1}, β_{3, 2}, β_{3, 3})}^{'} = {(2.0, 2.0, 2.0)}^{'}; ζ_{1, 1} = - 1.5, ζ_{1, 2} = - 0.5, ζ_{2, 1} = 0.5, ζ_{2, 2} = 0.5, ζ_{3, 1} = 0.5, ζ_{3, 2} = - 1.0; α_{1, 1} = (α_{1, 1, 1}, α_{1, 1, 2})' = (- 0.5, 0.5)', α_{1, 2} = (α_{1, 2, 1}, α_{1, 2, 2})' = (- 1.0, 1.0)', α_{2, 1} = (α_{2, 1, 1}, α_{2, 1, 2})' = (0.5, - 0.5)', α_{2, 2} = (α_{2, 2, 1}, α_{2, 2, 2})' = (1.0, - 1.0)', α_{3, 1} = (α_{3, 1, 1}, α_{3, 1, 2})' = (0.5, 0.5)', α_{3, 2} = (α_{3, 2, 1}, α_{3, 2, 2})' = (1.0, 1.0)'$ The parameters for the initial distribution of hidden states are $(p_{1}, p_{2}, p_{3}) = (0.3, 0.3, 0.4)$ .

A total of 1000 replicated datasets are generated from the abovementioned model under $(N, T) = (800, 8)$ and $ε_{i, t} \sim Γ (1, 2)$ . The hyperparameters are assigned in a similar manner as in Prior 1. Table A4 of the Online Appendix presents the estimation results. The performance of the parameter estimates is acceptable but worse than that reported in Tables 1 to 3. The main reason is that the number of parameters dramatically increases from 12 (in Simulation 1) to 30 (in Simulation 2). In particular, the number of regression coefficients $α$ in the transition model increases from 2 (in Simulation 1) to 12 (in Simulation 2). Hence, the expanded quantile HMM can accommodate unordered hidden states but may require a large sample size due to the increased complexity of the transition model.

In practice, we can determine whether the underlying states have a natural order on the basis of subject knowledge, experts’ experience, and/or common sense. For HMMs with ordered hidden states, the continuation-ratio logit model (4) is recommended to consider the ranking information of the hidden states, thereby facilitating a parsimonious transition model and simple interpretation. Otherwise, the multinomial logit model (10) is suggested to model the unordered hidden states.

5 Longitudinal analysis of cocaine use

We applied the proposed quantile HMM to analyze a longitudinal study concerning the prevention of cocaine use and its influential factors. The dataset was collected from the University of California, Los Angeles, Center for Advancing Longitudinal Drug Abuse Research. A total of 321 participants admitted in the West Los Angeles Veterans Affairs Medical Center between 1988 to 1989 were assessed at the baseline, 1 year after treatment, 2 years after treatment, and 12 years after treatment in 2002–2003; thus, T = 4. After deleting individuals with missing entries, N = 298 samples remained. The response of our primary interest, $Y_{i, t}$ , was measured as days of cocaine use per month, and the main goal was to investigate the potential influential factors of $Y_{i, t}$ . Song et al.¹² used a mean HMM with latent variables to analyze this dataset. They identified three cocaine use states, namely, addiction to, moderate dependence on, and swearing off cocaine, and found that the effects of patients’ treatment received and mental health on the mean of $Y_{i, t}$ vary across the three states. In this study, we applied the proposed quantile HMM to examine the potential covariate effects on the distribution of $Y_{i, t}$ , especially on the three quartiles of $Y_{i, t}$ . We selected the three quartiles ( $τ = 0.25, 0.5, 0.75$ ) because they are the most representative statistics to summarize the lower, central, and upper tendencies of the response variable, thereby enabling us to understand how potential predictors, such as treatment received and mental health, influence the response of patients whose cocaine use is at low, median, or high level given one specific state. Figure 1 presents the histograms of $Y_{i, t}$ at different time points, and it indicates that the distribution of $Y_{i, t}$ is highly skewed. The proposed quantile HMM accommodates the nonnormality of $Y_{i, t}$ and provides a complete picture of the associations between $Y_{i, t}$ and its influential factors.

We considered the following covariates in the conditional QR and continuation-logit transition models.

$X_{i, t, 1}$ : time per month in formal treatment with integer values in the range [0, 30];

$X_{i, t, 2}$ : depression, a positive continuous variable measuring the depression level;

$b_{i, t, 1}$ : employment, a binary variable taking the value of 1 if employed and 0 otherwise;

$b_{i, t, 2}$ : have someone to talk about the problem, a binary variable taking the value of 1 if a person has it and 0 otherwise.

Among these covariates, $X_{i, t, 1}$ and $X_{i, t, 2}$ are directly associated with cocaine use behavior and thus included in the conditional QR model, whereas $b_{i, t, 1}$ and $b_{i, t, 2}$ are typical individual characteristics and therefore included in the transition model. Considering that $Y_{i, t}$ has a restricted range of 0 to 30 days, we transformed $Y_{i, t}$ to ${\tilde{Y}}_{i, t} = \log \frac{Y_{i, t} + 0.5}{30 - Y_{i, t} + 0.5}$ such that ${\tilde{Y}}_{i, t}$ has an approximate range of $(- \infty, \infty)$ . Then, ${\tilde{Y}}_{i, t}$ and the continuous covariates $X_{i, t, 1}$ and $X_{i, t, 2}$ were standardized prior to analysis.

In accordance with a previous study,¹² we fixed S = 3 and regarded the three hidden states as severe, moderate, and minor cocaine addiction states. Then, the proposed model can be written as

\begin{array}{l} [Y_{i, t} | Z_{i, t} = s] = β_{s}^{'} X_{i, t} + ε_{i, t}, \\ logit (η_{i, t, u, s}) = ζ_{u, s} + α^{'} b_{i, t}, u = 1, \dots, S, s = 1, \dots, S - 1 \end{array}

where

X_{i, t} = {(1, X_{i, t, 1}, X_{i, t, 2})}^{'}, b_{i, t} = {(b_{i, t, 1}, b_{i, t, 2})}^{'}, β_{s} = {(β_{s, 1}, β_{s, 2}, β_{s, 3})}^{'}

, and

α = {(α_{1}, α_{2})}^{'}

. On the basis of the meaning of

Y_{i, t}

, the higher

Y_{i, t}

is, the worse the state is. Thus, we imposed the identifiability restriction

β_{τ, 1, 1} > β_{τ, 2, 1} > β_{τ, 3, 1}

in each MCMC iteration to avoid label switching and ensure that the states are ranked from the worst to the best. Two prior inputs were assigned in the same manner as that in the simulation study. We ran three parallel MCMC chains using different initial values under Prior 1 and calculated the EPSR values of the parameters to check the convergence of the MCMC algorithm.⁴⁶ Figure 3 shows that all of the EPSR values become stabilized and are less than 1.2 within 4000 iterations. Hence, we collected 10,000 posterior samples after discarding 4000 burn-in iterations to conduct posterior inference.

Figure 3.

Plot of EPSR values in the analysis of cocaine usage study under Prior 1.

Tables 4 and 5 present the estimation results under Priors 1 and 2, respectively. The results obtained with different priors are similar. For the conditional model, we derived the following observations. (i) Formal treatment ( $X_{i, t, 1}$ ) reduces the days of cocaine use per month, whereas depression ( $X_{i, t, 2}$ ) increases the days of cocaine use per month. This result agrees with the finding of Song et al.,¹² who initially formulated treatment and mental health as latent traits through their multiple surrogates and then investigated how the latent traits affected the mean of the response conditional on hidden states. (ii) The magnitudes of the covariate effects are the largest in the severe state (state 1) but gradually decrease from a moderate (state 2) to a minor cocaine addiction state (state 3). (iii) The covariate effects at the tails are slightly different from those at the median. In most of the cases, the magnitudes of the covariate effects slightly decline when τ varies from 0.5 to 0.25 and then to 0.75. This result provides a more complete picture of the associations between cocaine use and their influential factors than mean regression does. For the transition model, $b_{i, t, 1}$ and $b_{i, t, 2}$ have apparent negative effects on $logit (η_{i, t, u, s})$ ((4)) over the three quartiles, and the effects at the tails are slightly less pronounced than those at the median. This result implies that being employed and/or having someone to talk about the problem increases the likelihood of transitioning from a certain state to a better (less severe) one. Such subject-specific influential factors elicit heterogeneity in between-state transitions and cannot be discovered by existing homogeneous HMMs.

Table 4.

Parameter estimates in the cocaine usage study under Prior 1.

	$τ = 0.25$		$τ = 0.5$		$τ = 0.75$
Par	Est	Std	Est	Std	Est	Std
p _severe	0.481	0.030	0.479	0.031	0.478	0.031
p_moderate	0.280	0.029	0.281	0.031	0.279	0.029
p_minor	0.239	0.029	0.240	0.030	0.243	0.028
intercept _severe	0.988	0.293	3.668	0.345	6.385	0.296
treatment _severe	–1.049	0.312	–1.064	0.376	–1.061	0.316
depression _severe	1.147	0.294	1.172	0.358	1.144	0.298
intercept _moderate	–2.795	0.361	–0.085	0.426	2.651	0.355
treatment _moderate	–0.866	0.365	–0.895	0.442	–0.848	0.366
depression _moderate	0.939	0.378	0.975	0.450	0.923	0.366
intercept _minor	–6.007	0.223	–3.310	0.291	–0.636	0.224
treatment _minor	–0.764	0.222	–0.776	0.289	–0.758	0.218
depression _minor	0.606	0.234	0.619	0.307	0.607	0.231
${intercept}_{severe \to severe}$	0.418	0.152	0.411	0.158	0.412	0.151
${intercept}_{severe \to moderate}$	0.405	0.212	0.389	0.221	0.403	0.214
${intercept}_{moderate \to severe}$	–0.495	0.193	–0.497	0.191	–0.493	0.194
${intercept}_{moderate \to moderate}$	0.266	0.208	0.285	0.233	0.267	0.217
${intercept}_{minor \to severe}$	–1.806	0.234	–1.802	0.244	–1.818	0.241
${intercept}_{minor \to moderate}$	–0.740	0.199	–0.748	0.209	–0.743	0.187
$employment$	–0.780	0.127	–0.785	0.135	–0.783	0.138
$have someone to talk$	–0.764	0.150	–0.773	0.132	–0.744	0.135

Table 5.

Parameter estimates in the cocaine usage study under Prior 2.

	$τ = 0.25$		$τ = 0.5$		$τ = 0.75$
Par	Est	Std	Est	Std	Est	Std
p_severe	0.479	0.030	0.478	0.031	0.476	0.030
p_moderate	0.280	0.029	0.280	0.030	0.280	0.028
p_minor	0.241	0.029	0.242	0.030	0.244	0.028
intercept _severe	1.002	0.289	3.673	0.343	6.387	0.298
treatment _severe	–1.037	0.305	–1.071	0.369	–1.061	0.316
depression _severe	1.166	0.293	1.195	0.358	1.140	0.299
intercept _moderate	–2.754	0.358	–0.084	0.413	2.665	0.356
treatment _moderate	–0.854	0.363	–0.863	0.432	–0.829	0.369
depression _moderate	0.964	0.377	0.984	0.446	0.907	0.375
intercept _minor	–5.977	0.224	–3.256	0.287	–0.641	0.224
treatment _minor	–0.755	0.225	–0.782	0.285	–0.766	0.220
depression _minor	0.611	0.236	0.620	0.299	0.608	0.230
${intercept}_{severe \to severe}$	0.432	0.165	0.434	0.168	0.447	0.161
${intercept}_{severe \to moderate}$	0.423	0.226	0.404	0.227	0.436	0.222
${intercept}_{moderate \to severe}$	–0.493	0.199	–0.499	0.201	–0.478	0.194
${intercept}_{moderate \to moderate}$	0.282	0.231	0.256	0.243	0.296	0.230
${intercept}_{minor \to severe}$	–1.830	0.246	–1.833	0.261	–1.844	0.250
${intercept}_{minor \to moderate}$	–0.757	0.201	–0.731	0.204	–0.738	0.190
employment	–0.776	0.147	–0.808	0.138	–0.806	0.142
have someone to talk	–0.747	0.145	–0.775	0.147	–0.728	0.148

Furthermore, for any given subject i and time t, a transition matrix $P (i, t) = {[p_{i, t, u, s}]}_{3 \times 3}$ can be estimated based on (5). For example, when i = 1 and t = 2 , b_{1,2,1}=1, b_{1,2,2}=0, and

\hat{P} (1, 2) = [\begin{matrix} {\hat{p}}_{1, 2, 1, 1} & {\hat{p}}_{1, 2, 1, 2} & {\hat{p}}_{1, 2, 1, 3} \\ {\hat{p}}_{1, 2, 2, 1} & {\hat{p}}_{1, 2, 2, 2} & {\hat{p}}_{1, 2, 2, 3} \\ {\hat{p}}_{1, 2, 3, 1} & {\hat{p}}_{1, 2, 3, 2} & {\hat{p}}_{1, 2, 3, 3} \end{matrix}] = [\begin{matrix} 0.408 & 0.238 & 0.354 \\ 0.217 & 0.296 & 0.487 \\ 0.070 & 0.165 & 0.765 \end{matrix}]

Specifically, if the first patient’s baseline state is severe (u = 1), then after one-year treatment, the patient will stay in a severe state (s = 1) with a probability of 0.408, will transition to a moderate state (s = 2) with a probability of 0.238, and will jump to a minor state (s = 3) with a probability of 0.354. Likewise, if the patient’s baseline state is moderate (u = 2) or minor (u = 3), then after one-year treatment, the patient will transition to a severe (s = 1), moderate (s = 2), and minor (s = 3) state with probabilities 0.217 or 0.070, 0.296 or 0.165, and 0.487 or 0.765, respectively.

6 Discussion

We proposed a new quantile HMM to reveal the association structure of a longitudinal process and its dynamic heterogeneity simultaneously. The proposed model accommodates nonhomogeneous transitions between hidden states by allowing the transition probabilities to depend on exogenous covariates. Unlike mean regression that focuses on modeling the central tendency of the response variable, the proposed model provides a systematic and flexible method to analyze the entire conditional distribution of the response given the covariates and hidden states. We developed a Bayesian approach coupled with efficient MCMC methods to conduct statistical inference. Simulation studies demonstrated that the proposed method performs satisfactorily regardless of the error distributions and quantile levels.

The proposed method can be extended in several directions. First, we may consider nonparametric or semiparametric models, such as varying coefficient and additive models, to formulate the conditional QR and/or nonhomogeneous transition parts. Bayesian nonparametric techniques, such as Bayesian P-splines,^47,48 can be used to model the nonparametric functions. Second, model/variable selection is an important inference beyond estimation. Bayesian penalized likelihood methods, including but not limited to Bayesian Lasso⁴⁹ and spike-and-slab regression,^50,51 can be developed to determine the number of hidden states and select important variables in quantile HMMs. Third, Sriram et al.³⁶ developed the consistency of Bayesian QR by using the misspecified ALD in the case of QR on independent responses. However, whether this result would necessarily carry across to HMMs with longitudinal responses remains unclear. While our simulation studies provided empirical evidence on consistent estimates, we failed to provide theoretical guarantees of posterior consistency in the context of the proposed quantile HMM. This theoretical problem is worthy of investigation in the future. Lastly, the proposed method does not give valid posterior standard deviations. Yang et al.³⁰ pointed out that the posterior from the working likelihood based on ALD is not the conditional distribution of the parameter given the data; thus, the credible intervals obtained from the posterior do not generally have the right Bayesian confidence level. Our simulation confirmed this result and showed that the credible intervals (not reported) over-cover with almost 100% empirical coverage probabilities at the nominal level of 95%. Yang et al.³⁰ also proposed a correction method for standard deviations under the cross-sectional (single-state) QR setting. However, their correction method is not directly applicable to the present study due to the dynamic transition feature of the proposed model. How to overcome this limitation is an interesting direction for future research. Substantial efforts are required for the aforementioned extensions.

Footnotes

Acknowledgements

The authors would like to thank two anonymous reviewers, an associate editor, and the editor for constructive comments and helpful suggestions.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The work was supported by National Natural Science Foundation of China grants (11671268, 11871376), General Research Fund grants (14303017, 14302519) from the Research Grant Council of Hong Kong Special Administration Region, Yunnan Provincial Science and Technology Department Program 2018FH001-109, and Shanghai Pujiang Program 18PJ1409800.

ORCID iD

Xinyuan Song

References

Cappé

Moulines

Rydén

Inference in hidden Markov models. New York: Springer, 2005.

Maruotti

Mixed hidden Markov models for longitudinal data: an overview. Int Stat Rev 2011; 79: 427–454.

Bartolucci

Farcomeni

Pennoni

Latent Markov models for longitudinal data. Florida: Chapman & Hall/CRC, Taylor and Francis Group, 2013.

Albert

McFarland

Smith

, et al.

Time series for modelling counts from a relapsing-remitting disease: application to modelling disease activity in multiple sclerosis.

Stat Med 1994; 13: 453–466.

Vermunt

Langeheine

Böckenholt

Mixed hidden Markov models for longitudinal data: an overview. J Educ Behav Stat 1999; 24: 178–205.

Scott

James

Sugar

CA.

Hidden Markov models for longitudinal comparisons. J Am Stat Assoc 2005; 100: 359–369.

Altman

RM.

Mixed hidden Markov models. J Am Stat Assoc 2007; 102: 201–210.

Dannemann

Holzmann

Likelihood ratio testing for hidden Markov models under non-standard conditions. Scand J Stat 2007; 25: 309–321.

Bartolucci

Farcomeni

A multivariate extension of the dynamic logit model for longitudinal data based on a latent Markov heterogeneity structure. J Am Stat Assoc 2009; 104: 816–831.

10.

Chow

Grimm

Filteau

, et al.

Regime-switching bivariate dual change score model.

Multivar Behav Res 2013; 48: 463–502.

11.

Zhang

Rejeski

, et al.

Partially ordered mixed hidden Markov model for the disablement process of older adults.

J Am Stat Assoc 2013; 108: 370–384.

12.

Song

Xia

Zhu

HT.

Hidden Markov latent variable models with multivariate longitudinal data.

Biometrics 2017; 73: 313–323.

13.

Visser

Speekenbrink

depmixS4: an R-package for hidden Markov models. J Stat Softw 2010; 36: 1–21.

14.

Song

Kang

Ouyang

, et al. Bayesian analysis of semiparametric hidden Markov models with latent variables. Struct Equ Modeling 2018; 25: 1–20.

15.

Kang

Cai

Song

, et al.

Bayesian hidden Markov models for delineating the pathology of Alzheimer’s disease.

Stat Methods Med Res 2019; 28: 2112–2124.

16.

Koenker

Bassett

JG.

Regression quantiles. Econometrica 1978; 46: 33–50.

17.

Koenker

Quantile regression. Cambridge: Cambridge University Press, 2005.

18.

Farcomeni

Quantile regression for longitudinal data based on latent Markov subject-specific parameters. Stat Comput 2012; 22: 141–152.

19.

Marino

Tzavidis

Alfo

Quantile regression for longitudinal data: unobserved heterogeneity and informative missingness. Statistics 2015; 32: 1–40.

20.

Marino

Alfo

Latent drop-out based transitions in linear quantile hidden Markov models for longitudinal responses with attrition. Adv Data Anal Class 2015; 9: 483–502.

21.

Marino

Tzavidis

Alfo

Mixed hidden Markov quantile regression models for longitudinal data with possibly incomplete sequences.

Stat Methods Med Res 2018; 27: 2231–2246.

22.

Moyeed

RA.

Bayesian quantile regression. Stat Probab Lett 2001; 54: 437–447.

23.

Kottas

Gelfand

Bayesian semiparametric median regression modeling. J Am Stat Assoc 2001; 96: 1458–1468.

24.

Dunson

Watson

Taylor

Bayesian latent variable models for median regression on multiple outcomes.

Biometrics 2003; 59: 296–304.

25.

Stander

Bayesian analysis of a Tobit quantile regression model. J Econom 2007; 137: 260–276.

26.

Lin

Bayesian regularized quantile regression. Bayesian Anal 2010; 5: 533–556.

27.

Taddy

Kottas

A Bayesian nonparametric approach to inference for quantile regression. J Bus Econ Stat 2010; 28: 357–369.

28.

Al-Hamzawi

Benoit

Bayesian adaptive lasso quantile regression. Stat Model 2012; 12: 279–297.

29.

Wang

Feng

Song

XY.

Bayesian quantile structural equation models.

Struct Equ Model 2015; 23: 246–258.

30.

Yang

Wang

Posterior inference in Bayesian quantile regression with asymmetric Laplace likelihood. Int Stat Rev 2016; 84: 327–344.

31.

Santos

Bolfarine

On Bayesian quantile regression and outliers. arXiv:1601.07344, 2016.

32.

Feng

Wang

, et al. Bayesian regularized quantile structural equation models. J Multivar Anal 2017; 154: 234–248.

33.

Benoit

Van den Poel

bayesQR: a Bayesian approach to quantile regression. J Stat Softw 2017; 76: 1–32.

34.

Koutsourelis

Bayesian extreme quantile regression for hidden Markov models. PhD Thesis, Brunel University, London, 2012.

35.

Kozumi

Kobayashi

Gibbs sampling methods for Bayesian quantile regression. J Stat Comput Simul 2011; 81: 1565–1578.

36.

Sriram

Ramamoorthi

Ghosh

Posterior consistency of Bayesian quantile regression based on the misspecified asymmetric Laplace density. Bayesian Anal 2013; 8: 479–504.

37.

Celeux

Forbes

Robert

, et al. Deviance information criteria for missing data models. Bayesian Analysis 2006; 1: 651–674.

38.

Lee

Song

XY.

Bayesian model selection for mixture of structure equation models with an unknown number of components.

Brit J Math Stat Psychol 2003; 56: 145–165.

39.

Jasra

Holmes

Stephens

Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling. Stat Sci 2005; 20: 50–67.

40.

Spezia

Reversible jump and the label switching problem in hidden Markov models. J Stat Plan Infer 2009; 139: 2305–2315.

41.

Richardson

Green

PJ.

On Bayesian analysis of mixtures with an unknown number of components (with discussion). J R Stat Soc B 1997; 59: 731–792.

42.

Boys

Henderson

A comparison of reversible jump MCMC algorithms for DNA sequence segmentation using hidden Markov models. Comput Sci Stat 2001; 33: 35–49.

43.

Papastamoulis

Iliopoulos

Reversible jump MCMC in mixtures of normal distributions with the same component means. Comput Stat Data Anal 2009; 53: 900–911.

44.

Baum

Petrie

Soules

, et al. A maximization technique occurring in the statistical analysis of probabilistic function of Markov chains. Ann Math Stat 1970; 41: 164–171.

45.

Brooks

Gelman

General methods for monitoring convergence of iterative simulations. J Comput Graph Stat 1998; 7: 434–455.

46.

Gelman

Roberts

Gilks

WR.

Efficient metropolis jumping rules. In: Bernardo

Berger

Dawid

Smith

AFM

(eds) Bayesian Statistics. Oxford: Oxford University Press, 1996, pp.599–607.

47.

Fahrmeir

Raach

A Bayesian semiparametric latent variable model for mixed responses. Psychometrika 2007; 72: 327–346.

48.

Song

ZH.

Semiparametric latent variable models with Bayesian P-splines. J Comput Graph Stat 2010; 19: 590–608.

49.

Park

Casella

The Bayesian Lasso. J Am Stat Assoc 2008; 103: 681–686.

50.

Mitchell

Beauchamp

Bayesian variable selection in linear regression. J Am Stat Assoc 1988; 83: 1023–1032.

51.

Ishwaran

Rao

Spike and slab variable selection: frequentist and Bayesian strategies. Ann Stat 2005; 33: 730–773.