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

knitr::opts_chunk$set(echo = TRUE)

Libraries

#base glmm function
library(lme4)

#allows models with mildly informative priors
## improves convergned
library(blme)

# for checking convergence warnings from lme4
library(numDeriv)

library(afex)

# AIC calculation
library(bbmle)

# Estiamted SEs of BLUPs
library(arm)

# Confidence intervals for predicitons from lme4
library(merTools)

# Plotting
library(ggplot2)
library(cowplot)

# Plotting glmms
library(sjPlot)

Load data

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

load(file = "./data/ann_caps4age_ratio.RData")

Data subset

Remove Aceitillar b/c is does not have an age assigned; Aceitillar is analysed is sperate set of scripts

i.aceit <- which(ann_caps4age_ratio$site == "Aceitillar")

Below further species are removed because they are rare at Aceitillar; 2001-2002 data also removed

Model Aceitillar

Prep

ann_caps4age_ratio$year <- factor(ann_caps4age_ratio$year)
ann_caps4age_ratio$band <- factor(ann_caps4age_ratio$band)

summary(ann_caps4age_ratio[i.aceit,c("band","year","site","spp.code")],10)

Check average rates using GLM

An initial workup of data had problems with band coding for the Aceitillar that originated with a problem w/ Excel converting bands to dates. These GLMs were used as comparisons to the GLMMs run below to validate that something was wrong.

invlogit(coef(glm(age.AHY.01 ~ 1, data = ann_caps4age_ratio[i.aceit,],
              family = binomial)))

invlogit(coef(glm(age.AHY.01 ~ -1 + status.focals, data = ann_caps4age_ratio[i.aceit,],
              family = binomial)))

invlogit(coef(glm(age.AHY.01 ~ -1 + year, data = ann_caps4age_ratio[i.aceit,],
              family = binomial)))

invlogit(coef(glm(age.AHY.01 ~ -1 + year:status.focals, data = ann_caps4age_ratio[i.aceit,],
              family = binomial)))


invlogit(coef(glm(age.AHY.01 ~ -1 + spp.code, data = ann_caps4age_ratio[,],
              family = binomial)))

Subset just Aceitillar

ann_caps4age_ratio.aceit <- ann_caps4age_ratio[i.aceit,]

ID 2001 to 2002

Captures very rare

i.0102 <- which(ann_caps4age_ratio.aceit$year == "2001-2002")

length(i.0102)

ann_caps4age_ratio.aceit2 <- ann_caps4age_ratio.aceit[-i.0102,]

# 
# ann_caps4age_ratio.aceit2$year <- as.character(ann_caps4age_ratio.aceit2$year)
# ann_caps4age_ratio.aceit2$spp.code <- as.character(ann_caps4age_ratio.aceit2$spp.code)
# with(ann_caps4age_ratio.aceit2[-i.0102,],
#      table(spp.code,year))

How many cpatures per species

dcast(data = ann_caps4age_ratio.aceit2,
      formula = spp.code ~ age.AHY.01,
      fun.aggregate = length,
      value.var = "age.AHY.01")

Rare at Aceitillar

rare <- c("BFGR","HILC","NOMO","YFGR","GREP",
"COYE","PAWA","PRAW")

i.rare <- which(ann_caps4age_ratio.aceit2$spp.code %in% rare)

ann_caps4age_ratio.aceit3 <- ann_caps4age_ratio.aceit2[-i.rare,]

GLMM for Aceitillar

Null model

m.age.fixef.aceit.0 <- bglmer(age.AHY.01 ~ 1 +
               (1|band) +
               (1|year) +
               #(1|site) + #not applicable
               #(1|site:spp.code) #not applicable
               (1|spp.code), #site.age.cent not applicable
              data = ann_caps4age_ratio.aceit3,
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

Check results

#fixed effects
invlogit(fixef(m.age.fixef.aceit.0))

#random effects for each spp (BLUPs)
invlogit(ranef(m.age.fixef.aceit.0)[[2]])

fixed effeccts: status.focals

Model effect of status

#status.focals
m.age.fixef.aceit.stat <- update(m.age.fixef.aceit.0, . ~ . + status.focals)

#mean parameterization
m.age.fixef.aceit.means <- bglmer(age.AHY.01 ~ - 1 + status.focals+
               (1|band) +
               (1|year) +
               #(1|site) + #not applicable
               #(1|site:spp.code) #not applicable
               (1|spp.code), #site.age.cent not applicable
              data = ann_caps4age_ratio.aceit3,
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

Check output

invlogit(fixef(m.age.fixef.aceit.means))
invlogit(ranef(m.age.fixef.aceit.means)[[2]])

Inference

ICtab(m.age.fixef.aceit.0,
       m.age.fixef.aceit.stat,
       base = TRUE,
       logLik = TRUE,
       type = "AICc")

Calculate CIs

#merTools::predictInterval()

newdat.aceit <- expand.grid(
          band = m.age.fixef.aceit.stat@frame$band[1]
          ,year = unique(m.age.fixef.aceit.stat@frame$year)[1]
          #,site = unique(m.age.fixef.aceit.stat@frame$site)[1]
          #,site.age.cent = unique(m.age.fixef.aceit.stat@frame$site.age.cent)
          ,spp.code = unique(m.age.fixef.aceit.stat@frame$spp.code)[1]
          ,status.focals = unique(m.age.fixef.aceit.stat@frame$status.focals)
          )

names(m.age.fixef.aceit.stat@frame)
out.aceit <- predictInterval(m.age.fixef.aceit.stat, 
                       newdata =   newdat.aceit,
                which = "fixed", 
                level = 0.95,
                n.sims = 10,
                stat = "mean",
                include.resid.var = FALSE)

out.aceit <- cbind(newdat.aceit,out.aceit)

Graph

names(out.aceit) <- gsub("status.focals","Residency",names(out.aceit))
gg.aceit <- ggplot(data = out.aceit,
        aes(y = invlogit(fit),
            x = Residency,
            color = Residency,
            shape = Residency)) +
  geom_point(size = 3)+
  xlab("Residency") +
  ylab("Age ratio (AHY vs all ages)") +
  geom_errorbar(aes(ymax = invlogit(upr),
                  ymin = invlogit(lwr),
                  width = 0)) +
  ylim(0,1) +
  ylab("")

Plot graph of both Mencia and Aceitillar

plot_grid(gg.mencia, gg.aceit, labels = c('Mencia', 'Aceitillar'),
          rel_widths = c(2,1.5),
          label_x = c(0.125,0.125))

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