knitr::opts_chunk$set(echo = TRUE)
The script carries out an analysis of the abundance trends of each species using a regression style GLMM and extracting random slopes for each species. It loads a dataframe with the following columns
"year.num" "site" "spp.code" "hab1" "diet"
"year" "site.age.cent" "status.focals" "N" "i"
"site.age.init" "site.age"
load(file = "./data/ann_counts.RData") #loads "ann_counts"
levels(ann_counts$site) summary(net_hrs$year) summary(factor(ann_counts$year)) intersect(names(net_hrs), names(ann_counts))
Function that calls all libraries needed for analysis
source("./R/load_libraries.R") load_libraries()
names(ann_counts)
Aceitillar is run seperately using 1-way ANOVA style models in the next script in this series.
i.N.aceit <- which(ann_counts$site == "Aceitillar")
m.N.rand.slopes <- bglmer(N ~ 1 + (1|i) + (1|year) + (1|site) + (site.age.cent|spp.code), offset = log(net.hours), data = ann_counts[-i.N.aceit, ], family = poisson , glmerControl(optimizer = "Nelder_Mead") )
pois.blup.out <- data.frame(beta = ranef(m.N.rand.slopes)[["spp.code"]]$site.age.cent, se = se.coef(m.N.rand.slopes)["spp.code"][["spp.code"]][,"site.age.cent"]) pois.blup.out$Species <- row.names(pois.blup.out)
Set up names for labeling
focal.mig <- c("OVEN","BAWW","COYE","AMRE","CMWA", "BTBW","PAWA","PRAW", "MAWA","NOPA","PALM") #add to poisson focal.res <- c("HLCU","HILC", #fix - which name to use? "STOF", "RLTH","NOMO", "GRWA", "GTGT", #fix - which name to use? "BANA","BCPT", "YFGR","BFGR", "GABU", "BWVI", #32; Aceitillar / Cueva +misc "GREP") #18 just Aceitillar - ever yr focals <- c(focal.mig,focal.res)
Label residents vs. migrants
pois.blup.out$Species <- factor(pois.blup.out$Species) pois.blup.out$Status <- "Resident" pois.blup.out$Status[which(pois.blup.out$Species %in% focal.mig)] <- "Migrant"
To Do: -add endemic status -plot against observed abundance -NOTE: several spp. have very non-linear pattern! -sort by magnitude -add info on CIs
ggplot(dat = pois.blup.out, aes(y = exp(beta), x = Species, color = Status)) + 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_point(size = 4) + geom_hline(yintercept = exp(0)) + coord_flip() + facet_wrap(~Status,scales = "free") + ylab("Trend")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.