R/adjustReversible.R

In thoree/mut2: Reversible mutation matrices and pairwise likelihood ratios

#' 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))
  }