R-sig-mixed-models Digest, Vol 161, Issue 1 - R-SIG-mixed-models

Fri, May 1, 2020 3:20 AM #

Hi Carolina, your second model assumes homogenous tips variance across variables (random=~tips) - but in your prior you structure both tips and residual variances as 2x2 matrices. In the second model you have to use us/idh variance function for tips effects (e.g. us(trait):tips), and add rcov term to model similar structure for units (e.g. rcov=~us(trait):units).

Cheers,
Szymek



Dr hab. Szymon Drobniak

Institute of Environmental Sciences
Jagiellonian University, Krak?w, Poland

School of Biological, Environmental and Earth Sciences
University of New South Wales, Sydney, Australia

Google Scholar profile
szymekdrobniak.wordpress.com
szdrobniak.pl

Message: 1
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