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

knitr::opts_chunk$set(echo = TRUE)

Preliminaries

Libraries

Function that calls all libraries needed for analysis

source("./R/load_libraries.R")
load_libraries()

Load data

• Data has just focal species; rarities removed
• Within-year repatures removed

load(file = "./data/ann_counts") #loads "ann_counts"

Remove overlap of Youngest sites

The 2 youngest sites both have ages that equal 5 and 6 years: La Cueva's last 2 years, La Caoaba's 1st 2 years. Using "site" as a factor variable in ANOVA will confound site and age. This code subsets the data for these sites so thatall Cueva data is younger than Caoba data.

with(ann_counts, table(site,site.age))

summary(ann_counts[is.na(ann_counts$site.age),])

i.cueva.drop <- which(ann_counts$site == "La Cueva" & 
                      ann_counts$site.age > 4) #keep age 2, 3, 4
i.caoba.drop <- which(ann_counts$site == "La Caoba" & 
                      ann_counts$site.age < 7) #keep age 7, 8, 9

ann_counts_2 <- ann_counts[-c(i.cueva.drop,i.caoba.drop), ]

Remove "D" habitat

dim(ann_counts)
i.hab.D <- which(ann_counts$hab1 == "D" )


                 #ann_counts$spp.code %in% c("CGDO","GREP"))

length(i.hab.D)

ann_counts_3 <- ann_counts_2[-i.hab.D, ]

summary(ann_counts_3$site)
with(ann_counts_3, table(spp.code, site))

Model

Base ANOVA model - age ratio

Note that this model uses all species, if even they are rare at Aceitillar; for the Aceitillar only analysis really rare spp were removed.

m.Nvshab.ANOVA.0 <- bglmer(N ~ 1 +
               (1|i) +
               (1|year) +
               #(1|site) +
               (1|spp.code) ,#+
               #(1|site:spp.code) #this would be ideal but prob overkill
              data = ann_counts_3[,],
              family = poisson,
              glmerControl(optimizer = "Nelder_Mead")
              )

Run fixed effects

m.Nvshab.ANOVA.site          <- update(m.Nvshab.ANOVA.0, . ~ . + site)
m.Nvshab.ANOVA.status.focals <- update(m.Nvshab.ANOVA.0, . ~ . +         hab1)
m.Nvshab.ANOVA.add           <- update(m.Nvshab.ANOVA.0, . ~ . +  site + hab1)
m.Nvshab.ANOVA.full          <- update(m.Nvshab.ANOVA.0, . ~ . +       site*hab1)

#calcualte mean values for each year
m.Nvshab.ANOVA.means         <- update(m.Nvshab.ANOVA.0, . ~ . +  -1 + site:hab1)

fixef(m.Nvshab.ANOVA.means)

Inference

ICtab(m.Nvshab.ANOVA.0,
       m.Nvshab.ANOVA.site,
       m.Nvshab.ANOVA.status.focals,
       m.Nvshab.ANOVA.add,
       m.Nvshab.ANOVA.full,
        type = "AICc")

Contrast modeling

Create contrast matrix

contrst.mat <- rbind("Mig-linear" = c(-2,-1, 0, 1,  2, 0, 0, 0, 0, 0),
                     "Mig-quad"   = c(+2,-1,-2,-1, +2, 0, 0, 0, 0, 0),
                     "Mig-brk1"   = c(-2, 0, 1, 2, -1, 0, 0, 0, 0, 0),
                     "Mig-brk2"   = c(-1,-1,-1,-1,  5, 0, 0, 0, 0, 0),

                     "Res-linear" = c( 0, 0, 0, 0, 0,-2,-1, 0, 1, 2),
                     "Res-quad"   = c( 0, 0, 0, 0, 0,+2,-1,-2,-1,+2))

Run multcomp

mult.comp.working <- glht(m.Nvshab.ANOVA.means       
                          ,linfct =contrst.mat
                          ,alternative = c("two.sided"))


mult.comp.summary <- summary(mult.comp.working
                             ,test = adjusted(type = "none"))

mult.comp.summary

Plotting ANOVA

Calculate CIs

#merTools::predictInterval()
mod <- m.Nvshab.ANOVA.full
dat <- mod@frame
newdat <- expand.grid(
          i = 1
          #,band = dat$band[1]
          ,year = unique(dat$year)[1]
          #,site.age.cent = unique(dat$site.age.cent)
          ,spp.code = unique(dat$spp.code)[2]
          ,site = unique(dat$site)
          ,hab1 = unique(dat$hab1)
          )

names(dat)
names(newdat)
out <- predictInterval(mod, 
                       newdata =   newdat,
                which = "fixed", 
                level = 0.95,
                #n.sims = 10,
                stat = "median",
                include.resid.var = FALSE)

out <- cbind(newdat,out)

Plots

pd <- position_dodge(0.5)
ggplot(data = out,
        aes(y = exp(fit),
            x = site,
            color = hab1,
            group = hab1)) +
  geom_point(position = pd,
             size = 3,
             aes(shape = hab1)) +
  geom_line(position = pd)+
  xlab("Site") +
  ylab("Abundance") +
  geom_errorbar(position = pd,
                aes(ymax = exp(upr),
                  ymin = exp(lwr)),
              width = 0) 

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