thej022214 / corHMM

Fits a generalized form of the covarion model that allows different transition rate classes on different portions of a phylogeny by treating rate classes as “hidden” states in a Markov process.
11 stars 13 forks source link

dropStateMatPars bug #32

Closed Victaphanta closed 2 years ago

Victaphanta commented 2 years ago

Hi James, I am now observing odd behavior with this new syntax for dropStateMatPars.

I essentially want to take this matrix and remove some of the rates (i.e. make N/A), specifically rates 3,4,13,14,17,18

This is the original matrix. recon01b[["index.mat"]] 1 2 3 4 5 1 NA 5 9 13 17 2 1 NA 10 14 18 3 2 6 NA 15 19 4 3 7 11 NA 20 5 4 8 12 16 NA

Using this syntax

rate.matbr<-dropStateMatPars(StateMat = recon01b$index.mat, Pars = c(3,4,13,14,17,18))

I get this rate.matbr 1 2 3 4 5 1 NA 4 8 3 3 2 1 NA 9 3 3 3 2 5 NA 12 14 4 3 6 10 NA 15 5 3 7 11 13 NA

Instead of changing to N/A they all get allocated “3”

Not sure what I am doing wrong.

Thanks for your help.

Adnan

jboyko commented 2 years ago

Hi Adnan,

Thanks for pointing this out! This was a bug that we had fixed, but apparently has now reappeared after a recent version change. I've pushed the changes to the newest github version, so if you install the newest version the functions will be working. But, I'll also write them out here if you'd rather implement them in your code manually.

rate_mat <- matrix(c(
  NA, 5, 9, 13, 17,
  1, NA, 10, 14, 18,
  2, 6, NA, 15, 19,
  3, 7, 11, NA, 20,
  4, 8, 12, 16, NA), 5, 5, TRUE)

dropStateMatPars <- function(StateMat, Pars){
  for(i in Pars){
    StateMat[StateMat==i] <- 0
  }
  StateMat[StateMat == 0] <- NA
  pars <- unique(na.omit(as.vector(StateMat)))
  for(i in 1:length(pars)){
    StateMat[StateMat == pars[i]] <- i
  }
  return(StateMat)
}

equateStateMatPars <- function(StateMat, ParsList){
  if(!class(ParsList) == "list"){
    ParsList <- list(ParsList)
  }
  for(i in 1:length(ParsList)){
    min_par_i <- min(ParsList[[i]])
    StateMat[as.vector(StateMat) %in% ParsList[[i]]] <- min_par_i
  }
  StateMat[StateMat == 0] <- NA
  pars <- unique(na.omit(as.vector(StateMat)))
  for(i in 1:length(pars)){
    StateMat[StateMat == pars[i]] <- i
  }
  return(StateMat)
}

new_mat <- dropStateMatPars(StateMat = rate_mat, Pars = c(3,4,13,14,17,18))

new_mat_2 <- equateStateMatPars(new_mat, list(c(11,9,10), c(1,2,3,4)))