diff --git a/.RData b/.RData
index 0c42d36..0785a4d 100644
Binary files a/.RData and b/.RData differ
diff --git a/.Rhistory b/.Rhistory
index bd388fb..4448144 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,12 +1,41 @@
-mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
-mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
+mtext(text = "(d)", side = 1, line = -10, adj = 0.05)
+plot(NA,
+ylim = c(0,.25),
+xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
+xlab = "",
+ylab = "",
+bg='grey60',
+col = 'grey30',
+pch = 21,
+cex.lab = 1.5,
+bty = 'n',
+xaxt = 'n',
+yaxt = 'n')
+lines(inv.as(preds_tc_90_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+lines(inv.as(preds_tc_90_pg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#000000")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_pg$predictions[c(1,200)]), pch = 22, col = "#000000", bg = "#000000")
+lines(inv.as(preds_tc_90_pf$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 2, col = "#E69F00")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_pf$predictions[c(1,200)]), pch = 23, col = "#E69F00", bg = "#E69F00")
+lines(inv.as(preds_tc_90_sh$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 2, col = "#56B4E9")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_sh$predictions[c(1,200)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
+axis(side = 1)
+axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+text(x = -30, y = .18, labels = "90% CWD", cex = 1.3)
+mtext(text = "(e)", side = 1, line = -10, adj = 0.05)
 dev.off()
+par(pardefault)
+###################################################################################
+# Cheatgrass change over time
+# Figure 2
+###################################################################################
 ## new data
 model <- cg_change_lm
 model2 <- cg_change_lm_no_outliers
 # model <- model2
 compare$delta_cg <- compare$Cheatgrass - compare$Cheatgrass.Cover
-# compare2 <- compare[-c(5, 10), ]
+compare2 <- compare[-c(5, 10), ]
+compare_outliers <- compare[c(5,10), ]
 #create dataframe of new data to calculate predictions and envelope
 C_cg_change_10 <- data.frame(
 cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .15)),
@@ -63,102 +92,6 @@ yaxt = "n",
 bg='grey60',
 col = 'grey30',
 pch = 21,
-cex.lab = 1.5)
-abline(h = 0, lty = 4, lwd = 1.3)
-lines(cg_change[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 4, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$lo ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$hi ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-polygon(c(cg_change_no_out[[i]]$predictor, rev(cg_change_no_out[[i]]$predictor)), c(cg_change_no_out[[i]]$hi,
-rev(cg_change_no_out[[i]]$lo)),
-col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals
-# I also hate this: plot the residuals ~ delta_tc for each level of CWD
-points(
-if(i == 1){
-# part_resids[[i]][model$model[, "cwd_normal_cum"] < unname(quantile(compare$cwd_normal_cum, 0.3))] ~
-#    model$model[model$model[, "cwd_normal_cum"] < unname(quantile(compare$cwd_normal_cum, 0.3)), "Delta_tc"]
-compare[compare$cwd_normal_cum < unname(quantile(compare$cwd_normal_cum, 0.3)),]$delta_cg/100 ~
-compare[compare$cwd_normal_cum < unname(quantile(compare$cwd_normal_cum, 0.3)),]$Delta_tc
-}else if(i ==2){
-compare[findInterval(compare$cwd_normal_cum, unname(quantile(compare$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~
-compare[findInterval(compare$cwd_normal_cum, unname(quantile(compare$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc
-#   part_resids[[i]][findInterval(model$model[, "cwd_normal_cum"], unname(quantile(compare$cwd_normal_cum, c(0.3, 0.7)))) == 1] ~
-# model$model[findInterval(model$model[, "cwd_normal_cum"], unname(quantile(compare$cwd_normal_cum, c(0.3, 0.7)))) == 1, "Delta_tc"]
-}else if (i == 3){
-compare[compare$cwd_normal_cum > unname(quantile(compare$cwd_normal_cum, 0.7)),]$delta_cg/100 ~
-compare[compare$cwd_normal_cum > unname(quantile(compare$cwd_normal_cum, 0.7)),]$Delta_tc
-# part_resids[[i]][model$model[, "cwd_normal_cum"] > unname(quantile(compare$cwd_normal_cum, 0.7))] ~
-#   model$model[model$model[, "cwd_normal_cum"] > unname(quantile(compare$cwd_normal_cum, 0.7)), "Delta_tc"]
-},
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-if(i == 1){
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-# axis(side = 2, at = c(-0.1, 0, 0.1), labels = c(-10, 0, 10))
-}
-axis(side = 1, at = c(-40, -20, 0), labels = TRUE)
-axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25)
-text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
-}
-mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
-mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
-dev.off()
-compare_outliers <- compare[c(5,10), ]
-compare2 <- compare[-c(5, 10), ]
-compare_outliers[1]
-compare_outliers
-compare_outliers[1]
-compare_outliers
-compare_outliers[1, ]
-points(compare_outliers[1, ]$delta_cg,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$cwd_normal_cum,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-compare_outliers[1, ]$delta_cg
-compare_outliers[1, ]$cwd_normal_cum
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-compare_outliers[1, ]$delta_cg
-compare_outliers[1, ]$Delta_tc
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
-bg='white',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-#generate the figure
-tiff(filename="./outputs/Figure_2_cheatgrass_change_effect.tif",
-type="cairo",
-units="in",
-width = 4,
-height=3,
-pointsize=12,
-compression = "lzw",
-res=600)
-layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE))
-par(mar = c(0,1,0,0), oma = c(6,4,2,2), bty = 'n')
-panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
-for(i in c(1:3)){
-plot(NA,
-ylim = c(-0.06, 0.16),
-xlim = c(-70, 10 ),
-xlab = "",
-ylab = "",
-xaxt = "n",
-yaxt = "n",
-bg='grey60',
-col = 'grey30',
-pch = 21,
 cex = 0.8,
 cex.lab = 1.5)
 abline(h = 0, lty = 4, lwd = 1.3)
@@ -176,26 +109,31 @@ compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)
 bg='grey60',
 col = 'grey30',
 pch = 21,
-cex = 1.5)
+cex = 1.2)
 }else if(i ==2){
 points(compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~
 compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc,
 bg='grey60',
 col = 'grey30',
 pch = 21,
-cex = 1.5)
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
+cex = 1.2)
+points(compare_outliers[1, ]$delta_cg/100 ~ compare_outliers[1, ]$Delta_tc,
 bg='white',
 col = 'grey30',
 pch = 21,
-cex = 1.5)
+cex = 1.2)
 }else if (i == 3){
 points(compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$delta_cg/100 ~
 compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$Delta_tc,
 bg='grey60',
 col = 'grey30',
 pch = 21,
-cex = 1.5)
+cex = 1.2)
+points(compare_outliers[2, ]$delta_cg/100 ~ compare_outliers[2, ]$Delta_tc,
+bg='white',
+col = 'grey30',
+pch = 21,
+cex = 1.2)
 }
 if(i == 1){
 axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
@@ -208,305 +146,367 @@ text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
 mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
 mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
 dev.off()
-#generate the figure
-tiff(filename="./outputs/Figure_2_cheatgrass_change_effect.tif",
+#------------------------------------------------------------------------------
+# Effects plots with partial residuals for supplement
+#------------------------------------------------------------------------------
+library("effects")
+# Perennial grasses
+tiff(filename="./outputs/Supplementary figure effects cgrass.tiff",
 type="cairo",
 units="in",
-width = 4,
-height=3,
-pointsize=12,
+width = 7,
+height=5,
+pointsize=15,
 compression = "lzw",
 res=600)
-layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE))
-par(mar = c(0,1,0,0), oma = c(6,4,2,2), bty = 'n')
+#graphical params
+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), bty = 'n')
 panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
-for(i in c(1:3)){
+k <- 1
+#which  model to use?
+model <- pgrass_plot_noscale
+let <- letters
+vars <- c("Tree_cover", "AWC", "Delta_tc", "Delta_tc", "Delta_tc")
+x_labels <- c("Tree cover (%)", "Soil AWC (%)", "Change in live canopy (%)")
+for(i in 1:2){
+var <- vars[i]
+eff <- Effect(model, partial.residuals = TRUE, focal.predictors = var)
+line_data <- data.frame(y = inv.as(eff$fit),
+lower = inv.as(eff$lower),
+upper = inv.as(eff$upper),
+x = eff[["x"]][[var]]
+)
+x_new <- eff[["x"]][[var]]
+x_orig <- eff[["x.all"]][[var]]
+y_trans <- eff$fit
+fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+# that we can find the closest fitted value and add the residual for that
+p_resids <- inv.as(eff$residuals + fitted)
+# eff[["data"]][[var]]
+points_data <- data.frame(x = x_orig,
+y = p_resids)
 plot(NA,
-ylim = c(-0.06, 0.16),
-xlim = c(-70, 10 ),
+ylim = c(min(c(y, points_data$y)), max(c(y, points_data$y))),
+xlim = c(min(x_new), max(x_new)),
 xlab = "",
 ylab = "",
-xaxt = "n",
-yaxt = "n",
 bg='grey60',
 col = 'grey30',
 pch = 21,
-cex = 0.8,
-cex.lab = 1.5)
-abline(h = 0, lty = 4, lwd = 1.3)
-lines(cg_change[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 4, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$lo ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$hi ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-polygon(c(cg_change_no_out[[i]]$predictor, rev(cg_change_no_out[[i]]$predictor)), c(cg_change_no_out[[i]]$hi,
-rev(cg_change_no_out[[i]]$lo)),
-col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals
-# I also hate this: plot the residuals ~ delta_tc for each level of CWD
-if(i == 1){
-points(compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}else if(i ==2){
-points(compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~
-compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
-bg='white',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}else if (i == 3){
-points(compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}
-if(i == 1){
+cex.lab = 1.5,
+bty = 'n',
+xaxt = 'n',
+yaxt = 'n')
+lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+lines(line_data$upper ~ line_data$x)
+lines(line_data$lower ~ line_data$x)
+polygon(c(line_data$x, rev(line_data$x)),
+c(line_data$upper, rev(line_data$lower)),
+col = addTrans("#68EBC4",30), border = NA)
+points(points_data$y ~ points_data$x,
+pch = 21, col = "#1b9e77", bg = "#1b9e77")
+axis(side = 1)
 axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-# axis(side = 2, at = c(-0.1, 0, 0.1), labels = c(-10, 0, 10))
-}
-axis(side = 1, at = c(-40, -20, 0), labels = TRUE)
-axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25)
-text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
+mtext(text = var, 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)
 }
-mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
-mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
-dev.off()
-layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE))
-par(mar = c(0,1,0,0), oma = c(6,4,2,2), bty = 'n')
-panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
-for(i in c(1:3)){
+plot(NA)
+for(i in 1:3){
+eff <- Effect(model, partial.residuals = TRUE, focal.predictors = c("Delta_tc", "cwd_normal_cum"),
+xlevels = list(Delta_tc = 100),
+quantiles = c(0.1, 0.5, 0.9))
+cwd_levels <- eff$variables$cwd_normal_cum$levels
+line_data <- data.frame(y = inv.as(eff$fit[1:99 + (i-1) * 100]),
+lower = inv.as(eff$lower[1:99 + (i-1) * 100]),
+upper = inv.as(eff$upper[1:99 + (i-1) * 100]),
+x = eff[["x"]][["Delta_tc"]][1:99 + (i-1) * 100]
+)
+x_new <- eff[["x"]][["Delta_tc"]]
+x_orig_all <- eff[["x.all"]]
+which_x <- x_orig_all$cwd_normal_cum - cwd_levels[i] < 0.001
+x_orig <- x_orig_all[which_x, ][["Delta_tc"]]
+y_trans <- eff$fit[1:99 + (i-1) * 100]
+fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+# that we can find the closest fitted value and add the residual for that
+p_resids <- inv.as(eff$residuals[which_x] + fitted)
+points_data <- data.frame(x = x_orig,
+y = p_resids)
 plot(NA,
-ylim = c(-0.06, 0.16),
-xlim = c(-70, 10 ),
+ylim = c(0, max(c(y, points_data$y), na.rm = TRUE)),
+xlim = c(min(x_new), max(x_new)),
 xlab = "",
 ylab = "",
-xaxt = "n",
-yaxt = "n",
 bg='grey60',
 col = 'grey30',
 pch = 21,
-cex = 0.8,
-cex.lab = 1.5)
-abline(h = 0, lty = 4, lwd = 1.3)
-lines(cg_change[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 4, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$lo ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$hi ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-polygon(c(cg_change_no_out[[i]]$predictor, rev(cg_change_no_out[[i]]$predictor)), c(cg_change_no_out[[i]]$hi,
-rev(cg_change_no_out[[i]]$lo)),
-col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals
-# I also hate this: plot the residuals ~ delta_tc for each level of CWD
-if(i == 1){
-points(compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}else if(i ==2){
-points(compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~
-compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
-bg='white',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}else if (i == 3){
-points(compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}
-if(i == 1){
+cex.lab = 1.5,
+bty = 'n',
+xaxt = 'n',
+yaxt = 'n')
+lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+lines(line_data$upper ~ line_data$x)
+lines(line_data$lower ~ line_data$x)
+polygon(c(line_data$x, rev(line_data$x)),
+c(line_data$upper, rev(line_data$lower)),
+col = addTrans("#68EBC4",30), border = NA)
+points(points_data$y ~ points_data$x,
+pch = 21, col = "#1b9e77", bg = "#1b9e77")
+axis(side = 1)
 axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-# axis(side = 2, at = c(-0.1, 0, 0.1), labels = c(-10, 0, 10))
+mtext(text = var, 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)
 }
-axis(side = 1, at = c(-40, -20, 0), labels = TRUE)
-axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25)
-text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
+dev.off()
+shrub_plot
+summary(shrub_plot)
+all_plot_tc <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), data = plot_data)
+scatter.smooth(residuals(all_plot_tc) ~ predict(all_plot_tc))
+cheatgrass_plot_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc) * scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), data = plot_data)
+summary(cheatgrass_plot_tc, ddf = "Kenward-Roger")
+plot(allEffects(cheatgrass_plot_noscale_tc, partial.residuals = TRUE))
+cheatgrass_plot_noscale_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+AWC + (1|Cluster), data = plot_data)
+plot(allEffects(cheatgrass_plot_noscale_tc, partial.residuals = TRUE))
+cheatgrass_plot_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc) * scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), data = plot_data)
+summary(cheatgrass_plot_tc, ddf = "Kenward-Roger")
+AICc(cheatgrass_plot_tc)
+AICc(cheatgrass_plot)
+library("MuMIn")
+library("lme4")
+library("car")
+library("effects")
+library("plyr")
+library("lmerTest")
+library("vegan")
+library("reshape2")
+library("plyr")
+# import data created by data prep script
+plot_data <- read.csv("./clean data/plot_data.csv")
+plot_data$Delta_tc <- (plot_data$Tree_cover*100 - plot_data$Total_cover)
+plot_data$Delta_ba <- plot_data$Live_ba - plot_data$Live_ba_2005
+plot_data <- plot_data[plot_data$Cluster != "NPELECTRICEEL", ]
+# understory variables from Greenwood and Weisberg 2008
+greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv")
+#------------------------------------------------------------
+# Analysis
+#------------------------------------------------------------
+#---------------------------------------------------------------
+#linear mixed effects models for functional groups at plot level
+#-------------------------------------------------------------------
+## all understory vegetation
+all_plot <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), 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)
+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))
+## 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(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(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(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)
+scatter.smooth(residuals(pgrass_plot_tc) ~ predict(pgrass_plot_tc))
+pgrass_plot_tc <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), data = plot_data)
+scatter.smooth(residuals(pgrass_plot_tc) ~ predict(pgrass_plot_tc))
+summary(pgrass_plot_tc, ddf = "Kenward-Roger")
+plot(inv.as(predict(pforb_plot_tc)) ~ I(plot_data$Pforb/100))
+Q
+pforb_plot_tc <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), data = plot_data)
+pforb_plot_noscale_tc <- lmer(asin(sqrt(Pforb/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+AWC + (1|Cluster), data = plot_data)
+summary(pforb_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pforb_plot_tc)
+plot(inv.as(predict(pforb_plot_tc)) ~ I(plot_data$Pforb/100))
+abline(0,1)
+scatter.smooth(residuals(pforb_plot_tc) ~ predict(pforb_plot_tc))
+## Shrubs
+shrub_plot_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+scale(AWC) + (1|Cluster), data = plot_data)
+summary(shrub_plot_tc, ddf = "Kenward-Roger")
+plot(plot_data$Delta_pdc ~ plot_data$Delta_tc)
+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
+# These models were generated in plot_level_understory.R, and that
+# script needs to be run first and these models saved in the global environment.
+# Some other global variables are also needed, like the raw plot_data data frame.
+# Bad practices but I didn't think it was worth refactoring this thing.
+all <- all_plot_noscale_tc
+cg <- cheatgrass_plot_noscale
+pg <- pgrass_plot_noscale
+pf <- pforb_plot_noscale
+sh <- shrub_plot_noscale
+# 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))
 }
-mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
-mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
-bg='white',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-points(compare_outliers[1, ]$delta_cg ~ compare_outliers[1, ]$Delta_tc,
-bg='black',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-compare_outliers[1, ]$delta_cg
-compare_outliers[1, ]$Delta_tc
-delta_cg
-compare_outliers[1, ]$delta_cg/100
-#generate the figure
-tiff(filename="./outputs/Figure_2_cheatgrass_change_effect.tif",
+inv.as <- function(x){
+(sin(x))^2
+}
+# function to find the nearest y-value for an x data point, to add to partial residuals
+# for plotting. Stolen from effects package
+closest <- function(x, x0){
+apply(outer(x, x0, FUN = function(x, x0) abs(x - x0)), 1, which.min)
+}
+#####################################################################
+# All understory vegetation, just two panels
+# Figure 3
+#####################################################################
+vars <- c("Tree_cover", "Delta_tc")
+labels <- c("Tree cover", "Change in tree cover (%)")
+tiff(filename="./outputs/Figure_3_all_understory_effects.tiff",
 type="cairo",
 units="in",
-width = 4,
-height=3,
+width = 3,
+height=5,
 pointsize=12,
 compression = "lzw",
 res=600)
-layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE))
-par(mar = c(0,1,0,0), oma = c(6,4,2,2), bty = 'n')
-panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
-for(i in c(1:3)){
+par(mfrow = c(2,1), oma = c(2,3,0,0), mar = c(3,1,1,1), bty = 'n')
+for(i in c(1:2)){
+var <- vars[i]
+eff <- Effect(all_plot, partial.residuals = TRUE, focal.predictors = var)
+y <- eff$fit
+x <- eff$x[[var]]
+x.fit <- eff$x.all[[var]]
+fitted <- y[closest(x.fit, x)]
+resids <- inv.as(eff$residuals + fitted)
 plot(NA,
-ylim = c(-0.06, 0.16),
-xlim = c(-70, 10 ),
+xlim = c(min(x), max(x)),
+ylim = c(0, 50),
 xlab = "",
-ylab = "",
-xaxt = "n",
-yaxt = "n",
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 0.8,
-cex.lab = 1.5)
-abline(h = 0, lty = 4, lwd = 1.3)
-lines(cg_change[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 4, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$lo ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$hi ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-polygon(c(cg_change_no_out[[i]]$predictor, rev(cg_change_no_out[[i]]$predictor)), c(cg_change_no_out[[i]]$hi,
-rev(cg_change_no_out[[i]]$lo)),
-col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals
-# I also hate this: plot the residuals ~ delta_tc for each level of CWD
-if(i == 1){
-points(compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-}else if(i ==2){
-points(compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~
-compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-points(compare_outliers[1, ]$delta_cg/100 ~ compare_outliers[1, ]$Delta_tc,
-bg='black',
-col = 'grey30',
+ylab = "")
+points(resids*100 ~ x.fit,
 pch = 21,
-cex = 1.5)
-}else if (i == 3){
-points(compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$Delta_tc,
 bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
-points(compare_outliers[2, ]$delta_cg/100 ~ compare_outliers[2, ]$Delta_tc,
-bg='black',
-col = 'grey30',
-pch = 21,
-cex = 1.5)
+col = 'grey30')
+lines(inv.as(y)*100 ~ x, lwd = 2)
+lines(inv.as(eff$lower)*100 ~ x)
+lines(inv.as(eff$upper)*100 ~ x)
+polygon(c(x, rev(x)), c(inv.as(eff$upper)*100, rev(inv.as(eff$lower))*100),
+col = addTrans("#68EBC4",30), border = NA)
+mtext(side = 2, text = "Understory cover (%)", outer = TRUE, line = 1.7)
+mtext(side = 1, text = labels[i], line = 2)
 }
-if(i == 1){
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-# axis(side = 2, at = c(-0.1, 0, 0.1), labels = c(-10, 0, 10))
+dev.off()
+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
+# These models were generated in plot_level_understory.R, and that
+# script needs to be run first and these models saved in the global environment.
+# Some other global variables are also needed, like the raw plot_data data frame.
+# Bad practices but I didn't think it was worth refactoring this thing.
+all <- all_plot_noscale_tc
+cg <- cheatgrass_plot_noscale
+pg <- pgrass_plot_noscale
+pf <- pforb_plot_noscale
+sh <- shrub_plot_noscale
+# 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))
 }
-axis(side = 1, at = c(-40, -20, 0), labels = TRUE)
-axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25)
-text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
+inv.as <- function(x){
+(sin(x))^2
 }
-mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
-mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
-dev.off()
-#generate the figure
-tiff(filename="./outputs/Figure_2_cheatgrass_change_effect.tif",
+# function to find the nearest y-value for an x data point, to add to partial residuals
+# for plotting. Stolen from effects package
+closest <- function(x, x0){
+apply(outer(x, x0, FUN = function(x, x0) abs(x - x0)), 1, which.min)
+}
+#####################################################################
+# All understory vegetation, just two panels
+# Figure 3
+#####################################################################
+vars <- c("Tree_cover", "Delta_tc")
+labels <- c("Tree cover", "Change in tree cover (%)")
+tiff(filename="./outputs/Figure_3_all_understory_effects.tiff",
 type="cairo",
 units="in",
-width = 4,
-height=3,
+width = 3,
+height=5,
 pointsize=12,
 compression = "lzw",
 res=600)
-layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE))
-par(mar = c(0,1,0,0), oma = c(6,4,2,2), bty = 'n')
-panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
-for(i in c(1:3)){
+par(mfrow = c(2,1), oma = c(2,3,0,0), mar = c(3,1,1,1), bty = 'n')
+for(i in c(1:2)){
+var <- vars[i]
+eff <- Effect(all_plot_tc, partial.residuals = TRUE, focal.predictors = var)
+y <- eff$fit
+x <- eff$x[[var]]
+x.fit <- eff$x.all[[var]]
+fitted <- y[closest(x.fit, x)]
+resids <- inv.as(eff$residuals + fitted)
 plot(NA,
-ylim = c(-0.06, 0.16),
-xlim = c(-70, 10 ),
+xlim = c(min(x), max(x)),
+ylim = c(0, 50),
 xlab = "",
-ylab = "",
-xaxt = "n",
-yaxt = "n",
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 0.8,
-cex.lab = 1.5)
-abline(h = 0, lty = 4, lwd = 1.3)
-lines(cg_change[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 4, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$lo ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-lines(cg_change_no_out[[i]]$hi ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
-polygon(c(cg_change_no_out[[i]]$predictor, rev(cg_change_no_out[[i]]$predictor)), c(cg_change_no_out[[i]]$hi,
-rev(cg_change_no_out[[i]]$lo)),
-col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals
-# I also hate this: plot the residuals ~ delta_tc for each level of CWD
-if(i == 1){
-points(compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.2)
-}else if(i ==2){
-points(compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~
-compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc,
-bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.2)
-points(compare_outliers[1, ]$delta_cg/100 ~ compare_outliers[1, ]$Delta_tc,
-bg='white',
-col = 'grey30',
+ylab = "")
+points(resids*100 ~ x.fit,
 pch = 21,
-cex = 1.2)
-}else if (i == 3){
-points(compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$delta_cg/100 ~
-compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$Delta_tc,
 bg='grey60',
-col = 'grey30',
-pch = 21,
-cex = 1.2)
-points(compare_outliers[2, ]$delta_cg/100 ~ compare_outliers[2, ]$Delta_tc,
-bg='white',
-col = 'grey30',
-pch = 21,
-cex = 1.2)
+col = 'grey30')
+lines(inv.as(y)*100 ~ x, lwd = 2)
+lines(inv.as(eff$lower)*100 ~ x)
+lines(inv.as(eff$upper)*100 ~ x)
+polygon(c(x, rev(x)), c(inv.as(eff$upper)*100, rev(inv.as(eff$lower))*100),
+col = addTrans("#68EBC4",30), border = NA)
+mtext(side = 2, text = "Understory cover (%)", outer = TRUE, line = 1.7)
+mtext(side = 1, text = labels[i], line = 2)
 }
-if(i == 1){
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-# axis(side = 2, at = c(-0.1, 0, 0.1), labels = c(-10, 0, 10))
-}
-axis(side = 1, at = c(-40, -20, 0), labels = TRUE)
-axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25)
-text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
-}
-mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
-mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
 dev.off()
+plot(plot_data$Delta_pdc ~ plot_data$Dead_ba)
+plot(plot_data$Delta_tc ~ plot_data$Dead_ba)
+summary(lm(plot_data$Delta_tc ~ plot_data$Dead_ba))
+summary(lm(plot_data$Delta_pdc ~ plot_data$Dead_ba))
+summary(lm(plot_data$Delta_pdc ~ plot_data$Dead_rat))
+summary(lm(plot_data$Delta_tc ~ plot_data$Dead_rat))
diff --git a/.Rproj.user/EAEB7255/console06/INDEX001 b/.Rproj.user/EAEB7255/console06/INDEX001
index 22fbe44..32144ba 100644
--- a/.Rproj.user/EAEB7255/console06/INDEX001
+++ b/.Rproj.user/EAEB7255/console06/INDEX001
@@ -1 +1 @@
-[{"allow_restart":true,"alt_buffer":false,"autoclose":1,"buffered_output":"\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n","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
+[{"allow_restart":true,"alt_buffer":false,"autoclose":1,"buffered_output":"\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n","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 f53e558..70829f6 100644
--- a/.Rproj.user/EAEB7255/pcs/source-pane.pper
+++ b/.Rproj.user/EAEB7255/pcs/source-pane.pper
@@ -1,3 +1,3 @@
 {
-    "activeTab" : 4
+    "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 2291b13..0be2e78 100644
--- a/.Rproj.user/EAEB7255/pcs/windowlayoutstate.pper
+++ b/.Rproj.user/EAEB7255/pcs/windowlayoutstate.pper
@@ -1,14 +1,14 @@
 {
     "left" : {
-        "panelheight" : 779,
-        "splitterpos" : 324,
+        "panelheight" : 763,
+        "splitterpos" : 317,
         "topwindowstate" : "NORMAL",
-        "windowheight" : 817
+        "windowheight" : 801
     },
     "right" : {
-        "panelheight" : 779,
-        "splitterpos" : 487,
+        "panelheight" : 763,
+        "splitterpos" : 477,
         "topwindowstate" : "NORMAL",
-        "windowheight" : 817
+        "windowheight" : 801
     }
 }
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/prop/2A13C8B0 b/.Rproj.user/EAEB7255/sources/prop/2A13C8B0
index 01533ac..9c06b90 100644
--- a/.Rproj.user/EAEB7255/sources/prop/2A13C8B0
+++ b/.Rproj.user/EAEB7255/sources/prop/2A13C8B0
@@ -1,4 +1,4 @@
 {
-    "cursorPosition" : "44,0",
-    "scrollLine" : "32"
+    "cursorPosition" : "10,14",
+    "scrollLine" : "0"
 }
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC b/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC
index 90055be..44610cd 100644
--- a/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC
+++ b/.Rproj.user/EAEB7255/sources/prop/8DF0DCEC
@@ -1,4 +1,4 @@
 {
-    "cursorPosition" : "74,42",
-    "scrollLine" : "53"
+    "cursorPosition" : "3,46",
+    "scrollLine" : "0"
 }
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6 b/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6
index 9098667..7e12f07 100644
--- a/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6
+++ b/.Rproj.user/EAEB7255/sources/prop/A1FFFDE6
@@ -1,4 +1,4 @@
 {
-    "cursorPosition" : "97,3",
+    "cursorPosition" : "106,47",
     "scrollLine" : "0"
 }
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/prop/A30485D9 b/.Rproj.user/EAEB7255/sources/prop/A30485D9
new file mode 100644
index 0000000..44610cd
--- /dev/null
+++ b/.Rproj.user/EAEB7255/sources/prop/A30485D9
@@ -0,0 +1,4 @@
+{
+    "cursorPosition" : "3,46",
+    "scrollLine" : "0"
+}
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/prop/AA116D9C b/.Rproj.user/EAEB7255/sources/prop/AA116D9C
new file mode 100644
index 0000000..5b28f81
--- /dev/null
+++ b/.Rproj.user/EAEB7255/sources/prop/AA116D9C
@@ -0,0 +1,4 @@
+{
+    "cursorPosition" : "239,15",
+    "scrollLine" : "256"
+}
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/prop/D455C530 b/.Rproj.user/EAEB7255/sources/prop/D455C530
index 6b4e8eb..bb53532 100644
--- a/.Rproj.user/EAEB7255/sources/prop/D455C530
+++ b/.Rproj.user/EAEB7255/sources/prop/D455C530
@@ -1,5 +1,5 @@
 {
-    "cursorPosition" : "17,54",
+    "cursorPosition" : "13,32",
     "scrollLine" : "0",
     "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 a632a06..3912725 100644
--- a/.Rproj.user/EAEB7255/sources/prop/INDEX
+++ b/.Rproj.user/EAEB7255/sources/prop/INDEX
@@ -18,7 +18,9 @@ C%3A%2FUsers%2FSam%2FGoogle%20Drive%2FThesis%20related%2FUnderstory%2Funderstory
 ~%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%2Fplot_level_understory_tc.R="AA116D9C"
 ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_data_prep.R="A1FFFDE6"
 ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_effects_plots.R="8DF0DCEC"
+~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_effects_plots_tc.R="A30485D9"
 ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_electivity.R="D455C530"
 ~%2FResearch%2FMS%20Thesis%2FUnderstory%2Funderstory_electivity010819.R="8CD3B36B"
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/9E5C163F b/.Rproj.user/EAEB7255/sources/s-197F81BE/9E5C163F
similarity index 62%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/9E5C163F
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/9E5C163F
index 2087c3a..8b1b4ab 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/9E5C163F
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/9E5C163F
@@ -2,18 +2,18 @@
     "collab_server" : "",
     "contents" : "",
     "created" : 1547225685816.000,
-    "dirty" : true,
+    "dirty" : false,
     "encoding" : "UTF-8",
     "folds" : "88|48|165|0|\n",
     "hash" : "0",
     "id" : "9E5C163F",
-    "lastKnownWriteTime" : 1593721425,
-    "last_content_update" : 1593723732579,
-    "path" : "~/Research/MS Thesis/Understory/plot_level_understory.R",
-    "project_path" : "plot_level_understory.R",
+    "lastKnownWriteTime" : 1596053285,
+    "last_content_update" : 1596053285119,
+    "path" : "~/Research/MS Thesis/Understory/plot_level_understory_tc.R",
+    "project_path" : "plot_level_understory_tc.R",
     "properties" : {
-        "cursorPosition" : "44,0",
-        "scrollLine" : "32"
+        "cursorPosition" : "239,15",
+        "scrollLine" : "256"
     },
     "relative_order" : 1,
     "source_on_save" : false,
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/9E5C163F-contents b/.Rproj.user/EAEB7255/sources/s-197F81BE/9E5C163F-contents
similarity index 84%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/9E5C163F-contents
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/9E5C163F-contents
index 4fb1df5..d1cd503 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/9E5C163F-contents
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/9E5C163F-contents
@@ -18,7 +18,7 @@ library("plyr")
 
 # import data created by data prep script
 plot_data <- read.csv("./clean data/plot_data.csv")
-plot_data$Delta_tc <- plot_data$Tree_cover - plot_data$Total_cover
+plot_data$Delta_tc <- (plot_data$Tree_cover*100 - plot_data$Total_cover)
 plot_data$Delta_ba <- plot_data$Live_ba - plot_data$Live_ba_2005
 plot_data <- plot_data[plot_data$Cluster != "NPELECTRICEEL", ]
 
@@ -267,54 +267,78 @@ par(opar)
 
 ## all understory vegetation
 
-all_plot <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+all_plot_tc <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
                    scale(AWC) + (1|Cluster), 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))
+
+
+all_plot_noscale_tc <- lmer(asin(sqrt(All/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                   AWC + (1|Cluster), data = plot_data)
+
+summary(all_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(all_plot_tc)
+plot(allEffects(all_plot_tc, partial.residuals = TRUE))
+plot(all_plot_tc)
+AICc(all_plot_tc)
+scatter.smooth(residuals(all_plot_tc) ~ predict(all_plot_tc))
 
 ## Cheatgrass
 
-cheatgrass_plot <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+cheatgrass_plot_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc) * scale(cwd_normal_cum) +
                           scale(AWC) + (1|Cluster), 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))
+
+cheatgrass_plot_noscale_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                          AWC + (1|Cluster), data = plot_data)
+
+summary(cheatgrass_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(cheatgrass_plot_tc)
+plot(allEffects(cheatgrass_plot_noscale_tc, partial.residuals = TRUE))
+plot(cheatgrass_plot_tc)
+AICc(cheatgrass_plot_tc)
+scatter.smooth(residuals(cheatgrass_plot_tc) ~ predict(cheatgrass_plot_tc))
 
 ## Perr grass
 
-pgrass_plot <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+pgrass_plot_tc <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
                       scale(AWC) + (1|Cluster), data = plot_data)
-summary(pgrass_plot, ddf = "Kenward-Roger")
-r.squaredGLMM(pgrass_plot)
-plot(allEffects(pgrass_plot, partial.residuals = TRUE))
-plot(pgrass_plot)
+
+pgrass_plot_noscale_tc <- lmer(asin(sqrt(Pgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                      AWC + (1|Cluster), data = plot_data)
+
+summary(pgrass_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pgrass_plot_tc)
+plot(allEffects(pgrass_plot_tc, partial.residuals = TRUE))
+plot(pgrass_plot_tc)
+scatter.smooth(residuals(pgrass_plot_tc) ~ predict(pgrass_plot_tc))
+
 
 ## Perr forbs
 
-pforb_plot <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+pforb_plot_tc <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
                      scale(AWC) + (1|Cluster), 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))
+pforb_plot_noscale_tc <- lmer(asin(sqrt(Pforb/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                     AWC + (1|Cluster), data = plot_data)
+
+summary(pforb_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pforb_plot_tc)
+plot(allEffects(pforb_plot_tc, partial.residuals = TRUE))
+plot(inv.as(predict(pforb_plot_tc)) ~ I(plot_data$Pforb/100))
 abline(0,1)
-plot(pforb_plot)
+plot(pforb_plot_tc)
+scatter.smooth(residuals(pforb_plot_tc) ~ predict(pforb_plot_tc))
 
 
 ## Shrubs
-shrub_plot <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+shrub_plot_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_tc)*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))
+
+shrub_plot_noscale_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                     AWC + (1|Cluster), data = plot_data)
+
+
+summary(shrub_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(shrub_plot_tc)
+plot(allEffects(shrub_plot_tc, partial.residuals = TRUE))
+plot(inv.as(predict(shrub_plot_tc)) ~ I(plot_data$Shrub_cover_li))
 abline(0,1)
 
 #-----------------------------------------------
@@ -340,7 +364,7 @@ compare <- compare[-c(5, 10), ] #outlier to remove?
 
 cg_change_lm <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ Delta_tc*cwd_normal_cum, data= compare)
 cg_change_lm_no_outliers <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ 
-                                 Delta_tc*cwd_normal_cum, data= compare[-c(5,10), ])
+                                 Delta_tc + cwd_normal_cum, data= compare[-c(5,10), ])
 
 summary(cg_change_lm)
 plot(cg_change_lm)
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/BA4339 b/.Rproj.user/EAEB7255/sources/s-197F81BE/BA4339
similarity index 81%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/BA4339
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/BA4339
index 36faffd..125a586 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/BA4339
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/BA4339
@@ -7,12 +7,12 @@
     "folds" : "",
     "hash" : "0",
     "id" : "BA4339",
-    "lastKnownWriteTime" : 1579016181,
-    "last_content_update" : 1579016181197,
+    "lastKnownWriteTime" : 1596127268,
+    "last_content_update" : 1596127268334,
     "path" : "~/Research/MS Thesis/Understory/understory_electivity.R",
     "project_path" : "understory_electivity.R",
     "properties" : {
-        "cursorPosition" : "17,54",
+        "cursorPosition" : "13,32",
         "scrollLine" : "0",
         "tempName" : "Untitled1"
     },
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/BA4339-contents b/.Rproj.user/EAEB7255/sources/s-197F81BE/BA4339-contents
similarity index 66%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/BA4339-contents
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/BA4339-contents
index 27f7be7..0d52263 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/BA4339-contents
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/BA4339-contents
@@ -11,6 +11,7 @@ library(plyr)
 library(vioplot)
 library(reshape2)
 library(multcompView)
+source('addTrans.R', echo=FALSE)
 
 set.seed(16091315)
 
@@ -36,16 +37,6 @@ species$unique_quad <- paste0(species$Plot, species$Transect, species$Meter)
 
 count(species, vars = c("Plot", "Transect"))
 
-
-#little test to look at species relationships with CWD
-# species_aggregate <- aggregate(species[c("Cover")], by = list(species$Plot, species$Species), FUN = mean)
-# names(species_aggregate)[c(1,2)] <- c("Plot", "Species")
-# species_aggregate <- join(species_aggregate, plot_data[c("Plot", "cwd_normal_cum")], by = c("Plot"))
-# 
-# plot(log(Cover) ~ cwd_normal_cum, data = species_aggregate[species_aggregate$Species == "ARTTRIV", ],
-#      xlim = c(200, 450))
-# 
-
 #-----------------------------------------------------------------------
 # Calculate mean cover, cover by quadrat, cover by plot
 #-------------------------------------------------------------------
@@ -78,40 +69,21 @@ 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)
+
+#reduce the number of microsites
 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)))#))
+
+#standardize transect direction names
 ms[ms$Transect == "e", "Transect"] <- "E"
 ms[ms$Transect == "s", "Transect"] <- "S"
 ms[ms$Transect == "w", "Transect"] <- "W"
@@ -119,6 +91,7 @@ ms[ms$Transect == "n", "Transect"] <- "N"
 
 ms$Plot <- as.character(ms$Plot)
 
+#remove the "bonus" quadrats
 ms <- ms[(ms$Transect %in% c("N", "S", "E", "W")), ]
 
 ms[ms$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120"
@@ -126,44 +99,27 @@ ms[ms$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120"
 
 ms$unique_quad <- paste0(ms$Plot, ms$Transect, ms$Meter)
 
+
+
+#------------------------------------------------------------------------------
+# Electivity calculations
+#------------------------------------------------------------------------------
 #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)
+## select which plots have enough microsites of each type
+# for which the cover is enough to have 1% cover in each microsite type,
+# if the cover were uniformly distributed
+types <- c("All", "Perennial grass", "Cheatgrass", "Perennial forb", "Shrub")
+ntypes <- length(types)
 
-plots_to_use <- unique(ms$Plot)[plots_to_use_dead & plots_to_use_live]
+plots_to_use <- list(NA, NA, NA, NA, NA)
+names(plots_to_use) <- types
 
+n_plots <- 1
+type <- "Cheatgrass"  
 
-#------------------------------------------------------------------------------
-# 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){
@@ -175,51 +131,66 @@ for(i in 1:ntypes){
                                    DI = numeric(0),
                                    LI = numeric(0))
 }
+
 names(elect_results) <- types
 
 
-for(type in types){ # calculate electivity for each functional type
+for(type in types){
   
-  for (i in 1:length(plots_to_use)){
-    
+  for(j in 1:length(unique(ms$Plot))){
     if(type == "All"){
-      daub_cov <- daub[daub$Plot == as.character(plots_to_use[i]), ]
+      daub_cov <- daub[daub$Plot == unique(ms$Plot)[j], ]
       daub_cov <- daub_cov[daub_cov$Cover.type %in% c("Cheatgrass", "Other Ann Grass", 
-                        "Perennial grass", "Annual forb", "Perennial forb", "Shrub"), ]
+                                                      "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 <- daub[daub$Cover.type == type & daub$Plot == unique(ms$Plot)[j], ]
       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)
+    total_cov <- sum(daub_cov$Midpoint.value)
+    uniform_cov <- total_cov/20 #divy up cover evenly among quadrats
+    
+    table_ms <- table(daub_cov$ms)
+    
+    table_ms_uniform <- table_ms*uniform_cov #multiple average quadrat cover by microsite 
+                                             # amount to generate "expected" cover
     
-    cov <- aggregate(daub_cov[, "prop_cov"], by = list(daub_cov$ms), FUN = sum)
+    use_plot <- sum(table_ms_uniform >= 1) == 3 #are there greater than or equal to 1% cover
+                                               #in each microsite type?
     
-    prev <- table(daub_cov$ms)/20
+    plots_to_use[[type]][j] <- use_plot #save plot condition to global list
+    
+    if(use_plot){
+    
+      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 = unique(ms$Plot)[j],
+                          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 = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Live",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Live",2])
+      elect$DI = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Inter",2])
+      elect$LI = (cov[cov$Group.1 == "Live",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Live",2] + cov[cov$Group.1 == "Inter",2])
+      
+      elect_results[[type]] <- rbind(elect_results[[type]], elect) #I know this is slow
+    }
     
-    #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 = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Live",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Live",2])
-    elect$DI = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Inter",2])
-    elect$LI = (cov[cov$Group.1 == "Live",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Live",2] + cov[cov$Group.1 == "Inter",2])
     
-    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
@@ -237,6 +208,7 @@ for(i in 1:ntypes){
                                  DI = numeric(0),
                                  LI = numeric(0))
 }
+
 names(elect_means) <- types
 
 system.time(#takes about a half hour
@@ -246,40 +218,46 @@ system.time(#takes about a half hour
 for(type in types){
   
   plot_to_use_type <- subset(elect_results[[type]], !is.na(Dead))$Plot
+  nplots <- length(plot_to_use_type)
   
   for(j in 1:nit){
     
     #temporary storage for each iteration
-    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){
+    elect_rand <- data.frame(Plot = character(nplots),
+                             Dead = numeric(nplots),
+                             Live = numeric(nplots),
+                             Inter = numeric(nplots),
+                             DL = numeric(nplots),
+                             DI = numeric(nplots),
+                             LI = numeric(nplots), stringsAsFactors = FALSE)
+ 
+    for (k in 1:nplots){
+      # calculate electivity for a each plot using randomized cover values
+      plot_select <- plot_to_use_type[k]
+      
       if(type == "All"){
-        daub_cov <- daub[daub$Plot == as.character(plot_use), ]
+        daub_cov <- daub[daub$Plot == as.character(plot_select), ]
         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 <- daub[daub$Cover.type == type & daub$Plot == as.character(plot_select), ]
         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)
+      
+      #randomize the cover values by quadrat
       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
       
       #temporary df to store results from each plot
-      elect <- data.frame(Plot = plot_use,
+      elect <- data.frame(Plot = plot_select,
                           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")]) / 
@@ -290,14 +268,13 @@ for(type in types){
       elect$DI = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Inter",2])
       elect$LI = (cov[cov$Group.1 == "Live",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Live",2] + cov[cov$Group.1 == "Inter",2])
       
-      elect_rand <- rbind(elect_rand, elect) #add data from plot to temporary df
+      elect_rand[k, ] <- elect #add data from plot to temporary df
     }
     
     # add mean from electivity from this iteration to global df
     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
-    
   
   }
   
@@ -336,7 +313,7 @@ results_boots[[i]][2, ] <- pvals
 }
 
 #save and read results depending on what you need
-# saveRDS(results_boots, paste0("./outputs/results_boots_", n_plots, "dead.rds"))
+saveRDS(results_boots, paste0("./outputs/results_boots_", n_plots, "dead.rds"))
 
 # results_boots <- readRDS("./outputs/results_boots_1dead.rds")
 
@@ -344,8 +321,8 @@ results_boots[[i]][2, ] <- pvals
 #-----------------------------------------------------------------------------------------------
 # Summary table
 #-----------------------------------------------------------------------------------------------
-summary <- as.data.frame(matrix(nrow = 12, ncol = 8))
-for(i in 1:4){
+summary <- as.data.frame(matrix(nrow = 15, ncol = 8))
+for(i in 1:5){
   for(j in 1:6){
     summary[3*i - 2, j+2] <- results_boots[[i]][[j]][[1]]
     summary[3*i - 1, j+2] <- results_boots[[i]][[j]][[1]] - mean(elect_means[[i]][[j + 1]])
@@ -353,8 +330,8 @@ for(i in 1:4){
   }
 }
 
-summary[, 2] <- rep(c("Empirical", "Difference", "p-val"), times = 4)
-summary[, 1] <- rep(c("Perennial grass", "Cheatgrass", "Perennial forb", "Shrub"), each = 3)
+summary[, 2] <- rep(c("Empirical", "Difference", "p-val"), times = 5)
+summary[, 1] <- rep(types, each = 3)
 names(summary) <- c("FT", "var", "Dead", "Live", "Inter", "D-L", "D-I", "L-I")
 write.csv(summary, "./outputs/electivity summary table.csv")
 
@@ -363,7 +340,7 @@ write.csv(summary, "./outputs/electivity summary table.csv")
 # Figure 5
 #-----------------------------------------------------------------------------------------------
 
-tiff(filename="./outputs/Figure_5_electivity.tiff", 
+tiff(filename="./outputs/Figure_5_electivity_test.tiff", 
     type="cairo",
     units="in", 
     width = 4, 
@@ -376,7 +353,7 @@ par(mfrow = c(2,2),
     mar = c(2,1,1,2),
     oma = c(2,3,1,0))
 
-for(i in 1:ntypes){
+for(i in 2:5){
 
   melt_elect <- melt(elect_results[[i]][, 2:4])
   plot(NA,
@@ -393,22 +370,24 @@ for(i in 1:ntypes){
   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)
 
+  text(x = 3, y = 0.7, labels = paste0("n = ", nrow(melt_elect)/3))
   
-  if(i %in% c(1,2)){
+  if(i %in% c(2,3)){
     axis(1, at = c(1,2,3), labels = FALSE)
   }
   
-  if(i %in% c(3,4)){
+  if(i %in% c(4,5)){
   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 = paste0("(", letters[i-1], ")"), outer = FALSE, side = 3, at = 0.5, line = 0.3)
   mtext(text = "Ivlev's E", outer = TRUE, side = 2, line = 1.5)
   
 
diff --git a/.Rproj.user/EAEB7255/sources/s-197F81BE/C3F69618 b/.Rproj.user/EAEB7255/sources/s-197F81BE/C3F69618
new file mode 100644
index 0000000..f9cb123
--- /dev/null
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/C3F69618
@@ -0,0 +1,22 @@
+{
+    "collab_server" : "",
+    "contents" : "",
+    "created" : 1594816908657.000,
+    "dirty" : false,
+    "encoding" : "UTF-8",
+    "folds" : "",
+    "hash" : "3364106652",
+    "id" : "C3F69618",
+    "lastKnownWriteTime" : 1596053275,
+    "last_content_update" : 1596053275339,
+    "path" : "~/Research/MS Thesis/Understory/understory_effects_plots_tc.R",
+    "project_path" : "understory_effects_plots_tc.R",
+    "properties" : {
+        "cursorPosition" : "3,46",
+        "scrollLine" : "0"
+    },
+    "relative_order" : 4,
+    "source_on_save" : false,
+    "source_window" : "",
+    "type" : "r_source"
+}
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/40A653B0-contents b/.Rproj.user/EAEB7255/sources/s-197F81BE/C3F69618-contents
similarity index 80%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/40A653B0-contents
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/C3F69618-contents
index 374b55f..7e3bd73 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/40A653B0-contents
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/C3F69618-contents
@@ -20,11 +20,11 @@ source('addTrans.R', echo=FALSE)
 # Some other global variables are also needed, like the raw plot_data data frame.
 # Bad practices but I didn't think it was worth refactoring this thing. 
 
-all <- all_plot
-cg <- cheatgrass_plot
-pg <- pgrass_plot
-pf <- pforb_plot
-sh <- shrub_plot
+all <- all_plot_noscale_tc
+cg <- cheatgrass_plot_noscale
+pg <- pgrass_plot_noscale
+pf <- pforb_plot_noscale
+sh <- shrub_plot_noscale
 
 # 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
@@ -49,7 +49,7 @@ closest <- function(x, x0){
 # Figure 3
 #####################################################################
 vars <- c("Tree_cover", "Delta_tc")
-labels <- c("Tree cover", "Change in live canopy (%)")
+labels <- c("Tree cover", "Change in tree cover (%)")
 
 tiff(filename="./outputs/Figure_3_all_understory_effects.tiff", 
      type="cairo",
@@ -65,7 +65,7 @@ par(mfrow = c(2,1), oma = c(2,3,0,0), mar = c(3,1,1,1), bty = 'n')
 for(i in c(1:2)){
   var <- vars[i]
   
-  eff <- Effect(all_plot, partial.residuals = TRUE, focal.predictors = var)
+  eff <- Effect(all_plot_tc, partial.residuals = TRUE, focal.predictors = var)
   
   y <- eff$fit
   x <- eff$x[[var]]
@@ -80,7 +80,7 @@ for(i in c(1:2)){
        xlab = "",
        ylab = "")
   
-  points(resids*100 ~ eff$x.all[[var]], 
+  points(resids*100 ~ x.fit, 
          pch = 21, 
          bg='grey60',
          col = 'grey30')
@@ -169,7 +169,7 @@ calculate_preds <- function(model, C, predictor){
 
 # 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
+# this is the dumbest way to do this but I don't want to refactor 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")
@@ -195,10 +195,6 @@ preds_tc_90_pg <- calculate_preds(pg, C_tc_90, "Delta_tc")
 preds_tc_90_pf <- calculate_preds(pf, C_tc_90, "Delta_tc")
 preds_tc_90_sh <- calculate_preds(sh, C_tc_90, "Delta_tc")
 
-# 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
-
 opar <- par(no.readonly = TRUE)
 par(opar)
 
@@ -516,7 +512,7 @@ dev.off()
 library("effects")
 # Perennial grasses
 
-tiff(filename="./outputs/Supplementary figure effects pgrass.tiff", 
+tiff(filename="./outputs/Supplementary figure effects cgrass.tiff", 
      type="cairo",
      units="in", 
      width = 7, 
@@ -525,135 +521,141 @@ tiff(filename="./outputs/Supplementary figure effects pgrass.tiff",
      compression = "lzw",
      res=600)
 
-model <- cg
-
-eff <- Effect(focal.predictors = c("Tree_cover"), 
-              mod = cg, xlevels = 100, partial.residuals = TRUE)
-
-
-eff2 <- Effect(focal.predictors = c("PC1"), 
-               mod = mod3, xlevels = 100, transformation = list(link = log, inverse = exp))
-
-dat <- data.frame(y = exp(eff$fit),
-                  lower = exp(eff$lower),
-                  upper = exp(eff$upper),
-                  soil = eff[["x"]][["PC1"]])
-
+#graphical params
 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), bty = 'n')
 
 
-plot(NA,
-     ylim = c(0,.02),
-     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,200)]*100), inv.as(preds_tc_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-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,.02),
-     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,200)]), inv.as(preds_awc_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-mtext(text = "Soil AWC (%)", side = 1, line = 2.2)
-mtext(text = "(b)", side = 1, line = -10, adj = 0.05)
+panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
+k <- 1
 
+#which  model to use?
+model <- pgrass_plot_noscale
 
-plot.new()
-legend("topright", legend = c("Cheatgrass", "Per. Grass", "Per. Forb", "Shrub"), 
-       lty = 1, pch = c(21,22,23,24), lwd = 2, cex = 1.3, col = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"),
-       pt.bg = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"))
+let <- letters
 
-plot(NA,
-     ylim = c(0,.02),
-     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
-     xlab = "",
-     ylab = "",
-     bg='grey60',
-     col = 'grey30',
-     pch = 21,
-     cex.lab = 1.5,
-     bty = 'n',
-     xaxt = 'n',
-     yaxt = 'n')
-lines(inv.as(preds_tc_10_cg$predictions) ~ I(preds_tc_10_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_tc_10_cg$predictor[c(1,200)]), inv.as(preds_tc_10_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+vars <- c("Tree_cover", "AWC", "Delta_tc", "Delta_tc", "Delta_tc")
+x_labels <- c("Tree cover (%)", "Soil AWC (%)", "Change in live canopy (%)")
 
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-text(x = -30, y = .018, labels = "10% CWD", cex = 1.3)
-mtext(text = "(c)", side = 1, line = -10, adj = 0.05)
+for(i in 1:2){
 
-plot(NA,
-     ylim = c(0,.02),
-     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
-     xlab = "",
-     ylab = "",
-     bg='grey60',
-     col = 'grey30',
-     pch = 21,
-     cex.lab = 1.5,
-     bty = 'n',
-     xaxt = 'n',
-     yaxt = 'n')
-lines(inv.as(preds_tc_50_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_50_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  var <- vars[i]
+  
+  eff <- Effect(model, partial.residuals = TRUE, focal.predictors = var)
+  
+  line_data <- data.frame(y = inv.as(eff$fit),
+                    lower = inv.as(eff$lower),
+                    upper = inv.as(eff$upper),
+                    x = eff[["x"]][[var]]
+                    )
+  
+  x_new <- eff[["x"]][[var]]
+  x_orig <- eff[["x.all"]][[var]]
+  
+  y_trans <- eff$fit
+  
+  fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+  # that we can find the closest fitted value and add the residual for that 
+  
+  p_resids <- inv.as(eff$residuals + fitted)
+  # eff[["data"]][[var]]
+  
+  points_data <- data.frame(x = x_orig,
+                            y = p_resids)
+  
+  plot(NA,
+       ylim = c(min(c(y, points_data$y)), max(c(y, points_data$y))),
+       xlim = c(min(x_new), max(x_new)),
+       xlab = "",
+       ylab = "",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex.lab = 1.5,
+       bty = 'n',
+       xaxt = 'n',
+       yaxt = 'n')
+  
+  lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(line_data$upper ~ line_data$x)
+  lines(line_data$lower ~ line_data$x)
+  polygon(c(line_data$x, rev(line_data$x)), 
+          c(line_data$upper, rev(line_data$lower)),
+          col = addTrans("#68EBC4",30), border = NA)
+  
+  points(points_data$y ~ points_data$x,
+         pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  
+  axis(side = 1)
+  axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+  mtext(text = var, 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)
 
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-text(x = -30, y = .018, 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)
 
-plot(NA,
-     ylim = c(0,.02),
-     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
-     xlab = "",
-     ylab = "",
-     bg='grey60',
-     col = 'grey30',
-     pch = 21,
-     cex.lab = 1.5,
-     bty = 'n',
-     xaxt = 'n',
-     yaxt = 'n')
-lines(inv.as(preds_tc_90_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-text(x = -30, y = .018, labels = "90% CWD", cex = 1.3)
-mtext(text = "(e)", side = 1, line = -10, adj = 0.05)
+for(i in 1:3){
+  
+  eff <- Effect(model, partial.residuals = TRUE, focal.predictors = c("Delta_tc", "cwd_normal_cum"),
+                xlevels = list(Delta_tc = 100),
+                quantiles = c(0.1, 0.5, 0.9))
+  
+  cwd_levels <- eff$variables$cwd_normal_cum$levels
+  
+  line_data <- data.frame(y = inv.as(eff$fit[1:99 + (i-1) * 100]),
+                          lower = inv.as(eff$lower[1:99 + (i-1) * 100]),
+                          upper = inv.as(eff$upper[1:99 + (i-1) * 100]),
+                          x = eff[["x"]][["Delta_tc"]][1:99 + (i-1) * 100]
+  )
+  
+  x_new <- eff[["x"]][["Delta_tc"]]
+  x_orig_all <- eff[["x.all"]]
+  which_x <- x_orig_all$cwd_normal_cum - cwd_levels[i] < 0.001
+  x_orig <- x_orig_all[which_x, ][["Delta_tc"]]
+  
+  y_trans <- eff$fit[1:99 + (i-1) * 100]
+  
+  fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+  # that we can find the closest fitted value and add the residual for that 
+  
+  p_resids <- inv.as(eff$residuals[which_x] + fitted)
+  
+  points_data <- data.frame(x = x_orig,
+                            y = p_resids)
+  
+  plot(NA,
+       ylim = c(0, max(c(y, points_data$y), na.rm = TRUE)),
+       xlim = c(min(x_new), max(x_new)),
+       xlab = "",
+       ylab = "",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex.lab = 1.5,
+       bty = 'n',
+       xaxt = 'n',
+       yaxt = 'n')
+  
+  lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(line_data$upper ~ line_data$x)
+  lines(line_data$lower ~ line_data$x)
+  polygon(c(line_data$x, rev(line_data$x)), 
+          c(line_data$upper, rev(line_data$lower)),
+          col = addTrans("#68EBC4",30), border = NA)
+  
+  points(points_data$y ~ points_data$x,
+         pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  
+  axis(side = 1)
+  axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+  mtext(text = var, 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)
 
+}
 
 dev.off()
+
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/F6A93CBA b/.Rproj.user/EAEB7255/sources/s-197F81BE/F6A93CBA
similarity index 79%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/F6A93CBA
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/F6A93CBA
index 63d30e1..93eceb2 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/F6A93CBA
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/F6A93CBA
@@ -7,12 +7,12 @@
     "folds" : "",
     "hash" : "0",
     "id" : "F6A93CBA",
-    "lastKnownWriteTime" : 1593691764,
-    "last_content_update" : 1593691764080,
+    "lastKnownWriteTime" : 1594759880,
+    "last_content_update" : 1594759880143,
     "path" : "~/Research/MS Thesis/Understory/understory_data_prep.R",
     "project_path" : "understory_data_prep.R",
     "properties" : {
-        "cursorPosition" : "97,3",
+        "cursorPosition" : "106,47",
         "scrollLine" : "0"
     },
     "relative_order" : 4,
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/F6A93CBA-contents b/.Rproj.user/EAEB7255/sources/s-197F81BE/F6A93CBA-contents
similarity index 95%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/F6A93CBA-contents
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/F6A93CBA-contents
index 29c89e9..0407209 100644
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/F6A93CBA-contents
+++ b/.Rproj.user/EAEB7255/sources/s-197F81BE/F6A93CBA-contents
@@ -28,6 +28,9 @@ 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")
 
+# individual tree data from Greenwood and Weisberg 2008
+greenwood_trees <- read.csv("./raw data/Individual_TreeData_AllData_SF_edits.csv")
+
 # 2005 understory data, from Greenwood and Weisberg, 2008, 
 # Density-dependent tree mortality in pinyon-juniper woodlands, FEM
 greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv")
@@ -39,8 +42,11 @@ tcover <- read.csv("./raw data/tree_and_shrub_cover_020815.csv")
 # 2005 tree cover and other variables from Greenwood and Weisberg 2008
 greenwood_tree_cover <- read.csv("./raw data/Structure_vs_environmental_all_data_SF_edits.csv")
   
-# individual tree data from Greenwood and Weisberg 2008
-greenwood_trees <- read.csv("./raw data/Individual_TreeData_AllData_SF_edits.csv")
+# calculate basal area of ingrowths
+# in_ba <- aggregate(trees[trees$Code == "P", "BA_m2_per_ha"], 
+#           by = list(trees[trees$Code == "P", "Plot"]), FUN = sum)
+# names(in_ba) <- c("Plot", "In_ba")
+# in_ba_data <-  join(in_ba, plot_data, by = c("Plot"))
 
 # soil data
 source("./calculate_awc.R")
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/lock_file b/.Rproj.user/EAEB7255/sources/s-197F81BE/lock_file
similarity index 100%
rename from .Rproj.user/EAEB7255/sources/s-9C9BC584/lock_file
rename to .Rproj.user/EAEB7255/sources/s-197F81BE/lock_file
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/40A653B0 b/.Rproj.user/EAEB7255/sources/s-9C9BC584/40A653B0
deleted file mode 100644
index e6e9694..0000000
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/40A653B0
+++ /dev/null
@@ -1,22 +0,0 @@
-{
-    "collab_server" : "",
-    "contents" : "",
-    "created" : 1565286118499.000,
-    "dirty" : false,
-    "encoding" : "UTF-8",
-    "folds" : "",
-    "hash" : "0",
-    "id" : "40A653B0",
-    "lastKnownWriteTime" : 1593814330,
-    "last_content_update" : 1593814330391,
-    "path" : "~/Research/MS Thesis/Understory/understory_effects_plots.R",
-    "project_path" : "understory_effects_plots.R",
-    "properties" : {
-        "cursorPosition" : "74,42",
-        "scrollLine" : "53"
-    },
-    "relative_order" : 7,
-    "source_on_save" : false,
-    "source_window" : "",
-    "type" : "r_source"
-}
\ No newline at end of file
diff --git a/.Rproj.user/EAEB7255/sources/s-9C9BC584/E4AFA6C6 b/.Rproj.user/EAEB7255/sources/s-9C9BC584/E4AFA6C6
deleted file mode 100644
index 42099de..0000000
--- a/.Rproj.user/EAEB7255/sources/s-9C9BC584/E4AFA6C6
+++ /dev/null
@@ -1,29 +0,0 @@
-{
-    "collab_server" : "",
-    "contents" : "",
-    "created" : 1593692468100.000,
-    "dirty" : false,
-    "encoding" : "",
-    "folds" : "",
-    "hash" : "0",
-    "id" : "E4AFA6C6",
-    "lastKnownWriteTime" : 7308604914264011873,
-    "last_content_update" : 1593692468100,
-    "path" : null,
-    "project_path" : null,
-    "properties" : {
-        "cacheKey" : "E6FB805F",
-        "caption" : "plot_data",
-        "contentUrl" : "grid_resource/gridviewer.html?env=&obj=plot_data&cache_key=E6FB805F",
-        "displayedObservations" : 102,
-        "environment" : "",
-        "object" : "plot_data",
-        "preview" : 0,
-        "totalObservations" : 102,
-        "variables" : 30
-    },
-    "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/s-9C9BC584/E4AFA6C6-contents b/.Rproj.user/EAEB7255/sources/s-9C9BC584/E4AFA6C6-contents
deleted file mode 100644
index e69de29..0000000
diff --git a/.Rproj.user/EAEB7255/viewer-cache/E6FB805F.Rdata b/.Rproj.user/EAEB7255/viewer-cache/E6FB805F.Rdata
deleted file mode 100644
index 453f526..0000000
Binary files a/.Rproj.user/EAEB7255/viewer-cache/E6FB805F.Rdata and /dev/null differ
diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths
index 5c772b5..ae3664b 100644
--- a/.Rproj.user/shared/notebooks/paths
+++ b/.Rproj.user/shared/notebooks/paths
@@ -1,4 +1,6 @@
 C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory.R="3F501BF8"
+C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory_tc.R="B649CF02"
 C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_data_prep.R="318C3941"
 C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_effects_plots.R="9317131D"
+C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_effects_plots_tc.R="B77DF8FD"
 C:/Users/Sam/Documents/Research/MS Thesis/Understory/understory_electivity.R="A106E7C2"
diff --git a/outputs/Figure_2_cheatgrass_change_effect.tif b/outputs/Figure_2_cheatgrass_change_effect.tif
index e5078b5..d3a3715 100644
Binary files a/outputs/Figure_2_cheatgrass_change_effect.tif and b/outputs/Figure_2_cheatgrass_change_effect.tif differ
diff --git a/outputs/Figure_3_all_understory_effects.tiff b/outputs/Figure_3_all_understory_effects.tiff
index be22226..52cc91c 100644
Binary files a/outputs/Figure_3_all_understory_effects.tiff and b/outputs/Figure_3_all_understory_effects.tiff differ
diff --git a/outputs/Figure_4_understory_effects.tiff b/outputs/Figure_4_understory_effects.tiff
index 1c28e7a..3f390ab 100644
Binary files a/outputs/Figure_4_understory_effects.tiff and b/outputs/Figure_4_understory_effects.tiff differ
diff --git a/outputs/Figure_5_electivity_test.tiff b/outputs/Figure_5_electivity_test.tiff
new file mode 100644
index 0000000..69e5d1c
Binary files /dev/null and b/outputs/Figure_5_electivity_test.tiff differ
diff --git a/outputs/Supplementary figure effects cgrass.tiff b/outputs/Supplementary figure effects cgrass.tiff
new file mode 100644
index 0000000..0c6c3fe
Binary files /dev/null and b/outputs/Supplementary figure effects cgrass.tiff differ
diff --git a/outputs/all_understory_effects.tiff b/outputs/all_understory_effects.tiff
deleted file mode 100644
index 276d4f8..0000000
Binary files a/outputs/all_understory_effects.tiff and /dev/null differ
diff --git a/outputs/cheatgrass_change_effect_no_outlier.tif b/outputs/cheatgrass_change_effect_no_outlier.tif
deleted file mode 100644
index 6cb20f4..0000000
Binary files a/outputs/cheatgrass_change_effect_no_outlier.tif and /dev/null differ
diff --git a/outputs/cheatgrass_change_effect_with_outlier.tif b/outputs/cheatgrass_change_effect_with_outlier.tif
deleted file mode 100644
index 2f0cc06..0000000
Binary files a/outputs/cheatgrass_change_effect_with_outlier.tif and /dev/null differ
diff --git a/outputs/electivity summary table.csv b/outputs/electivity summary table.csv
new file mode 100644
index 0000000..d2a1e42
--- /dev/null
+++ b/outputs/electivity summary table.csv	
@@ -0,0 +1,16 @@
+"","FT","var","Dead","Live","Inter","D-L","D-I","L-I"
+"1","All","Empirical",-0.0743033885607182,-0.0278109494523768,-0.0555980269498611,-0.507713667149215,-0.583200773715881,-0.120598140688693
+"2","All","Difference",0.0595523857436357,0.00182766977481981,-0.0402022583005446,0.0410063082768805,0.0750098341371376,0.0315387471759683
+"3","All","p-val",0.189,0.484,0.04,0.201,0.024,0.206
+"4","Perennial grass","Empirical",-0.192674459715908,0.0159977449726142,-0.0956047526050771,-0.565051037512717,-0.580468770773183,-0.0578091868907063
+"5","Perennial grass","Difference",-0.00991214048926803,0.0633071495429758,-0.0715729859102648,-0.0168567040219004,0.0766639090453757,0.0947762978648043
+"6","Perennial grass","p-val",0.462,0.058,0.011,0.416,0.07,0.029
+"7","Cheatgrass","Empirical",-0.304836422090473,-0.282639621711735,-0.0945149761162297,-0.311111111111111,-0.381768825584231,-0.300064487499948
+"8","Cheatgrass","Difference",0.0213814705753234,-0.119879042994129,0.0163599657013229,0.0600085773710695,0.155554894091391,-0.108947866057829
+"9","Cheatgrass","p-val",0.453,0.175,0.451,0.369,0.164,0.263
+"10","Perennial forb","Empirical",-0.297111840871393,-0.351160282239902,-0.00331648580873326,-0.305092592592593,-0.517196267784503,-0.376805046805047
+"11","Perennial forb","Difference",-0.0483139050742002,-0.269849324449275,0.0567734664670934,0.125049831123907,0.0421834555196033,-0.19936984983386
+"12","Perennial forb","p-val",0.367,0.004,0.24,0.172,0.349,0.024
+"13","Shrub","Empirical",-0.433014652475664,-0.247885213803435,-0.0515697173985787,-0.408156755525177,-0.643507702854854,-0.325430330683954
+"14","Shrub","Difference",0.0532953063016339,-0.0928465788199659,0.035714791494117,0.172941632384489,0.0605932038748218,-0.14982426902294
+"15","Shrub","p-val",0.342,0.164,0.328,0.088,0.263,0.114
diff --git a/outputs/electivity.tiff b/outputs/electivity.tiff
deleted file mode 100644
index d0b9424..0000000
Binary files a/outputs/electivity.tiff and /dev/null differ
diff --git a/outputs/results_boots_1dead.rds b/outputs/results_boots_1dead.rds
index f675689..8baf38a 100644
Binary files a/outputs/results_boots_1dead.rds and b/outputs/results_boots_1dead.rds differ
diff --git a/outputs/understory_effects.tiff b/outputs/understory_effects.tiff
deleted file mode 100644
index f5c9f5e..0000000
Binary files a/outputs/understory_effects.tiff and /dev/null differ
diff --git a/plot_level_understory.R b/plot_level_understory.R
index 4fb1df5..d1cd503 100644
--- a/plot_level_understory.R
+++ b/plot_level_understory.R
@@ -18,7 +18,7 @@ library("plyr")
 
 # import data created by data prep script
 plot_data <- read.csv("./clean data/plot_data.csv")
-plot_data$Delta_tc <- plot_data$Tree_cover - plot_data$Total_cover
+plot_data$Delta_tc <- (plot_data$Tree_cover*100 - plot_data$Total_cover)
 plot_data$Delta_ba <- plot_data$Live_ba - plot_data$Live_ba_2005
 plot_data <- plot_data[plot_data$Cluster != "NPELECTRICEEL", ]
 
@@ -267,54 +267,78 @@ par(opar)
 
 ## all understory vegetation
 
-all_plot <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+all_plot_tc <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
                    scale(AWC) + (1|Cluster), 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))
+
+
+all_plot_noscale_tc <- lmer(asin(sqrt(All/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                   AWC + (1|Cluster), data = plot_data)
+
+summary(all_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(all_plot_tc)
+plot(allEffects(all_plot_tc, partial.residuals = TRUE))
+plot(all_plot_tc)
+AICc(all_plot_tc)
+scatter.smooth(residuals(all_plot_tc) ~ predict(all_plot_tc))
 
 ## Cheatgrass
 
-cheatgrass_plot <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+cheatgrass_plot_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc) * scale(cwd_normal_cum) +
                           scale(AWC) + (1|Cluster), 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))
+
+cheatgrass_plot_noscale_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                          AWC + (1|Cluster), data = plot_data)
+
+summary(cheatgrass_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(cheatgrass_plot_tc)
+plot(allEffects(cheatgrass_plot_noscale_tc, partial.residuals = TRUE))
+plot(cheatgrass_plot_tc)
+AICc(cheatgrass_plot_tc)
+scatter.smooth(residuals(cheatgrass_plot_tc) ~ predict(cheatgrass_plot_tc))
 
 ## Perr grass
 
-pgrass_plot <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+pgrass_plot_tc <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
                       scale(AWC) + (1|Cluster), data = plot_data)
-summary(pgrass_plot, ddf = "Kenward-Roger")
-r.squaredGLMM(pgrass_plot)
-plot(allEffects(pgrass_plot, partial.residuals = TRUE))
-plot(pgrass_plot)
+
+pgrass_plot_noscale_tc <- lmer(asin(sqrt(Pgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                      AWC + (1|Cluster), data = plot_data)
+
+summary(pgrass_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pgrass_plot_tc)
+plot(allEffects(pgrass_plot_tc, partial.residuals = TRUE))
+plot(pgrass_plot_tc)
+scatter.smooth(residuals(pgrass_plot_tc) ~ predict(pgrass_plot_tc))
+
 
 ## Perr forbs
 
-pforb_plot <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+pforb_plot_tc <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
                      scale(AWC) + (1|Cluster), 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))
+pforb_plot_noscale_tc <- lmer(asin(sqrt(Pforb/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                     AWC + (1|Cluster), data = plot_data)
+
+summary(pforb_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pforb_plot_tc)
+plot(allEffects(pforb_plot_tc, partial.residuals = TRUE))
+plot(inv.as(predict(pforb_plot_tc)) ~ I(plot_data$Pforb/100))
 abline(0,1)
-plot(pforb_plot)
+plot(pforb_plot_tc)
+scatter.smooth(residuals(pforb_plot_tc) ~ predict(pforb_plot_tc))
 
 
 ## Shrubs
-shrub_plot <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+shrub_plot_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_tc)*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))
+
+shrub_plot_noscale_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                     AWC + (1|Cluster), data = plot_data)
+
+
+summary(shrub_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(shrub_plot_tc)
+plot(allEffects(shrub_plot_tc, partial.residuals = TRUE))
+plot(inv.as(predict(shrub_plot_tc)) ~ I(plot_data$Shrub_cover_li))
 abline(0,1)
 
 #-----------------------------------------------
@@ -340,7 +364,7 @@ compare <- compare[-c(5, 10), ] #outlier to remove?
 
 cg_change_lm <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ Delta_tc*cwd_normal_cum, data= compare)
 cg_change_lm_no_outliers <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ 
-                                 Delta_tc*cwd_normal_cum, data= compare[-c(5,10), ])
+                                 Delta_tc + cwd_normal_cum, data= compare[-c(5,10), ])
 
 summary(cg_change_lm)
 plot(cg_change_lm)
diff --git a/plot_level_understory_tc.R b/plot_level_understory_tc.R
new file mode 100644
index 0000000..d1cd503
--- /dev/null
+++ b/plot_level_understory_tc.R
@@ -0,0 +1,423 @@
+# Plot-level understory analysis
+# author: Sam Flake
+# email: sflake@gmail.com
+# Description: this script runs all of the plot-level analysis for the study,
+# not including the electivity analysis. It takes the processed data from
+# understory_data_prep.R and does statistical analysis, as well as creating 
+# Figure 1.
+
+library("MuMIn")
+library("lme4")
+library("car")
+library("effects")
+library("plyr")
+library("lmerTest")
+library("vegan")
+library("reshape2")
+library("plyr")
+
+# import data created by data prep script
+plot_data <- read.csv("./clean data/plot_data.csv")
+plot_data$Delta_tc <- (plot_data$Tree_cover*100 - plot_data$Total_cover)
+plot_data$Delta_ba <- plot_data$Live_ba - plot_data$Live_ba_2005
+plot_data <- plot_data[plot_data$Cluster != "NPELECTRICEEL", ]
+
+# understory variables from Greenwood and Weisberg 2008
+greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv")
+
+#------------------------------------------------------------
+# Analysis
+#------------------------------------------------------------
+#---------------------------------------------------------------
+#linear mixed effects models for functional groups at plot level
+#-------------------------------------------------------------------
+
+## all understory vegetation
+
+all_plot <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_pdc)*scale(cwd_normal_cum) +
+                          scale(AWC) + (1|Cluster), 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)
+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))
+
+## 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(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(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(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)
+
+#-----------------------------------------------
+# Comparison between 2005 and 2015 samples
+#-----------------------------------------------
+compare <- join(plot_data, greenwood_under, by = "Plot", type = "inner")
+
+#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 <- compare[-c(5, 10), ] #outlier to remove?
+
+cg_change_lm <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ Delta_tc*cwd_normal_cum, data= compare)
+
+summary(cg_change_lm)
+plot(cg_change_lm)
+plot(allEffects(cg_change_lm, partial.residuals = TRUE))
+
+
+#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 ~ cwd_normal_cum + Delta_tc, data= compare)
+summary(pg_change_lm)
+plot(allEffects(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_tc*cwd_normal_cum, data= compare)
+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_tc*cwd_normal_cum, data= compare)
+summary(s_change_lm)
+
+######################
+# Plot of comparisons between 2005 and 2015
+# Figure 1
+######################
+opar <- par(no.readonly = TRUE)
+
+tiff(filename="./outputs/Figure_1_change_in_cover.tif", 
+    type="cairo",
+    units="in", 
+    width = 4, 
+    height=4, 
+    pointsize=12, 
+    compression = "lzw", 
+    res=600)
+
+layout(matrix(c(1,2,3,4), nrow = 2, ncol = 2, byrow = TRUE))
+
+par(oma = c(2,2,0,0), 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$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',
+     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')
+  
+
+
+mtext("2005 Cover (%)", side = 1, outer = TRUE)
+mtext("2015 Cover (%)", side = 2, outer = TRUE)
+dev.off()
+
+par(opar)
+
+
+
+
+
+#----------------------------------------------------------------------------------------
+# Redo the analysis with change in tree cover
+#----------------------------------------------------------------------------------------
+#---------------------------------------------------------------
+#linear mixed effects models for functional groups at plot level
+#-------------------------------------------------------------------
+
+## all understory vegetation
+
+all_plot_tc <- lmer(asin(sqrt(All/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+                   scale(AWC) + (1|Cluster), data = plot_data)
+
+
+all_plot_noscale_tc <- lmer(asin(sqrt(All/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                   AWC + (1|Cluster), data = plot_data)
+
+summary(all_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(all_plot_tc)
+plot(allEffects(all_plot_tc, partial.residuals = TRUE))
+plot(all_plot_tc)
+AICc(all_plot_tc)
+scatter.smooth(residuals(all_plot_tc) ~ predict(all_plot_tc))
+
+## Cheatgrass
+
+cheatgrass_plot_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc) * scale(cwd_normal_cum) +
+                          scale(AWC) + (1|Cluster), data = plot_data)
+
+cheatgrass_plot_noscale_tc <- lmer(asin(sqrt(Cheatgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                          AWC + (1|Cluster), data = plot_data)
+
+summary(cheatgrass_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(cheatgrass_plot_tc)
+plot(allEffects(cheatgrass_plot_noscale_tc, partial.residuals = TRUE))
+plot(cheatgrass_plot_tc)
+AICc(cheatgrass_plot_tc)
+scatter.smooth(residuals(cheatgrass_plot_tc) ~ predict(cheatgrass_plot_tc))
+
+## Perr grass
+
+pgrass_plot_tc <- lmer(asin(sqrt(Pgrass/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+                      scale(AWC) + (1|Cluster), data = plot_data)
+
+pgrass_plot_noscale_tc <- lmer(asin(sqrt(Pgrass/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                      AWC + (1|Cluster), data = plot_data)
+
+summary(pgrass_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pgrass_plot_tc)
+plot(allEffects(pgrass_plot_tc, partial.residuals = TRUE))
+plot(pgrass_plot_tc)
+scatter.smooth(residuals(pgrass_plot_tc) ~ predict(pgrass_plot_tc))
+
+
+## Perr forbs
+
+pforb_plot_tc <- lmer(asin(sqrt(Pforb/100)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+                     scale(AWC) + (1|Cluster), data = plot_data)
+pforb_plot_noscale_tc <- lmer(asin(sqrt(Pforb/100)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                     AWC + (1|Cluster), data = plot_data)
+
+summary(pforb_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(pforb_plot_tc)
+plot(allEffects(pforb_plot_tc, partial.residuals = TRUE))
+plot(inv.as(predict(pforb_plot_tc)) ~ I(plot_data$Pforb/100))
+abline(0,1)
+plot(pforb_plot_tc)
+scatter.smooth(residuals(pforb_plot_tc) ~ predict(pforb_plot_tc))
+
+
+## Shrubs
+shrub_plot_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ scale(Tree_cover) + scale(Delta_tc)*scale(cwd_normal_cum) +
+                     scale(AWC) + (1|Cluster), data = plot_data)
+
+shrub_plot_noscale_tc <- lmer(asin(sqrt(Shrub_cover_li)) ~ Tree_cover + Delta_tc*cwd_normal_cum +
+                     AWC + (1|Cluster), data = plot_data)
+
+
+summary(shrub_plot_tc, ddf = "Kenward-Roger")
+r.squaredGLMM(shrub_plot_tc)
+plot(allEffects(shrub_plot_tc, partial.residuals = TRUE))
+plot(inv.as(predict(shrub_plot_tc)) ~ I(plot_data$Shrub_cover_li))
+abline(0,1)
+
+#-----------------------------------------------
+# Comparison between 2005 and 2015 samples
+#-----------------------------------------------
+compare <- join(plot_data, greenwood_under, by = "Plot", type = "inner")
+
+#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 <- compare[-c(5, 10), ] #outlier to remove?
+
+cg_change_lm <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ Delta_tc*cwd_normal_cum, data= compare)
+cg_change_lm_no_outliers <- lm(I((Cheatgrass - Cheatgrass.Cover)/100) ~ 
+                                 Delta_tc + cwd_normal_cum, data= compare[-c(5,10), ])
+
+summary(cg_change_lm)
+plot(cg_change_lm)
+plot(allEffects(cg_change_lm, partial.residuals = TRUE))
+
+summary(cg_change_lm_no_outliers)
+plot(cg_change_lm)
+plot(allEffects(cg_change_lm_no_outliers, partial.residuals = TRUE))
+
+#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
+
+#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_tc*cwd_normal_cum, data= compare)
+summary(pg_change_lm)
+plot(allEffects(pg_change_lm, partial.residuals = TRUE))
+
+## 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_tc*cwd_normal_cum, data= compare)
+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_tc*cwd_normal_cum, data= compare)
+summary(s_change_lm)
diff --git a/understory_data_prep.R b/understory_data_prep.R
index 29c89e9..0407209 100644
--- a/understory_data_prep.R
+++ b/understory_data_prep.R
@@ -28,6 +28,9 @@ 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")
 
+# individual tree data from Greenwood and Weisberg 2008
+greenwood_trees <- read.csv("./raw data/Individual_TreeData_AllData_SF_edits.csv")
+
 # 2005 understory data, from Greenwood and Weisberg, 2008, 
 # Density-dependent tree mortality in pinyon-juniper woodlands, FEM
 greenwood_under <- read.csv("./raw data/Greenwood_Understory_Variables_SF_edits.csv")
@@ -39,8 +42,11 @@ tcover <- read.csv("./raw data/tree_and_shrub_cover_020815.csv")
 # 2005 tree cover and other variables from Greenwood and Weisberg 2008
 greenwood_tree_cover <- read.csv("./raw data/Structure_vs_environmental_all_data_SF_edits.csv")
   
-# individual tree data from Greenwood and Weisberg 2008
-greenwood_trees <- read.csv("./raw data/Individual_TreeData_AllData_SF_edits.csv")
+# calculate basal area of ingrowths
+# in_ba <- aggregate(trees[trees$Code == "P", "BA_m2_per_ha"], 
+#           by = list(trees[trees$Code == "P", "Plot"]), FUN = sum)
+# names(in_ba) <- c("Plot", "In_ba")
+# in_ba_data <-  join(in_ba, plot_data, by = c("Plot"))
 
 # soil data
 source("./calculate_awc.R")
diff --git a/understory_effects_plots.R b/understory_effects_plots.R
index 374b55f..7e3bd73 100644
--- a/understory_effects_plots.R
+++ b/understory_effects_plots.R
@@ -20,11 +20,11 @@ source('addTrans.R', echo=FALSE)
 # Some other global variables are also needed, like the raw plot_data data frame.
 # Bad practices but I didn't think it was worth refactoring this thing. 
 
-all <- all_plot
-cg <- cheatgrass_plot
-pg <- pgrass_plot
-pf <- pforb_plot
-sh <- shrub_plot
+all <- all_plot_noscale_tc
+cg <- cheatgrass_plot_noscale
+pg <- pgrass_plot_noscale
+pf <- pforb_plot_noscale
+sh <- shrub_plot_noscale
 
 # 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
@@ -49,7 +49,7 @@ closest <- function(x, x0){
 # Figure 3
 #####################################################################
 vars <- c("Tree_cover", "Delta_tc")
-labels <- c("Tree cover", "Change in live canopy (%)")
+labels <- c("Tree cover", "Change in tree cover (%)")
 
 tiff(filename="./outputs/Figure_3_all_understory_effects.tiff", 
      type="cairo",
@@ -65,7 +65,7 @@ par(mfrow = c(2,1), oma = c(2,3,0,0), mar = c(3,1,1,1), bty = 'n')
 for(i in c(1:2)){
   var <- vars[i]
   
-  eff <- Effect(all_plot, partial.residuals = TRUE, focal.predictors = var)
+  eff <- Effect(all_plot_tc, partial.residuals = TRUE, focal.predictors = var)
   
   y <- eff$fit
   x <- eff$x[[var]]
@@ -80,7 +80,7 @@ for(i in c(1:2)){
        xlab = "",
        ylab = "")
   
-  points(resids*100 ~ eff$x.all[[var]], 
+  points(resids*100 ~ x.fit, 
          pch = 21, 
          bg='grey60',
          col = 'grey30')
@@ -169,7 +169,7 @@ calculate_preds <- function(model, C, predictor){
 
 # 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
+# this is the dumbest way to do this but I don't want to refactor 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")
@@ -195,10 +195,6 @@ preds_tc_90_pg <- calculate_preds(pg, C_tc_90, "Delta_tc")
 preds_tc_90_pf <- calculate_preds(pf, C_tc_90, "Delta_tc")
 preds_tc_90_sh <- calculate_preds(sh, C_tc_90, "Delta_tc")
 
-# 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
-
 opar <- par(no.readonly = TRUE)
 par(opar)
 
@@ -516,7 +512,7 @@ dev.off()
 library("effects")
 # Perennial grasses
 
-tiff(filename="./outputs/Supplementary figure effects pgrass.tiff", 
+tiff(filename="./outputs/Supplementary figure effects cgrass.tiff", 
      type="cairo",
      units="in", 
      width = 7, 
@@ -525,135 +521,141 @@ tiff(filename="./outputs/Supplementary figure effects pgrass.tiff",
      compression = "lzw",
      res=600)
 
-model <- cg
-
-eff <- Effect(focal.predictors = c("Tree_cover"), 
-              mod = cg, xlevels = 100, partial.residuals = TRUE)
-
-
-eff2 <- Effect(focal.predictors = c("PC1"), 
-               mod = mod3, xlevels = 100, transformation = list(link = log, inverse = exp))
-
-dat <- data.frame(y = exp(eff$fit),
-                  lower = exp(eff$lower),
-                  upper = exp(eff$upper),
-                  soil = eff[["x"]][["PC1"]])
-
+#graphical params
 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), bty = 'n')
 
 
-plot(NA,
-     ylim = c(0,.02),
-     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,200)]*100), inv.as(preds_tc_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-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,.02),
-     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,200)]), inv.as(preds_awc_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-mtext(text = "Soil AWC (%)", side = 1, line = 2.2)
-mtext(text = "(b)", side = 1, line = -10, adj = 0.05)
+panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
+k <- 1
 
+#which  model to use?
+model <- pgrass_plot_noscale
 
-plot.new()
-legend("topright", legend = c("Cheatgrass", "Per. Grass", "Per. Forb", "Shrub"), 
-       lty = 1, pch = c(21,22,23,24), lwd = 2, cex = 1.3, col = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"),
-       pt.bg = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"))
+let <- letters
 
-plot(NA,
-     ylim = c(0,.02),
-     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
-     xlab = "",
-     ylab = "",
-     bg='grey60',
-     col = 'grey30',
-     pch = 21,
-     cex.lab = 1.5,
-     bty = 'n',
-     xaxt = 'n',
-     yaxt = 'n')
-lines(inv.as(preds_tc_10_cg$predictions) ~ I(preds_tc_10_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_tc_10_cg$predictor[c(1,200)]), inv.as(preds_tc_10_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+vars <- c("Tree_cover", "AWC", "Delta_tc", "Delta_tc", "Delta_tc")
+x_labels <- c("Tree cover (%)", "Soil AWC (%)", "Change in live canopy (%)")
 
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-text(x = -30, y = .018, labels = "10% CWD", cex = 1.3)
-mtext(text = "(c)", side = 1, line = -10, adj = 0.05)
+for(i in 1:2){
 
-plot(NA,
-     ylim = c(0,.02),
-     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
-     xlab = "",
-     ylab = "",
-     bg='grey60',
-     col = 'grey30',
-     pch = 21,
-     cex.lab = 1.5,
-     bty = 'n',
-     xaxt = 'n',
-     yaxt = 'n')
-lines(inv.as(preds_tc_50_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_50_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  var <- vars[i]
+  
+  eff <- Effect(model, partial.residuals = TRUE, focal.predictors = var)
+  
+  line_data <- data.frame(y = inv.as(eff$fit),
+                    lower = inv.as(eff$lower),
+                    upper = inv.as(eff$upper),
+                    x = eff[["x"]][[var]]
+                    )
+  
+  x_new <- eff[["x"]][[var]]
+  x_orig <- eff[["x.all"]][[var]]
+  
+  y_trans <- eff$fit
+  
+  fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+  # that we can find the closest fitted value and add the residual for that 
+  
+  p_resids <- inv.as(eff$residuals + fitted)
+  # eff[["data"]][[var]]
+  
+  points_data <- data.frame(x = x_orig,
+                            y = p_resids)
+  
+  plot(NA,
+       ylim = c(min(c(y, points_data$y)), max(c(y, points_data$y))),
+       xlim = c(min(x_new), max(x_new)),
+       xlab = "",
+       ylab = "",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex.lab = 1.5,
+       bty = 'n',
+       xaxt = 'n',
+       yaxt = 'n')
+  
+  lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(line_data$upper ~ line_data$x)
+  lines(line_data$lower ~ line_data$x)
+  polygon(c(line_data$x, rev(line_data$x)), 
+          c(line_data$upper, rev(line_data$lower)),
+          col = addTrans("#68EBC4",30), border = NA)
+  
+  points(points_data$y ~ points_data$x,
+         pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  
+  axis(side = 1)
+  axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+  mtext(text = var, 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)
 
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-text(x = -30, y = .018, 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)
 
-plot(NA,
-     ylim = c(0,.02),
-     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
-     xlab = "",
-     ylab = "",
-     bg='grey60',
-     col = 'grey30',
-     pch = 21,
-     cex.lab = 1.5,
-     bty = 'n',
-     xaxt = 'n',
-     yaxt = 'n')
-lines(inv.as(preds_tc_90_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
-points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
-
-axis(side = 1)
-axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
-text(x = -30, y = .018, labels = "90% CWD", cex = 1.3)
-mtext(text = "(e)", side = 1, line = -10, adj = 0.05)
+for(i in 1:3){
+  
+  eff <- Effect(model, partial.residuals = TRUE, focal.predictors = c("Delta_tc", "cwd_normal_cum"),
+                xlevels = list(Delta_tc = 100),
+                quantiles = c(0.1, 0.5, 0.9))
+  
+  cwd_levels <- eff$variables$cwd_normal_cum$levels
+  
+  line_data <- data.frame(y = inv.as(eff$fit[1:99 + (i-1) * 100]),
+                          lower = inv.as(eff$lower[1:99 + (i-1) * 100]),
+                          upper = inv.as(eff$upper[1:99 + (i-1) * 100]),
+                          x = eff[["x"]][["Delta_tc"]][1:99 + (i-1) * 100]
+  )
+  
+  x_new <- eff[["x"]][["Delta_tc"]]
+  x_orig_all <- eff[["x.all"]]
+  which_x <- x_orig_all$cwd_normal_cum - cwd_levels[i] < 0.001
+  x_orig <- x_orig_all[which_x, ][["Delta_tc"]]
+  
+  y_trans <- eff$fit[1:99 + (i-1) * 100]
+  
+  fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+  # that we can find the closest fitted value and add the residual for that 
+  
+  p_resids <- inv.as(eff$residuals[which_x] + fitted)
+  
+  points_data <- data.frame(x = x_orig,
+                            y = p_resids)
+  
+  plot(NA,
+       ylim = c(0, max(c(y, points_data$y), na.rm = TRUE)),
+       xlim = c(min(x_new), max(x_new)),
+       xlab = "",
+       ylab = "",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex.lab = 1.5,
+       bty = 'n',
+       xaxt = 'n',
+       yaxt = 'n')
+  
+  lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(line_data$upper ~ line_data$x)
+  lines(line_data$lower ~ line_data$x)
+  polygon(c(line_data$x, rev(line_data$x)), 
+          c(line_data$upper, rev(line_data$lower)),
+          col = addTrans("#68EBC4",30), border = NA)
+  
+  points(points_data$y ~ points_data$x,
+         pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  
+  axis(side = 1)
+  axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+  mtext(text = var, 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)
 
+}
 
 dev.off()
+
diff --git a/understory_effects_plots_tc.R b/understory_effects_plots_tc.R
new file mode 100644
index 0000000..7e3bd73
--- /dev/null
+++ b/understory_effects_plots_tc.R
@@ -0,0 +1,661 @@
+# Plot-level understory analysis
+# author: Sam Flake
+# email: sflake@gmail.com
+# Description: this script generates Figures 2, 3, and 4. It relies on the effects
+# package or on code borrowed from that package. You'll have to run plot_level_understory.R
+# first, because this script needs a lot of global variables generated during that
+# analysis. Outputs some .tiff files. 
+
+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
+# These models were generated in plot_level_understory.R, and that
+# script needs to be run first and these models saved in the global environment.
+# Some other global variables are also needed, like the raw plot_data data frame.
+# Bad practices but I didn't think it was worth refactoring this thing. 
+
+all <- all_plot_noscale_tc
+cg <- cheatgrass_plot_noscale
+pg <- pgrass_plot_noscale
+pf <- pforb_plot_noscale
+sh <- shrub_plot_noscale
+
+# 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
+}
+
+# function to find the nearest y-value for an x data point, to add to partial residuals
+# for plotting. Stolen from effects package
+closest <- function(x, x0){
+  apply(outer(x, x0, FUN = function(x, x0) abs(x - x0)), 1, which.min)
+}
+
+#####################################################################
+# All understory vegetation, just two panels
+# Figure 3
+#####################################################################
+vars <- c("Tree_cover", "Delta_tc")
+labels <- c("Tree cover", "Change in tree cover (%)")
+
+tiff(filename="./outputs/Figure_3_all_understory_effects.tiff", 
+     type="cairo",
+     units="in", 
+     width = 3, 
+     height=5, 
+     pointsize=12,
+     compression = "lzw", 
+     res=600)
+
+par(mfrow = c(2,1), oma = c(2,3,0,0), mar = c(3,1,1,1), bty = 'n')
+
+for(i in c(1:2)){
+  var <- vars[i]
+  
+  eff <- Effect(all_plot_tc, partial.residuals = TRUE, focal.predictors = var)
+  
+  y <- eff$fit
+  x <- eff$x[[var]]
+  x.fit <- eff$x.all[[var]]
+  
+  fitted <- y[closest(x.fit, x)]
+  resids <- inv.as(eff$residuals + fitted)
+  
+  plot(NA,
+       xlim = c(min(x), max(x)),
+       ylim = c(0, 50),
+       xlab = "",
+       ylab = "")
+  
+  points(resids*100 ~ x.fit, 
+         pch = 21, 
+         bg='grey60',
+         col = 'grey30')
+  
+  lines(inv.as(y)*100 ~ x, lwd = 2)
+  lines(inv.as(eff$lower)*100 ~ x)
+  lines(inv.as(eff$upper)*100 ~ x)
+  polygon(c(x, rev(x)), c(inv.as(eff$upper)*100, rev(inv.as(eff$lower))*100),
+          col = addTrans("#68EBC4",30), border = NA)
+
+  mtext(side = 2, text = "Understory cover (%)", outer = TRUE, line = 1.7)
+  mtext(side = 1, text = labels[i], line = 2)
+}
+
+dev.off()
+
+
+####################################################################
+## Multipanel figure to compare between functional types
+# Figure 4
+####################################################################
+# 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_tc_mean <-  mean(plot_data$Delta_tc)
+Delta_tc_sd <- sd(plot_data$Delta_tc)
+
+#--------------------------------------------------------------------------
+## Calcuate partial residuals, effects
+#--------------------------------------------------------------------------
+
+C_tc <- data.frame(
+  Tree_cover = seq(min(plot_data$Tree_cover), max(plot_data$Tree_cover), length.out = 200),
+  cwd_normal_cum = mean(plot_data$cwd_normal_cum),
+  Delta_tc = mean(plot_data$Delta_tc),
+  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_tc = mean(plot_data$Delta_tc),
+  AWC = seq(min(plot_data$AWC), max(plot_data$AWC), length.out = 200)
+)
+
+C_tc_10 <- data.frame(
+  Tree_cover = mean(plot_data$Tree_cover),
+  cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .1)),
+  Delta_tc = seq(min(plot_data$Delta_tc), max(plot_data$Delta_tc), length.out = 200),
+  AWC = mean(plot_data$AWC)
+)
+
+C_tc_50 <- data.frame(
+  Tree_cover = mean(plot_data$Tree_cover),
+  cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .5)),
+  Delta_tc = seq(min(plot_data$Delta_tc), max(plot_data$Delta_tc), length.out = 200),
+  AWC = mean(plot_data$AWC)
+)
+
+C_tc_90 <- data.frame(
+  Tree_cover = mean(plot_data$Tree_cover),
+  cwd_normal_cum = unname(quantile(plot_data$cwd_normal_cum, .9)),
+  Delta_tc = seq(min(plot_data$Delta_tc), max(plot_data$Delta_tc), length.out = 200),
+  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 refactor 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_tc_10_cg <- calculate_preds(cg, C_tc_10, "Delta_tc")
+preds_tc_10_pg <- calculate_preds(pg, C_tc_10, "Delta_tc")
+preds_tc_10_pf <- calculate_preds(pf, C_tc_10, "Delta_tc")
+preds_tc_10_sh <- calculate_preds(sh, C_tc_10, "Delta_tc")
+
+preds_tc_50_cg <- calculate_preds(cg, C_tc_50, "Delta_tc")
+preds_tc_50_pg <- calculate_preds(pg, C_tc_50, "Delta_tc")
+preds_tc_50_pf <- calculate_preds(pf, C_tc_50, "Delta_tc")
+preds_tc_50_sh <- calculate_preds(sh, C_tc_50, "Delta_tc")
+
+preds_tc_90_cg <- calculate_preds(cg, C_tc_90, "Delta_tc")
+preds_tc_90_pg <- calculate_preds(pg, C_tc_90, "Delta_tc")
+preds_tc_90_pf <- calculate_preds(pf, C_tc_90, "Delta_tc")
+preds_tc_90_sh <- calculate_preds(sh, C_tc_90, "Delta_tc")
+
+opar <- par(no.readonly = TRUE)
+par(opar)
+
+tiff(filename="./outputs/Figure_4_understory_effects.tiff", 
+    type="cairo",
+    units="in", 
+    width = 7, 
+    height=5, 
+    pointsize=15,
+    compression = "lzw",
+    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), 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,200)]*100), inv.as(preds_tc_cg$predictions[c(1,200)]), 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,200)]*100), inv.as(preds_tc_pg$predictions[c(1,200)]), 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,200)]*100), inv.as(preds_tc_pf$predictions[c(1,200)]), 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,200)]*100), inv.as(preds_tc_sh$predictions[c(1,200)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
+
+axis(side = 1)
+axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+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,200)]), inv.as(preds_awc_cg$predictions[c(1,200)]), 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,200)]), inv.as(preds_awc_pg$predictions[c(1,200)]), pch = 22, 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,200)]), inv.as(preds_awc_pf$predictions[c(1,200)]), 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,200)]), inv.as(preds_awc_sh$predictions[c(1,200)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
+
+axis(side = 1)
+axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+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 = 1, pch = c(21,22,23,24), lwd = 2, cex = 1.3, col = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"),
+       pt.bg = c("#1b9e77", "#000000", "#E69F00", "#56B4E9"))
+
+plot(NA,
+     ylim = c(0,.25),
+     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
+     xlab = "",
+     ylab = "",
+     bg='grey60',
+     col = 'grey30',
+     pch = 21,
+     cex.lab = 1.5,
+     bty = 'n',
+     xaxt = 'n',
+     yaxt = 'n')
+lines(inv.as(preds_tc_10_cg$predictions) ~ I(preds_tc_10_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
+points(I(preds_tc_10_cg$predictor[c(1,200)]), inv.as(preds_tc_10_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+lines(inv.as(preds_tc_10_pg$predictions) ~ I(preds_tc_10_cg$predictor), lwd = 2, lty = 1, col = "#000000")
+points(I(preds_tc_10_cg$predictor[c(1,200)]), inv.as(preds_tc_10_pg$predictions[c(1,200)]), pch = 22, col = "#000000", bg = "#000000")
+lines(inv.as(preds_tc_10_pf$predictions) ~ I(preds_tc_10_cg$predictor), lwd = 2, lty = 2, col = "#E69F00")
+points(I(preds_tc_10_cg$predictor[c(1,200)]), inv.as(preds_tc_10_pf$predictions[c(1,200)]), pch = 23, col = "#E69F00", bg = "#E69F00")
+lines(inv.as(preds_tc_10_sh$predictions) ~ I(preds_tc_10_cg$predictor), lwd = 2, lty = 2, col = "#56B4E9")
+points(I(preds_tc_10_cg$predictor[c(1,200)]), inv.as(preds_tc_10_sh$predictions[c(1,200)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
+
+axis(side = 1)
+axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+text(x = -30, y = .18, labels = "10% CWD", cex = 1.3)
+mtext(text = "(c)", side = 1, line = -10, adj = 0.05)
+
+plot(NA,
+     ylim = c(0,.25),
+     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
+     xlab = "",
+     ylab = "",
+     bg='grey60',
+     col = 'grey30',
+     pch = 21,
+     cex.lab = 1.5,
+     bty = 'n',
+     xaxt = 'n',
+     yaxt = 'n')
+lines(inv.as(preds_tc_50_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_50_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+lines(inv.as(preds_tc_50_pg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#000000")
+points(I(preds_tc_50_cg$predictor[c(1,200)]), inv.as(preds_tc_50_pg$predictions[c(1,200)]), pch = 22, col = "#000000", bg = "#000000")
+lines(inv.as(preds_tc_50_pf$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 2, col = "#E69F00")
+points(I(preds_tc_50_cg$predictor[c(1,200)]), inv.as(preds_tc_50_pf$predictions[c(1,200)]), pch = 23, col = "#E69F00", bg = "#E69F00")
+lines(inv.as(preds_tc_50_sh$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 2, col = "#56B4E9")
+points(I(preds_tc_50_cg$predictor[c(1,200)]), inv.as(preds_tc_50_sh$predictions[c(1,200)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
+
+axis(side = 1)
+axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+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,.25),
+     xlim = c(min(plot_data$Delta_tc), max(plot_data$Delta_tc)),
+     xlab = "",
+     ylab = "",
+     bg='grey60',
+     col = 'grey30',
+     pch = 21,
+     cex.lab = 1.5,
+     bty = 'n',
+     xaxt = 'n',
+     yaxt = 'n')
+lines(inv.as(preds_tc_90_cg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#1b9e77")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_cg$predictions[c(1,200)]), pch = 21, col = "#1b9e77", bg = "#1b9e77")
+lines(inv.as(preds_tc_90_pg$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 1, col = "#000000")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_pg$predictions[c(1,200)]), pch = 22, col = "#000000", bg = "#000000")
+lines(inv.as(preds_tc_90_pf$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 2, col = "#E69F00")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_pf$predictions[c(1,200)]), pch = 23, col = "#E69F00", bg = "#E69F00")
+lines(inv.as(preds_tc_90_sh$predictions) ~ I(preds_tc_90_cg$predictor), lwd = 2, lty = 2, col = "#56B4E9")
+points(I(preds_tc_90_cg$predictor[c(1,200)]), inv.as(preds_tc_90_sh$predictions[c(1,200)]), pch = 24, col = "#56B4E9", bg = "#56B4E9")
+
+axis(side = 1)
+axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+text(x = -30, y = .18, labels = "90% CWD", cex = 1.3)
+mtext(text = "(e)", side = 1, line = -10, adj = 0.05)
+
+
+dev.off()
+
+par(pardefault)
+
+
+###################################################################################
+# Cheatgrass change over time
+# Figure 2
+###################################################################################
+## new data
+model <- cg_change_lm
+model2 <- cg_change_lm_no_outliers
+
+# model <- model2
+
+compare$delta_cg <- compare$Cheatgrass - compare$Cheatgrass.Cover
+
+
+compare2 <- compare[-c(5, 10), ]
+compare_outliers <- compare[c(5,10), ]
+
+#create dataframe of new data to calculate predictions and envelope
+C_cg_change_10 <- data.frame(
+  cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .15)),
+  Delta_tc = seq(-50,10, length.out = 200)
+)
+
+C_cg_change_50 <- data.frame(
+  cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .5)),
+  Delta_tc = seq(-50,10, length.out = 200)
+)
+
+
+C_cg_change_90 <- data.frame(
+  cwd_normal_cum = unname(quantile(compare$cwd_normal_cum, .85)),
+  Delta_tc = seq(-50,10, length.out = 200)
+)
+
+C_cg <- list(C_cg_change_10, C_cg_change_50, C_cg_change_90)
+
+#function to calculate predictions and envelope
+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)
+}
+
+#calculate predictions and prediction envelope for each level of CWD
+cg_change_10 <- preds_cg_change(cg_change_lm, C_cg_change_10, "Delta_tc")
+cg_change_50 <- preds_cg_change(cg_change_lm, C_cg_change_50, "Delta_tc")
+cg_change_90 <- preds_cg_change(cg_change_lm, C_cg_change_90, "Delta_tc")
+
+cg_change_10_no_out <- preds_cg_change(cg_change_lm_no_outliers, C_cg_change_10, "Delta_tc")
+cg_change_50_no_out <- preds_cg_change(cg_change_lm_no_outliers, C_cg_change_50, "Delta_tc")
+cg_change_90_no_out <- preds_cg_change(cg_change_lm_no_outliers, C_cg_change_90, "Delta_tc")
+
+cg_change <- list(cg_change_10, cg_change_50, cg_change_90)
+cg_change_no_out <- list(cg_change_10_no_out, cg_change_50_no_out, cg_change_90_no_out)
+
+#generate the figure
+tiff(filename="./outputs/Figure_2_cheatgrass_change_effect.tif", 
+    type="cairo",
+    units="in", 
+    width = 4, 
+    height=3, 
+    pointsize=12,
+    compression = "lzw", 
+    res=600)
+
+layout(matrix(c(1,2,3), nrow=1, ncol=3, byrow = TRUE))
+par(mar = c(0,1,0,0), oma = c(6,4,2,2), bty = 'n')
+
+panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
+
+for(i in c(1:3)){
+  plot(NA,
+       ylim = c(-0.06, 0.16),
+       xlim = c(-70, 10 ),
+       xlab = "",
+       ylab = "",
+       xaxt = "n",
+       yaxt = "n",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex = 0.8,
+       cex.lab = 1.5)
+  abline(h = 0, lty = 4, lwd = 1.3)
+  
+  lines(cg_change[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 4, col = "#1b9e77")
+  lines(cg_change_no_out[[i]]$predictions ~ cg_change[[i]]$predictor, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(cg_change_no_out[[i]]$lo ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
+  lines(cg_change_no_out[[i]]$hi ~ cg_change[[i]]$predictor, lwd = 1, lty = 1, col = "#1b9e77")
+  polygon(c(cg_change_no_out[[i]]$predictor, rev(cg_change_no_out[[i]]$predictor)), c(cg_change_no_out[[i]]$hi, 
+                                                                    rev(cg_change_no_out[[i]]$lo)),
+          col = addTrans("#68EBC4",30), border = NA) #fills in the area between high and low confidence intervals
+  
+  # I also hate this: plot the residuals ~ delta_tc for each level of CWD
+      if(i == 1){
+        points(compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$delta_cg/100 ~ 
+          compare2[compare2$cwd_normal_cum < unname(quantile(compare2$cwd_normal_cum, 0.3)),]$Delta_tc,
+          bg='grey60',
+          col = 'grey30',
+          pch = 21,
+          cex = 1.2)
+      }else if(i ==2){
+        points(compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$delta_cg/100 ~ 
+          compare2[findInterval(compare2$cwd_normal_cum, unname(quantile(compare2$cwd_normal_cum, c(0.3, 0.7)))) == 1,]$Delta_tc,
+          bg='grey60',
+          col = 'grey30',
+          pch = 21,
+          cex = 1.2)
+        points(compare_outliers[1, ]$delta_cg/100 ~ compare_outliers[1, ]$Delta_tc, 
+               bg='white',
+               col = 'grey30',
+               pch = 21,
+               cex = 1.2)
+      }else if (i == 3){
+        points(compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$delta_cg/100 ~ 
+          compare2[compare2$cwd_normal_cum > unname(quantile(compare2$cwd_normal_cum, 0.7)),]$Delta_tc,
+          bg='grey60',
+          col = 'grey30',
+          pch = 21,
+          cex = 1.2)
+        points(compare_outliers[2, ]$delta_cg/100 ~ compare_outliers[2, ]$Delta_tc, 
+               bg='white',
+               col = 'grey30',
+               pch = 21,
+               cex = 1.2)
+      }
+  
+ if(i == 1){
+   axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+   # axis(side = 2, at = c(-0.1, 0, 0.1), labels = c(-10, 0, 10))
+   }
+ axis(side = 1, at = c(-40, -20, 0), labels = TRUE)
+ axis(side = 1, at = c(-50, -30, -10, 10), labels = FALSE, tcl = -.25)
+ text(x = -30, y = .07, labels = panel_labels[i], cex = 1.2)
+
+}
+   
+  mtext(side = 1, text = "Change in live canopy (%)", outer = TRUE, line = 2.6)
+  mtext(side = 2, text = "Change in cheatgrass cover (%)", outer = TRUE, line = 2.3)
+
+dev.off()
+
+
+
+#------------------------------------------------------------------------------
+# Effects plots with partial residuals for supplement
+#------------------------------------------------------------------------------
+library("effects")
+# Perennial grasses
+
+tiff(filename="./outputs/Supplementary figure effects cgrass.tiff", 
+     type="cairo",
+     units="in", 
+     width = 7, 
+     height=5, 
+     pointsize=15,
+     compression = "lzw",
+     res=600)
+
+#graphical params
+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), bty = 'n')
+
+
+panel_labels <- c("10% CWD", "50% CWD", "90% CWD")
+k <- 1
+
+#which  model to use?
+model <- pgrass_plot_noscale
+
+let <- letters
+
+vars <- c("Tree_cover", "AWC", "Delta_tc", "Delta_tc", "Delta_tc")
+x_labels <- c("Tree cover (%)", "Soil AWC (%)", "Change in live canopy (%)")
+
+for(i in 1:2){
+
+  var <- vars[i]
+  
+  eff <- Effect(model, partial.residuals = TRUE, focal.predictors = var)
+  
+  line_data <- data.frame(y = inv.as(eff$fit),
+                    lower = inv.as(eff$lower),
+                    upper = inv.as(eff$upper),
+                    x = eff[["x"]][[var]]
+                    )
+  
+  x_new <- eff[["x"]][[var]]
+  x_orig <- eff[["x.all"]][[var]]
+  
+  y_trans <- eff$fit
+  
+  fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+  # that we can find the closest fitted value and add the residual for that 
+  
+  p_resids <- inv.as(eff$residuals + fitted)
+  # eff[["data"]][[var]]
+  
+  points_data <- data.frame(x = x_orig,
+                            y = p_resids)
+  
+  plot(NA,
+       ylim = c(min(c(y, points_data$y)), max(c(y, points_data$y))),
+       xlim = c(min(x_new), max(x_new)),
+       xlab = "",
+       ylab = "",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex.lab = 1.5,
+       bty = 'n',
+       xaxt = 'n',
+       yaxt = 'n')
+  
+  lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(line_data$upper ~ line_data$x)
+  lines(line_data$lower ~ line_data$x)
+  polygon(c(line_data$x, rev(line_data$x)), 
+          c(line_data$upper, rev(line_data$lower)),
+          col = addTrans("#68EBC4",30), border = NA)
+  
+  points(points_data$y ~ points_data$x,
+         pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  
+  axis(side = 1)
+  axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+  mtext(text = var, 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)
+
+for(i in 1:3){
+  
+  eff <- Effect(model, partial.residuals = TRUE, focal.predictors = c("Delta_tc", "cwd_normal_cum"),
+                xlevels = list(Delta_tc = 100),
+                quantiles = c(0.1, 0.5, 0.9))
+  
+  cwd_levels <- eff$variables$cwd_normal_cum$levels
+  
+  line_data <- data.frame(y = inv.as(eff$fit[1:99 + (i-1) * 100]),
+                          lower = inv.as(eff$lower[1:99 + (i-1) * 100]),
+                          upper = inv.as(eff$upper[1:99 + (i-1) * 100]),
+                          x = eff[["x"]][["Delta_tc"]][1:99 + (i-1) * 100]
+  )
+  
+  x_new <- eff[["x"]][["Delta_tc"]]
+  x_orig_all <- eff[["x.all"]]
+  which_x <- x_orig_all$cwd_normal_cum - cwd_levels[i] < 0.001
+  x_orig <- x_orig_all[which_x, ][["Delta_tc"]]
+  
+  y_trans <- eff$fit[1:99 + (i-1) * 100]
+  
+  fitted <- y_trans[closest(x_orig, x_new)] #match original data to the nearest new x value, so
+  # that we can find the closest fitted value and add the residual for that 
+  
+  p_resids <- inv.as(eff$residuals[which_x] + fitted)
+  
+  points_data <- data.frame(x = x_orig,
+                            y = p_resids)
+  
+  plot(NA,
+       ylim = c(0, max(c(y, points_data$y), na.rm = TRUE)),
+       xlim = c(min(x_new), max(x_new)),
+       xlab = "",
+       ylab = "",
+       bg='grey60',
+       col = 'grey30',
+       pch = 21,
+       cex.lab = 1.5,
+       bty = 'n',
+       xaxt = 'n',
+       yaxt = 'n')
+  
+  lines(line_data$y ~ line_data$x, lwd = 2, lty = 1, col = "#1b9e77")
+  lines(line_data$upper ~ line_data$x)
+  lines(line_data$lower ~ line_data$x)
+  polygon(c(line_data$x, rev(line_data$x)), 
+          c(line_data$upper, rev(line_data$lower)),
+          col = addTrans("#68EBC4",30), border = NA)
+  
+  points(points_data$y ~ points_data$x,
+         pch = 21, col = "#1b9e77", bg = "#1b9e77")
+  
+  axis(side = 1)
+  axis(side = 2, at = axTicks(2), labels = axTicks(2)*100)
+  mtext(text = var, 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)
+
+}
+
+dev.off()
+
diff --git a/understory_electivity.R b/understory_electivity.R
index 27f7be7..0d52263 100644
--- a/understory_electivity.R
+++ b/understory_electivity.R
@@ -11,6 +11,7 @@ library(plyr)
 library(vioplot)
 library(reshape2)
 library(multcompView)
+source('addTrans.R', echo=FALSE)
 
 set.seed(16091315)
 
@@ -36,16 +37,6 @@ species$unique_quad <- paste0(species$Plot, species$Transect, species$Meter)
 
 count(species, vars = c("Plot", "Transect"))
 
-
-#little test to look at species relationships with CWD
-# species_aggregate <- aggregate(species[c("Cover")], by = list(species$Plot, species$Species), FUN = mean)
-# names(species_aggregate)[c(1,2)] <- c("Plot", "Species")
-# species_aggregate <- join(species_aggregate, plot_data[c("Plot", "cwd_normal_cum")], by = c("Plot"))
-# 
-# plot(log(Cover) ~ cwd_normal_cum, data = species_aggregate[species_aggregate$Species == "ARTTRIV", ],
-#      xlim = c(200, 450))
-# 
-
 #-----------------------------------------------------------------------
 # Calculate mean cover, cover by quadrat, cover by plot
 #-------------------------------------------------------------------
@@ -78,40 +69,21 @@ 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)
+
+#reduce the number of microsites
 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)))#))
+
+#standardize transect direction names
 ms[ms$Transect == "e", "Transect"] <- "E"
 ms[ms$Transect == "s", "Transect"] <- "S"
 ms[ms$Transect == "w", "Transect"] <- "W"
@@ -119,6 +91,7 @@ ms[ms$Transect == "n", "Transect"] <- "N"
 
 ms$Plot <- as.character(ms$Plot)
 
+#remove the "bonus" quadrats
 ms <- ms[(ms$Transect %in% c("N", "S", "E", "W")), ]
 
 ms[ms$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120"
@@ -126,44 +99,27 @@ ms[ms$Plot == "NPElectricEel120", "Plot"] <- "NPELECTRICEEL120"
 
 ms$unique_quad <- paste0(ms$Plot, ms$Transect, ms$Meter)
 
+
+
+#------------------------------------------------------------------------------
+# Electivity calculations
+#------------------------------------------------------------------------------
 #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)
+## select which plots have enough microsites of each type
+# for which the cover is enough to have 1% cover in each microsite type,
+# if the cover were uniformly distributed
+types <- c("All", "Perennial grass", "Cheatgrass", "Perennial forb", "Shrub")
+ntypes <- length(types)
 
-plots_to_use <- unique(ms$Plot)[plots_to_use_dead & plots_to_use_live]
+plots_to_use <- list(NA, NA, NA, NA, NA)
+names(plots_to_use) <- types
 
+n_plots <- 1
+type <- "Cheatgrass"  
 
-#------------------------------------------------------------------------------
-# 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){
@@ -175,51 +131,66 @@ for(i in 1:ntypes){
                                    DI = numeric(0),
                                    LI = numeric(0))
 }
+
 names(elect_results) <- types
 
 
-for(type in types){ # calculate electivity for each functional type
+for(type in types){
   
-  for (i in 1:length(plots_to_use)){
-    
+  for(j in 1:length(unique(ms$Plot))){
     if(type == "All"){
-      daub_cov <- daub[daub$Plot == as.character(plots_to_use[i]), ]
+      daub_cov <- daub[daub$Plot == unique(ms$Plot)[j], ]
       daub_cov <- daub_cov[daub_cov$Cover.type %in% c("Cheatgrass", "Other Ann Grass", 
-                        "Perennial grass", "Annual forb", "Perennial forb", "Shrub"), ]
+                                                      "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 <- daub[daub$Cover.type == type & daub$Plot == unique(ms$Plot)[j], ]
       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)
+    total_cov <- sum(daub_cov$Midpoint.value)
+    uniform_cov <- total_cov/20 #divy up cover evenly among quadrats
+    
+    table_ms <- table(daub_cov$ms)
+    
+    table_ms_uniform <- table_ms*uniform_cov #multiple average quadrat cover by microsite 
+                                             # amount to generate "expected" cover
     
-    cov <- aggregate(daub_cov[, "prop_cov"], by = list(daub_cov$ms), FUN = sum)
+    use_plot <- sum(table_ms_uniform >= 1) == 3 #are there greater than or equal to 1% cover
+                                               #in each microsite type?
     
-    prev <- table(daub_cov$ms)/20
+    plots_to_use[[type]][j] <- use_plot #save plot condition to global list
+    
+    if(use_plot){
+    
+      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 = unique(ms$Plot)[j],
+                          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 = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Live",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Live",2])
+      elect$DI = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Inter",2])
+      elect$LI = (cov[cov$Group.1 == "Live",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Live",2] + cov[cov$Group.1 == "Inter",2])
+      
+      elect_results[[type]] <- rbind(elect_results[[type]], elect) #I know this is slow
+    }
     
-    #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 = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Live",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Live",2])
-    elect$DI = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Inter",2])
-    elect$LI = (cov[cov$Group.1 == "Live",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Live",2] + cov[cov$Group.1 == "Inter",2])
     
-    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
@@ -237,6 +208,7 @@ for(i in 1:ntypes){
                                  DI = numeric(0),
                                  LI = numeric(0))
 }
+
 names(elect_means) <- types
 
 system.time(#takes about a half hour
@@ -246,40 +218,46 @@ system.time(#takes about a half hour
 for(type in types){
   
   plot_to_use_type <- subset(elect_results[[type]], !is.na(Dead))$Plot
+  nplots <- length(plot_to_use_type)
   
   for(j in 1:nit){
     
     #temporary storage for each iteration
-    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){
+    elect_rand <- data.frame(Plot = character(nplots),
+                             Dead = numeric(nplots),
+                             Live = numeric(nplots),
+                             Inter = numeric(nplots),
+                             DL = numeric(nplots),
+                             DI = numeric(nplots),
+                             LI = numeric(nplots), stringsAsFactors = FALSE)
+ 
+    for (k in 1:nplots){
+      # calculate electivity for a each plot using randomized cover values
+      plot_select <- plot_to_use_type[k]
+      
       if(type == "All"){
-        daub_cov <- daub[daub$Plot == as.character(plot_use), ]
+        daub_cov <- daub[daub$Plot == as.character(plot_select), ]
         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 <- daub[daub$Cover.type == type & daub$Plot == as.character(plot_select), ]
         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)
+      
+      #randomize the cover values by quadrat
       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
       
       #temporary df to store results from each plot
-      elect <- data.frame(Plot = plot_use,
+      elect <- data.frame(Plot = plot_select,
                           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")]) / 
@@ -290,14 +268,13 @@ for(type in types){
       elect$DI = (cov[cov$Group.1 == "Dead",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Dead",2] + cov[cov$Group.1 == "Inter",2])
       elect$LI = (cov[cov$Group.1 == "Live",2] - cov[cov$Group.1 == "Inter",2]) / (cov[cov$Group.1 == "Live",2] + cov[cov$Group.1 == "Inter",2])
       
-      elect_rand <- rbind(elect_rand, elect) #add data from plot to temporary df
+      elect_rand[k, ] <- elect #add data from plot to temporary df
     }
     
     # add mean from electivity from this iteration to global df
     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
-    
   
   }
   
@@ -336,7 +313,7 @@ results_boots[[i]][2, ] <- pvals
 }
 
 #save and read results depending on what you need
-# saveRDS(results_boots, paste0("./outputs/results_boots_", n_plots, "dead.rds"))
+saveRDS(results_boots, paste0("./outputs/results_boots_", n_plots, "dead.rds"))
 
 # results_boots <- readRDS("./outputs/results_boots_1dead.rds")
 
@@ -344,8 +321,8 @@ results_boots[[i]][2, ] <- pvals
 #-----------------------------------------------------------------------------------------------
 # Summary table
 #-----------------------------------------------------------------------------------------------
-summary <- as.data.frame(matrix(nrow = 12, ncol = 8))
-for(i in 1:4){
+summary <- as.data.frame(matrix(nrow = 15, ncol = 8))
+for(i in 1:5){
   for(j in 1:6){
     summary[3*i - 2, j+2] <- results_boots[[i]][[j]][[1]]
     summary[3*i - 1, j+2] <- results_boots[[i]][[j]][[1]] - mean(elect_means[[i]][[j + 1]])
@@ -353,8 +330,8 @@ for(i in 1:4){
   }
 }
 
-summary[, 2] <- rep(c("Empirical", "Difference", "p-val"), times = 4)
-summary[, 1] <- rep(c("Perennial grass", "Cheatgrass", "Perennial forb", "Shrub"), each = 3)
+summary[, 2] <- rep(c("Empirical", "Difference", "p-val"), times = 5)
+summary[, 1] <- rep(types, each = 3)
 names(summary) <- c("FT", "var", "Dead", "Live", "Inter", "D-L", "D-I", "L-I")
 write.csv(summary, "./outputs/electivity summary table.csv")
 
@@ -363,7 +340,7 @@ write.csv(summary, "./outputs/electivity summary table.csv")
 # Figure 5
 #-----------------------------------------------------------------------------------------------
 
-tiff(filename="./outputs/Figure_5_electivity.tiff", 
+tiff(filename="./outputs/Figure_5_electivity_test.tiff", 
     type="cairo",
     units="in", 
     width = 4, 
@@ -376,7 +353,7 @@ par(mfrow = c(2,2),
     mar = c(2,1,1,2),
     oma = c(2,3,1,0))
 
-for(i in 1:ntypes){
+for(i in 2:5){
 
   melt_elect <- melt(elect_results[[i]][, 2:4])
   plot(NA,
@@ -393,22 +370,24 @@ for(i in 1:ntypes){
   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)
 
+  text(x = 3, y = 0.7, labels = paste0("n = ", nrow(melt_elect)/3))
   
-  if(i %in% c(1,2)){
+  if(i %in% c(2,3)){
     axis(1, at = c(1,2,3), labels = FALSE)
   }
   
-  if(i %in% c(3,4)){
+  if(i %in% c(4,5)){
   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 = paste0("(", letters[i-1], ")"), outer = FALSE, side = 3, at = 0.5, line = 0.3)
   mtext(text = "Ivlev's E", outer = TRUE, side = 2, line = 1.5)