使用 R 包中的“ace”函数和“所有速率不同”模型时返回的无意义可能性：“猿”

问题描述

我在使用“ace”函数的“所有速率不同”(ARD) 模型时遇到了一些问题。该函数似乎返回不同的错误/警告，具体取决于计算速率矩阵指数的指定方法。最值得注意的是，它对于 >1 或

library(ape)
library(phytools)

#function to make a perfectly asymmetrical tree
make.asym.tree <- function(n = 32,sum.of.edges){
    taxa <- paste("t",1:n,sep = "")
    taxa <- as.list(taxa)
    o <- n - 1
    bl <- sum.of.edges/((o^2+3*o)/2)
    tree <- paste(taxa[[1]],":",bl,sep = "") 
    for(i in 2:n){
        tree <- paste("(",tree,",taxa[[i]],(bl * (i-1)),")",sep = "")
        if(i != n){
            tree <- paste(tree,sep = "")
        }
    }
tree <- paste(tree,";",sep = "")
tree <- read.newick(text = tree)
return(tree)
}

#make asymmetrical tree
asym <- make.asym.tree(64,sum.of.edges = 14)

#get tip data - this data was originally simulated using 'rTraitdisc()'
tip_data <- c(t1 = 3,t2 = 3,t3 = 3,t4 = 3,t5 = 3,t6 = 3,t7 = 3,t8 = 3,t9 = 3,t10 = 3,t11 = 3,t12 = 3,t13 = 3,t14 = 3,t15 = 3,t16 = 3,t17 = 3,t18 = 3,t19 = 3,t20 = 4,t21 = 3,t22 = 3,t23 = 3,t24 = 1,t25 = 3,t26 = 3,t27 = 3,t28 = 3,t29 = 3,t30 = 3,t31 = 3,t32 = 3,t33 = 3,t34 = 2,t35 = 3,t36 = 3,t37 = 3,t38 = 3,t39 = 3,t40 = 3,t41 = 4,t42 = 3,t43 = 3,t44 = 3,t45 = 3,t46 = 3,t47 = 3,t48 = 3,t49 = 3,t50 = 1,t51 = 3,t52 = 3,t53 = 3,t54 = 3,t55 = 3,t56 = 3,t57 = 1,t58 = 3,t59 = 3,t60 = 3,t61 = 3,t62 = 2,t63 = 3,t64 = 3)
tip_data <- as.factor(tip_data)

您可以使用以下方法可视化提示状态：

#plot data
cols<-setNames(c("red","blue","cyan","orange"),levels(tip_data))
plottree(asym,fsize=0.7,ftype="i",lwd=1)
tiplabels(pie=to.matrix(tip_data[asym$tip.label],levels(tip_data)),piecol=cols,cex=0.3)

首先我尝试了默认设置：

ARD_results <- ace(tip_data,asym,model = "ARD",type = "discrete")

这会产生大量警告，说明：“在 nlminb(rep(ip,length.out = np),function(p) dev(p),... : 强制丢弃的虚部"

为每个状态返回的似然是“复杂”类，有些大于 1，或者小于 0。

ARD_results$lik.anc

我认为这可能是使用 'matexpo' 函数的问题，所以我遵循了帮助中的建议，并使用了 'expm' 函数。

ARD_results2 <- ace(tip_data,type = "discrete",use.expm = T,use.eigen = F)

这会产生不同的警告消息：“In sqrt(diag(solve(h))) : NaNs产生” 这一次，状态的似然是“数字”类，但有些仍然大于 1，或者小于 0。

ARD_results2$lik.anc

最后，我尝试将 'use.expm' 和 'use.eigen' 都设置为 false。

ARD_results3 <- ace(tip_data,use.expm = F,use.eigen = F)

这仍然会产生警告信息：“In sqrt(diag(solve(h))) : NaNs产生” 但是，状态的可能性现在属于“数字”类，并且应该在 0 到 1 之间。

ARD_results3$lik.anc

但我不确定这种方法是否有效。 'use.expm' 和 'use.eigen' 都设置为 false，如何计算速率矩阵的指数？

如果有人解决了将“use.expm”或“use.eigen”设置为 true 的问题，我将不胜感激。否则，如果有人能向我解释将这两个设置为 false 的后果，我将不胜感激。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）