使用MCMCglmm对每个物种的多个条目进行系统发育校正

问题描述 投票:0回答:1

我有一个210行(67种)的数据集(data),11个连续变量(从PC1到PC11)和1个分类变量(具有6个级别的簇)和一个系统树(phylo)。这是我的数据集的一个子集,可以更清楚地显示格式:

Names   PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10    PC11    species Clusters
Accipiter_gentilis1 -1.34E-01   1.98E-02    -1.72E-02   4.00E-02    -1.93E-03   2.45E-02    -1.28E-04   1.03E-02    9.09E-03    -5.33E-03   -1.30E-02   Accipiter_gentilis  "Soaring"
Accipiter_gentilis2 -1.26E-01   1.22E-02    1.66E-02    9.22E-03    1.15E-02    1.68E-02    -1.50E-02   1.02E-02    7.93E-03    -8.47E-03   -1.02E-02   Accipiter_gentilis  "Soaring"
Accipiter_nisus1    -1.81E-01   5.76E-02    -1.82E-02   6.15E-02    -9.25E-03   3.40E-02    -1.77E-02   5.45E-03    7.01E-03    -2.07E-02   -8.78E-03   Accipiter_nisus "Soaring"
Accipiter_nisus2    -2.00E-01   7.05E-02    -1.12E-02   5.94E-02    3.49E-03    3.10E-02    -1.58E-03   -1.55E-03   6.92E-03    -3.54E-02   -1.80E-02   Accipiter_nisus "Soaring"
Accipiter_nisus3    -8.14E-02   -3.39E-04   -8.88E-03   4.25E-02    -5.48E-04   -8.51E-03   5.07E-03    4.56E-03    1.97E-02    -1.46E-02   -1.43E-02   Accipiter_nisus "Soaring"
Accipiter_nisus4    -2.06E-01   7.05E-02    -2.17E-02   6.38E-02    -1.61E-02   2.80E-02    8.70E-03    -5.96E-03   6.15E-03    -5.29E-02   -2.05E-02   Accipiter_nisus "Soaring"
Actitis_hypoleucos1 2.27E-02    -2.74E-03   4.79E-02    -2.30E-02   -2.76E-02   -2.36E-02   1.70E-02    2.43E-03    3.82E-03    1.15E-02    -9.87E-03   Actitis_hypoleucos  "Continuous flapping"
Actitis_hypoleucos2 6.67E-02    -1.05E-02   5.12E-02    -2.65E-02   -3.21E-02   -2.61E-02   3.21E-03    7.46E-03    7.29E-03    4.70E-03    -1.37E-02   Actitis_hypoleucos  "Continuous flapping"
Aix_sponsa1 -3.70E-02   -1.41E-02   1.13E-02    3.16E-02    2.32E-02    -1.70E-02   2.32E-02    1.91E-03    2.91E-02    -7.71E-03   7.40E-03    Aix_sponsa  "Continuous flapping"
Aix_sponsa2 1.03E-02    -4.08E-02   -6.62E-03   1.19E-02    2.83E-02    -1.49E-02   3.78E-02    6.98E-03    2.91E-02    -4.32E-03   2.54E-03    Aix_sponsa  "Continuous flapping"
Aix_sponsa3 1.19E-02    -3.48E-02   -1.53E-02   3.76E-03    2.17E-02    1.47E-02    8.84E-03    2.39E-02    -9.20E-03   -1.78E-02   8.76E-04    Aix_sponsa  "Continuous flapping"
Aix_sponsa4 -3.37E-02   -1.75E-02   -8.06E-03   3.64E-02    -5.50E-03   1.03E-02    2.37E-02    3.33E-03    -1.04E-03   -2.00E-02   5.89E-03    Aix_sponsa  "Continuous flapping"
Aix_sponsa5 -2.30E-02   -9.59E-03   1.06E-02    3.01E-02    7.10E-03    -1.23E-02   2.08E-02    1.17E-02    1.59E-03    2.83E-03    8.75E-03    Aix_sponsa  "Continuous flapping"
Aix_sponsa6 -1.70E-02   -2.98E-02   -1.96E-02   1.76E-02    1.23E-02    4.92E-03    5.45E-03    1.99E-02    -6.43E-03   -9.63E-04   1.99E-03    Aix_sponsa  "Continuous flapping"
Aix_sponsa7 2.57E-02    -4.22E-02   -1.60E-02   1.75E-02    3.41E-03    5.80E-03    2.89E-02    6.10E-03    7.12E-03    2.75E-03    4.99E-03    Aix_sponsa  "Continuous flapping"
Aix_sponsa8 4.09E-02    -4.46E-02   -4.49E-03   2.24E-02    2.37E-03    -5.90E-03   2.78E-02    -8.26E-04   1.17E-02    -5.71E-03   -1.77E-03   Aix_sponsa  "Continuous flapping"

我正在尝试查看我的连续变量是否可以在校正系统发育时用于预测我的类别。理想情况下,我将使用多元系统发育校正的判别分析,但似乎不存在这样的事情(是吗?)。因此,我正在尝试做下一个最好的事情,并使用MCMCglmm包实施多元系统发育分析。我尝试遵循此网站(https://github.com/JonBrommer/Multivariate-Mixed-Models-in-R/wiki/MCMCglmm-examples#organisational-level-4)上所做的操作,该网站还提供了要使用的数据集和树的链接,并获得了以下信息:

Ainv.phylo<-inverseA(phylo,nodes="TIPS",scale=F)$Ainv
phen.var<-apply(data[,1:11],2,function(m) var(m,na.rm=T))
prior.phylo<-list(R=list(V=diag(11)*phen.var, nu=0.002),
                      G=list(G1=list(V=diag(11)*phen.var, nu=0.002)))
model<-MCMCglmm(cbind(PC1,PC2,PC3,PC4,PC5,PC6,PC7,PC8,PC9,PC10,PC11)~Clusters,
                       random=~us(Clusters):species,
                       rcov=~us(Clusters):units,
                       ginverse=list(tips=Ainv.phylo),
                       family=c("gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian"),
                       data=data,
                       prior=prior.phylo)
summary(model)

但是,我不断收到此错误:

Error in priorformat(if (NOpriorG) { : 
  V is the wrong dimension for some prior$G/prior$R elements

我真的很努力地思考如何创建先验和应该使用的维数(我很抱歉-使用矩阵从来都不是我的强项),因此任何帮助或建议都将不胜感激。

谢谢,

卡罗琳娜

r matrix matrix-multiplication mcmc phylogeny
1个回答
0
投票

我真的不是专家,但是之前已经运行了一些MCMCglmm模型,并且您的模型结构和先前的结构看起来还不错,所以我认为可能是Ainv.phylo或phen.var的尺寸正在改变尺寸先前的V值(给您错误消息)。

由于您希望先前的V矩阵为11x11,所以我将尝试使用简单的V = diag(11)结构运行模型,看看是否可行,然后慢慢添加phen.var,一点一点地算出哪个命令是弄乱了您的11x11矩阵?

© www.soinside.com 2019 - 2024. All rights reserved.