SPI6m_Prod_final.R

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

setwd("/Users/marketa/Library/Mobile Documents/com~apple~CloudDocs/School/Thesis/SPI")
SPI6m <- read.csv("SPI6m.csv", header = TRUE)
names(SPI6m) <- 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
SPI6m_all <- stack(lapply(SPI6m, as.character))
#set column "values" as numeric value
SPI6m_all$values <- as.numeric(SPI6m_all$values)
#ad an ID column to be able to merge data frames
SPI6m_all$ID <- seq.int(nrow(SPI6m_all))
#create SD and mean variables
SPI6m_sd <- sd(SPI6m_all$values)*sqrt((length(SPI6m_all$values)-1)/(length(SPI6m_all$values)))
SPI6m_mean <- mean(SPI6m_all$values)                                    
#create z-score column
SPI6m_all$zscore <- (SPI6m_all$values-SPI6m_mean)/SPI6m_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_SPI6m <- merge(SPI6m_all, Productivity_all, by = "ID")
colnames(Prod_SPI6m) <- c("ID", "SPI6m", "year","z_score", "Prod", "year2")
prod.lm = lm(Productivity_all$values~SPI6m_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_SPI6m$year)
final_colors <- color_palette_function(numcolors)

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

p <- ggplot(Prod_SPI6m) + 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_SPI6m, 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 ~ 6-month SPI index", subtitle = "Density of points, r2 = 0.663", x = "SPI6m 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"))