Update to comments

Tidied up code and wrote better comments

.Rhistory

@@ -1,102 +1,3 @@
-axis(side = 1)
-axis(side = 2)
-text(x = -30, y = .18, 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,.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_90_cg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_pdc_90_cg$predictor[c(1,102)]), inv.as(preds_pdc_90_cg$predictions[c(1,102)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-lines(inv.as(preds_pdc_90_pg$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 1, col = "#000000")
-points(I(preds_pdc_90_cg$predictor[c(1,102)]), inv.as(preds_pdc_90_pg$predictions[c(1,102)]), pch = 21, col = "#000000", bg = "#000000")
-lines(inv.as(preds_pdc_90_pf$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#E69F00")
-points(I(preds_pdc_90_cg$predictor[c(1,102)]), inv.as(preds_pdc_90_pf$predictions[c(1,102)]), pch = 23, col = "#E69F00", bg = "#E69F00")
-lines(inv.as(preds_pdc_90_sh$predictions) ~ I(preds_pdc_90_cg$predictor), lwd = 2, lty = 2, col = "#56B4E9")
-points(I(preds_pdc_90_cg$predictor[c(1,102)]), inv.as(preds_pdc_90_sh$predictions[c(1,102)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
-axis(side = 1)
-axis(side = 2)
-text(x = -30, y = .18, labels = "90% CWD", cex = 1.3)
-mtext(text = "(e)", side = 1, line = -10, adj = 0.05)
-dev.off()
-tiff(filename="understory_effects_test.tiff",
-type="cairo",
-units="in",
-width = 7,
-height=5,
-pointsize=15,
-res=600)
-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,.12),
-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 = 2, col = "#1b9e77")
-points(I(preds_tc_cg$predictor[c(1,102)]*100), inv.as(preds_tc_cg$predictions[c(1,102)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-lines(inv.as(preds_tc_pg$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 2, col = "#000000")
-points(I(preds_tc_cg$predictor[c(1,102)]*100), inv.as(preds_tc_pg$predictions[c(1,102)]), pch = 22, col = "#000000", bg = "#000000")
-lines(inv.as(preds_tc_pf$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 1, col = "#E69F00")
-points(I(preds_tc_cg$predictor[c(1,102)]*100), inv.as(preds_tc_pf$predictions[c(1,102)]), pch = 23, col = "#E69F00", bg = "#E69F00")
-lines(inv.as(preds_tc_sh$predictions) ~ I(preds_tc_cg$predictor*100), lwd = 2, lty = 1, col = "#56B4E9")
-points(I(preds_tc_cg$predictor[c(1,102)]*100), inv.as(preds_tc_sh$predictions[c(1,102)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
-axis(side = 1)
-axis(side = 2)
-mtext(text = "Tree cover (%)", side = 1, line = 2.2)
-mtext(text = "Understory 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 = 2, col = "#1b9e77")
-points(I(preds_awc_cg$predictor[c(1,102)]), inv.as(preds_awc_cg$predictions[c(1,102)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-lines(inv.as(preds_awc_pg$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 2, col = "#000000")
-points(I(preds_awc_cg$predictor[c(1,102)]), inv.as(preds_awc_pg$predictions[c(1,102)]), pch = 21, col = "#000000", bg = "#000000")
-lines(inv.as(preds_awc_pf$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 2, col = "#E69F00")
-points(I(preds_awc_cg$predictor[c(1,102)]), inv.as(preds_awc_pf$predictions[c(1,102)]), pch = 23, col = "#E69F00", bg = "#E69F00")
-lines(inv.as(preds_awc_sh$predictions) ~ I(preds_awc_cg$predictor), lwd = 2, lty = 2, col = "#56B4E9")
-points(I(preds_awc_cg$predictor[c(1,102)]), inv.as(preds_awc_sh$predictions[c(1,102)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
-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", "Per. Grass", "Per. Forb", "Shrub"),
-lty = c(1), lwd = 2, cex = 1.3, col = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"))
-plot(NA,
-ylim = c(0,.25),
-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',
@@ -510,3 +411,102 @@ axis(side = 2)
 text(x = -30, y = .18, labels = "90% CWD", cex = 1.3)
 mtext(text = "(e)", side = 1, line = -10, adj = 0.05)
 dev.off()
+## 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(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)
+tiff(filename="./outputs/electivity.tiff",
+type="cairo",
+units="in",
+width = 4,
+height=4,
+pointsize=10,
+res=600)
+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(NA,
+ylim = c(-1,1),
+xlim = c(.7, 3.3),
+xaxt = "n")
+pvals <- results_boots[[i]][2, ]
+abline(h = 0)
+vioplot(elect_means[[i]][[2]], add = TRUE, col = addTrans("blue", 30), drawRect = FALSE)
+vioplot(elect_means[[i]][[3]], at = 2, add = TRUE, col = addTrans("blue", 30), drawRect = FALSE)
+vioplot(elect_means[[i]][[4]], at = 3, add = TRUE, col = addTrans("blue", 30), drawRect = FALSE)
+points(melt_elect$value ~ I(as.numeric(melt_elect$variable)+ runif(nrow(melt_elect), -.1, .1 )),
+pch = 21,
+bg = "grey")
+means <- aggregate(melt_elect$value, by = list(melt_elect$variable), FUN = function(x){mean(x, na.rm = TRUE)})
+segments(x0 = c(0.85, 1.85, 2.85), y0 = means$x, x1 = c(1.15, 2.15, 3.15), 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", "Inter"))
+}
+# text(x = c(1,2,3), y = 0.9, labels = pvals[4:6])
+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 = "Ivlev's E", outer = TRUE, side = 2, line = 1.5)
+}
+dev.off()
+library(plyr)
+library(vioplot)
+library(reshape2)
+library(multcompView)
+tiff(filename="./outputs/electivity.tiff",
+type="cairo",
+units="in",
+width = 4,
+height=4,
+pointsize=10,
+res=600)
+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(NA,
+ylim = c(-1,1),
+xlim = c(.7, 3.3),
+xaxt = "n")
+pvals <- results_boots[[i]][2, ]
+abline(h = 0)
+vioplot(elect_means[[i]][[2]], add = TRUE, col = addTrans("blue", 30), drawRect = FALSE)
+vioplot(elect_means[[i]][[3]], at = 2, add = TRUE, col = addTrans("blue", 30), drawRect = FALSE)
+vioplot(elect_means[[i]][[4]], at = 3, add = TRUE, col = addTrans("blue", 30), drawRect = FALSE)
+points(melt_elect$value ~ I(as.numeric(melt_elect$variable)+ runif(nrow(melt_elect), -.1, .1 )),
+pch = 21,
+bg = "grey")
+means <- aggregate(melt_elect$value, by = list(melt_elect$variable), FUN = function(x){mean(x, na.rm = TRUE)})
+segments(x0 = c(0.85, 1.85, 2.85), y0 = means$x, x1 = c(1.15, 2.15, 3.15), 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", "Inter"))
+}
+# text(x = c(1,2,3), y = 0.9, labels = pvals[4:6])
+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 = "Ivlev's E", outer = TRUE, side = 2, line = 1.5)
+}
+dev.off()