knitr::opts_chunk$set(echo = TRUE)
Function that calls all libraries needed for analysis
source("./R/load_libraries.R") load_libraries()
Dataframe contains the following variables for each site and year
"dominance" also reported in the original MS but not calculate as of 2/12/18
load(file = "./data/community_dat.RData")
# library(ggpubr) # # ggplot(data = community_dat, # aes(y = spp.rich, # x = site.age)) + # geom_point() + # geom_smooth(method = "lm")
i.aceit <- which(community_dat$site == "Aceitillar") #community_dat$i <- 1:dim(community_dat)[1] #null model m.div.0 <- blmer(exp(simpson) ~ 1 + #(1|i) + (1|site) + (1|year), #family = "poisson", data = community_dat[-i.aceit,]) #site model m.div.site <- update(m.div.0, . ~ . + site.age.cent) #results anova(m.div.0, m.div.site)
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), ]
#null - w/ ranef for site m.anova.div.0 <- blmer(exp(simpson) ~ 1 + #(1|i) + (1|site) + (1|year), #family = "poisson", data = community_dat_2[,]) # site m.anova.div.site <- blmer(exp(simpson) ~ site + #(1|i) + #(1|site) + (1|year), #family = "poisson", data = community_dat_2[,]) # means m.anova.div.means <- blmer(exp(simpson) ~ -1 + site + #(1|i) + #(1|site) + (1|year), # family = "poisson", data = community_dat_2[,])
ICtab(m.anova.div.0, m.anova.div.site, type = "AICc", weights = T, delta = T)
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.div <- glht(m.anova.div.means ,linfct =contrst.mat.comm ,alternative = c("two.sided")) mult.comp.summary.div <- summary(mult.comp.working.div ,test = adjusted(type = "none"))
Mult comp output
mult.comp.summary.div
#merTools::predictInterval() source('~/1_R/git/git-aviary/DRmencia/R/build_CIs_regrxn_pois_comm.R') div.rgrxn.preds <- buidCIs_regrxn_pois_comm(m.div.site,mod.type = "normal")
Determine mean of site.age variable to de-center it while plotting
covar.mean <- mean(community_dat$site.age.init, na.rm = TRUE)
source('~/1_R/git/git-aviary/DRmencia/R/buidCIs_ANOVA_pois_comm.R') div.ANOVA.preds <- buidCIs_ANOVA_pois_comm(mod = m.anova.div.site, mod.type = "normal")
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') div.rgrxn.preds <- back_transform_reals(div.rgrxn.preds,fnxn = "none") regrxnplot_pois_div <- plot_regrxn_pois_comm(dat = div.rgrxn.preds, ANOVA.dat = div.ANOVA.preds, mod = m.div.site, mean.of.covar = covar.mean, ylab. = "exp(Simpson's Diversity)", ymin. = 5, ymax = 20)
i.aceit <- which(div.ANOVA.preds$site == "Aceitillar") i.aceit.mod <- which(m.anova.div.site@frame$site == "Aceitillar") div.ANOVA.preds <- back_transform_reals(div.ANOVA.preds,fnxn = "none") source('~/1_R/git/git-aviary/DRmencia/R/plot_ANOVA_pois_comm.R') anovaplot_pois_div <- plot_ANOVA_pois_comm(div.ANOVA.preds[i.aceit,], mod = m.anova.div.site, mod.dat.i = i.aceit.mod, ylab. = "", point.sz = 4, ymin. = 5, ymax = 20)
save.image("./plots/regrxnplot_pois_div.RData") save.image("./plots/anovaplot_pois_div.RData")
cowplot::plot_grid(regrxnplot_pois_div, anovaplot_pois_div, rel_widths = c(8,1))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.