mHMMbayes: Multilevel Hidden Markov Models Using Bayesian Estimation

Description

An implementation of the multilevel (also known as mixed or random effects) hidden Markov model using Bayesian estimation in R. The multilevel hidden Markov model (HMM) is a generalization of the well-known hidden Markov model, for the latter see Rabiner (1989) doi:10.1109/5.18626. The multilevel HMM is tailored to accommodate (intense) longitudinal data of multiple individuals simultaneously, see e.g., de Haan-Rietdijk et al. doi:10.1080/00273171.2017.1370364. Using a multilevel framework, we allow for heterogeneity in the model parameters (transition probability matrix and conditional distribution), while estimating one overall HMM. The model can be fitted on multivariate data with either a categorical, normal, or Poisson distribution, and include individual level covariates (allowing for e.g., group comparisons on model parameters). Parameters are estimated using Bayesian estimation utilizing the forward-backward recursion within a hybrid Metropolis within Gibbs sampler. Missing data (NA) in the dependent variables is accommodated assuming MAR. The package also includes various visualization options, a function to simulate data, and a function to obtain the most likely hidden state sequence for each individual using the Viterbi algorithm.

Author(s)

Maintainer: Emmeke Aarts e.aarts@uu.nl

Authors:

• Sebastian Mildiner Moraga

See Also

Useful links:

• https://CRAN.R-project.org/package=mHMMbayes

• Report bugs at https://github.com/emmekeaarts/mHMMbayes/issues

Transforming a set of Multinomial logit regression intercepts to probabilities

Description

int_to_prob transforms a set of Multinomial logit regression intercepts to the corresponding state transition or categorical emission observation probabilities. Note that the first state or category is assumed to be the reference category, hence no intercept is to specified for the first state or category.

Usage

int_to_prob(int_matrix)

Arguments

Details

Designed to ease the specification of informative hyper-prior values for the mean intercepts of the transition probability matrix gamma and categorical emission distribution(s) of the multilevel hidden Markov model through the functions prior_gamma and prior_emiss_cat. No check is performed on correct specifications of the dimensions.

Value

int_to_prob returns a matrix containing probabilities with each row summing to one, with the number of columns equal to the number of states / categories and the number of rows equal to the number rows specified in the input matrix.

See Also

prob_to_int for transforming a set of probabilities to a set of Multinomial logit regression intercepts, prior_gamma and prior_emiss_cat for specifying informative hyper-priors for the the multilevel hidden Markov model and mHMM to fit a multilevel hidden Markov model.

Examples

# example for transition probability matrix gamma with 3 states
m <- 3
gamma_int <- matrix(c(-1, -1,
                       3,  0,
                       0,  2), ncol = m-1, nrow = m, byrow = TRUE)
gamma_prob <- int_to_prob(gamma_int)
gamma_prob

Multilevel hidden Markov model using Bayesian estimation

Description

mHMM fits a multilevel (also known as mixed or random effects) hidden Markov model (HMM) to intense longitudinal data with categorical, continuous (i.e., normally distributed), or count (i.e., Poisson distributed) observations of multiple subjects using Bayesian estimation, and creates an object of class mHMM. By using a multilevel framework, we allow for heterogeneity in the model parameters between subjects, while estimating one overall HMM. The function includes the possibility to add covariates at level 2 (i.e., at the subject level) and have varying observation lengths over subjects. For a short description of the package see mHMMbayes. See vignette("tutorial-mhmm") for an introduction to multilevel hidden Markov models and the package, and see vignette("estimation-mhmm") for an overview of the used estimation algorithms.

Usage

mHMM(
  s_data,
  data_distr = "categorical",
  gen,
  xx = NULL,
  start_val,
  mcmc,
  return_path = FALSE,
  show_progress = TRUE,
  gamma_hyp_prior = NULL,
  emiss_hyp_prior = NULL,
  gamma_sampler = NULL,
  emiss_sampler = NULL
)

Arguments

Details

Covariates specified in xx can either be dichotomous or continuous variables. Dichotomous variables have to be coded as 0/1 variables. Categorical or factor variables can as yet not be used as predictor covariates. The user can however break up the categorical variable in multiple dummy variables (i.e., dichotomous variables), which can be used simultaneously in the analysis. Continuous predictors are automatically centered. That is, the mean value of the covariate is subtracted from all values of the covariate such that the new mean equals zero. This is done such that the presented probabilities in the output (i.e., for the population transition probability matrix and population emission probabilities) correspond to the predicted probabilities at the average value of the covariate(s).

Value

mHMM returns an object of class mHMM, which has print and summary methods to see the results. The object contains always the following components:

PD_subj

A list containing one list per subject with the elements trans_prob, cat_emiss, cont_emiss, or count_emiss in case of categorical, continuous, and count observations, respectively, and log_likl, providing the subject parameter estimates over the iterations of the MCMC sampler. trans_prob relates to the transition probabilities gamma, cat_emiss to the categorical emission distribution (emission probabilities), cont_emiss to the continuous emission distributions (subsequently the the emission means and the (fixed over subjects) emission standard deviation), count_emiss to the count emission distribution means in the natural (real-positive) scale, and log_likl to the log likelihood over the MCMC iterations. Iterations are contained in the rows, the parameters in the columns.

gamma_prob_bar

A matrix containing the group level parameter estimates of the transition probabilities over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows, and the columns contain the group level parameter estimates. If covariates were included in the analysis, the group level probabilities represent the predicted probability given that the covariate is at the average value for continuous covariates, or given that the covariate equals zero for dichotomous covariates.

gamma_int_bar

A matrix containing the group level intercepts of the Multinomial logistic regression modeling the transition probabilities over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows, and the columns contain the group level intercepts.

gamma_cov_bar

A matrix containing the group level regression coefficients of the Multinomial logistic regression predicting the transition probabilities over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows, and the columns contain the group level regression coefficients.

gamma_V_int_bar

A matrix containing the variance components for the subject-level intercepts (between subject variances) of the multinomial logistic regression modeling the transition probabilities over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows, and the columns contain the variance components for the subject level intercepts. Note that only the intercept variances (and not the co-variances) are returned.

gamma_int_subj

A list containing one matrix per subject denoting the subject level intercepts of the Multinomial logistic regression modeling the transition probabilities over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows, and the columns contain the subject level intercepts.

gamma_naccept

A matrix containing the number of accepted draws at the subject level RW Metropolis step for each set of parameters of the transition probabilities. The subjects are contained in the rows, and the columns contain the sets of parameters.

input

Overview of used input specifications: the distribution type of the observations data_distr, the number of states m, the number of used dependent variables n_dep, (in case of categorical observations) the number of output categories for each of the dependent variables q_emiss, the number of iterations J and the specified burn in period burn_in of the hybrid Metropolis within Gibbs sampler, the number of subjects n_subj, the observation length for each subject n_vary, and the column names of the dependent variables dep_labels.

sample_path

A list containing one matrix per subject with the sampled hidden state sequence over the hybrid Metropolis within Gibbs sampler. The time points of the dataset are contained in the rows, and the sampled paths over the iterations are contained in the columns. Only returned if return_path = TRUE.

Additionally, in case of categorical observations, the mHMM return object contains:

emiss_prob_bar

A list containing one matrix per dependent variable, denoting the group level emission probabilities of each dependent variable over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level emission probabilities. If covariates were included in the analysis, the group level probabilities represent the predicted probability given that the covariate is at the average value for continuous covariates, or given that the covariate equals zero for dichotomous covariates.

emiss_int_bar

A list containing one matrix per dependent variable, denoting the group level intercepts of each dependent variable of the Multinomial logistic regression modeling the probabilities of the emission distribution over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level intercepts.

emiss_cov_bar

A list containing one matrix per dependent variable, denoting the group level regression coefficients of the Multinomial logistic regression predicting the emission probabilities within each of the dependent variables over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level regression coefficients.

emiss_V_int_bar

A list containing one matrix per dependent variable, denoting the variance components for the subject-level intercepts (between subject variances) of the multinomial logistic regression modeling the categorical emission probabilities over the iterations of the hybrid Metropolis within Gibbs sampler. The iterations of the sampler are contained in the rows, and the columns contain the variance components for the subject level intercepts. Note that only the intercept variances (and not the co-variances) are returned.

emiss_int_subj

A list containing one list per subject denoting the subject level intercepts of each dependent variable of the Multinomial logistic regression modeling the probabilities of the emission distribution over the iterations of the hybrid Metropolis within Gibbs sampler. Each lower level list contains one matrix per dependent variable, in which iterations of the sampler are contained in the rows, and the columns contain the subject level intercepts.

emiss_naccept

A list containing one matrix per dependent variable with the number of accepted draws at the subject level RW Metropolis step for each set of parameters of the emission distribution. The subjects are contained in the rows, and the columns of the matrix contain the sets of parameters.

In case of continuous observations, the mHMM return object contains:

emiss_mu_bar

A list containing one matrix per dependent variable, denoting the group level means of the Normal emission distribution of each dependent variable over the iterations of the Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level emission means. If covariates were included in the analysis, the group level means represent the predicted mean given that the covariate is at the average value for continuous covariates, or given that the covariate equals zero for dichotomous covariates.

emiss_varmu_bar

A list containing one matrix per dependent variable, denoting the variance between the subject level means of the Normal emission distributions. over the iterations of the Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level variance in the mean.

emiss_sd_bar

A list containing one matrix per dependent variable, denoting the (fixed over subjects) standard deviation of the Normal emission distributions over the iterations of the Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level emission variances.

emiss_cov_bar

A list containing one matrix per dependent variable, denoting the group level regression coefficients predicting the emission means within each of the dependent variables over the iterations of the Gibbs sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group level regression coefficients.

emiss_naccept

A list containing one matrix per dependent variable with the number of accepted draws at the subject level RW Metropolis step for each set of parameters of the emission distribution. The subjects are contained in the rows, and the columns of the matrix contain the sets of parameters.

input

Overview of used input specifications: the number of states m, the number of used dependent variables n_dep, the number of iterations J and the specified burn in period burn_in of the hybrid Metropolis within Gibbs sampler, the number of subjects n_subj, the observation length for each subject n_vary, and the column names of the dependent variables dep_labels.

In case of count observations, the mHMM return object contains:

emiss_mu_bar

A list containing one matrix per dependent variable, denoting the group-level means of the Poisson emission distribution of each dependent variable over the MCMC iterations. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group-level Poisson emission means in the positive scale. If covariates were included in the analysis, the group-level means represent the predicted mean counts given that the covariate is at the average value for continuous covariates, or given that the covariate equals zero for dichotomous covariates.

emiss_varmu_bar

A list containing one matrix per dependent variable, denoting the variance between the subject-level means of the Poisson emission distribution(s) over the iterations of the MCMC sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group-level variance(s) between subject-specific means (in the positive scale).

emiss_logmu_bar

A list containing one matrix per dependent variable, denoting the group-level logmeans of the Poisson emission distribution of each dependent variable over the MCMC iterations. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group-level Poisson emission logmeans in the logarithmic scale. If covariates were included in the analysis, the group-level means represent the predicted logmean counts given that the covariate is at the average value for continuous covariates, or given that the covariate equals zero for dichotomous covariates.

emiss_logvarmu_bar

A list containing one matrix per dependent variable, denoting the logvariance between the subject-level logmeans of the Poisson emission distribution(s) over the iterations of the MCMC sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group-level logvariance(s) between subject-specific means (in the logarithmic scale).

emiss_cov_bar

A list containing one matrix per dependent variable, denoting the group level regression coefficients predicting the emission means in the logarithmic scale within each of the dependent variables over the iterations of the MCMC sampler. The iterations of the sampler are contained in the rows of the matrix, and the columns contain the group-level regression coefficients in the logarithmic scale.

input

Overview of used input specifications: the number of states m, the number of used dependent variables n_dep, the number of iterations J and the specified burn in period burn_in of the hybrid Metropolis within Gibbs sampler, the number of subjects n_subj, the observation length for each subject n_vary, and the column names of the dependent variables dep_labels.

References

Rabiner LR (1989). “A tutorial on hidden Markov models and selected applications in speech recognition.” Proceedings of the IEEE, 77(2), 257–286.

Scott SL (2002). “Bayesian methods for hidden Markov models: Recursive computing in the 21st century.” Journal of the American Statistical Association, 97(457), 337–351.

Altman RM (2007). “Mixed hidden Markov models: an extension of the hidden Markov model to the longitudinal data setting.” Journal of the American Statistical Association, 102(477), 201–210.

Rossi PE, Allenby GM, McCulloch R (2012). Bayesian statistics and marketing. John Wiley & Sons.

Zucchini W, MacDonald IL, Langrock R (2017). Hidden Markov models for time series: an introduction using R. Chapman and Hall/CRC.

See Also

sim_mHMM for simulating multilevel hidden Markov data, vit_mHMM for obtaining the most likely hidden state sequence for each subject using the Viterbi algorithm, obtain_gamma and obtain_emiss for obtaining the transition or emission distribution probabilities of a fitted model at the group or subject level, and plot.mHMM for plotting the posterior densities of a fitted model.

Examples

###### Example on package (categorical) example data, see ?nonverbal

# specifying general model properties:
m <- 2
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
start_TM <- diag(.8, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
start_EM <- list(matrix(c(0.05, 0.90, 0.05,
                          0.90, 0.05, 0.05), byrow = TRUE,
                        nrow = m, ncol = q_emiss[1]), # vocalizing patient
                 matrix(c(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), byrow = TRUE,
                        nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                 matrix(c(0.1, 0.9,
                          0.1, 0.9), byrow = TRUE, nrow = m,
                        ncol = q_emiss[4])) # looking therapist

# Run a model without covariate(s):
# 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_2st <- mHMM(s_data = nonverbal,
                gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                start_val = c(list(start_TM), start_EM),
                mcmc = list(J = 11, burn_in = 5))

out_2st
summary(out_2st)

# plot the posterior densities for the transition and emission probabilities
plot(out_2st, component = "gamma", col =c("darkslategray3", "goldenrod"))

# Run a model including a covariate (see ?nonverbal_cov) to predict the
# emission distribution for each of the 4 dependent variables:

n_subj <- 10
xx_emiss <- rep(list(matrix(c(rep(1, n_subj),nonverbal_cov$std_CDI_change),
                            ncol = 2, nrow = n_subj)), n_dep)
xx <- c(list(matrix(1, ncol = 1, nrow = n_subj)), xx_emiss)
out_2st_c <- mHMM(s_data = nonverbal, xx = xx,
                  gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                  start_val = c(list(start_TM), start_EM),
                  mcmc = list(J = 11, burn_in = 5))


###### Example on categorical simulated data
# Simulate data for 10 subjects with each 100 observations:
n_t <- 100
n <- 10
m <- 2
n_dep <- 1
q_emiss <- 3
gamma <- matrix(c(0.8, 0.2,
                  0.3, 0.7), ncol = m, byrow = TRUE)
emiss_distr <- list(matrix(c(0.5, 0.5, 0.0,
                        0.1, 0.1, 0.8), nrow = m, ncol = q_emiss, byrow = TRUE))
data1 <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                  gamma = gamma, emiss_distr = emiss_distr, var_gamma = .5, var_emiss = .5)

# Specify remaining required analysis input (for the example, we use simulation
# input as starting values):
n_dep <- 1
q_emiss <- 3

# Run the model on the simulated data:
out_2st_sim <- mHMM(s_data = data1$obs,
                 gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                 start_val = c(list(gamma), emiss_distr),
                 mcmc = list(J = 11, burn_in = 5))

###### Example on continuous simulated data
 # simulating multivariate continuous data
n_t     <- 100
n       <- 10
m       <- 3
n_dep   <- 2

gamma   <- matrix(c(0.8, 0.1, 0.1,
                    0.2, 0.7, 0.1,
                    0.2, 0.2, 0.6), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c( 50, 10,
                              100, 10,
                              150, 10), nrow = m, byrow = TRUE),
                    matrix(c(5, 2,
                             10, 5,
                             20, 3), nrow = m, byrow = TRUE))

data_cont <- sim_mHMM(n_t = n_t, n = n, data_distr = 'continuous',
                      gen = list(m = m, n_dep = n_dep),
                      gamma = gamma, emiss_distr = emiss_distr,
                      var_gamma = .1, var_emiss = c(5^2, 0.2^2))

# Specify hyper-prior for the continuous emission distribution
manual_prior_emiss <- prior_emiss_cont(
                        gen = list(m = m, n_dep = n_dep),
                        emiss_mu0 = list(matrix(c(30, 70, 170), nrow = 1),
                                         matrix(c(7, 8, 18), nrow = 1)),
                        emiss_K0 = list(1, 1),
                        emiss_V =  list(rep(5^2, m), rep(0.5^2, m)),
                        emiss_nu = list(1, 1),
                        emiss_a0 = list(rep(1.5, m), rep(1, m)),
                        emiss_b0 = list(rep(20, m), rep(4, m)))

# Run the model on the simulated data:
# 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_cont_sim <- mHMM(s_data = data_cont$obs,
                         data_distr = 'continuous',
                         gen = list(m = m, n_dep = n_dep),
                         start_val = c(list(gamma), emiss_distr),
                         emiss_hyp_prior = manual_prior_emiss,
                         mcmc = list(J = 11, burn_in = 5))




###### Example on multivariate count data
# Simulate data with one covariate for each count dependent variable

n_t     <- 200     # Number of observations on the dependent variable
m       <- 3        # Number of hidden states
n_dep   <- 2        # Number of dependent variables
n_subj  <- 30        # Number of subjects

gamma   <- matrix(c(0.9, 0.05, 0.05,
                    0.2, 0.7, 0.1,
                    0.2,0.3, 0.5), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c(20,
                             10,
                             5), nrow = m, byrow = TRUE),
                    matrix(c(50,
                             3,
                             20), nrow = m, byrow = TRUE))

# Define between subject variance to use on the simulating function:
# here, the variance is varied over states within the dependent variable.
var_emiss <- list(matrix(c(5.0, 3.0, 1.5), nrow = m),
                  matrix(c(5.0, 5.0, 5.0), nrow = m))

# Simulate count data:
data_count <- sim_mHMM(n_t = n_t,
                       n = n_subj,
                       data_distr = "count",
                       gen = list(m = m, n_dep = n_dep),
                       gamma = gamma,
                       emiss_distr = emiss_distr,
                       var_gamma = 0.1,
                       var_emiss = var_emiss,
                       return_ind_par = TRUE)


# Transition probabilities
start_gamma <- diag(0.8, m)
start_gamma[lower.tri(start_gamma) | upper.tri(start_gamma)] <-
                                            (1 - diag(start_gamma)) / (m - 1)

# Emission distribution
start_emiss <- list(matrix(c(20,10, 5), nrow = m, byrow = TRUE),
                    matrix(c(50, 3,20), nrow = m, byrow = TRUE))

# Specify hyper-prior for the count emission distribution
manual_prior_emiss <- prior_emiss_count(
  gen = list(m = m, n_dep = n_dep),
  emiss_mu0 = list(matrix(c(20, 10, 5), byrow = TRUE, ncol = m),
                   matrix(c(50, 3, 20), byrow = TRUE, ncol = m)),
  emiss_K0  = rep(list(0.1),n_dep),
  emiss_nu  = rep(list(0.1),n_dep),
  emiss_V   = rep(list(rep(10, m)),n_dep)
)

# Run model
# 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_count_sim <- mHMM(s_data = data_count$obs,
                          data_distr = 'count',
                          gen = list(m = m, n_dep = n_dep),
                          start_val = c(list(start_gamma), start_emiss),
                          emiss_hyp_prior = manual_prior_emiss,
                          mcmc = list(J = 11, burn_in = 5),
                          show_progress = TRUE)

summary(out_3st_count_sim)

Nonverbal communication of patients and therapist

Description

A dataset containing the nonverbal communication of 10 patient-therapist couples, recorded for 15 minutes at a frequency of 1 observation per second (= 900 observations per couple).

Usage

nonverbal

Format

A matrix with 10 * 900 rows and 5 variables:

id

id variable of patient - therapist couple to distinguish which observation belongs to which couple

p_verbalizing

verbalizing behavior of the patient, consisting of 1 = not verbalizing, 2 = verbalizing, 3 = back channeling

p_looking

looking behavior of the patient, consisting of 1 = not looking at therapist, 2 = looking at therapist

t_verbalizing

verbalizing behavior of the therapist, consisting of 1 = not verbalizing, 2 = verbalizing, 3 = back channeling

t_looking

looking behavior of the therapist, consisting of 1 = not looking at patient, 2 = looking at patient

Predictors of nonverbal communication

Description

A dataset containing predictors of nonverbal communication of 10 patient-therapist couples.

Usage

nonverbal_cov

Format

A matrix with 10 rows and 3 variables:

diagnosis

Diagnosis of the patient, consisting of 0 = depression, 1 = anxiety

std_CDI_change

Change in measure for depression (CDI) before and after therapy, standardized scale

std_SCA_change

Change in measure for anxiety (SCARED) before and after therapy, standardized scale

Obtain the emission distribution probabilities for a fitted multilevel HMM

Description

obtain_emiss obtains the emission distribution for an object containing a fitted multilevel hidden Markov model, either at the group level, i.e., representing the average emission distribution over all subjects, or at the subject level, returning the emission distribution for each subject.

Usage

obtain_emiss(object, level = "group", burn_in = NULL)

Arguments

Value

obtain_emiss creates an object of the class mHMM_emiss. Depending on the specification at the input variable level, the output is either a list of matrices with the emission distribution at the group level (if level = "group") for each dependent variable, or a list of lists, where for each dependent variable a list is returned with the number of elements equal to the number of subjects analyzed (if level = 'subject'). In the latter scenario, each matrix in the lower level list represents the subject specific emission distribution for a specific dependent variable.

See Also

mHMM for fitting the multilevel hidden Markov model.

Examples

###### Example on package data, see ?nonverbal

# specifying general model properties:
m <- 2
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
start_TM <- diag(.8, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
start_EM <- list(matrix(c(0.05, 0.90, 0.05,
                          0.90, 0.05, 0.05), byrow = TRUE,
                        nrow = m, ncol = q_emiss[1]), # vocalizing patient
                 matrix(c(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), byrow = TRUE,
                        nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                 matrix(c(0.1, 0.9,
                          0.1, 0.9), byrow = TRUE, nrow = m,
                        ncol = q_emiss[4])) # looking therapist

# Run a model without covariate(s):
out_2st <- mHMM(s_data = nonverbal,
                gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                start_val = c(list(start_TM), start_EM),
                mcmc = list(J = 11, burn_in = 5))

out_2st
summary(out_2st)

# obtaining the emission probabilities at the group and subject level
obtain_emiss(out_2st, level = "group")
obtain_emiss(out_2st, level = "subject")

Obtain the transition probabilities gamma for a fitted multilevel HMM

Description

obtain_gamma obtains the transition probability matrix (TPM) for an object containing a fitted multilevel hidden Markov model. The TPM can be obtained either at the group level, i.e., representing the average transition probability matrix over all subjects, or at the subject level, returning the transition probability matrices for each subject.

Usage

obtain_gamma(object, level = "group", burn_in = NULL)

Arguments

Value

obtain_gamma creates an object of the class mHMM_gamma. This object can be directly plotted using the function plot.mHMM_gamma(), or simply plot(). Depending on the specification at the input variable level, the output is either a matrix with the transition probabilities at the group level (if level = "group"), or a list of matrices (with the number of elements equal to the number of subjects analyzed, if level = 'subject'), where each matrix in the list represents a subject specific transition probability matrix.

See Also

mHMM for fitting the multilevel hidden Markov model, and plot.mHMM_gamma for plotting the obtained transition probabilities.

Examples

###### Example on package data

# specifying general model properties:
m <- 2
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
start_TM <- diag(.8, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
start_EM <- list(matrix(c(0.05, 0.90, 0.05,
                          0.90, 0.05, 0.05), byrow = TRUE,
                        nrow = m, ncol = q_emiss[1]), # vocalizing patient
                 matrix(c(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), byrow = TRUE,
                        nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                 matrix(c(0.1, 0.9,
                          0.1, 0.9), byrow = TRUE, nrow = m,
                        ncol = q_emiss[4])) # looking therapist

# Run a model without covariate(s):
out_2st <- mHMM(s_data = nonverbal,
                gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                start_val = c(list(start_TM), start_EM),
                mcmc = list(J = 11, burn_in = 5))

out_2st
summary(out_2st)

# obtaining the transition probabilities at the group and subject level
obtain_gamma(out_2st, level = "group")
obtain_gamma(out_2st, level = "subject")

Proposal distribution settings RW Metropolis sampler for mHMM categorical emission distribution(s)

Description

pd_RW_emiss_cat provides a framework to manually specify the settings of the proposal distribution of the random walk (RW) Metropolis sampler of categorical emission distribution(s) of the multilevel hidden Markov model, and creates on object of the class mHMM_pdRW_emiss. The RW metropolis sampler is used for sampling the subject level parameter estimates relating to the emission distributions of the dependent variables k, that is, the Multinomial logistic regression intercepts.

Usage

pd_RW_emiss_cat(gen, emiss_int_mle0, emiss_scalar, emiss_w)

Arguments

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_emiss_cat is not used), all elements in emiss_int_mle0 set to 0, emiss_scalar set to 2.93 / sqrt(q_emiss[k] - 1), and emiss_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 emission distribution(s) can be tracked using the output parameter emiss_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_emiss_cat returns an object of class mHMM_pdRW_emiss, containing settings of the proposal distribution of the random walk (RW) Metropolis sampler on the categorical emission distribution(s) 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:

gen

A list containing the elements m, n_dep, and q_emiss, used for checking equivalent general model properties specified under pd_RW_emiss_cat and mHMM.

emiss_int_mle0

A list containing n_dep elements, where each element is a matrix containing the starting values for the maximum likelihood (ML) estimates of the Multinomial logit regression intercepts of the emission distribution(s).

emiss_scalar

A list containing n_dep elements denoting the scale factor s of the proposal distribution.

emiss_w

A list containing n_dep elements 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
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying manual values for RW metropolis sampler on emission distributions
emiss_int_mle0 <- list(matrix(c( 2,  0,
                                -2, -2,
                                 0, -1), byrow = TRUE, nrow = m, ncol = q_emiss[1] - 1),
                       matrix(c( 2,
                                 2,
                                 2), byrow = TRUE, nrow = m, ncol = q_emiss[2] - 1),
                       matrix(c(-2, -2,
                                 2,  0,
                                 0, -1), byrow = TRUE, nrow = m, ncol = q_emiss[3] - 1),
                       matrix(c( 2,
                                 2,
                                 2), byrow = TRUE, nrow = m, ncol = q_emiss[4] - 1))
emiss_scalar <- list(c(2), c(3), c(2), c(3))
emiss_w <- rep(list(c(0.2)), n_dep)
manual_emiss_sampler <- pd_RW_emiss_cat(gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                                        emiss_int_mle0 = emiss_int_mle0,
                                        emiss_scalar = emiss_scalar,
                                        emiss_w = emiss_w)

# specifying starting values
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_RWemiss <- mHMM(s_data = nonverbal,
                         gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                         start_val = c(list(start_TM), start_EM),
                         emiss_sampler = manual_emiss_sampler,
                         mcmc = list(J = 11, burn_in = 5))

out_3st_RWemiss
summary(out_3st_RWemiss)

# checking acceptance rate (for illustrative purposes, in the example,
# J is too low for getting a fair indication)
div_J <- function(x, J) x / J
J_it <- 11 - 1 # accept/reject starts at iteration 2 of MCMC algorithm
RW_emiss_accept <- sapply(out_3st_RWemiss$emiss_naccept, div_J, J_it, simplify = FALSE)

# average acceptance rate over all subjects per parameter
# rows represent each of the n_dep dependent variables, columns represent the m states
t(sapply(RW_emiss_accept, apply, MARGIN = 2, mean))

Proposal distribution settings RW Metropolis sampler for mHMM Poisson-lognormal emission distribution(s)

Description

pd_RW_emiss_count provides a framework to manually specify the settings of the proposal distribution of the random walk (RW) Metropolis sampler of Poisson emission distribution(s) of the multilevel hidden Markov model, and creates on object of the class mHMM_pdRW_emiss. The RW metropolis sampler is used for sampling the subject level parameter estimates relating to the emission distributions of the dependent variables k, that is, the Poisson parameters lambda.

Usage

pd_RW_emiss_count(gen, emiss_scalar, emiss_w)

Arguments

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_emiss_count is not used), emiss_scalar set to 2.38, and emiss_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 emission distribution(s) can be tracked using the output parameter emiss_naccept. An acceptance rate of about 45% is considered optimal when a single parameter is being updated (Gelman, Carlin, Stern & Rubin, 2014).

Value

pd_RW_emiss_count returns an object of class mHMM_pdRW_emiss, containing settings of the proposal distribution of the random walk (RW) Metropolis sampler on the categorical emission distribution(s) 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:

gen

A list containing the elements m and n_dep, used for checking equivalent general model properties specified under pd_RW_emiss_count and mHMM.

emiss_scalar

A list containing n_dep elements denoting the scale factor s of the proposal distribution.

emiss_w

A list containing n_dep elements 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 simulated data
# specifying general model properties:
n_t     <- 200     # Number of observations on the dependent variable
m       <- 3        # Number of hidden states
n_dep   <- 2        # Number of dependent variables
n_subj  <- 30        # Number of subjects

gamma   <- matrix(c(0.9, 0.05, 0.05,
                    0.2, 0.7, 0.1,
                    0.2,0.3, 0.5), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c(20,
                             10,
                             5), nrow = m, byrow = TRUE),
                    matrix(c(50,
                             3,
                             20), nrow = m, byrow = TRUE))

# Define between subject variance to use on the simulating function:
# here, the variance is varied over states within the dependent variable.
var_emiss <- list(matrix(c(5.0, 3.0, 1.5), nrow = m),
                  matrix(c(5.0, 5.0, 5.0), nrow = m))

# Simulate count data:
data_count <- sim_mHMM(n_t = n_t,
                       n = n_subj,
                       data_distr = "count",
                       gen = list(m = m, n_dep = n_dep),
                       gamma = gamma,
                       emiss_distr = emiss_distr,
                       var_gamma = 0.1,
                       var_emiss = var_emiss,
                       return_ind_par = TRUE)


# Transition probabilities
start_gamma <- diag(0.8, m)
start_gamma[lower.tri(start_gamma) | upper.tri(start_gamma)] <- (1 - diag(start_gamma)) / (m - 1)

# Emission distribution
start_emiss <- list(matrix(c(20,10, 5), nrow = m, byrow = TRUE),
                    matrix(c(50, 3,20), nrow = m, byrow = TRUE))

# Specify hyper-prior for the count emission distribution
manual_prior_emiss <- prior_emiss_count(
  gen = list(m = m, n_dep = n_dep),
  emiss_mu0 = list(matrix(c(20, 10, 5), byrow = TRUE, ncol = m),
                   matrix(c(50, 3, 20), byrow = TRUE, ncol = m)),
  emiss_K0  = rep(list(0.1),n_dep),
  emiss_nu  = rep(list(0.1),n_dep),
  emiss_V   = rep(list(rep(10, m)),n_dep)
)

# Specify the desired values for the sampler
manual_emiss_sampler <- pd_RW_emiss_count(gen = list(m = m, n_dep = n_dep),
                                        emiss_scalar = rep(list(2.38),n_dep),
                                        emiss_w = rep(list(0.1), n_dep))

# Run model
# 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_count_RWemiss <- mHMM(s_data = data_count$obs,
                          data_distr = 'count',
                          gen = list(m = m, n_dep = n_dep),
                          start_val = c(list(start_gamma), start_emiss),
                          emiss_hyp_prior = manual_prior_emiss,
                          emiss_sampler = manual_emiss_sampler,
                          mcmc = list(J = 11, burn_in = 5),
                          show_progress = TRUE)

# Examine acceptance rates over dependent variable, individual, and states:
lapply(out_3st_count_RWemiss$emiss_naccept, function(e) e/out_3st_count_RWemiss$input$J)

# Finally, take the average acceptance rate by dependent variable and state:
lapply(out_3st_count_RWemiss$emiss_naccept, function(e) colMeans(e/out_3st_count_RWemiss$input$J))

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

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 and mHMM.

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)

Plotting the posterior densities for a fitted multilevel HMM

Description

plot.mHMM plots the posterior densities for a fitted multilevel hidden Markov model for the group and subject level parameters simultaneously. The plotted posterior densities are either for the transition probability matrix gamma, or for the emission distribution probabilities (categorical data) or means and standard deviation (continuous data).

Usage

## S3 method for class 'mHMM'
plot(
  x,
  component = "gamma",
  dep = 1,
  col,
  dep_lab,
  cat_lab,
  lwd1 = 2,
  lwd2 = 1,
  lty1 = 1,
  lty2 = 3,
  legend_cex,
  burn_in,
  ...
)

Arguments

Details

Note that the standard deviation of the (variable and) state specific Normal emission distribution in case of continuous data is fixed over subjects. Hence, for the standard deviation, only the posterior distribution at the group level is plotted.

Value

plot.mHMM returns a plot of the posterior densities. Depending on whether (component = "gamma") or (component = "emiss"), the plotted posterior densities are either for the transition probability matrix gamma or for the emission distribution probabilities, respectively.

See Also

mHMM for fitting the multilevel hidden Markov model, creating the object mHMM.

Examples

###### Example on package example data, see ?nonverbal
# First run the function mHMM on example data

# specifying general model properties:
m <- 2
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
start_TM <- diag(.8, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
start_EM <- list(matrix(c(0.05, 0.90, 0.05, 0.90, 0.05, 0.05), byrow = TRUE,
                        nrow = m, ncol = q_emiss[1]), # vocalizing patient
                 matrix(c(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), byrow = TRUE,
                        nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                 matrix(c(0.1, 0.9, 0.1, 0.9), byrow = TRUE, nrow = m,
                        ncol = q_emiss[4])) # looking therapist

# Run a model without covariate(s):
out_2st <- mHMM(s_data = nonverbal, gen = list(m = m, n_dep = n_dep,
                q_emiss = q_emiss), start_val = c(list(start_TM), start_EM),
                mcmc = list(J = 11, burn_in = 5))

## plot the posterior densities for gamma
plot(out_2st, component = "gamma")

Plotting the transition probabilities gamma for a fitted multilevel HMM

Description

plot.mHMM_gamma plots the transition probability matrix for a fitted multilevel hidden Markov model, by means of an alluvial plot (also known as Sankey diagram or riverplot) using the R package alluvial. The plotted transition probability matrix either represents the probabilities at the group level, i.e., representing the average transition probability matrix over all subjects, or at the subject level. In case of the latter, the user has to specify for which subject the transition probability matrix should be plotted.

Usage

## S3 method for class 'mHMM_gamma'
plot(x, subj_nr = NULL, cex = 0.8, col, hide, ...)

Arguments

Value

plot.mHMM_gamma returns a plot of the transition probability matrix. Depending on whether the input object represents the transition probabilities at the group level or the subject specific transition probability matrices, the returned plot represents either the group transition probability matrix, or the transition probability matrix for a given subject, specified by subject_nr.

See Also

mHMM for fitting the multilevel hidden Markov model, creating the object mHMM, and obtain_gamma to obtain the transition probabilities gamma for a fitted multilevel HMM, creating the object mHMM_gamma.

Examples

#' ###### Example on package data, see ?nonverbal
# specifying general model properties:
m <- 2
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
start_TM <- diag(.8, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
start_EM <- list(matrix(c(0.05, 0.90, 0.05,
                          0.90, 0.05, 0.05), byrow = TRUE,
                        nrow = m, ncol = q_emiss[1]), # vocalizing patient
                 matrix(c(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), byrow = TRUE,
                        nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                 matrix(c(0.1, 0.9,
                          0.1, 0.9), byrow = TRUE, nrow = m,
                        ncol = q_emiss[4])) # looking therapist

# Run a model without covariate(s):
out_2st <- mHMM(s_data = nonverbal,
                gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                start_val = c(list(start_TM), start_EM),
                mcmc = list(J = 11, burn_in = 5))

out_2st
summary(out_2st)

# obtaining the transition probabilities at the group and subject level
est_gamma_group <- obtain_gamma(out_2st, level = "group")

# plot the obtained transition probabilities
plot(est_gamma_group, col = rep(c("green", "blue"), each = m))

Specifying informative hyper-prior on the categorical emission distribution(s) of the multilevel hidden Markov model

Description

prior_emiss_cat provides a framework to manually specify an informative hyper-prior on the categorical emission distribution(s). prior_emiss_cat creates an object of class mHMM_prior_emiss used by the function mHMM, and additionally attaches the class cat to signal use for categorical observations. Note that the hyper-prior distribution on the categorical emission probabilities are on the intercepts (and, if subject level covariates are used, regression coefficients) of the Multinomial logit model used to accommodate the multilevel framework of the data, instead of on the probabilities directly. The set of hyper-prior distributions consists of a multivariate Normal hyper-prior distribution on the vector of means (i.e., intercepts and regression coefficients), and an Inverse Wishart hyper-prior distribution on the covariance matrix.

Usage

prior_emiss_cat(
  gen,
  emiss_mu0,
  emiss_K0 = NULL,
  emiss_V = NULL,
  emiss_nu = NULL,
  n_xx_emiss = NULL
)

Arguments

Details

Estimation of the mHMM proceeds within a Bayesian context, hence a hyper-prior distribution has to be defined for the group level parameters. Default, non-informative priors are used unless specified otherwise by the user. For each dependent variable, each row of the categorical emission probability matrix (i.e., the probability to observe each category (columns) within each of the states (rows)) has its own set of Multinomial logit intercepts, which are assumed to follow a multivariate normal distribution. Hence, the hyper-prior distributions for the intercepts consists of a multivariate Normal hyper-prior distribution on the vector of means, and an Inverse Wishart hyper-prior distribution on the covariance matrix. Note that only the general model properties (number of states m, number of dependent variables n_dep, and number of observed categories for each dependent variable q_emiss) and values of the hypothesized hyper-prior mean values of the Multinomial logit intercepts have to be specified by the user, default values are available for all other hyper-prior distribution parameters.

Given that the hyper-priors are specified on the intercepts of the Multinomial logit model intercepts instead of on the categorical emission probabilities directly, specifying a hyper-prior can seem rather daunting. However, see the function prob_to_int and int_to_prob for translating probabilities to a set of Multinomial logit intercepts and vice versa.

Note that emiss_K0, emiss_nu and emiss_V are assumed equal over the states. When the hyper-prior values for emiss_K0, emiss_nu and emiss_V are not manually specified, the default values are as follows. emiss_K0 set to 1, emiss_nu set to 3 + q_emiss[k] - 1, and the diagonal of gamma_V (i.e., the variance) set to 3 + q_emiss[k] - 1 and the off-diagonal elements (i.e., the covariance) set to 0. In addition, when no manual values for the hyper-prior on the categorical emission distribution are specified at all (that is, the function prior_emiss_cat is not used), all elements of the matrices contained in emiss_mu0 are set to 0 in the function mHMM.

Note that in case covariates are specified, the hyper-prior parameter values of the inverse Wishart distribution on the covariance matrix remain unchanged, as the estimates of the regression coefficients for the covariates are fixed over subjects.

Value

prior_emiss_cat returns an object of class mHMM_prior_emiss, containing informative hyper-prior values for the categorical emission distribution(s) of the multilevel hidden Markov model. The object is specifically created and formatted for use by the function mHMM, and thoroughly checked for correct input dimensions. The object contains the following components:

gen

A list containing the elements m, n_dep, and q_emiss, used for checking equivalent general model properties specified under prior_emiss_cat and mHMM.

emiss_mu0

A list of lists containing the hypothesized hyper-prior mean values of the intercepts of the Multinomial logit model on the categorical emission probabilities.

emiss_K0

A list containing n_dep elements denoting the number of hypothetical prior subjects on which the set of hyper-prior mean intercepts specified in emiss_mu0 are based.

emiss_nu

A list containing n_dep elements denoting the degrees of freedom of the hyper-prior Inverse Wishart distribution on the covariance of the Multinomial logit intercepts.

emiss_V

A list containing n_dep elements containing the variance-covariance of the hyper-prior Inverse Wishart distribution on the covariance of the Multinomial logit intercepts.

n_xx_emiss

A numeric vector denoting the number of (level 2) covariates used to predict the emission distribution of each of the dependent variables. When no covariates are used, n_xx_emiss equals NULL.

See Also

prior_gamma for manually specifying an informative hyper-prior on the transition probability matrix gamma, prob_to_int for transforming a set of probabilities to a set of Multinomial logit regression intercepts, and mHMM for fitting a multilevel hidden Markov model.

Examples

###### Example using package example data, see ?nonverbal
# specifying general model properties:
m <- 3
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# hypothesized mean emission probabilities
prior_prob_emiss_cat <- list(matrix(c(0.10, 0.80, 0.10,
                                      0.80, 0.10, 0.10,
                                      0.40, 0.40, 0.20), byrow = TRUE,
                                    nrow = m, ncol = q_emiss[1]), # vocalizing patient,
                                    # prior belief: state 1 - much talking, state 2 -
                                    # no talking, state 3 - mixed
                             matrix(c(0.30, 0.70,
                                      0.30, 0.70,
                                      0.30, 0.70), byrow = TRUE, nrow = m,
                                    ncol = q_emiss[2]), # looking patient
                                    # prior belief: all 3 states show frequent looking
                                    # behavior
                             matrix(c(0.80, 0.10, 0.10,
                                      0.10, 0.80, 0.10,
                                      0.40, 0.40, 0.20), byrow = TRUE,
                                    nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                                    # prior belief: state 1 - no talking, state 2 -
                                    # frequent talking, state 3 - mixed
                             matrix(c(0.30, 0.70,
                                      0.30, 0.70,
                                      0.30, 0.70), byrow = TRUE, nrow = m,
                                    ncol = q_emiss[4])) # looking therapist
                                    # prior belief: all 3 states show frequent looking
                                    # behavior

# using the function prob_to_int to obtain intercept values for the above specified
# categorical emission distributions
prior_int_emiss <- sapply(prior_prob_emiss_cat, prob_to_int)
emiss_mu0 <- rep(list(vector(mode = "list", length = m)), n_dep)
for(k in 1:n_dep){
  for(i in 1:m){
  emiss_mu0[[k]][[i]] <- matrix(prior_int_emiss[[k]][i,], nrow = 1)
  }
}

emiss_K0 <- rep(list(c(1)), n_dep)
emiss_nu <- list(c(5), c(4), c(5), c(4))
emiss_V <- list(diag(5, q_emiss[1] - 1),
                diag(4, q_emiss[2] - 1),
                diag(5, q_emiss[3] - 1),
                diag(4, q_emiss[4] - 1))

manual_prior_emiss <- prior_emiss_cat(gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                                  emiss_mu0 = emiss_mu0, emiss_K0 = emiss_K0,
                                  emiss_nu = emiss_nu, emiss_V = emiss_V)


# using the informative hyper-prior in a model

# specifying starting values
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_infemiss <- mHMM(s_data = nonverbal,
                    gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                    start_val = c(list(start_TM), start_EM),
                    emiss_hyp_prior = manual_prior_emiss,
                    mcmc = list(J = 11, burn_in = 5))

out_3st_infemiss
summary(out_3st_infemiss)

Specifying informative hyper-prior on the continuous emission distribution(s) of the multilevel hidden Markov model

Description

prior_emiss_cont provides a framework to manually specify an informative hyper-prior on the Normal (i.e., Gaussian) emission distributions. prior_emiss_cont creates an object of class mHMM_prior_emiss used by the function mHMM, and additionally attaches the class cont to signal use for continuous observations. The set of hyper-prior distributions consists of a Normal-Inverse-Gamma distribution (i.e., assuming both unknown population mean and variance between subject level means) on the vector of means (i.e., intercepts and regression coefficients), and an Inverse gamma distribution (i.e., assuming a known mean) on each of the (fixed over subjects) emission variances.

Usage

prior_emiss_cont(
  gen,
  emiss_mu0,
  emiss_K0,
  emiss_V,
  emiss_nu,
  emiss_a0,
  emiss_b0,
  n_xx_emiss = NULL
)

Arguments

Details

Estimation of the mHMM proceeds within a Bayesian context, hence a hyper-prior distribution has to be defined for the group level parameters. To avoid problems with 'label switching' when dealing with continuous emission distribution(s) (i.e., switching of the labels of the hidden states while sampling from the MCMC), the user is forced to specify hyper-prior parameter values when using continuous emission distributions (i.e., default, non-informative priors are not available for continuous emission distributions).

Note that emiss_K0 and emiss_nu are assumed equal over the states. Also note that in case covariates are specified, the hyper-prior parameter values of the inverse Wishart distribution on the covariance matrix remain unchanged, as the estimates of the regression coefficients for the covariates are fixed over subjects.

Value

prior_emiss_cont returns an object of class mHMM_prior_emiss, containing informative hyper-prior values for the continuous emission distribution(s) of the multilevel hidden Markov model. The object is specifically created and formatted for use by the function mHMM, and thoroughly checked for correct input dimensions. The object contains the following components:

gen

A list containing the elements m, and n_dep, used for checking equivalent general model properties specified under prior_emiss_cont and mHMM.

emiss_mu0

A lists containing the hypothesized hyper-prior means of the Normal distribution on the continuous emission probabilities.

emiss_K0

A list containing n_dep elements denoting the number of hypothetical prior subjects on which the set of hyper-prior means specified in emiss_mu0 are based.

emiss_V

A list containing n_dep elements containing the variance of the Inverse Gamma hyper-prior distribution on the between subject variance of the emission distribution means.

emiss_nu

A list containing n_dep elements denoting the degrees of freedom of the Inverse Gamma hyper-prior distribution on the between subject variance of the emission distribution means.

emiss_a0

A list containing n_dep elements denoting the shape values of the Inverse Gamma hyper-prior on each of the (fixed over subjects) emission standard deviation^2 of the Normal emission distributions.

emiss_b0

A list containing n_dep elements denoting the scale values of the Inverse Gamma hyper-prior on each of the (fixed over subjects) emission standard deviation^2 of the Normal emission distributions.

n_xx_emiss

A numeric vector denoting the number of (level 2) covariates used to predict the emission distribution of each of the dependent variables. When no covariates are used, n_xx_emiss equals NULL.

See Also

prior_gamma for manually specifying an informative hyper-prior on the transition probability matrix gamma, and mHMM for fitting a multilevel hidden Markov model.

Examples

###### Example using simulated data
# specifying general model properties:
m <- 3
n_dep <- 2

# hypothesized hyper-prior values for the continuous emission distribution
manual_prior_emiss <- prior_emiss_cont(
                        gen = list(m = m, n_dep = n_dep),
                        emiss_mu0 = list(matrix(c(30, 70, 170), nrow = 1),
                                         matrix(c(7, 8, 18), nrow = 1)),
                        emiss_K0 = list(1, 1),
                         emiss_V =  list(rep(5^2, m), rep(0.5^2, m)),
                        emiss_nu = list(1, 1),
                        emiss_a0 = list(rep(1.5, m), rep(1, m)),
                        emiss_b0 = list(rep(20, m), rep(4, m)))

# to use the informative priors in a model, simulate multivariate continuous data
n_t     <- 100
n       <- 10

gamma   <- matrix(c(0.8, 0.1, 0.1,
                    0.2, 0.7, 0.1,
                    0.2, 0.2, 0.6), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c( 50, 10,
                              100, 10,
                              150, 10), nrow = m, byrow = TRUE),
                    matrix(c(5, 2,
                             10, 5,
                             20, 3), nrow = m, byrow = TRUE))

data_cont <- sim_mHMM(n_t = n_t, n = n, data_distr = 'continuous',
                      gen = list(m = m, n_dep = n_dep),
                      gamma = gamma, emiss_distr = emiss_distr,
                      var_gamma = .1, var_emiss = c(5^2, 0.2^2))

# using the informative hyper-prior in a model
# 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_cont_sim_infemiss <- mHMM(s_data = data_cont$obs,
                    data_distr = "continuous",
                    gen = list(m = m, n_dep = n_dep),
                    start_val = c(list(gamma), emiss_distr),
                    emiss_hyp_prior = manual_prior_emiss,
                    mcmc = list(J = 11, burn_in = 5))

out_3st_cont_sim_infemiss
summary(out_3st_cont_sim_infemiss)

Specifying informative hyper-priors on the count emission distribution(s) of the multilevel hidden Markov model

Description

prior_emiss_count provides a framework to manually specify an informative hyper-prior on the count emission distributions. prior_emiss_count creates an object of class mHMM_prior_emiss used by the function mHMM, and additionally attaches the class count to signal use for count observations. The set of hyper-prior distributions consists of a lognormal-Inverse-Gamma distribution (i.e., assuming both unknown population mean and variance between subject level means) on the vector of Poisson means (i.e., intercepts and regression coefficients).

Usage

prior_emiss_count(
  gen,
  emiss_mu0,
  emiss_K0,
  emiss_V,
  emiss_nu,
  n_xx_emiss = NULL,
  log_scale = FALSE
)

Arguments

Details

Estimation of the mHMM proceeds within a Bayesian context, hence a hyper-prior distribution has to be defined for the group level parameters. To avoid problems with 'label switching' when dealing with continuous emission distribution(s) (i.e., switching of the labels of the hidden states while sampling from the MCMC), the user is forced to specify hyper-prior parameter values when using count emission distributions (i.e., default, non-informative priors are not available for count emission distributions).

Note that emiss_K0 and emiss_nu are assumed equal over the states. Also note that in case covariates are specified, the hyper-prior parameter values of the inverse Wishart distribution on the covariance matrix remain unchanged, as the estimates of the regression coefficients for the covariates are fixed over subjects.

Also note that for simplicity the hyper-prior means and variances of the lognormal distribution, emiss_mu0 and emiss_V, by default have to be specified in the natural (positive real numbers) scale and not in the logarithmic scale. prior_emiss_count returns the corresponding values of the parameters on the logarithmic scale. If the user wants to manually specify these values on the logarithmic scale, please set the argument log_scale to TRUE in prior_emiss_count. If covariates are used to predict the emission distribution, then the logarithmic scale should be used for the inputs emiss_mu0 and emiss_V, and set log_scale = TRUE. To aid the user in transforming the variance to the logarithmic scale, the function 'var_to_logvar()' can be used.

Value

prior_emiss_count returns an object of class mHMM_prior_emiss, containing informative hyper-prior values for the continuous emission distribution(s) of the multilevel hidden Markov model. The object is specifically created and formatted for use by the function mHMM, and thoroughly checked for correct input dimensions. The object contains the following components:

gen

A list containing the elements m, and n_dep, used for checking equivalent general model properties specified under prior_emiss_count and mHMM.

emiss_mu0

A lists containing the hypothesized hyper-prior means of the the Poisson distribution used to model the count emissions.

emiss_K0

A list containing n_dep elements denoting the number of hypothetical prior subjects on which the set of hyper-prior means specified in emiss_mu0 are based.

emiss_V

A list containing n_dep elements containing the variance of the Inverse Gamma hyper-prior distribution on the between subject variance of the emission distribution means.

emiss_nu

A list containing n_dep elements denoting the degrees of freedom of the Inverse Gamma hyper-prior distribution on the between subject variance of the emission distribution means.

n_xx_emiss

A numeric vector denoting the number of (level 2) covariates used to predict the emission distribution of each of the dependent variables. When no covariates are used, n_xx_emiss equals NULL.

See Also

prior_gamma for manually specifying an informative hyper-prior on the transition probability matrix gamma, and mHMM for fitting a multilevel hidden Markov model.

Examples

###### Example using simulated data
# specifying general model properties:
m <- 3
n_dep <- 2

# hypothesized hyper-prior values for the count emission distribution
manual_prior_emiss <- prior_emiss_count(
                        gen = list(m = m, n_dep = n_dep),
                        emiss_mu0 = list(matrix(c(30, 70, 170), nrow = 1),
                                         matrix(c(7, 8, 18), nrow = 1)),
                        emiss_K0 = list(1, 1),
                        emiss_V =  list(rep(16, m), rep(4, m)),
                        emiss_nu = list(0.1, 0.1))

# to use the informative priors in a model, simulate multivariate count data
n_t     <- 100
n       <- 10

# Specify group-level transition and emission means
gamma   <- matrix(c(0.8, 0.1, 0.1,
                    0.2, 0.7, 0.1,
                    0.2, 0.2, 0.6), ncol = m, byrow = TRUE)
emiss_distr <- list(matrix(log(c( 50,
                              100,
                              150)), nrow = m, byrow = TRUE),
                    matrix(log(c(5,
                             10,
                             20)), nrow = m, byrow = TRUE))
# Simulate data
data_count <- sim_mHMM(n_t = n_t, n = n, data_distr = 'count',
                       gen = list(m = m, n_dep = n_dep),
                       gamma = gamma, emiss_distr = emiss_distr,
                       var_gamma = .1, var_emiss = c(.05, 0.01), log_scale = TRUE)

# Specify starting values
start_gamma <- gamma
start_emiss <- list(matrix(c(50,
                             100,
                             150), nrow = m, byrow = TRUE),
                    matrix(c(5,
                             10,
                             20), nrow = m, byrow = TRUE))


# using the informative hyper-prior in a model
# 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_count_sim_infemiss <- mHMM(s_data = data_count$obs,
                    data_distr = "count",
                    gen = list(m = m, n_dep = n_dep),
                    start_val = c(list(start_gamma), start_emiss),
                    emiss_hyp_prior = manual_prior_emiss,
                    mcmc = list(J = 11, burn_in = 5))

out_3st_count_sim_infemiss
summary(out_3st_count_sim_infemiss)

Specifying informative hyper-prior on the transition probability matrix gamma of the multilevel hidden Markov model

Description

prior_gamma provides a framework to manually specify an informative hyper-prior on the transition probability matrix gamma, and creates an object of class mHMM_prior_gamma used by the function mHMM. Note that the hyper-prior distribution on the transition probabilities are on the intercepts (and, if subject level covariates are used, regression coefficients) of the Multinomial logit model used to accommodate the multilevel framework of the data, instead of on the probabilities directly. The set of hyper-prior distributions consists of a multivariate Normal hyper-prior distribution on the vector of means (i.e., intercepts and regression coefficients), and an Inverse Wishart hyper-prior distribution on the covariance matrix.

Usage

prior_gamma(
  m,
  gamma_mu0,
  gamma_K0 = NULL,
  gamma_nu = NULL,
  gamma_V = NULL,
  n_xx_gamma = NULL
)

Arguments

Details

Estimation of the mHMM proceeds within a Bayesian context, hence a hyper-prior distribution has to be defined for the group level parameters. Default, non-informative priors are used unless specified otherwise by the user. Each row of the transition probability matrix has its own set of Multinomial logit intercepts, which are assumed to follow a multivariate normal distribution. Hence, the hyper-prior distributions for the intercepts consists of a multivariate Normal hyper-prior distribution on the vector of means, and an Inverse Wishart hyper-prior distribution on the covariance matrix. Note that only the number of states m and values of the hypothesized hyper-prior mean values of the Multinomial logit intercepts have to be specified by the user, default values are available for all other hyper-prior distribution parameters.

Given that the hyper-priors are specified on the intercepts of the Multinomial logit model intercepts instead of on the probabilities of the transition probability matrix gamma directly, specifying a hyper-prior can seem rather daunting. However, see the function prob_to_int and int_to_prob for translating probabilities to a set of Multinomial logit intercepts and vice versa.

Note that gamma_K0, gamma_nu and gamma_V are assumed equal over the states. When the hyper-prior values for gamma_K0, gamma_nu and gamma_V are not manually specified, the default values are as follows. gamma_K0 set to 1, gamma_nu set to 3 + m - 1, and the diagonal of gamma_V (i.e., the variance) set to 3 + m - 1 and the off-diagonal elements (i.e., the covariance) set to 0. In addition, when no manual values for the hyper-prior on gamma are specified at all (that is, the function prior_gamma is not used), all elements of the matrices contained in gamma_mu0 are set to 0 in the function mHMM.

Note that in case covariates are specified, the hyper-prior parameter values of the inverse Wishart distribution on the covariance matrix remain unchanged, as the estimates of the regression coefficients for the covariates are fixed over subjects.

Value

prior_gamma returns an object of class mHMM_prior_gamma, containing informative hyper-prior values for 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 thoroughly 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 prior_gamma and mHMM.

gamma_mu0

A list containing the hypothesized hyper-prior mean values of the intercepts of the Multinomial logit model on the transition probability matrix gamma.

gamma_K0

A numeric vector denoting the number of hypothetical prior subjects on which the set of hyper-prior mean intercepts specified in gamma_mu0 are based.

gamma_nu

A numeric vector denoting the degrees of freedom of the hyper-prior Inverse Wishart distribution on the covariance of the Multinomial logit intercepts.

gamma_V

A matrix containing the variance-covariance of the hyper-prior Inverse Wishart distribution on the covariance of the Multinomial logit intercepts.

n_xx_gamma

A numeric vector denoting the number of (level 2) covariates used to predict the transition probability matrix gamma. When no covariates are used, n_xx_gamma equals NULL.

See Also

prior_emiss_cat for manually specifying an informative hyper-prior on the categorical emission distribution(s), prob_to_int for transforming a set of probabilities to a set of Multinomial logit regression intercepts, and mHMM for fitting a multilevel hidden Markov model.

Examples

###### Example using package example data, see ?nonverbal
# specifying general model properties:
m <- 3
# representing a prior belief that switching to state 3 does not occur often and
# state 3 has a relative short duration
prior_prob_gamma <- matrix(c(0.70, 0.25, 0.05,
                             0.25, 0.70, 0.05,
                             0.30, 0.30, 0.40), nrow = m, ncol = m, byrow = TRUE)

# using the function prob_to_int to obtain intercept values for the above specified
# transition probability matrix gamma
prior_int_gamma <- prob_to_int(prior_prob_gamma)
gamma_mu0 <- list(matrix(prior_int_gamma[1,], nrow = 1, ncol = m-1),
                  matrix(prior_int_gamma[2,], nrow = 1, ncol = m-1),
                  matrix(prior_int_gamma[3,], nrow = 1, ncol = m-1))

gamma_K0 <- 1
gamma_nu <- 5
gamma_V <- diag(5, m - 1)

manual_prior_gamma <- prior_gamma(m = m, gamma_mu0 = gamma_mu0,
                                  gamma_K0 = gamma_K0, gamma_nu = gamma_nu,
                                  gamma_V = gamma_V)


# using the informative hyper-prior in a model
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
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_infgamma <- 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_hyp_prior = manual_prior_gamma,
                    mcmc = list(J = 11, burn_in = 5))

out_3st_infgamma
summary(out_3st_infgamma)

Transforming a set of probabilities to Multinomial logit regression intercepts

Description

prob_to_int transforms a set of state transition or categorical emission observation probabilities to the corresponding Multinomial logit regression intercepts. Note that the first category is assumed to be the reference category, hence no intercept is returned for the first state or category.

Usage

prob_to_int(prob_matrix)

Arguments

Details

Designed to ease the specification of informative hyper-prior values for the mean intercepts of the transition probability matrix gamma and categorical emission distribution(s) of the multilevel hidden Markov model through the functions prior_gamma and prior_emiss_cat. No check is performed on correct specifications of the dimensions.

Value

prob_to_int returns a matrix containing Multinomial logit regression intercepts, with the number of columns equal to (number of states or categories - 1) and the number of rows equal to the number rows specified in the input matrix. The first state / category is assumed to be the reference category, hence no intercept is returned for this first category.

See Also

int_to_prob for transforming a set of Multinomial logit regression intercepts to a probabilities, prior_gamma and prior_emiss_cat for specifying informative hyper-priors for the the multilevel hidden Markov model and mHMM to fit a multilevel hidden Markov model.

Examples

# example for transition probability matrix gamma with 3 states
m <- 3
gamma_prob <- matrix(c(0.6, 0.2, 0.2,
                       0.1, 0.8, 0.1,
                       0.1, 0.1, 0.8), ncol = m, nrow = m, byrow = TRUE)
gamma_int <- prob_to_int(gamma_prob)
gamma_int

Simulate data using a multilevel hidden Markov model

Description

sim_mHMM simulates data for multiple subjects, for which the data have either categorical or continuous (i.e., normally distributed) observations that follow a hidden Markov model (HMM) with a multilevel structure. The multilevel structure implies that each subject is allowed to have its own set of parameters, and that the parameters at the subject level (level 1) are tied together by a population distribution at level 2 for each of the corresponding parameters. The shape of the population distribution for each of the parameters is a normal distribution. In addition to (natural and/or unexplained) heterogeneity between subjects, the subjects parameters can also depend on a covariate.

Usage

sim_mHMM(
  n_t,
  n,
  data_distr = "categorical",
  gen,
  gamma,
  emiss_distr,
  start_state = NULL,
  xx_vec = NULL,
  beta = NULL,
  var_gamma = 0.1,
  var_emiss = NULL,
  return_ind_par = FALSE,
  m,
  n_dep,
  q_emiss,
  log_scale = FALSE
)

Arguments

Details

In simulating the data, having a multilevel structure means that the parameters for each subject are sampled from the population level distribution of the corresponding parameter. The user specifies the population distribution for each parameter: the average population transition probability matrix and its variance, and the average population emission distribution and its variance. For now, the variance of the mean population parameters is assumed fixed for all components of the transition probability matrix and for all components of the emission distribution.

One can simulate multivariate data. That is, the hidden states depend on more than 1 observed variable simultaneously. The distributions of multiple dependent variables for multivariate data are assumed to be independent, and all distributions for one dataset have to be of the same type (either categorical or continuous).

Note that the subject specific initial state distributions (i.e., the probability of each of the states at the first time point) needed to simulate the data are obtained from the stationary distributions of the subject specific transition probability matrices gamma.

beta: As the first element in each row of gamma is used as reference category in the Multinomial logistic regression, the first matrix in the list beta used to predict transition probability matrix gamma has a number of rows equal to m and the number of columns equal to m - 1. The first element in the first row corresponds to the probability of switching from state one to state two. The second element in the first row corresponds to the probability of switching from state one to state three, and so on. The last element in the first row corresponds to the probability of switching from state one to the last state. The same principle holds for the second matrix in the list beta used to predict categorical emission distribution(s) emiss_distr: the first element in the first row corresponds to the probability of observing category two in state one. The second element in the first row corresponds to the probability of observing category three is state one, and so on. The last element in the first row corresponds to the probability of observing the last category in state one.

Note that when simulating count data (data_distr = 'count'), by default the emission distribution parameters emiss_distr and var_emiss are specified in the natural (positive real numbers) scale (as opposite to the logarithmic scale). If the user wants to manually specify these values on the logarithmic scale, please set the argument log_scale = TRUE. Also note that if covariates are used to predict a count emission distribution, then the logarithmic scale should be used for the inputs emiss_distr, beta, and var_emiss, and also set log_scale = TRUE.

Value

The following components are returned by the function sim_mHMM:

states

A matrix containing the simulated hidden state sequences, with one row per hidden state per subject. The first column indicates subject id number. The second column contains the simulated hidden state sequence, consecutively for all subjects. Hence, the id number is repeated over the rows (with the number of repeats equal to the length of the simulated hidden state sequence T for each subject).

obs

A matrix containing the simulated observed outputs, with one row per simulated observation per subject. The first column indicates subject id number. The second column contains the simulated observation sequence, consecutively for all subjects. Hence, the id number is repeated over rows (with the number of repeats equal to the length of the simulated observation sequence T for each subject).

gamma

A list containing n elements with the simulated subject specific transition probability matrices gamma. Only returned if return_ind_par is set to TRUE.

emiss_distr

A list containing n elements with the simulated subject specific emission probability matrices emiss_distr. Only returned if return_ind_par is set to TRUE.

See Also

mHMM for analyzing multilevel hidden Markov data.

Examples

## Examples on univariate categorical data
# simulating data for 10 subjects with each 100 categorical observations
n_t     <- 100
n       <- 10
m       <- 3
n_dep   <- 1
q_emiss <- 4
gamma   <- matrix(c(0.8, 0.1, 0.1,
                    0.2, 0.7, 0.1,
                    0.2, 0.2, 0.6), ncol = m, byrow = TRUE)
emiss_distr <- list(matrix(c(0.5, 0.5, 0.0, 0.0,
                             0.1, 0.1, 0.8, 0.0,
                             0.0, 0.0, 0.1, 0.9), nrow = m, ncol = q_emiss, byrow = TRUE))
data1 <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                  gamma = gamma, emiss_distr = emiss_distr, var_gamma = 1, var_emiss = 1)
head(data1$obs)
head(data1$states)

# including a covariate to predict (only) the transition probability matrix gamma
beta      <- rep(list(NULL), 2)
beta[[1]] <- matrix(c(0.5, 1.0,
                     -0.5, 0.5,
                      0.0, 1.0), byrow = TRUE, ncol = 2)
xx_vec      <- rep(list(NULL),2)
xx_vec[[1]] <-  c(rep(0,5), rep(1,5))
data2 <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                  gamma = gamma, emiss_distr = emiss_distr, beta = beta, xx_vec = xx_vec,
                  var_gamma = 1, var_emiss = 1)


# simulating subject specific transition probability matrices and emission distributions only
n_t <- 0
n <- 5
m <- 3
n_dep   <- 1
q_emiss <- 4
gamma <- matrix(c(0.8, 0.1, 0.1,
                  0.2, 0.7, 0.1,
                  0.2, 0.2, 0.6), ncol = m, byrow = TRUE)
emiss_distr <- list(matrix(c(0.5, 0.5, 0.0, 0.0,
                             0.1, 0.1, 0.8, 0.0,
                             0.0, 0.0, 0.1, 0.9), nrow = m, ncol = q_emiss, byrow = TRUE))
data3 <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                  gamma = gamma, emiss_distr = emiss_distr, var_gamma = 1, var_emiss = 1)
data3

data4 <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
                  gamma = gamma, emiss_distr = emiss_distr, var_gamma = .5, var_emiss = .5)
data4


## Example on multivariate continuous data
# simulating multivariate continuous data
n_t     <- 100
n       <- 10
m       <- 3
n_dep   <- 2

gamma   <- matrix(c(0.8, 0.1, 0.1,
                    0.2, 0.7, 0.1,
                    0.2, 0.2, 0.6), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c( 50, 10,
                              100, 10,
                              150, 10), nrow = m, byrow = TRUE),
                    matrix(c(5, 2,
                             10, 5,
                             20, 3), nrow = m, byrow = TRUE))

data_cont <- sim_mHMM(n_t = n_t, n = n, data_distr = 'continuous',
                      gen = list(m = m, n_dep = n_dep),
                      gamma = gamma, emiss_distr = emiss_distr,
                      var_gamma = .5, var_emiss = c(5^2, 0.2^2))

head(data_cont$states)
head(data_cont$obs)


## Example on multivariate count data without covariates
n_t     <- 200     # Number of observations on the dependent variable
m       <- 3        # Number of hidden states
n_dep   <- 2        # Number of dependent variables
n_subj  <- 30        # Number of subjects

gamma   <- matrix(c(0.9, 0.05, 0.05,
                    0.2, 0.7, 0.1,
                    0.2,0.3, 0.5), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c(20,
                             10,
                             5), nrow = m, byrow = TRUE),
                    matrix(c(50,
                             3,
                             20), nrow = m, byrow = TRUE))

# Define between subject variance to use on the simulating function:
# here, the variance is varied over states within the dependent variable.
var_emiss <- list(matrix(c(5.0, 3.0, 1.5), nrow = m),
                  matrix(c(5.0, 5.0, 5.0), nrow = m))

# Simulate count data:
data_count <- sim_mHMM(n_t = n_t,
                       n = n_subj,
                       data_distr = "count",
                       gen = list(m = m, n_dep = n_dep),
                       gamma = gamma,
                       emiss_distr = emiss_distr,
                       var_gamma = 0.1,
                       var_emiss = var_emiss,
                       return_ind_par = TRUE)


## Example on multivariate count data with covariates
# Simulate data with one covariate for each count dependent variable

n_t     <- 200     # Number of observations on the dependent variable
m       <- 3        # Number of hidden states
n_dep   <- 3        # Number of dependent variables
n_subj  <- 30        # Number of subjects

gamma   <- matrix(c(0.9, 0.05, 0.05,
                    0.2, 0.7, 0.1,
                    0.2,0.3, 0.5), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c(20,
                             10,
                             5), nrow = m, byrow = TRUE),
                    matrix(c(15,
                             2,
                             5), nrow = m, byrow = TRUE),
                    matrix(c(50,
                             3,
                             20), nrow = m, byrow = TRUE))

# Since we have covariates, we specify inputs in the log scale:
emiss_distr_log <- lapply(emiss_distr, function(e) log(e))

# Define list of vectors of covariate values
set.seed(42)
xx_vec <- c(list(NULL),rep(list(rnorm(n_subj,mean = 0, sd = 0.1)),3))

# Define object beta with regression coefficients for the three dependent variables
beta      <- rep(list(NULL), n_dep+1)
beta[[2]] <- matrix(c(1,-1,0), byrow = TRUE, ncol = 1)
beta[[3]] <- matrix(c(2,0,-2), byrow = TRUE, ncol = 1)
beta[[4]] <- matrix(c(-1,3,1), byrow = TRUE, ncol = 1)

# Calculate logvar to use on the simulating function:
logvar <- var_to_logvar(gen = list(m = m, n_dep = n_dep),
                        emiss_mu = emiss_distr,
                        var_emiss =  list(rep(4, m),
                                          rep(2, m),
                                          rep(5, m)),
                        byrow = FALSE)

# Put logvar in the right format:
logvar <- lapply(logvar, function(q) matrix(q, nrow = m))

# Simulate count data:
data_count <- sim_mHMM(n_t = n_t,
                       n = n_subj,
                       data_distr = "count",
                       gen = list(m = m, n_dep = n_dep),
                       gamma = gamma,
                       emiss_distr = emiss_distr_log,
                       xx_vec = xx_vec,
                       beta = beta,
                       var_gamma = 0.1,
                       var_emiss = logvar,
                       return_ind_par = TRUE,
                       log_scale = TRUE)

head(data_count$states)
head(data_count$obs)

Transform the between-subject variance in the positive scale to the logvariance in the logarithmic scale

Description

var_to_logvar returns the desired between-subject logvariance in the logarithmic scale corresponding to the between-subject variance in the (real) positive scale specified as input. The logvariance is used as the dispersion parameter in the lognormal distribution adopted as prior for the subject-specific Poisson means in the implementation of mHMM and sim_mHMM of mHMMbayes for count distributed data. It takes as inputs the group-level Poisson means emiss_mu and the desired between-subject variance var_emiss, both in the (real) positive scale.

Usage

var_to_logvar(gen, emiss_mu, var_emiss, byrow = TRUE)

Arguments

Value

var_to_logvar Returns a list of n_dep numeric vectors of m elements denoting the state i-specific logvariance (between-subject variance in the logarithmic scale) for the k-dependent variable used as dispersion parameter in the lognormal prior for the Poisson emission distribution.

See Also

mHMM for fitting the multilevel hidden Markov model, and sim_mHMM for simulating data from a multilevel hidden Markov model.

Examples

###### Example on package data


## Example: specifying count priors manually using both positive and log scale:

# Define general parameters:
m       <- 3        # Number of hidden states
n_dep   <- 2        # Number of dependent variables

# Specify priors manually on the positive scale for emuss_mu0 and emiss_V:
manual_prior_emiss_1 <- prior_emiss_count(
  gen = list(m = m, n_dep = n_dep),
  emiss_mu0 = list(matrix(c(30, 70, 170), nrow = 1),
                   matrix(c(7, 8, 18), nrow = 1)),
  emiss_K0 = list(1, 1),
  emiss_V =  list(c(16,25,32),
                  rep(4, m)),
  emiss_nu = list(0.1, 0.1),
  log_scale = FALSE)

# Define logmu and logvar:
logmu <-  list(matrix(log(c(30, 70, 170)), nrow = 1),
               matrix(log(c(7, 8, 18)), nrow = 1))

logvar <- var_to_logvar(gen = list(m = m, n_dep = n_dep),
                        emiss_mu = list(matrix(c(30, 70, 170), nrow = 1),
                                        matrix(c(7, 8, 18), nrow = 1)),
                        var_emiss = list(c(16,25,32),
                                         rep(4, m)),
                        byrow = TRUE)

# Specify priors manually on the log scale for emiss_mu0 and emiss_V:
manual_prior_emiss_2 <- prior_emiss_count(
  gen = list(m = m, n_dep = n_dep),
  emiss_mu0 = logmu,
  emiss_K0 = list(1, 1),
  emiss_V =  logvar,
  emiss_nu = list(0.1, 0.1),
  log_scale = TRUE)

# Check whether they are identical:
identical(manual_prior_emiss_1, manual_prior_emiss_2)

Obtain hidden state sequence for each subject using the Viterbi algorithm

Description

vit_mHMM obtains the most likely state sequence (for each subject) from an object of class mHMM (generated by the function mHMM()), using (an extended version of) the Viterbi algorithm. This is also known as global decoding. Optionally, the state probabilities themselves at each point in time can also be obtained.

Usage

vit_mHMM(object, s_data, burn_in = NULL, return_state_prob = FALSE)

Arguments

Details

Note that local decoding is also possible, by inferring the most frequent state at each point in time for each subject from the sampled state path at each iteration of the MCMC algorithm. This information is contained in the output object return_path of the function mHMM().

Value

The function vit_mHMM returns a matrix containing the most likely state at each point in time. The first column indicates subject id number. Hence, the id number is repeated over rows equal to the number of observations for that subject. The second column contains the most likely state at each point in time. If requested, the subsequent columns contain the state probabilities at each point in time for each subject.

References

Viterbi A (1967). “Error bounds for convolutional codes and an asymptotically optimum decoding algorithm.” IEEE transactions on Information Theory, 13(2), 260–269.

Rabiner LR (1989). “A tutorial on hidden Markov models and selected applications in speech recognition.” Proceedings of the IEEE, 77(2), 257–286.

See Also

mHMM for analyzing multilevel hidden Markov data and obtaining the input needed for vit_mHMM, and sim_mHMM for simulating multilevel hidden Markov data.

Examples

###### Example on package example categorical data, see ?nonverbal
# First fit the multilevel HMM on the example data

# specifying general model properties:
m <- 2
n_dep <- 4
q_emiss <- c(3, 2, 3, 2)

# specifying starting values
start_TM <- diag(.8, m)
start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
start_EM <- list(matrix(c(0.05, 0.90, 0.05, 0.90, 0.05, 0.05), byrow = TRUE,
                        nrow = m, ncol = q_emiss[1]), # vocalizing patient
                 matrix(c(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), byrow = TRUE,
                        nrow = m, ncol = q_emiss[3]), # vocalizing therapist
                 matrix(c(0.1, 0.9, 0.1, 0.9), byrow = TRUE, nrow = m,
                        ncol = q_emiss[4])) # looking therapist

# Fit the multilevel HMM model:
# 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_2st <- mHMM(s_data = nonverbal, gen = list(m = m, n_dep = n_dep,
                q_emiss = q_emiss), start_val = c(list(start_TM), start_EM),
                mcmc = list(J = 3, burn_in = 1))

###### obtain the most likely state sequence with the Viterbi algorithm
states1 <- vit_mHMM(s_data = nonverbal, object = out_2st)


###### Example on simulated multivariate continuous data
# simulating multivariate continuous data
n_t     <- 100
n       <- 10
m       <- 3
n_dep   <- 2

gamma   <- matrix(c(0.8, 0.1, 0.1,
                    0.2, 0.7, 0.1,
                    0.2, 0.2, 0.6), ncol = m, byrow = TRUE)

emiss_distr <- list(matrix(c( 50, 10,
                              100, 10,
                              150, 10), nrow = m, byrow = TRUE),
                    matrix(c(5, 2,
                             10, 5,
                             20, 3), nrow = m, byrow = TRUE))

data_cont <- sim_mHMM(n_t = n_t, n = n, data_distr = 'continuous',
                      gen = list(m = m, n_dep = n_dep), gamma = gamma,
                      emiss_distr = emiss_distr, var_gamma = .1,
                      var_emiss = c(5^2, 0.2^2))

# Specify hyper-prior for the continuous emission distribution
manual_prior_emiss <- prior_emiss_cont(
                        gen = list(m = m, n_dep = n_dep),
                        emiss_mu0 = list(matrix(c(30, 70, 170), nrow = 1),
                                         matrix(c(7, 8, 18), nrow = 1)),
                        emiss_K0 = list(1, 1),
                        emiss_V =  list(rep(5^2, m), rep(0.5^2, m)),
                        emiss_nu = list(1, 1),
                        emiss_a0 = list(rep(1.5, m), rep(1, m)),
                        emiss_b0 = list(rep(20, m), rep(4, m)))

# Run the model on the simulated data:
# 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_cont_sim <- mHMM(s_data = data_cont$obs,
                    data_distr = "continuous",
                    gen = list(m = m, n_dep = n_dep),
                    start_val = c(list(gamma), emiss_distr),
                    emiss_hyp_prior = manual_prior_emiss,
                    mcmc = list(J = 11, burn_in = 5))

###### obtain the most likely state sequence with the Viterbi algorithm,
# including state probabilities at each time point.
states2 <- vit_mHMM(s_data = data_cont$obs, object = out_3st_cont_sim,
                    return_state_prob = TRUE)