########## Code for Nelson et al. in submission ######## ## Forb diversity globally is harmed by nutrient enrichment but can be rescued by large mammalian herbivory #### ######### by Rebecca Nelson ############# ########### created on 7/4/22 ########## ####### last updated: 7/4/24 ######## ###################################### ### All data files associated with code are available in the Nutrient Network Dropbox Data Repository ### rm(list = ls()) ######## load required packages ######## #library(tidyverse) #require(lme4) #require(nlme) #library(sjPlot) #for plotting lmer and glmer mods #library(sjmisc) #library(effects) #library(sjstats) #use for r2 functions #library(Rmisc) #library(ggpubr) #require(effects) #library(gridExtra) ######################### ### Read in data #### #comb <- read.csv("comb-by-plot-2022-05-02.csv", # na.strings=c("NULL", "NA"), # header=T, strip.white=T) #clim <- read.csv("site-worldclim-9-April-2021.csv", # na.strings=c("NULL", "NA"), # header=T, strip.white=T) #cru <- read.csv("CRU-annual_2018-07-06.csv", # na.strings=c("NULL", "NA"), # header=T, strip.white=T) #biomass <- read.csv("full-biomass-2022-05-02.csv", # na.strings=c("NULL", "NA"), # header=T, strip.white=T) #cover <- read.csv("full-cover-2022-05-02.csv", # na.strings=c("NULL", "NA"), # header=T, strip.white=T) #site <- read.csv("sites-2022-10-27.csv", # na.strings=c("NULL", "NA"), # header=T, strip.white=T) ### Organize data ##### ### biomass variable control <- biomass %>% filter(live == 1) %>% filter(year_trt == 1) %>% filter(trt == "Control") %>% filter(site_code %in% c("hero.uk", "kiny.au", "pinj.au", "rook.uk", "sereng.tz" , "spin.us" , "tyso.us" , "ahth.is" ,"amlr.is" , "bayr.de" , "bldr.us" , "kbs.us" , "pape.de" ,"thth.is" , "tmlr.is" , "ethamc.au" , "bogong.au" , "comp.pt" ,"glcr.us" , "tejon_foot.us" , "tejon_south.us" ,"burren.ie" ,"hnvr.us" ,"lancaster.uk" , "neba.jp" , "trel.us" ,"ethass.au" ,"kibber.in" ,"kilp.fi" , "mcla.us" , "sage.us" ,"sgs.us" ,"sier.us" ,"wupatki.us" , "yarra.au" , "bnbt.us" ,"burrawan.au" ,"cbgb.us" ,"elliot.us" ,"koffler.ca" , "msla.us" ,"saana.fi" ,"saline.us" ,"barta.us" ,"cdpt.us" , "cowi.ca" ,"gilb.za" ,"kidman.au" ,"look.us" ,"msla_2.us" , "msla_3.us" ,"nilla.au" ,"smith.us" ,"arch.us" ,"azi.cn" , "bari.ar" ,"frue.ch" ,"hall.us" , "hustai.mn" ,"mtca.au" , "ping.au" ,"valm.ch" ,"bnch.us" , "cdcr.us" ,"chilcas.ar", "hart.us" ,"hopl.us" ,"konz.us" , "potrok.ar" ,"sava.us" , "veluwe.nl" , "shps.us" ,"lubb.us" , "unc.us" ,"marc.ar" )) %>% group_by(plot, trt, site_code, block) %>% dplyr::summarise(Control_Live_Biomass = sum(mass)) fence <- biomass %>% filter(live == 1) %>% filter(year_trt == 1) %>% filter(trt == "Fence") %>% group_by(plot, trt, site_code, block) %>% dplyr::summarise(Fence_Live_Biomass = sum(mass)) live <- left_join(fence, control, by = c("site_code", "block")) intensity <- live %>% group_by(block, site_code) %>% dplyr::summarise(Difference = (Control_Live_Biomass - Fence_Live_Biomass)) %>% group_by(site_code) %>% dplyr::summarise(mean(Difference)) intensity$Intensity <- intensity$`mean(Difference)` #live <- biomass %>% filter(live == 1) %>% filter(year_trt == 1) %>% filter(trt %in% c( "Fence", "Control")) %>% group_by(plot, trt, site_code, block) %>% dplyr::summarise(Live_Biomass = sum(mass)) ### create columns to test for interactions ####### ### NPK x Fence exp cover <- cover %>% mutate( Fence = case_when( trt == "Fence" ~ 1, trt == "NPK+Fence" ~ 1, trt == "Control" ~ 0, trt == "NPK" ~ 0, )) cover$Fence <- as.factor(cover$Fence) cover <- cover %>% mutate( NPK = case_when( trt == "Fence" ~ 0, trt == "NPK+Fence" ~ 1, trt == "Control" ~ 0, trt == "NPK" ~ 1, )) cover$NPK <- as.factor(cover$NPK) ### Fert exp cover <- cover %>% mutate( N = case_when( trt == "Control" ~ 0, trt == "K" ~ 0, trt == "P" ~ 0, trt == "N" ~ 1, trt == "PK" ~ 0, trt == "NK" ~ 1, trt == "NP" ~ 1, trt == "NPK" ~ 1, )) cover$N <- as.factor(cover$N) cover <- cover %>% mutate( P = case_when( trt == "Control" ~ 0, trt == "K" ~ 0, trt == "P" ~ 1, trt == "N" ~ 0, trt == "PK" ~ 1, trt == "NK" ~ 0, trt == "NP" ~ 1, trt == "NPK" ~ 1, )) cover$P <- as.factor(cover$P) cover <- cover %>% mutate( K = case_when( trt == "Control" ~ 0, trt == "K" ~ 1, trt == "P" ~ 0, trt == "N" ~ 0, trt == "PK" ~ 1, trt == "NK" ~ 1, trt == "NP" ~ 0, trt == "NPK" ~ 1, )) cover$K <- as.factor(cover$K) ## remove oddball site cover <- cover %>% filter(site_code != "cereep.fr") obs.cover <- cover[cover$year_trt == 0,] cover <- cover[cover$year_trt > 0, ] #### Q1-2 Forbs ############# ## 1. organize variables ###### #cover includes all trt #using obs.cover calculates variables for non-pretreatment years ## calculate total cover for each plot in a given a year cover <- cover %>% group_by(plot, year, trt, site_code) %>% mutate(Total_Plot_Cover = sum(max_cover)) obs.cover <- obs.cover %>% group_by(plot, year, trt, site_code) %>% mutate(Total_Plot_Cover = sum(max_cover)) #Total_Plot_Cover = total cover for a plot ##calculate normalized cover for each species in a given sample cover <- cover %>% group_by(plot, year, trt, site_code, Taxon) %>% mutate(Normalized_Cover = max_cover/Total_Plot_Cover) obs.cover <- obs.cover %>% group_by(plot, year, trt, site_code, Taxon) %>% mutate(Normalized_Cover = max_cover/Total_Plot_Cover) #remove non plant cover: cover <- cover %>% filter(Taxon != "OTHER LITTER") %>% filter(Taxon != "GROUND") obs.cover <- obs.cover %>% filter(Taxon != "OTHER LITTER") %>% filter(Taxon != "GROUND") ### calculate proportion of sp richness from key families ##### ##using pre-trt data obs.cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(year_trt) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #1576 total species obs.cover %>% filter(Family %in% c("Compositae")) %>% group_by(year_trt) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #439 total species obs.cover %>% filter(Family %in% c("Fabaceae")) %>% group_by(year_trt) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #199 total species obs.cover %>% filter(Family %in% c("Geraniaceae")) %>% group_by(year_trt) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #25 total species obs.cover %>% filter(Family %in% c("Apiaceae")) %>% group_by(year_trt) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #59 total species obs.cover %>% filter(Family %in% c("Polygonaceae")) %>% group_by(year_trt) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #35 total species #list of all families present fam_list <- obs.cover %>% filter(functional_group %in% c("FORB", "LEGUME")) unique(fam_list$Family) fam_table <- obs.cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(year_trt, Family) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) write.csv(fam_table,"~fam_table.csv", row.names = FALSE) ## calculate response variables #### #Forb Species Richness sp_richness <- cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(block, plot, year, trt, year_trt, site_code, Fence, N, K, P, NPK) %>% dplyr::summarise(Forb_Sp_Richness = length(unique(Taxon))) obs.sp_richness <- obs.cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(block, plot, year, trt, year_trt, site_code, Fence, N, K, P, NPK) %>% dplyr::summarise(Forb_Sp_Richness.0 = length(unique(Taxon))) #Forb Family Richness fam_richness <- cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(block, plot, year, trt, year_trt, site_code, Fence, N, K, P, NPK) %>% filter(Family != "NA") %>% dplyr::summarise(Forb_Fam_Richness= length(unique(Family))) obs.fam_richness <- obs.cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(block, plot, year, trt, year_trt, site_code, Fence, N, K, P, NPK) %>% filter(Family != "NA") %>% dplyr::summarise(Forb_Fam_Richness.0 = length(unique(Family))) ##Normalized Forb Cover (total forb & legume cover per plot normalized) forb_cover <- cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, Fence, N, K, P, NPK) %>% dplyr::summarise(Forb_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Forb_Cover = Forb_Cover/Total_Plot_Cover) obs.forb_cover <- obs.cover %>% filter(functional_group %in% c("FORB", "LEGUME")) %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, Fence, NPK, N, K, P) %>% dplyr::summarise(Forb_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Forb_Cover.0 = Forb_Cover/Total_Plot_Cover) ### Grass Richness (total grass & gramminoid richness) grass_sp_richness <- cover %>% filter(functional_group %in% c("GRASS", "GRAMMINOID")) %>% group_by(block, plot, year, trt, year_trt, site_code, Fence, N, K, P, NPK) %>% dplyr::summarise(Grass_Richness = length(unique(Taxon))) obs.grass_sp_richness <- obs.cover %>% filter(functional_group %in% c("GRASS", "GRAMMINOID")) %>% group_by(block, plot, year, trt, year_trt, site_code, Fence, N, K, P, NPK) %>% dplyr::summarise(Grass_Richness.0 = length(unique(Taxon))) ##Normalized Grass Cover (total grass & gramminoid cover per plot normalized) grass_cover <- cover %>% filter(functional_group %in% c("GRASS", "GRAMMINOID")) %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, Fence, N, K, P, NPK) %>% dplyr::summarise(Grass_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Grass_Cover = Grass_Cover/Total_Plot_Cover) obs.grass_cover <- obs.cover %>% filter(functional_group %in% c("GRASS", "GRAMMINOID")) %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, Fence, N, K, P, NPK) %>% dplyr::summarise(Grass_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Grass_Cover.0 = Grass_Cover/Total_Plot_Cover) ## calculate LRRs ########## #Merge in pre-treatment for each data set sp.rich <- merge(sp_richness, obs.sp_richness[c("site_code", "plot", "Forb_Sp_Richness.0")], by=c("site_code", "plot"), all.x=T) fam.rich <- merge(fam_richness, obs.fam_richness[c("site_code", "plot", "Forb_Fam_Richness.0")], by=c("site_code", "plot"), all.x=T) forb.cover <- merge(forb_cover, obs.forb_cover[c("site_code", "plot", "Normalized_Forb_Cover.0")], by=c("site_code", "plot"), all.x=T) grass.rich <- merge(grass_sp_richness, obs.grass_sp_richness[c("site_code", "plot", "Grass_Richness.0")], by=c("site_code", "plot"), all.x=T) grass.cover <- merge(grass_cover, obs.grass_cover[c("site_code", "plot", "Normalized_Grass_Cover.0")], by=c("site_code", "plot"), all.x=T) #Calculate LRR effect using pre-treatment data in same plot sp.rich$lrr2 <- with(sp.rich, log(Forb_Sp_Richness/Forb_Sp_Richness.0)) fam.rich$lrr2 <- with(fam.rich, log(Forb_Fam_Richness/Forb_Fam_Richness.0)) forb.cover$lrr2 <- with(forb.cover, log(Normalized_Forb_Cover/Normalized_Forb_Cover.0)) grass.cover$lrr2 <- with(grass.cover, log(Normalized_Grass_Cover/Normalized_Grass_Cover.0)) grass.rich$lrr2 <- with(grass.rich, log(Grass_Richness/Grass_Richness.0)) ### 1. Analysis ########## ## LRRs show change from pre-trt to trt data for each trt_year #### Part 1 and 2 ######### ### Example: family richness ##### fam.rich <- left_join(fam.rich, intensity, by = "site_code") #fam.rich.PAR <-left_join(fam.rich, PAR_clean, by = c("block", "plot","year", "trt", "year_trt", "site_code")) fam_richness_nut <- fam.rich %>% filter(trt %in% c("Control", "K", "NK", "NPK", "P", "PK", "N", "NP")) fam_richness_fence <- fam.rich %>% filter(trt %in% c("Control", "Fence","NPK", "NPK+Fence")) mod_fam_rich.2 <- lme(lrr2 ~ N*P*K, data = fam_richness_nut, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) ### with climate added in ##CRU fam.rich <- left_join (fam.rich, cru, by=c("site_code", "year")) #note 2017 onward doesn't have climate data fam.rich$PET.s <- scale(fam.rich$PET, center = TRUE, scale = TRUE) mod_fam_rich.clim <- lme(lrr2 ~ NPK*Fence*(PET.s + Intensity), data = fam.rich, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) plot(fitted(mod_fam_rich.clim),resid(mod_fam_rich.clim)) qqnorm(resid(mod_fam_rich.clim)) summary(mod_fam_rich.clim) ## repeat for each response variable #### make figure and table 1 ####### plot_models(list(mod_sp_rich.clim, mod_grass_rich.clim), axis.labels = c("NPKμ x Fence x Intensity", "NPKμ x Fence x PET", "Fence x Intensity", "Fence x PET", "NPKμ x Intensity", "NPKμ x PET","NPKμ x Fence", "Intensity", "PET", "Fence", "NPKμ"), show.values=FALSE, show.p=TRUE, colors = c("forestgreen", "purple4"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.4, .4) plot_models(list(mod_forb_cover_clim, mod_grass_cover.clim), axis.labels = c("NPKμ x Fence x Intensity", "NPKμ x Fence x PET", "Fence x Intensity", "Fence x PET", "NPKμ x Intensity", "NPKμ x PET","NPKμ x Fence", "Intensity", "PET", "Fence", "NPKμ"), show.values=FALSE, show.p=TRUE, colors = c("forestgreen", "purple4"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.4, .4) tab_model(mod_sp_rich.clim, mod_fam_rich.clim, mod_forb_cover_clim, mod_grass_rich.clim, mod_grass_cover.clim, show.stat = TRUE, dv.labels = c("Forb Species Richness", "Forb Family Richness", "Forb Cover", "Grass Species Richness", "Grass Cover"), collapse.ci = TRUE, p.style = "numeric_stars" , file = "Table_1.doc") tab_model(mod_sp_rich.clim, mod_fam_rich.clim, mod_forb_cover_clim, mod_grass_rich.clim, mod_grass_cover.clim, show.stat = TRUE, dv.labels = c("Forb Species Richness", "Forb Family Richness", "Forb Cover", "Grass Species Richness", "Grass Cover"), collapse.ci = TRUE, p.style = "numeric_stars") ## interaction plots library(effects) eff.p1 <- effect("NPK*Fence*Intensity", mod_forb_cover_clim, KR=T) plot(eff.p1) #### figure 2 and table 2 ######### plot_models(list(mod_sp_rich.2, mod_grass_rich.2), axis.labels = c("N x P x Kμ", "P x Kμ", "N x Kμ", "N x P", "Kμ", "P", "N"), show.values=FALSE, show.p=TRUE, colors = c("forestgreen", "purple4"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.3, .3) plot_models(list(mod_forb_cover.2, mod_grass_cover.2), axis.labels = c("N x P x Kμ", "P x Kμ", "N x Kμ", "N x P", "Kμ", "P", "N"), show.values=FALSE, show.p=TRUE, colors = c("forestgreen", "purple4"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.3, .3) plot_models(list(mod_forb_cover.2.drop, mod_grass_cover.2.drop), axis.labels = c("Kμ", "P", "N"), show.values=FALSE, show.p=TRUE, colors = c("forestgreen", "purple4"), vline.color = "gray", dot.size = 5, line.size = 2 ) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.3, .3) tab_model(mod_sp_rich.2, mod_fam_rich.2, mod_forb_cover.2.drop, mod_grass_rich.2, mod_grass_cover.2.drop, show.stat = TRUE, dv.labels = c("Forb Species Richness", "Forb Family Richness", "Forb Cover", "Grass Species Richness", "Grass Cover"), collapse.ci = TRUE, p.style = "numeric_stars", file = "Table_2.doc") tab_model(mod_sp_rich.2, mod_fam_rich.2, mod_forb_cover.2.drop, mod_grass_rich.2, mod_grass_cover.2.drop, show.stat = TRUE, dv.labels = c("Forb Species Richness", "Forb Family Richness", "Forb Cover", "Grass Species Richness", "Grass Cover"), collapse.ci = TRUE, p.style = "numeric_stars") ## supplementary table of full interactions for cover tab_model(mod_forb_cover.2, mod_grass_cover.2, show.stat = TRUE, dv.labels = c("Forb Cover", "Grass Cover", "Grass Cover"), collapse.ci = TRUE, p.style = "numeric_stars", file = "Sup_Table_4.doc") ### Q3 family-level analysis: ######################################## ##### Compositae ################ ### Species Richness within family comp_rich <- cover %>% filter(Family == ("Compositae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness = length(unique(Taxon))) comp_rich.0 <- obs.cover %>% filter(Family == ("Compositae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness.0 = length(unique(Taxon))) # Total normalized cover at the family level comp_fam.0 <- obs.cover %>% filter(Family == "Compositae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover.0 = Fam_Cover/Total_Plot_Cover) comp_fam <- cover %>% filter(Family == "Compositae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover = Fam_Cover/Total_Plot_Cover) #### Fabaceae ################# ### Species Richness within family fab_rich <- cover %>% filter(Family == ("Fabaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarise(Sp_Richness = length(unique(Taxon))) fab_rich.0 <- obs.cover %>% filter(Family == ("Fabaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarise(Sp_Richness.0 = length(unique(Taxon))) # Total normalized cover at the family level fab_fam.0 <- obs.cover %>% filter(Family == "Fabaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover.0 = Fam_Cover/Total_Plot_Cover) # Total normalized cover at the family level fab_fam <- cover %>% filter(Family == "Fabaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover = Fam_Cover/Total_Plot_Cover) #### Geraniaceae ##### ### Species Richness within family ger_rich <- cover %>% filter(Family == ("Geraniaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness = length(unique(Taxon))) ger_rich.0 <- obs.cover %>% filter(Family == ("Geraniaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness.0 = length(unique(Taxon))) # Total normalized cover at the family level ger_fam.0 <- obs.cover %>% filter(Family == "Geraniaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover.0 = Fam_Cover/Total_Plot_Cover) # Total normalized cover at the family level ger_fam <- cover %>% filter(Family == "Geraniaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover = Fam_Cover/Total_Plot_Cover) #### Apiaceae ######### ### Species Richness within family api_rich <- cover %>% filter(Family == ("Apiaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness = length(unique(Taxon))) api_rich.0 <- obs.cover %>% filter(Family == ("Apiaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness.0 = length(unique(Taxon))) # Total normalized cover at the family level api_fam.0 <- obs.cover %>% filter(Family == "Apiaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover.0 = Fam_Cover/Total_Plot_Cover) # Total normalized cover at the family level api_fam <- cover %>% filter(Family == "Apiaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover = Fam_Cover/Total_Plot_Cover) ##### Polygonaceae ############# ### Species Richness within family pol_rich <- cover %>% filter(Family == ("Polygonaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness = length(unique(Taxon))) pol_rich.0 <- obs.cover %>% filter(Family == ("Polygonaceae")) %>% group_by(block, plot, year, trt, year_trt, site_code, NPK, Fence, N, P, K) %>% dplyr::summarize(Sp_Richness.0 = length(unique(Taxon))) # Total normalized cover at the family level pol_fam.0 <- obs.cover %>% filter(Family == "Polygonaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover.0 = Fam_Cover/Total_Plot_Cover) # Total normalized cover at the family level pol_fam <- cover %>% filter(Family == "Polygonaceae") %>% group_by(block, plot, year, trt, year_trt, site_code, Total_Plot_Cover, NPK, Fence, N, P, K) %>% dplyr::summarise(Fam_Cover = sum(max_cover)) %>% dplyr::mutate(Normalized_Fam_Cover = Fam_Cover/Total_Plot_Cover) ## calculate LRRs ########## #Merge in pre-treatment for each data set comp.rich <- merge(comp_rich, comp_rich.0[c("site_code", "plot", "Sp_Richness.0")], by=c("site_code", "plot"), all.x=T) comp.fam <- merge(comp_fam, comp_fam.0[c("site_code", "plot", "Normalized_Fam_Cover.0")], by=c("site_code", "plot"), all.x=T) fab.rich <- merge(fab_rich, fab_rich.0[c("site_code", "plot", "Sp_Richness.0")], by=c("site_code", "plot"), all.x=T) fab.fam <- merge(fab_fam, fab_fam.0[c("site_code", "plot", "Normalized_Fam_Cover.0")], by=c("site_code", "plot"), all.x=T) ger.rich <- merge(ger_rich, ger_rich.0[c("site_code", "plot", "Sp_Richness.0")], by=c("site_code", "plot"), all.x=T) ger.fam <- merge(ger_fam, ger_fam.0[c("site_code", "plot", "Normalized_Fam_Cover.0")], by=c("site_code", "plot"), all.x=T) api.rich <- merge(api_rich, api_rich.0[c("site_code", "plot", "Sp_Richness.0")], by=c("site_code", "plot"), all.x=T) api.fam <- merge(api_fam, api_fam.0[c("site_code", "plot", "Normalized_Fam_Cover.0")], by=c("site_code", "plot"), all.x=T) pol.rich <- merge(pol_rich, pol_rich.0[c("site_code", "plot", "Sp_Richness.0")], by=c("site_code", "plot"), all.x=T) pol.fam <- merge(pol_fam, pol_fam.0[c("site_code", "plot", "Normalized_Fam_Cover.0")], by=c("site_code", "plot"), all.x=T) #Calculate LRR effect using pre-treatment data in same plot comp.rich$lrr2 <- with(comp.rich, log(Sp_Richness/Sp_Richness.0)) comp.fam$lrr2 <- with(comp.fam, log(Normalized_Fam_Cover/Normalized_Fam_Cover.0)) fab.rich$lrr2 <- with(fab.rich, log(Sp_Richness/Sp_Richness.0)) fab.fam$lrr2 <- with(fab.fam, log(Normalized_Fam_Cover/Normalized_Fam_Cover.0)) ger.rich$lrr2 <- with(ger.rich, log(Sp_Richness/Sp_Richness.0)) ger.fam$lrr2 <- with(ger.fam, log(Normalized_Fam_Cover/Normalized_Fam_Cover.0)) api.rich$lrr2 <- with(api.rich, log(Sp_Richness/Sp_Richness.0)) api.fam$lrr2 <- with(api.fam, log(Normalized_Fam_Cover/Normalized_Fam_Cover.0)) pol.rich$lrr2 <- with(pol.rich, log(Sp_Richness/Sp_Richness.0)) pol.fam$lrr2 <- with(pol.fam, log(Normalized_Fam_Cover/Normalized_Fam_Cover.0)) ##### 2. updated stats with LRR ########## #### composites ########### ## richness comp.rich <- left_join(comp.rich, intensity, by="site_code") comp_rich_nut <- comp.rich %>% filter(trt %in% c("Control", "K", "NK", "NPK", "P", "PK", "N", "NP")) comp_rich_fence <- comp.rich %>% filter(trt %in% c("Control", "Fence","NPK", "NPK+Fence")) mod_comp_rich <- lme(lrr2 ~ N*P*K, data = comp_rich_nut, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) mod_comp_rich_fence <- lme(lrr2 ~ NPK*Fence*Intensity, data = comp_rich_fence, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) ## family level cover comp.fam <- left_join(comp.fam, intensity, by="site_code") comp_fam_nut <- comp.fam %>% filter(trt %in% c("Control", "K", "NK", "NPK", "P", "PK", "N", "NP")) comp_fam_fence <- comp.fam %>% filter(trt %in% c("Control", "Fence","NPK", "NPK+Fence")) mod_comp_cover <- lme(lrr2 ~ N+P+K, data = comp_fam_nut, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) mod_comp_cover.2 <- lme(lrr2 ~ N*P*K, data = comp_fam_nut, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) mod_comp_cover_fence <- lme(lrr2 ~ NPK*Fence*Intensity, data = comp_fam_fence, random = ~1|site_code/block, correlation = corAR1(form = ~year_trt|site_code/block/plot), na.action=na.exclude) # repeat for other families ### Figure 3 and Table 3-4 ######### #Figure 3 plot_models(list(mod_comp_rich_fence, mod_fab_rich_fence), axis.labels = c("NPKμ x Fence x Intensity", "Fence x Intensity", "NPKμ x Intensity", "NPKμ x Fence", "Intensity", "Fence", "NPKμ"), show.values=FALSE, show.p=TRUE, colors = c( "mediumslateblue", "goldenrod2"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.7, .7) plot_models(list(mod_comp_cover_fence, mod_fab_cover_fence), axis.labels = c("NPKμ x Fence x Intensity", "Fence x Intensity", "NPKμ x Intensity", "NPKμ x Fence", "Intensity", "Fence", "NPKμ"), show.values=FALSE, show.p=TRUE, colors = c( "mediumslateblue", "goldenrod2"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.7, .7) plot_models(mod_comp_rich_fence, mod_fab_rich_fence, mod_ger_rich_fence, mod_api_rich_fence, mod_pol_rich_fence, axis.labels = c("NPKμ x Fence x Intensity", "Fence x Intensity", "NPKμ x Intensity", "NPKμ x Fence", "Intensity", "Fence", "NPKμ"), show.values=FALSE, show.p=TRUE, colors = c("sienna", "magenta4", "hotpink1", "mediumslateblue", "goldenrod2"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.7, .7) plot_models(mod_comp_cover_fence, mod_fab_cover_fence, mod_ger_cover_fence, mod_api_cover_fence, mod_pol_cover_fence, axis.labels = c("NPKμ x Fence x Intensity", "Fence x Intensity", "NPKμ x Intensity", "NPKμ x Fence", "Intensity", "Fence", "NPKμ"), show.values=FALSE, show.p=TRUE, colors = c("sienna", "magenta4", "hotpink1", "mediumslateblue", "goldenrod2"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.7, .7) plot_models(mod_comp_rich, mod_fab_rich, mod_ger_rich, mod_api_rich, mod_pol_rich, axis.labels = c("N x P x Kμ", "P x Kμ", "N x Kμ", "N x P", "Kμ", "P", "N"), show.values=FALSE, show.p=TRUE, colors = c("sienna", "magenta4", "hotpink1", "mediumslateblue", "goldenrod2"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.7, .7) plot_models(mod_comp_cover, mod_fab_cover, mod_ger_cover, mod_api_cover, mod_pol_cover, axis.labels = c( "Kμ", "P", "N"), show.values=FALSE, show.p=TRUE, colors = c("sienna", "magenta4","hotpink1", "mediumslateblue", "goldenrod2"), vline.color = "gray", dot.size = 5, line.size = 2) + set_theme( base = theme_classic(), axis.title.size = 1.8, axis.textsize = 1.8) + ylim(-.7, .7) tab_model(mod_comp_rich_fence, mod_comp_cover_fence, mod_fab_rich_fence, mod_fab_cover_fence,mod_ger_rich_fence, mod_ger_cover_fence, mod_api_rich_fence, mod_api_cover_fence, mod_pol_rich_fence, mod_pol_cover_fence, show.stat = TRUE, dv.labels = c("Asteraceae Richness", "Asteraceae Cover", "Fabaceae Richness", "Fabaceae Cover", "Geraniaceae Richness", "Geraniaceae Cover", "Apiaceae Richness", "Apiaceae Cover", "Polygonaceae Richness", "Polygonaceae Cover"), collapse.ci = TRUE, p.style = "numeric_stars", file = "Table_3.doc") tab_model(mod_comp_rich, mod_comp_cover, mod_fab_rich, mod_fab_cover,mod_ger_rich, mod_ger_cover, mod_api_rich, mod_api_cover, mod_pol_rich, mod_pol_cover, show.stat = TRUE, dv.labels = c("Asteraceae Richness", "Asteraceae Cover", "Fabaceae Richness", "Fabaceae Cover", "Geraniaceae Richness", "Geraniaceae Cover", "Apiaceae Richness", "Apiaceae Cover", "Polygonaceae Richness", "Polygonaceae Cover"), collapse.ci = TRUE, p.style = "numeric_stars", file = "Table_4.doc") tab_model(mod_comp_rich, mod_comp_cover, mod_fab_rich, mod_fab_cover,mod_ger_rich, mod_ger_cover, mod_api_rich, mod_api_cover, mod_pol_rich, mod_pol_cover, show.stat = TRUE, dv.labels = c("Asteraceae Richness", "Asteraceae Cover", "Fabaceae Richness", "Fabaceae Cover", "Geraniaceae Richness", "Geraniaceae Cover", "Apiaceae Richness", "Apiaceae Cover", "Polygonaceae Richness", "Polygonaceae Cover"), collapse.ci = TRUE, p.style = "numeric_stars") tab_model(mod_comp_rich_fence, mod_comp_cover_fence, mod_fab_rich_fence, mod_fab_cover_fence,mod_ger_rich_fence, mod_ger_cover_fence, mod_api_rich_fence, mod_api_cover_fence, mod_pol_rich_fence, mod_pol_cover_fence, show.stat = TRUE, dv.labels = c("Asteraceae Richness", "Asteraceae Cover", "Fabaceae Richness", "Fabaceae Cover", "Geraniaceae Richness", "Geraniaceae Cover", "Apiaceae Richness", "Apiaceae Cover", "Polygonaceae Richness", "Polygonaceae Cover"), collapse.ci = TRUE, p.style = "numeric_stars") tab_model(mod_comp_cover.2, mod_fab_cover.2, mod.ger.cover.2, mod_api_cover.2, mod_pol_cover.2, show.stat = TRUE, dv.labels = c("Asteraceae Cover", "Fabaceae Cover", "Geraniaceae Cover", "Apiaceae Cover", "Polygonaceae Cover"), collapse.ci = TRUE, p.style = "numeric_stars", file = "Sup_Table_5.doc") library(effects) eff.p1 <- effect("NPK*Fence*Intensity", mod_comp_cover_fence, KR=T) plot(eff.p1) eff.p <- effect("NPK*Fence*Intensity", mod_comp_rich_fence, KR=T) plot(eff.p) eff.p <- effect("NPK*Fence*Intensity", mod_pol_rich_fence, KR=T) plot(eff.p)