Created
November 6, 2018 21:41
-
-
Save MarketaP/aefad3e07e1aab5ae5ca2b722b92eece to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(car) | |
library(MASS) | |
library(ggplot2) | |
library(viridisLite) | |
setwd("C:/Users/mpodebradska2/Desktop/School/Thesis/Indices_analysis/VegDRI") | |
VegDRI <- read.csv("VegDRI.csv", header = TRUE) | |
names(VegDRI) <- c( "2009", "2010", "2011", "2012", "2013","2014", "2015", "2016") | |
Productivity <- read.csv("prod_anomaly.csv", header = TRUE) | |
Productivity <- Productivity[-480,-c(1,2,3,4,5,6,7,8,9)] | |
names(Productivity) <- c("2009", "2010", "2011", "2012", "2013","2014", "2015", "2016") | |
#stack all columns into one and save headers as a factor | |
VegDRI_all <- stack(lapply(VegDRI, as.character)) | |
#set column "values" as numeric value | |
VegDRI_all$values <- as.numeric(VegDRI_all$values) | |
#ad an ID column to be able to merge data frames | |
VegDRI_all$ID <- seq.int(nrow(VegDRI_all)) | |
#create SD and mean variables | |
VegDRI_sd <- sd(VegDRI_all$values)*sqrt((length(VegDRI_all$values)-1)/(length(VegDRI_all$values))) | |
VegDRI_mean <- mean(VegDRI_all$values) | |
#create z-score column | |
VegDRI_all$zscore <- (VegDRI_all$values-VegDRI_mean)/VegDRI_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_VegDRI <- merge(VegDRI_all, Productivity_all, by = "ID") | |
colnames(Prod_VegDRI) <- c("ID", "VegDRI", "year","z_score", "Prod", "year2") | |
prod.lm = lm(Productivity_all$values~VegDRI_all$zscore) | |
summary(prod.lm)$r.squared | |
prod.lm1 = lm(Productivity_all$values~VegDRI_all$values) | |
summary(prod.lm1)$r.squared | |
#Setting the colors for the plot | |
color_palette_function <- colorRampPalette(colors = c("darkred", "orange", "yellow", "lightgreen", "darkblue"), space = "Lab") | |
numcolors <- nlevels(Prod_VegDRI$year) | |
final_colors <- color_palette_function(numcolors) | |
#colorful plot based on a factor=year | |
plot(Prod_VegDRI$z_score, Prod_VegDRI$Prod, pch = 18, col = (final_colors[Prod_VegDRI$year]), main = "Productivity anomaly ~ VegDRI\n r2 = 0.581", xlab = "VegDRI index (z-score)", ylab = "Productivity anomaly [kg ha-1 year-1]") | |
legend(x="bottomright", legend = paste(levels(Prod_VegDRI$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~VegDRI_all$zscore), col = "red") | |
#colorful plot based on a factor=year | |
plot(Prod_VegDRI$VegDRI, Prod_VegDRI$Prod, pch = 18, col = (final_colors[Prod_VegDRI$year]), main = "Productivity anomaly ~ VegDRI\n r2 = 0.581", xlab = "VegDRI index", ylab = "Productivity anomaly [kg ha-1 year-1]") | |
legend(x="bottomright", legend = paste(levels(Prod_VegDRI$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~VegDRI_all$values), 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_VegDRI$density <- get_density(Prod_VegDRI$z_score, Prod_VegDRI$Prod) | |
p <- ggplot(Prod_VegDRI) + 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_VegDRI, 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 ~ VegDRI index", subtitle = "Density of points, r2 = 0.581", x = "VegDRI 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")) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment