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

knitr::opts_chunk$set(echo = TRUE)

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_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" &
                   ann_caps4age_ratio$spp.code %in% c("BANA",
                                                      "STOF",
                                                      "RLTH",
                                                      "NOMO",

                                                      "PAWA",
                                                      "OVEN",
                                                      "BCPT",
                                                      "BWVI"))

Explore data

Number of individuals per age class

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

Number of individuals per sex

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

Recode sex variable to be binary

"sex" column current coded

• U = unknown / can't be sexed
• F = female
• M = male

Recode to be NA, F, M

ann_caps4age_ratio$sex.MF <- NA
ann_caps4age_ratio$sex.MF[which(ann_caps4age_ratio$sex == "M")] <- "M"
ann_caps4age_ratio$sex.MF[which(ann_caps4age_ratio$sex == "F")] <- "F"

ann_caps4age_ratio$sex.MF <- factor(ann_caps4age_ratio$sex.MF)

Random slopes model - age ratio

Model

Extract blups from model w/o age

summary(ann_caps4age_ratio$age.AHY.01) #age2 is 2nd level
summary(ann_caps4age_ratio$status.focals)

ann_caps4age_ratio$status.focals <- factor(ann_caps4age_ratio$status.focals,
                                           levels = c("res","mig"))

m.age.rand.slopes <- bglmer(age.AHY.01 ~ 1 +
               (1|band) + 
               (1|year) +
               (1|site) +
               #(1|site:spp.code) #this would be ideal but cause converge problems
               (site.age.cent|spp.code),
              data = ann_caps4age_ratio[-i.aceit,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

summary(m.age.rand.slopes)

Extract blups from model WITH SEX

Estimate sep. slopes for each sex. THis model fitting performs very poorly and was not used in the final MS>

with(ann_caps4age_ratio[-i.aceit], 
     table(age.AHY.01,sex.MF,spp.code))
m.age.rand.slopes.with.sex <- bglmer(age.AHY.01 ~ 1 +
               (1|band) + 
               (1|year) +
               (1|site) +
               (site.age.cent:sex.MF|spp.code),
              data = ann_caps4age_ratio[-i.aceit,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

summary(m.age.rand.slopes.with.sex)

Extract blups and calculate SEs

extract BLUPs for model w/o age

out.blups <- data.frame(beta = ranef(m.age.rand.slopes)[["spp.code"]]$site.age.cent,
se  = se.coef(m.age.rand.slopes)["spp.code"][["spp.code"]][,"site.age.cent"])

row.names(out.blups) <- gsub("[ ][LE][al]","",row.names(out.blups))
ids <- stringi::stri_split(row.names(out.blups),regex= " ", simplify = T)

out.blups$Species <- ids[,1]
out.blups$Site <- ids[,1]

Extract blups from model WITH age

mod <- m.age.rand.slopes.with.sex

age.blups.with.sex <- data.frame(
Species = rownames(ranef(mod)[["spp.code"]]),
Sex = c(rep("Female",dim(ranef(mod)[["spp.code"]])[1]),
        rep("Male",dim(ranef(mod)[["spp.code"]])[1])),
beta = c(ranef(mod)[["spp.code"]][,2],
              ranef(mod)[["spp.code"]][,3]),
se  = c(se.coef(mod)["spp.code"][["spp.code"]][,2],
        se.coef(mod)["spp.code"][["spp.code"]][,3])
)

age.blups.with.sex$Age <- as.character(age.blups.with.sex$Age)
age.blups.with.sex$Age <- ifelse(age.blups.with.sex$Age == "age1","HY","AHY")

Plot random slopes

From model m.age.rand.slopes

focal.mig <- c("OVEN","BAWW","COYE","AMRE","CMWA",
               "BTBW","PAWA","PRAW")

out.blups$Species <- factor(out.blups$Species,
                      levels = focals)
out.blups$Status <- "Resident"
out.blups$Status[which(out.blups$Species %in% focal.mig)] <- "Migrant"

gg.mencia.age.blups <- ggplot(dat = out.blups,
       aes(y = exp(beta),
           x = Species,
           color = Status,
           shape  = Status)) +
  geom_point(size = 4) +
   geom_errorbar(aes(ymax = exp(beta+1*se),
                    ymin = exp(beta-1*se)),
                width = 0,
                size = 2) +
  geom_errorbar(aes(ymax = exp(beta+1.96*se),
                    ymin = exp(beta-1.96*se)),
                width = 0) +
  geom_hline(yintercept = exp(0)) +
  coord_flip() +
  facet_wrap(~Status,scales = "free") +
  ylab("Trend % Male") +
  ylim(0.90,1.075) +
  theme(legend.position = "bottom", 
    legend.direction = "horizontal",
    legend.justification =  "center",
    legend.title=element_blank()
    )

Model overal trend by status (mig vs. resident)

post-pasture sites

Aceitillar models seperatealy b/c it does not have a value for "site.age"

Null model

Similar to above but required tweaking to converged. Current version should be identical but w/ different optimizer

summary(ann_caps4age_ratio$status.focals)

m.age.fixef.0 <- bglmer(as.numeric(age.AHY.01)-1 ~ 1 +
               (1|band) +
               (1|year) +
               (1|site) +
               #(1|site:spp.code) #this would be ideal but prob overkill
               (site.age.cent|spp.code),
              data = ann_caps4age_ratio[-i.aceit,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

Model fixed effects

#status fixed effect
## (almost?) always converged
m.age.fixef.stat     <- update(m.age.rand.slopes, . ~ . + status.focals)

#site.age.cent fixed effect
## throws warning
m.age.fixef.site.age <- update(m.age.rand.slopes, . ~ . +           site.age.cent)

#convergence issues
##fails: with bobyqa under most/all conditions
##works: w/ Nelder_Mean when only intial captures are used AND when all recaps used
### even with Nelder_Mean does throw warning
m.age.fixef.add <- update(m.age.rand.slopes, . ~ . + status.focals + site.age.cent)

Full model

convergence issues

• fails with bobyqa under ALL conditions
• works: w/ Nelder_Mean when only intial captures are used
• afex::all_fit() used to find best optimizr

m.age.fixef.full <- update(m.age.rand.slopes, . ~ . + status.focals*site.age.cent)

#full models throws warnings; run afex::all_fit()
m.age.fixef.full.allfit <- all_fit(m.age.fixef.full)

#list of output
str(m.age.fixef.full.allfit,1)
summary(m.age.fixef.full.allfit$Nelder_Mead.)
summary(m.age.fixef.full.allfit$optimx.nlminb)
summary(m.age.fixef.full.allfit[[4]])
summary(m.age.fixef.full.allfit$nloptwrap.NLOPT_LN_NELDERMEAD)#
summary(m.age.fixef.full.allfit$nloptwrap.NLOPT_LN_BOBYQA)#
summary(m.age.fixef.full.allfit$nmkbw.) #gradient good: 0.00455954

Check convergence warnings

See https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b615d8633c01d51.html

Also stackoverflow posts by Ben Bolker, eg https://stackoverflow.com/questions/21344555/convergence-error-for-development-version-of-lme4

use output of afex::all_fit to check different optimizers.

Check singularity

"The definition of singularity is that some of the constrained parameters of the random effects theta parameters are on the boundary (equal to zero, or very very close to zero, say <10e−6):"

tt <- getME(m.age.fixef.stat,"theta")
ll <- getME(m.age.fixef.stat,"lower")
min(tt[ll==0])

Check gradients

Warning was of the form "Model failed to converge with max|grad| = 0.00176624 (tol = 0.001, component 1)"

Worst was for m.age.fixef.full

mod2check <- m.age.fixef.full.allfit$nmkbw. 
summary(mod2check)
derivs1 <- mod2check@optinfo$derivs
sc_grad1 <- with(derivs1,solve(Hessian,gradient))
max(abs(sc_grad1)) #0.0007450714

max(pmin(abs(sc_grad1),abs(derivs1$gradient)))


dd <- update(mod2check,devFunOnly=TRUE)
pars <- unlist(getME(mod2check,c("theta","fixef")))
grad2 <- grad(dd,pars)
hess2 <- hessian(dd, pars)
sc_grad2 <- solve(hess2, grad2)
max(pmin(abs(sc_grad2),abs(grad2)))

Inference

LRT

Averaged accross species the age ratio of residents was constant (Fig. 4a) but

summary(m.age.fixef.full.allfit$nmkbw.)




anova(m.age.fixef.full.allfit$nmkbw.,
      m.age.fixef.add)

AICtab

ICtab(m.age.rand.slopes,
       m.age.fixef.stat,

       m.age.fixef.site.age, #conv warn
       m.age.fixef.add,      #fails
       m.age.fixef.full.allfit$nmkbw.,
       base = TRUE,
       logLik = TRUE,
       type = "AICc")

Plotting

Calculate CIs

#merTools::predictInterval()
m.age.fixef.full <- m.age.fixef.full.allfit$nmkbw.
newdat <- expand.grid(
          band = m.age.fixef.full@frame$band[1]
          ,year = unique(m.age.fixef.full@frame$year)[1]
          ,site = unique(m.age.fixef.full@frame$site)[1]
          ,site.age.cent = unique(m.age.fixef.full@frame$site.age.cent)
          ,spp.code = unique(m.age.fixef.full@frame$spp.code)[1]
          ,status.focals = unique(m.age.fixef.full@frame$status.focals)
          )

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

out <- cbind(newdat,out)

Determine mean of site.age variable to de-center it while plotting

covar.mean <- mean(ann_caps4age_ratio$site.age.init,na.rm = TRUE)

Plots

gg.mencia.mean.age <- ggplot(data = out,
        aes(y = invlogit(fit),
            x = site.age.cent+covar.mean,
            color = status.focals)) +
  geom_line(aes(linetype = status.focals),
            size = 1) +
  xlab("Site age") +
  ylab("Age ratio (AHY vs all ages)") +
  geom_ribbon(aes(ymax = invlogit(upr),
                  ymin = invlogit(lwr),
                  fill = status.focals),
              alpha = 0.125,
              linetype = 0) +
  ylim(0,1) + theme(legend.position="none") +
  theme(legend.position = "bottom", #legend.background = element_rect(color = "black", 
    #fill = "grey90", size = 1, linetype = "solid"), 
    legend.direction = "horizontal",
    legend.title=element_blank())

Build multi-paneled plot

cowplot::plot_grid(gg.mencia.mean.age,
                   gg.mencia.age.blups,
                   nrow = 1,rel_widths = c(1,1.75),
                   labels = c("a)","b)"))

save.image(file = "./models/WORKSPACE_age_ratio_regression.RData")

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