Analyse transect data from Shultz et al (2016) PeerJ

data-plus-effects-original-scale.png

data-plus-effects-sqrt-scale.png

model-effects.png

transect-model.R

## Analyse full transect data of J.A. Shultz et al (2016) PeerJ 4:e1980; DOI 10.7717/peerj.1980
## https://doi.org/10.7717/peerj.1980

## packages
library('readr')
library('dplyr')
library('ggplot2')
library('tidyr')
library('lme4')
library("ggstance")

## ggplot2 theme
theme_set(theme_bw())

## read in the data
transects <- read_csv("transectcounts.csv")

## transects in long form
lt <- gather(transects, species, count, - ssws, - site, - transect)
lt

## convert to factors for modelling
lt <- mutate(lt,
             ssws = factor(ssws, levels = c('before', 'after')),
             site = factor(site),
             transect = factor(transect))

## subset down to the two species of interest
spp <- filter(lt, species %in% c('sunflower.star', 'green.urchin'))

## update the species names
spp <- mutate(spp, species = case_when(species == 'sunflower.star' ~ 'Sunflower Star',
                                       species == 'green.urchin'   ~ 'Green Urchin'),
              species = factor(species, levels = c('Sunflower Star', 'Green Urchin')))

## Model fitting - sunflower
m.sun <- glmer(count ~ ssws + (1 | site / transect), data = spp, family = poisson,
               subset = species == 'Sunflower Star')
m0.sun <- update(m.sun, . ~ . - ssws)

## clear support for an effect
anova(m.sun, m0.sun)

## model fitting - green urchin
m.gur <- glmer(count ~ ssws + (1 | site / transect), data = spp, family = poisson,
               subset = species == 'Green Urchin')
m0.gur <- update(m.gur, . ~ . - ssws)

## clear support for an effect
anova(m.gur, m0.gur)

## Let's assume these are good models -
##   should really do a lot of diagnostics...

## Extract data from the model for plotting
ilink <- family(m.sun)$linkinv
effs.s <- tibble(species = 'Sunflower Star', ssws = c('before', 'after'), fitted = fixef(m.sun),
                 se = diag(vcov(m.sun)))
effs.g <- tibble(species = 'Green Urchin', ssws = c('before', 'after'), fitted = fixef(m.gur),
                 se = diag(vcov(m.gur)))

effs <- bind_rows(effs.s, effs.g)
effs <- mutate(effs, est = ilink(fitted), upper = ilink(fitted + (2*se)),
               lower = ilink(fitted - (2*se)),
               species = factor(species, levels = c('Sunflower Star','Green Urchin')),
               ssws = factor(ssws, levels = c('before', 'after')))

## base plot
base.plt <- ggplot() +
    facet_wrap(~ species, scales = 'free_y') +
    labs(x = 'Sea star mortality event')

## plot of effects
eff.plt <- base.plt +
    geom_pointrange(aes(y = est, x = ssws, ymax = upper, ymin = lower), data = effs) +
    labs(y = 'Expected count')
eff.plt

ggsave('model-effects.png', eff.plt, width = 7, height = 3.5)

## plot of effects on data
data.plt <- base.plt +
    geom_point(aes(y = count, x = ssws), data = spp, alpha = 0.2,
               position = position_jitter(width = 0.2)) +
    geom_pointrange(aes(y = est, x = ssws, ymax = upper, ymin = lower, colour = ssws),
                    data = effs) +
    theme(legend.position = 'none') +
    labs(y = 'Number of individuals')
data.plt

ggsave('data-plus-effects-original-scale.png', data.plt, width = 7, height = 3.5)

## Data on sqrt scale
data.plt + scale_y_sqrt(minor_breaks = NULL) +
    expand_limits(y = 0)

ggsave('data-plus-effects-sqrt-scale.png', width = 7, height = 3.5)

transectcounts.csv