R/assign.covergence.R

In compbiomed/ASSIGN: Adaptive Signature Selection and InteGratioN (ASSIGN)

#' Check the convergence of the MCMC chain
#'
#' The assign.convergence checks the convergence of the MCMC chain of the model
#' parameters generated by the Gibbs sampling algorithm.
#'
#' To compute the convergence of the gth gene in B, set whichGene=g,
#' whichSample=NA, whichPath=NA.
#'
#' To compute the convergence of the gth gene in the kth pathway within the
#' signature matrix (S), set whichGene=g, whichSample=NA, whichPath=NA.
#'
#' To compute the convergence of the kth pathway in the jth test sample within
#' the pathway activation matrix (A), set whichGene=NA, whichSample=n,
#' whichPath=k.
#'
#' @param test The list object returned from the assign.mcmc function. The list
#' components are the MCMC chains of the B, S, Delta, beta, gamma, and sigma.
#' @param burn_in The number of burn-in iterations. These iterations are
#' discarded when computing the posterior means of the model parameters. The
#' default is 0.
#' @param iter The number of total iterations. The default is 2000.
#' @param parameter A character string indicating which model parameter is to be
#' checked for convergence. This must be one of "B", "S", "Delta", "beta",
#' "kappa", "gamma", and "sigma".
#' @param whichGene A numerical value indicating which gene is to be checked
#' for convergence. The value has to be in the range between 1 and G.
#' @param whichSample A numerical value indicating which test sample is to be
#' checked for convergence. The value has to be in the range between 1 and N.
#' @param whichPath A numerical value indicating which pathway is to be checked
#' for convergence. The value has to be in the range between 1 and K.
#' @return The assign.convergence function returns the a vector of the
#' estimated values from each Gibbs sampling iteration of the model parameter
#' to be checked, and a trace plot of this parameter.
#' @author Ying Shen
#' @examples
#'
#' \dontrun{
#' # check the 10th gene in the 1st pathway for the convergence
#' trace.plot <- assign.convergence(test=mcmc.chain, burn_in=0, iter=2000,
#'                                  parameter="S", whichGene=10, whichSample=NA,
#'                                  whichPath=1)
#' }
#'
#' @export assign.convergence
assign.convergence <- function(test, burn_in=0, iter=2000,
                               parameter=c("B", "S", "Delta", "beta", "kappa",
                                           "gamma", "sigma"), whichGene,
                               whichSample, whichPath) {
  if (parameter == "B") {
    B_pos <- test$B_mcmc[burn_in:iter, whichGene]
    graphics::plot(B_pos, type = "l", xlab = "Iteration", ylab = "Posterior B",
                   main = "Convergence of posterior B")
    return(B_pos)
  }

  if (parameter == "S") {
    S_pos <- test$S_mcmc[burn_in:iter, whichGene, whichPath]
    graphics::plot(S_pos, type = "l", xlab = "Iteration",
                   ylab = paste("Posterior S of pathway ", whichPath, " gene ",
                              whichGene, sep = ""),
                   main = "Convergence of posterior S")
    return(S_pos)
  }

  if (parameter == "Delta") {
    Delta_pos <- test$Delta_mcmc[burn_in:iter, whichGene, whichPath]
    graphics::plot(Delta_pos, type = "l", xlab = "Iteration",
                   ylab = paste("Posterior Delta of pathway ", whichPath,
                              " gene ", whichGene, sep = ""),
                   main = "Convergence of posterior Delta")
    return(Delta_pos)
  }

  if (parameter == "beta") {
    beta_pos <- test$beta_mcmc[burn_in:iter, whichPath, whichSample]
    graphics::plot(beta_pos, type = "l", xlab = "Iteration",
                   ylab = paste("Posterior beta of sample ", whichSample,
                                " pathway ", whichPath, sep = ""),
                   main = "Convergence of posterior beta")
    return(beta_pos)
  }

  if (parameter == "kappa") {
    kappa_pos <- test$kappa_mcmc[burn_in:iter, whichPath, whichSample]
    graphics::plot(kappa_pos, type = "l", xlab = "Iteration",
                   ylab = paste("Posterior kappa of sample ", whichSample,
                                " pathway ", whichPath, sep = ""),
                   main = "Convergence of posterior kappa")
    return(kappa_pos)
  }

  if (parameter == "gamma") {
    gamma_pos <- test$gamma_mcmc[burn_in:iter, whichPath, whichSample]
    graphics::plot(gamma_pos, type = "l", xlab = "Iteration",
                   ylab = paste("Posterior gamma of sample ", whichSample,
                                " pathway ", whichPath, sep = ""),
                   main = "Convergence of posterior gamma")
    return(gamma_pos)
  }

  if (parameter == "sigma") {
    sigma_pos <- 1 / test$tau2_mcmc[burn_in:iter, whichGene]
    graphics::plot(sigma_pos, type = "l", xlab = "Iteration",
                   ylab = paste("Posterior sigma^2 of gene ",
                                whichGene, sep = ""),
                   main = "Convergence of posterior sigma^2")
    return(sigma_pos)
  }
}