There are several methods to estimate the parameters of a
hidden Markov model (HMM). To be complete, we discuss three methods
typically used when assumed that the data is generated by one
(uni-level) model. That is, the Maximum Likelihood approach, the Baum
Welch algorithm utilizing the forward backward probabilities, and the
Bayesian estimation method. Note that all these methods assume that the
number of states is known from the context of the application, i.e.,
specified by the user. The issue of determining the number of states is
discussed in the vignette “tutorial-mhmm”.
After discussing the simplified case of estimating the parameters where
the data consists of only one observed sequence (or, with multiple
sequences, assuming that all data is generated by one identical model),
we proceed with elaborating on estimating the parameters of a multilevel
hidden Markov model.
ML estimation can be used to estimate the parameters of the HMM. The relevant likelihood function has a convenient form: In equation , P(ot) denotes a diagonal matrix with the state-dependent conditional probabilities of observing Ot = o as entries, δ denotes the distribution of the initial probabilities πi, Γ denotes the transition probability matrix, and 1′ is a column vector consisting of m (i.e., the number of distinct states) elements which all have the value one. See the vignette “tutorial-mhmm” for an explanation of these quantities. Direct maximization of the log-likelihood poses no problems even for very long sequences, provided measures are taken to avoid numerical underflow.
The EM algorithm (Dempster, Laird, and Rubin 1977), in
this context also known as the Baum-Welch algorithm (Baum et al. 1970; Rabiner 1989), can also be used to
maximize the log-likelihood function. Here, the unobserved latent states
are treated as missing data, and quantities known as the forward and the
backward probabilities are used to obtain the ‘complete-data
log-likelihood’ of the HMM parameters: the log-likelihood based on both
the observed event sequence and the unobserved, or “missing”, latent
states. The forward probabilities αt(i)
denote the joint probability of the observed event sequence from time
point 1 to t and state S at time point t being i: The name “forward probabilities”
derives from the fact that when computing the forward probabilities
αt, one
evaluates the sequence of hidden states in the chronological order
(i.e., forward in time) until time point t. The backward probabilities βt(i)
denote the conditional probability of the observed event sequence after
time point t until the end, so
from t + 1, t + 2, …, T,
given that state S at time
point t equals i:
When computing the backward probabilities βt, one
evaluates the sequence of hidden states in the reversed order, i.e.,
from ST, ST − 1, …, St + 1.
The forward and backward probabilities together cover the complete event
sequence from t = 1 to T, and combined give the joint
probability of the complete event sequence and state S at time point t being i: We refer to Cappé (2005) for a
discussion on the advantages of combining forward and backward
probability information in the EM algorithm over direct maximization of
the likelihood for the HMM.
A third approach is to use Bayesian estimation to infer the
parameters of the HMM. We refer to e.g., Gelman
et al. (2014), Lynch (2007), Rossi, Allenby, and McCulloch (2012) for an in-depth exposition of
Bayesian statistics. In general terms, the difference between
frequentist and Bayesian estimation is the following. In frequentist
estimation, we view the parameters as fixed entities in the population,
which are subject only to sampling fluctuation (as quantified in the
standard error of the estimate). In Bayesian estimation, however, we
assume that each parameter follows a given distribution. The general
shape of this distribution is determined beforehand, using a prior
distribution. This prior distribution not only determines the shape of
the parameter distribution, but also allows for giving some information
to the model with respect to the most likely values of the parameter
that is estimated. To arrive at the final distribution for the parameter
- which is called the posterior distribution -, the prior distribution
is combined with the likelihood function of the data using Bayes’
theorem. Here, the likelihood function provides us with the probability
of the data given the parameters. While any aspect of these
distributions may be of interest, the emphasis is usually on the mean
(or median) of the posterior distribution, which serves as the point
estimate of the parameter of interest (analogous to the frequentist
parameter estimates). In the event that one has no or vague expectations
about the possible parameter values, one can specify “non-informative”
priors (e.g., uniform distribution). That is, one can choose parameters
of the prior distributions, so called hyper-parameters, such that the
parameters may assume a wide range of possible values. Non-informative
priors therefore express a lack of knowledge.
In the implemented hidden Markov model, both the transitions from state
i at time point t to any of the other states at time
point t + 1 and the observed
outcomes within state i follow
a categorical distribution, with parameter sets Γi
(i.e., the probabilities in row i of the transition probability
matrix Γ) and θi
(i.e., the state i dependent
probabilities of observing an act). A convenient (conjugate) prior
distribution on the parameters of the categorical distribution is a
(symmetric) Dirichlet distribution. We assume that the rows of Γ and the state-dependent
probabilities θi are
independent. That is, where the probability distribution of the current
state St
is given by the row of the transition probability matrix Γ corresponding to the
previous state in the hidden state sequence St − 1. The
probability distribution of St given by
Γ holds for states
after the first time point, i.e., t starts at 2 as there is no previous
state in the hidden state sequence for state S at t = 1. The probability of the first
state in the hidden state sequence S1 is given by the
initial probabilities of the states πi. The
probability distribution of the observed event Ot is given by
state-dependent probabilities θi
corresponding to the current state St. The
hyper-parameter a10 of the prior
Dirichlet distribution on Γi is a
vector with length equal to the number of states m, and the hyper-parameter a20 of the prior
Dirichlet distribution on θi is a
vector with length equal to the number of categorical outcomes q. Note that in this model, the
hyper-parameter values are assumed invariant over the states i. The initial probabilities of the
states πi
are assumed to coincide with the stationary distribution of Γ and are therefore not
independent (to-be-estimated) parameters.
Given these distributions, our goal is to construct the joint posterior
distribution of the hidden state sequence and the parameter estimates,
given the observed event sequence and the hyper-parameters by drawing
samples from the posterior distribution. By applying a Gibbs sampler, we
can iteratively sample from the appropriate conditional posterior
distributions of St, Γi and
θi,
given the remaining parameters in the model. In short, the Gibbs sampler
iterates between the following two steps: first the hidden state
sequence S1, S2, …, ST
is sampled, given, the observed event sequence O1, O2, …, OT,
and the current values of the parameters Γ and θi.
Subsequently, the remaining parameters in the model (Γi and
θi) are
updated by sampling them conditional on the sampled hidden state
sequence S1, S2, …, ST
and observed event sequence O1, O2, …, OT.
Sampling the hidden state sequence of the HMM by means of the Gibbs
sampler can be performed in various ways. Here, we use the approach
outlined by Scott (2002). That is, we use the
forward-backward Gibbs sampler, in which first the forward probabilities
αt(i)
(i.e., the joint probability of state S = i at time point t and the observed event sequence
from time point 1 to t) as
given in equation are obtained, after which the hidden state sequence is
sampled in a backward run (i.e., drawing ST, ST − 1, …, S1)
using the corresponding forward probabilities αT : 1.
The forward-backward Gibbs sampler produces sampled values that rapidly
represent the complete area of the posterior distribution, and produces
useful quantities as byproducts, such as the log-likelihood of the
observed data given the current draws of the parameters in each
iteration (Scott
2002). In the section “” , we provide a more detailed
description of how the Gibbs sampler proceeds for the HMM.
As it generally takes a number of iterations before the Gibbs sampler
converges to the appropriate region of the posterior distribution, the
initial iterations are usually discarded as a ‘burn-in’ period. The
remaining sampled values of Γi and
θi
provide the posterior distributions of their respective
parameters.
A problem that can arise when using Bayesian estimation in this context
is “label switching”, i.e., as the hidden states of the HMM have no a
priori ordering or interpretation, their labels (i.e., which state
represents what) can switch over the iterations of the Gibbs sampler,
without affecting the likelihood of the model (see e.g., Scott 2002; Jasra, Holmes, and Stephens 2005). As a
result, the marginal posterior distributions of the parameters are
impossible to interpret because they represent the distribution of
multiple states. Sometimes, using reasonable starting values (i.e., the
user-specified parameter values of the “zero-th” iteration used to start
the MCMC sampler) suffices to prevent label switching. Otherwise,
possible solutions are to set constraints on the parameters of the
state-dependent distribution, or use (weakly) informative priors on the
state-dependent distributions (Scott 2002). Hence, before making
inferences from the obtained marginal distributions, one should first
assess if the problem of label switching is present (e.g., by using
plots of the sampled parameter values of the state-dependent
distributions over the iterations), and if necessary, take steps to
prevent the problem of label switching. In our own experience, the use
of reasonable starting values always sufficed to prevent label
switching.
Both EM and Bayesian Gibbs sampling are viable inferential procedures
for HMMs, but for more complex HMMs such as multilevel HMMs, the
Bayesian estimation method has several advantages (e.g., lower
computational cost, and less computation time) over the EM algorithm. We
refer to Rydén (2008) for a comparison on frequentist
(i.e., the EM algorithm) and Bayesian approaches.
Bayesian estimation is particularly suited to model multilevel
models. In the multilevel model, we have a multi-layered structure in
the parameters. For the HMM, we have subject level parameters at the
first level pertaining to the observations within a subject, and group
level parameters at the second level that describe the mean and
variation within the group, as inferred from the sample of subjects. To
illustrate the multilevel model, suppose that we have K subjects for which we have each
H observations on their number
of cups of coffee consumed per day y, i.e., subject k ∈ {1, 2, …, K} and
observation h ∈ {1, 2, …, H}. Hence, at
the first level, we have daily observations on coffee consumption within
subjects: y11,
y12, …, y1H, y21, y22, …, y2H, yK1, yK2, …, yKH.
Using a multilevel model, the observations of each subject are
distributed according to the same distribution Q, but each subject has its own
parameter set θk. That is: In
addition, the subject-specific parameter sets θk are
realizations of a common group level distribution W with parameter set Λ: That is, in the multilevel model,
the subject level model parameters that pertain to the observations
within a subject are assumed to be random draws from a given
distribution, and, as such, are denoted as “random”, independent of the
used estimation method (i.e., Bayesian or classical frequentist
estimation). This multi-layered structure fits naturally into a Bayesian
paradigm since in Bayesian estimation, model parameters are by
definition viewed as random. That is, parameters follow a given
distribution, where the prior distribution expresses the prior
expectations with respect to the most likely values of the model
parameters. In the multilevel model, the prior expectations of the
subject level model parameter values are reflected in the group level
distribution. Hence, in Bayesian estimation, the prior distribution for
the subject level parameters, is given by the group level distribution.
The group level distribution provides information on the location (e.g.,
mean) of the subject level (i.e., subject-specific) parameters, and on
the variation in the subject level parameters. As the Normal
distribution is a flexible distribution with parameters that easily
relate to this interpretation, the group level distribution is often
taken to be a normal distribution.
To illustrate the notion of the group level (prior) distribution,
suppose we assume a Poisson distribution for the observations on daily
coffee consumption within each subject k, and a Normal group level
distribution on the Poisson mean. In this case, the set of
hyper-parameters (i.e., the parameters of the group level distribution,
here the mean (Λμ) and variance
(Λσ2)
of the Normal distribution) on the Poisson mean denote the group mean
number of cups of coffee consumed per day over subjects, and the
variation in the mean number of cups of coffee consumed per day between
subjects.
Finally, in fitting the multilevel model using Bayesian estimation, a
prior distribution is placed on each of these hyper-parameters. Prior
distributions on hyper-parameters are referred to as hyper-priors and
allow the hyper-parameters to have a distribution instead of being
fixed. That is, as the parameters that characterize the group level
prior distribution (i.e., the hyper-parameters) are now also quantities
of interest (i.e., to-be-estimated), they are viewed as random in
Bayesian estimation methods. The randomness in the hyper-parameters is
thus specific to the Bayesian estimation method of the multilevel model,
in contrast to the randomness in the subject level parameters.
To continue our example, the hyper-prior on the mean of the Normal prior
distribution for the subject level mean of cups of coffee consumed daily
denote our prior belief on the mean number of cups of coffee consumed
per day in the group. The hyper-prior on the variance of the Normal
prior distribution for the subject level mean of cups of coffee consumed
per day denote our prior belief on how much this mean number of cups of
coffee varies over subjects. Often, the hyper-prior distribution and its
values are chosen to be vague (i.e., not informative), like a uniform
distribution: See e.g., Gelman et al. (2014), Lynch
(2007), Rossi, Allenby, and McCulloch (2012) for an in-depth exposition of
various multilevel Bayesian models and e.g. Snijders and Bosker (2011), Hox,
Moerbeek, and Schoot (2017), Goldstein (2011) for coverage of the
classical, frequentist approach to multilevel (also called hierarchical
or random effects) models.
The present implemented multilevel model pertains only to data comprised
of (multivariate) categorical observations, and, possibly, time
invariant covariates (i.e., for each covariate, we have one value per
subject). As such, data is comprised of Od, k, t
observations on the categorical outcome(s) for categorical outcome
variable d = 1, 2, …, D for subject
k = 1, 2, …, K at
time point t = 1, 2, …, T. In
addition, we have a matrix X that consists of k covariate vectors with length
p × 1, Xk = (Xk1, Xk2, …, Xkp).
As yet, the explanation of the estimation procedure in this vignette is
restricted to the univariate case for simplicity. That is, to the
instance that we only have one observed categorical outcome variable per
subject, and our outcome data is comprised of Okt
observations. However, the explanation extents quite naturally to the
multivariate case.
Given these observations, we construct a multilevel model for each of
the parameters in the HMM with q observable categorical outcomes,
and m hidden states, possibly
predicted by p covariates.
Using the multilevel framework, each parameter is assumed to follow a
given distribution, and the parameter value of a given subject
represents a draw from this common (i.e., group level) distribution.
Hence, in the multilevel Bayesian HMM, the parameters are: the
subject-specific transition probability matrix Γk with
transition probabilities γkij
and the subject-specific state-dependent probability distribution
denoting the subject-specific probabilities θki
of categorical outcomes within state i. The initial probabilities of the
states πk, j
are not estimated as πk is assumed to
be the stationary distribution of Γk.
Subsequently, the parameter values of the subjects are assumed to be
realizations of a model component and state specific (multivariate)
Normal distribution. We discuss the multilevel model for the two
components of the HMM (Γk and
θki
separately. Table provides an overview of the used symbols in the
multilevel models related to the two components of the HMM. We use the
subscript 0 to denote values of the hyper-prior distribution
parameters.
In the standard (non-multilevel) Bayesian HMM estimation, we
specified a Dirichlet prior distribution on the state-dependent
probabilities θi. To
provide a flexible model that allows for the inclusion of random effects
and (time invariant) covariates, we follow Altman
(2007) and extend the
subject-specific state-dependent probabilities θki
to a Multinomial logit (MNL) model. Hence, we utilize a linear predictor
function to estimate the probability of observing categorical outcome
l within state i. The state i specific linear predictor function
at the subject level consists of q − 1 random intercepts (i.e., each
subject has its own intercept). That is, each categorical outcome l has its own intercept, with the
exception of the first categorical outcome in the set for which the
intercept is omitted for reasons of model identification (i.e., not all
probabilities can be estimated freely as within subject k and state i, the probabilities need to add up
to 1). By making the intercepts random (i.e., each subject has its own
intercept), we accommodate heterogeneity between subjects in their state
conditional probabilities. Hence, in the MNL model for θki,
subject k’s probabilities of
observing categorical outcome l ∈ {1, 2, …, q} within a
state i ∈ {1, 2, …, m}, θkil,
are modeled using m batches of
q − 1 random intercepts, α(O)ki = (α(O)ki2, α(O)ki3, …, α(O)kiq).
That is, where K, m, and q are the number of subjects,
states, and categorical outcomes, respectively. The numerator is set
equal to one for l = 1, making
the first categorical outcome in the set the baseline category in every
state.
At the group level, these subject-level intercepts are (possibly) partly
determined by covariates that differentiate between subjects. Thus, in
addition to the subject level random intercepts α(O)ki,
we have m matrices of p * (q − 1) fixed
regression coefficients, β(O)i,
where p denotes the number of
used covariates. The columns of β(O)i
are β(O)il = (β(O)il1, β(O)il2, …, β(O)ilp)
to model the random intercepts for state i and categorical outcome l given p covariates. Combining both terms,
each batch of random intercepts α(O)ki
(i.e., the batch of q − 1
intercepts for the state i
conditional probabilities of a categorical outcome for subject k) come from a state i specific population level
multivariate Normal distribution, with mean vector ᾱ(O)i + Xk⊤β(O)il
that has length q − 1, and
covariance Φi that denotes
the covariance between the q − 1 state i specific intercepts over subjects
and models the dependence of the probabilities of categorical outcomes
within state i (i.e., we
specify a state specific multivariate Normal prior distribution on the
subject-specific α(O)ki
parameters). A convenient hyper-prior on the hyper-parameters of the
group level prior distribution is a multivariate Normal distribution for
the mean vector ᾱ(O)i
and the fixed regression coefficients β(O)il,
and an Inverse Wishart distribution for the covariance Φi (see e.g., Gelman et al.
2014). That is, where the probability distribution of the
observed categorical outcomes Okt is
given by the subject k
specific state-dependent probabilities θki
corresponding to the current state Skt
(where t indicates the time
point, see Table ). The parameters α(O)0,
β(O)0,
and K0 denote the
values of the parameters of the hyper-prior on the group (mean) vector
ᾱ(O)i
and β(O)i,
respectively. Here, α(O)0
and β(O)0
represent a vector of means and K0 denotes the number of
observations (i.e., the number of hypothetical prior subjects) on which
the prior mean vector α(O)0
and β(O)0
are based, i.e., K0
determines the weight of the prior on ᾱ(O)i
and β(O)i.
The parameters Φ0
and df0,
respectively, denote the values of the covariance and the degrees of
freedom of the hyper-prior Inverse Wishart distribution on the
population variance Φi of the
subject-specific random intercepts α(O)ki.
Note that we chose the values of the parameters of the hyper-prior
distributions that result in uninformative hyper-prior distributions, as
such the values of the parameters of the hyper-priors are assumed
invariant over the states i.
Similar to the state-dependent probabilities θki,
we extend each set of state i
specific state transition probabilities γkij
to a MNL model to allow for the inclusion of random effects and (time
invariant) covariates. Hence, we use a linear predictor function to
estimate the probability to transition from behavioral state i to state j. The linear predictor function
consists of m − 1 random
intercepts to allow for heterogeneity between subjects in their
probabilities to switch between states. That is, within row i of the transition probability
matrix Γk,
each state j has its own
intercept, where the intercept that relates to transitioning to the
first state in the set is omitted for reasons of model identification
(i.e., not all probabilities can be estimated freely as the
row-probabilities need to add up to 1). Hence, each subject’s
probability to transition from behavioral state i ∈ {1, 2, …, m} to state
j ∈ {1, 2, …, m} is
modeled using m batches of
m − 1 random intercepts, α(S)ki = (α(S)k13, …, α(S)k1m, α(S)k23, …, α(S)k2m,
…, α(S)km2,
…, α(S)km(m − 1)).
That is, where K and m are again the number of subjects
in the dataset, and the distinct number of states, respectively. The
numerator is set equal to 1 for j = 1, making the first state of
every row of the transition probability matrix Γk the
baseline category.
At the group level, these subject-level intercepts are (possibly) partly
determined by covariates that differentiate between subjects. Thus, in
addition to the subject level random intercepts α(S)ki,
we have m matrices of p * (q − 1) fixed
regression coefficients, β(S)i,
where p denotes the number of
used covariates. The columns of β(S)i
are β(S)ij = (β(S)ij1, β(S)ij2, …, β(S)ijp)
to model the random intercepts denoting the probability of transitioning
from behavioral state i to
state j given p covariates. Combining both terms,
each batch of random intercepts α(S)ki
come from a state i specific
population level multivariate Normal distribution, with mean vector
ᾱ(S)i + Xk⊤β(S)ij
that has length q − 1, and
covariance Ψii
that denotes the covariance between the q − 1 state i specific intercepts over subjects,
and models the dependency between the probabilities of states within
random intercept batch α(S)ki
(i.e., we specify a state specific multivariate Normal prior
distribution on the subject-specific α(S)ki
parameters). A convenient hyper-prior on the hyper-parameters of the
group level prior distribution is a multivariate Normal distribution for
the mean vector ᾱ(S)i
and the fixed regression coefficients β(S)il
and an Inverse Wishart distribution for the covariance Ψi. That is,
where the subject-specific probability distribution of the current state
Skt is
given by the row of the transition probability matrix Γk
corresponding to the previous state in the hidden state sequence Sk, t − 1.
The probability distribution of Skt
given by Γk holds
for states after the first time point, i.e., t starts at 2 as there is no
previous state in the hidden state sequence for state Skt at
t = 1. The probability of the
first state in the hidden state sequence Sk, 1 is given
by the initial probabilities of the states πk, j.
The parameters α(S)0,
β(S)0,
and K0 denote the
values of the parameters of the hyper-prior on the group (mean) vector
ᾱ(S)i
and β(S)i,
respectively. Here, α(S)0
and β(S)0
represent a vector of means and K0 denotes the number of
observations (i.e., the number of hypothetical prior subjects) on which
the prior mean vector α(S)0
and β(S)0
are based. The parameters Ψ0 and df0,
respectively, denote values of the covariance and the degrees of freedom
of the hyper-prior Inverse Wishart distribution on the group variance
Ψi of the
subject-specific random intercepts α(S)ki.
Given the above distributions, our goal is to construct the joint
posterior distribution of the parameters - i.e., the subject-specific
hidden state sequences, the subject level (i.e., subject-specific)
parameters and the group level parameter estimates - given the
observations (i.e., the observed event sequences for all k subjects that are analyzed
simultaneously as one group, and the hyper-prior parameter values) by
drawing samples from the posterior distribution. We follow a MCMC
sampler algorithm to iteratively sample from the appropriate conditional
posterior distributions of α(O)ki,
α(S)ki,
ᾱ(O)i,
β(O)i,
Φi, ᾱ(S)i,
β(S)i,
and Ψi
given the remaining parameters in the model (see below). The conditional
posterior distributions of all parameters are provided in the Section
“”.
In Bayesian estimation, it is preferable to use the natural conjugate
prior as prior distribution, as this conveniently results in a closed
form expression of the (conditional) posterior distribution(s), making
Gibbs sampling possible. However, as the non-conjugate Normal prior
provides a much more intuitive interpretation of the prior group level
distribution compared to using the natural conjugate prior of the MNL
model, and since the asymptotic Normal approximation is excellent for
the MNL likelihood (Rossi, Allenby, and McCulloch 2012), we
opt for the former and do not use the conjugate prior of the MNL model.
Therefore, we cannot use a Gibbs sampler to update the parameters of the
subject-specific state-dependent distributions and the subject-specific
transition probabilities, α(O)ki
and α(S)ki,
respectively. Instead, we use a combination of the Gibbs sampler and the
Metropolis algorithm, i.e., a Hybrid Metropolis within Gibbs sampler.
That is, we use a Metropolis sampler to update α(O)ki
and α(S)ki,
and we use a Gibbs sampler to update all other model parameters. There
are various types of Metropolis algorithms, and each type involves
specific choices. Simulation studies showed that, in line with Rossi, Allenby, and McCulloch (2012), the Random Walk (RW) Metropolis
sampler outperformed the Independence Metropolis sampler in terms of
efficiency for estimating the parameters of the (multilevel) HMM, we
chose to use the RW Metropolis sampler to update the parameters of the
subject-specific state-dependent distributions (α(O)ki)
and subject-specific state transition probabilities (α(S)ki)
in our Hybrid Metropolis within Gibbs sampler.
The Hybrid Metropolis within Gibbs sampler for the multilevel HMM
proceeds in a similar fashion as the Gibbs sampler for the HMM: first
the hidden state sequences are sampled (for each subject separately),
after which the (subject level and group level) parameters are sampled
given the observed event sequence (for each subject, Okt),
the sampled hidden state sequences (of each subject, Skt),
and the current values of the remaining parameters in the model. We
provide a stepwise walkthrough of the hybrid Metropolis within Gibbs
sampler for the multilevel HMM below.
The Hybrid Metropolis within Gibbs sampler used to fit the multilevel HMM proceeds as described below. We use the subscript c to denote the current (i.e., updated using a combination of the value of the hyper-prior and the data) parameters of the conditional posterior distributions.
RW Metropolis sampler is repeated for each subject k ∈ {1, 2, …, K}. \end{itemize}
These steps are repeated for a large number of iterations, and, after
discarding the first iterations as a “burn-in” period, the sampled
parameter estimates provide the empirical posterior distribution of the
model parameters.
Regarding the acceptance rate of the RW Metropolis sampler for the
subject-specific sets of intercepts α(O)ki
and α(S)ki
(i.e., related to the subject-specific state-dependent probabilities
θki
and the subject-specific state transition probability matrix Γk,
respectively), an acceptance rate of ∼ 23% is considered optimal when many
parameters are being updated at once (Gelman et al. 2014). Within the R
package mHMMbayes, the number of accepted draws of a model are stored in
emiss_naccept and gamma_naccept for the conditional distributions and
the transition probabilities, respectively.
To obtain the scale parameter Σα(O)ki
and Σα(S)ki
of the proposal distributions of the RW Metropolis sampler for α(O)ki
and α(S)ki,
respectively, we followed the method outlined by Rossi, Allenby, and McCulloch (2012), which has several advantages as
discussed below.
The general challenge of the RW Metropolis sampler is that it has to be
“tuned” by choosing the scale of the symmetric proposal distribution
(e.g., the variance or covariance of a Normal or multivariate Normal
proposal distribution, respectively). The scale of the proposal
distribution is composed of a covariance matrix Σ, which is then tuned by
multiplying it by a scaling factor s2. Hence we denote the
scale of the proposal distribution by s2Σ. The scale
s2Σ has to
be set such that the drawn parameter estimates cover the entire area of
the posterior distribution (i.e., the scale Σ should not be set too narrow
because then only candidate parameters in close proximity of the current
parameter will be drawn), but remains reasonably efficient (i.e., the
scale Σ should not be set too
wide because then many candidate parameters will be rejected resulting
in a slowly progressing chain of drawn parameter estimates).
There are various options for the covariance matrix Σ. Often, the covariance matrix
Σ is set such that it
resembles the covariance matrix of the actual posterior distribution. To
capture the curvature of each subject’s conditional posterior
distribution, the scale of the RW Metropolis proposal distribution
should be customized to each subject. This also facilitates the
possibility to let the amount of information available within the data
of a subject for a parameter determine to which degree the group level
distribution dominates the estimation of the subject-specific
parameters. Hence, to approximate the conditional posterior distribution
of each subject, the covariance matrix is set to be a combination of the
covariance matrix obtained from the subject data and the group level
covariance matrix Φi or Ψi. To estimate
the covariance matrix from the subject data, which is only used for the
proposal distribution of the RW Metropolis sampler, we simply use a
Maximum Likelihood Estimate (MLE), as this quantity is only used for the
purpose of scaling the proposal distribution and is not part of the
estimated parameter values that constitute the posterior distribution.
The MLE estimate of the covariance matrix is obtained by maximizing the
likelihood of the Multinomial Logit (MNL) model on the data, and
retrieving the Hessian matrix Hki
(i.e., the second order partial derivatives of the likelihood function
with respect to the parameters). The covariance matrix of the parameters
is the inverse of the Hessian matrix, Hki−1.
Hence, the covariance matrices for α(O)ki
and α(S)ki,
are defined by Σα(O)ki = (Hα(O)ki + Φi−1)−1
and Σα(S)ki = (Hα(S)ki + Ψi−1)−1,
respectively. For α(O)ki,
the data on which the Hessian is obtained is the frequency with which a
categorical outcome is observed in state i of subject k. For α(S)ki,
this data is the frequency with which state i transitions to another state
within subject k. Hence, the
subject-specific covariance matrix (i.e., the inverse of the Hessian
matrix) is based on the sampled hidden state sequence. Therefore, the
MLE estimates of the subject-specific covariance matrices that are used
for the RW Metropolis proposal distributions have to be obtained in each
iteration, as the sampled hidden state sequence changes in each
iteration.
A potential problem with maximizing the log-likelihood of each subject’s
data, is that a certain state might not be sampled for a subject. To
circumvent this problem, we modify the subject likelihood function by
adding a so-called regularizing likelihood function that has a defined
maximum to the subject-level likelihood function. We maximize the
resulting pooled likelihood function in order to obtain the MLE
estimates. Here, we use the likelihood function of the combined data
over all subjects that are considered to be part of one group as the
regularizing likelihood function. The pooled likelihood function is
scaled by 1 − w × subject-level
likelihood + w × overall
likelihoodn.obsk/N.obs,
so that the overall likelihood function does not dominate the
subject-level likelihood function, and where n.obsk
is the number of data observations for subject k and N.obs is
the total number of data observations over all subjects in a
group.
Now that we defined the covariance matrix Σ for the scale of the RW Metropolis
sampler proposal distribution, we have to define the scalar factor s2 to obtain the scale
s2Σ of the
proposal distribution. As in Rossi, Allenby, and
McCulloch (2012)}, we adopt the
scaling proposal of Roberts, Rosenthal, et al.
(2001), and set scaling to $s^2 = (2.93 / \sqrt{\text{n.param}})^2$,
where n.param
is the number of parameters to be estimated for α(S)ki
or α(O)ki
in the RW Metropolis sampler, which equals m − 1 in case of α(S)ki
(where m denotes the number of
states) and q − 1 in case of
α(O)ki
(where q denotes the number of
categorical outcomes).
In summary, the scale parameter s2Σα(O)ki
and s2Σα(S)ki
of the proposal distributions of the RW Metropolis sampler for α(O)ki
and α(S)ki
are defined as: where Hα(O)ki
is the Hessian of the kth
subject’s data of the frequency with which a categorical outcome is
observed within state i
evaluated at the MLE of the pooled likelihood, Hα(S)ki
is the Hessian of the kth
subject’s data of the frequency with which state i transitions to another state
evaluated at the MLE of the pooled likelihood, and Φi−1
and Ψi−1
are the inverses of the group level covariance matrices. This provides
us with m pairs of scale
parameters that closely resemble the scale of the subject-level
conditional posterior distribution, and that 1) are automatically tuned
(i.e., we do not require experimentation to determine s2 to tune the covariance
matrix), 2) allow the amount of information available within the data of
a specific subject to determine the degree to which the group level
distribution dominates the estimation of that subject’s level
parameters, and 3) do not require each state to be sampled in the hidden
state sequence as not each subject-level likelihood is required to have
a maximum.
For the random intercepts α(O)ki
and α(S)ki,
related to the subject-specific state-dependent probabilities of
observing a categorical outcome θki
and the subject-specific state transition probability matrix Γk,
respectively, the choice of prior distributions does not result in
closed form expressions of the full conditional posterior distributions.
That is, for the subject-specific sets of intercepts α(O)ki
related to the subject-specific state-dependent probabilities of
observing a categorical outcome within state i, the full conditional posterior
distribution when we assess a standard multivariate normal prior is: and
the likelihood is the product of the probabilities of the observed
outcomes Okt = l ∈ {1, 2, …, q}
within sampled states S = i in subject k over time points t: where the product is restricted
to the set of time points that coincide with the sampled state S for subject k at time point t being i, and q is the number of categorical
outcomes. The numerator is set equal to one for l = 1, making the first categorical
outcome in the set the baseline category in every state.
For the subject-specific sets of intercepts α(S)ki
related to the state-transition probabilities to transition from state
i to any of the other states
j ∈ {1, 2, …, m}, the full conditional posterior
distribution when we assess a standard multivariate normal prior is: and
the likelihood is the product of the probabilities of the observed
transitions from state i in
the previous time point t − 1
to any of the other states Skt = j
over time points t in subject
k:
where the product is restricted to the set of time points that coincide
with the sampled state S in
the previous time point t − 1
being i for subject k, and m is the number of states. The
numerator is set equal to 1 for j = 1, making the first state of
every row of the transition probability matrix Γk the
baseline category.
As the conditional posterior distributions for α(O)ki
and α(S)ki
do not result in a closed form expression of a know distribution, we
cannot directly sample values of α(O)ki
and α(S)ki
from their conditional posterior distributions with pre-defined
equations on how to obtain the current (i.e., updated using a
combination of the value of the hyper-prior and the data) parameters of
the conditional posterior distributions. Instead, new values for α(O)ki
and α(S)ki
are sampled using a RW Metropolis sampler, as described above.
}