# R/obtain_gamma.R In mHMMbayes: Multilevel Hidden Markov Models Using Bayesian Estimation

#### Documented in obtain_gamma

```#' Obtain the transition probabilities gamma for a fitted multilevel HMM
#'
#' \code{obtain_gamma} obtains the transition probability matrix for a
#' fitted multilevel hidden Markov model, for either the group level, i.e.,
#'   representing the average transition probability matrix over all subjects,
#'   or at the subject level, returning the transition probability matrices for
#'   each subject.
#'
#' @param object An object of class \code{mHMM}, generated by the function
#' @param level String specifying if the returned transition probability matrix
#'   gamma should be at the group level (\code{level = "group"}), i.e.,
#'   representing the average transition probability matrix 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_gamma} returns the object \code{est_gamma} of the class
#'   \code{mHMM_gamma}. This object can be directly plotted using the function
#'   \code{plot.mHMM_gamma()}, or simply \code{plot()}. Depending on the
#'   specification at the input variable \code{level}, \code{est_gamma} is
#'   either a matrix with the transition probabilities at the group level (if
#'   \code{level = "group"}), or a list of matrices (with the number of elements
#'   equal to the number of subjects analyzed, if \code{level = 'subject'}),
#'   where each matrix in the list represents a subject specific transition
#'   probability matrix.
#'
#' @seealso \code{\link{mHMM}} for fitting the multilevel hidden Markov
#' model, creating the object \code{mHMM}, and \code{\link{plot.mHMM_gamma}} for
#' plotting the obtained transition probabilities.
#'
#'
#'
#' @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), # vocalizing patient
#'                  matrix(c(0.1, 0.9,
#'                           0.1, 0.9), byrow = TRUE, nrow = m,
#'                         ncol = q_emiss), # looking patient
#'                  matrix(c(0.90, 0.05, 0.05,
#'                           0.05, 0.90, 0.05), byrow = TRUE,
#'                         nrow = m, ncol = q_emiss), # vocalizing therapist
#'                  matrix(c(0.1, 0.9,
#'                           0.1, 0.9), byrow = TRUE, nrow = m,
#'                         ncol = q_emiss)) # 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
#' obtain_gamma(out_2st, level = "group")
#' obtain_gamma(out_2st, level = "subject")
#'
#' }
#'
#' @export
#'
obtain_gamma <- 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
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
est <- matrix(, ncol = m, nrow = m)
colnames(est) <- paste("To state", 1:m)
rownames(est) <- paste("From state", 1:m)
if (level == "group"){
est[] <- matrix(round(apply(object\$gamma_prob_bar[((burn_in + 1): J),], 2, median),3),
byrow = TRUE, ncol = m, nrow = m)
est_gamma <- est
}
if (level == "subject"){
est_gamma <- rep(list(est), n_subj)
names(est_gamma) <- paste("Subject", 1:n_subj)
for(i in 1:n_subj){
est_gamma[[i]][] <- matrix(round(apply(object\$PD_subj[[i]][burn_in:J, (1 + sum(m * q_emiss)) : (sum(m * q_emiss) + m*m)], 2, median), 3),
byrow = TRUE, ncol = m, nrow = m)
}
}
class(est_gamma) <- append(class(est_gamma), "mHMM_gamma")
return(est_gamma)
}
```

## 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.