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

knitr::opts_chunk$set(echo = TRUE)

Libraries

Function that calls all libraries needed for analysis

source("./R/load_libraries.R")
load_libraries()

Load data

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

Plot data

library(ggpubr)

ggplot(data = community_dat,
            aes(y = exp(simpson),
            x = site.age)) +
  geom_point() +
  geom_smooth(method = "lm")

Model

i.aceit <- which(community_dat$site == "Aceitillar")
community_dat$i <- 1:dim(community_dat)[1]
m.simpson.0 <- blmer(exp(simpson) ~ 1 + 
                      #  (1|i) +
                       (1|site) +
                       (1|year),
                     #family = "poisson",
                     data = community_dat[-i.aceit,])

m.simpson.site <- update(m.simpson.0, . ~ . + site.age.cent)

anova(m.simpson.0,
      m.simpson.site)

Plotting Species Richness

Calculate CIs

#merTools::predictInterval()
mod <- m.simpson.site

names(mod@frame)

newdat <- expand.grid(
          #i = mod@frame$i[1]
          site = unique(mod@frame$site)[1]
          ,year = unique(mod@frame$year)[1]
          ,site.age.cent = unique(mod@frame$site.age.cent)
          )


out <- predictInterval(mod, 
                       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_caps4sex_ratio$site.age.init,
                   na.rm = TRUE)

Plots

gg.sppdiv<- ggplot(data = out,
        aes(y = fit,
            x = site.age.cent+covar.mean)) +
  geom_line(aes(),
            size = 1) +
  xlab("Site age") +
  ylab("exp(simpson)") +
  geom_ribbon(aes(ymax = upr,
                  ymin = lwr),
              alpha = 0.125,
              linetype = 0)  + theme(legend.position="none") +
  geom_point(data = mod@frame,
             aes(y = mod@frame[,1],
                 x = site.age.cent+covar.mean))

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