R/vit_mHMM.R

In mHMMbayes: Multilevel Hidden Markov Models Using Bayesian Estimation

#' Obtain hidden state sequence for each subject using the Viterbi
#' algorithm
#'
#' \code{vit_mHMM} obtains the most likely state sequence (for each subject)
#' from an object of class \code{mHMM} (generated by the function
#' \code{mHMM()}), using (an extended version of) the Viterbi algorithm. This is
#' also known as global decoding. Optionally, the state probabilities
#' themselves at each point in time can also be obtained.
#'
#' Note that local decoding is also possible, by inferring the most frequent
#' state at each point in time for each subject from the sampled state path at
#' each iteration of the MCMC algorithm. This information is contained in the
#' output object \code{return_path} of the function \code{mHMM()}.
#'
#' @param object An object of class \code{mHMM}, generated
#'   by the function \code{\link{mHMM}}.
#' @inheritParams mHMM
#' @param burn_in The number of iterations to be discarded from the MCMC
#'   algorithm when inferring the transition probability matrix gamma and the
#'   emission distribution of (each of) the dependent variable(s) for each
#'   subject from \code{s_data}. If omitted, defaults to \code{NULL} and
#'   \code{burn_in} specified at the function \code{mHMM} will be used.
#' @param return_state_prob A logical scaler. Should the function, in addition
#'   to the most likely state sequence for each subject, also return the state
#'   probabilities at each point in time for each subject
#'   (\code{return_state_prob = TRUE}) or not (\code{return_state_prob =
#'   FALSE}). Defaults to \code{return_state_prob = FALSE}.
#'
#' @return The function \code{vit_mHMM} returns a matrix containing the most
#'   likely state at each point in time. The first column indicates subject id
#'   number. Hence, the id number is repeated over rows equal to the number of
#'   observations for that subject. The second column contains the most likely
#'   state at each point in time. If requested, the subsequent columns contain
#'   the state probabilities at each point in time for each subject.
#'
#' @examples
#' ###### Example on package example categorical data, see ?nonverbal
#' # First fit the multilevel HMM on the 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
#'
#' # Fit the multilevel HMM model:
#' # Note that for reasons of running time, J is set at a ridiculous low value.
#' # One would typically use a number of iterations J of at least 1000,
#' # and a burn_in of 200.
#' out_2st <- mHMM(s_data = nonverbal, gen = list(m = m, n_dep = n_dep,
#'                 q_emiss = q_emiss), start_val = c(list(start_TM), start_EM),
#'                 mcmc = list(J = 3, burn_in = 1))
#'
#' ###### obtain the most likely state sequence with the Viterbi algorithm
#' states1 <- vit_mHMM(s_data = nonverbal, object = out_2st)
#' }
#'
#' ###### Example on simulated multivariate continuous data
#' # simulating multivariate continuous data
#' n_t     <- 100
#' n       <- 10
#' m       <- 3
#' n_dep   <- 2
#'
#' gamma   <- matrix(c(0.8, 0.1, 0.1,
#'                     0.2, 0.7, 0.1,
#'                     0.2, 0.2, 0.6), ncol = m, byrow = TRUE)
#'
#' emiss_distr <- list(matrix(c( 50, 10,
#'                               100, 10,
#'                               150, 10), nrow = m, byrow = TRUE),
#'                     matrix(c(5, 2,
#'                              10, 5,
#'                              20, 3), nrow = m, byrow = TRUE))
#'
#' data_cont <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep), data_distr = 'continuous',
#'                   gamma = gamma, emiss_distr = emiss_distr, var_gamma = .1, var_emiss = c(.5, 0.01))
#'
#' # 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(100, m), rep(25, m)),
#'                         emiss_nu = list(1, 1),
#'                         emiss_a0 = list(rep(1, m), rep(1, m)),
#'                         emiss_b0 = list(rep(1, m), rep(1, m)))
#'
#' # Run the model on the simulated data:
#' # Note that for reasons of running time, J is set at a ridiculous low value.
#' # One would typically use a number of iterations J of at least 1000,
#' # and a burn_in of 200.
#' out_3st_cont_sim <- mHMM(s_data = data_cont$obs,
#'                     data_distr = "continuous",
#'                     gen = list(m = m, n_dep = n_dep),
#'                     start_val = c(list(gamma), emiss_distr),
#'                     emiss_hyp_prior = manual_prior_emiss,
#'                     mcmc = list(J = 11, burn_in = 5))
#'
#' ###### obtain the most likely state sequence with the Viterbi algorithm,
#' # including state probabilities at each time point.
#' states2 <- vit_mHMM(s_data = data_cont$obs, object = out_3st_cont_sim,
#'                     return_state_prob = TRUE)
#'
#'
#'
#' @seealso \code{\link{mHMM}} for analyzing multilevel hidden Markov data
#'   and obtaining the input needed for \code{vit_mHMM}, and
#'   \code{\link{sim_mHMM}} for simulating multilevel hidden Markov data.
#'
#' @references
#'  \insertRef{viterbi1967}{mHMMbayes}
#'
#'  \insertRef{rabiner1989}{mHMMbayes}
#'
#' @export
#'
#'
vit_mHMM <- function(object, s_data, burn_in = NULL, return_state_prob = FALSE){
  message("Please note that the output format is changed from wide to long format to facilitate aditionally returning the state probabilities, see the section 'Value' in the help file for more information.")
  if (!is.mHMM(object)){
    stop("The input object used should be from the class mHMM, obtained by using the function mHMM.")
  }
  if(sum(objects(object$PD_subj[[1]]) %in% "log_likl") != 1){
    stop("The input object is created using an earlier version of the mHMMbayes package. Please re-run the function mHMM with the current package version, or post-process the object using the earlier version of the package.")
  }
  input      <- object$input
  data_distr <- input$data_distr
  id         <- unique(s_data[,1])
  n_subj     <- length(id)
  if(length(object$PD_subj) != n_subj){
    stop("s_data used should be from the same subjects used for creating the object in mHMM.
         The number of subjects in the datasets are not the same.")
  }
  n_vary     <- table(s_data[,1])
  max_n      <- max(n_vary)
  state_seq  <- matrix(,ncol = n_subj, nrow = max_n)
  probs      <- vector(mode = "list", length = n_subj)
  n_dep      <- input$n_dep
  m          <- input$m
  if(is.null(burn_in)){
    burn_in  <- input$burn_in
  }
  J          <- input$J
  if (burn_in >= (J-1)){
    stop(paste("The specified burn in period should be at least 2 points smaller
               compared to the number of iterations J, J =", J))
  }
  if(data_distr == "categorical"){
    q_emiss    <- input$q_emiss
    int_est_emiss <- rep(list(lapply(q_emiss-1, dif_matrix, rows = m)), n_subj)
    est_emiss  <- rep(list(lapply(q_emiss, dif_matrix, rows = m)), n_subj)
    for(s in 1:n_subj){
      for(j in 1:n_dep){
        int_est_emiss[[s]][[j]][] <- matrix(apply(object$emiss_int_subj[[s]][[j]][burn_in:J, ], 2, median),
                                            byrow = TRUE, ncol = q_emiss[j]-1, nrow = m)
        est_emiss[[s]][[j]][] <- int_to_prob(int_est_emiss[[s]][[j]])
      }
    }
  } else if(data_distr == "continuous"){
    est_emiss  <- rep(list(rep(list(matrix(,nrow = m, ncol = 2)),n_dep)), n_subj)
    for(s in 1:n_subj){
      for(q in 1:n_dep){
        est_emiss[[s]][[q]][] <- matrix(round(c(apply(object$PD_subj[[s]]$cont_emiss[((burn_in + 1): J),((q-1) * m + 1):(q * m)], 2, median),
                                                apply(object$PD_subj[[s]]$cont_emiss[((burn_in + 1): J), (n_dep * m + (q-1) * m + 1):(n_dep * m + q * m)], 2, median)),3),
                                        ncol = 2, nrow = m)
      }
    }
  }

  est_gamma <- obtain_gamma(object, level = "subject")
  for(s in 1:n_subj){
    emiss <- est_emiss[[s]]
    gamma    <- est_gamma[[s]]
    if(data_distr == "categorical"){
      probs[[s]]    <- cat_mult_fw_r_to_cpp(x = as.matrix(s_data[s_data[,1] == id[s],][,-1], ncol = n_dep),
                                       m = m, emiss = emiss, gamma = gamma, n_dep = n_dep, delta=NULL)[[1]]
    } else if(data_distr == "continuous"){
      probs[[s]]    <- cont_mult_fw_r_to_cpp(x = as.matrix(s_data[s_data[,1] == id[s],][,-1], ncol = n_dep),
                            m = m, emiss = emiss, n_dep = n_dep, gamma = gamma)[[1]]
    }
    state_seq[1:n_vary[s], s] <- apply(probs[[s]], 2, which.max)
  }
  state_seq <- as.data.frame(state_seq)
  colnames(state_seq) <- id
  state_seq_stacked <- utils::stack(state_seq)[,c(2,1)]
  colnames(state_seq_stacked) <- c("subj", "state")
  if(return_state_prob == TRUE){
    probs <- sapply(probs, t, simplify = FALSE)
    probs <- do.call(rbind, probs)
    colnames(probs) <- paste0("pr_state_", 1:m)
    state_seq_stacked <- cbind(state_seq_stacked, probs)
  }
  return(state_seq = state_seq_stacked)
}