knitr::opts_chunk$set(echo = TRUE)
load(file = "./data/ann_counts") #loads "ann_counts"
Function that calls all libraries needed for analysis
source("./R/load_libraries.R") load_libraries()
names(ann_counts)
i.N.aceit <- which(ann_counts$site == "Aceitillar")
Aceitillar models seperatealy b/c it does not have a value for "site.age"
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") )
#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)
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]]
#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]]
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")
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)
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")
#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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.