[R-sig-ME] Creating a phylogenetically corrected multivariate linear model using MCMCglmm

Thu Apr 30 10:30:55 CEST 2020

Hi all,

I'm trying to use the package MCMCglmm to run a multivariate linear model that is phylogenetically corrected. Here is a subset of my data (there are 210 entries in total for 67 species and 6 clusters):

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"

I created a univariate model taking inspiration from this link:

https://github.com/JonBrommer/Multivariate-Mixed-Models-in-R/wiki/MCMCglmm-examples#organisational-level-4https://github.com/JonBrommer/Multivariate-Mixed-Models-in-R/wiki/MCMCglmm-examples#organisational-level-4

And when I try to run it, it works (code below, phylo refers to my tree):

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)

However, as soon as I try to make a multivariate model, I get an error:

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

I do want to use all 11 PCs in the model so it's not encouraging that I can't seem to get it to work with just 2 of them...  Does anyone have any ideas of what could have gone wrong?  I would like to create other models with the same structure using different data and trees so it would be good to understand what's going on and how to create a prior properly for next time.

Thanks,

Carolina

	[[alternative HTML version deleted]]