R/plot.mHMM_gamma.R

In mHMMbayes: Multilevel Hidden Markov Models Using Bayesian Estimation

#' Plotting the transition probabilities gamma for a fitted multilevel HMM
#'
#' \code{plot.mHMM_gamma} plots the transition probability matrix for a fitted
#' multilevel hidden Markov model, by means of an alluvial plot (also known as
#' Sankey diagram or riverplot) using the R package \code{alluvial}. The plotted
#' transition probability matrix either represents the probabilities at the
#' group level, i.e., representing the average transition probability matrix
#' over all subjects, or at the subject level. In case of the latter, the user
#' has to specify for which subject the transition probability matrix should be
#' plotted.
#'
#' @param x An object of class \code{mHMM_gamma}, generated by the function
#'   \code{\link{obtain_gamma}}.
#' @param subj_nr An integer specifying for which specific subject the
#'   transition probability matrix should be plotted. Only required if the input
#'   object represents the subject specific transition probability matrices.
#' @param cex An integer specifying scaling of fonts of category labels. When
#'   not specified, defaults to \code{cex = 0.8}.
#' @param col An optional vector with length \code{m} * \code{m} (i.e., where
#'   \code{m} denotes the number of hidden states) specifying the used colors in
#'   the alluvial plot.
#' @param hide An optional logical vector with  length \code{m} * \code{m}
#'   (i.e., where \code{m} denotes the number of hidden states) specifying
#'   whether particular stripes should be plotted. When not specified, omits
#'   the lines representing a value of exactly zero.
#' @param ... Arguments to be passed to alluvial (see
#'   \code{\link[alluvial]{alluvial}})
#'
#' @return \code{plot.mHMM_gamma} returns a plot of the transition probability
#'   matrix. Depending on whether the input object represents the transition
#'   probabilities at the group level or the subject specific transition
#'   probability matrices, the returned plot represents either the group
#'   transition probability matrix, or the transition probability matrix for a
#'   given subject, specified by \code{subject_nr}.
#'
#' @seealso \code{\link{mHMM}} for fitting the multilevel hidden Markov
#'   model, creating the object \code{mHMM}, and \code{\link{obtain_gamma}} to
#'   obtain the transition probabilities gamma for a fitted multilevel HMM,
#'   creating the object \code{mHMM_gamma}.
#'
#'
#' @examples
#' \donttest{
#' #' ###### Example on package data, see ?nonverbal
#' # specifying general model properties:
#' m <- 2
#' n_dep <- 4
#' q_emiss <- c(3, 2, 3, 2)
#'
#' # specifying starting values
#' start_TM <- diag(.8, m)
#' start_TM[lower.tri(start_TM) | upper.tri(start_TM)] <- .2
#' start_EM <- list(matrix(c(0.05, 0.90, 0.05,
#'                           0.90, 0.05, 0.05), byrow = TRUE,
#'                         nrow = m, ncol = q_emiss[1]), # vocalizing patient
#'                  matrix(c(0.1, 0.9,
#'                           0.1, 0.9), byrow = TRUE, nrow = m,
#'                         ncol = q_emiss[2]), # looking patient
#'                  matrix(c(0.90, 0.05, 0.05,
#'                           0.05, 0.90, 0.05), byrow = TRUE,
#'                         nrow = m, ncol = q_emiss[3]), # vocalizing therapist
#'                  matrix(c(0.1, 0.9,
#'                           0.1, 0.9), byrow = TRUE, nrow = m,
#'                         ncol = q_emiss[4])) # looking therapist
#'
#' # Run a model without covariate(s):
#' out_2st <- mHMM(s_data = nonverbal,
#'                 gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
#'                 start_val = c(list(start_TM), start_EM),
#'                 mcmc = list(J = 11, burn_in = 5))
#'
#' out_2st
#' summary(out_2st)
#'
#' # obtaining the transition probabilities at the group and subject level
#' est_gamma_group <- obtain_gamma(out_2st, level = "group")
#'
#' # plot the obtained transition probabilities
#' plot(est_gamma_group, col = rep(c("green", "blue"), each = m))
#'
#' }
#' \dontshow{
#' ###### Example on simulated data
#' # Simulate data for 10 subjects with each 100 observations:
#' n_t <- 100
#' n <- 10
#' m <- 2
#' n_dep <- 1
#' q_emiss <- 3
#' gamma <- matrix(c(0.8, 0.2,
#'                   0.3, 0.7), ncol = m, byrow = TRUE)
#' emiss_distr <- list(matrix(c(0.5, 0.5, 0.0,
#'                         0.1, 0.1, 0.8), nrow = m, ncol = q_emiss, byrow = TRUE))
#' data1 <- sim_mHMM(n_t = n_t, n = n, gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss), gamma = gamma,
#'                   emiss_distr = emiss_distr, var_gamma = .5, var_emiss = .5)
#'
#' # Specify remaining required analysis input (for the example, we use simulation
#' # input as starting values):
#' n_dep <- 1
#' q_emiss <- 3
#'
#' # Run the model on the simulated data:
#' out_2st_sim <- mHMM(s_data = data1$obs,
#'                  gen = list(m = m, n_dep = n_dep, q_emiss = q_emiss),
#'                  start_val = c(list(gamma), emiss_distr),
#'                  mcmc = list(J = 11, burn_in = 5))
#'
#' # obtaining the transition probabilities at the group and subject level
#' est_gamma_group_sim <- obtain_gamma(out_2st_sim, level = "group")
#'
#' # plot the obtained transition probabilities
#' plot(est_gamma_group_sim, col = rep(c("green", "blue"), each = m))
#'
#' }
#'
#' @export


plot.mHMM_gamma <- function(x, subj_nr = NULL, cex = 0.8, col, hide, ...){
  if (!requireNamespace("alluvial", quietly = TRUE)) {
    stop("Package \"alluvial\" needed for this function to work. Please install it.",
         call. = FALSE)
  }
  if (!requireNamespace("grDevices", quietly = TRUE)) {
    stop("Package \"grDevices\" needed for this function to work. Please install it.",
         call. = FALSE)
  }
  if (!is.mHMM_gamma(x)){
    stop("The input object x should be from the class mHMM_gamma, obtained with the function obtain_gamma.")
  }
  old_par <- graphics::par(no.readonly =TRUE)
  on.exit(graphics::par(old_par))
  if (is.list(x)){
    if (is.null(subj_nr)){
      stop("When the input object x represents the subject specific transition
           probability matrices, the subject for which the probabilities should
           be plotted needs to be specified with the input variable -subj_nr-.")
    }
    m <- dim(x[[subj_nr]])[1]
    From <- paste("State", rep(1:m, each = m))
    To <-  paste("State", rep(1:m, m))
    trans <- as.vector(t(x[[subj_nr]]))
    foo <- data.frame(From, To, trans)
    if(missing(col)){
      col <- c(rep(grDevices::rainbow(m), eac = m))
    }
    if (missing(hide)){
      hide <- foo$trans == 0
    }
    alluvial::alluvial(foo[,1:2], freq=foo$trans,
                       cex = cex,
                       col = col,
                       hide = hide, ...)
  } else {
    if(!is.null(subj_nr)){
      warning("The subject number can only be specified when plotting the subject level transition probabilities. Currently, the group level transition probabilities are plotted.")
    }
    m <- dim(x)[1]
    From <- paste("State", rep(1:m, each = m))
    To <-  paste("State", rep(1:m, m))
    trans <- as.vector(t(x))
    foo <- data.frame(From, To, trans)
    if(missing(col)){
      col <- c(rep(grDevices::rainbow(m), eac = m))
    }
    if (missing(hide)){
      hide <- foo$trans == 0
    }
    alluvial::alluvial(foo[,1:2], freq=foo$trans,
                     cex = cex,
                     col = col,
                     hide =  hide, ...)
  }
}