Skip to content

Instantly share code, notes, and snippets.

@willpearse
Created October 30, 2017 14:12
Show Gist options
  • Save willpearse/54614ea7269681c021773ccf3e9d8e52 to your computer and use it in GitHub Desktop.
Save willpearse/54614ea7269681c021773ccf3e9d8e52 to your computer and use it in GitHub Desktop.
Simple demos of PGLMM with traits and environment. YMMV...
# 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?
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment