R/obtain_emiss.R

Defines functions obtain_emiss

Documented in obtain_emiss

#' Obtain the emission distribution probabilities for a fitted multilevel HMM
#'
#' \code{obtain_emiss} obtains the emission distribution probabilities (also
#' known as conditional probabilities) for a fitted multilevel hidden Markov
#' model, for either the group level, i.e., representing the average emission
#' distribution probabilities over all subjects, or at the subject level,
#' returning the emission distribution probabilities for each subject.
#'
#' @param object An object of class \code{mHMM}, generated by the function
#'   \code{\link{mHMM}}.
#' @param level String specifying if the returned emission distribution
#'   probabilities should be at the group level (\code{level = "group"}), i.e.,
#'   representing the average emission distribution probabilities over all
#'   subjects, or at the subject level (\code{level = "subject"}).
#' @param burn_in An integer which specifies the number of iterations to discard
#'   when obtaining the model parameter summary statistics. When left
#'   unspecified (\code{burn_in = NULL}), the burn in period specified when
#'   creating the \code{mHMM} object with the function \code{\link{mHMM}}
#'   will be used.
#'
#' @return \code{obtain_emiss} returns the object \code{est_emiss}. Depending on
#'   the specification at the input variable \code{level}, \code{est_emiss} is
#'   either a list of matrices with the emission distribution probabilities at
#'   the group level (if \code{level = "group"}) for each dependent variable, or
#'   a list of lists, where for each dependent variable a list is returned with
#'   the number of elements equal to the number of subjects analyzed, if
#'   \code{level = 'subject'}). In the latter scenario, each matrix in the lower
#'   level list represents the subject specific emission distribution
#'   probabilities for a specific dependent variable.
#'
#' @seealso \code{\link{mHMM}} for fitting the multilevel hidden Markov
#' model, creating the object \code{mHMM}.
#'
#'
#' @examples
#' ###### Example on package 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):
#' out_2st <- mHMM(s_data = nonverbal,
#'                 gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
#'                 start_val = c(list(start_TM), start_EM),
#'                 mcmc = list(J = 11, burn_in = 5))
#'
#' out_2st
#' summary(out_2st)
#'
#' # obtaining the emission probabilities at the group and subject level
#' obtain_emiss(out_2st, level = "group")
#' obtain_emiss(out_2st, level = "subject")
#'
#' }
#'
#' @export
#'
obtain_emiss <- function(object, level = "group", burn_in = NULL){
  if (!is.mHMM(object)){
    stop("The input object used should be from the class mHMM, obtained by using the function mHMM.")
  }
  if (level != "group" & level != "subject"){
    stop("The specification at the input variable -level- should be set to either group or subject")
  }
  input   <- object$input
  dep_labels <- input$dep_labels
  n_subj  <- input$n_subj
  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))
  }
  m       <- input$m
  q_emiss <- input$q_emiss
  n_dep   <- input$n_dep
  if (level == "group"){
    est <- lapply(q_emiss, dif_matrix, rows = m)
    for(j in 1:n_dep){
      colnames(est[[j]]) <- paste("Category", 1:q_emiss[j])
      rownames(est[[j]]) <- paste("State", 1:m)
    }
    names(est) <- dep_labels
    for(j in 1: n_dep){
      est[[j]][] <- matrix(round(apply(object$emiss_prob_bar[[j]][((burn_in + 1): J),], 2, median),3),
                      byrow = TRUE, ncol = q_emiss[j], nrow = m)
    }
    est_emiss <- est
  }
  if (level == "subject"){
    est_emiss <- vector("list", n_dep)
    names(est_emiss) <- dep_labels
    for(j in 1:n_dep){
      est <- matrix(,ncol = q_emiss[j], nrow = m)
      colnames(est) <- paste("Category", 1:q_emiss[j])
      rownames(est) <- paste("State", 1:m)
      est_emiss[[j]] <- rep(list(est), n_subj)
      names(est_emiss[[j]]) <- paste("Subject", 1:n_subj)
    }
    start <- c(0, q_emiss * m)
    for(i in 1:n_subj){
      for(j in 1:n_dep){
        est_emiss[[j]][[i]][] <- matrix(round(apply(object$PD_subj[[i]][burn_in:J, (sum(start[1:j]) + 1) : sum(start[1:(j+1)])], 2, median), 3),
                                   byrow = TRUE, ncol = q_emiss[j], nrow = m)
      }
    }
  }
  return(est_emiss)
}

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.