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

knitr::opts_chunk$set(echo = TRUE)

Preliminaries

Load cleaned, merged and scrubbed data

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

Load libraries

Function that calls all libraries needed for analysis

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

Look at data

names(ann_counts)

Index for Removing aceitillar

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

1) 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"

Model random slopes

Full random slopes model

m.N.rand.slopes <- bglmer(N ~ 1 +
               (1|i) + 
               (1|year) +
               (1|site) +
               #(1|site:spp.code) #this would be ideal but prob overkill
               (site.age.cent|spp.code),
              data = ann_counts[-i.N.aceit, ],
              family = poisson ,
              glmerControl(optimizer = "Nelder_Mead")
              )

Null model

Model fixed effects

#status fixed effect
m.N.fixef.stat     <- update(m.N.rand.slopes, . ~ . + status.focals)

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


anova(m.N.fixef.stat,
      m.N.fixef.site.age)

Addive model + status.focals + site.age.cent

Convergence issues

#convergence issues
## trows gradient warning
m.N.fixef.add <- update(m.N.rand.slopes, . ~ . + status.focals + site.age.cent)




m.N.fixef.add.ALL <- afex::all_fit(m.N.fixef.add)
str(m.N.fixef.add.ALL,1)
summary(m.N.fixef.add.ALL$bobyqa.)
summary(m.N.fixef.add.ALL[[2]])
summary(m.N.fixef.add.ALL[[3]])
summary(m.N.fixef.add.ALL[[4]])
summary(m.N.fixef.add.ALL[[5]])
summary(m.N.fixef.add.ALL[[6]])
summary(m.N.fixef.add.ALL[[7]])


m.N.fixef.add.bobyqa <- update(m.N.rand.slopes, . ~ . + status.focals + site.age.cent, glmerControl(optimizer = "bobyqa"))


## check gradient
# mod2check <- m.N.fixef.add.ALL[[7]]
# summary(mod2check)
# derivs1 <- mod2check@optinfo$derivs
# sc_grad1 <- with(derivs1,solve(Hessian,gradient))
# max(abs(sc_grad1)) #0.00227505
# 
# 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)))     


m.N.fixef.add <- m.N.fixef.add.ALL[[7]]

Full model

#convergence issues
##fails with bobyqa under ALL conditions
##fails ...
##works: w/ Nelder_Mean when only intial captures are used 
##       
m.N.fixef.full <- update(m.N.rand.slopes, . ~ . + status.focals*site.age.cent)

m.N.fixef.full.ALL <- afex::all_fit(m.N.fixef.full)


summary(m.N.fixef.full.ALL$bobyqa.)
summary(m.N.fixef.full.ALL[[2]])
summary(m.N.fixef.full.ALL[[3]])
summary(m.N.fixef.full.ALL[[4]])
summary(m.N.fixef.full.ALL[[5]])
summary(m.N.fixef.full.ALL[[6]])
summary(m.N.fixef.full.ALL[[7]])


## check gradient
mod2check <- m.N.fixef.full.ALL[[7]]
summary(mod2check)
derivs1 <- mod2check@optinfo$derivs
sc_grad1 <- with(derivs1,solve(Hessian,gradient))
max(abs(sc_grad1)) # 0.0002914756

m.N.fixef.full <- m.N.fixef.full.ALL[[7]]

Inference

AIC

ICtab(m.N.rand.slopes,
       m.N.fixef.stat,
       m.N.fixef.site.age, #conv warn
       m.N.fixef.add,      #fails
       m.N.fixef.full,
       base = TRUE,
       logLik = TRUE,
       type = "AICc")

Likelihood ratio test

anova(m.N.fixef.full,
      m.N.fixef.add)

anova(m.N.fixef.add,
      m.N.fixef.site.age)

anova(m.N.fixef.site.age,
      m.N.fixef.stat)

library(lmerTest)
lmerTest::anova(m.N.fixef.full)

anova(m.N.fixef.full,
      m.N.fixef.add)

Save model model

ind.spp.abund.trends <- list(m.N.rand.slopes =m.N.rand.slopes,
                             m.N.fixef.stat = m.N.fixef.stat,
                           m.N.fixef.site.age = m.N.fixef.site.age,
                           m.N.fixef.add= m.N.fixef.add,
                           m.N.fixef.full = m.N.fixef.full)

save(ind.spp.abund.trends, 
           file = "./models/ind_spp_abund_trends.RData")

Plotting

Calculate CIs

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


out <- predictInterval(m.N.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_caps$site.N.init,na.rm = TRUE)

Plots

gg.mencia <- ggplot(data = out,
        aes(y = invlogit(fit),
            x = site.age.cent+7.557782,
            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")

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