ESI_productivity_final.R

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

setwd("/Users/marketa/Library/Mobile Documents/com~apple~CloudDocs/School/Thesis/ESI")
ESI <- read.csv("ESI.csv", header = TRUE)
names(ESI) <- c("2001", "2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012", "2013","2014", "2015", "2016")
Productivity <- read.csv("prod_anomaly.csv", header = TRUE)
Productivity <- Productivity[,-1]
Productivity <- Productivity[-c(594, 595, 601, 616),]
names(Productivity) <- c("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
ESI_all <- stack(lapply(ESI, as.character))
#set column "values" as numeric value
ESI_all$values <- as.numeric(ESI_all$values)
#ad an ID column to be able to merge data frames
ESI_all$ID <- seq.int(nrow(ESI_all))
#create SD and mean variables
ESI_sd <- sd(ESI_all$values)*sqrt((length(ESI_all$values)-1)/(length(ESI_all$values)))
ESI_mean <- mean(ESI_all$values)                                    
#create z-score column
ESI_all$zscore <- (ESI_all$values-ESI_mean)/ESI_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_ESI <- merge(ESI_all, Productivity_all, by = "ID")
colnames(Prod_ESI) <- c("ID", "ESI", "year","z_score", "Prod", "year2")
prod.lm = lm(Productivity_all$values~ESI_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_ESI$year)
final_colors <- color_palette_function(numcolors)

#colorful plot based on a factor=year
plot(Prod_ESI$z_score, Prod_ESI$Prod, pch = 18, col = (final_colors[Prod_ESI$year]), main = "Productivity anomaly ~ 3-month ESI\n r2 = 0.661", xlab = "3-month ESI index (z-score)", ylab = "Productivity anomaly [kg ha-1 year-1]", xlim = c(-4,3))
legend(x="bottomright", legend = paste(levels(Prod_ESI$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~ESI_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_ESI$density <- get_density(Prod_ESI$z_score, Prod_ESI$Prod)

p <- ggplot(Prod_ESI) + geom_point(aes(z_score, Prod, color = density)) + scale_color_gradientn(colors = terrain.colors(10) , breaks = c(2e-05, 1.3e-03), labels = c("Low", "High"))
p <- p + geom_smooth(data = Prod_ESI, 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 ESI index", subtitle = "Density of points, r2 = 0.661", x = "3-month ESI 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"))