In aliceyiwang/mvabund: Statistical Methods for Analysing Multivariate Abundance Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(knitr)

This vignette takes you through a case study where offsets are used during analysis with mvabund. We recommend reading the "Analysing discrete data" chapter from the Ecostats textbook.

Under some circumstances, we may tempted to standardise or average our abundance counts across replicates before analysing them. Often this is to account for differences in sampling intensity (different sized plots, different numbers of replicates), but we caution users from doing so! Standardisation/averaging counts can change the distribution of your data, and affect the mean-variance relationship. Instead you should use an offset, a term specially designed to account for known sources of variation among replicates.

Lets take a look at a case study that uses an offset in mvabund

Anthony wants to evaluate how well earthworms numbers were doing after some bush regeneration efforts. Here are some worm counts from pitfall traps across each of 10 sites (where C=control, R=bush regen). Notice there were only 4 pitfall traps at one site!

worms <- data.frame(matrix(nrow = 3, ncol = 10) )
rownames(worms) <- c("Treatment", "Count", "# pitfalls")
names(worms) <- NULL

worms[1,] <- c("C", "R", "R", "R", "C", "R", "R", "R", "R", "R")
worms[2,] <- c(0,3,1,3,1,2,12,1,18,0)
worms[3,] <- c(5,5,5,5,5,5,5,4,5,5)

knitr::kable(worms, "simple")

Anthony would like to model the number of worms (Haplotaxida) from the reveg dataset as a function of treatment (bush regeneration) while accounting variation in sampling effort.

Lets load mvabund and the dataset first:

library(mvabund)
library(ecostats)

data(reveg) 
attach(reveg)

skimr::skim(abund$Haplotaxida) # Great function to get an overview of the data

Now we will use the manyglm function. We specified family = "negative.binomial" as count data often follows a negative binomial distribution. We have also included an offset term to deal with the fact that sampling effort varied across sites – with 5 pitfall traps at most sites but just 4 at one site. Offsets are most commonly used in a log-linear model, where an offset is a predictor known to have an exactly proportional effect on the response. The offset in Anthony's model for is log(pitfalls). By including this term, our model changes from predicting mean abundance to modelling mean abundance per pitfall trap. This is what you would be trying to do if you averaged the counts across pitfalls prior to analysing them, but using an offset achieves this without messing with the data and its mean-variance relationship.

worms_offset <-  manyglm(Haplotaxida~treatment+offset(log(pitfalls)), family="negative.binomial", data=abund)

anova(worms_offset)

aliceyiwang/mvabund documentation built on March 13, 2024, 1:58 a.m.