R/plot.mHMM.R

Defines functions plot.mHMM

Documented in plot.mHMM

#' Plotting the posterior densities for a fitted multilevel HMM
#'
#' \code{plot.mHMM} plots the posterior densities for a fitted multilevel hidden
#' Markov model for the group and subject level parameters simultaneously. The
#' plotted posterior densities are either for the transition probability matrix
#' gamma, or for the emission distribution probabilities.
#'
#' @param x Object of class \code{mHMM}, generated by the function
#'   \code{\link{mHMM}}.
#' @param component String specifying if the displayed posterior densities
#'   should be for the transition probability matrix gamma (\code{component =
#'   "gamma"}), or for the emission distribution probabilities (\code{component
#'   = "emiss"}). In case of the latter and the model is based on multiple
#'   dependent variables, the user has to indicate for which dependent variable
#'   the posterior densities have to be plotted, see \code{dep}.
#' @param dep Integer specifying for which dependent variable the posterior
#'   densities should be plotted. Only required if one wishes to plot the
#'   emission distribution probabilities and the model is based on multiple
#'   dependent variables.
#' @param col Vector of colors for the posterior density lines. If one is
#'   plotting the posterior densities for gamma, the vector has length \code{m}
#'   (i.e., number of hidden states). If one is plotting the posterior densities
#'   for the emission probabilities, the vector has length \code{q_emiss[k]}
#'   (i.e., the number of outcome categories for the dependent variable
#'   \code{k}).
#' @param cat_lab Optional vector of strings when plotting the posterior
#'   densities of the emission probabilities, denoting the labels of the
#'   categorical outcome values. Automatically generated when not provided.
#' @param dep_lab Optional string when plotting the posterior
#'   densities of the emission probabilities with length 1, denoting the label
#'   for the dependent variable plotted. Automatically obtained from the input
#'   object \code{x} when not specified.
#' @param lwd1 Positive number indicating the line width of the posterior
#'   density at the group level.
#' @param lwd2 Positive number indicating the line width of the posterior
#'   density at the subject level.
#' @param lty1 Positive number indicating the line type of the posterior
#'   density at the group level.
#' @param lty2 Positive number indicating the line type of the posterior
#'   density at the subject level.
#' @param burn_in An integer which specifies the number of iterations to discard
#'   when obtaining the model parameter summary statistics. When left
#'   unspecified, the burn in period specified at creating the \code{mHMM}
#'   object with the function \code{\link{mHMM}} will be used.
#' @param legend_cex A numerical value giving the amount by which plotting text
#'   and symbols in the legend should be magnified relative to the default.
#' @param ... Arguments to be passed to methods (see \code{\link[graphics]{par}})
#'
#' @return \code{plot.mHMM} returns a plot of the posterior densities. Depending
#'   on whether (\code{component = "gamma"}) or (\code{component = "emiss"}),
#'   the plotted posterior densities are either for the transition probability
#'   matrix gamma or for the emission distribution probabilities, respectively.
#'
#' @seealso \code{\link{mHMM}} for fitting the multilevel hidden Markov
#'   model, creating the object \code{mHMM}.
#'
#' @examples
#' ###### First run the function mHMM on the nonverbal 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))
#'
#' ## plot the posterior densities for gamma
#' plot(out_2st, component = "gamma")
#' }
#'
#' @export
#'
plot.mHMM <- function(x, component = "gamma", dep = 1, col, cat_lab,
                      dep_lab, lwd1 = 2, lwd2 = 1, lty1 = 1, lty2 = 3,
                      legend_cex, burn_in, ...){
  if (!is.mHMM(x)){
    stop("The input object x should be from the class mHMM, obtained with the function mHMM.")
  }
  if (component != "gamma" & component != "emiss"){
    stop("The input specified under component should be a string, restrectid to state either gamma or emiss.")
  }
  object <- x
  input   <- x$input
  n_subj  <- input$n_subj
  if (missing(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))
  }
  old_par <- graphics::par(no.readonly =TRUE)
  on.exit(graphics::par(old_par))
  m       <- input$m
  q_emiss <- input$q_emiss
  n_dep   <- input$n_dep

  if(component == "gamma"){
    if (missing(col)){
      state_col <- grDevices::rainbow(m)
    } else {
      state_col <- col
    }
    if(m > 3){
      graphics::par(mfrow = c(2,ceiling(m/2)), mar = c(4,2,3,1) + 0.1, mgp = c(2,1,0))
    } else {
      graphics::par(mfrow = c(1,m), mar = c(4,2,3,1) + 0.1, mgp = c(2,1,0))
    }
    for(i in 1:m){
      max <- 0
      for(j in 1:m){
        new <- max(stats::density(object$gamma_prob_bar[burn_in:J, m * (i-1) + j])$y)
        if(new > max){max <- new}
      }
      graphics::plot.default(x = 1, ylim = c(0, max), xlim = c(0,1), type = "n", cex = .8,  main =
             paste("From state", i, "to state ..."), yaxt = "n", ylab = "",
           xlab = "Transition probability", ...)
      graphics::title(ylab="Density", line=.5)
      for(j in 1:m){
        graphics::lines(stats::density(object$gamma_prob_bar[burn_in:J,m * (i-1) + j]),
              type = "l", col = state_col[j], lwd = lwd1, lty = lty1)
        for(s in 1:n_subj){
          graphics::lines(stats::density(object$PD_subj[[s]][burn_in:J,(sum(q_emiss * m) + m * (i-1) + j)]),
                type = "l", col = state_col[j], lwd = lwd2, lty = lty2)
        }
      }
      graphics::legend("topright", col = state_col, legend = paste("To state", 1:m),
             bty = 'n', lty = 1, lwd = 2, cex = .8)
    }
  } else if (component == "emiss"){
    if (missing(cat_lab)){
      cat_lab <- paste("Category", 1:q_emiss[dep])
    }
    if (missing(dep_lab)){
      dep_lab <- input$dep_labels[dep]
    }
    start <- c(0, q_emiss * m)
    start2 <- c(0, seq(from = (q_emiss[dep]-1) * 2, to = (q_emiss[dep]-1) * 2 * m, by = (q_emiss[dep]-1) * 2))
    if (missing(col)){
      cat_col <- grDevices::rainbow(q_emiss[dep])
    } else {
      cat_col <- col
    }
    if(m > 3){
      graphics::par(mfrow = c(2,ceiling(m/2)), mar = c(4,2,3,1) + 0.1, mgp = c(2,1,0))
    } else {
      graphics::par(mfrow = c(1,m), mar = c(4,2,3,1) + 0.1, mgp = c(2,1,0))
    }
    for(i in 1:m){
      # determining the scale of the y axis
      max <- 0
      for(q in 1:q_emiss[dep]){
        new <- max(stats::density(object$emiss_prob_bar[[dep]][burn_in:J,q_emiss[dep] * (i-1) + q])$y)
        if(new > max){max <- new}
      }
      # set plotting area
      graphics::plot.default(x = 1, ylim = c(0, max), xlim = c(0,1), type = "n",
           main = paste(dep_lab, ", state", i),
           yaxt = "n", ylab = "", xlab = "Conditional probability", ...)
      graphics::title(ylab="Density", line=.5)
      for(q in 1:q_emiss[dep]){
        # add density curve for population level posterior distribution
        graphics::lines(stats::density(object$emiss_prob_bar[[dep]][burn_in:J,q_emiss[dep] * (i-1) + q]),
              type = "l", col = cat_col[q], lwd = lwd1, lty = lty1)
        # add density curves for subject posterior distributions
        for(s in 1:10){
          graphics::lines(stats::density(object$PD_subj[[s]][burn_in:J,(sum(start[1:dep])
                                                     + (i-1)*q_emiss[dep] + q)]),
                type = "l", col = cat_col[q], lwd = lwd2, lty = lty2)
        }
      }
      graphics::legend("topright", col = cat_col, legend = cat_lab, bty = 'n', lty = 1, lwd = 2, cex = .7)
    }
  }
}

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.