Creating a phylogenetically corrected multivariate linear model using MCMCglmm - R-SIG-mixed-models

Fri, May 1, 2020 3:28 AM #

Hi Carolina,

Looks like you have used the prior from Jon's 'multivariate
phylogenetic example' section but haven't updated the random and rcov
terms in your MCMCglmm model accordingly, which should be something
like:

random=~us(trait):tips, rcov=~us(trait):units

Hope that helps...

Cheers

Tom

Date: Thu, 30 Apr 2020 08:30:55 +0000
From: Carolina Karoullas
        <carolina.karoullas at postgrad.manchester.ac.uk>
To: "r-sig-mixed-models at r-project.org"
        <r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Creating a phylogenetically corrected multivariate
        linear model using MCMCglmm
Message-ID:
        <56C3CB6634CED94898F1E5FC8E71A09CB602BF47 at MBXP03.ds.man.ac.uk>
Content-Type: text/plain; charset="utf-8"

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