R/hilda_diagnosis.R

In HiLDA: Conducting statistical inference on comparing the mutational exposures of mutational signatures by using hierarchical latent Dirichlet allocation

#' Output the maximum potential scale reduction statistic of all parameters
#' estimated
#'
#' @param jagsOutput the output jags file generated by the jags function from
#' the R2jags package.
#'
#' @return a number for the Rhat statistic.
#'
#' @importFrom methods is
#' 
#' @examples
#' 
#' inputFile <- system.file("extdata/hildaLocal.rdata", package="HiLDA")
#' hildaLocal <- readRDS(inputFile)
#' hildaRhat(hildaLocal)
#'
#' @export
#

hildaRhat <- function(jagsOutput) {
    if(is(jagsOutput) != "rjags") {
        stop("Not an output object from running the HiLDA tests.")
    }
    return(max(jagsOutput$BUGSoutput$summary[, "Rhat"]))
}


#' Extract the posterior distributions of the mean differences in muational
#' exposures
#'
#' @param jagsOutput the output jags file generated by the jags function from
#' the R2jags package.
#'
#' @importFrom methods is
#'
#' @examples
#' 
#' inputFile <- system.file("extdata/hildaLocal.rdata", package="HiLDA")
#' hildaLocal <- readRDS(inputFile)
#' hildaLocalResult(hildaLocal)
#'
#' @importFrom stringr str_detect
#' @return a data frame that contains the posterior distributions of difference.
#' @export
#'
#

hildaLocalResult <- function(jagsOutput) {
    if(is(jagsOutput) != "rjags") {
        stop("Not an output object from running the HiLDA tests.")
    }
    rownames <- rownames(jagsOutput$BUGSoutput$summary)
    beta.index <- stringr::str_detect(rownames, "beta")
    return(jagsOutput$BUGSoutput$summary[beta.index, ])
}

#' Compute the Bayes factor
#'
#' @param jagsOutput the output jags file generated by the jags function from
#' the R2jags package.
#' @param pM1 the probability of sampling the null (default: 0.5)
#' @return a number for the Bayes factor
#'
#' @importFrom methods is
#'
#' @examples
#' 
#' load(system.file("extdata/sample.rdata", package="HiLDA"))
#' hildaGlobal <- hildaTest(inputG=G, numSig=3, refGroup=1:4, nIter=1000,
#' localTest=TRUE)
#' hildaGlobalResult(hildaGlobal)
#'
#' @export
#

hildaGlobalResult <- function(jagsOutput, pM1=0.5) {
    if(is(jagsOutput) != "rjags") {
        stop("Not an output object from running the HiLDA tests.")
    }
    
    if(pM1 <= 0 | pM1 >= 1) {
        stop("Please input a fraction between 0 and 1.")
    }
    
    freq <- table(jagsOutput$BUGSoutput$sims.list$pM2)
    if (length(freq) == 1) {
        stop("It got stuck in the model ", as.numeric(names(freq)) + 1)
    }
    return(as.vector(freq)[2]/as.vector(freq)[1] * (pM1)/(1 - pM1))
}