knitr::opts_chunk$set(echo = TRUE)
Function that calls all libraries needed for analysis
source("./R/load_libraries.R") load_libraries()
load(file = "./data/ann_caps4age_ratio.RData")
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"))
dcast(data = ann_caps4age_ratio, formula = spp.code ~ age.AHY.01, fun.aggregate = length, value.var = "age.AHY.01")
dcast(data = ann_caps4age_ratio, formula = spp.code ~ sex.MF, fun.aggregate = length, value.var = "age.AHY.01")
"sex" column current coded
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)
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)
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)
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]
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")
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() )
Aceitillar models seperatealy b/c it does not have a value for "site.age"
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") )
#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)
convergence issues
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
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.
"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])
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)))
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)
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")
#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())
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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.