Nothing
R/mHMM.R
In mHMMbayes: Multilevel Hidden Markov Models Using Bayesian Estimation
Defines functions
mHMM
Documented in
mHMM
#' Multilevel hidden Markov model using Bayesian estimation
#'
#' \code{mHMM} fits a multilevel (also known as mixed or random effects) hidden
#' Markov model (HMM) to intense longitudinal data with categorical 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 \link{mHMMbayes}. See
#' \code{vignette("tutorial-mhmm")} for an introduction to multilevel hidden
#' Markov models and the package, and see \code{vignette("estimation-mhmm")} for
#' an overview of the used estimation algorithms.
#'
#' Covariates specified in \code{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).
#'
#' If covariates are specified and the user wants to set the values for the
#' parameters of the hyper-prior distributions manually, the specification of
#' the elements in the arguments of the hyper-prior parameter values of the
#' normal distribution on the means change as follows: the number of rows in the
#' matrices \code{gamma_mu0} and \code{emiss_mu0} are equal to 1 + the number of
#' covariates used to predict the transition probability matrix for
#' \code{gamma_mu0} and the emission distribution for \code{emiss_mu0} (i.e.,
#' the first row correspond to the hyper-prior mean values of the intercepts,
#' the subsequent rows correspond to the hyper-prior mean values of the
#' regression coefficients connected to each of the covariates), and
#' \code{gamma_K0} and \code{emiss_K0} are now a matrix with the number of
#' hypothetical prior subjects on which the vectors of means (i.e., the rows in
#' \code{gamma_mu0} or \code{emiss_mu0}) are based on the diagonal, and
#' off-diagonal elements equal to 0. Note that the hyper-prior parameter values
#' of the inverse Wishart distribution on the covariance matrix remains
#' unchanged, as the estimates of the regression coefficients for the covariates
#' are fixed over subjects.
#'
#' @param s_data A matrix containing the observations to be modelled, where the
#' rows represent the observations over time. In \code{s_data}, 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 subsequent
#' columns contain the dependent variable(s). Note that the dependent
#' variables have to be numeric, i.e., they cannot be a (set of) factor
#' variable(s). The total number of rows are equal to the sum over the number
#' of observations of each subject, and the number of columns are equal to the
#' number of dependent variables (\code{n_dep}) + 1. The number of
#' observations can vary over subjects.
#' @param gen List containing the following elements denoting the general model
#' properties:
#' \itemize{\item{\code{m}: numeric vector with length 1 denoting the number
#' of hidden states}
#' \item{\code{n_dep}: numeric vector with length 1 denoting the
#' number of dependent variables}
#' \item{\code{q_emiss}: numeric vector with length \code{n_dep} denoting the
#' number of observed categories for the categorical emission distribution
#' for each of the dependent variables.}}
#' @param xx An optional list of (level 2) covariates to predict the transition
#' matrix and/or the emission probabilities. Level 2 covariate(s) means that
#' there is one observation per subject of each covariate. The first element
#' in the list \code{xx} is used to predict the transition matrix. Subsequent
#' elements in the list are used to predict the emission distribution of (each
#' of) the dependent variable(s). Each element in the list is a matrix, with
#' the number of rows equal to the number of subjects. The first column of
#' each matrix represents the intercept, that is, a column only consisting of
#' ones. Subsequent columns correspond to covariates used to predict the
#' transition matrix / emission distribution. See \emph{Details} for more
#' information on the use of covariates.
#'
#' If \code{xx} is omitted completely, \code{xx} defaults to \code{NULL},
#' resembling no covariates. Specific elements in the list can also be left
#' empty (i.e., set to \code{NULL}) to signify that either the transition
#' probability matrix or a specific emission distribution is not predicted by
#' covariates.
#' @param start_val List containing the start values for the transition
#' probability matrix gamma and the emission distribution(s). The first
#' element of the list contains a \code{m} by \code{m} matrix with the start
#' values for gamma. The subsequent elements contain \code{m} by
#' \code{q_emiss[k]} matrices for the start values for each of the \code{k}
#' emission distribution(s). Note that \code{start_val} should not contain
#' nested lists (i.e., lists within lists).
#' @param mcmc List of Markov chain Monte Carlo (MCMC) arguments, containing the
#' following elements:
#' \itemize{\item{\code{J}: numeric vector with length 1 denoting the number
#' of iterations of the MCMC algorithm}
#' \item{\code{burn_in}: numeric vector with length 1 denoting the
#' burn-in period for the MCMC algorithm.}}
#' @param return_path A logical scalar. Should the sampled state sequence
#' obtained at each iteration and for each subject be returned by the function
#' (\code{sample_path = TRUE}) or not (\code{sample_path = FALSE}). Note that
#' the sampled state sequence is quite a large object, hence the default
#' setting is \code{sample_path = FALSE}. Can be used for local decoding
#' purposes.
#' @param print_iter A logical scaler. Should the function print the progress of
#' the algorithm by returning the current iteration number at every 10th
#' iteration (\code{print_iter = TRUE}) or not (\code{print_iter = FALSE}).
#' Defaults to \code{print_iter = FALSE}.
#' @param gamma_hyp_prior An optional list containing user specified parameters
#' of the hyper-prior distribution on the multivariate normal distribution
#' of the intercepts (and regression coefficients given that covariates are
#' used) of the multinomial regression model of the transition probability
#' matrix gamma. The hyper-prior of the mean intercepts is a multivariate
#' Normal distribution, the hyper-prior of the covariance matrix between the
#' set of (state specific) intercepts is an Inverse Wishart distribution.
#'
#' Hence, the list \code{gamma_hyp_prior} contains the following elements:
#' \itemize{\item{\code{gamma_mu0}: a list containing m matrices; one matrix
#' for each row of the transition probability matrix gamma. Each matrix
#' contains the hypothesized mean values of the intercepts. Hence, each matrix
#' consists of one row (when not including covariates in the model) and
#' \code{m} - 1 columns}
#' \item{\code{gamma_K0}: numeric vector with length 1 denoting the number of
#' hypothetical prior subjects on which the vector of means \code{gamma_mu0} is
#' based}
#' \item{\code{gamma_nu}: numeric vector with length 1 denoting the degrees of
#' freedom of the Inverse Wishart distribution}
#' \item{\code{gamma_V}: matrix of \code{m} - 1 by \code{m} - 1 containing the
#' hypothesized variance-covariance matrix between the set of intercepts.}}
#' Note that \code{gamma_K0}, \code{gamma_nu} and \code{gamma_V} are assumed
#' equal over the states. The mean values of the intercepts (and regression
#' coefficients of the covariates) denoted by \code{gamma_mu0} are allowed to
#' vary over the states.
#'
#' The default values for the hyper-prior on gamma are: all elements of the
#' matrices contained in \code{gamma_mu0} set to 0, \code{gamma_K0} set to 1,
#' \code{gamma_nu} set to 3 + m - 1, and the diagonal of \code{gamma_V} (i.e.,
#' the variance) set to 3 + m - 1 and the off-diagonal elements (i.e., the
#' covariance) set to 0.
#'
#' See \emph{Details} below if covariates are used for changes in the settings
#' of the arguments of \code{gamma_hyp_prior}.
#' @param emiss_hyp_prior An optional list containing user specified parameters
#' of the hyper-prior distribution on the multivariate normal distribution
#' of the intercepts (and regression coefficients given that covariates are
#' used) of the multinomial regression model of the emission distribution.
#' The hyper-prior for the mean intercepts is a multivariate
#' Normal distribution, the hyper-prior for the covariance matrix between the
#' set of (state specific) intercepts is an Inverse Wishart distribution.
#'
#' Hence, the list \code{emiss_hyp_prior} contains the following elements:
#' \itemize{\item{\code{emiss_mu0}: a list of lists: \code{emiss_mu0} contains
#' \code{ndep} lists, i.e., one list for each dependent variable \code{k}. Each
#' of these lists contains m matrices; one matrix for each set of emission
#' probabilities within a state. The matrices contain the hypothesized mean
#' values of the intercepts. Hence, each matrix consists of one row (when not
#' including covariates in the model) and \code{q_emiss[k]} - 1 columns}
#' \item{\code{emiss_K0}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables, where each element is a numeric vector
#' with length 1 denoting the number of hypothetical prior subjects on which
#' the vector of means \code{emiss_mu0} is based}
#' \item{\code{emiss_nu}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables, where each element is a numeric vector
#' with length 1 denoting the degrees of freedom of the Inverse Wishart
#' distribution}
#' \item{\code{emiss_V}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables \code{k}, where each element is a matrix
#' of \code{q_emiss[k]} - 1 by \code{q_emiss[k]} - 1 containing the
#' hypothesized variance-covariance matrix between the set of intercepts.}}
#' Note that \code{emiss_K0}, \code{emiss_nu} and \code{emiss_V} are assumed
#' equal over the states. The mean values of the intercepts (and regression
#' coefficients of the covariates) denoted by \code{emiss_mu0} are allowed to
#' vary over the states.
#'
#' The default values for the hyper-prior on the emission distribution(s) are:
#' all elements of the matrices contained in \code{emiss_K0} set to 0,
#' \code{emiss_K0} set to 1, \code{emiss_nu} set to 3 + \code{q_emiss[k]} - 1,
#' and the diagonal of \code{gamma_V} (i.e., the variance) set to 3 +
#' \code{q_emiss[k]} - 1 and the off-diagonal elements (i.e., the covariance)
#' set to 0.
#'
#' See \emph{Details} below if covariates are used for changes in the settings
#' of the arguments of \code{emiss_hyp_prior}.
#' @param gamma_sampler An optional list containing user specified settings for
#' the proposal distribution of the random walk (RW) Metropolis sampler for
#' the subject level parameter estimates of the intercepts modeling the
#' transition probability matrix. The list \code{gamma_sampler} contains the
#' following elements:
#' \itemize{\item{\code{gamma_int_mle0}: a numeric vector with length \code{m}
#' - 1 denoting the start values for the maximum likelihood estimates of the
#' intercepts for the transition probability matrix gamma, based on the pooled
#' data (data over all subjects)}
#' \item{\code{gamma_scalar}: a numeric vector with length 1 denoting the scale
#' factor \code{s}. That is, The scale of the proposal distribution is composed
#' of a covariance matrix Sigma, which is then tuned by multiplying it by a
#' scaling factor \code{s}^2}
#' \item{\code{gamma_w}: a numeric vector with length 1 denoting the weight for
#' the overall log likelihood (i.e., log likelihood based on the pooled data
#' over all subjects) in the fractional likelihood.}}
#' Default settings are: all elements in \code{gamma_int_mle0} set to 0,
#' \code{gamma_scalar} set to 2.93 / sqrt(\code{m} - 1), and \code{gamma_w} set to
#' 0.1. See the section \emph{Scaling the proposal distribution of the RW
#' Metropolis sampler} in \code{vignette("estimation-mhmm")} for details.
#' @param emiss_sampler An optional list containing user specified settings for
#' the proposal distribution of the random walk (RW) Metropolis sampler for
#' the subject level parameter estimates of the intercepts modeling the
#' emission distributions of the dependent variables \code{k}. The list
#' \code{emiss_sampler} contains the following elements:
#' \itemize{\item{\code{emiss_int_mle0}: a list containing \code{ndep} elements
#' corresponding to each of the dependent variables \code{k}, where each
#' element is a a numeric vector with length \code{q_emiss[k]} - 1 denoting the
#' start values for the maximum likelihood estimates of the intercepts for
#' the emission distribution, based on the pooled data (data over all
#' subjects)}
#' \item{\code{emiss_scalar}: a list containing \code{ndep} elements
#' corresponding to each of the dependent variables, where each element is a
#' numeric vector with length 1 denoting the scale factor \code{s}. That is,
#' The scale of the proposal distribution is composed of a covariance matrix
#' Sigma, which is then tuned by multiplying it by a scaling factor \code{s}^2}
#' \item{\code{emiss_w}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables, where each element is a numeric vector
#' with length 1 denoting the weight for the overall log likelihood (i.e., log
#' likelihood based on the pooled data over all subjects) in the fractional
#' likelihood.}}
#' Default settings are: all elements in \code{emiss_int_mle0} set to 0,
#' \code{emiss_scalar} set to 2.93 / sqrt(\code{q_emiss[k]} - 1), and
#' \code{emiss_w} set to 0.1. See the section \emph{Scaling the proposal
#' distribution of the RW Metropolis sampler} in
#' \code{vignette("estimation-mhmm")} for details.
#'
#' @return \code{mHMM} returns an object of class \code{mHMM}, which has
#' \code{print} and \code{summary} methods to see the results.
#' The object contains the following components:
#' \describe{
#' \item{\code{PD_subj}}{A list containing one matrix per subject with the
#' subject level parameter estimates and the log likelihood 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 (parameter) estimates of subsequently the emission probabilities, the
#' transition probabilities and the log likelihood.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{input}}{Overview of used input specifications: the number of
#' states \code{m}, the number of used dependent variables \code{n_dep}, the
#' number of output categories for each of the dependent variables
#' \code{q_emiss}, the number of iterations \code{J} and the specified burn in
#' period \code{burn_in} of the hybrid Metropolis within Gibbs sampler, the
#' number of subjects \code{n_subj}, the observation length for each subject
#' \code{n_vary}, and the column names of the dependent variables
#' \code{dep_labels}.}
#' \item{\code{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 \code{return_path = TRUE}. }
#' }
#'
#' @seealso \code{\link{sim_mHMM}} for simulating multilevel hidden Markov data,
#' \code{\link{vit_mHMM}} for obtaining the most likely hidden state sequence
#' for each subject using the Viterbi algorithm, \code{\link{obtain_gamma}}
#' and \code{\link{obtain_emiss}} for obtaining the transition or emission
#' distribution probabilities of a fitted model at the group or subject level,
#' and \code{\link{plot.mHMM}} for plotting the posterior densities of a
#' fitted model.
#'
#' @references
#' \insertRef{rabiner1989}{mHMMbayes}
#'
#' \insertRef{scott2002}{mHMMbayes}
#'
#' \insertRef{altman2007}{mHMMbayes}
#'
#' \insertRef{rossi2012}{mHMMbayes}
#'
#' \insertRef{zucchini2017}{mHMMbayes}
#'
#' @examples
#' ###### Example on package example data
#' \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. Here, the covariate (standardized CDI
#' # change) predicts the emission distribution for each of the 4 dependent
#' # variables:
#'
#' n_subj <- 10
#' xx <- rep(list(matrix(1, ncol = 1, nrow = n_subj)), (n_dep + 1))
#' for(i in 2:(n_dep + 1)){
#' xx[[i]] <- cbind(xx[[i]], nonverbal_cov$std_CDI_change)
#' }
#' 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 simulated data
#' # Simulate data for 10 subjects with each 100 observations:
#' n_t <- 100
#' n <- 10
#' m <- 2
#' q_emiss <- 3
#' gamma <- matrix(c(0.8, 0.2,
#' 0.3, 0.7), ncol = m, byrow = TRUE)
#' emiss_distr <- 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, m = m, 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 = list(gamma, emiss_distr),
#' mcmc = list(J = 11, burn_in = 5))
#'
#'
#' @export
#'
#'
mHMM <- function(s_data, gen, xx = NULL, start_val, mcmc, return_path = FALSE, print_iter = FALSE,
gamma_hyp_prior = NULL, emiss_hyp_prior = NULL,
gamma_sampler = NULL, emiss_sampler = NULL){
# Initialize data -----------------------------------
# dependent variable(s), sample size, dimensions gamma and conditional distribuiton
n_dep <- gen$n_dep
dep_labels <- colnames(s_data[,2:(n_dep+1)])
id <- unique(s_data[,1])
n_subj <- length(id)
subj_data <- rep(list(NULL), n_subj)
if(sum(sapply(s_data, is.factor)) > 0 ){
stop("Your data contains factorial variables, which cannot be used as input in the function mHMM. All variables have to be numerical.")
}
for(s in 1:n_subj){
subj_data[[s]]$y <- as.matrix(s_data[s_data[,1] == id[s],][,-1], ncol = n_dep)
}
ypooled <- n <- NULL
n_vary <- numeric(n_subj)
m <- gen$m
q_emiss <- gen$q_emiss
emiss_int_mle <- rep(list(NULL), n_dep)
emiss_mhess <- rep(list(NULL), n_dep)
for(q in 1:n_dep){
emiss_int_mle[[q]] <- matrix(, m, (q_emiss[q] - 1))
emiss_mhess[[q]] <- matrix(, (q_emiss[q] - 1) * m, (q_emiss[q] - 1))
}
for(s in 1:n_subj){
ypooled <- rbind(ypooled, subj_data[[s]]$y)
n <- dim(subj_data[[s]]$y)[1]
n_vary[s] <- n
subj_data[[s]] <- c(subj_data[[s]], n = n, list(gamma_converge = numeric(m), gamma_int_mle = matrix(, m, (m - 1)),
gamma_mhess = matrix(, (m - 1) * m, (m - 1)), emiss_converge =
rep(list(numeric(m)), n_dep), emiss_int_mle = emiss_int_mle, emiss_mhess = emiss_mhess))
}
n_total <- dim(ypooled)[1]
# covariates
n_dep1 <- 1 + n_dep
nx <- numeric(n_dep1)
if (is.null(xx)){
xx <- rep(list(matrix(1, ncol = 1, nrow = n_subj)), n_dep1)
nx[] <- 1
} else {
if(!is.list(xx) | length(xx) != n_dep1){
stop("If xx is specified, xx should be a list, with the number of elements equal to the number of dependent variables + 1")
}
for(i in 1:n_dep1){
if (is.null(xx[[i]])){
xx[[i]] <- matrix(1, ncol = 1, nrow = n_subj)
nx[i] <- 1
} else {
nx[i] <- ncol(xx[[i]])
if (sum(xx[[i]][,1] != 1)){
stop("If xx is specified, the first column in each element of xx has to represent the intercept. That is, a column that only consists of the value 1")
}
if(nx[i] > 1){
for(j in 2:nx[i]){
if(is.factor(xx[[i]][,j])){
stop("Factors currently cannot be used as covariates, see help file for alternatives")
}
if((length(unique(xx[[i]][,j])) == 2) & (sum(xx[[i]][,j] != 0 & xx[[i]][,j] !=1) > 0)){
stop("Dichotomous covariates in xx need to be coded as 0 / 1 variables. That is, only conisting of the values 0 and 1")
}
if(length(unique(xx[[i]][,j])) > 2){
xx[[i]][,j] <- xx[[i]][,j] - mean(xx[[i]][,j])
}
}
}
}
}
}
# Initialize mcmc argumetns
J <- mcmc$J
burn_in <- mcmc$burn_in
# Initalize priors and hyper priors --------------------------------
# Initialize gamma sampler
if(is.null(gamma_sampler)) {
gamma_int_mle0 <- rep(0, m - 1)
gamma_scalar <- 2.93 / sqrt(m - 1)
gamma_w <- .1
} else {
gamma_int_mle0 <- gamma_sampler$gamma_int_mle0
gamma_scalar <- gamma_sampler$gamma_scalar
gamma_w <- gamma_sampler$gamma_w
}
# Initialize emiss sampler
if(is.null(emiss_sampler)){
emiss_int_mle0 <- rep(list(NULL), n_dep)
emiss_scalar <- rep(list(NULL), n_dep)
for(q in 1:n_dep){
emiss_int_mle0[[q]] <- rep(0, q_emiss[q] - 1)
emiss_scalar[[q]] <- 2.93 / sqrt(q_emiss[q] - 1)
}
emiss_w <- .1
} else {
emiss_int_mle0 <- emiss_sampler$emiss_int_mle0
emiss_scalar <- emiss_sampler$emiss_scalar
emiss_w <- emiss_sampler$emiss_w
}
# Initialize Gamma hyper prior
if(is.null(gamma_hyp_prior)){
gamma_mu0 <- rep(list(matrix(0,nrow = nx[1], ncol = m - 1)), m)
gamma_K0 <- diag(1, nx[1])
gamma_nu <- 3 + m - 1
gamma_V <- gamma_nu * diag(m - 1)
} else {
###### BUILD in a warning / check if gamma_mu0 is a matrix when given, with nrows equal to the number of covariates
gamma_mu0 <- gamma_hyp_prior$gamma_mu0
gamma_K0 <- gamma_hyp_prior$gamma_K0
gamma_nu <- gamma_hyp_prior$gamma_nu
gamma_V <- gamma_hyp_prior$gamma_V
}
# Initialize Pr hyper prior
# emiss_mu0: for each dependent variable, emiss_mu0 is a list, with one element for each state.
# Each element is a matrix, with number of rows equal to the number of covariates (with the intercept being one cov),
# and the number of columns equal to q_emiss[q] - 1.
emiss_mu0 <- rep(list(vector("list", m)), n_dep)
emiss_nu <- rep(list(NULL), n_dep)
emiss_V <- rep(list(NULL), n_dep)
emiss_K0 <- rep(list(NULL), n_dep)
if(is.null(emiss_hyp_prior)){
for(q in 1:n_dep){
for(i in 1:m){
emiss_mu0[[q]][[i]] <- matrix(0, ncol = q_emiss[q] - 1, nrow = nx[1 + q])
}
emiss_nu[[q]] <- 3 + q_emiss[q] - 1
emiss_V[[q]] <- emiss_nu[[q]] * diag(q_emiss[q] - 1)
emiss_K0[[q]] <- diag(1, nx[1 + q])
}
} else {
for(q in 1:n_dep){
# emiss_hyp_prior[[q]]$emiss_mu0 has to contain a list with lenght equal to m, and each list contains matrix with number of rows equal to number of covariates for that dep. var.
# stil build in a CHECK, with warning / stop / switch to default prior
emiss_mu0[[q]] <- emiss_hyp_prior[[q]]$emiss_mu0
emiss_nu[[q]] <- emiss_hyp_prior[[q]]$emiss_nu
emiss_V[[q]] <- emiss_hyp_prior[[q]]$emiss_V
emiss_K0[[q]] <- diag(emiss_hyp_prior$emiss_K0, nx[1 + q])
}
}
# Define objects used to store data in mcmc algorithm, not returned ----------------------------
# overall
c <- llk <- numeric(1)
sample_path <- lapply(n_vary, dif_matrix, cols = J)
trans <- rep(list(vector("list", m)), n_subj)
# gamma
gamma_mle_pooled <-gamma_int_mle_pooled <- gamma_mhess_pooled <- gamma_pooled_ll <- vector("list", m)
gamma_c_int <- rep(list(matrix(, n_subj, (m-1))), m)
gamma_mu_int_bar <- gamma_V_int <- vector("list", m)
gamma_mu_prob_bar <- rep(list(numeric(m)), m)
gamma_naccept <- matrix(0, n_subj, m)
# emiss
cond_y <- lapply(rep(n_dep, n_subj), nested_list, m = m)
emiss_mle_pooled <- emiss_int_mle_pooled <- emiss_mhess_pooled <- emiss_pooled_ll <- rep(list(vector("list", n_dep)), m)
emiss_c_int <- rep(list(lapply(q_emiss - 1, dif_matrix, rows = n_subj)), m)
emiss_mu_int_bar <- emiss_V_int <- rep(list(vector("list", n_dep)), m)
emiss_mu_prob_bar <- rep(list(lapply(q_emiss, dif_vector)), m)
emiss_naccept <- rep(list(matrix(0, n_subj, m)), n_dep)
# Define objects that are returned from mcmc algorithm ----------------------------
# Define object for subject specific posterior density, put start values on first row
if(length(start_val) != n_dep + 1){
stop("The number of elements in the list start_val should be equal to 1 + the number of dependent variables,
and should not contain nested lists (i.e., lists within lists)")
}
PD <- matrix(, nrow = J, ncol = sum(m * q_emiss) + m * m + 1)
PD_emiss_names <- paste("q", 1, "_emiss", rep(1:q_emiss[1], m), "_S", rep(1:m, each = q_emiss[1]), sep = "")
if(n_dep > 1){
for(q in 2:n_dep){
PD_emiss_names <- c(PD_emiss_names, paste("q", q, "_emiss", rep(1:q_emiss[q], m), "_S", rep(1:m, each = q_emiss[q]), sep = ""))
}
}
colnames(PD) <- c(PD_emiss_names, paste("S", rep(1:m, each = m), "toS", rep(1:m, m), sep = ""), "LL")
PD[1, ((sum(m * q_emiss) + 1)) :((sum(m * q_emiss) + m * m))] <- unlist(sapply(start_val, t))[1:(m*m)]
PD[1, 1:((sum(m * q_emiss)))] <- unlist(sapply(start_val, t))[(m*m + 1): (m*m + sum(m * q_emiss))]
PD_subj <- rep(list(PD), n_subj)
# Define object for population posterior density (probabilities and regression coefficients parameterization )
gamma_prob_bar <- matrix(, nrow = J, ncol = (m * m))
colnames(gamma_prob_bar) <- paste("S", rep(1:m, each = m), "toS", rep(1:m, m), sep = "")
gamma_prob_bar[1,] <- PD[1,(sum(m*q_emiss) + 1):(sum(m * q_emiss) + m * m)]
emiss_prob_bar <- lapply(q_emiss * m, dif_matrix, rows = J)
names(emiss_prob_bar) <- dep_labels
for(q in 1:n_dep){
colnames(emiss_prob_bar[[q]]) <- paste("Emiss", rep(1:q_emiss[q], m), "_S", rep(1:m, each = q_emiss[q]), sep = "")
start <- c(0, q_emiss * m)
emiss_prob_bar[[q]][1,] <- PD[1,(sum(start[1:q]) + 1):(sum(start[1:q]) + (m * q_emiss[q]))]
}
gamma_int_bar <- matrix(, nrow = J, ncol = ((m-1) * m))
colnames(gamma_int_bar) <- paste("int_S", rep(1:m, each = m-1), "toS", rep(2:m, m), sep = "")
gamma_int_bar[1,] <- as.vector(prob_to_int(matrix(gamma_prob_bar[1,], byrow = TRUE, ncol = m, nrow = m)))
if(nx[1] > 1){
gamma_cov_bar <- matrix(, nrow = J, ncol = ((m-1) * m) * (nx[1] - 1))
colnames(gamma_cov_bar) <- paste( paste("cov", 1 : (nx[1] - 1), "_", sep = ""), "S", rep(1:m, each = (m-1) * (nx[1] - 1)), "toS", rep(2:m, m * (nx[1] - 1)), sep = "")
gamma_cov_bar[1,] <- 0
} else{
gamma_cov_bar <- "No covariates where used to predict the transition probability matrix"
}
emiss_int_bar <- lapply((q_emiss-1) * m, dif_matrix, rows = J)
names(emiss_int_bar) <- dep_labels
for(q in 1:n_dep){
colnames(emiss_int_bar[[q]]) <- paste("int_Emiss", rep(2:q_emiss[q], m), "_S", rep(1:m, each = q_emiss[q] - 1), sep = "")
emiss_int_bar[[q]][1,] <- as.vector(prob_to_int(matrix(emiss_prob_bar[[q]][1,], byrow = TRUE, ncol = q_emiss[q], nrow = m)))
}
if(sum(nx[-1]) > n_dep){
emiss_cov_bar <- lapply((q_emiss-1) * m * (nx[-1] - 1 ), dif_matrix, rows = J)
names(emiss_cov_bar) <- dep_labels
for(q in 1:n_dep){
if(nx[1 + q] > 1){
colnames(emiss_cov_bar[[q]]) <- paste( paste("cov", 1 : (nx[1 + q] - 1), "_", sep = ""), "emiss", rep(2:q_emiss[q], m * (nx[1 + q] - 1)), "_S", rep(1:m, each = (q_emiss[q] - 1) * (nx[1 + q] - 1)), sep = "")
emiss_cov_bar[[q]][1,] <- 0
} else {
emiss_cov_bar[[q]] <- "No covariates where used to predict the emission probabilities for this outcome"
}
}
} else{
emiss_cov_bar <- "No covariates where used to predict the emission probabilities"
}
# Define object for subject specific posterior density (regression coefficients parameterization )
gamma_int_subj <- rep(list(gamma_int_bar), n_subj)
emiss_int_subj <- rep(list(emiss_int_bar), n_subj)
# Put starting values in place for fist run forward algorithm
emiss_sep <- vector("list", n_dep)
for(q in 1:n_dep){
start <- c(0, q_emiss * m)
emiss_sep[[q]] <- matrix(PD[1,(sum(start[1:q]) + 1):(sum(start[1:q]) + (m * q_emiss[q]))], byrow = TRUE, ncol = q_emiss[q], nrow = m)
}
emiss <- rep(list(emiss_sep), n_subj)
gamma <- rep(list(matrix(PD[1,(sum(m*q_emiss) + 1):(sum(m * q_emiss) + m * m)], byrow = TRUE, ncol = m)), n_subj)
delta <- rep(list(solve(t(diag(m) - gamma[[1]] + 1), rep(1, m))), n_subj)
# Start analysis --------------------------------------------
# Run the MCMC algorithm
itime <- proc.time()[3]
for (iter in 2 : J){
# For each subject, obtain sampled state sequence with subject individual parameters ----------
for(s in 1:n_subj){
# Run forward algorithm, obtain subject specific forward proababilities and log likelihood
forward <- cat_Mult_HMM_fw(x = subj_data[[s]]$y, m = m, emiss = emiss[[s]], gamma = gamma[[s]], n_dep = n_dep, delta=NULL)
alpha <- forward$forward_p
c <- max(forward$la[, subj_data[[s]]$n])
llk <- c + log(sum(exp(forward$la[, subj_data[[s]]$n] - c)))
PD_subj[[s]][iter, sum(m * q_emiss) + m * m + 1] <- llk
# Using the forward probabilites, sample the state sequence in a backward manner.
# In addition, saves state transitions in trans, and conditional observations within states in cond_y
trans[[s]] <- vector("list", m)
sample_path[[s]][n_vary[[s]], iter] <- sample(1:m, 1, prob = c(alpha[, n_vary[[s]]]))
for(t in (subj_data[[s]]$n - 1):1){
sample_path[[s]][t,iter] <- sample(1:m, 1, prob = (alpha[, t] * gamma[[s]][,sample_path[[s]][t + 1, iter]]))
trans[[s]][[sample_path[[s]][t,iter]]] <- c(trans[[s]][[sample_path[[s]][t, iter]]], sample_path[[s]][t + 1, iter])
}
for (i in 1:m){
trans[[s]][[i]] <- c(trans[[s]][[i]], 1:m)
trans[[s]][[i]] <- rev(trans[[s]][[i]])
for(q in 1:n_dep){
cond_y[[s]][[i]][[q]] <- c(subj_data[[s]]$y[sample_path[[s]][, iter] == i, q], 1:q_emiss[q])
}
}
}
# The remainder of the mcmc algorithm is state specific
for(i in 1:m){
# Obtain MLE of the covariance matrices and log likelihood of gamma and emiss at subject and population level -----------------
# used to scale the propasal distribution of the RW Metropolis sampler
# population level, transition matrix
trans_pooled <- factor(c(unlist(sapply(trans, "[[", i)), c(1:m)))
gamma_mle_pooled[[i]] <- optim(gamma_int_mle0, llmnl_int, Obs = trans_pooled, n_cat = m, method = "BFGS", hessian = TRUE, control = list(fnscale = -1))
gamma_int_mle_pooled[[i]] <- gamma_mle_pooled[[i]]$par
gamma_pooled_ll[[i]] <- gamma_mle_pooled[[i]]$value
# population level, conditional probabilities, seperate for each dependent variable
for(q in 1:n_dep){
cond_y_pooled <- numeric()
### MOET OOK ECHT BETER KUNNEN, eerst # cond_y_pooled <- unlist(sapply(cond_y, "[[", m))
for(s in 1:n_subj){
cond_y_pooled <- c(cond_y_pooled, cond_y[[s]][[i]][[q]])
}
emiss_mle_pooled[[i]][[q]] <- optim(emiss_int_mle0[[q]], llmnl_int, Obs = c(cond_y_pooled, c(1:q_emiss[q])), n_cat = q_emiss[q],
method = "BFGS", hessian = TRUE, control = list(fnscale = -1))
emiss_int_mle_pooled[[i]][[q]] <- emiss_mle_pooled[[i]][[q]]$par
emiss_pooled_ll[[i]][[q]] <- emiss_mle_pooled[[i]][[q]]$value
}
# subject level, transition matrix
for (s in 1:n_subj){
wgt <- subj_data[[s]]$n / n_total
gamma_out <- optim(gamma_int_mle_pooled[[i]], llmnl_int_frac, Obs = c(trans[[s]][[i]], c(1:m)), n_cat = m,
pooled_likel = gamma_pooled_ll[[i]], w = gamma_w, wgt = wgt, method="BFGS", hessian = TRUE, control = list(fnscale = -1))
if(gamma_out$convergence == 0){
subj_data[[s]]$gamma_converge[i] <- 1
subj_data[[s]]$gamma_mhess[(1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1)), ] <- mnlHess_int(int = gamma_out$par, Obs = c(trans[[s]][[i]], c(1:m)), n_cat = m)
subj_data[[s]]$gamma_int_mle[i,] <- gamma_out$par
} else {
subj_data[[s]]$gamma_converge[i] <- 0
subj_data[[s]]$gamma_mhess[(1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1)),] <- diag(m-1)
subj_data[[s]]$gamma_int_mle[i,] <- rep(0, m - 1)
}
# if this is first iteration, use MLE for current values RW metropolis sampler
if (iter == 2){
gamma_c_int[[i]][s,] <- gamma_out$par
}
# subject level, conditional probabilities, seperate for each dependent variable
for(q in 1:n_dep){
emiss_out <- optim(emiss_int_mle_pooled[[i]][[q]], llmnl_int_frac, Obs = c(cond_y[[s]][[i]][[q]], c(1:q_emiss[q])),
n_cat = q_emiss[q], pooled_likel = emiss_pooled_ll[[i]][[q]], w = emiss_w, wgt = wgt, method = "BFGS", hessian = TRUE, control = list(fnscale = -1))
if(emiss_out$convergence == 0){
subj_data[[s]]$emiss_converge[[q]][i] <- 1
subj_data[[s]]$emiss_mhess[[q]][(1 + (i - 1) * (q_emiss[q] - 1)):((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1)), ] <- mnlHess_int(int = emiss_out$par, Obs = c(cond_y[[s]][[i]][[q]], c(1:q_emiss[q])), n_cat = q_emiss[q])
subj_data[[s]]$emiss_int_mle[[q]][i,] <- emiss_out$par
} else {
subj_data[[s]]$emiss_converge[[q]][i] <- 0
subj_data[[s]]$emiss_mhess[[q]][(1 + (i - 1) * (q_emiss[q] - 1)):((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1)), ] <- diag(q_emiss[q] - 1)
subj_data[[s]]$emiss_int_mle[[q]][i,] <- rep(0, q_emiss[q] - 1)
}
# if this is first iteration, use MLE for current values RW metropolis sampler
if (iter == 2){
emiss_c_int[[i]][[q]][s,] <- emiss_out$par
}
}
}
# Sample pouplaton values for gamma and conditional probabilities using Gibbs sampler -----------
# gamma_mu0_n and gamma_mu_int_bar are matrices, with the number of rows equal to the number of covariates, and ncol equal to number of intercepts estimated
gamma_mu0_n <- solve(t(xx[[1]]) %*% xx[[1]] + gamma_K0) %*% (t(xx[[1]]) %*% gamma_c_int[[i]] + gamma_K0 %*% gamma_mu0[[i]])
gamma_V_n <- gamma_V + t(gamma_c_int[[i]] - xx[[1]] %*% gamma_mu0_n) %*% (gamma_c_int[[i]] - xx[[1]] %*% gamma_mu0_n) + t(gamma_mu0_n - gamma_mu0[[i]]) %*% gamma_K0 %*% (gamma_mu0_n - gamma_mu0[[i]])
gamma_V_int[[i]] <- solve(rwish(S = solve(gamma_V_n), v = gamma_nu + n_subj))
gamma_mu_int_bar[[i]] <- gamma_mu0_n + solve(chol(t(xx[[1]]) %*% xx[[1]] + gamma_K0)) %*% matrix(rnorm((m - 1) * nx[1]), nrow = nx[1]) %*% t(solve(chol(solve(gamma_V_int[[i]]))))
gamma_exp_int <- matrix(exp(c(0, gamma_mu_int_bar[[i]][1,] )), nrow = 1)
gamma_mu_prob_bar[[i]] <- gamma_exp_int / as.vector(gamma_exp_int %*% c(rep(1,(m))))
for(q in 1:n_dep){
emiss_mu0_n <- solve(t(xx[[1 + q]]) %*% xx[[1 + q]] + emiss_K0[[q]]) %*% (t(xx[[1 + q]]) %*% emiss_c_int[[i]][[q]] + emiss_K0[[q]] %*% emiss_mu0[[q]][[i]])
emiss_V_n <- emiss_V[[q]] + t(emiss_c_int[[i]][[q]] - xx[[1 + q]] %*% emiss_mu0_n) %*% (emiss_c_int[[i]][[q]] - xx[[1 + q]] %*% emiss_mu0_n) + t(emiss_mu0_n - emiss_mu0[[q]][[i]]) %*% emiss_K0[[q]] %*% (emiss_mu0_n - emiss_mu0[[q]][[i]])
emiss_V_int[[i]][[q]] <- solve(rwish(S = solve(emiss_V_n), v = emiss_nu[[q]] + n_subj))
emiss_mu_int_bar[[i]][[q]] <- emiss_mu0_n + solve(chol(t(xx[[1 + q]]) %*% xx[[1 + q]] + emiss_K0[[q]])) %*% matrix(rnorm((q_emiss[q] - 1) * nx[1 + q]), nrow = nx[1 + q]) %*% t(solve(chol(solve(emiss_V_int[[i]][[q]]))))
emiss_exp_int <- matrix(exp(c(0, emiss_mu_int_bar[[i]][[q]][1, ])), nrow = 1)
emiss_mu_prob_bar[[i]][[q]] <- emiss_exp_int / as.vector(emiss_exp_int %*% c(rep(1, (q_emiss[q]))))
}
# Sample subject values for gamma and conditional probabilities using RW Metropolis sampler -----------
for (s in 1:n_subj){
gamma_candcov_comb <- chol2inv(chol(subj_data[[s]]$gamma_mhess[(1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1)), ] + chol2inv(chol(gamma_V_int[[i]]))))
gamma_RWout <- mnl_RW_once(int1 = gamma_c_int[[i]][s,], Obs = trans[[s]][[i]], n_cat = m, mu_int_bar1 = c(t(gamma_mu_int_bar[[i]]) %*% xx[[1]][s,]), V_int1 = gamma_V_int[[i]], scalar = gamma_scalar, candcov1 = gamma_candcov_comb)
gamma[[s]][i,] <- PD_subj[[s]][iter, c((sum(m * q_emiss) + 1 + (i - 1) * m):(sum(m * q_emiss) + (i - 1) * m + m))] <- gamma_RWout$prob
gamma_naccept[s, i] <- gamma_naccept[s, i] + gamma_RWout$accept
gamma_c_int[[i]][s,] <- gamma_RWout$draw_int
gamma_int_subj[[s]][iter, (1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1))] <- gamma_c_int[[i]][s,]
start <- c(0, q_emiss * m)
for(q in 1:n_dep){
emiss_candcov_comb <- chol2inv(chol(subj_data[[s]]$emiss_mhess[[q]][(1 + (i - 1) * (q_emiss[q] - 1)):((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1)), ] + chol2inv(chol(emiss_V_int[[i]][[q]]))))
emiss_RWout <- mnl_RW_once(int1 = emiss_c_int[[i]][[q]][s,], Obs = cond_y[[s]][[i]][[q]], n_cat = q_emiss[q], mu_int_bar1 = c(t(emiss_mu_int_bar[[i]][[q]]) %*% xx[[1 + q]][s,]), V_int1 = emiss_V_int[[i]][[q]], scalar = emiss_scalar[[q]], candcov1 = emiss_candcov_comb)
emiss[[s]][[q]][i,] <- PD_subj[[s]][iter, (sum(start[1:q]) + 1 + (i - 1) * q_emiss[q]):(sum(start[1:q]) + (i - 1) * q_emiss[q] + q_emiss[q])] <- emiss_RWout$prob
emiss_naccept[[q]][s, i] <- emiss_naccept[[q]][s, i] + emiss_RWout$accept
emiss_c_int[[i]][[q]][s,] <- emiss_RWout$draw_int
emiss_int_subj[[s]][[q]][iter, (1 + (i - 1) * (q_emiss[q] - 1)) : ((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1))] <- emiss_c_int[[i]][[q]][s,]
}
if(i == m){
delta[[s]] <- solve(t(diag(m) - gamma[[s]] + 1), rep(1, m))
}
}
}
# End of 1 MCMC iteration, save output values --------
gamma_int_bar[iter, ] <- unlist(lapply(gamma_mu_int_bar, "[",1,))
if(nx[1] > 1){
gamma_cov_bar[iter, ] <- unlist(lapply(gamma_mu_int_bar, "[",-1,))
}
gamma_prob_bar[iter,] <- unlist(gamma_mu_prob_bar)
for(q in 1:n_dep){
emiss_int_bar[[q]][iter, ] <- as.vector(unlist(lapply(
lapply(emiss_mu_int_bar, "[[", q), "[",1,)
))
if(nx[1+q] > 1){
emiss_cov_bar[[q]][iter, ] <- as.vector(unlist(lapply(
lapply(emiss_mu_int_bar, "[[", q), "[",-1,)
))
}
emiss_prob_bar[[q]][iter,] <- as.vector(unlist(sapply(emiss_mu_prob_bar, "[[", q)))
}
if(is.whole(iter/10) & print_iter == TRUE){
if(iter == 10){
cat("Iteration:", "\n")
}
cat(c(iter), "\n")
}
}
# End of function, return output values --------
ctime = proc.time()[3]
message(paste("Total time elapsed (hh:mm:ss):", hms(ctime-itime)))
if(return_path == TRUE){
out <- list(input = list(m = m, n_dep = n_dep, q_emiss = q_emiss, J = J,
burn_in = burn_in, n_subj = n_subj, n_vary = n_vary, dep_labels = dep_labels),
PD_subj = PD_subj, gamma_int_subj = gamma_int_subj, emiss_int_subj = emiss_int_subj,
gamma_int_bar = gamma_int_bar, gamma_cov_bar = gamma_cov_bar, emiss_int_bar = emiss_int_bar,
emiss_cov_bar = emiss_cov_bar, gamma_prob_bar = gamma_prob_bar,
emiss_prob_bar = emiss_prob_bar, gamma_naccept = gamma_naccept, emiss_naccept = emiss_naccept,
sample_path = sample_path)
} else {
out <- list(input = list(m = m, n_dep = n_dep, q_emiss = q_emiss, J = J,
burn_in = burn_in, n_subj = n_subj, n_vary = n_vary, dep_labels = dep_labels),
PD_subj = PD_subj, gamma_int_subj = gamma_int_subj, emiss_int_subj = emiss_int_subj,
gamma_int_bar = gamma_int_bar, gamma_cov_bar = gamma_cov_bar, emiss_int_bar = emiss_int_bar,
emiss_cov_bar = emiss_cov_bar, gamma_prob_bar = gamma_prob_bar,
emiss_prob_bar = emiss_prob_bar, gamma_naccept = gamma_naccept, emiss_naccept = emiss_naccept)
}
class(out) <- append(class(out), "mHMM")
return(out)
}
Try the mHMMbayes package in your browser
Any scripts or data that you put into this service are public.
mHMMbayes documentation built on Oct. 30, 2019, 5:05 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions mHMM
Documented in mHMM
#' Multilevel hidden Markov model using Bayesian estimation
#'
#' \code{mHMM} fits a multilevel (also known as mixed or random effects) hidden
#' Markov model (HMM) to intense longitudinal data with categorical 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 \link{mHMMbayes}. See
#' \code{vignette("tutorial-mhmm")} for an introduction to multilevel hidden
#' Markov models and the package, and see \code{vignette("estimation-mhmm")} for
#' an overview of the used estimation algorithms.
#'
#' Covariates specified in \code{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).
#'
#' If covariates are specified and the user wants to set the values for the
#' parameters of the hyper-prior distributions manually, the specification of
#' the elements in the arguments of the hyper-prior parameter values of the
#' normal distribution on the means change as follows: the number of rows in the
#' matrices \code{gamma_mu0} and \code{emiss_mu0} are equal to 1 + the number of
#' covariates used to predict the transition probability matrix for
#' \code{gamma_mu0} and the emission distribution for \code{emiss_mu0} (i.e.,
#' the first row correspond to the hyper-prior mean values of the intercepts,
#' the subsequent rows correspond to the hyper-prior mean values of the
#' regression coefficients connected to each of the covariates), and
#' \code{gamma_K0} and \code{emiss_K0} are now a matrix with the number of
#' hypothetical prior subjects on which the vectors of means (i.e., the rows in
#' \code{gamma_mu0} or \code{emiss_mu0}) are based on the diagonal, and
#' off-diagonal elements equal to 0. Note that the hyper-prior parameter values
#' of the inverse Wishart distribution on the covariance matrix remains
#' unchanged, as the estimates of the regression coefficients for the covariates
#' are fixed over subjects.
#'
#' @param s_data A matrix containing the observations to be modelled, where the
#' rows represent the observations over time. In \code{s_data}, 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 subsequent
#' columns contain the dependent variable(s). Note that the dependent
#' variables have to be numeric, i.e., they cannot be a (set of) factor
#' variable(s). The total number of rows are equal to the sum over the number
#' of observations of each subject, and the number of columns are equal to the
#' number of dependent variables (\code{n_dep}) + 1. The number of
#' observations can vary over subjects.
#' @param gen List containing the following elements denoting the general model
#' properties:
#' \itemize{\item{\code{m}: numeric vector with length 1 denoting the number
#' of hidden states}
#' \item{\code{n_dep}: numeric vector with length 1 denoting the
#' number of dependent variables}
#' \item{\code{q_emiss}: numeric vector with length \code{n_dep} denoting the
#' number of observed categories for the categorical emission distribution
#' for each of the dependent variables.}}
#' @param xx An optional list of (level 2) covariates to predict the transition
#' matrix and/or the emission probabilities. Level 2 covariate(s) means that
#' there is one observation per subject of each covariate. The first element
#' in the list \code{xx} is used to predict the transition matrix. Subsequent
#' elements in the list are used to predict the emission distribution of (each
#' of) the dependent variable(s). Each element in the list is a matrix, with
#' the number of rows equal to the number of subjects. The first column of
#' each matrix represents the intercept, that is, a column only consisting of
#' ones. Subsequent columns correspond to covariates used to predict the
#' transition matrix / emission distribution. See \emph{Details} for more
#' information on the use of covariates.
#'
#' If \code{xx} is omitted completely, \code{xx} defaults to \code{NULL},
#' resembling no covariates. Specific elements in the list can also be left
#' empty (i.e., set to \code{NULL}) to signify that either the transition
#' probability matrix or a specific emission distribution is not predicted by
#' covariates.
#' @param start_val List containing the start values for the transition
#' probability matrix gamma and the emission distribution(s). The first
#' element of the list contains a \code{m} by \code{m} matrix with the start
#' values for gamma. The subsequent elements contain \code{m} by
#' \code{q_emiss[k]} matrices for the start values for each of the \code{k}
#' emission distribution(s). Note that \code{start_val} should not contain
#' nested lists (i.e., lists within lists).
#' @param mcmc List of Markov chain Monte Carlo (MCMC) arguments, containing the
#' following elements:
#' \itemize{\item{\code{J}: numeric vector with length 1 denoting the number
#' of iterations of the MCMC algorithm}
#' \item{\code{burn_in}: numeric vector with length 1 denoting the
#' burn-in period for the MCMC algorithm.}}
#' @param return_path A logical scalar. Should the sampled state sequence
#' obtained at each iteration and for each subject be returned by the function
#' (\code{sample_path = TRUE}) or not (\code{sample_path = FALSE}). Note that
#' the sampled state sequence is quite a large object, hence the default
#' setting is \code{sample_path = FALSE}. Can be used for local decoding
#' purposes.
#' @param print_iter A logical scaler. Should the function print the progress of
#' the algorithm by returning the current iteration number at every 10th
#' iteration (\code{print_iter = TRUE}) or not (\code{print_iter = FALSE}).
#' Defaults to \code{print_iter = FALSE}.
#' @param gamma_hyp_prior An optional list containing user specified parameters
#' of the hyper-prior distribution on the multivariate normal distribution
#' of the intercepts (and regression coefficients given that covariates are
#' used) of the multinomial regression model of the transition probability
#' matrix gamma. The hyper-prior of the mean intercepts is a multivariate
#' Normal distribution, the hyper-prior of the covariance matrix between the
#' set of (state specific) intercepts is an Inverse Wishart distribution.
#'
#' Hence, the list \code{gamma_hyp_prior} contains the following elements:
#' \itemize{\item{\code{gamma_mu0}: a list containing m matrices; one matrix
#' for each row of the transition probability matrix gamma. Each matrix
#' contains the hypothesized mean values of the intercepts. Hence, each matrix
#' consists of one row (when not including covariates in the model) and
#' \code{m} - 1 columns}
#' \item{\code{gamma_K0}: numeric vector with length 1 denoting the number of
#' hypothetical prior subjects on which the vector of means \code{gamma_mu0} is
#' based}
#' \item{\code{gamma_nu}: numeric vector with length 1 denoting the degrees of
#' freedom of the Inverse Wishart distribution}
#' \item{\code{gamma_V}: matrix of \code{m} - 1 by \code{m} - 1 containing the
#' hypothesized variance-covariance matrix between the set of intercepts.}}
#' Note that \code{gamma_K0}, \code{gamma_nu} and \code{gamma_V} are assumed
#' equal over the states. The mean values of the intercepts (and regression
#' coefficients of the covariates) denoted by \code{gamma_mu0} are allowed to
#' vary over the states.
#'
#' The default values for the hyper-prior on gamma are: all elements of the
#' matrices contained in \code{gamma_mu0} set to 0, \code{gamma_K0} set to 1,
#' \code{gamma_nu} set to 3 + m - 1, and the diagonal of \code{gamma_V} (i.e.,
#' the variance) set to 3 + m - 1 and the off-diagonal elements (i.e., the
#' covariance) set to 0.
#'
#' See \emph{Details} below if covariates are used for changes in the settings
#' of the arguments of \code{gamma_hyp_prior}.
#' @param emiss_hyp_prior An optional list containing user specified parameters
#' of the hyper-prior distribution on the multivariate normal distribution
#' of the intercepts (and regression coefficients given that covariates are
#' used) of the multinomial regression model of the emission distribution.
#' The hyper-prior for the mean intercepts is a multivariate
#' Normal distribution, the hyper-prior for the covariance matrix between the
#' set of (state specific) intercepts is an Inverse Wishart distribution.
#'
#' Hence, the list \code{emiss_hyp_prior} contains the following elements:
#' \itemize{\item{\code{emiss_mu0}: a list of lists: \code{emiss_mu0} contains
#' \code{ndep} lists, i.e., one list for each dependent variable \code{k}. Each
#' of these lists contains m matrices; one matrix for each set of emission
#' probabilities within a state. The matrices contain the hypothesized mean
#' values of the intercepts. Hence, each matrix consists of one row (when not
#' including covariates in the model) and \code{q_emiss[k]} - 1 columns}
#' \item{\code{emiss_K0}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables, where each element is a numeric vector
#' with length 1 denoting the number of hypothetical prior subjects on which
#' the vector of means \code{emiss_mu0} is based}
#' \item{\code{emiss_nu}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables, where each element is a numeric vector
#' with length 1 denoting the degrees of freedom of the Inverse Wishart
#' distribution}
#' \item{\code{emiss_V}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables \code{k}, where each element is a matrix
#' of \code{q_emiss[k]} - 1 by \code{q_emiss[k]} - 1 containing the
#' hypothesized variance-covariance matrix between the set of intercepts.}}
#' Note that \code{emiss_K0}, \code{emiss_nu} and \code{emiss_V} are assumed
#' equal over the states. The mean values of the intercepts (and regression
#' coefficients of the covariates) denoted by \code{emiss_mu0} are allowed to
#' vary over the states.
#'
#' The default values for the hyper-prior on the emission distribution(s) are:
#' all elements of the matrices contained in \code{emiss_K0} set to 0,
#' \code{emiss_K0} set to 1, \code{emiss_nu} set to 3 + \code{q_emiss[k]} - 1,
#' and the diagonal of \code{gamma_V} (i.e., the variance) set to 3 +
#' \code{q_emiss[k]} - 1 and the off-diagonal elements (i.e., the covariance)
#' set to 0.
#'
#' See \emph{Details} below if covariates are used for changes in the settings
#' of the arguments of \code{emiss_hyp_prior}.
#' @param gamma_sampler An optional list containing user specified settings for
#' the proposal distribution of the random walk (RW) Metropolis sampler for
#' the subject level parameter estimates of the intercepts modeling the
#' transition probability matrix. The list \code{gamma_sampler} contains the
#' following elements:
#' \itemize{\item{\code{gamma_int_mle0}: a numeric vector with length \code{m}
#' - 1 denoting the start values for the maximum likelihood estimates of the
#' intercepts for the transition probability matrix gamma, based on the pooled
#' data (data over all subjects)}
#' \item{\code{gamma_scalar}: a numeric vector with length 1 denoting the scale
#' factor \code{s}. That is, The scale of the proposal distribution is composed
#' of a covariance matrix Sigma, which is then tuned by multiplying it by a
#' scaling factor \code{s}^2}
#' \item{\code{gamma_w}: a numeric vector with length 1 denoting the weight for
#' the overall log likelihood (i.e., log likelihood based on the pooled data
#' over all subjects) in the fractional likelihood.}}
#' Default settings are: all elements in \code{gamma_int_mle0} set to 0,
#' \code{gamma_scalar} set to 2.93 / sqrt(\code{m} - 1), and \code{gamma_w} set to
#' 0.1. See the section \emph{Scaling the proposal distribution of the RW
#' Metropolis sampler} in \code{vignette("estimation-mhmm")} for details.
#' @param emiss_sampler An optional list containing user specified settings for
#' the proposal distribution of the random walk (RW) Metropolis sampler for
#' the subject level parameter estimates of the intercepts modeling the
#' emission distributions of the dependent variables \code{k}. The list
#' \code{emiss_sampler} contains the following elements:
#' \itemize{\item{\code{emiss_int_mle0}: a list containing \code{ndep} elements
#' corresponding to each of the dependent variables \code{k}, where each
#' element is a a numeric vector with length \code{q_emiss[k]} - 1 denoting the
#' start values for the maximum likelihood estimates of the intercepts for
#' the emission distribution, based on the pooled data (data over all
#' subjects)}
#' \item{\code{emiss_scalar}: a list containing \code{ndep} elements
#' corresponding to each of the dependent variables, where each element is a
#' numeric vector with length 1 denoting the scale factor \code{s}. That is,
#' The scale of the proposal distribution is composed of a covariance matrix
#' Sigma, which is then tuned by multiplying it by a scaling factor \code{s}^2}
#' \item{\code{emiss_w}: a list containing \code{ndep} elements corresponding
#' to each of the dependent variables, where each element is a numeric vector
#' with length 1 denoting the weight for the overall log likelihood (i.e., log
#' likelihood based on the pooled data over all subjects) in the fractional
#' likelihood.}}
#' Default settings are: all elements in \code{emiss_int_mle0} set to 0,
#' \code{emiss_scalar} set to 2.93 / sqrt(\code{q_emiss[k]} - 1), and
#' \code{emiss_w} set to 0.1. See the section \emph{Scaling the proposal
#' distribution of the RW Metropolis sampler} in
#' \code{vignette("estimation-mhmm")} for details.
#'
#' @return \code{mHMM} returns an object of class \code{mHMM}, which has
#' \code{print} and \code{summary} methods to see the results.
#' The object contains the following components:
#' \describe{
#' \item{\code{PD_subj}}{A list containing one matrix per subject with the
#' subject level parameter estimates and the log likelihood 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 (parameter) estimates of subsequently the emission probabilities, the
#' transition probabilities and the log likelihood.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{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.}
#' \item{\code{input}}{Overview of used input specifications: the number of
#' states \code{m}, the number of used dependent variables \code{n_dep}, the
#' number of output categories for each of the dependent variables
#' \code{q_emiss}, the number of iterations \code{J} and the specified burn in
#' period \code{burn_in} of the hybrid Metropolis within Gibbs sampler, the
#' number of subjects \code{n_subj}, the observation length for each subject
#' \code{n_vary}, and the column names of the dependent variables
#' \code{dep_labels}.}
#' \item{\code{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 \code{return_path = TRUE}. }
#' }
#'
#' @seealso \code{\link{sim_mHMM}} for simulating multilevel hidden Markov data,
#' \code{\link{vit_mHMM}} for obtaining the most likely hidden state sequence
#' for each subject using the Viterbi algorithm, \code{\link{obtain_gamma}}
#' and \code{\link{obtain_emiss}} for obtaining the transition or emission
#' distribution probabilities of a fitted model at the group or subject level,
#' and \code{\link{plot.mHMM}} for plotting the posterior densities of a
#' fitted model.
#'
#' @references
#' \insertRef{rabiner1989}{mHMMbayes}
#'
#' \insertRef{scott2002}{mHMMbayes}
#'
#' \insertRef{altman2007}{mHMMbayes}
#'
#' \insertRef{rossi2012}{mHMMbayes}
#'
#' \insertRef{zucchini2017}{mHMMbayes}
#'
#' @examples
#' ###### Example on package example data
#' \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. Here, the covariate (standardized CDI
#' # change) predicts the emission distribution for each of the 4 dependent
#' # variables:
#'
#' n_subj <- 10
#' xx <- rep(list(matrix(1, ncol = 1, nrow = n_subj)), (n_dep + 1))
#' for(i in 2:(n_dep + 1)){
#' xx[[i]] <- cbind(xx[[i]], nonverbal_cov$std_CDI_change)
#' }
#' 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 simulated data
#' # Simulate data for 10 subjects with each 100 observations:
#' n_t <- 100
#' n <- 10
#' m <- 2
#' q_emiss <- 3
#' gamma <- matrix(c(0.8, 0.2,
#' 0.3, 0.7), ncol = m, byrow = TRUE)
#' emiss_distr <- 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, m = m, 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 = list(gamma, emiss_distr),
#' mcmc = list(J = 11, burn_in = 5))
#'
#'
#' @export
#'
#'
mHMM <- function(s_data, gen, xx = NULL, start_val, mcmc, return_path = FALSE, print_iter = FALSE,
gamma_hyp_prior = NULL, emiss_hyp_prior = NULL,
gamma_sampler = NULL, emiss_sampler = NULL){
# Initialize data -----------------------------------
# dependent variable(s), sample size, dimensions gamma and conditional distribuiton
n_dep <- gen$n_dep
dep_labels <- colnames(s_data[,2:(n_dep+1)])
id <- unique(s_data[,1])
n_subj <- length(id)
subj_data <- rep(list(NULL), n_subj)
if(sum(sapply(s_data, is.factor)) > 0 ){
stop("Your data contains factorial variables, which cannot be used as input in the function mHMM. All variables have to be numerical.")
}
for(s in 1:n_subj){
subj_data[[s]]$y <- as.matrix(s_data[s_data[,1] == id[s],][,-1], ncol = n_dep)
}
ypooled <- n <- NULL
n_vary <- numeric(n_subj)
m <- gen$m
q_emiss <- gen$q_emiss
emiss_int_mle <- rep(list(NULL), n_dep)
emiss_mhess <- rep(list(NULL), n_dep)
for(q in 1:n_dep){
emiss_int_mle[[q]] <- matrix(, m, (q_emiss[q] - 1))
emiss_mhess[[q]] <- matrix(, (q_emiss[q] - 1) * m, (q_emiss[q] - 1))
}
for(s in 1:n_subj){
ypooled <- rbind(ypooled, subj_data[[s]]$y)
n <- dim(subj_data[[s]]$y)[1]
n_vary[s] <- n
subj_data[[s]] <- c(subj_data[[s]], n = n, list(gamma_converge = numeric(m), gamma_int_mle = matrix(, m, (m - 1)),
gamma_mhess = matrix(, (m - 1) * m, (m - 1)), emiss_converge =
rep(list(numeric(m)), n_dep), emiss_int_mle = emiss_int_mle, emiss_mhess = emiss_mhess))
}
n_total <- dim(ypooled)[1]
# covariates
n_dep1 <- 1 + n_dep
nx <- numeric(n_dep1)
if (is.null(xx)){
xx <- rep(list(matrix(1, ncol = 1, nrow = n_subj)), n_dep1)
nx[] <- 1
} else {
if(!is.list(xx) | length(xx) != n_dep1){
stop("If xx is specified, xx should be a list, with the number of elements equal to the number of dependent variables + 1")
}
for(i in 1:n_dep1){
if (is.null(xx[[i]])){
xx[[i]] <- matrix(1, ncol = 1, nrow = n_subj)
nx[i] <- 1
} else {
nx[i] <- ncol(xx[[i]])
if (sum(xx[[i]][,1] != 1)){
stop("If xx is specified, the first column in each element of xx has to represent the intercept. That is, a column that only consists of the value 1")
}
if(nx[i] > 1){
for(j in 2:nx[i]){
if(is.factor(xx[[i]][,j])){
stop("Factors currently cannot be used as covariates, see help file for alternatives")
}
if((length(unique(xx[[i]][,j])) == 2) & (sum(xx[[i]][,j] != 0 & xx[[i]][,j] !=1) > 0)){
stop("Dichotomous covariates in xx need to be coded as 0 / 1 variables. That is, only conisting of the values 0 and 1")
}
if(length(unique(xx[[i]][,j])) > 2){
xx[[i]][,j] <- xx[[i]][,j] - mean(xx[[i]][,j])
}
}
}
}
}
}
# Initialize mcmc argumetns
J <- mcmc$J
burn_in <- mcmc$burn_in
# Initalize priors and hyper priors --------------------------------
# Initialize gamma sampler
if(is.null(gamma_sampler)) {
gamma_int_mle0 <- rep(0, m - 1)
gamma_scalar <- 2.93 / sqrt(m - 1)
gamma_w <- .1
} else {
gamma_int_mle0 <- gamma_sampler$gamma_int_mle0
gamma_scalar <- gamma_sampler$gamma_scalar
gamma_w <- gamma_sampler$gamma_w
}
# Initialize emiss sampler
if(is.null(emiss_sampler)){
emiss_int_mle0 <- rep(list(NULL), n_dep)
emiss_scalar <- rep(list(NULL), n_dep)
for(q in 1:n_dep){
emiss_int_mle0[[q]] <- rep(0, q_emiss[q] - 1)
emiss_scalar[[q]] <- 2.93 / sqrt(q_emiss[q] - 1)
}
emiss_w <- .1
} else {
emiss_int_mle0 <- emiss_sampler$emiss_int_mle0
emiss_scalar <- emiss_sampler$emiss_scalar
emiss_w <- emiss_sampler$emiss_w
}
# Initialize Gamma hyper prior
if(is.null(gamma_hyp_prior)){
gamma_mu0 <- rep(list(matrix(0,nrow = nx[1], ncol = m - 1)), m)
gamma_K0 <- diag(1, nx[1])
gamma_nu <- 3 + m - 1
gamma_V <- gamma_nu * diag(m - 1)
} else {
###### BUILD in a warning / check if gamma_mu0 is a matrix when given, with nrows equal to the number of covariates
gamma_mu0 <- gamma_hyp_prior$gamma_mu0
gamma_K0 <- gamma_hyp_prior$gamma_K0
gamma_nu <- gamma_hyp_prior$gamma_nu
gamma_V <- gamma_hyp_prior$gamma_V
}
# Initialize Pr hyper prior
# emiss_mu0: for each dependent variable, emiss_mu0 is a list, with one element for each state.
# Each element is a matrix, with number of rows equal to the number of covariates (with the intercept being one cov),
# and the number of columns equal to q_emiss[q] - 1.
emiss_mu0 <- rep(list(vector("list", m)), n_dep)
emiss_nu <- rep(list(NULL), n_dep)
emiss_V <- rep(list(NULL), n_dep)
emiss_K0 <- rep(list(NULL), n_dep)
if(is.null(emiss_hyp_prior)){
for(q in 1:n_dep){
for(i in 1:m){
emiss_mu0[[q]][[i]] <- matrix(0, ncol = q_emiss[q] - 1, nrow = nx[1 + q])
}
emiss_nu[[q]] <- 3 + q_emiss[q] - 1
emiss_V[[q]] <- emiss_nu[[q]] * diag(q_emiss[q] - 1)
emiss_K0[[q]] <- diag(1, nx[1 + q])
}
} else {
for(q in 1:n_dep){
# emiss_hyp_prior[[q]]$emiss_mu0 has to contain a list with lenght equal to m, and each list contains matrix with number of rows equal to number of covariates for that dep. var.
# stil build in a CHECK, with warning / stop / switch to default prior
emiss_mu0[[q]] <- emiss_hyp_prior[[q]]$emiss_mu0
emiss_nu[[q]] <- emiss_hyp_prior[[q]]$emiss_nu
emiss_V[[q]] <- emiss_hyp_prior[[q]]$emiss_V
emiss_K0[[q]] <- diag(emiss_hyp_prior$emiss_K0, nx[1 + q])
}
}
# Define objects used to store data in mcmc algorithm, not returned ----------------------------
# overall
c <- llk <- numeric(1)
sample_path <- lapply(n_vary, dif_matrix, cols = J)
trans <- rep(list(vector("list", m)), n_subj)
# gamma
gamma_mle_pooled <-gamma_int_mle_pooled <- gamma_mhess_pooled <- gamma_pooled_ll <- vector("list", m)
gamma_c_int <- rep(list(matrix(, n_subj, (m-1))), m)
gamma_mu_int_bar <- gamma_V_int <- vector("list", m)
gamma_mu_prob_bar <- rep(list(numeric(m)), m)
gamma_naccept <- matrix(0, n_subj, m)
# emiss
cond_y <- lapply(rep(n_dep, n_subj), nested_list, m = m)
emiss_mle_pooled <- emiss_int_mle_pooled <- emiss_mhess_pooled <- emiss_pooled_ll <- rep(list(vector("list", n_dep)), m)
emiss_c_int <- rep(list(lapply(q_emiss - 1, dif_matrix, rows = n_subj)), m)
emiss_mu_int_bar <- emiss_V_int <- rep(list(vector("list", n_dep)), m)
emiss_mu_prob_bar <- rep(list(lapply(q_emiss, dif_vector)), m)
emiss_naccept <- rep(list(matrix(0, n_subj, m)), n_dep)
# Define objects that are returned from mcmc algorithm ----------------------------
# Define object for subject specific posterior density, put start values on first row
if(length(start_val) != n_dep + 1){
stop("The number of elements in the list start_val should be equal to 1 + the number of dependent variables,
and should not contain nested lists (i.e., lists within lists)")
}
PD <- matrix(, nrow = J, ncol = sum(m * q_emiss) + m * m + 1)
PD_emiss_names <- paste("q", 1, "_emiss", rep(1:q_emiss[1], m), "_S", rep(1:m, each = q_emiss[1]), sep = "")
if(n_dep > 1){
for(q in 2:n_dep){
PD_emiss_names <- c(PD_emiss_names, paste("q", q, "_emiss", rep(1:q_emiss[q], m), "_S", rep(1:m, each = q_emiss[q]), sep = ""))
}
}
colnames(PD) <- c(PD_emiss_names, paste("S", rep(1:m, each = m), "toS", rep(1:m, m), sep = ""), "LL")
PD[1, ((sum(m * q_emiss) + 1)) :((sum(m * q_emiss) + m * m))] <- unlist(sapply(start_val, t))[1:(m*m)]
PD[1, 1:((sum(m * q_emiss)))] <- unlist(sapply(start_val, t))[(m*m + 1): (m*m + sum(m * q_emiss))]
PD_subj <- rep(list(PD), n_subj)
# Define object for population posterior density (probabilities and regression coefficients parameterization )
gamma_prob_bar <- matrix(, nrow = J, ncol = (m * m))
colnames(gamma_prob_bar) <- paste("S", rep(1:m, each = m), "toS", rep(1:m, m), sep = "")
gamma_prob_bar[1,] <- PD[1,(sum(m*q_emiss) + 1):(sum(m * q_emiss) + m * m)]
emiss_prob_bar <- lapply(q_emiss * m, dif_matrix, rows = J)
names(emiss_prob_bar) <- dep_labels
for(q in 1:n_dep){
colnames(emiss_prob_bar[[q]]) <- paste("Emiss", rep(1:q_emiss[q], m), "_S", rep(1:m, each = q_emiss[q]), sep = "")
start <- c(0, q_emiss * m)
emiss_prob_bar[[q]][1,] <- PD[1,(sum(start[1:q]) + 1):(sum(start[1:q]) + (m * q_emiss[q]))]
}
gamma_int_bar <- matrix(, nrow = J, ncol = ((m-1) * m))
colnames(gamma_int_bar) <- paste("int_S", rep(1:m, each = m-1), "toS", rep(2:m, m), sep = "")
gamma_int_bar[1,] <- as.vector(prob_to_int(matrix(gamma_prob_bar[1,], byrow = TRUE, ncol = m, nrow = m)))
if(nx[1] > 1){
gamma_cov_bar <- matrix(, nrow = J, ncol = ((m-1) * m) * (nx[1] - 1))
colnames(gamma_cov_bar) <- paste( paste("cov", 1 : (nx[1] - 1), "_", sep = ""), "S", rep(1:m, each = (m-1) * (nx[1] - 1)), "toS", rep(2:m, m * (nx[1] - 1)), sep = "")
gamma_cov_bar[1,] <- 0
} else{
gamma_cov_bar <- "No covariates where used to predict the transition probability matrix"
}
emiss_int_bar <- lapply((q_emiss-1) * m, dif_matrix, rows = J)
names(emiss_int_bar) <- dep_labels
for(q in 1:n_dep){
colnames(emiss_int_bar[[q]]) <- paste("int_Emiss", rep(2:q_emiss[q], m), "_S", rep(1:m, each = q_emiss[q] - 1), sep = "")
emiss_int_bar[[q]][1,] <- as.vector(prob_to_int(matrix(emiss_prob_bar[[q]][1,], byrow = TRUE, ncol = q_emiss[q], nrow = m)))
}
if(sum(nx[-1]) > n_dep){
emiss_cov_bar <- lapply((q_emiss-1) * m * (nx[-1] - 1 ), dif_matrix, rows = J)
names(emiss_cov_bar) <- dep_labels
for(q in 1:n_dep){
if(nx[1 + q] > 1){
colnames(emiss_cov_bar[[q]]) <- paste( paste("cov", 1 : (nx[1 + q] - 1), "_", sep = ""), "emiss", rep(2:q_emiss[q], m * (nx[1 + q] - 1)), "_S", rep(1:m, each = (q_emiss[q] - 1) * (nx[1 + q] - 1)), sep = "")
emiss_cov_bar[[q]][1,] <- 0
} else {
emiss_cov_bar[[q]] <- "No covariates where used to predict the emission probabilities for this outcome"
}
}
} else{
emiss_cov_bar <- "No covariates where used to predict the emission probabilities"
}
# Define object for subject specific posterior density (regression coefficients parameterization )
gamma_int_subj <- rep(list(gamma_int_bar), n_subj)
emiss_int_subj <- rep(list(emiss_int_bar), n_subj)
# Put starting values in place for fist run forward algorithm
emiss_sep <- vector("list", n_dep)
for(q in 1:n_dep){
start <- c(0, q_emiss * m)
emiss_sep[[q]] <- matrix(PD[1,(sum(start[1:q]) + 1):(sum(start[1:q]) + (m * q_emiss[q]))], byrow = TRUE, ncol = q_emiss[q], nrow = m)
}
emiss <- rep(list(emiss_sep), n_subj)
gamma <- rep(list(matrix(PD[1,(sum(m*q_emiss) + 1):(sum(m * q_emiss) + m * m)], byrow = TRUE, ncol = m)), n_subj)
delta <- rep(list(solve(t(diag(m) - gamma[[1]] + 1), rep(1, m))), n_subj)
# Start analysis --------------------------------------------
# Run the MCMC algorithm
itime <- proc.time()[3]
for (iter in 2 : J){
# For each subject, obtain sampled state sequence with subject individual parameters ----------
for(s in 1:n_subj){
# Run forward algorithm, obtain subject specific forward proababilities and log likelihood
forward <- cat_Mult_HMM_fw(x = subj_data[[s]]$y, m = m, emiss = emiss[[s]], gamma = gamma[[s]], n_dep = n_dep, delta=NULL)
alpha <- forward$forward_p
c <- max(forward$la[, subj_data[[s]]$n])
llk <- c + log(sum(exp(forward$la[, subj_data[[s]]$n] - c)))
PD_subj[[s]][iter, sum(m * q_emiss) + m * m + 1] <- llk
# Using the forward probabilites, sample the state sequence in a backward manner.
# In addition, saves state transitions in trans, and conditional observations within states in cond_y
trans[[s]] <- vector("list", m)
sample_path[[s]][n_vary[[s]], iter] <- sample(1:m, 1, prob = c(alpha[, n_vary[[s]]]))
for(t in (subj_data[[s]]$n - 1):1){
sample_path[[s]][t,iter] <- sample(1:m, 1, prob = (alpha[, t] * gamma[[s]][,sample_path[[s]][t + 1, iter]]))
trans[[s]][[sample_path[[s]][t,iter]]] <- c(trans[[s]][[sample_path[[s]][t, iter]]], sample_path[[s]][t + 1, iter])
}
for (i in 1:m){
trans[[s]][[i]] <- c(trans[[s]][[i]], 1:m)
trans[[s]][[i]] <- rev(trans[[s]][[i]])
for(q in 1:n_dep){
cond_y[[s]][[i]][[q]] <- c(subj_data[[s]]$y[sample_path[[s]][, iter] == i, q], 1:q_emiss[q])
}
}
}
# The remainder of the mcmc algorithm is state specific
for(i in 1:m){
# Obtain MLE of the covariance matrices and log likelihood of gamma and emiss at subject and population level -----------------
# used to scale the propasal distribution of the RW Metropolis sampler
# population level, transition matrix
trans_pooled <- factor(c(unlist(sapply(trans, "[[", i)), c(1:m)))
gamma_mle_pooled[[i]] <- optim(gamma_int_mle0, llmnl_int, Obs = trans_pooled, n_cat = m, method = "BFGS", hessian = TRUE, control = list(fnscale = -1))
gamma_int_mle_pooled[[i]] <- gamma_mle_pooled[[i]]$par
gamma_pooled_ll[[i]] <- gamma_mle_pooled[[i]]$value
# population level, conditional probabilities, seperate for each dependent variable
for(q in 1:n_dep){
cond_y_pooled <- numeric()
### MOET OOK ECHT BETER KUNNEN, eerst # cond_y_pooled <- unlist(sapply(cond_y, "[[", m))
for(s in 1:n_subj){
cond_y_pooled <- c(cond_y_pooled, cond_y[[s]][[i]][[q]])
}
emiss_mle_pooled[[i]][[q]] <- optim(emiss_int_mle0[[q]], llmnl_int, Obs = c(cond_y_pooled, c(1:q_emiss[q])), n_cat = q_emiss[q],
method = "BFGS", hessian = TRUE, control = list(fnscale = -1))
emiss_int_mle_pooled[[i]][[q]] <- emiss_mle_pooled[[i]][[q]]$par
emiss_pooled_ll[[i]][[q]] <- emiss_mle_pooled[[i]][[q]]$value
}
# subject level, transition matrix
for (s in 1:n_subj){
wgt <- subj_data[[s]]$n / n_total
gamma_out <- optim(gamma_int_mle_pooled[[i]], llmnl_int_frac, Obs = c(trans[[s]][[i]], c(1:m)), n_cat = m,
pooled_likel = gamma_pooled_ll[[i]], w = gamma_w, wgt = wgt, method="BFGS", hessian = TRUE, control = list(fnscale = -1))
if(gamma_out$convergence == 0){
subj_data[[s]]$gamma_converge[i] <- 1
subj_data[[s]]$gamma_mhess[(1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1)), ] <- mnlHess_int(int = gamma_out$par, Obs = c(trans[[s]][[i]], c(1:m)), n_cat = m)
subj_data[[s]]$gamma_int_mle[i,] <- gamma_out$par
} else {
subj_data[[s]]$gamma_converge[i] <- 0
subj_data[[s]]$gamma_mhess[(1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1)),] <- diag(m-1)
subj_data[[s]]$gamma_int_mle[i,] <- rep(0, m - 1)
}
# if this is first iteration, use MLE for current values RW metropolis sampler
if (iter == 2){
gamma_c_int[[i]][s,] <- gamma_out$par
}
# subject level, conditional probabilities, seperate for each dependent variable
for(q in 1:n_dep){
emiss_out <- optim(emiss_int_mle_pooled[[i]][[q]], llmnl_int_frac, Obs = c(cond_y[[s]][[i]][[q]], c(1:q_emiss[q])),
n_cat = q_emiss[q], pooled_likel = emiss_pooled_ll[[i]][[q]], w = emiss_w, wgt = wgt, method = "BFGS", hessian = TRUE, control = list(fnscale = -1))
if(emiss_out$convergence == 0){
subj_data[[s]]$emiss_converge[[q]][i] <- 1
subj_data[[s]]$emiss_mhess[[q]][(1 + (i - 1) * (q_emiss[q] - 1)):((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1)), ] <- mnlHess_int(int = emiss_out$par, Obs = c(cond_y[[s]][[i]][[q]], c(1:q_emiss[q])), n_cat = q_emiss[q])
subj_data[[s]]$emiss_int_mle[[q]][i,] <- emiss_out$par
} else {
subj_data[[s]]$emiss_converge[[q]][i] <- 0
subj_data[[s]]$emiss_mhess[[q]][(1 + (i - 1) * (q_emiss[q] - 1)):((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1)), ] <- diag(q_emiss[q] - 1)
subj_data[[s]]$emiss_int_mle[[q]][i,] <- rep(0, q_emiss[q] - 1)
}
# if this is first iteration, use MLE for current values RW metropolis sampler
if (iter == 2){
emiss_c_int[[i]][[q]][s,] <- emiss_out$par
}
}
}
# Sample pouplaton values for gamma and conditional probabilities using Gibbs sampler -----------
# gamma_mu0_n and gamma_mu_int_bar are matrices, with the number of rows equal to the number of covariates, and ncol equal to number of intercepts estimated
gamma_mu0_n <- solve(t(xx[[1]]) %*% xx[[1]] + gamma_K0) %*% (t(xx[[1]]) %*% gamma_c_int[[i]] + gamma_K0 %*% gamma_mu0[[i]])
gamma_V_n <- gamma_V + t(gamma_c_int[[i]] - xx[[1]] %*% gamma_mu0_n) %*% (gamma_c_int[[i]] - xx[[1]] %*% gamma_mu0_n) + t(gamma_mu0_n - gamma_mu0[[i]]) %*% gamma_K0 %*% (gamma_mu0_n - gamma_mu0[[i]])
gamma_V_int[[i]] <- solve(rwish(S = solve(gamma_V_n), v = gamma_nu + n_subj))
gamma_mu_int_bar[[i]] <- gamma_mu0_n + solve(chol(t(xx[[1]]) %*% xx[[1]] + gamma_K0)) %*% matrix(rnorm((m - 1) * nx[1]), nrow = nx[1]) %*% t(solve(chol(solve(gamma_V_int[[i]]))))
gamma_exp_int <- matrix(exp(c(0, gamma_mu_int_bar[[i]][1,] )), nrow = 1)
gamma_mu_prob_bar[[i]] <- gamma_exp_int / as.vector(gamma_exp_int %*% c(rep(1,(m))))
for(q in 1:n_dep){
emiss_mu0_n <- solve(t(xx[[1 + q]]) %*% xx[[1 + q]] + emiss_K0[[q]]) %*% (t(xx[[1 + q]]) %*% emiss_c_int[[i]][[q]] + emiss_K0[[q]] %*% emiss_mu0[[q]][[i]])
emiss_V_n <- emiss_V[[q]] + t(emiss_c_int[[i]][[q]] - xx[[1 + q]] %*% emiss_mu0_n) %*% (emiss_c_int[[i]][[q]] - xx[[1 + q]] %*% emiss_mu0_n) + t(emiss_mu0_n - emiss_mu0[[q]][[i]]) %*% emiss_K0[[q]] %*% (emiss_mu0_n - emiss_mu0[[q]][[i]])
emiss_V_int[[i]][[q]] <- solve(rwish(S = solve(emiss_V_n), v = emiss_nu[[q]] + n_subj))
emiss_mu_int_bar[[i]][[q]] <- emiss_mu0_n + solve(chol(t(xx[[1 + q]]) %*% xx[[1 + q]] + emiss_K0[[q]])) %*% matrix(rnorm((q_emiss[q] - 1) * nx[1 + q]), nrow = nx[1 + q]) %*% t(solve(chol(solve(emiss_V_int[[i]][[q]]))))
emiss_exp_int <- matrix(exp(c(0, emiss_mu_int_bar[[i]][[q]][1, ])), nrow = 1)
emiss_mu_prob_bar[[i]][[q]] <- emiss_exp_int / as.vector(emiss_exp_int %*% c(rep(1, (q_emiss[q]))))
}
# Sample subject values for gamma and conditional probabilities using RW Metropolis sampler -----------
for (s in 1:n_subj){
gamma_candcov_comb <- chol2inv(chol(subj_data[[s]]$gamma_mhess[(1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1)), ] + chol2inv(chol(gamma_V_int[[i]]))))
gamma_RWout <- mnl_RW_once(int1 = gamma_c_int[[i]][s,], Obs = trans[[s]][[i]], n_cat = m, mu_int_bar1 = c(t(gamma_mu_int_bar[[i]]) %*% xx[[1]][s,]), V_int1 = gamma_V_int[[i]], scalar = gamma_scalar, candcov1 = gamma_candcov_comb)
gamma[[s]][i,] <- PD_subj[[s]][iter, c((sum(m * q_emiss) + 1 + (i - 1) * m):(sum(m * q_emiss) + (i - 1) * m + m))] <- gamma_RWout$prob
gamma_naccept[s, i] <- gamma_naccept[s, i] + gamma_RWout$accept
gamma_c_int[[i]][s,] <- gamma_RWout$draw_int
gamma_int_subj[[s]][iter, (1 + (i - 1) * (m - 1)):((m - 1) + (i - 1) * (m - 1))] <- gamma_c_int[[i]][s,]
start <- c(0, q_emiss * m)
for(q in 1:n_dep){
emiss_candcov_comb <- chol2inv(chol(subj_data[[s]]$emiss_mhess[[q]][(1 + (i - 1) * (q_emiss[q] - 1)):((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1)), ] + chol2inv(chol(emiss_V_int[[i]][[q]]))))
emiss_RWout <- mnl_RW_once(int1 = emiss_c_int[[i]][[q]][s,], Obs = cond_y[[s]][[i]][[q]], n_cat = q_emiss[q], mu_int_bar1 = c(t(emiss_mu_int_bar[[i]][[q]]) %*% xx[[1 + q]][s,]), V_int1 = emiss_V_int[[i]][[q]], scalar = emiss_scalar[[q]], candcov1 = emiss_candcov_comb)
emiss[[s]][[q]][i,] <- PD_subj[[s]][iter, (sum(start[1:q]) + 1 + (i - 1) * q_emiss[q]):(sum(start[1:q]) + (i - 1) * q_emiss[q] + q_emiss[q])] <- emiss_RWout$prob
emiss_naccept[[q]][s, i] <- emiss_naccept[[q]][s, i] + emiss_RWout$accept
emiss_c_int[[i]][[q]][s,] <- emiss_RWout$draw_int
emiss_int_subj[[s]][[q]][iter, (1 + (i - 1) * (q_emiss[q] - 1)) : ((q_emiss[q] - 1) + (i - 1) * (q_emiss[q] - 1))] <- emiss_c_int[[i]][[q]][s,]
}
if(i == m){
delta[[s]] <- solve(t(diag(m) - gamma[[s]] + 1), rep(1, m))
}
}
}
# End of 1 MCMC iteration, save output values --------
gamma_int_bar[iter, ] <- unlist(lapply(gamma_mu_int_bar, "[",1,))
if(nx[1] > 1){
gamma_cov_bar[iter, ] <- unlist(lapply(gamma_mu_int_bar, "[",-1,))
}
gamma_prob_bar[iter,] <- unlist(gamma_mu_prob_bar)
for(q in 1:n_dep){
emiss_int_bar[[q]][iter, ] <- as.vector(unlist(lapply(
lapply(emiss_mu_int_bar, "[[", q), "[",1,)
))
if(nx[1+q] > 1){
emiss_cov_bar[[q]][iter, ] <- as.vector(unlist(lapply(
lapply(emiss_mu_int_bar, "[[", q), "[",-1,)
))
}
emiss_prob_bar[[q]][iter,] <- as.vector(unlist(sapply(emiss_mu_prob_bar, "[[", q)))
}
if(is.whole(iter/10) & print_iter == TRUE){
if(iter == 10){
cat("Iteration:", "\n")
}
cat(c(iter), "\n")
}
}
# End of function, return output values --------
ctime = proc.time()[3]
message(paste("Total time elapsed (hh:mm:ss):", hms(ctime-itime)))
if(return_path == TRUE){
out <- list(input = list(m = m, n_dep = n_dep, q_emiss = q_emiss, J = J,
burn_in = burn_in, n_subj = n_subj, n_vary = n_vary, dep_labels = dep_labels),
PD_subj = PD_subj, gamma_int_subj = gamma_int_subj, emiss_int_subj = emiss_int_subj,
gamma_int_bar = gamma_int_bar, gamma_cov_bar = gamma_cov_bar, emiss_int_bar = emiss_int_bar,
emiss_cov_bar = emiss_cov_bar, gamma_prob_bar = gamma_prob_bar,
emiss_prob_bar = emiss_prob_bar, gamma_naccept = gamma_naccept, emiss_naccept = emiss_naccept,
sample_path = sample_path)
} else {
out <- list(input = list(m = m, n_dep = n_dep, q_emiss = q_emiss, J = J,
burn_in = burn_in, n_subj = n_subj, n_vary = n_vary, dep_labels = dep_labels),
PD_subj = PD_subj, gamma_int_subj = gamma_int_subj, emiss_int_subj = emiss_int_subj,
gamma_int_bar = gamma_int_bar, gamma_cov_bar = gamma_cov_bar, emiss_int_bar = emiss_int_bar,
emiss_cov_bar = emiss_cov_bar, gamma_prob_bar = gamma_prob_bar,
emiss_prob_bar = emiss_prob_bar, gamma_naccept = gamma_naccept, emiss_naccept = emiss_naccept)
}
class(out) <- append(class(out), "mHMM")
return(out)
}
Try the mHMMbayes package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.