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

knitr::opts_chunk$set(echo = TRUE)

Load data

Data has just focal species; rarities removed
Within-year repatures removed

load(file = "./data/ann_caps.RData")

Fix missing mig/res status for rare spp

ann_caps[which(is.na(ann_caps$status.focals) == TRUE), c("spp.code")]

ann_caps$status.focals[which(ann_caps$spp.code == "LOST")] <- NA

ann_caps$status.focals[which(ann_caps$spp.code == "WOTH")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "PALM")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "MAWA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "INBU")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "TEWA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "NOPA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "LOWA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "BWWA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "HOWA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "GRCA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "RBGR")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "SUTA")] <- "mig"
ann_caps$status.focals[which(ann_caps$spp.code == "YRWA")] <- "mig"


ann_caps$status.focals[which(ann_caps$spp.code == "YBCU")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "SSHA")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "SBMU")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "SBAN")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "RTSO")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "NOPO")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "HITR")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "GRAK")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "WFQD")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "ZEND")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "HIWO")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "MACU")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "AMKE")] <- "res"

ann_caps$status.focals[which(ann_caps$spp.code == "HIPE")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "KWQD")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "RUQD")] <- "res"

ann_caps$status.focals[which(ann_caps$spp.code == "BWVI")] <- "res"

ann_caps$status.focals[which(ann_caps$spp.code == "HISP")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "COGD")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "ANPI")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "BLGR")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "BLPW")] <- "res"
ann_caps$status.focals[which(ann_caps$spp.code == "BUOW")] <- "res"

summary(factor(ann_caps$status.focals))

Label endemics

unique(ann_caps$spp.code[which(ann_caps$status.focals == "res")] )

library(wrapr)
endems <- wrapr::qc(BBTO,BCPT,FBVI,GTGT,HHTA,HIEM,HILC,
HIOR,HIPE,HIPK,HISP,HITR,HIWO,LATH,NBTO,PALM,
WCHT,WFQD)

ann_caps$endemicYN <- "not.endem"

ann_caps$endemicYN[which(ann_caps$spp.code %in% endems)] <- "endem"

i.cueva.drop <- which(ann_caps$site == "La Cueva" & 
                      ann_caps$site.age > 4) #keep age 2, 3, 4
i.caoba.drop <- which(ann_caps$site == "La Caoba" & 
                      ann_caps$site.age < 7) #keep age 7, 8, 9

ann_caps_2 <- ann_caps[-c(i.cueva.drop,i.caoba.drop), ]

Order factors

summary(ann_counts$site)
ann_caps$site <- factor(ann_caps$site,
                        levels =c("La Cueva","La Caoba",
                                  "Morelia","El Corral","Aceitillar"))

ann_caps_2$site <- factor(ann_caps_2$site,
                        levels =c("La Cueva","La Caoba",
                                  "Morelia","El Corral","Aceitillar"))

Reshape to abundance - by residency status

tot_caps_by_status <- dcast(data = ann_caps,
                            formula = year + year.num + site.age + 
                              site  ~ status.focals,
                            #value.var = "spp.code",
                            fun.aggregate = length)

Merge with net hours

net_hours <- read.csv("./data/net_hours.csv")
tot_caps_by_status2 <- merge(tot_caps_by_status,
                             net_hours, all = T)

dim(tot_caps_by_status)
dim(tot_caps_by_status2)

Plot as counts

i.aceit2 <- which(tot_caps_by_status2$site == "Aceitillar")
ggplot(data = tot_caps_by_status2[,],
       aes(y = mig/(res+mig),
           x = site)) + #site.age
  geom_point(aes(color = year)) +
  geom_smooth(method = "lm", se= FALSE)

Model as counts

Model migrants vs. Residents

use ann_caps not ann_caps2 b/c regression can accomodate the overlap in ages of the 2 youngest sites

i.aceit <- which(ann_caps$site == "Aceitillar")
summary(ann_caps[i.aceit,"site.age.cent"])

ann_caps$status.focals <- factor(ann_caps$status.focals,
                                 levels = c("res","mig"))
summary(ann_caps$status.focals)

m.mig.reg.0 <- bglmer(status.focals ~ 1 +
               (1|band) + 
               (1|year) +
               (1|site),
              data = ann_caps[-i.aceit,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )
m.mig.reg.x <- update(m.mig.reg.x, .~ . + site.age.cent)
m.mig.reg.x2 <- update(m.mig.reg.x, .~ . + I(site.age.cent^2))
summary(m.mig.reg.x)
summary(m.mig.reg.x2)

anova(m.mig.reg.0,
      m.mig.reg.x)
anova(m.mig.reg.x,
      m.mig.reg.x2)

The overall proportion of captures which were migrants increased as pastures aged (slope = 0.01, SE = 0.03) but the trend was not significant (Likelihood ratio test: chi2 = 0.25, df = (4,5), p = 0.61). Capture rates of migrants decline from the second oldest to the oldest pasture but there was no signficant quadratic effect in this rergession (Likelihood ratio test: chi2 = 0, df = (5,6), p = 1). The proportion of migrants in the mature forest was lower than any other site (0.15), resulting in a signifcant negative quadratic trend accross the chronosequence (z = -5.25, p < 0.001).

Anova models w/ all sites

note: don't use species as a random effect! want overall rate. (probably wont work anyway)

ann_caps_2$status.focals <- factor(ann_caps_2$status.focals,
                                   levels = c("res","mig"))
summary(ann_caps_2$status.focals)

m.mig.anova <- bglmer(status.focals ~ site +
               (1|band) + 
               (1|year) ,
              data = ann_caps_2[,], #ann_caps_2 set up so no overlap of young sites
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )
summary(m.mig.anova)
m.mig.anova.means <- bglmer(status.focals ~ -1 + site +
               (1|band) + 
               (1|year) ,
              data = ann_caps_2[,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

summary(m.mig.anova.means)

m.mig.anova.means
invlogit(fixef(m.mig.anova.means))

Contrast modeling

Create contrast matrix

contrst.mat.ratios <- rbind("Linear" = c(-2,-1, 0, 1, 2),
                          "Quad"   = c(+2,-1,-2,-1,+2))

Run multcomp

mult.comp.working.mig <- glht(m.mig.anova.means       
                          ,linfct =contrst.mat.ratios
                          ,alternative = c("two.sided"))


mult.comp.summary.mig <- summary(mult.comp.working.mig
                             ,test = adjusted(type = "none"))

Mult comp output

mult.comp.summary.mig

Model endemics

i.res <- which(ann_caps$status.focals == "res")

ann_caps$endemicYN <- factor(ann_caps$endemicYN,
                                 levels = c("not.endem","endem"))
summary(factor(ann_caps$endemicYN))

m.endem <- bglmer(endemicYN ~ site.age.cent +
               (1|band) + 
               (1|year) +
               (1|site),
              data = ann_caps[i.res,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

summary(m.endem)

The propotion of individuals that were endemics increased on average as pastures aged but the change was not significant (beta = 0.06, SE = 0.04, z = 1.4, p = 0.15). Across the entire chronosequence there was a significant linear trend in the frequency of endemic individuals (z = 5.28, p < 0.001). THough their community composition was different as disucssed above the frequency of endemic indivdiuals was similar in the oldest pasture (25%) and the mature forest (23%).

ann_caps_2$endemicYN <- factor(ann_caps_2$endemicYN,
                                 levels = c("not.endem","endem"))
m.endem.anova <- bglmer(endemicYN ~ site +
               #(1|band) + 
               (1|year) ,
              data = ann_caps_2[,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )
summary(m.endem.anova)

m.endem.anova.means <- bglmer(endemicYN ~ -1 + site +
               #(1|band) + 
               (1|year) ,
              data = ann_caps_2[,],
              family = binomial,
              glmerControl(optimizer = "Nelder_Mead")
              )

summary(m.endem.anova.means)

Contrast modeling

Run multcomp

mult.comp.working.endem <- glht(m.endem.anova.means       
                          ,linfct =contrst.mat.ratios
                          ,alternative = c("two.sided"))


mult.comp.summary.endem <- summary(mult.comp.working.endem
                             ,test = adjusted(type = "none"))

Mult comp output

mult.comp.summary.endem

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

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

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

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

Load data

Fix missing mig/res status for rare spp

Label endemics

Order factors

Reshape to abundance - by residency status

Merge with net hours

Plot as counts

Model as counts

Model migrants vs. Residents

Anova models w/ all sites

Contrast modeling

Model endemics

Contrast modeling

R Package Documentation

Browse R Packages

We want your feedback!

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

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

Load data

Fix missing mig/res status for rare spp

Label endemics

Order factors

Reshape to abundance - by residency status

Merge with net hours

Plot as counts

Model as counts

Model migrants vs. Residents

Anova models w/ all sites

Contrast modeling

Model endemics

Contrast modeling

R Package Documentation

Browse R Packages

We want your feedback!

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