使用每个物种的多个条目的系统发育模型答案

【问题标题】：Phylogenetic model using multiple entries for each species使用每个物种的多个条目的系统发育模型
【发布时间】：2018-07-02 08:50:08
【问题描述】：

我对系统发育回归模型比较陌生。过去，当我的树中每个物种只有 1 个条目时，我使用 PGLS。现在我有一个包含 9 个物种的数千条记录的数据集，我想运行一个系统发育模型。我阅读了最常见软件包（例如 caper）的教程，但我不确定如何构建模型。

当我尝试为 caper 创建对象时，即使用：

obj <- comparative.data(phy = Tree, data = Data, names.col = species, vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE)

我收到消息：

row.names<-.data.frame(*tmp*, value = value) 中的错误：不允许重复的“row.names” 另外：警告信息：设置'row.names'时的非唯一值：'Species1'、'Species2'、'Species3'、'Species4'、'Species5'、'Species6'、'Species7'、'Species8'、'Species9'

我知道我可以通过应用 MCMCglmm 模型来解决这个问题，但我不熟悉贝叶斯模型。

提前感谢您的帮助。

【问题讨论】：

标签： mixed-models phylogeny

【解决方案1】：

这确实不适用于来自caper 的简单 PGLS，因为它无法将个人视为随机效应。我建议你使用MCMCglmm，它理解起来并不复杂，并且允许你将个体作为随机效应。您可以从包的作者here 或here 或更多处理包的某些特定方面（即树不确定性）here 的替代文档中找到优秀的文档。

真的很简短地让你开始：

## Your comparative data
comp_data <- comparative.data(phy = my_tree, data =my_data,
      names.col = species, vcv = TRUE)

请注意，您可以有一个如下所示的样本列：

   taxa        var1 var2 specimen
1     A  0.08730689    a    spec1
2     B  0.47092692    a    spec1
3     C -0.26302706    b    spec1
4     D  0.95807782    b    spec1
5     E  2.71590217    b    spec1
6     A -0.40752058    a    spec2
7     B -1.37192856    a    spec2
8     C  0.30634567    b    spec2
9     D -0.49828379    b    spec2
10    E  1.42722363    b    spec2

然后您可以设置您的公式（类似于简单的lm 公式）：

## Your formula
my_formula <- variable1 ~ variable2

还有你的 MCMC 设置：

## Setting the prior list (see the MCMCglmm course notes for details)
prior <- list(R = list(V=1, nu=0.002),
              G = list(G1 = list(V=1, nu=0.002)))

## Setting the MCMC parameters
## Number of interactions
nitt <- 12000

## Length of burnin
burnin <- 2000

## Amount of thinning
thin <- 5

然后您应该能够运行默认的MCMCglmm：

## Extracting the comparative data
mcmc_data <- comp_data$data

## As MCMCglmm requires a column named animal for it to identify it as a phylo
## model we include an extra column with the species names in it.
mcmc_data <- cbind(animal = rownames(mcmc_data), mcmc_data)
mcmc_tree <- comp_data$phy

## The MCMCglmmm
mod_mcmc <- MCMCglmm(fixed = my_formula, 
                     random = ~ animal + specimen, 
                     family = "gaussian",
                     pedigree = mcmc_tree, 
                     data = mcmc_data,
                     nitt = nitt,
                     burnin = burnin,
                     thin = thin,
                     prior = prior)

【讨论】：

非常感谢。会试一试，让你知道我是否可以处理。
将个体指定为随机因素有什么好处？我的理解是，它所做的只是考虑过度分散（正如我们在广义混合效应模型中所做的那样，正如 Zuur 等人在他们的一本书中所建议的那样）。我对么？提前谢谢！
是否有一个函数（可能在MCMCglmm）允许在每个物种存在多个观察值时“组装”比较数据？在@thomas-guillerme 建议的代码中，comp_data 仍然是使用来自caper 的comparative.data 创建的，这给出了@marco-gamba 描述的错误。干杯
您可以使用mulTree 包中的函数as.mulTree，该函数创建处理多个样本的对象，并且可以通过mulTree 传递给MCMCglmm。有关详细信息，请参阅 mulTree 文档示例。