R/assign.covergence.R

In wevanjohnson/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 convergency of the gth gene in B, set whichGene=g,
#' whichSample=NA, whichPath=NA.
#' 
#' To compute the convergency of the gth gene in the kth pathway within the
#' signature matrix (S), set whichGene=g, whichSample=NA, whichPath=NA.
#' 
#' To compute the convergency 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 charater string indicating which model parameter is to be
#' checked for convergency. 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 convergency
#' 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]
    plot(B_pos, type="l", xlab="Iteration", ylab="Posterior B", main="Convergency of posterior B")
    return(B_pos)
  }
  
  if (parameter == "S"){
    S_pos <- test$S_mcmc[burn_in:iter, whichGene, whichPath]
    plot(S_pos, type="l", xlab="Iteration", ylab=paste("Posterior S of pathway ", whichPath, " gene ", whichGene, sep=""), main="Convergency of posterior S")
    return(S_pos)
  }
  
  if (parameter == "Delta"){
    Delta_pos <- test$Delta_mcmc[burn_in:iter, whichGene, whichPath]
    plot(Delta_pos, type="l", xlab="Iteration", ylab=paste("Posterior Delta of pathway ", whichPath, " gene ", whichGene, sep=""), main="Convergency of posterior Delta")
    return(Delta_pos)
  }
  
  
  if (parameter == "beta"){
    beta_pos <- test$beta_mcmc[burn_in:iter,whichPath,whichSample]
    plot(beta_pos, type="l", xlab="Iteration", ylab=paste("Posterior beta of sample ", whichSample, " pathway ", whichPath, sep=""), main="Convergency of posterior beta")
    return(beta_pos)
  }
  
  if (parameter == "kappa"){
    kappa_pos <- test$kappa_mcmc[burn_in:iter,whichPath,whichSample]
    plot(kappa_pos, type="l", xlab="Iteration", ylab=paste("Posterior kappa of sample ", whichSample, " pathway ", whichPath, sep=""), main="Convergency of posterior kappa")
    return(kappa_pos)
  }
  
  if (parameter == "gamma"){
    gamma_pos <- test$gamma_mcmc[burn_in:iter,whichPath,whichSample]
    plot(gamma_pos, type="l", xlab="Iteration", ylab=paste("Posterior gamma of sample ", whichSample, " pathway ", whichPath, sep=""), main="Convergency of posterior gamma")
    return(gamma_pos)
  }
  
  if (parameter == "sigma"){
    sigma_pos <- 1/test$tau2_mcmc[burn_in:iter, whichGene]
    plot(sigma_pos, type="l",xlab="Iteration", ylab=paste("Posterior sigma^2 of gene ", whichGene, sep=""), main="Convergency of posterior sigma^2")
    return(sigma_pos)
  }
}