#' Adjusts reversible mutation matrix
#'
#' The Metropolis - Hastings conversions may give a mutation matrix \code{balancedMutmat}
#' with too small expected mutation rate (gamma). The \code{balancedMutmat} matrix is adjusted
#' to have expected mutation rate a equal to that of the original mutation matrix \code{mutmat}.
#'
#' @param mutmat Original, non balanced, mutation matrix.
#' @param balancedMutmat Balanced, mutation matrix.
#' @param method Character. MH, PM or BA conversions.
#' @param afreq A vector with allele frequencies.
#' of the same length as the size of \code{mutmat}.
#' @param check Logical. Checks if \code{balancedMutmat} is reversible.
#'
#' @return Adjusted mutation matrix.
#'
#' @details If \code{balancedMutmat == NULL}, \code{mutmat} is first balanced.
#' The adjusted balanced matrix is
#'
#' \code{alpha * balancedMutmat + (1-alpha) * I}
#'
#' where
#'
#' \code{alpha} is the ratio of the (expected mutation) rates of the original matrix,
#' \code{mutmat} to the balanced version \code{balancedMutmat} and \code{I} is the identity matrix.
#'
#' @seealso [makeReversible()]
#'
#' @author Thore Egeland.
#'
#' @export
#'
#' @examples
#' library(pedmut)
#' afreq = c(0.1, 0.3, 0.4, 0.2)
#' names(afreq) = 1:4
#' mutmat = mutationMatrix("onestep", rate = 0.02, alleles = 1:4)
#' adj = adjustReversible(mutmat, balancedMutmat = NULL,
#' method = "BA", afreq = afreq, check = TRUE)
#' attr(mutmat, "rate") - attr(adj, "rate")
adjustReversible = function(mutmat, balancedMutmat, method = "MH",
afreq = NULL, check = TRUE){
if (is.null(afreq))
afreq = attr(mutmat, "afreq")
if(is.null(balancedMutmat))
balancedMutmat = makeReversible(mutmat, method = method,
afreq = afreq, check = check)
if(check){
if(!isReversible(balancedMutmat, afreq))
stop("Second argument needs to a balanced mutation matrix")
}
RM = expectedMutationRate(mutmat, afreq)
RP = expectedMutationRate(balancedMutmat, afreq)
if(is.null(RM) | is.null(RP) | (RP == 0))
stop("Non-negative (expected mutation) rates of matrices are needed")
alpha = RM/RP
newM = alpha * balancedMutmat + (1-alpha) * diag(rep(1,length(afreq)))
pedmut:::newMutationMatrix(newM, afreq = afreq, model = "custom",
rate = expectedMutationRate(newM, afreq))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.