Stochastic character mapping with precursor model

[eliotmiller]

Wondering if it's possible or would be desirable to run stochastic character mapping with a precursor model. It currently does not work for me.

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")
model <- Precur_res.corHMM$solution
simmap <- makeSimmap(tree=phy, data=data, model=model, rate.cat=2, nSim=1, nCores=1)

returns

Error in `[<-`(`*tmp*`, i, state_index, value = 1) : 
  subscript out of bounds

Am I possibly just doing this wrong?

[jboyko]

Yes, the data your inputing is not the same between the model fit and the stochastic map function. You'll also need to specify the root prior again when doing makeSimmap. The correct code would like this: simmap <- makeSimmap(tree=phy, data=Precur_Dat, model=model, rate.cat=2, nSim=1, nCores=1, root.p = "maddfitz") Best, James

[eliotmiller]

Thanks! This was a problem in my worked example here, I wasn't making the same mistake in the actual code. My problem in the actual code was failing to reference the $solution object specifically for the model argument. I mistook errors in both cases to mean it wasn't working with the precursor model, but they were two separate mistakes I'd made. Thanks again!