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

Precursor model no longer running #41

Closed eliotmiller closed 2 years ago

eliotmiller commented 2 years ago

Code that previously worked for me to run precursor models is no longer working. I'm not able to run the code in the vignette to see if the example works there, as the data file in the vignette directory appears corrupt. Suspect the example would fail too, possibly a good place to start. Am getting this error when I run the model.

Error in liks.down[root, ] <- liks.down[root, ] * root.p : number of items to replace is not a multiple of replacement length

jboyko commented 2 years ago

Hi,

The code still works, but the "yang" root.prior does not because the 2,R1 state is not allowed into the model. You can either specify a root.prior, use "maddfitz", or use "flat". Example code below.

require(corHMM)
require(MASS)
data(primates)
phy <- multi2di(primates[[1]])
data <- primates[[2]]

Precur_Dat <- data[,c(1,2)]

Precur_LegendAndMat <- getStateMat4Dat(Precur_Dat)
Precur_R1 <- Precur_LegendAndMat$rate.mat 
Precur_R1 <- dropStateMatPars(Precur_R1, c(1, 2)) 
Precur_R2 <- equateStateMatPars(Precur_LegendAndMat$rate.mat, c(1,2))
RateClassMat <- equateStateMatPars(getRateCatMat(2), c(1,2))

Precur_FullMat <- getFullMat(list(Precur_R1, Precur_R2), RateClassMat) 
Precur_FullMat[c(4, 2), c(2, 4)] <- 0
Precur_FullMat

Precur_res.corHMM <- corHMM(phy = phy, data = Precur_Dat, rate.cat = 2, rate.mat = Precur_FullMat, root.p = "maddfitz")

Best, James

eliotmiller commented 2 years ago

That does it! Thanks!