In brouwern/DRmencia: Data and code to replicate analyses for Latta, Brouwer, et al. 2018. PeerJ.

knitr::opts_chunk$set(echo = TRUE)

Introduction

GLMM and MANOVA analyses of leaf litter insects. See seperate script for flying insects from sticky traps. Annotation in that script is generally more thorough than this one.

Load Leaf flyingl data

load("./data/flying.RData")
library(reshape2)
library(vegan)
library(lme4)
library(blme)
library(ggplot2)
library(cowplot)
library(merTools)
source('~/1_R/git/git-aviary/DRmencia/R/back_transform_reals.R')

Reshape wide data to community matrix

flying.comm.matrix <- dcast(data = flying,
                            formula = site + sample ~ taxa,
                            value.var = "N",
                            fun.aggregate = sum)

#remove any rows that are all zero; 
flying.comm.matrix <- na.omit(flying.comm.matrix)

Reshape wide to total abundance

Total number of arthropods in a sample, ignoring taxa

flying.tot.N <- dcast(data = flying,
                            formula = site + sample ~ .,
                            value.var = "N",
                            fun.aggregate = sum)
names(flying.tot.N)[3] <- "tot.N"


#reset factor levels
flying.tot.N$site2 <- factor(flying.tot.N$site,
                        levels =c("La Cueva",
                                  "La Caoba",
                                  "Morelia",
                                  "El Corral",
                                "Aceitillar"))

Set ages for regression

Not age for the broadleaf forest at Aceitillar

flying.tot.N$site.age.init <- NA
flying.tot.N$site.age.init[which(flying.tot.N$site == "La Cueva")] <- 2
flying.tot.N$site.age.init[which(flying.tot.N$site == "La Caoba")] <- 5
flying.tot.N$site.age.init[which(flying.tot.N$site == "Morelia")] <- 10 
flying.tot.N$site.age.init[which(flying.tot.N$site == "El Corral")] <- 20

Plot

Plot total abundance

Plot means

#mean SE/CIs on the fly
ggplot(data = flying.tot.N,
       aes(y = tot.N,
           x = site)) +
  #ylim(25,85) +
  stat_summary(fun.data  = "mean_cl_boot",size = 2)

Plot means w/ age as continous

pd <- position_dodge(0.5)
i.max <- which(flying.tot.N$tot.N > 100 |
               flying.tot.N$site == "Broadleaf" |
               is.na(flying.tot.N$tot.N) == TRUE)
ggplot(data = flying.tot.N[-i.max,],
       aes(y = tot.N,
           x = site.age.init)) +
  ylim(25,75) +
  geom_smooth(method=lm)+
  #geom_point(position = pd) +
  stat_summary(fun.data  = "mean_cl_boot",size = 2)

Calculate percentages

flying2 <- merge(flying,flying.tot.N)

dim(flying2)

dim(flying)

Plot percentages

flying2$site2 <- flying2$site
flying2$site2 <- gsub("La ","",flying2$site2)
flying2$site2 <- gsub("El ","",flying2$site2)
flying2$site2 <- gsub("Broadleaf","Aceitillar",flying2$site2)

summary(factor(flying2$site))
summary(factor(flying2$site2))
flying2$site2 <- factor(flying2$site2,
                        levels =c("Cueva",
                                  "Caoba",
                                  "Morelia",
                                  "Corral",
                                "Aceitillar"))
summary(flying$site)

ggplot(data = flying2,
       aes(y = N/tot.N,
           x = site2)) +
  stat_summary(fun.data  = "mean_cl_boot") +
  facet_wrap(~taxa,
             scale = "free") +
  xlab("Site") +
  ylab("Percentage composition") +
  ggtitle("Flying arthropods")

Model

flying.tot.N$i <- 1:dim(flying.tot.N)[1]
flying.tot.N$year <- flying.tot.N$site.age.init
i.aceit <- which(flying.tot.N$site == "Broadleaf")
flying.regress0 <- bglmer(tot.N ~ 1 +
                (1|site) +
                (1|i),
              family = poisson,
              data = flying.tot.N[-i.aceit,])

flying.regress1 <- update(flying.regress0, 
                          . ~ . + site.age.init)

anova(flying.regress0,
      flying.regress1)

i.max2 <- which(flying.tot.N$tot.N > 100 |
               #flying.tot.N$site == "Broadleaf" |
               is.na(flying.tot.N$tot.N) == TRUE)

#don't use means parameterization; doesn't play nice
#with downstream plotting
flying.anova <- bglmer(tot.N ~  site +
                #(1|site) +
                (1|i),
              family = poisson,
              data = flying.tot.N[,])

BUild models which pool all pasture sites

#make factor that designates Aceit vs. pasture
flying.tot.N$Aceit.YN <- ifelse(flying.tot.N$site == "Broadleaf","Aceit.YES",
                                "Aceit.NO")

#model
flying.pasture.vs.broadleaf.0 <- bglmer(tot.N ~ 1 +
                                         (1|site) +
                (1|i),
              family = poisson,
              data = flying.tot.N[-i.max2,])
flying.pasture.vs.broadleaf.1 <- update(flying.pasture.vs.broadleaf.0,
                                        .~. + Aceit.YN)

#test 
anova(flying.pasture.vs.broadleaf.0,
      flying.pasture.vs.broadleaf.1)

Calcualte CIs: determine by what factor they are different

#ci.out <-  confint(flying.pasture.vs.broadleaf.1,method = "boot")

                   2.5 %    97.5 %

.sig01 0.19080235 0.3162773 .sig02 0.06657935 0.3454228 (Intercept) 3.43297028 3.8569679 Aceit.YNAceit.YES 0.02347341 0.9440533

Plot

# data for prediction
newdat.flying.reg <- expand.grid(
     site.age.init = seq(min(flying.regress1@frame$site.age.init),
                       max(flying.regress1@frame$site.age.init))
    ,site = unique(flying.regress1@frame$site)[1]
    ,i = flying.regress1@frame$i[1]
    )

newdat.flying.anova <- expand.grid(
      site = unique(flying.anova@frame$site)
      ,i = flying.anova@frame$i[1]
    )

## predictions
flying.rgrxn.pred <- merTools::predictInterval(flying.regress1,
                                   newdata =   newdat.flying.reg,
                                   which = "fixed",
                                   level = 0.95,
                                   #n.sims = 10,
                                   stat = "median",
                                   include.resid.var = FALSE)

flying.rgrxn.pred2 <- back_transform_reals(flying.rgrxn.pred,fnxn = "exp")

flying.rgrxn.pred3 <- cbind(newdat.flying.reg
                            ,flying.rgrxn.pred2)

flying.anova.pred <- merTools::predictInterval(flying.anova,
                                   newdata =   newdat.flying.anova,
                                   which = "fixed",
                                   level = 0.95,
                                   #n.sims = 10,
                                   stat = "median",
                                   include.resid.var = FALSE)

flying.anova.pred2 <- back_transform_reals(flying.anova.pred,fnxn = "exp")
flying.anova.pred3 <- cbind(newdat.flying.anova,flying.anova.pred2)

ymin. <- 15
ymax.<- 100

gg.out.fly.reg <-  ggplot(data = flying.rgrxn.pred3,
         aes(y = fit.real,
             x = site.age.init,
             color = site,
             shape = site)) +
    geom_line(aes()) +
   ylim(ymin.,ymax.) +
    xlab("Site age") +
    ylab("Flying arthropod abundance") +
   scale_shape_manual(values=c(4,0, 1, 2, 5)) +

    geom_ribbon(aes(ymax = upr.real,
                    ymin = lwr.real),
                alpha = 0.05,
                linetype = 0)  +
    #ylim(ymin., ymax.) +
    theme(legend.position="none")

  ## add raw data
  gg.out.fly.reg1 <- gg.out.fly.reg + geom_point(data = flying.regress1@frame,
                      aes(y = flying.regress1@frame[,1],
                          x = site.age.init),
                      size = 3)

  ## add ANOVA dat
  nudge. <- 0.5

  flying.anova.pred3$site.age.init <- c(2,5,10,20, NA)
  gg.out.fly.reg2 <- gg.out.fly.reg1 + geom_point(data = flying.anova.pred3,
                                aes(y = fit.real,
                                    x =  site.age.init+nudge.),
                                size = 5,shape =16,
                                position = pd) +
                 geom_errorbar(data = flying.anova.pred3,
                               aes(x = site.age.init+nudge.,
                                  ymax = upr.real,
                                   ymin = lwr.real),
                               width = 0)

flying.anova.pred3$site <- gsub("Broadleaf",
                                "Aceitillar",
                                flying.anova.pred3$site)

flying.anova@frame$site <- gsub("Broadleaf",
                                "Aceitillar",
                                flying.anova@frame$site)


i.broad <- which(flying.anova.pred3$site == "Aceitillar")
i.broad2 <- which(flying.anova@frame$site == "Aceitillar")

y.mean <- flying.anova.pred3[i.broad,"fit.real"]
ymin.CI <- flying.anova.pred3[i.broad,"lwr.real"]
ymax.CI <- flying.anova.pred3[i.broad,"upr.real"]

gg.anova.fly <-  ggplot(data = flying.anova.pred3[i.broad,],
                    aes(y = fit.real,
                        x = site)) +
    geom_line(aes(),
              size = 1) +
    ylim(ymin.,ymax.) +
     xlab("") +
    # ylab(ylab.) +
  annotate("pointrange", x = 1.25, 
           y = y.mean,
           ymin = ymin.CI, 
           ymax = ymax.CI,
  size = 1.5,
  fatten = 2) +
    # geom_errorbar(aes(ymax = upr.real,
    #                 ymin = lwr.real),
    #             width = 0)  +
    # geom_point(size = 3) +
    theme(legend.position="none")# +
    #ylim(ymin., ymax.)

  ## add raw data
#  pd <- position_dodge(0.25)
gg.anova.fly2 <- gg.anova.fly + geom_jitter(data = flying.anova@frame[i.broad2,],
                                #position = pd,
                                width  = 0.125,
                                shape = 10,
                                aes(y = flying.anova@frame[i.broad2,1],
                                    x = site)) +
    theme(axis.line.y=element_blank(),
          #axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks=element_blank(),
          #axis.title.x=element_blank(),
          axis.title.y=element_blank())

Check plots

cowplot::plot_grid(gg.out.fly.reg2,
                   gg.anova.fly2,
                   rel_widths = c(8,1))

Save plots

save(gg.out.fly.reg2, file = "gg_flying_insects_regression.RData")
save(gg.anova.fly2, file = "gg_flying_insects_anova.RData")

Analyis

Manova

with(na.omit(flying.comm.matrix), table(site, sample))
i.drop <- which(flying.comm.matrix$sample > 10)

manova.dat <- flying.comm.matrix[-i.drop,]
manova.out <- manova(cbind(Coleoptera,Diptera,Homoptera,Hymenoptera,
             Lepidoptera,Orthoptera,Spiders,
             Unidentified) ~ site, data = manova.dat )
summary(manova.out)
summary.aov(manova.out)

      Df Pillai approx F num Df den Df   Pr(>F)

site 4 1.1174 1.9867 32 164 0.002933 ** Residuals 45

Save workspac

save.image(file = "WORKSPACE_leaf_litter_insect_analysis.RData")

brouwern/DRmencia documentation built on May 6, 2019, 12:24 p.m.