Как создать многомерную линейную модель с коррекцией филогенеза c?

Question

Я пытаюсь использовать пакет MCMCglmm для запуска многомерной линейной модели, которая филогенетически исправлена. Вот подмножество моих данных (в общей сложности 210 записей для 67 видов и 6 кластеров):

Names PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 tips 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"

Когда я пытаюсь запустить одномерную модель (код ниже, phylo относится к моему дереву ) это работает:

Ainv<-inverseA(phylo,nodes="TIPS",scale=F)$Ainv
p.var=var(data[,c("PC1")]) 
prior1<-list(R=list(V=(p.var),nu=0.002),G=list(G1=list(V=(p.var),nu=0.002)))
m7.phylo<-MCMCglmm(PC1~Clusters,
                   random=~tips,
                   family=rep("gaussian",1),
                   ginverse=list(tips=Ainv),
                   data=data,
                   prior=prior1)

Однако, как только я пытаюсь сделать многомерную модель, я получаю ошибку:

Ainv<-inverseA(phylo,nodes="TIPS",scale=F)$Ainv
p.var=var(data[,c("PC1","PC2")]) 
prior1<-list(R=list(V=(diag(2)*p.var),nu=0.002),G=list(G1=list(V=(diag(2)*p.var),nu=0.002)))
m7.phylo<-MCMCglmm(cbind(PC1,PC2)~Clusters,
                   random=~tips,
                   family=rep("gaussian",2),
                   ginverse=list(tips=Ainv),
                   data=data,
                   prior=prior1)

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

Есть ли у кого-нибудь какие-либо идеи о том, что могло бы пойти не так?

Спасибо,

Каролина

Answer 1

Итак, этот код работает:

Ainv<-inverseA(phylo,nodes="TIPS",scale=F)$Ainv
p.var=var(data[,c(1:11)]) 
prior2<-list(R=list(V=p.var,nu=0.002),G=list(G1=list(V=p.var,nu=0.002)))
model2<-MCMCglmm(cbind(PC1,PC2,PC3,PC4,PC5,PC6,PC7,PC8,PC9,PC10,PC11)~Clusters,
                   random = ~us(trait):tips,
                   rcov = ~idh(trait):units,
                   family=rep("gaussian",11),
                   ginverse=list(tips=Ainv),
                   data=data,
                   prior=prior2)

И он, кажется, делает то, что я хочу, но я действительно не уверен, что знаю достаточно, чтобы сказать это наверняка! Любой вклад все равно будет оценен!

Редактировать: поговорив об этом с другими людьми, похоже, это соответствует решению, которое они дали, поэтому я буду считать его правильным.