EDDI_productivity.R

library(car)
library(MASS)
library(ggplot2)
library(viridisLite)

setwd("/Users/marketa/Library/Mobile Documents/com~apple~CloudDocs/School/Thesis/EDDI")
EDDI <- read.csv("EDDI_3mn_all.csv", header = TRUE)
EDDI <- EDDI[,-1]
EDDI <- EDDI*-1
names(EDDI) <- c("2000", "2001", "2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012", "2013","2014", "2015", "2016")
Productivity <- read.csv("Prod_anomaly.csv", header = TRUE)
names(Productivity) <- c("2000", "2001", "2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012", "2013","2014", "2015", "2016")

#stack all columns into one and save headers as a factor
EDDI_all <- stack(lapply(EDDI, as.character))
#set column "values" as numeric value
EDDI_all$values <- as.numeric(EDDI_all$values)
#ad an ID column to be able to merge data frames
EDDI_all$ID <- seq.int(nrow(EDDI_all))
#create SD and mean variables
EDDI_sd <- sd(EDDI_all$values)*sqrt((length(EDDI_all$values)-1)/(length(EDDI_all$values)))
EDDI_mean <- mean(EDDI_all$values)                                    
#create z-score column
EDDI_all$zscore <- (EDDI_all$values-EDDI_mean)/EDDI_sd

Productivity_all <- stack(lapply(Productivity, as.character))
Productivity_all$values <- as.numeric(Productivity_all$values)
Productivity_all$ID <- seq.int(nrow(Productivity_all))

#merge columns togther so you can do linear model
Prod_EDDI <- merge(EDDI_all, Productivity_all, by = "ID")
colnames(Prod_EDDI) <- c("ID", "EDDI", "year","z_score", "Prod", "year2")
prod.lm = lm(Productivity_all$values~EDDI_all$zscore)
summary(prod.lm)$r.squared

#Setting the colors for the plot
color_palette_function <- colorRampPalette(colors = c("darkred", "orange", "yellow", "lightgreen", "darkblue"), space = "Lab")
numcolors <- nlevels(Prod_EDDI$year)
final_colors <- color_palette_function(numcolors)

#colorful plot based on a factor=year
plot(Prod_EDDI$z_score, Prod_EDDI$Prod, pch = 18, col = (final_colors[Prod_EDDI$year]), main = "Productivity anomaly ~ 3-month EDDI\n r2 = 0.383", xlab = "3-month EDDI index (z-score)", ylab = "Productivity anomaly [kg ha-1 year-1]", xlim = c(-4,3))
legend(x="bottomright", legend = paste(levels(Prod_EDDI$year)), col = final_colors, pch = 19, cex = .5 )
#inserting lines in the plot
abline(h=0, v=0, col = "darkgray", lty = 2)
abline(lm(Productivity_all$values~EDDI_all$zscore), col = "red")

#calculating a density of points in a 100-part grid
theme_set(theme_bw(base_size = 16))

get_density <- function(x, y, n = 100){
  dens <- MASS::kde2d(x = x, y = y, n = n)
  ix <- findInterval(x, dens$x)
  iy <- findInterval(y, dens$y)
  ii <- cbind(ix, iy)
  return(dens$z[ii])
}

#creating a density plot
Prod_EDDI$density <- get_density(Prod_EDDI$z_score, Prod_EDDI$Prod)

p <- ggplot(Prod_EDDI) + geom_point(aes(z_score, Prod, color = density)) + scale_color_gradientn(colors = terrain.colors(10) , breaks = c(2e-05, 1.15e-03), labels = c("Low", "High"))
p <- p + geom_smooth(data = Prod_EDDI, aes(x=z_score, y=Prod), method = "lm", formula = y~x)

#title, subtitle, labels, title of the legend
p <- p + xlim(-4, 3)
p <- p + labs(title = "Productivity anomaly ~ 3-month EDDI index", subtitle = "Density of points, r2 = 0.383", x = "3-month EDDI index (z-score)", y = "Productivity anomaly [kg ha-1 yr-1]", color = "Density")

#black background
p + theme(plot.title = element_text(size = 15), plot.subtitle = element_text(color = "darkgrey"), panel.background = element_rect(fill = "black", colour = "black", size = 0.5, linetype = "solid"))