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

Dataframe contains the following variables for each site and year

• simpson's index
• spp.rich
• eveness

"dominance" also reported in the original MS but not calculate as of 2/12/18

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

Plot data

# library(ggpubr)
# 
# ggplot(data = community_dat,
#             aes(y = spp.rich,
#             x = site.age)) +
#   geom_point() +
#   geom_smooth(method = "lm")

Model

Model Pasture sites

i.aceit <- which(community_dat$site == "Aceitillar")
community_dat$i <- 1:dim(community_dat)[1]

#null model
m.spp.rich.0 <- bglmer(spp.rich ~ 1 + 
                        (1|i) +
                       (1|site) +
                       (1|year),
                     family = "poisson",
                     data = community_dat[-i.aceit,])

#site model
m.spp.rich.site <- update(m.spp.rich.0, . ~ . + site.age.cent)

#results
anova(m.spp.rich.0,
      m.spp.rich.site)

Model Aceitillar alone

This was used in previous verions but has been replaced by output from 1-way ANOVA style model.

m.aceit.spp.rich.0 <- bglmer(spp.rich ~ 1 + 
                        (1|i) +
                        #(1|site) +
                        (1|year),
                     family = "poisson",
                     data = community_dat[i.aceit,])

1-way ANOVA style model

Remove overlap

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


community_dat_2 <- community_dat[-c(i.cueva.drop.comm,
                                    i.caoba.drop.comm), ]

Run models

#null - w/ ranef for site
m.anova.spp.rich.0 <- bglmer(spp.rich ~ 1 + 
                       (1|i) +
                       (1|site) +
                       (1|year),
                     family = "poisson",
                     data = community_dat_2[,])

# site
m.anova.spp.rich.site <- bglmer(spp.rich ~ site + 
                       (1|i) +
                       #(1|site) +
                       (1|year),
                     family = "poisson",
                     data = community_dat_2[,])
# means
m.anova.spp.rich.means <- bglmer(spp.rich ~ -1 + site + 
                       (1|i) +
                       #(1|site) +
                       (1|year),
                     family = "poisson",
                     data = community_dat_2[,])

Inference

Model comparisons

ICtab(m.anova.spp.rich.0,
      m.anova.spp.rich.site,
      type = "AICc",
      weights = T,
      delta = T)

Contrast modeling

Create contrast matrix

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

Run multcomp

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


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

Mult comp output

mult.comp.summary.spp.rich

Plotting Species Richness

Plot pasture sites

Calculate CIs for pasture sites

#merTools::predictInterval()
source('~/1_R/git/git-aviary/DRmencia/R/build_CIs_regrxn_pois_comm.R')
spp.rich.rgrxn.preds <- buidCIs_regrxn_pois_comm(m.spp.rich.site)

Determine mean of site.age variable to de-center it while plotting

covar.mean <- mean(community_dat$site.age.init,
                   na.rm = TRUE)

Calculate CIs for ANOVA

source('~/1_R/git/git-aviary/DRmencia/R/buidCIs_ANOVA_pois_comm.R')
spp.rich.ANOVA.preds <- buidCIs_ANOVA_pois_comm(mod = m.anova.spp.rich.site)

Build regression Plot

source('~/1_R/git/git-aviary/DRmencia/R/back_transform_reals.R')
source('~/1_R/git/git-aviary/DRmencia/R/plot_regrxn_pois_comm.R')
spp.rich.rgrxn.preds <- back_transform_reals(spp.rich.rgrxn.preds)
spp.rich.ANOVA.preds <- back_transform_reals(spp.rich.ANOVA.preds)


regrxnplot_pois_spprich <- plot_regrxn_pois_comm(dat = spp.rich.rgrxn.preds,
                      ANOVA.dat = spp.rich.ANOVA.preds, 
                      mod = m.spp.rich.site,
                      mean.of.covar = covar.mean,
                      ylab. = "Species Richness",
                      ymin. = 0,
                      ymax. = 35)

Build Aceitillar plot

i.aceit <- which(spp.rich.ANOVA.preds$site == "Aceitillar")
i.aceit.mod <- which(m.anova.spp.rich.site@frame$site == "Aceitillar")
source('~/1_R/git/git-aviary/DRmencia/R/plot_ANOVA_pois_comm.R')
anovaplot_pois_spprich <- plot_ANOVA_pois_comm(spp.rich.ANOVA.preds[i.aceit,],
                                               mod = m.anova.spp.rich.site, 
                                               mod.dat.i = i.aceit.mod,
                                               ylab. = "",
                                               point.sz = 4,
                                               ymin. = 0,
                                               ymax. = 35,
                                               xlab. = "")

Save constituent plots

save.image("./plots/regrxnplot_pois_spprich.RData")
save.image("./plots/anovaplot_pois_spprich.RData")

Plot both Regression and ANOVA together

cowplot::plot_grid(regrxnplot_pois_spprich,
                   anovaplot_pois_spprich,
                   rel_widths = c(8,1))

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