rm(list=ls())# wipes slate clean library(mgcv) library(lme4) library(lmerTest) library(ggplot2) library(dplyr) library( geosphere ) library( stringr) setwd("C:/Users/Ragenaky/Desktop/Thesis chapter 3/Data/Master Sheets") #Calculate a standard error stderr <- function(x, ...) sd(x, na.rm = TRUE) / sqrt(length(is.na(x == FALSE)) ) ### Install this When you start for Multiplots!!!##### # # ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) # - cols: Number of columns in layout # - layout: A matrix specifying the layout. If present, 'cols' is ignored. # # If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), # then plot 1 will go in the upper left, 2 will go in the upper right, and # 3 will go all the way across the bottom. # multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { library(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } #Fig 4a Open_Data_IC_RC_WD<-read.csv("Open_Data_IC_RC_WD.CSV") # Mean rainfall model1 <- lm( log( Control_Mean + 1) ~ Mean_RF, data = Open_Data_IC_RC_WD ) anova(model1) summary(model1) Open_Data_IC_WDA <- Open_Data_IC_RC_WD Open_Data_IC_WDA$rainCat <- round(Open_Data_IC_WDA$ Mean_RF / 1.5) * 1.5 summaryRain <- Open_Data_IC_WDA %>% group_by( rainCat ) %>% summarise( meanN = mean(log( Control_Mean + 1), na.rm = TRUE), SE = stderr(log( Control_Mean + 1), na.rm = TRUE) ) fig4a <- ggplot( summaryRain,aes(x = rainCat, y = meanN) ) + geom_point(size = 1) + geom_errorbar(aes( ymin = meanN - SE, ymax = meanN + SE), width = 0.5, size = 0.25 ) + theme_bw() + theme( panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = 'black', size = 0.25), axis.ticks = element_line(colour = "black", size = 0.25), axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), legend.position="none", axis.title.x=element_text( size = 12 ), axis.title.y=element_text( size = 12 ) ) + labs( x = "Mean rainfall (mm)", y = "Log Weed density") + theme(axis.text.x = element_text(angle = 90)) fig4a # Precipitation seasonality model2 <- lm( log( Control_Mean + 1) ~ RFCV, data = Open_Data_IC_RC_WD ) anova(model2) summary(model2) Open_Data_IC_WDA <- Open_Data_IC_RC_WD Open_Data_IC_WDA$RFCVCat <- round(Open_Data_IC_WDA$ RFCV / 1.5) * 1.5 summaryRFCV <- Open_Data_IC_WDA %>% group_by( RFCVCat ) %>% summarise( meanN = mean(log( Control_Mean + 1), na.rm = TRUE), SE = stderr(log( Control_Mean + 1), na.rm = TRUE) ) fig4b <- ggplot(summaryRFCV, aes(x = RFCVCat, y = meanN) ) + geom_point(size = 1) + geom_errorbar(aes( ymin = meanN - SE, ymax = meanN + SE), width = 0.5, size = 0.25 ) + theme_bw() + theme( panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = 'black', size = 0.25), axis.ticks = element_line(colour = "black", size = 0.25), axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), legend.position="none", axis.title.x=element_text( size = 12 ), axis.title.y=element_text( size = 12 ) ) + labs( x = "Precipitation seasonality (CV)", y = "Log Weed density") + theme(axis.text.x = element_text(angle = 90)) fig4b # ------------------------------ # altitude model3 <- lm( log( Control_Mean + 1) ~ Alt, data = Open_Data_IC_RC_WD) anova(model3) summary(model3) Open_Data_IC_WDA <- Open_Data_IC_RC_WD Open_Data_IC_WDA$altCat <- round(Open_Data_IC_WDA$Alt / 100) * 100 summaryAlt <- Open_Data_IC_WDA %>% group_by( altCat ) %>% summarise( meanN = mean(log( Control_Mean + 1), na.rm = TRUE), SE = stderr(log( Control_Mean + 1), na.rm = TRUE) ) fig4c <- ggplot(summaryAlt, aes(x = altCat, y = meanN) ) + geom_point(size = 1) + geom_errorbar(aes( ymin = meanN - SE, ymax = meanN + SE), width = 0.5, size = 0.25 ) + theme_bw() + theme( panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = 'black', size = 0.25), axis.ticks = element_line(colour = "black", size = 0.25), axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), legend.position="none", axis.title.x=element_text( size = 12 ), axis.title.y=element_text( size = 12 ) ) + labs( x = "Altitude (m)", y = "Log Weed density") + theme(axis.text.x = element_text(angle = 90)) fig4c # Mean temperature model4 <- lm( log( Control_Mean + 1) ~ Mean_TA, data = Open_Data_IC_RC_WD) anova(model4) summary(model4) Open_Data_IC_WDA <- Open_Data_IC_RC_WD Open_Data_IC_WDA$tempCat <- round(Open_Data_IC_WDA$Mean_TA / 1) * 1 summaryTemp <- Open_Data_IC_WDA %>% group_by( tempCat ) %>% summarise( meanN = mean (log( Control_Mean + 1), na.rm = TRUE), SE = stderr(log( Control_Mean + 1), na.rm = TRUE) ) fig4d <- ggplot(summaryTemp, aes(x = tempCat, y = meanN) ) + geom_point(size = 1) + geom_errorbar(aes( ymin = meanN - SE, ymax = meanN + SE), width = 0.5, size = 0.25 ) + theme_bw() + theme( panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = 'black', size = 0.25), axis.ticks = element_line(colour = "black", size = 0.25), axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm"), size = 10), legend.position="none", axis.title.x=element_text( size = 12 ), axis.title.y=element_text( size = 12 ) ) + labs( x = "Mean Temperature (\u00B0C)", y = "Log Weed Density") + theme(axis.text.x = element_text(angle = 90)) fig4d multiplot(fig4a + labs( tag = "A"), fig4b+ labs( tag = "B"), fig4c+ labs( tag = "C"), fig4d+ labs( tag = "D"), cols = 2)