[R-sig-ME] Phylogenetic component using lme4ord

[Gustaf Granath]

Hi all,
Im testing different ways to quantify/incorporate the phylogenetic 
component in the analyzes of trait data (e.g. phylo signal/lambda). I 
noticed that lme4ord has progressed and started to play around with it. 
Im particularly interested in using data with with-species samples 
(PGLMM). Im not sure I understand the fitting/output from lme4ord though.

# First load packages and data
require(phytools)
require(geiger)
require(MCMCglmm)
data(bird.families)

# this data has high phylogentic signal so we lower it a bit
library(geiger)
set.seed(350)
bird.families.rs <- rescale(bird.families, 'lambda', lambda=0.5)
phylo.effect <- fastBM(bird.families.rs)
phylosig(bird.families.rs, phylo.effect, method='lambda', test = TRUE)
# lambda ~0.7

# Lets add within-species variation
phylo.effect.rep <- data.frame(y = c(phylo.effect,
                                     rnorm(length(phylo.effect), 0, 1) + 
phylo.effect,
                                     rnorm(length(phylo.effect), 0, 1) + 
phylo.effect),
                                 species = rep(names(phylo.effect), 3))

# Anlyze data using lme4ord (install from github)
# I think phyloEdge() was earlier pagelLamda()
library(lme4ord)
formula <- y ~ 1 + (1| species) + phyloEdge(1 | species, phylo = 
bird.families.rs)
mod <- strucGlmer(formula, family = gaussian(), data = phylo.effect.rep, 
optMaxit = 1000000)
summary(mod)
# sumary() does not give random effects
getReTrm(mod)
# species.phyloEdge component VERY small.
# It seems like only 2 components are estimated.
# I thought there would be 3: phylo, species (between species variation 
not captured by the BM model), within-species (residual).

# And now compare with MCMCglmm
# WARNING slow!!
library(MCMCglmm)
Ainv <- inverseA(bird.families.rs, scale = TRUE)$Ainv
prior = list(R = list(V = 1, nu = 0.02),
              G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 
1000),
                       G2 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 
1000)))
phylo.effect.rep$species.ide <- phylo.effect.rep$species

m1.phy <- MCMCglmm(y ~ 1, random =~ species + species.ide,
                     data = phylo.effect.rep, ginverse = list(species = 
Ainv),
                     nitt = 700000, thin=100, prior = prior, burnin = 
50000, pr=TRUE)

# without species.ide
prior = list(R = list(V = 1, nu = 0.02),
              G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 
1000)))
m2.phy <- MCMCglmm(y ~ 1, random =~ species,
                     data = phylo.effect.rep, ginverse = list(species = 
Ainv),
                     nitt = 500000, thin=100, prior = prior, burnin = 
50000, pr=TRUE)

summary(m1.phy)
# 3 components estimated
summary(m2.phy)
# in m2.phy species and species.ide are merged into one component.
# Thus Lambda cannot be estimated for m2.phy.

# Lambda
# Here we assume within-species variance is measurement error so not 
included in the denominator
lambda.1 <- 
m1.phy$VCV[,"species"]/(m1.phy$VCV[,"species"]+m1.phy$VCV[,"species.ide"])
mean(lambda.1)
# ~0.78. Not too bad and longer chains (more samples) may improve this
# however the posterior of species.ide does not look good! But it seems 
to give an OK result.

Question: can the MCMCglmm results be replicated using lme4ord? What 
model is actually fitted by lme4ord?

Hope someone can help.

Cheers,

Gustaf


-- 
Gustaf Granath (PhD)
Post doc
Swedish University of Agricultural Sciences
Department of Ecology