###### Example on package (categorical) example data, see ?nonverbal
# \donttest{
# 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
# \donttest{
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)
# }
Run the code above in your browser using DataLab