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

knitr::opts_chunk$set(echo = TRUE)

fit spp/site as fixed effect

# 
# m.age.spp.code <- bglmer(as.numeric(age.AHY.01)-1 ~ -1 + spp.code +
#                (1|band) +
#                (1|year) , #+
#                #(1|site) +
#               # (1|spp.code), # +
#               # (1|spp_site),
#               data = ann_caps2[,],
#               family = binomial#,
#               #glmerControl(optimizer = "Nelder_Mead")
#               )

# m.age.spp.codeMENAS <- bglmer(as.numeric(age.AHY.01)-1 ~ -1 + spp.code:site +
#                (1|band) +
#                (1|year) , #+
#                #(1|site) +
#               # (1|spp.code), # +
#               # (1|spp_site),
#               data = ann_caps2[,],
#               family = binomial#,
#               #glmerControl(optimizer = "Nelder_Mead")
#               )

fit spp/site as random effects

# m.age.site.fixed <- bglmer(as.numeric(age.AHY.01)-1 ~ 1 +
#                (1|band) +
#                (1|year) +
#                #(1|site) +
#                (1+ site|spp.code),
#               data = ann_caps2[,],
#               family = binomial#,
#               #glmerControl(optimizer = "Nelder_Mead")
#               )

#create unique spp-site combo
ann_caps2$spp_site <- with(ann_caps2,paste(spp.code, site))

#cente site age variable
ann_caps2$site.age.cent <- scale(ann_caps2$site.age,scale = F)


m.age.site.fixed.alt <- bglmer(as.numeric(age.AHY.01)-1 ~ 1 +
               (1|band) +
               (1|year) +
               (1|site) +
               (1|spp.code) +
               (1|spp_site),
              data = ann_caps2[,],
              family = binomial#,
              #glmerControl(optimizer = "Nelder_Mead")
              )


#se.coef(m.age.site.fixed)["spp.code"]
# mod.dat <- m.age.site.fixed.alt@frame
# 
# dim(mod.dat)
# 
# unique(spp_site)
# 
# i <- match(unique(mod.dat$band),mod.dat$band)
# 
# newdat <- expand.grid(year = unique(mod.dat$year),
#                      spp.code = unique(mod.dat$spp.code),
#                      site = unique(mod.dat$site))
# 
# predictInterval(m.age.site.fixed,
#                 as.matrix(newdat),
#                 which = "random", 
#                 level = 0.95,n.sims =100)

bootstrap

m <- m.age.site.fixed.alt

#ranef(m.age.site.fixed.alt)[["spp.code"]][,1]

summ.blups <- function(x){
  out <- c(ranef(x)[["spp_site"]][,1],
         ranef(x)[["spp.code"]][,1],
         ranef(x)[["site"]][,1])
  return(out)
  print(Sys.time())
  }


bb$t

## run of 10
#load(file = "./data/bb.RData")

library(lme4)
bb2 <- bootMer(m.age.site.fixed.alt,
              FUN=summ.blups,
              nsim=250,
              #use.u=FALSE,
              parallel = "multicore")

save(bb2, file = "./data/bb2.RData")

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