pd_RW_gamma {mHMMbayes} | R Documentation |
Proposal distribution settings RW Metropolis sampler for mHMM transition probability matrix gamma
Description
pd_RW_gamma
provides a framework to manually specify the settings of
the proposal distribution of the random walk (RW) Metropolis sampler of the
transition probability matrix gamma of the multilevel hidden Markov model,
and creates on object of the class mHMM_pdRW_gamma
. The RW metropolis
sampler is used for sampling the subject level parameter estimates relating
to the transition probability matrix gamma, that is, the Multinomial logistic
regression intercepts.
Usage
pd_RW_gamma(m, gamma_int_mle0, gamma_scalar, gamma_w)
Arguments
m |
Numeric vector with length 1 denoting the number of hidden states. |
gamma_int_mle0 |
A matrix with |
gamma_scalar |
A numeric vector with length 1 denoting the scale factor
|
gamma_w |
A numeric vector with length 1 denoting the weight for the overall log likelihood (i.e., log likelihood based on the pooled data over all subjects) in the fractional likelihood. |
Details
When no manual values for the settings of the proposal distribution of the
random walk (RW) Metropolis sampler are specified at all (that is, the
function pd_RW_gamma
is not used), all elements in
gamma_int_mle0
set to 0, gamma_scalar
set to 2.93 /
sqrt(m
- 1), and gamma_w
set to 0.1. See the section
Scaling the proposal distribution of the RW Metropolis sampler in
vignette("estimation-mhmm")
for details.
Within the function mHMM
, the acceptance rate of the RW metropolis
sampler relating to the transition probability matrix gamma can be tracked
using the output parameter gamma_naccept
. An acceptance rate of about
23% is considered optimal when many parameters are being updated at once
(Gelman, Carlin, Stern & Rubin, 2014).
Value
pd_RW_gamma
returns an object of class
mHMM_pdRW_gamma
, containing settings of the proposal distribution of
the random walk (RW) Metropolis sampler on the transition probability
matrix gamma of the multilevel hidden Markov model. The object is
specifically created and formatted for use by the function mHMM
, and
checked for correct input dimensions. The object contains the following
components:
m
Numeric vector denoting the number of hidden states, used for checking equivalent general model properties specified under
pd_RW_gamma
andmHMM
.gamma_int_mle0
A matrix containing the starting values for the maximum likelihood (ML) estimates of the Multinomial logit regression intercepts of the transition probability matrix gamma.
gamma_scalar
A numeric vector with length 1 denoting the scale factor
s
of the proposal distribution.gamma_w
A numeric vector with length 1 denoting denoting the weight for the overall log likelihood in the fractional likelihood.
References
Gelman A, Carlin JB, Stern HS, Rubin DB (2014). Bayesian Data Analysis vol. 2. Taylor & Francis.
Rossi PE, Allenby GM, McCulloch R (2012). Bayesian statistics and marketing. John Wiley & Sons.
Examples
###### Example using package example data, see ?nonverbal
# specifying general model properties:
m <- 3
# specifying manual values for RW metropolis sampler on gamma
gamma_int_mle0 <- matrix(c( -2, -2,
2, 0,
0, 3), byrow = TRUE, nrow = m, ncol = m - 1)
gamma_scalar <- c(2)
gamma_w <- c(0.2)
manual_gamma_sampler <- pd_RW_gamma(m = m, gamma_int_mle0 = gamma_int_mle0,
gamma_scalar = gamma_scalar,
gamma_w = gamma_w)
# specifying starting values
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)
start_TM <- diag(.7, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .1
start_EM <- list(matrix(c(0.05, 0.90, 0.05,
0.90, 0.05, 0.05,
0.55, 0.45, 0.05), byrow = TRUE,
nrow = m, ncol = q_emiss[1]), # vocalizing patient
matrix(c(0.1, 0.9,
0.1, 0.9,
0.1, 0.9), byrow = TRUE, nrow = m,
ncol = q_emiss[2]), # looking patient
matrix(c(0.90, 0.05, 0.05,
0.05, 0.90, 0.05,
0.55, 0.45, 0.05), byrow = TRUE,
nrow = m, ncol = q_emiss[3]), # vocalizing therapist
matrix(c(0.1, 0.9,
0.1, 0.9,
0.1, 0.9), byrow = TRUE, nrow = m,
ncol = q_emiss[4])) # looking therapist
# Note that for reasons of running time, J is set at a ridiculous low value.
# One would typically use a number of iterations J of at least 1000,
# and a burn_in of 200.
out_3st_RWgamma <- mHMM(s_data = nonverbal,
gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
start_val = c(list(start_TM), start_EM),
gamma_sampler = manual_gamma_sampler,
mcmc = list(J = 11, burn_in = 5))
out_3st_RWgamma
summary(out_3st_RWgamma)
# checking acceptance rate (for illustrative purposes, in the example,
# J is too low for getting a fair indication)
J_it <- 11 - 1 # accept/reject starts at iteration 2 of MCMC algorithm
out_3st_RWgamma$gamma_naccept / J_it
# average acceptance rate over all subjects per parameter
apply(out_3st_RWgamma$gamma_naccept / J_it, 2, mean)