Simple demos of PGLMM with traits and environment. YMMV...

playing-with-pglmm.r

# Headers
library(pez)
data(laja)

#############################
# SETUP #####################
#############################
    
# Make a comparative community object from a phylogeny and some
# site-species community data - some invertebrates found in some
# rivers
# - for the sake of "simplicity", let's work with presence/absence
#   data only
pres.abs <- river.sites
for(i in seq_len(ncol(pres.abs)))
    pres.abs[pres.abs[,i] != 0,i] <- 1
c.data <- comparative.comm(invert.tree, pres.abs)
c.data

# Convert from comparative.comm to long-format
l.data <- as.data.frame(c.data)

# Create a (scaled for ease of fitting and comparison)
# Variance-CoVariance matrix that represents whether species are
# closely-related
.standardise <- function(x) x <- x/(det(x)^(1/nrow(x)))
vcv <- .standardise(vcv(tree(c.data))) # Standardise species' phylogenetic structure

# IMPORTANT PIECES OF ADVICE:

# * It's important that species and sites are factors. as.data.frame
#   can get confused when sites are numbers (as here)
# * Non-conformable arguments means things aren't the same length, or
#   not of the right type (see above). I appreciate this isn't a
#   helpful error message, but it's also not wrong... :p

#############################
# Are closely-related #######
# species found together? ###
#############################
# Setup random effect terms
rnd.spp <- list(1, sp = l.data$species, covar = diag(length(species(c.data)))) # no phylogeny
rnd.attract <- list(1, sp = l.data$species, covar = vcv) # close relatives similar

# Fit model and examine
simple.model <- communityPGLMM(presence ~ 1, data=l.data, family="binomial",
                               sp=l.data$species, site=factor(l.data$site),
                               random.effects=list(rnd.spp,rnd.attract), REML=FALSE)

summary(simple.model)
# ...the intercept is significant, great...

# Check "significance" of the random effect terms
communityPGLMM.binary.LRT(simple.model, 1)
communityPGLMM.binary.LRT(simple.model, 2)
# ... and they're not, which is a shame, but hey ho

#############################
# Controlling for site- #####
# specific-variation ########
#############################
rnd.site <- list(1, site = factor(l.data$site), covar = diag(length(sites(c.data))))
site.model <- communityPGLMM(presence ~ 1, data=l.data, family="binomial",
                               sp=l.data$species, site=factor(l.data$site),
                               random.effects=list(rnd.spp,rnd.attract,rnd.site), REML=FALSE)
summary(site.model)
communityPGLMM.binary.LRT(site.model, 1)
communityPGLMM.binary.LRT(site.model, 2)
communityPGLMM.binary.LRT(site.model, 3)
#...effects of site and species, but nothing for phylogeny :-(

#############################
# Incorporating species' ####
# traits ####################
#############################
# Add in the trait data
c.data <- comparative.comm(invert.tree, pres.abs, invert.traits)
l.data <- as.data.frame(c.data)

# Fit model
trait.model <- communityPGLMM(presence ~ length + fish.pref, data=l.data, family="binomial",
                               sp=l.data$species, site=factor(l.data$site),
                               random.effects=list(rnd.spp,rnd.attract,rnd.site), REML=FALSE)

summary(trait.model)
#Assessing the significance of fixed effects is a bit easier

#############################
# Incorporating environemnt #
#############################
# Add in the environmental data
c.data <- comparative.comm(invert.tree, pres.abs, invert.traits, river.env)
l.data <- as.data.frame(c.data)

# Fit model
env.model <- communityPGLMM(presence ~ fish.pref + temperature, data=l.data, family="binomial",
                               sp=l.data$species, site=factor(l.data$site),
                               random.effects=list(rnd.spp,rnd.attract,rnd.site), REML=FALSE)

summary(env.model)
#... you could do other variables and interactions as fixed effects if
#you wished. It would be nice if something were significant!... :-(

fixed.interactions <- communityPGLMM(presence ~ fish.pref * temperature, data=l.data, family="binomial",
                               sp=l.data$species, site=factor(l.data$site),
                               random.effects=list(rnd.spp,rnd.attract,rnd.site), REML=FALSE)
#...which might be interesting to contrast with the below

#############################
# Interactions  #############
#############################
# For traits
rnd.attract.trt <- list(l.data$length, sp = l.data$species, covar = vcv) # close relatives similar
trait.model <- communityPGLMM(presence ~ length + fish.pref, data=l.data, family="binomial",
                              sp=l.data$species, site=factor(l.data$site),
                              random.effects=list(rnd.spp,rnd.attract,rnd.site,rnd.attract.trt),
                              REML=FALSE)


# For traits
rnd.spp.trt <- list(l.data$length, sp = l.data$species, covar = diag(length(species(c.data)))) # close relatives similar
rnd.attract.trt <- list(l.data$length, sp = l.data$species, covar = vcv) # close relatives similar
trait.phy.model <- communityPGLMM(presence ~ length + fish.pref, data=l.data, family="binomial",
                              sp=l.data$species, site=factor(l.data$site),
                              random.effects=list(rnd.spp,rnd.attract,rnd.site,rnd.attract.trt),
                              REML=FALSE)
#...it's basically just the same as the pure-species model, but with
#   the trait "length" in the first argument slot. That's all there is
#   to it - if you want to have interactions any other way, just plonk
#   them in that first slot. Does that make sense?