diff --git a/.Rproj.user/EAEB7255/console06/98DD95A6 b/.Rproj.user/EAEB7255/console06/98DD95A6 new file mode 100644 index 0000000..994544e --- /dev/null +++ b/.Rproj.user/EAEB7255/console06/98DD95A6 @@ -0,0 +1,4 @@ +Microsoft Windows [Version 6.1.7601][?25l +Copyright (c) 2009 Microsoft Corporation. All rights reserved. + +C:\Users\Sam\Documents\Research\MS Thesis\Understory>[?25h \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/console06/INDEX001 b/.Rproj.user/EAEB7255/console06/INDEX001 index 0637a08..7b88d0e 100644 --- a/.Rproj.user/EAEB7255/console06/INDEX001 +++ b/.Rproj.user/EAEB7255/console06/INDEX001 @@ -1 +1 @@ -[] \ No newline at end of file +[{"allow_restart":true,"alt_buffer":false,"autoclose":1,"buffered_output":"","caption":"Terminal 1","channel_id":"7526","channel_mode":1,"child_procs":false,"cols":93,"cwd":"","exit_code":15,"handle":"98DD95A6","interaction_mode":2,"max_output_lines":1000,"restarted":false,"rows":48,"shell_type":3,"show_on_output":false,"terminal_sequence":1,"title":"","track_env":false,"zombie":false}] \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/pcs/source-pane.pper b/.Rproj.user/EAEB7255/pcs/source-pane.pper index 1743e40..9efb51b 100644 --- a/.Rproj.user/EAEB7255/pcs/source-pane.pper +++ b/.Rproj.user/EAEB7255/pcs/source-pane.pper @@ -1,3 +1,3 @@ { - "activeTab" : 0 + "activeTab" : -1 } \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/pcs/windowlayoutstate.pper b/.Rproj.user/EAEB7255/pcs/windowlayoutstate.pper index 3745081..48df712 100644 --- a/.Rproj.user/EAEB7255/pcs/windowlayoutstate.pper +++ b/.Rproj.user/EAEB7255/pcs/windowlayoutstate.pper @@ -2,7 +2,7 @@ "left" : { "panelheight" : 779, "splitterpos" : 326, - "topwindowstate" : "NORMAL", + "topwindowstate" : "HIDE", "windowheight" : 817 }, "right" : { diff --git a/.Rproj.user/EAEB7255/sources/per/t/679490EB b/.Rproj.user/EAEB7255/sources/per/t/679490EB deleted file mode 100644 index 8ce3d6b..0000000 --- a/.Rproj.user/EAEB7255/sources/per/t/679490EB +++ /dev/null @@ -1,29 +0,0 @@ -{ - "collab_server" : "", - "contents" : "", - "created" : 1536263599868.000, - "dirty" : false, - "encoding" : "", - "folds" : "", - "hash" : "0", - "id" : "679490EB", - "lastKnownWriteTime" : 7308604914264011873, - "last_content_update" : 1536263599868, - "path" : null, - "project_path" : null, - "properties" : { - "cacheKey" : "29B3C11D", - "caption" : "ms", - "contentUrl" : "grid_resource/gridviewer.html?env=&obj=ms&cache_key=29B3C11D", - "displayedObservations" : "1020", - "environment" : "", - "object" : "ms", - "preview" : "0", - "totalObservations" : "1020", - "variables" : "13" - }, - "relative_order" : 5, - "source_on_save" : false, - "source_window" : "", - "type" : "r_dataframe" -} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/per/t/6DD0CFE5 b/.Rproj.user/EAEB7255/sources/per/t/6DD0CFE5 deleted file mode 100644 index 9fb9528..0000000 --- a/.Rproj.user/EAEB7255/sources/per/t/6DD0CFE5 +++ /dev/null @@ -1,29 +0,0 @@ -{ - "collab_server" : "", - "contents" : "", - "created" : 1547069427357.000, - "dirty" : false, - "encoding" : "", - "folds" : "", - "hash" : "0", - "id" : "6DD0CFE5", - "lastKnownWriteTime" : 1546898750, - "last_content_update" : 1547069427357, - "path" : null, - "project_path" : null, - "properties" : { - "cacheKey" : "6F672F8B", - "caption" : "plot_daub_cover", - "contentUrl" : "grid_resource/gridviewer.html?env=&obj=plot_daub_cover&cache_key=6F672F8B", - "displayedObservations" : 102, - "environment" : "", - "object" : "plot_daub_cover", - "preview" : 0, - "totalObservations" : 102, - "variables" : 13 - }, - "relative_order" : 9, - "source_on_save" : false, - "source_window" : "", - "type" : "r_dataframe" -} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/per/t/7EBDE282 b/.Rproj.user/EAEB7255/sources/per/t/7EBDE282 deleted file mode 100644 index 8eeb248..0000000 --- a/.Rproj.user/EAEB7255/sources/per/t/7EBDE282 +++ /dev/null @@ -1,29 +0,0 @@ -{ - "collab_server" : "", - "contents" : "", - "created" : 1536262763531.000, - "dirty" : false, - "encoding" : "", - "folds" : "", - "hash" : "0", - "id" : "7EBDE282", - "lastKnownWriteTime" : 7308604914264011873, - "last_content_update" : 1536262763531, - "path" : null, - "project_path" : null, - "properties" : { - "cacheKey" : "71129F7A", - "caption" : "spp_cov", - "contentUrl" : "grid_resource/gridviewer.html?env=&obj=spp_cov&cache_key=71129F7A", - "displayedObservations" : 2002, - "environment" : "", - "object" : "spp_cov", - "preview" : 0, - "totalObservations" : 2002, - "variables" : 5 - }, - "relative_order" : 4, - "source_on_save" : false, - "source_window" : "", - "type" : "r_dataframe" -} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/per/t/EE7BA85C b/.Rproj.user/EAEB7255/sources/per/t/EE7BA85C deleted file mode 100644 index 127c787..0000000 --- a/.Rproj.user/EAEB7255/sources/per/t/EE7BA85C +++ /dev/null @@ -1,23 +0,0 @@ -{ - "collab_server" : "", - "contents" : "", - "created" : 1546972070942.000, - "dirty" : false, - "encoding" : "UTF-8", - "folds" : "", - "hash" : "495975393", - "id" : "EE7BA85C", - "lastKnownWriteTime" : 1547068489, - "last_content_update" : 1547068489536, - "path" : "~/Research/MS Thesis/Understory/understory_electivity.R", - "project_path" : "understory_electivity.R", - "properties" : { - "cursorPosition" : "6,18", - "scrollLine" : "0", - "tempName" : "Untitled1" - }, - "relative_order" : 8, - "source_on_save" : false, - "source_window" : "", - "type" : "r_source" -} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/25B4B1DC b/.Rproj.user/EAEB7255/sources/prop/25B4B1DC new file mode 100644 index 0000000..652161e --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/prop/25B4B1DC @@ -0,0 +1,4 @@ +{ + "cursorPosition" : "36,29", + "scrollLine" : "27" +} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/2A13C8B0 b/.Rproj.user/EAEB7255/sources/prop/2A13C8B0 index 241bcf3..3dff34d 100644 --- a/.Rproj.user/EAEB7255/sources/prop/2A13C8B0 +++ b/.Rproj.user/EAEB7255/sources/prop/2A13C8B0 @@ -1,4 +1,4 @@ { - "cursorPosition" : "180,52", - "scrollLine" : "162" + "cursorPosition" : "298,47", + "scrollLine" : "286" } \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/2C2B9CE b/.Rproj.user/EAEB7255/sources/prop/2C2B9CE new file mode 100644 index 0000000..7a73a41 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/prop/2C2B9CE @@ -0,0 +1,2 @@ +{ +} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/5E4621F b/.Rproj.user/EAEB7255/sources/prop/5E4621F new file mode 100644 index 0000000..652161e --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/prop/5E4621F @@ -0,0 +1,4 @@ +{ + "cursorPosition" : "36,29", + "scrollLine" : "27" +} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC b/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC new file mode 100644 index 0000000..1c98d2e --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC @@ -0,0 +1,4 @@ +{ + "cursorPosition" : "13,21", + "scrollLine" : "4" +} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6 b/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6 index fb1ea21..1529b1f 100644 --- a/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6 +++ b/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6 @@ -1,4 +1,4 @@ { - "cursorPosition" : "94,50", - "scrollLine" : "20" + "cursorPosition" : "23,28", + "scrollLine" : "17" } \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/D455C530 b/.Rproj.user/EAEB7255/sources/prop/D455C530 index 245dc77..b709900 100644 --- a/.Rproj.user/EAEB7255/sources/prop/D455C530 +++ b/.Rproj.user/EAEB7255/sources/prop/D455C530 @@ -1,5 +1,5 @@ { - "cursorPosition" : "6,18", - "scrollLine" : "0", + "cursorPosition" : "359,2", + "scrollLine" : "348", "tempName" : "Untitled1" } \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/prop/INDEX b/.Rproj.user/EAEB7255/sources/prop/INDEX index 9399f56..a632a06 100644 --- a/.Rproj.user/EAEB7255/sources/prop/INDEX +++ b/.Rproj.user/EAEB7255/sources/prop/INDEX @@ -15,7 +15,10 @@ C%3A%2FUsers%2FSam%2FGoogle%20Drive%2FThesis%20related%2FUnderstory%2Funderstory C%3A%2FUsers%2FSam%2FGoogle%20Drive%2FThesis%20related%2FUnderstory%2Funderstory_effects_plots.R="5FD4B92F" C%3A%2FUsers%2FSam%2FGoogle%20Drive%2FThesis%20related%2FUnderstory%2Funderstory_electivity.R="16C87132" ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Fplot_level_understory.R="2A13C8B0" +~%2FResearch%2FMS%20Thesis%2FUnderstory%2Fplot_level_understory_082418.R="25B4B1DC" ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Fplot_level_understory_090618.R="B04D5E6E" +~%2FResearch%2FMS%20Thesis%2FUnderstory%2Fplot_level_understory_112817.R="5E4621F" ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_data_prep.R="A1FFFDE6" +~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_effects_plots.R="8DF0DCEC" ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_electivity.R="D455C530" ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_electivity010819.R="8CD3B36B" diff --git a/.Rproj.user/EAEB7255/sources/per/t/679490EB-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/17AED2FF-contents similarity index 100% rename from .Rproj.user/EAEB7255/sources/per/t/679490EB-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/17AED2FF-contents diff --git a/.Rproj.user/EAEB7255/sources/per/t/EE7BA85C-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/488CD1EA-contents similarity index 99% rename from .Rproj.user/EAEB7255/sources/per/t/EE7BA85C-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/488CD1EA-contents index 5fc890e..d70bd7e 100644 --- a/.Rproj.user/EAEB7255/sources/per/t/EE7BA85C-contents +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/488CD1EA-contents @@ -143,10 +143,10 @@ plots_to_use <- unique(ms$Plot)[plots_to_use_dead & plots_to_use_live] #------------------------------------------------------------------------------ #perform analysis for the four FTs types <- c("All", "Perennial grass", "Cheatgrass", "Perennial forb", "Shrub") - +ntypes <- length(types) #initialize a list of dataframes to catch the electivities at each plot/FT elect_results <- list() -for(i in 1:5){ +for(i in 1:ntypes){ elect_results[[i]] <- data.frame(Plot = numeric(0), Dead = numeric(0), Live = numeric(0), @@ -206,7 +206,7 @@ nit <- 999 #number of monte carlo draws #initialize a list of dataframes to catch the electivities at each plot/FT elect_means <- list() -for(i in 1:5){ +for(i in 1:ntypes){ elect_means[[i]] <- data.frame(Iter = numeric(0), Dead = numeric(0), Live = numeric(0), @@ -279,7 +279,7 @@ for(type in types){ results_boots <- list() -for(i in 1:5){ +for(i in 1:ntypes){ results_boots[[i]] <- data.frame(Dead = numeric(2), Live = numeric(2), Inter = numeric(2), @@ -291,7 +291,7 @@ for(i in 1:5){ names(results_boots) <- types #calculate p-values from distribution of bootstrapped mean electivities -for(i in 1:5){ +for(i in 1:ntypes){ emp_mean <- apply(elect_results[[i]][, 2:7], 2, FUN = function(x){mean(x, na.rm = TRUE)}) pvals <- c(sum(elect_means[[i]][2] > emp_mean[1]), @@ -327,7 +327,7 @@ par(mfrow = c(3,2), mar = c(2,1,1,2), oma = c(2,3,1,0)) -for(i in 1:5){ +for(i in 1:ntypes){ melt_elect <- melt(elect_results[[i]][, 2:4]) plot(melt_elect$value ~ I(as.numeric(melt_elect$variable)+ runif(nrow(melt_elect), -.1, .1 )), diff --git a/.Rproj.user/EAEB7255/sources/per/t/6DD0CFE5-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/49D1D8FF-contents similarity index 100% rename from .Rproj.user/EAEB7255/sources/per/t/6DD0CFE5-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/49D1D8FF-contents diff --git a/.Rproj.user/EAEB7255/sources/per/t/7EBDE282-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/71D46E86-contents similarity index 100% rename from .Rproj.user/EAEB7255/sources/per/t/7EBDE282-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/71D46E86-contents diff --git a/.Rproj.user/EAEB7255/sources/per/t/F6290AA8-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/9CCF2D15-contents similarity index 100% rename from .Rproj.user/EAEB7255/sources/per/t/F6290AA8-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/9CCF2D15-contents diff --git a/.Rproj.user/EAEB7255/sources/per/t/BBCED885 b/.Rproj.user/EAEB7255/sources/s-A74A5E15/9E5C163F similarity index 56% rename from .Rproj.user/EAEB7255/sources/per/t/BBCED885 rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/9E5C163F index 8dd655d..8551ebe 100644 --- a/.Rproj.user/EAEB7255/sources/per/t/BBCED885 +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/9E5C163F @@ -1,21 +1,21 @@ { "collab_server" : "", "contents" : "", - "created" : 1547060059118.000, + "created" : 1547225685816.000, "dirty" : false, "encoding" : "UTF-8", "folds" : "", - "hash" : "482206355", - "id" : "BBCED885", - "lastKnownWriteTime" : 1547087115, - "last_content_update" : 1547087115203, + "hash" : "1518548566", + "id" : "9E5C163F", + "lastKnownWriteTime" : 1551912261, + "last_content_update" : 1551912261877, "path" : "~/Research/MS Thesis/Understory/plot_level_understory.R", "project_path" : "plot_level_understory.R", "properties" : { - "cursorPosition" : "180,52", - "scrollLine" : "162" + "cursorPosition" : "298,47", + "scrollLine" : "286" }, - "relative_order" : 8, + "relative_order" : 1, "source_on_save" : false, "source_window" : "", "type" : "r_source" diff --git a/.Rproj.user/EAEB7255/sources/per/t/BBCED885-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/9E5C163F-contents similarity index 86% rename from .Rproj.user/EAEB7255/sources/per/t/BBCED885-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/9E5C163F-contents index f691c4c..1fc9220 100644 --- a/.Rproj.user/EAEB7255/sources/per/t/BBCED885-contents +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/9E5C163F-contents @@ -11,8 +11,8 @@ library("plyr") setwd("C:/Users/Sam/Documents/Research/MS Thesis/Understory") plot_data <- read.csv("./clean data/plot_data.csv") +greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv") -# process other files and import #------------------------------------------------------------ # Analysis @@ -22,7 +22,7 @@ plot_data <- read.csv("./clean data/plot_data.csv") #exploration hist(asin(sqrt(plot_data$All/100))) hist(logit(plot_data$All/100 + min(plot_data[plot_data$All != 0, "All"]))) -hist(log(plot_data$All/100 + .01)) +hist(log(plot_data$All/100)) hist(asin(sqrt(plot_data$Cheatgrass/100))) hist(logit(plot_data$Cheatgrass/100 + min(plot_data[plot_data$Cheatgrass != 0, "Cheatgrass"]))) @@ -48,13 +48,32 @@ plot(plot_data$Cheatgrass ~ plot_data$cwd_normal_cum) #--------------------------------------------------------------- #linear mixed effects models for functional groups at plot level #------------------------------------------------------------------- + +## all + +all_plot <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +all_plot <- lmer(logit(All/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +summary(lm(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC), data = plot_data)) + +summary(all_plot, ddf = "Kenward-Roger") +r.squaredGLMM(all_plot) +plot(allEffects(all_plot, partial.residuals = TRUE)) +plot(all_plot) +AICc(all_plot) +scatter.smooth(residuals(all_plot) ~ predict(all_plot)) + ## Cheatgrass cheatgrass_plot <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC) + (1|Cluster), data = plot_data) -cheatgrass_plot <- lmer(logit(Cheatgrass/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + - scale(AWC) + (1|Cluster), data = plot_data) +# cheatgrass_plot <- lmer(logit(Cheatgrass/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + # scale(AWC) + (1|Cluster), data = plot_data) summary(lm(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC), data = plot_data)) @@ -86,35 +105,40 @@ plot(pgrass_plot) pforb_plot <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC) + (1|Cluster), data = plot_data) -summary(lm(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + - scale(AWC), data = plot_data)) - summary(pforb_plot, ddf = "Kenward-Roger") r.squaredGLMM(pforb_plot) plot(allEffects(pforb_plot, partial.residuals = TRUE)) +plot(inv.as(predict(pforb_plot)) ~ I(plot_data$Pforb/100)) +abline(0,1) plot(pforb_plot) ## Shrubs shrub_plot <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC) + (1|Cluster), data = plot_data) - -summary(lm(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + - scale(AWC), data = plot_data)) +# shrub_plot <- lmer(log(Shrub_cover_li+.01) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + +# scale(AWC) + (1|Cluster), data = plot_data) summary(shrub_plot, ddf = "Kenward-Roger") r.squaredGLMM(shrub_plot) plot(allEffects(shrub_plot, partial.residuals = TRUE)) +plot(inv.as(predict(shrub_plot)) ~ I(plot_data$Shrub_cover_li)) +abline(0,1) plot(shrub_plot) +### Plot regression coefficients +library(coefplot2) +coefplot2(pgrass_plot) +coefplot2(cheatgrass_plot, add = TRUE, spacing = 10) + + + #----------------------------------------------- # Comparison of greenwood and 2015 #----------------------------------------------- -compare <- join(plot_daub_cover, greenwood_under, by = "Plot", type = "inner") -compare <- join(compare, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot", type = "inner") -compare <- join(compare, plot_data[, c("Plot", "Dead_ba", "Delta_pdc", "cwd_normal_cum")], by = "Plot") +compare <- join(plot_data, greenwood_under, by = "Plot", type = "inner") #cheatgrass plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover) @@ -159,7 +183,7 @@ mean(compare$Pgrass-compare$Perrenial.Grass.Cover) t.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") wilcox.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") -pg_change_lm <- lm(Pgrass - Perrenial.Grass.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +pg_change_lm <- lm(Pgrass - Perrenial.Grass.Cover ~ Delta_pdc*cwd_normal_cum, data= compare) summary(pg_change_lm) ## Perr forb @@ -172,7 +196,7 @@ mean(I(compare$Aforb + compare$Pforb) - compare$Forb.Cover) t.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") wilcox.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") -pf_change_lm <- lm(Aforb + Pforb - Forb.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +pf_change_lm <- lm(Aforb + Pforb - Forb.Cover ~ Delta_pdc*cwd_normal_cum, data= compare) summary(pf_change_lm) ## Shrub @@ -187,7 +211,7 @@ t.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="les wilcox.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="less") hist(compare$Shrub - I(compare$Shrub.Cover.Total)) -s_change_lm <- lm(Shrub - Shrub.Cover.Total ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +s_change_lm <- lm(Shrub - Shrub.Cover.Total ~ Delta_pdc*cwd_normal_cum, data= compare) summary(s_change_lm) ###################### @@ -199,7 +223,7 @@ png(filename="change_in_cover.png", type="cairo", units="in", width = 4, - height=6, + height=4, pointsize=12, res=160) @@ -223,6 +247,22 @@ plot(compare$Pgrass ~ compare$Perrenial.Grass.Cover, xlim = c(0,35), ylim = c(0, axis(side = 2, at = c(0,5,15,25,35), tcl = -.24, labels = FALSE) abline(0,1, lty = 2, col = 'grey30') +par(mar = c(2,2,2,1)) +plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover, xlim = c(0, 20), ylim = c(0, 20), + xlab = "", ylab = "", main = "Cheatgrass", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 1, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 2, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + + par(mar = c(2,2,2,1)) plot(compare$Pforb ~ compare$Forb.Cover, xlim = c(0, 12), ylim = c(0, 12), main = "Perennial Forb", xlab = "", ylab = "", xaxt = 'n', yaxt = 'n', @@ -255,21 +295,7 @@ plot(compare$Shrub ~ compare$Shrub.Cover.Total, xlim = c(0, 30), ylim = c(0,30), abline(0,1, lty = 2, col = 'grey30') -par(mar = c(2,2,2,1)) -plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover, xlim = c(0, 20), ylim = c(0, 20), - xlab = "", ylab = "", main = "Cheatgrass", xaxt = 'n', yaxt = 'n', - bg='grey60', - col = 'grey30', - pch = 21, - cex = 1.5, - cex.axis = .8, - cex.main = 1.3) - axis(side = 1, at = c(0,5,10,15,20), labels = TRUE) - axis(side = 1, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) - axis(side = 2, at = c(0,5,10,15,20), labels = TRUE) - axis(side = 2, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) - abline(0,1, lty = 2, col = 'grey30') - + mtext("2005 Cover (%)", side = 1, outer = TRUE) mtext("2015 Cover (%)", side = 2, outer = TRUE) dev.off() diff --git a/.Rproj.user/EAEB7255/sources/per/t/F6290AA8 b/.Rproj.user/EAEB7255/sources/s-A74A5E15/A8857C5B similarity index 67% rename from .Rproj.user/EAEB7255/sources/per/t/F6290AA8 rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/A8857C5B index be4d563..ff63e6e 100644 --- a/.Rproj.user/EAEB7255/sources/per/t/F6290AA8 +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/A8857C5B @@ -1,28 +1,28 @@ { "collab_server" : "", "contents" : "", - "created" : 1536264339225.000, + "created" : 1551911591824.000, "dirty" : false, "encoding" : "", "folds" : "", "hash" : "0", - "id" : "F6290AA8", - "lastKnownWriteTime" : 7308604914264011873, - "last_content_update" : 1536264339225, + "id" : "A8857C5B", + "lastKnownWriteTime" : 1547476559, + "last_content_update" : 1551911591824, "path" : null, "project_path" : null, "properties" : { - "cacheKey" : "EDD72D8A", + "cacheKey" : "14D3D425", "caption" : "plot_data", - "contentUrl" : "grid_resource/gridviewer.html?env=&obj=plot_data&cache_key=EDD72D8A", + "contentUrl" : "grid_resource/gridviewer.html?env=&obj=plot_data&cache_key=14D3D425", "displayedObservations" : 102, "environment" : "", "object" : "plot_data", "preview" : 0, "totalObservations" : 102, - "variables" : 23 + "variables" : 25 }, - "relative_order" : 7, + "relative_order" : 4, "source_on_save" : false, "source_window" : "", "type" : "r_dataframe" diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/A8857C5B-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/A8857C5B-contents new file mode 100644 index 0000000..e69de29 diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/AD22B2CF-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/AD22B2CF-contents new file mode 100644 index 0000000..d7aa228 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/AD22B2CF-contents @@ -0,0 +1,439 @@ +## Understory electivity +library(plyr) +library(vioplot) +library(reshape2) +library(multcompView) + +set.seed(16091315) + +daub <- read.csv("./raw data/daub_cover.csv", stringsAsFactors = FALSE) + +#some data proofing + daub <- daub[(daub$Transect %in% c("N", "E", "S", "W")), ] + daub$Plot <- as.character(daub$Plot) + daub[daub$Plot == "NPElectricEel", "Plot"] <- "NPELECTRICEEL" + daub[daub$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120" + daub[daub$Plot == "NPElectricEel240", "Plot"] <- "NPELECTRICEEL240" + daub[daub$Plot == "NPElectricEel360", "Plot"] <- "NPELECTRICEEL360" +daub$Midpoint.value <- as.numeric(as.character((daub$Midpoint.value))) +daub$unique_quad <- paste0(daub$Plot, daub$Transect, daub$Meter) +daub[daub$Cover.type == "Perennial forb ", "Cover.type"] <- "Perennial forb" +daub[daub$Cover.type == "Shrub ", "Cover.type"] <- "Shrub" + +species<- read.csv("./raw data/spp_cover2.csv") +species <- species[(species$Transect %in% c("N", "E", "S", "W")), ] +species$unique_quad <- paste0(species$Plot, species$Transect, species$Meter) + +count(species, vars = c("Plot", "Transect")) + +#----------------------------------------------------------------------- +# Calculate mean cover, cover by quadrat, cover by plot +#------------------------------------------------------------------- +# Mean overall cover +length(unique(daub$unique_quad)) +spp_abund <- aggregate(species$Cover / 2040, by = list(species$Species), FUN = sum) +spp_abund <- spp_abund[order(spp_abund[, 1], decreasing = FALSE), ] + + + +#percent of quadrats occupied +quadrat_abund <- as.data.frame(table(species[species$Cover != 0, ]$Species)) +quadrat_abund <- quadrat_abund[order(quadrat_abund[, 1], decreasing = FALSE), ] +quadrat_abund[, 2] <- quadrat_abund[, 2] * 100 / 2040 + + +count <- table(species$Species, species$Plot) +is.not.zero <- function(x){ + if((x) == 0){return(0)} + else{return(1)} +} + +presence <- apply(count, c(1,2), FUN = is.not.zero) +plot_presence <- as.data.frame(rowSums(presence))/102*100 + +spp_abund <- cbind(spp_abund, quadrat_abund[-1, 2], plot_presence[-1, 1]) + +write.csv(spp_abund, "./outputs/species_summary.csv") + + +## generate cover data for different subsets +n_quads_occ <- length(unique(species$unique_quad)) #number of quadrat-by-group records + +spp_cov <- data.frame(unique_quad = unique(species$unique_quad), + poasec = numeric(n_quads_occ), + othergrass = numeric(n_quads_occ), + phlhoo = numeric(n_quads_occ), + otherforb = numeric(n_quads_occ)) + + +for(i in 1:n_quads_occ){ + spp_cov$poasec[i] <- ifelse(species[species$unique_quad == spp_cov$unique_quad[i], ]$Species == "POASEC", + species[species$unique_quad == spp_cov$unique_quad[i] & species$Species == "POASEC", ]$Cover, 0) + + spp_cov$othergrass[i] <- ifelse(species[species$unique_quad == spp_cov$unique_quad[i], ]$Species %in% c("PSESPI", "LEYCIN", "ELYELY", "STITHU"), + sum(species[species$unique_quad == spp_cov$unique_quad[i] + & species$Species %in% c("PSESPI", "LEYCIN", "ELYELY", "STITHU"), ]$Cover), 0) + + spp_cov$phlhoo[i] <- ifelse(species[species$unique_quad == spp_cov$unique_quad[i], ]$Species == "PHLHOO", + species[species$unique_quad == spp_cov$unique_quad[i] & species$Species == "PHLHOO", ]$Cover, 0) + + spp_cov$otherforb[i] <- ifelse(species[species$unique_quad == spp_cov$unique_quad[i], ]$Species %in% c("CREACA", "CREACU", "STEACA", "ERICAE", "LEPPUN", "ERIUMB", "AREACU","LESKIN","GALAPA","ASTPUR","BALSAG","CRYFLA","ERIMIC","MINKIN","ASTOOP","OPUPOL","ARAHOL","ANTDIM","BOEHOL"), + sum(species[species$unique_quad == spp_cov$unique_quad[i] + & species$Species %in% c("CREACA", "CREACU", "STEACA", "ERICAE", "LEPPUN", "ERIUMB", "AREACU","LESKIN","GALAPA","ASTPUR","BALSAG","CRYFLA","ERIMIC","MINKIN","ASTOOP","OPUPOL","ARAHOL","ANTDIM","BOEHOL"), ]$Cover), 0) +} + +#---------------------------------------------------------------------------- + +# import and recode microsite data +ms <- read.csv("./raw data/microsite.csv") +ms$Microsite <- toupper(ms$Microsite) +ms$ms <- ifelse(ms$Microsite %in% c("PI", "JI", "CI", "PO", "JO", "CO"), "Live", + # ifelse(ms$Microsite %in% c("PO", "JO", "CO"), "Live Outer", + ifelse(ms$Microsite %in% c("PI(S)", "PO(S)", "JI(S)", "JO(S)", "CI(S)", "CO(S)", "LOG"), "Dead", + # ifelse(ms$Microsite =="LOG", "Log", + ifelse(ms$Microsite == "I", "Inter", NA)))#)) +ms[ms$Transect == "e", "Transect"] <- "E" +ms[ms$Transect == "s", "Transect"] <- "S" +ms[ms$Transect == "w", "Transect"] <- "W" +ms[ms$Transect == "n", "Transect"] <- "N" + +ms$Plot <- as.character(ms$Plot) + +ms <- ms[(ms$Transect %in% c("N", "S", "E", "W")), ] + +ms[ms$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120" + + +ms$unique_quad <- paste0(ms$Plot, ms$Transect, ms$Meter) + +#get total cover for each plot +total_cover <- aggregate(daub$Midpoint.value, by = list(daub$Plot, daub$Cover.type), FUN = sum) + +## select which plots have enough microsites of each type +# this is a real friggin mess but it pulls out the plots which have +# at least n_plots dead and at least n_plots live ms +plots_to_use_dead <- NA +plots_to_use_live <- NA +n_plots <- 1 + +for(i in 1:length(unique(ms$Plot))){ + table <- table(ms[ms$Plot == unique(ms$Plot)[i], ]$ms) + + if("Dead" %in% names(table)){ + plots_to_use_dead[i] <- table["Dead"] + } else{ + plots_to_use_dead[i] <- 0 + } + + if("Live" %in% names(table)){ + plots_to_use_live[i] <- table["Live"] + } else{ + plots_to_use_live[i] <- 0 + } +} + +plots_to_use_dead <- (plots_to_use_dead >= n_plots) +plots_to_use_live <- (plots_to_use_live >= n_plots) + +plots_to_use <- unique(ms$Plot)[plots_to_use_dead & plots_to_use_live] + + +#------------------------------------------------------------------------------ +# Calculate test and randomized distributions +#------------------------------------------------------------------------------ +#perform analysis for the four FTs +types <- c("Perennial grass", "Cheatgrass", "Perennial forb", "Shrub") +ntypes <- length(types) +#initialize a list of dataframes to catch the electivities at each plot/FT +elect_results <- list() +for(i in 1:ntypes){ + elect_results[[i]] <- data.frame(Plot = numeric(0), + Dead = numeric(0), + Live = numeric(0), + Inter = numeric(0), + DL = numeric(0), + DI = numeric(0), + LI = numeric(0)) +} +names(elect_results) <- types + + +for(type in types){ + + for (i in 1:length(plots_to_use)){ + + if(type == "All"){ + daub_cov <- daub[daub$Plot == as.character(plots_to_use[i]), ] + daub_cov <- daub_cov[daub_cov$Cover.type %in% c("Cheatgrass", "Other Ann Grass", + "Perennial grass", "Annual forb", "Perennial forb", "Shrub"), ] + daub_cov_2 <- aggregate(daub_cov$Midpoint.value, by = list(daub_cov$unique_quad), FUN = sum) + names(daub_cov_2) <- c("unique_quad", "Midpoint.value") + daub_cov <- join(daub_cov_2, ms[, c("unique_quad", "ms")], by = "unique_quad", type = "inner") + }else{ + daub_cov <- daub[daub$Cover.type == type & daub$Plot == as.character(plots_to_use[i]), ] + daub_cov <- join(daub_cov, ms[, c("unique_quad", "ms")], by = "unique_quad", type = "inner") + + } + + daub_cov$prop_cov <- daub_cov$Midpoint.value / sum(daub_cov$Midpoint.value) + + cov <- aggregate(daub_cov[, "prop_cov"], by = list(daub_cov$ms), FUN = sum) + + prev <- table(daub_cov$ms)/20 + + #this is a mess because tables are hard to use. There's definitely a better way to do this + elect <- data.frame(Plot = plots_to_use[i], + Dead = (cov[cov$Group.1 == "Dead",2] - prev[which(names(prev) == "Dead")]) / + (cov[cov$Group.1 == "Dead", 2] + prev[which(names(prev) == "Dead")]), + Live = (cov[cov$Group.1 == "Live",2] - prev[which(names(prev) == "Live")]) / + (cov[cov$Group.1 == "Live",2] + prev[which(names(prev) == "Live")]), + Inter = (cov[cov$Group.1 == "Inter",2] - prev[which(names(prev) == "Inter")]) / + (cov[cov$Group.1 == "Inter",2] + prev[which(names(prev) == "Inter")])) + elect$DL = elect$Dead - elect$Live + elect$DI = elect$Dead - elect$Inter + elect$LI = elect$Live - elect$Inter + + elect_results[[type]] <- rbind(elect_results[[type]], elect) + } +} + +saveRDS(elect_results, paste0("./outputs/results_elect_", n_plots, "dead.rds")) +# elect_results <- readRDS("./outputs/results_elect_1dead.rds") +#------------------------------------------------------------------------------ +#Calculate monte carlo means +#------------------------------------------------------------------------------ +nit <- 999 #number of monte carlo draws + +#initialize a list of dataframes to catch the electivities at each plot/FT +elect_means <- list() +for(i in 1:ntypes){ + elect_means[[i]] <- data.frame(Iter = numeric(0), + Dead = numeric(0), + Live = numeric(0), + Inter = numeric(0), + DL = numeric(0), + DI = numeric(0), + LI = numeric(0)) +} +names(elect_means) <- types + +system.time(#takes about a half hour + +for(type in types){ + + plot_to_use_type <- subset(elect_results[[type]], !is.na(Dead))$Plot + + for(j in 1:nit){ + + elect_rand <- data.frame(Plot = character(0), + Dead = numeric(0), + Live = numeric(0), + Inter = numeric(0), + DL = numeric(0), + DI = numeric(0), + DL = numeric(0)) + + for (plot_use in plot_to_use_type){ + if(type == "All"){ + daub_cov <- daub[daub$Plot == as.character(plot_use), ] + daub_cov <- daub_cov[daub_cov$Cover.type %in% c("Cheatgrass", "Other Ann Grass", + "Perennial grass", "Annual forb", "Perennial forb", "Shrub"), ] + daub_cov_2 <- aggregate(daub_cov$Midpoint.value, by = list(daub_cov$unique_quad), FUN = sum) + names(daub_cov_2) <- c("unique_quad", "Midpoint.value") + daub_cov <- join(daub_cov_2, ms[, c("unique_quad", "ms")], by = "unique_quad", type = "inner") + }else{ + daub_cov <- daub[daub$Cover.type == type & daub$Plot == as.character(plot_use), ] + daub_cov <- join(daub_cov, ms[, c("unique_quad", "ms")], by = "unique_quad", type = "inner") + + } + + daub_cov$prop_cov <- daub_cov$Midpoint.value / sum(daub_cov$Midpoint.value) + daub_cov$prop_cov_ran <- sample(daub_cov$prop_cov, size = 20, replace = FALSE) + + cov <- aggregate(daub_cov[, "prop_cov_ran"], by = list(daub_cov$ms), FUN = sum) + prev <- table(daub_cov$ms)/20 + + elect <- data.frame(Plot = plot_use, + Dead = (cov[cov$Group.1 == "Dead",2] - prev[which(names(prev) == "Dead")]) / + (cov[cov$Group.1 == "Dead", 2] + prev[which(names(prev) == "Dead")]), + Live = (cov[cov$Group.1 == "Live",2] - prev[which(names(prev) == "Live")]) / + (cov[cov$Group.1 == "Live",2] + prev[which(names(prev) == "Live")]), + Inter = (cov[cov$Group.1 == "Inter",2] - prev[which(names(prev) == "Inter")]) / + (cov[cov$Group.1 == "Inter",2] + prev[which(names(prev) == "Inter")])) + elect$DL = elect$Dead - elect$Live + elect$DI = elect$Dead - elect$Inter + elect$LI = elect$Live - elect$Inter + + elect_rand <- rbind(elect_rand, elect) + } + + means <- apply(elect_rand[, c(2:7)], 2, FUN = function(x){mean(x, na.rm = TRUE)}) + elect_means[[type]][j, 1] <- j + elect_means[[type]][j, 2:7] <- means + + + } + +} +) + + +results_boots <- list() +for(i in 1:ntypes){ + results_boots[[i]] <- data.frame(Dead = numeric(2), + Live = numeric(2), + Inter = numeric(2), + DL = numeric(2), + DI = numeric(2), + LI = numeric(2)) + rownames(results_boots[[i]]) <- c("Empirical", "N_boots_greater") +} +names(results_boots) <- types + +#calculate p-values from distribution of bootstrapped mean electivities +for(i in 1:ntypes){ +emp_mean <- apply(elect_results[[i]][, 2:7], 2, FUN = function(x){mean(x, na.rm = TRUE)}) + +pvals <- c(sum(elect_means[[i]][2] > emp_mean[1]), + sum(elect_means[[i]][3] > emp_mean[2]), + sum(elect_means[[i]][4] > emp_mean[3]), + sum(elect_means[[i]][5] > emp_mean[4]), + sum(elect_means[[i]][6] > emp_mean[5]), + sum(elect_means[[i]][7] > emp_mean[6])) + +#flip around large p-values to make them small p-values for the left tail +pvals <- ifelse(pvals > (nit + 1)/2, (nit + 1 - pvals)/(nit + 1), pvals/(nit + 1)) + +results_boots[[i]][1, ] <- emp_mean +results_boots[[i]][2, ] <- pvals +} + +saveRDS(results_boots, paste0("./outputs/results_boots_", n_plots, "dead.rds")) + +# results_boots <- readRDS("./outputs/results_boots_1dead.rds") +#----------------------------------------------------------------------------------------------- +# Plot of electivity for each FT +#----------------------------------------------------------------------------------------------- + +png(filename="electivity.png", + type="cairo", + units="in", + width = 4, + height=4, + pointsize=10, + res=160) + +par(mfrow = c(2,2), + mar = c(2,1,1,2), + oma = c(2,3,1,0)) + +for(i in 1:ntypes){ + + melt_elect <- melt(elect_results[[i]][, 2:4]) + plot(melt_elect$value ~ I(as.numeric(melt_elect$variable)+ runif(nrow(melt_elect), -.1, .1 )), + ylim = c(-1,1), + xlim = c(.7, 3.3), + pch = 21, + bg = "grey", + xaxt = "n") + pvals <- results_boots[[i]][2, ] + abline(h = 0) + + means <- aggregate(melt_elect$value, by = list(melt_elect$variable), FUN = function(x){mean(x, na.rm = TRUE)}) + + segments(x0 = c(0.9, 1.9, 2.9), y0 = means$x, x1 = c(1.1, 2.1, 3.1), lwd = 3) + + + # segments(x0 = 1, x1 = 2, y0 = melt_elect[melt_elect$variable == "Dead", "value"], + # y1 = melt_elect[melt_elect$variable == "Live", "value"]) + # + # segments(x0 = 2, x1 = 3, y0 = melt_elect[melt_elect$variable == "Live", "value"], + # y1 = melt_elect[melt_elect$variable == "Inter", "value"]) + + if(i %in% c(1,2)){ + axis(1, at = c(1,2,3), labels = FALSE) + } + + if(i %in% c(3,4)){ + axis(1, at = c(1,2,3), labels = c("Dead", "Live", "Interspace")) + } + + text(x = c(1,2,3), y = 0.9, labels = pvals[2:4]) + mtext(text = types[i], outer = FALSE, side = 3, line = 0.3) + mtext(text = paste0("(", letters[i], ")"), outer = FALSE, side = 3, at = 0.5, line = 0.3) + mtext(text = "Electivity", outer = TRUE, side = 2, line = 1.5) + + +} + +dev.off() +# +# #----------------------------------------------------------------------------------------------- +# # For whole study area +# #----------------------------------------------------------------------------------------------- +# all_merged <- join(daub, ms, by = "unique_quad", type = "inner") +# +# daub$unique_quad +# ms$unique_quad +# +# +# all_plots_prev <- table(all_merged$ms) +# all_plots_prev <- all_plots_prev / sum(all_plots_prev) +# +# all_plots_cover <- aggregate(all_merged$Midpoint.value, by = list(all_merged$Cover.type, all_merged$ms), FUN = sum) +# all_plots_cover <- all_plots_cover[all_plots_cover$Group.1 %in% c("Perennial grass", "Perennial forb ", "Shrub ", "Cheatgrass"), ] +# +# for (i in 1:4){ # standardized cover so it adds to 1 for each functional type +# types <- c("Cheatgrass", "Perennial forb ", "Perennial grass", "Shrub ") +# all_plots_cover[all_plots_cover$Group.1 == types[i], "x"] <- all_plots_cover[all_plots_cover$Group.1 == types[i], "x"] / +# sum(all_plots_cover[all_plots_cover$Group.1 == types[i], "x"]) +# } +# +# +# +# all_plots_cover <- dcast(all_plots_cover, Group.1 ~ Group.2) #reshapes the table +# +# all_elect <- data.frame(type = c("Cheatgrass", "Perennial forb", "Perennial grass", "Shrub"), +# Dead = numeric(4), +# Inter = numeric(4), +# LiveInner = numeric(4), +# LiveOuter = numeric(4), +# Log = numeric(4)) +# +# for (i in 2:6){ #calculate electivity +# for (j in 1:4){ +# all_elect[j, i] <- (all_plots_cover[j,i] - all_plots_prev[i - 1]) / (all_plots_cover[j,i] + all_plots_prev[i - 1]) +# } +# } +# +# +# #-------------------------------------------------------------------- +# +# spp_merged <- join(spp_cov, ms, by = "unique_quad", type = "inner") +# head(spp_merged) +# +# +# spp_plots_prev <- table(spp_merged$ms) +# spp_plots_prev <- spp_plots_prev / sum(spp_plots_prev) +# +# spp_merged$poasec <- spp_merged$poasec / sum(spp_merged$poasec) +# spp_merged$othergrass <- spp_merged$othergrass / sum(spp_merged$othergrass) +# spp_merged$phlhoo <- spp_merged$phlhoo / sum(spp_merged$phlhoo) +# spp_merged$otherforb <- spp_merged$otherforb / sum(spp_merged$otherforb) +# +# spp_ag <- aggregate(spp_merged[, c("poasec", "othergrass", "phlhoo", "otherforb")], by = list(spp_merged$ms), FUN = sum) +# +# +# spp_elect <- data.frame(type = c("POASEC", "Other Perr Grass", "PHLHOO", "Other Forb"), +# Dead = numeric(4), +# Inter = numeric(4), +# LiveInner = numeric(4), +# LiveOuter = numeric(4), +# Log = numeric(4)) +# +# for (i in 1:5){ #calculate electivity +# for (j in 1:4){ +# spp_elect[j, i+1] <- (unname(spp_ag[i,j + 1]) - unname(all_plots_prev[i])) / (unname(spp_ag[i,j + 1]) + unname(all_plots_prev[i])) +# } +# } +# diff --git a/.Rproj.user/EAEB7255/sources/per/t/8F5BFB6A-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/B0ABC46C-contents similarity index 100% rename from .Rproj.user/EAEB7255/sources/per/t/8F5BFB6A-contents rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/B0ABC46C-contents diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/C58850AC-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/C58850AC-contents new file mode 100644 index 0000000..1630074 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/C58850AC-contents @@ -0,0 +1,95 @@ +##Understory data formatting +## +library("reshape2") +library("plyr") +library("effects") + + +daub <- read.csv("./raw data/daub_cover.csv", stringsAsFactors = FALSE) + +#some data proofing +daub <- daub[(daub$Transect %in% c("N", "E", "S", "W")), ] +daub$Plot <- as.character(daub$Plot) +daub[daub$Plot == "NPElectricEel", "Plot"] <- "NPELECTRICEEL" +daub[daub$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120" +daub[daub$Plot == "NPElectricEel240", "Plot"] <- "NPELECTRICEEL240" +daub[daub$Plot == "NPElectricEel360", "Plot"] <- "NPELECTRICEEL360" +daub$Midpoint.value <- as.numeric(as.character((daub$Midpoint.value))) +daub$unique_quad <- paste0(daub$Plot, daub$Transect, daub$Meter) +daub[daub$Cover.type == "Perennial forb ", "Cover.type"] <- "Perennial forb" +daub[daub$Cover.type == "Shrub ", "Cover.type"] <- "Shrub" + +trees <- read.csv("./raw data/trees_updated_with_logs_041716.csv") +tree_pdc <- read.csv("./raw data/all_trees_with_delta_and_ENN_041916.csv") +greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv") +tcover <- read.csv("./raw data/tree_and_shrub_cover_020815.csv") + names(tcover)[3] <- "Shrub_cover_li" +source("./calculate_awc.R") +soil <- calculate_awc(soil_raw = "./raw data/soils_missing_imputed_012016.csv") +clim <- read.csv("./raw data/ALL_climate_variables.csv") + names(clim)[2] <- "Plot" +other_vars <- read.csv("./raw data/all_vars_EXPORT.csv") + + +#remaking data into plot average cover +plot_daub_cover <- dcast(daub, Plot ~ Cover.type, value.var = "Midpoint.value", fun.aggregate = mean) +names(plot_daub_cover) <- c("Plot", "Aforb", "Bg", "Cheatgrass", "Crust", "Gravel", "Litter", + "Agrass", "Pforb", "Pgrass", "Rock", "Shrub") +plot_daub_cover$All <- plot_daub_cover$Aforb + plot_daub_cover$Cheatgrass + plot_daub_cover$Agrass + plot_daub_cover$Pforb + + plot_daub_cover$Pgrass + plot_daub_cover$Shrub + +#------------------------------------------------------------- +#Create tree data plot-level + +#impute meanas for missing values +for (i in which(sapply(trees, is.numeric))) { + for (j in which(is.na(trees[, i]))) { + trees[j, i] <- mean(trees[trees[, "Plot"] == trees[j, "Plot"], i], na.rm = TRUE) + } +} + +plot_trees <- data.frame(Plot = character(102), + Live_ba = numeric(102), + Dead_ba = numeric(102), + Total_ba = numeric(102), + Dead_rat = numeric(102), + Delta_pdc = numeric(102)) + +plot_trees$Plot <- as.character(unique(plot_daub_cover$Plot)) + +for(i in (1:nrow(plot_trees))){ + plot_trees$Live_ba[i] <- sum(trees[trees$Live == "Y" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Dead_ba[i] <- sum(trees[trees$Live == "N" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Total_ba[i] <- sum(trees[trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, na.rm = TRUE) + plot_trees$Dead_rat[i] <- plot_trees$Dead_ba[i] / plot_trees$Total_ba[i] + plot_trees$Delta_pdc[i] <- weighted.mean(tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "Delta_pdc"], + tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "BA_cm"]) +} + + +#----------------------------------------------- +# Plot-level cover for everything and also predictor variables + +#fill in any NAs + +plot_data <- join(plot_daub_cover, plot_trees, by = c("Plot")) +plot_data <- join(plot_data, clim[, c("Plot", "cwd_normal_cum")], by = "Plot") +plot_data <- join(plot_data, soil[, c("Plot", "AWC")], by = "Plot") +plot_data <- join(plot_data, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot") +plot_data <- join(plot_data, other_vars[, c("Plot", "Cluster", "Avg_depth")], by = "Plot") + + +for (i in which(sapply(plot_data, is.numeric))) { + for (j in which(is.na(plot_data[, i]))) { + plot_data[j, i] <- mean(plot_data[, i], na.rm = TRUE) + } +} + +levels(plot_data$Cluster) <- c(levels(plot_data$Cluster), "NPELECTRICEEL") +plot_data[(plot_data$Plot %in% c("NPELECTRICEEL", "NPELECTRICEEL120", + "NPELECTRICEEL240", "NPELECTRICEEL360")), "Cluster"] <- "NPELECTRICEEL" + + +write.csv(plot_data, "./clean data/plot_data.csv") \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/C9AC37D6-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/C9AC37D6-contents new file mode 100644 index 0000000..321e94b --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/C9AC37D6-contents @@ -0,0 +1,404 @@ +library("MuMIn") +library("lme4") +library("car") +library("glmmADMB") +library("effects") +library("plyr") +library("lmerTest") +library("vegan") +library("reshape2") +library("plyr") + +setwd("C:/Users/Sam/Google Drive/Projects/MS Thesis/Understory") + +#-------------------------------------------------------- +##Understory data formatting +#-------------------------------------------------------- + + +daub <- read.csv("./Raw data/daub_cover.csv") + daub$Plot <- as.character(daub$Plot) + daub[daub$Plot == "SS2556120", "Plot"] <- "SS256120" + daub$Midpoint.value <- as.numeric(as.character((daub$Midpoint.value))) +trees <- read.csv("./Raw data/trees_updated_with_logs_041716.csv") +tree_pdc <- read.csv("./Raw data/all_trees_with_delta_and_ENN_041916.csv") +greenwood_under <- read.csv("./Raw data/Greenwood_Understory_Variables_SF_edits.csv") +tcover <- read.csv("./Raw data/tree_and_shrub_cover_020815.csv") +names(tcover)[3] <- "Shrub_cover_li" +source("calculate_awc.R") +soil <- calculate_awc("./Raw data/soils_missing_imputed_012016.csv") +clim <- read.csv("./Raw data/ALL_climate_variables.csv") +names(clim)[2] <- "Plot" +other_vars <- read.csv("./Raw data/all_vars_EXPORT.csv") +map_mat <- read.csv("./Raw data/map_mat_normals.csv") +names(map_mat)[2] <- "Plot" + + +#remaking data into plot average cover +plot_daub_cover <- dcast(daub, Plot ~ Cover.type, value.var = "Midpoint.value", fun.aggregate = mean) +names(plot_daub_cover) <- c("Plot", "Aforb", "Bg", "Cheatgrass", "Crust", "Gravel", "Litter", + "Agrass", "Pforb", "Pgrass", "Rock", "Shrub") + + +#------------------------------------------------------------- +#Create tree data plot-level + +#impute meanas for missing values +for (i in which(sapply(trees, is.numeric))) { + for (j in which(is.na(trees[, i]))) { + trees[j, i] <- mean(trees[trees[, "Plot"] == trees[j, "Plot"], i], na.rm = TRUE) + } +} + +plot_trees <- data.frame(Plot = character(102), + Live_ba = numeric(102), + Dead_ba = numeric(102), + Total_ba = numeric(102), + Dead_rat = numeric(102), + Delta_pdc = numeric(102)) + +plot_trees$Plot <- as.character(unique(plot_daub_cover$Plot)) + +for(i in (1:nrow(plot_trees))){ + plot_trees$Live_ba[i] <- sum(trees[trees$Live == "Y" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Dead_ba[i] <- sum(trees[trees$Live == "N" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Total_ba[i] <- sum(trees[trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, na.rm = TRUE) + plot_trees$Dead_rat[i] <- plot_trees$Dead_ba[i] / plot_trees$Total_ba[i] + plot_trees$Delta_pdc[i] <- weighted.mean(tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "Delta_pdc"], + tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "BA_cm"]) +} + +#----------------------------------------------- +# Plot-level cover for everything and also predictor variables + +#fill in any NAs + +plot_data <- join(plot_daub_cover, plot_trees, by = c("Plot")) +plot_data <- join(plot_data, clim[, c("Plot", "cwd_normal_cum", "avg_ann_ppt")], by = "Plot") +plot_data <- join(plot_data, soil[, c("Plot", "AWC")], by = "Plot") +plot_data <- join(plot_data, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot") +plot_data <- join(plot_data, other_vars[, c("Plot", "Cluster", "Avg_depth", "Elev", "Slope", "Swness")], by = "Plot") +plot_data <- join(plot_data, map_mat[, c("Plot", "tmean", "ppt")], by = "Plot") + +for (i in which(sapply(plot_data, is.numeric))) { + for (j in which(is.na(plot_data[, i]))) { + plot_data[j, i] <- mean(plot_data[, i], na.rm = TRUE) + } +} + +#------------------------------------------------------------ +# Analysis +#------------------------------------------------------------ +# +# plot_data <- plot_data[!(plot_data$Plot %in% c("NPELECTRICEEL", "NPELECTRICEEL120", +# "NPELECTRICEEL240", "NPELECTRICEEL360")), ] +# plot_data$Cluster <- droplevels(plot_data$Cluster) + +#exploration +hist(asin(sqrt(plot_data$Cheatgrass/100))) +hist(log(plot_data$Cheatgrass/100 + .01)) + +nrow(plot_data[plot_data$Cheatgrass != 0, ]) + +hist(asin(sqrt(plot_data$Pgrass/100))) +nrow(plot_data[plot_data$Pgrass != 0, ]) + +hist(asin(sqrt(plot_data$Pforb/100))) +nrow(plot_data[plot_data$Pforb != 0, ]) + +hist(asin(sqrt(plot_data$Shrub_cover_li))) +nrow(plot_data[plot_data$Shrub_cover_li != 0, ]) + +plot(plot_data$Cheatgrass ~ plot_data$cwd_normal_cum) + + +#--------------------------------------------------------------- +#linear mixed effects models for functional groups at plot level +#------------------------------------------------------------------- + +## Cheatgrass + +cheatgrass_plot <- lmer(log(Cheatgrass/100 + .01) ~ scale(Tree_cover) + scale(Dead_ba)*scale(cwd_normal_cum) + + scale(AWC) + scale(Avg_depth) + scale(Dead_ba)*scale(cwd_normal_cum) + (1|Cluster), data = plot_data, + na.action = na.omit) +summary(cheatgrass_plot, ddf = "Kenward-Roger") +r.squaredGLMM(cheatgrass_plot) +plot(allEffects(cheatgrass_plot, partial.residuals = TRUE)) +plot(cheatgrass_plot) +AICc(cheatgrass_plot) +scatter.smooth(residuals(cheatgrass_plot) ~ predict(cheatgrass_plot)) + +cgrass_set <- dredge(cheatgrass_plot) +head(cgrass_set) + +cg_model <- lmer(log(Cheatgrass/100 + .01) ~ scale(Delta_pdc)*scale(cwd_normal_cum) + (1|Cluster), data = plot_data, + na.action = na.fail) +summary(cg_model, ddf = "Kenward-Roger") +r.squaredGLMM(cg_model) + +## Perr grass + +pgrass_plot <- lmer(log(Pgrass/100+.01) ~ scale(Tree_cover) + scale(Dead_ba)*scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + scale(Avg_depth) + scale(Dead_ba)*scale(cwd_normal_cum) + scale(Dead_rat) + (1|Cluster), data = plot_data, + na.action = na.omit) +summary(pgrass_plot, ddf = "Kenward-Roger") +r.squaredGLMM(pgrass_plot) +plot(allEffects(pgrass_plot, partial.residuals = TRUE)) +plot(pgrass_plot) + +pgrass_set <- dredge(pgrass_plot) +head(pgrass_set) + +pg_model <- lmer(log(Pgrass/100) ~ (1|Cluster), data = plot_data, na.action = na.fail) +summary(pg_model) +r.squaredGLMM(pg_model) + +## Perr forbs + +pforb_plot <- lmer(log(Pforb/100+.01) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + scale(Avg_depth) + scale(Dead_ba)*scale(cwd_normal_cum) + (1|Cluster), data = plot_data, + na.action = na.fail) +summary(pforb_plot, ddf = "Kenward-Roger") +r.squaredGLMM(pforb_plot) +plot(allEffects(pforb_plot, partial.residuals = TRUE)) +plot(pforb_plot) + +pforb_set <- dredge(pforb_plot) +head(pforb_set) + +pf_model <- lmer(log(Pforb/100 +.01) ~ (1|Cluster), data = plot_data, na.action = na.fail) +summary(pf_model) +r.squaredGLMM(pf_model) + +## Shrubs +shrub_plot <- lmer(log(Shrub_cover_li + .01) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + scale(Avg_depth) + scale(Dead_ba)*scale(cwd_normal_cum) + (1|Cluster), data = plot_data, + na.action = na.fail) +summary(shrub_plot, ddf = "Kenward-Roger") +r.squaredGLMM(shrub_plot) +plot(allEffects(shrub_plot, partial.residuals = TRUE)) +plot(shrub_plot) + +shrub_set <- dredge(shrub_plot) +head(shrub_set) + +shrub_model<- lmer(log(Shrub_cover_li + .01) ~ scale(Tree_cover) + (1|Cluster), data = plot_data, + na.action = na.fail) +summary(shrub_model) +r.squaredGLMM(shrub_model) +plot(allEffects(shrub_model, partial.residuals = TRUE)) + +#----------------------------------------------- +# Comparison of greenwood and 2015 +#----------------------------------------------- +compare <- join(plot_daub_cover, greenwood_under, by = "Plot", type = "inner") +compare <- join(compare, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot", type = "inner") +compare <- join(compare, plot_data[, c("Plot", "Dead_ba", "Delta_pdc", "cwd_normal_cum")], by = "Plot") + +#cheatgrass +plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover) +summary(lm(compare$Cheatgrass ~ compare$Cheatgrass.Cover)) +cor.test(compare$Cheatgrass, compare$Cheatgrass.Cover, + method = c("spearman"), exact = FALSE) +boxplot(compare$Cheatgrass, compare$Cheatgrass.Cover) +mean(compare$Cheatgrass-compare$Cheatgrass.Cover) +wilcox.test(compare$Cheatgrass, compare$Cheatgrass.Cover, paired = TRUE, alternative="greater") +t.test(compare$Cheatgrass, compare$Cheatgrass.Cover, paired = TRUE, alternative="greater") + +cor.test(compare$Cheatgrass - compare$Cheatgrass.Cover, compare$cwd_normal_cum, method = "pearson", + alternative = "two.sided", exact = FALSE) + +#linear model of log change ~ delta pdc and cwd +compare$log_cg_change <- log(I(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100 - min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100) + .01)) +cg_change_lm <- lm(log_cg_change ~ Delta_pdc*cwd_normal_cum, data= compare) + +summary(cg_change_lm) +plot(cg_change_lm) +plot(allEffects(cg_change_lm)) + + +#seeing how far out these outliers are +sd(compare$log_cg_change[-c(5,10)]) +mean(compare$log_cg_change[-c(5,10)]) +(compare$log_cg_change[c(5)] - 1.72) / .18 +(compare$log_cg_change[c(10)] - 1.72) / .18 + + +summary(cg_change_lm) +plot(cg_change_lm) +plot(allEffects(cg_change_lm, partial.residuals = TRUE)) + +#Perr grass +plot(compare$Pgrass ~ compare$Perrenial.Grass.Cover) +summary(lm(compare$Pgrass ~ compare$Perrenial.Grass.Cover + 0)) +cor.test(compare$Pgrass, compare$Perrenial.Grass.Cover, + method = c("spearman"), exact = FALSE) +boxplot(log(compare$Pgrass/compare$Perrenial.Grass.Cover)) +mean(compare$Pgrass-compare$Perrenial.Grass.Cover) +t.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") +wilcox.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") + +pg_change_lm <- lm(Pgrass - Perrenial.Grass.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +summary(pg_change_lm) + +## Perr forb +plot(I(compare$Aforb + compare$Pforb) ~ compare$Forb.Cover) +summary(lm(I(compare$Aforb + compare$Pforb) ~ compare$Forb.Cover + 0)) +cor.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, + method = c("spearman"), exact = FALSE) +boxplot(I(compare$Aforb + compare$Pforb) - compare$Forb.Cover) +mean(I(compare$Aforb + compare$Pforb) - compare$Forb.Cover) +t.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") +wilcox.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") + +pf_change_lm <- lm(Aforb + Pforb - Forb.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +summary(pf_change_lm) + +## Shrub +plot(compare$Shrub ~ I(compare$Shrub.Cover.Total)) + abline(0,1) +summary(lm(compare$Shrub ~ I(compare$Shrub.Cover.Total) + 0)) +cor.test(compare$Shrub, compare$Shrub.Cover.Total, + method = c("spearman"), exact = FALSE) +boxplot(compare$Shrub - compare$Shrub.Cover.Total) +mean(compare$Shrub - compare$Shrub.Cover.Total) +t.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="less") +wilcox.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="less") +hist(compare$Shrub - I(compare$Shrub.Cover.Total)) + +s_change_lm <- lm(Shrub - Shrub.Cover.Total ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +summary(s_change_lm) + +###################### +# Plot of comparisons between 2005 and 2015 +###################### +opar <- par(no.readonly = TRUE) + +png(filename="change_in_cover.png", + type="cairo", + units="in", + width = 4, + height=4, + pointsize=12, + res=160) + +layout(matrix(c(1,2,3,4), nrow = 2, ncol = 2, byrow = TRUE)) + +par(oma = c(2,2,0,0), family = "serif", bty = 'n') + + +par(mar = c(2,2,2,1)) +plot(compare$Pgrass ~ compare$Perrenial.Grass.Cover, xlim = c(0,35), ylim = c(0, 35), + xlab = "", ylab = "", main = "Perennial Grass", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,10,20,30), labels = TRUE) + axis(side = 1, at = c(0,5,15,25,35), tcl = -.24, labels = FALSE) + axis(side = 2, at = c(0,10,20,30), labels = TRUE) + axis(side = 2, at = c(0,5,15,25,35), tcl = -.24, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + +par(mar = c(2,2,2,1)) +plot(compare$Pforb ~ compare$Forb.Cover, xlim = c(0, 12), ylim = c(0, 12), + main = "Perennial Forb", xlab = "", ylab = "", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,4,8,12), labels = TRUE) + axis(side = 1, at = c(0,2,6,10), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,4,8,12), labels = TRUE) + axis(side = 2, at = c(0,2,6,10), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + + +par(mar = c(2,2,2,1)) +plot(compare$Shrub ~ compare$Shrub.Cover.Total, xlim = c(0, 30), ylim = c(0,30), + xlab = "", ylab = "", main = "Shrub",xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,10,20,30), labels = TRUE) + axis(side = 1, at = c(0,5,15,25), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,10,20,30), labels = TRUE) + axis(side = 2, at = c(0,5,15,25), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + + +par(mar = c(2,2,2,1)) +plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover, xlim = c(0, 20), ylim = c(0, 20), + xlab = "", ylab = "", main = "Cheatgrass", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 1, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 2, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + +mtext("2005 Cover (%)", side = 1, outer = TRUE) +mtext("2015 Cover (%)", side = 2, outer = TRUE) +dev.off() + +par(opar) + +#---------------------------------------- +# Ordination analysis +#---------------------------------------- +raw_cover <- read.csv("./Raw data/spp_cover2.csv") #spp_cover2 has unknowns filled with "unknown" +raw_cover <- raw_cover[raw_cover$Transect %in% c("E", "N", "S", "W"), ] #take out the "bonus" quadrats + +spp_summary <- read.csv("./Raw data/species_summary.csv") +spp_priority <- spp_summary[order(spp_summary[, "x"], decreasing = TRUE), ] + +species_list <- unique(raw_cover$Species) +plot_list <- unique(raw_cover$Plot)[order(unique(raw_cover$Plot))] + +spp_matrix <- matrix(ncol = length(species_list), nrow = length(plot_list)) +row.names(spp_matrix) <- plot_list +colnames(spp_matrix) <- species_list + +env_data <- plot_data[, c("cwd_normal_cum", "Delta_pdc", "AWC", "Live_ba", + "Dead_ba", "Avg_depth", "Elev", "Slope", "Swness", + "tmean", "ppt")] + + +for(i in 1:nrow(spp_matrix)){ + for(j in 1:ncol(spp_matrix)){ + + present <- species_list[j] %in% raw_cover[raw_cover$Plot == plot_list[i], "Species"] + if(present){ + plot_cover <- raw_cover[raw_cover$Plot == plot_list[i] & raw_cover$Species == species_list[j], ] + spp_matrix[i, j] <- sum(plot_cover$Cover) / 20 / 100 + } else{ + spp_matrix[i, j] <- 0 + } + + + } +} + + +nmds_test <- metaMDS(spp_matrix, try = 50) +nmds_test +plot(nmds_test, display = "species") +# text(nmds_test, display = "species") +orditorp(nmds_test, display = "species", priority = spp_priority$Group.1, + col = "forestgreen", pch = 2, cex = 1, air = 0.3) + +plot(envfit(nmds_test, env_data, labels = list(vectors = names(env_data)))) diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/DC1CB2D9-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/DC1CB2D9-contents new file mode 100644 index 0000000..6b8e6c7 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/DC1CB2D9-contents @@ -0,0 +1,397 @@ +library("MuMIn") +library("lme4") +library("car") +library("glmmADMB") +library("effects") +library("plyr") +library("lmerTest") +library("vegan") +library("reshape2") +library("plyr") + +setwd("C:/Users/Sam/Google Drive/Projects/MS Thesis/Understory") + +#-------------------------------------------------------- +##Understory data formatting +#-------------------------------------------------------- + + +daub <- read.csv("./Raw data/daub_cover.csv") + daub$Plot <- as.character(daub$Plot) + daub[daub$Plot == "SS2556120", "Plot"] <- "SS256120" + daub$Midpoint.value <- as.numeric(as.character((daub$Midpoint.value))) +trees <- read.csv("./Raw data/trees_updated_with_logs_041716.csv") +tree_pdc <- read.csv("./Raw data/all_trees_with_delta_and_ENN_041916.csv") +greenwood_under <- read.csv("./Raw data/Greenwood_Understory_Variables_SF_edits.csv") +tcover <- read.csv("./Raw data/tree_and_shrub_cover_020815.csv") +names(tcover)[3] <- "Shrub_cover_li" +source("calculate_awc.R") +soil <- calculate_awc("./Raw data/soils_missing_imputed_012016.csv") +clim <- read.csv("./Raw data/ALL_climate_variables.csv") +names(clim)[2] <- "Plot" +other_vars <- read.csv("./Raw data/all_vars_EXPORT.csv") +map_mat <- read.csv("./Raw data/map_mat_normals.csv") +names(map_mat)[2] <- "Plot" + + +#remaking data into plot average cover +plot_daub_cover <- dcast(daub, Plot ~ Cover.type, value.var = "Midpoint.value", fun.aggregate = mean) +names(plot_daub_cover) <- c("Plot", "Aforb", "Bg", "Cheatgrass", "Crust", "Gravel", "Litter", + "Agrass", "Pforb", "Pgrass", "Rock", "Shrub") + + +#------------------------------------------------------------- +#Create tree data plot-level + +#impute meanas for missing values +for (i in which(sapply(trees, is.numeric))) { + for (j in which(is.na(trees[, i]))) { + trees[j, i] <- mean(trees[trees[, "Plot"] == trees[j, "Plot"], i], na.rm = TRUE) + } +} + +plot_trees <- data.frame(Plot = character(102), + Live_ba = numeric(102), + Dead_ba = numeric(102), + Total_ba = numeric(102), + Dead_rat = numeric(102), + Delta_pdc = numeric(102)) + +plot_trees$Plot <- as.character(unique(plot_daub_cover$Plot)) + +for(i in (1:nrow(plot_trees))){ + plot_trees$Live_ba[i] <- sum(trees[trees$Live == "Y" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Dead_ba[i] <- sum(trees[trees$Live == "N" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Total_ba[i] <- sum(trees[trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, na.rm = TRUE) + plot_trees$Dead_rat[i] <- plot_trees$Dead_ba[i] / plot_trees$Total_ba[i] + plot_trees$Delta_pdc[i] <- weighted.mean(tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "Delta_pdc"], + tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "BA_cm"]) +} + +#----------------------------------------------- +# Plot-level cover for everything and also predictor variables + +#fill in any NAs + +plot_data <- join(plot_daub_cover, plot_trees, by = c("Plot")) +plot_data <- join(plot_data, clim[, c("Plot", "cwd_normal_cum", "avg_ann_ppt")], by = "Plot") +plot_data <- join(plot_data, soil[, c("Plot", "AWC")], by = "Plot") +plot_data <- join(plot_data, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot") +plot_data <- join(plot_data, other_vars[, c("Plot", "Cluster", "Avg_depth", "Elev", "Slope", "Swness")], by = "Plot") +plot_data <- join(plot_data, map_mat[, c("Plot", "tmean", "ppt")], by = "Plot") + +for (i in which(sapply(plot_data, is.numeric))) { + for (j in which(is.na(plot_data[, i]))) { + plot_data[j, i] <- mean(plot_data[, i], na.rm = TRUE) + } +} + +#------------------------------------------------------------ +# Analysis +#------------------------------------------------------------ +# +levels(plot_data$Cluster) <- c(levels(plot_data$Cluster), "NPELECTRICEEL") +plot_data[(plot_data$Plot %in% c("NPELECTRICEEL", "NPELECTRICEEL120", + "NPELECTRICEEL240", "NPELECTRICEEL360")), "Cluster"] <- "NPELECTRICEEL" + + +#exploration +hist(asin(sqrt(plot_data$Cheatgrass/100))) +hist(logit(plot_data$Cheatgrass/100 + min(plot_data[plot_data$Cheatgrass != 0, "Cheatgrass"]))) +hist(log(plot_data$Cheatgrass/100 + .01)) + +nrow(plot_data[plot_data$Cheatgrass != 0, ]) + +hist(asin(sqrt(plot_data$Pgrass/100))) +hist(logit(plot_data$Pgrass/100)) +nrow(plot_data[plot_data$Pgrass != 0, ]) + +hist(asin(sqrt(plot_data$Pforb/100))) +hist(logit(plot_data$Pforb/100)) +nrow(plot_data[plot_data$Pforb != 0, ]) + +hist(asin(sqrt(plot_data$Shrub_cover_li))) +hist(logit(plot_data$Shrub_cover_li)) +nrow(plot_data[plot_data$Shrub_cover_li != 0, ]) + +plot(plot_data$Cheatgrass ~ plot_data$cwd_normal_cum) + + +#--------------------------------------------------------------- +#linear mixed effects models for functional groups at plot level +#------------------------------------------------------------------- +## Cheatgrass + +cheatgrass_plot <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +cheatgrass_plot <- lmer(logit(Cheatgrass/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +summary(lm(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC), data = plot_data)) + +summary(cheatgrass_plot, ddf = "Kenward-Roger") +r.squaredGLMM(cheatgrass_plot) +plot(allEffects(cheatgrass_plot, partial.residuals = TRUE)) +plot(cheatgrass_plot) +AICc(cheatgrass_plot) +scatter.smooth(residuals(cheatgrass_plot) ~ predict(cheatgrass_plot)) + +leveneTest(residuals(cheatgrass_plot) ~ plot_data$Cheatgrass) + +## Perr grass + +pgrass_plot <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +summary(lm(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC), data = plot_data)) + +summary(pgrass_plot, ddf = "Kenward-Roger") +r.squaredGLMM(pgrass_plot) +plot(allEffects(pgrass_plot, partial.residuals = TRUE)) +plot(pgrass_plot) + +## Perr forbs + +pforb_plot <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +summary(lm(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC), data = plot_data)) + + +summary(pforb_plot, ddf = "Kenward-Roger") +r.squaredGLMM(pforb_plot) +plot(allEffects(pforb_plot, partial.residuals = TRUE)) +plot(pforb_plot) + + +## Shrubs +shrub_plot <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +summary(lm(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC), data = plot_data)) + +summary(shrub_plot, ddf = "Kenward-Roger") +r.squaredGLMM(shrub_plot) +plot(allEffects(shrub_plot, partial.residuals = TRUE)) +plot(shrub_plot) + + +#----------------------------------------------- +# Comparison of greenwood and 2015 +#----------------------------------------------- +compare <- join(plot_daub_cover, greenwood_under, by = "Plot", type = "inner") +compare <- join(compare, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot", type = "inner") +compare <- join(compare, plot_data[, c("Plot", "Dead_ba", "Delta_pdc", "cwd_normal_cum")], by = "Plot") + +#cheatgrass +plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover) +summary(lm(compare$Cheatgrass ~ compare$Cheatgrass.Cover)) +cor.test(compare$Cheatgrass, compare$Cheatgrass.Cover, + method = c("spearman"), exact = FALSE) +boxplot(compare$Cheatgrass, compare$Cheatgrass.Cover) +mean(compare$Cheatgrass-compare$Cheatgrass.Cover) +wilcox.test(compare$Cheatgrass, compare$Cheatgrass.Cover, paired = TRUE, alternative="greater") +t.test(compare$Cheatgrass, compare$Cheatgrass.Cover, paired = TRUE, alternative="greater") + +cor.test(compare$Cheatgrass - compare$Cheatgrass.Cover, compare$cwd_normal_cum, method = "pearson", + alternative = "two.sided", exact = FALSE) + +#linear model of log change ~ delta pdc and cwd +compare$log_cg_change <- log(I(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100 - min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100) + .01)) +cg_change_lm <- lm(log_cg_change ~ Delta_pdc*cwd_normal_cum, data= compare) + +summary(cg_change_lm) +plot(cg_change_lm) +plot(allEffects(cg_change_lm)) + + +#seeing how far out these outliers are +sd(compare$log_cg_change[-c(5,10)]) +mean(compare$log_cg_change[-c(5,10)]) +(compare$log_cg_change[c(5)] - 1.72) / .18 +(compare$log_cg_change[c(10)] - 1.72) / .18 + + +summary(cg_change_lm) +plot(cg_change_lm) +plot(allEffects(cg_change_lm, partial.residuals = TRUE)) + +#Perr grass +plot(compare$Pgrass ~ compare$Perrenial.Grass.Cover) +summary(lm(compare$Pgrass ~ compare$Perrenial.Grass.Cover + 0)) +cor.test(compare$Pgrass, compare$Perrenial.Grass.Cover, + method = c("spearman"), exact = FALSE) +boxplot(log(compare$Pgrass/compare$Perrenial.Grass.Cover)) +mean(compare$Pgrass-compare$Perrenial.Grass.Cover) +t.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") +wilcox.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") + +pg_change_lm <- lm(Pgrass - Perrenial.Grass.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +summary(pg_change_lm) + +## Perr forb +plot(I(compare$Aforb + compare$Pforb) ~ compare$Forb.Cover) +summary(lm(I(compare$Aforb + compare$Pforb) ~ compare$Forb.Cover + 0)) +cor.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, + method = c("spearman"), exact = FALSE) +boxplot(I(compare$Aforb + compare$Pforb) - compare$Forb.Cover) +mean(I(compare$Aforb + compare$Pforb) - compare$Forb.Cover) +t.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") +wilcox.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") + +pf_change_lm <- lm(Aforb + Pforb - Forb.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +summary(pf_change_lm) + +## Shrub +plot(compare$Shrub ~ I(compare$Shrub.Cover.Total)) + abline(0,1) +summary(lm(compare$Shrub ~ I(compare$Shrub.Cover.Total) + 0)) +cor.test(compare$Shrub, compare$Shrub.Cover.Total, + method = c("spearman"), exact = FALSE) +boxplot(compare$Shrub - compare$Shrub.Cover.Total) +mean(compare$Shrub - compare$Shrub.Cover.Total) +t.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="less") +wilcox.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="less") +hist(compare$Shrub - I(compare$Shrub.Cover.Total)) + +s_change_lm <- lm(Shrub - Shrub.Cover.Total ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +summary(s_change_lm) + +###################### +# Plot of comparisons between 2005 and 2015 +###################### +opar <- par(no.readonly = TRUE) + +png(filename="change_in_cover.png", + type="cairo", + units="in", + width = 4, + height=4, + pointsize=12, + res=160) + +layout(matrix(c(1,2,3,4), nrow = 2, ncol = 2, byrow = TRUE)) + +par(oma = c(2,2,0,0), family = "serif", bty = 'n') + + +par(mar = c(2,2,2,1)) +plot(compare$Pgrass ~ compare$Perrenial.Grass.Cover, xlim = c(0,35), ylim = c(0, 35), + xlab = "", ylab = "", main = "Perennial Grass", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,10,20,30), labels = TRUE) + axis(side = 1, at = c(0,5,15,25,35), tcl = -.24, labels = FALSE) + axis(side = 2, at = c(0,10,20,30), labels = TRUE) + axis(side = 2, at = c(0,5,15,25,35), tcl = -.24, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + +par(mar = c(2,2,2,1)) +plot(compare$Pforb ~ compare$Forb.Cover, xlim = c(0, 12), ylim = c(0, 12), + main = "Perennial Forb", xlab = "", ylab = "", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,4,8,12), labels = TRUE) + axis(side = 1, at = c(0,2,6,10), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,4,8,12), labels = TRUE) + axis(side = 2, at = c(0,2,6,10), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + + +par(mar = c(2,2,2,1)) +plot(compare$Shrub ~ compare$Shrub.Cover.Total, xlim = c(0, 30), ylim = c(0,30), + xlab = "", ylab = "", main = "Shrub",xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,10,20,30), labels = TRUE) + axis(side = 1, at = c(0,5,15,25), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,10,20,30), labels = TRUE) + axis(side = 2, at = c(0,5,15,25), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + + +par(mar = c(2,2,2,1)) +plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover, xlim = c(0, 20), ylim = c(0, 20), + xlab = "", ylab = "", main = "Cheatgrass", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 1, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 2, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + +mtext("2005 Cover (%)", side = 1, outer = TRUE) +mtext("2015 Cover (%)", side = 2, outer = TRUE) +dev.off() + +par(opar) + +#---------------------------------------- +# Ordination analysis +#---------------------------------------- +raw_cover <- read.csv("./Raw data/spp_cover2.csv") #spp_cover2 has unknowns filled with "unknown" +raw_cover <- raw_cover[raw_cover$Transect %in% c("E", "N", "S", "W"), ] #take out the "bonus" quadrats + +spp_summary <- read.csv("./Raw data/species_summary.csv") +spp_priority <- spp_summary[order(spp_summary[, "x"], decreasing = TRUE), ] + +species_list <- unique(raw_cover$Species) +plot_list <- unique(raw_cover$Plot)[order(unique(raw_cover$Plot))] + +spp_matrix <- matrix(ncol = length(species_list), nrow = length(plot_list)) +row.names(spp_matrix) <- plot_list +colnames(spp_matrix) <- species_list + +env_data <- plot_data[, c("cwd_normal_cum", "Delta_pdc", "AWC", "Live_ba", + "Dead_ba", "Avg_depth", "Elev", "Slope", "Swness", + "tmean", "ppt")] + + +for(i in 1:nrow(spp_matrix)){ + for(j in 1:ncol(spp_matrix)){ + + present <- species_list[j] %in% raw_cover[raw_cover$Plot == plot_list[i], "Species"] + if(present){ + plot_cover <- raw_cover[raw_cover$Plot == plot_list[i] & raw_cover$Species == species_list[j], ] + spp_matrix[i, j] <- sum(plot_cover$Cover) / 20 / 100 + } else{ + spp_matrix[i, j] <- 0 + } + + + } +} + + +nmds_test <- metaMDS(spp_matrix, try = 50) +nmds_test +plot(nmds_test, display = "species") +# text(nmds_test, display = "species") +orditorp(nmds_test, display = "species", priority = spp_priority$Group.1, + col = "forestgreen", pch = 2, cex = 1, air = 0.3) + +plot(envfit(nmds_test, env_data, labels = list(vectors = names(env_data)))) diff --git a/.Rproj.user/EAEB7255/sources/per/t/8F5BFB6A b/.Rproj.user/EAEB7255/sources/s-A74A5E15/E666EB0A similarity index 63% rename from .Rproj.user/EAEB7255/sources/per/t/8F5BFB6A rename to .Rproj.user/EAEB7255/sources/s-A74A5E15/E666EB0A index 929d756..8f99c7a 100644 --- a/.Rproj.user/EAEB7255/sources/per/t/8F5BFB6A +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/E666EB0A @@ -1,21 +1,21 @@ { "collab_server" : "", "contents" : "", - "created" : 1536264244143.000, + "created" : 1551911553391.000, "dirty" : false, "encoding" : "UTF-8", "folds" : "", - "hash" : "1779826032", - "id" : "8F5BFB6A", + "hash" : "0", + "id" : "E666EB0A", "lastKnownWriteTime" : 1547069446, - "last_content_update" : 1547069446252, + "last_content_update" : 1547069446, "path" : "~/Research/MS Thesis/Understory/understory_data_prep.R", "project_path" : "understory_data_prep.R", "properties" : { - "cursorPosition" : "94,50", - "scrollLine" : "20" + "cursorPosition" : "23,28", + "scrollLine" : "17" }, - "relative_order" : 6, + "relative_order" : 3, "source_on_save" : false, "source_window" : "", "type" : "r_source" diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/E666EB0A-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/E666EB0A-contents new file mode 100644 index 0000000..1630074 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/E666EB0A-contents @@ -0,0 +1,95 @@ +##Understory data formatting +## +library("reshape2") +library("plyr") +library("effects") + + +daub <- read.csv("./raw data/daub_cover.csv", stringsAsFactors = FALSE) + +#some data proofing +daub <- daub[(daub$Transect %in% c("N", "E", "S", "W")), ] +daub$Plot <- as.character(daub$Plot) +daub[daub$Plot == "NPElectricEel", "Plot"] <- "NPELECTRICEEL" +daub[daub$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120" +daub[daub$Plot == "NPElectricEel240", "Plot"] <- "NPELECTRICEEL240" +daub[daub$Plot == "NPElectricEel360", "Plot"] <- "NPELECTRICEEL360" +daub$Midpoint.value <- as.numeric(as.character((daub$Midpoint.value))) +daub$unique_quad <- paste0(daub$Plot, daub$Transect, daub$Meter) +daub[daub$Cover.type == "Perennial forb ", "Cover.type"] <- "Perennial forb" +daub[daub$Cover.type == "Shrub ", "Cover.type"] <- "Shrub" + +trees <- read.csv("./raw data/trees_updated_with_logs_041716.csv") +tree_pdc <- read.csv("./raw data/all_trees_with_delta_and_ENN_041916.csv") +greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv") +tcover <- read.csv("./raw data/tree_and_shrub_cover_020815.csv") + names(tcover)[3] <- "Shrub_cover_li" +source("./calculate_awc.R") +soil <- calculate_awc(soil_raw = "./raw data/soils_missing_imputed_012016.csv") +clim <- read.csv("./raw data/ALL_climate_variables.csv") + names(clim)[2] <- "Plot" +other_vars <- read.csv("./raw data/all_vars_EXPORT.csv") + + +#remaking data into plot average cover +plot_daub_cover <- dcast(daub, Plot ~ Cover.type, value.var = "Midpoint.value", fun.aggregate = mean) +names(plot_daub_cover) <- c("Plot", "Aforb", "Bg", "Cheatgrass", "Crust", "Gravel", "Litter", + "Agrass", "Pforb", "Pgrass", "Rock", "Shrub") +plot_daub_cover$All <- plot_daub_cover$Aforb + plot_daub_cover$Cheatgrass + plot_daub_cover$Agrass + plot_daub_cover$Pforb + + plot_daub_cover$Pgrass + plot_daub_cover$Shrub + +#------------------------------------------------------------- +#Create tree data plot-level + +#impute meanas for missing values +for (i in which(sapply(trees, is.numeric))) { + for (j in which(is.na(trees[, i]))) { + trees[j, i] <- mean(trees[trees[, "Plot"] == trees[j, "Plot"], i], na.rm = TRUE) + } +} + +plot_trees <- data.frame(Plot = character(102), + Live_ba = numeric(102), + Dead_ba = numeric(102), + Total_ba = numeric(102), + Dead_rat = numeric(102), + Delta_pdc = numeric(102)) + +plot_trees$Plot <- as.character(unique(plot_daub_cover$Plot)) + +for(i in (1:nrow(plot_trees))){ + plot_trees$Live_ba[i] <- sum(trees[trees$Live == "Y" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Dead_ba[i] <- sum(trees[trees$Live == "N" & trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, + na.rm = TRUE) + plot_trees$Total_ba[i] <- sum(trees[trees$Plot == as.character(plot_trees$Plot[i]), ]$BA_cm, na.rm = TRUE) + plot_trees$Dead_rat[i] <- plot_trees$Dead_ba[i] / plot_trees$Total_ba[i] + plot_trees$Delta_pdc[i] <- weighted.mean(tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "Delta_pdc"], + tree_pdc[tree_pdc$Plot == as.character(plot_trees$Plot[i]), "BA_cm"]) +} + + +#----------------------------------------------- +# Plot-level cover for everything and also predictor variables + +#fill in any NAs + +plot_data <- join(plot_daub_cover, plot_trees, by = c("Plot")) +plot_data <- join(plot_data, clim[, c("Plot", "cwd_normal_cum")], by = "Plot") +plot_data <- join(plot_data, soil[, c("Plot", "AWC")], by = "Plot") +plot_data <- join(plot_data, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot") +plot_data <- join(plot_data, other_vars[, c("Plot", "Cluster", "Avg_depth")], by = "Plot") + + +for (i in which(sapply(plot_data, is.numeric))) { + for (j in which(is.na(plot_data[, i]))) { + plot_data[j, i] <- mean(plot_data[, i], na.rm = TRUE) + } +} + +levels(plot_data$Cluster) <- c(levels(plot_data$Cluster), "NPELECTRICEEL") +plot_data[(plot_data$Plot %in% c("NPELECTRICEEL", "NPELECTRICEEL120", + "NPELECTRICEEL240", "NPELECTRICEEL360")), "Cluster"] <- "NPELECTRICEEL" + + +write.csv(plot_data, "./clean data/plot_data.csv") \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/F8B392EE b/.Rproj.user/EAEB7255/sources/s-A74A5E15/F8B392EE new file mode 100644 index 0000000..a70d475 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/F8B392EE @@ -0,0 +1,22 @@ +{ + "collab_server" : "", + "contents" : "", + "created" : 1547230037453.000, + "dirty" : false, + "encoding" : "UTF-8", + "folds" : "", + "hash" : "0", + "id" : "F8B392EE", + "lastKnownWriteTime" : 1547472910, + "last_content_update" : 1547472910192, + "path" : "~/Research/MS Thesis/Understory/understory_effects_plots.R", + "project_path" : "understory_effects_plots.R", + "properties" : { + "cursorPosition" : "13,21", + "scrollLine" : "4" + }, + "relative_order" : 2, + "source_on_save" : false, + "source_window" : "", + "type" : "r_source" +} \ No newline at end of file diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/F8B392EE-contents b/.Rproj.user/EAEB7255/sources/s-A74A5E15/F8B392EE-contents new file mode 100644 index 0000000..7d8e0c1 --- /dev/null +++ b/.Rproj.user/EAEB7255/sources/s-A74A5E15/F8B392EE-contents @@ -0,0 +1,447 @@ +### Code to make effects plots + +# TODO scheck what models are fit, fit them in another script and import them to here + +library(visreg) +library(car) +library(effects) +library(lme4) + +pardefault <- par(no.readonly = T) +source('addTrans.R', echo=FALSE) + +### set some parameters to be use throughout +cg <- cheatgrass_plot +pg <- pgrass_plot +pf <- pforb_plot +sh <- shrub_plot + +# function to calculate standard errors for a given vcov matrix and row of x values, to get +# point estimates of prediction error for a given prediction. See http://www.ats.ucla.edu/stat/r/faq/deltamethod.htm +# This gets used for the confidence interval of the effects plot + +calc.se <- function(x){ + sqrt(as.matrix(x) %*% vcov(model) %*% t(x)) +} + +inv.as <- function(x){ + (sin(x))^2 +} + +############################################## +## Creating new predictions +############################################## +model <- cg + + +# Calculate residuals from predicted values and observed values +resid <- residuals(model) + +# for unscaling variables -- these things come from a full model in the model set, you'll have to change +# the number or the path depending upon what kind of model object and the order of your models, etc. +# You can also just get them from the raw data. +tc_mean <- mean(plot_data$Tree_cover) +tc_sd <- sd(plot_data$Tree_cover) + +cwd_mean <- mean(plot_data$cwd_normal_cum) +cwd_sd <- sd(plot_data$cwd_normal_cum) + +Delta_pdc_mean <- mean(plot_data$Delta_pdc) +Delta_pdc_sd <- sd(plot_data$Delta_pdc) + +################################ +## Calcuate partial residuals, effects, for avg_BA +################################ + +model <- cg + +C_tc <- data.frame( + Tree_cover = seq(min(plot_data$Tree_cover), max(plot_data$Tree_cover), length.out = 102), + cwd_normal_cum = mean(plot_data$cwd_normal_cum), + Delta_pdc = mean(plot_data$Delta_pdc), + AWC = mean(plot_data$AWC) +) + +C_awc <- data.frame( + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = mean(plot_data$cwd_normal_cum), + Delta_pdc = mean(plot_data$Delta_pdc), + AWC = seq(min(plot_data$AWC), max(plot_data$AWC), length.out = 102) +) + +C_pdc_10 <- data.frame( + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .1)), + Delta_pdc = seq(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc), length.out = 102), + AWC = mean(plot_data$AWC) +) + +C_pdc_50 <- data.frame( + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .5)), + Delta_pdc = seq(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc), length.out = 102), + AWC = mean(plot_data$AWC) +) + +C_pdc_90 <- data.frame( + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .9)), + Delta_pdc = seq(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc), length.out = 102), + AWC = mean(plot_data$AWC) +) + + +# model <- cg +# C <- C_tc +# predictor <- "Tree_cover" +calculate_preds <- function(model, C, predictor){ + preds <- data.frame( + predictions = predict(model, re.form = NA, newdata = C, type = "response"), + predictor = C[, predictor] + ) + + return(preds) + +} + +# Use the function to get all the stuff we need +# (predictions, partial residuals, and SEs) +# this is the dumbest way to do this but I don't want to rewrite it +preds_tc_cg <- calculate_preds(cg, C_tc, "Tree_cover") +preds_tc_pg <- calculate_preds(pg, C_tc, "Tree_cover") +preds_tc_pf <- calculate_preds(pf, C_tc, "Tree_cover") +preds_tc_sh <- calculate_preds(sh, C_tc, "Tree_cover") + +preds_awc_cg <- calculate_preds(cg, C_awc, "AWC") +preds_awc_pg <- calculate_preds(pg, C_awc, "AWC") +preds_awc_pf <- calculate_preds(pf, C_awc, "AWC") +preds_awc_sh <- calculate_preds(sh, C_awc, "AWC") + +preds_pdc_10_cg <- calculate_preds(cg, C_pdc_10, "Delta_pdc") +preds_pdc_10_pg <- calculate_preds(pg, C_pdc_10, "Delta_pdc") +preds_pdc_10_pf <- calculate_preds(pf, C_pdc_10, "Delta_pdc") +preds_pdc_10_sh <- calculate_preds(sh, C_pdc_10, "Delta_pdc") + +preds_pdc_50_cg <- calculate_preds(cg, C_pdc_50, "Delta_pdc") +preds_pdc_50_pg <- calculate_preds(pg, C_pdc_50, "Delta_pdc") +preds_pdc_50_pf <- calculate_preds(pf, C_pdc_50, "Delta_pdc") +preds_pdc_50_sh <- calculate_preds(sh, C_pdc_50, "Delta_pdc") + +preds_pdc_90_cg <- calculate_preds(cg, C_pdc_90, "Delta_pdc") +preds_pdc_90_pg <- calculate_preds(pg, C_pdc_90, "Delta_pdc") +preds_pdc_90_pf <- calculate_preds(pf, C_pdc_90, "Delta_pdc") +preds_pdc_90_sh <- calculate_preds(sh, C_pdc_90, "Delta_pdc") + +# create plots of log-partial-residuals y-axis versus natural x-scale +# this is also dumb and I could do it in a loop except I don't want to +# put all the predictions in a list +setwd("C:/Users/Sam/Documents/Research/MS Thesis/Understory/outputs") #where the plots go + +opar <- par(no.readonly = TRUE) +par(opar) + +png(filename="understory_effects.png", + type="cairo", + units="in", + width = 7, + height=5, + pointsize=15, + res=160) + +layout(matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3, byrow = TRUE)) + +par(oma = c(2,4,0,0), mar = c(3,1,1,1), family = "serif", bty = 'n') + + +plot(NA, + ylim = c(0,.1), + xlim = c(I(min(plot_data$Tree_cover)*100), I(max(plot_data$Tree_cover)*100)), + xlab = "", + ylab = "", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5, + bty = 'n', + xaxt = 'n', + yaxt = 'n') +lines(inv.as(preds_tc_cg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_tc_pg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_tc_pf$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_tc_sh$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 4, col = "#e7298a") +axis(side = 1) +axis(side = 2) +mtext(text = "Tree Cover (%)", side = 1, line = 2.2) +mtext(text = "Predicted Cover", side = 2, outer = TRUE, line = 2, cex = 1) +mtext(text = "(a)", side = 1, line = -10, adj = 0.05) + + +plot(NA, + ylim = c(0,.1), + xlim = c(min(plot_data$AWC), max(plot_data$AWC)), + xlab = "", + ylab = "", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5, + bty = 'n', + xaxt = 'n', + yaxt = 'n') +lines(inv.as(preds_awc_cg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_awc_pg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_awc_pf$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_awc_sh$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +axis(side = 1) +axis(side = 2) +mtext(text = "Soil AWC (%)", side = 1, line = 2.2) +mtext(text = "(b)", side = 1, line = -10, adj = 0.05) + + +plot.new() +legend("topright", legend = c("Cheatgrass", "Perr. Grass", "Perr. Forb", "Shrub"), + lty = c(1,2,3,4), lwd = 2, cex = 1.3, col = c("#1b9e77", "#d95f02", "#7570b3", "#e7298a")) + +plot(NA, + ylim = c(0,.2), + xlim = c(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc)), + xlab = "", + ylab = "", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5, + bty = 'n', + xaxt = 'n', + yaxt = 'n') +lines(inv.as(preds_pdc_10_cg$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_pdc_10_pg$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_pdc_10_pf$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_pdc_10_sh$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +axis(side = 1) +axis(side = 2) +text(x = -30, y = .09, labels = "10% CWD", cex = 1.3) +mtext(text = "(c)", side = 1, line = -10, adj = 0.05) + +plot(NA, + ylim = c(0,.2), + xlim = c(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc)), + xlab = "", + ylab = "", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5, + bty = 'n', + xaxt = 'n', + yaxt = 'n') +lines(inv.as(preds_pdc_50_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_pdc_50_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_pdc_50_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_pdc_50_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +axis(side = 1) +axis(side = 2) +text(x = -30, y = .09, labels = "50% CWD", cex = 1.3) +mtext(text = "Change in Live Canopy (%)", side = 1, line = 2.2) +mtext(text = "(d)", side = 1, line = -10, adj = 0.05) + + +plot(NA, + ylim = c(0,.3), + xlim = c(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc)), + xlab = "", + ylab = "", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5, + bty = 'n', + xaxt = 'n', + yaxt = 'n') +lines(inv.as(preds_pdc_90_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_pdc_90_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_pdc_90_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_pdc_90_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +axis(side = 1) +axis(side = 2) +text(x = -30, y = .09, labels = "90% CWD", cex = 1.3) +mtext(text = "(e)", side = 1, line = -10, adj = 0.05) + + +dev.off() + +par(pardefault) + + +#---------------------------------------------------------------------------------- +# Cheatgrass change over time +#---------------------------------------------------------------------------------- +## new data +model <- + +C_cg_change_10 <- data.frame( + cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .1)), + Delta_pdc = seq(-60,15, length.out = 100) +) + + +C_cg_change_50 <- data.frame( + cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .5)), + Delta_pdc = seq(-60,15, length.out = 100) +) + + +C_cg_change_90 <- data.frame( + cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .9)), + Delta_pdc = seq(-60,15, length.out = 100) +) + +#calculate partial residuals for different levels of cwd + +part_resids_10 <- residuals(model) + model$model[, "Delta_pdc"] * coef(model)["Delta_pdc"] + + unname(quantile(compare$cwd_normal_cum, .1)) * coef(model)["cwd_normal_cum"] + + model$model[, "Delta_pdc"] * unname(quantile(compare$cwd_normal_cum, .1)) * coef(model)["Delta_pdc:cwd_normal_cum"] + + coef(model)["(Intercept)"] + +part_resids_50 <- residuals(model) + model$model[, "Delta_pdc"] * coef(model)["Delta_pdc"] + + unname(quantile(compare$cwd_normal_cum, .5)) * coef(model)["cwd_normal_cum"] + + model$model[, "Delta_pdc"] * unname(quantile(compare$cwd_normal_cum, .5)) * coef(model)["Delta_pdc:cwd_normal_cum"] + + coef(model)["(Intercept)"] + +part_resids_90 <- residuals(model) + model$model[, "Delta_pdc"] * coef(model)["Delta_pdc"] + + unname(quantile(compare$cwd_normal_cum, .9)) * coef(model)["cwd_normal_cum"] + + model$model[, "Delta_pdc"] * unname(quantile(compare$cwd_normal_cum, .9)) * coef(model)["Delta_pdc:cwd_normal_cum"] + + coef(model)["(Intercept)"] + + +#function to calculate predictions +preds_cg_change <- function(model, C, predictor){ + preds <- data.frame( + predictions = predict(model, se = TRUE, newdata = C, type = "response")$fit, + lo = predict(model, se = TRUE, newdata = C)$fit - 1.96*predict(model, se = TRUE, newdata = C)$se.fit, + hi = predict(model, se = TRUE, newdata = C)$fit + 1.96*predict(model, se = TRUE, newdata = C)$se.fit, + predictor = C[, predictor] + ) + + return(preds) +} + +cg_change_10 <- preds_cg_change(cg_change_lm, C_cg_change_10, "Delta_pdc") +cg_change_50 <- preds_cg_change(cg_change_lm, C_cg_change_50, "Delta_pdc") +cg_change_90 <- preds_cg_change(cg_change_lm, C_cg_change_90, "Delta_pdc") + + + +setwd("C:\\Users\\Sam\\Google Drive\\Thesis related\\Understory\\outputs") #where the plots go + +png(filename="cheatgrass_change_effect.png", + type="cairo", + units="in", + width = 5, + height=3, + pointsize=15, + res=160) + +layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE)) +par(mar = c(0,1,0,0), oma = c(3,4,2,2), family = "serif", bty = 'n') + + plot(NA, + ylim = c(-.03,.06), + xlim = c(-50,10), + xlab = "", + ylab = "", + xaxt = "n", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5) + abline(h = 0, lty = 4, lwd = 1.3) + lines((inv.as(cg_change_10$predictions) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 2, lty = 1, col = "#1b9e77") + lines((inv.as(cg_change_10$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + lines((inv.as(cg_change_10$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + polygon(c(cg_change_10$predictor, rev(cg_change_10$predictor)), c((exp(cg_change_10$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)), + rev((exp(cg_change_10$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)))), + col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals + + points((exp(part_resids_10[model$model[, "cwd_normal_cum"] < unname(quantile(compare$cwd_normal_cum, .3))]) + + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100))~ + model$model[model$model[, "cwd_normal_cum"] < unname(quantile(compare$cwd_normal_cum, .3)), "Delta_pdc"], + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5) + axis(side = 1, at = c(-40, -20, 0), labels = TRUE, mgp=c(1.5,-0.5,-1.5)) + axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25, mgp=c(1.5,-0.5,-1.5)) + text(x = -20, y = .06, labels = "10% CWD") + + + plot(NA, + ylim = c(-.03,.06), + xlim = c(-50,10), + xlab = "", + ylab = "", + yaxt = "n", + xaxt = "n", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5) + abline(h = 0, lty = 4, lwd = 1.3) + lines((exp(cg_change_50$predictions) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 2, lty = 1, col = "#1b9e77") + lines((exp(cg_change_50$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + lines((exp(cg_change_50$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + polygon(c(cg_change_50$predictor, rev(cg_change_50$predictor)), c((exp(cg_change_50$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)), + rev((exp(cg_change_50$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)))), + col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals + + + points((exp(part_resids_50[model$model[, "cwd_normal_cum"] > unname(quantile(compare$cwd_normal_cum, .3)) & + model$model[, "cwd_normal_cum"] < unname(quantile(compare$cwd_normal_cum, .7))]) + + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ + model$model[model$model[, "cwd_normal_cum"] > unname(quantile(compare$cwd_normal_cum, .3)) & + model$model[, "cwd_normal_cum"] < unname(quantile(compare$cwd_normal_cum, .7)), "Delta_pdc"], + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5) + axis(side = 2, labels = FALSE) + axis(side = 1, at = c(-40, -20, 0), labels = TRUE, mgp=c(1.5,-0.5,-1.5)) + axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25, mgp=c(1.5,-0.5,-1.5)) + text(x = -20, y = .06, labels = "50% CWD") + + + plot(NA, + ylim = c(-.03,.06), + xlim = c(-50,10), + xlab = "", + ylab = "", + yaxt = "n", + xaxt = "n", + bg='grey60', + col = 'grey30', + pch = 21, + cex.lab = 1.5) + abline(h = 0, lty = 4, lwd = 1.3) + lines((exp(cg_change_90$predictions) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 2, lty = 1, col = "#1b9e77") + lines((exp(cg_change_90$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + lines((exp(cg_change_90$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + polygon(c(cg_change_90$predictor, rev(cg_change_90$predictor)), c((exp(cg_change_90$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)), + rev((exp(cg_change_90$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)))), + col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals + + points((exp(part_resids_90[model$model[, "cwd_normal_cum"] > unname(quantile(compare$cwd_normal_cum, .7))]) + + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ + model$model[model$model[, "cwd_normal_cum"] > unname(quantile(compare$cwd_normal_cum, .7)), "Delta_pdc"], + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5) + axis(side = 2, labels = FALSE) + axis(side = 1, at = c(-40, -20, 0), labels = TRUE, mgp=c(1.5,-0.5,-1.5)) + axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25, mgp=c(1.5,-0.5,-1.5)) + text(x = -20, y = .06, labels = "90% CWD") + + mtext(side = 1, text = "Change in Live Crown", outer = TRUE, line = 1.7) + mtext(side = 2, text = "Pred. Change Cheatgrass Cover", outer = TRUE, line = 2.3) + +dev.off() diff --git a/.Rproj.user/EAEB7255/sources/s-A74A5E15/lock_file b/.Rproj.user/EAEB7255/sources/s-A74A5E15/lock_file new file mode 100644 index 0000000..e69de29 diff --git a/.Rproj.user/EAEB7255/viewer-cache/29B3C11D.Rdata b/.Rproj.user/EAEB7255/viewer-cache/29B3C11D.Rdata deleted file mode 100644 index 4cca1bd..0000000 Binary files a/.Rproj.user/EAEB7255/viewer-cache/29B3C11D.Rdata and /dev/null differ diff --git a/.Rproj.user/EAEB7255/viewer-cache/6F672F8B.Rdata b/.Rproj.user/EAEB7255/viewer-cache/6F672F8B.Rdata deleted file mode 100644 index ea663f3..0000000 Binary files a/.Rproj.user/EAEB7255/viewer-cache/6F672F8B.Rdata and /dev/null differ diff --git a/.Rproj.user/EAEB7255/viewer-cache/71129F7A.Rdata b/.Rproj.user/EAEB7255/viewer-cache/71129F7A.Rdata deleted file mode 100644 index 196938b..0000000 Binary files a/.Rproj.user/EAEB7255/viewer-cache/71129F7A.Rdata and /dev/null differ diff --git a/.Rproj.user/EAEB7255/viewer-cache/EDD72D8A.Rdata b/.Rproj.user/EAEB7255/viewer-cache/EDD72D8A.Rdata deleted file mode 100644 index ad8bdda..0000000 Binary files a/.Rproj.user/EAEB7255/viewer-cache/EDD72D8A.Rdata and /dev/null differ diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths index fab7c01..fef631f 100644 --- a/.Rproj.user/shared/notebooks/paths +++ b/.Rproj.user/shared/notebooks/paths @@ -1,5 +1,8 @@ C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory.R="3F501BF8" +C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory_082418.R="86948D8" C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory_090618.R="1D3078D6" +C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory_112817.R="3740AC65" +C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_data_prep.R="318C3941" C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_electivity.R="A106E7C2" C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_electivity010819.R="EAB749D7" C:/Users/Sam/Google Drive/Teaching/Spring 2018/Lab 11 - Population/beetle_script.R="BCD3989" diff --git a/Rplots.pdf b/Rplots.pdf index 5b5f423..437f719 100644 Binary files a/Rplots.pdf and b/Rplots.pdf differ diff --git a/change_in_cover.png b/change_in_cover.png index cd347af..7ec24cf 100644 Binary files a/change_in_cover.png and b/change_in_cover.png differ diff --git a/electivity.png b/electivity.png index d35fe86..682e696 100644 Binary files a/electivity.png and b/electivity.png differ diff --git a/outputs/change_in_cover.png b/outputs/change_in_cover.png deleted file mode 100644 index cd347af..0000000 Binary files a/outputs/change_in_cover.png and /dev/null differ diff --git a/outputs/cheatgrass_change_effect.png b/outputs/cheatgrass_change_effect.png deleted file mode 100644 index 6daaa00..0000000 Binary files a/outputs/cheatgrass_change_effect.png and /dev/null differ diff --git a/outputs/understory_effects.png b/outputs/understory_effects.png index 5403808..9075e64 100644 Binary files a/outputs/understory_effects.png and b/outputs/understory_effects.png differ diff --git a/plot_level_understory.R b/plot_level_understory.R index f691c4c..1fc9220 100644 --- a/plot_level_understory.R +++ b/plot_level_understory.R @@ -11,8 +11,8 @@ library("plyr") setwd("C:/Users/Sam/Documents/Research/MS Thesis/Understory") plot_data <- read.csv("./clean data/plot_data.csv") +greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv") -# process other files and import #------------------------------------------------------------ # Analysis @@ -22,7 +22,7 @@ plot_data <- read.csv("./clean data/plot_data.csv") #exploration hist(asin(sqrt(plot_data$All/100))) hist(logit(plot_data$All/100 + min(plot_data[plot_data$All != 0, "All"]))) -hist(log(plot_data$All/100 + .01)) +hist(log(plot_data$All/100)) hist(asin(sqrt(plot_data$Cheatgrass/100))) hist(logit(plot_data$Cheatgrass/100 + min(plot_data[plot_data$Cheatgrass != 0, "Cheatgrass"]))) @@ -48,13 +48,32 @@ plot(plot_data$Cheatgrass ~ plot_data$cwd_normal_cum) #--------------------------------------------------------------- #linear mixed effects models for functional groups at plot level #------------------------------------------------------------------- + +## all + +all_plot <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +all_plot <- lmer(logit(All/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC) + (1|Cluster), data = plot_data) + +summary(lm(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + scale(AWC), data = plot_data)) + +summary(all_plot, ddf = "Kenward-Roger") +r.squaredGLMM(all_plot) +plot(allEffects(all_plot, partial.residuals = TRUE)) +plot(all_plot) +AICc(all_plot) +scatter.smooth(residuals(all_plot) ~ predict(all_plot)) + ## Cheatgrass cheatgrass_plot <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC) + (1|Cluster), data = plot_data) -cheatgrass_plot <- lmer(logit(Cheatgrass/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + - scale(AWC) + (1|Cluster), data = plot_data) +# cheatgrass_plot <- lmer(logit(Cheatgrass/100) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + + # scale(AWC) + (1|Cluster), data = plot_data) summary(lm(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC), data = plot_data)) @@ -86,35 +105,40 @@ plot(pgrass_plot) pforb_plot <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC) + (1|Cluster), data = plot_data) -summary(lm(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + - scale(AWC), data = plot_data)) - summary(pforb_plot, ddf = "Kenward-Roger") r.squaredGLMM(pforb_plot) plot(allEffects(pforb_plot, partial.residuals = TRUE)) +plot(inv.as(predict(pforb_plot)) ~ I(plot_data$Pforb/100)) +abline(0,1) plot(pforb_plot) ## Shrubs shrub_plot <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + scale(AWC) + (1|Cluster), data = plot_data) - -summary(lm(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + - scale(AWC), data = plot_data)) +# shrub_plot <- lmer(log(Shrub_cover_li+.01) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) + +# scale(AWC) + (1|Cluster), data = plot_data) summary(shrub_plot, ddf = "Kenward-Roger") r.squaredGLMM(shrub_plot) plot(allEffects(shrub_plot, partial.residuals = TRUE)) +plot(inv.as(predict(shrub_plot)) ~ I(plot_data$Shrub_cover_li)) +abline(0,1) plot(shrub_plot) +### Plot regression coefficients +library(coefplot2) +coefplot2(pgrass_plot) +coefplot2(cheatgrass_plot, add = TRUE, spacing = 10) + + + #----------------------------------------------- # Comparison of greenwood and 2015 #----------------------------------------------- -compare <- join(plot_daub_cover, greenwood_under, by = "Plot", type = "inner") -compare <- join(compare, tcover[, c("Plot", "Tree_cover", "Shrub_cover_li")], by = "Plot", type = "inner") -compare <- join(compare, plot_data[, c("Plot", "Dead_ba", "Delta_pdc", "cwd_normal_cum")], by = "Plot") +compare <- join(plot_data, greenwood_under, by = "Plot", type = "inner") #cheatgrass plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover) @@ -159,7 +183,7 @@ mean(compare$Pgrass-compare$Perrenial.Grass.Cover) t.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") wilcox.test(compare$Pgrass, compare$Perrenial.Grass.Cover, paired = TRUE, alternative="less") -pg_change_lm <- lm(Pgrass - Perrenial.Grass.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +pg_change_lm <- lm(Pgrass - Perrenial.Grass.Cover ~ Delta_pdc*cwd_normal_cum, data= compare) summary(pg_change_lm) ## Perr forb @@ -172,7 +196,7 @@ mean(I(compare$Aforb + compare$Pforb) - compare$Forb.Cover) t.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") wilcox.test(I(compare$Aforb + compare$Pforb), compare$Forb.Cover, paired = TRUE, alternative="less") -pf_change_lm <- lm(Aforb + Pforb - Forb.Cover ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +pf_change_lm <- lm(Aforb + Pforb - Forb.Cover ~ Delta_pdc*cwd_normal_cum, data= compare) summary(pf_change_lm) ## Shrub @@ -187,7 +211,7 @@ t.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="les wilcox.test(compare$Shrub, compare$Shrub.Cover.Total, paired = TRUE, alternative="less") hist(compare$Shrub - I(compare$Shrub.Cover.Total)) -s_change_lm <- lm(Shrub - Shrub.Cover.Total ~ Delta_pdc*cwd_normal_cum, data= compare[-5,]) +s_change_lm <- lm(Shrub - Shrub.Cover.Total ~ Delta_pdc*cwd_normal_cum, data= compare) summary(s_change_lm) ###################### @@ -199,7 +223,7 @@ png(filename="change_in_cover.png", type="cairo", units="in", width = 4, - height=6, + height=4, pointsize=12, res=160) @@ -223,6 +247,22 @@ plot(compare$Pgrass ~ compare$Perrenial.Grass.Cover, xlim = c(0,35), ylim = c(0, axis(side = 2, at = c(0,5,15,25,35), tcl = -.24, labels = FALSE) abline(0,1, lty = 2, col = 'grey30') +par(mar = c(2,2,2,1)) +plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover, xlim = c(0, 20), ylim = c(0, 20), + xlab = "", ylab = "", main = "Cheatgrass", xaxt = 'n', yaxt = 'n', + bg='grey60', + col = 'grey30', + pch = 21, + cex = 1.5, + cex.axis = .8, + cex.main = 1.3) + axis(side = 1, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 1, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + axis(side = 2, at = c(0,5,10,15,20), labels = TRUE) + axis(side = 2, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) + abline(0,1, lty = 2, col = 'grey30') + + par(mar = c(2,2,2,1)) plot(compare$Pforb ~ compare$Forb.Cover, xlim = c(0, 12), ylim = c(0, 12), main = "Perennial Forb", xlab = "", ylab = "", xaxt = 'n', yaxt = 'n', @@ -255,21 +295,7 @@ plot(compare$Shrub ~ compare$Shrub.Cover.Total, xlim = c(0, 30), ylim = c(0,30), abline(0,1, lty = 2, col = 'grey30') -par(mar = c(2,2,2,1)) -plot(compare$Cheatgrass ~ compare$Cheatgrass.Cover, xlim = c(0, 20), ylim = c(0, 20), - xlab = "", ylab = "", main = "Cheatgrass", xaxt = 'n', yaxt = 'n', - bg='grey60', - col = 'grey30', - pch = 21, - cex = 1.5, - cex.axis = .8, - cex.main = 1.3) - axis(side = 1, at = c(0,5,10,15,20), labels = TRUE) - axis(side = 1, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) - axis(side = 2, at = c(0,5,10,15,20), labels = TRUE) - axis(side = 2, at = c(0,2.5,7.5,12.5,17.5), tcl = -.25, labels = FALSE) - abline(0,1, lty = 2, col = 'grey30') - + mtext("2005 Cover (%)", side = 1, outer = TRUE) mtext("2015 Cover (%)", side = 2, outer = TRUE) dev.off() diff --git a/understory_effects_plots.R b/understory_effects_plots.R index 0aa6948..7d8e0c1 100644 --- a/understory_effects_plots.R +++ b/understory_effects_plots.R @@ -11,10 +11,10 @@ pardefault <- par(no.readonly = T) source('addTrans.R', echo=FALSE) ### set some parameters to be use throughout -cg <- cheatgrass_ag -pg <- pgrass_ag -pf <- pforb_ag -sh <- shrub_ag +cg <- cheatgrass_plot +pg <- pgrass_plot +pf <- pforb_plot +sh <- shrub_plot # function to calculate standard errors for a given vcov matrix and row of x values, to get # point estimates of prediction error for a given prediction. See http://www.ats.ucla.edu/stat/r/faq/deltamethod.htm @@ -24,9 +24,15 @@ calc.se <- function(x){ sqrt(as.matrix(x) %*% vcov(model) %*% t(x)) } +inv.as <- function(x){ + (sin(x))^2 +} + ############################################## ## Creating new predictions ############################################## +model <- cg + # Calculate residuals from predicted values and observed values resid <- residuals(model) @@ -34,14 +40,14 @@ resid <- residuals(model) # for unscaling variables -- these things come from a full model in the model set, you'll have to change # the number or the path depending upon what kind of model object and the order of your models, etc. # You can also just get them from the raw data. -tc_mean <- mean(ag_data$Tree_cover) -tc_sd <- sd(ag_data$Tree_cover) +tc_mean <- mean(plot_data$Tree_cover) +tc_sd <- sd(plot_data$Tree_cover) -cwd_mean <- mean(ag_data$cwd_normal_cum) -cwd_sd <- sd(ag_data$cwd_normal_cum) +cwd_mean <- mean(plot_data$cwd_normal_cum) +cwd_sd <- sd(plot_data$cwd_normal_cum) -Delta_pdc_mean <- mean(ag_data$Delta_pdc) -Delta_pdc_sd <- sd(ag_data$Delta_pdc) +Delta_pdc_mean <- mean(plot_data$Delta_pdc) +Delta_pdc_sd <- sd(plot_data$Delta_pdc) ################################ ## Calcuate partial residuals, effects, for avg_BA @@ -50,50 +56,47 @@ Delta_pdc_sd <- sd(ag_data$Delta_pdc) model <- cg C_tc <- data.frame( - Tree_cover = seq(min(ag_data$Tree_cover), max(ag_data$Tree_cover), length.out = 100), - cwd_normal_cum = mean(ag_data$cwd_normal_cum), - Delta_pdc = mean(ag_data$Delta_pdc), - AWC = mean(ag_data$AWC) + Tree_cover = seq(min(plot_data$Tree_cover), max(plot_data$Tree_cover), length.out = 102), + cwd_normal_cum = mean(plot_data$cwd_normal_cum), + Delta_pdc = mean(plot_data$Delta_pdc), + AWC = mean(plot_data$AWC) ) C_awc <- data.frame( - Tree_cover = mean(ag_data$Tree_cover), - cwd_normal_cum = mean(ag_data$cwd_normal_cum), - Delta_pdc = mean(ag_data$Delta_pdc), - AWC = seq(min(ag_data$AWC), max(ag_data$AWC), length.out = 100) + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = mean(plot_data$cwd_normal_cum), + Delta_pdc = mean(plot_data$Delta_pdc), + AWC = seq(min(plot_data$AWC), max(plot_data$AWC), length.out = 102) ) C_pdc_10 <- data.frame( - Tree_cover = mean(ag_data$Tree_cover), - cwd_normal_cum = unname(quantile(ag_data$cwd_normal_cum, .1)), - Delta_pdc = seq(min(ag_data$Delta_pdc), max(ag_data$Delta_pdc), length.out = 100), - AWC = mean(ag_data$AWC) + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .1)), + Delta_pdc = seq(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc), length.out = 102), + AWC = mean(plot_data$AWC) ) C_pdc_50 <- data.frame( - Tree_cover = mean(ag_data$Tree_cover), - cwd_normal_cum = unname(quantile(ag_data$cwd_normal_cum, .5)), - Delta_pdc = seq(min(ag_data$Delta_pdc), max(ag_data$Delta_pdc), length.out = 100), - AWC = mean(ag_data$AWC) + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .5)), + Delta_pdc = seq(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc), length.out = 102), + AWC = mean(plot_data$AWC) ) C_pdc_90 <- data.frame( - Tree_cover = mean(ag_data$Tree_cover), - cwd_normal_cum = unname(quantile(ag_data$cwd_normal_cum, .9)), - Delta_pdc = seq(min(ag_data$Delta_pdc), max(ag_data$Delta_pdc), length.out = 100), - AWC = mean(ag_data$AWC) + Tree_cover = mean(plot_data$Tree_cover), + cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .9)), + Delta_pdc = seq(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc), length.out = 102), + AWC = mean(plot_data$AWC) ) - -model <- cg -C <- C_tc -predictor <- "Tree_cover" +# model <- cg +# C <- C_tc +# predictor <- "Tree_cover" calculate_preds <- function(model, C, predictor){ preds <- data.frame( - predictions = predict(model, se = TRUE, newdata = C, type = "response")$fit, - lo = predict(model, se = TRUE, newdata = C)$fit - 1.96*predict(model, se = TRUE, newdata = C)$se.fit, - hi = predict(model, se = TRUE, newdata = C)$fit + 1.96*predict(model, se = TRUE, newdata = C)$se.fit, + predictions = predict(model, re.form = NA, newdata = C, type = "response"), predictor = C[, predictor] ) @@ -101,7 +104,9 @@ calculate_preds <- function(model, C, predictor){ } -# Use the function to get all the stuff we need (predictions, partial residuals, and SEs) +# Use the function to get all the stuff we need +# (predictions, partial residuals, and SEs) +# this is the dumbest way to do this but I don't want to rewrite it preds_tc_cg <- calculate_preds(cg, C_tc, "Tree_cover") preds_tc_pg <- calculate_preds(pg, C_tc, "Tree_cover") preds_tc_pf <- calculate_preds(pf, C_tc, "Tree_cover") @@ -128,8 +133,9 @@ preds_pdc_90_pf <- calculate_preds(pf, C_pdc_90, "Delta_pdc") preds_pdc_90_sh <- calculate_preds(sh, C_pdc_90, "Delta_pdc") # create plots of log-partial-residuals y-axis versus natural x-scale - -setwd("C:\\Users\\Sam\\Google Drive\\Thesis related\\Understory\\outputs") #where the plots go +# this is also dumb and I could do it in a loop except I don't want to +# put all the predictions in a list +setwd("C:/Users/Sam/Documents/Research/MS Thesis/Understory/outputs") #where the plots go opar <- par(no.readonly = TRUE) par(opar) @@ -149,7 +155,7 @@ par(oma = c(2,4,0,0), mar = c(3,1,1,1), family = "serif", bty = 'n') plot(NA, ylim = c(0,.1), - xlim = c(I(min(ag_data$Tree_cover)*100), I(max(ag_data$Tree_cover)*100)), + xlim = c(I(min(plot_data$Tree_cover)*100), I(max(plot_data$Tree_cover)*100)), xlab = "", ylab = "", bg='grey60', @@ -159,10 +165,10 @@ plot(NA, bty = 'n', xaxt = 'n', yaxt = 'n') -lines(exp(preds_tc_cg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 1, col = "#1b9e77") -lines(exp(preds_tc_pg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 2, col = "#d95f02") -lines(exp(preds_tc_pf$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 3, col = "#7570b3") -lines(exp(preds_tc_sh$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 4, col = "#e7298a") +lines(inv.as(preds_tc_cg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_tc_pg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_tc_pf$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_tc_sh$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 4, col = "#e7298a") axis(side = 1) axis(side = 2) mtext(text = "Tree Cover (%)", side = 1, line = 2.2) @@ -172,7 +178,7 @@ mtext(text = "(a)", side = 1, line = -10, adj = 0.05) plot(NA, ylim = c(0,.1), - xlim = c(min(ag_data$AWC), max(ag_data$AWC)), + xlim = c(min(plot_data$AWC), max(plot_data$AWC)), xlab = "", ylab = "", bg='grey60', @@ -182,10 +188,10 @@ plot(NA, bty = 'n', xaxt = 'n', yaxt = 'n') -lines(exp(preds_awc_cg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") -lines(exp(preds_awc_pg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") -lines(exp(preds_awc_pf$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") -lines(exp(preds_awc_sh$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +lines(inv.as(preds_awc_cg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_awc_pg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_awc_pf$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_awc_sh$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") axis(side = 1) axis(side = 2) mtext(text = "Soil AWC (%)", side = 1, line = 2.2) @@ -197,8 +203,8 @@ legend("topright", legend = c("Cheatgrass", "Perr. Grass", "Perr. Forb", "Shrub" lty = c(1,2,3,4), lwd = 2, cex = 1.3, col = c("#1b9e77", "#d95f02", "#7570b3", "#e7298a")) plot(NA, - ylim = c(0,.1), - xlim = c(min(ag_data$Delta_pdc), max(ag_data$Delta_pdc)), + ylim = c(0,.2), + xlim = c(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc)), xlab = "", ylab = "", bg='grey60', @@ -208,18 +214,18 @@ plot(NA, bty = 'n', xaxt = 'n', yaxt = 'n') -lines(exp(preds_pdc_10_cg$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") -lines(exp(preds_pdc_10_pg$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") -lines(exp(preds_pdc_10_pf$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") -lines(exp(preds_pdc_10_sh$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +lines(inv.as(preds_pdc_10_cg$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_pdc_10_pg$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_pdc_10_pf$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_pdc_10_sh$predictions) ~ I(preds_pdc_10_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") axis(side = 1) axis(side = 2) text(x = -30, y = .09, labels = "10% CWD", cex = 1.3) mtext(text = "(c)", side = 1, line = -10, adj = 0.05) plot(NA, - ylim = c(0,.1), - xlim = c(min(ag_data$Delta_pdc), max(ag_data$Delta_pdc)), + ylim = c(0,.2), + xlim = c(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc)), xlab = "", ylab = "", bg='grey60', @@ -229,10 +235,10 @@ plot(NA, bty = 'n', xaxt = 'n', yaxt = 'n') -lines(exp(preds_pdc_50_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") -lines(exp(preds_pdc_50_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") -lines(exp(preds_pdc_50_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") -lines(exp(preds_pdc_50_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +lines(inv.as(preds_pdc_50_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_pdc_50_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_pdc_50_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_pdc_50_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") axis(side = 1) axis(side = 2) text(x = -30, y = .09, labels = "50% CWD", cex = 1.3) @@ -241,8 +247,8 @@ mtext(text = "(d)", side = 1, line = -10, adj = 0.05) plot(NA, - ylim = c(0,.1), - xlim = c(min(ag_data$Delta_pdc), max(ag_data$Delta_pdc)), + ylim = c(0,.3), + xlim = c(min(plot_data$Delta_pdc), max(plot_data$Delta_pdc)), xlab = "", ylab = "", bg='grey60', @@ -252,10 +258,10 @@ plot(NA, bty = 'n', xaxt = 'n', yaxt = 'n') -lines(exp(preds_pdc_90_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") -lines(exp(preds_pdc_90_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") -lines(exp(preds_pdc_90_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") -lines(exp(preds_pdc_90_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") +lines(inv.as(preds_pdc_90_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77") +lines(inv.as(preds_pdc_90_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#d95f02") +lines(inv.as(preds_pdc_90_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 3, col = "#7570b3") +lines(inv.as(preds_pdc_90_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 4, col = "#e7298a") axis(side = 1) axis(side = 2) text(x = -30, y = .09, labels = "90% CWD", cex = 1.3) @@ -350,9 +356,9 @@ par(mar = c(0,1,0,0), oma = c(3,4,2,2), family = "serif", bty = 'n') pch = 21, cex.lab = 1.5) abline(h = 0, lty = 4, lwd = 1.3) - lines((exp(cg_change_10$predictions) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 2, lty = 1, col = "#1b9e77") - lines((exp(cg_change_10$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") - lines((exp(cg_change_10$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + lines((inv.as(cg_change_10$predictions) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 2, lty = 1, col = "#1b9e77") + lines((inv.as(cg_change_10$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") + lines((inv.as(cg_change_10$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)) ~ I(cg_change_10$predictor), lwd = 1, lty = 1, col = "#1b9e77") polygon(c(cg_change_10$predictor, rev(cg_change_10$predictor)), c((exp(cg_change_10$hi) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)), rev((exp(cg_change_10$lo) + min(compare$Cheatgrass/100 - compare$Cheatgrass.Cover/100)))), col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals diff --git a/understory_electivity.R b/understory_electivity.R index 5fc890e..d7aa228 100644 --- a/understory_electivity.R +++ b/understory_electivity.R @@ -142,11 +142,11 @@ plots_to_use <- unique(ms$Plot)[plots_to_use_dead & plots_to_use_live] # Calculate test and randomized distributions #------------------------------------------------------------------------------ #perform analysis for the four FTs -types <- c("All", "Perennial grass", "Cheatgrass", "Perennial forb", "Shrub") - +types <- c("Perennial grass", "Cheatgrass", "Perennial forb", "Shrub") +ntypes <- length(types) #initialize a list of dataframes to catch the electivities at each plot/FT elect_results <- list() -for(i in 1:5){ +for(i in 1:ntypes){ elect_results[[i]] <- data.frame(Plot = numeric(0), Dead = numeric(0), Live = numeric(0), @@ -206,7 +206,7 @@ nit <- 999 #number of monte carlo draws #initialize a list of dataframes to catch the electivities at each plot/FT elect_means <- list() -for(i in 1:5){ +for(i in 1:ntypes){ elect_means[[i]] <- data.frame(Iter = numeric(0), Dead = numeric(0), Live = numeric(0), @@ -279,7 +279,7 @@ for(type in types){ results_boots <- list() -for(i in 1:5){ +for(i in 1:ntypes){ results_boots[[i]] <- data.frame(Dead = numeric(2), Live = numeric(2), Inter = numeric(2), @@ -291,7 +291,7 @@ for(i in 1:5){ names(results_boots) <- types #calculate p-values from distribution of bootstrapped mean electivities -for(i in 1:5){ +for(i in 1:ntypes){ emp_mean <- apply(elect_results[[i]][, 2:7], 2, FUN = function(x){mean(x, na.rm = TRUE)}) pvals <- c(sum(elect_means[[i]][2] > emp_mean[1]), @@ -319,15 +319,15 @@ png(filename="electivity.png", type="cairo", units="in", width = 4, - height=6, - pointsize=12, + height=4, + pointsize=10, res=160) -par(mfrow = c(3,2), +par(mfrow = c(2,2), mar = c(2,1,1,2), oma = c(2,3,1,0)) -for(i in 1:5){ +for(i in 1:ntypes){ melt_elect <- melt(elect_results[[i]][, 2:4]) plot(melt_elect$value ~ I(as.numeric(melt_elect$variable)+ runif(nrow(melt_elect), -.1, .1 )), @@ -359,8 +359,8 @@ for(i in 1:5){ } text(x = c(1,2,3), y = 0.9, labels = pvals[2:4]) - mtext(text = types[i], outer = FALSE, side = 3, line = 0) - mtext(text = paste0("(", letters[i], ")"), outer = FALSE, side = 3, at = 0.5, line = 0) + mtext(text = types[i], outer = FALSE, side = 3, line = 0.3) + mtext(text = paste0("(", letters[i], ")"), outer = FALSE, side = 3, at = 0.5, line = 0.3) mtext(text = "Electivity", outer = TRUE, side = 2, line = 1.5)