citation() getwd() meta<-read.csv("salinity.C.flux.090521.csv") head(meta) daily<-read.csv("salinity.C.flux.daily.060321.csv") head(daily) monthly<-read.csv("salinity.C.flux.monthly.060321.csv") head(monthly) #Subset Data exp.meta<-meta[meta$Type %in% c("manipulation"),] exp.meta gradient.meta<-meta[meta$Type %in% c("salinity gradient"),] gradient.meta library(ISLR) library(mgcv) library(voxel) library(tidyverse) library(gridExtra) library(car) library(mgcViz) library("grid") library("ggplotify") library(ggplot2) library(dplyr) library(tidymv) library(boot) library(ggplot2) theme_set(theme_bw()) library(sf) library(rnaturalearth) library(rnaturalearthdata) world <- ne_countries(scale = "medium", returnclass = "sf") class(world) par(mar=c(5,5,5,5)) ggplot(data = world) + geom_sf() ggplot(data = world) + geom_sf() + coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) + geom_point(data=meta, aes(x=Longitude, y=Latitude), cex=3, col="blue") ggplot(data = world) + geom_sf() + xlab("Longitude") + ylab("Latitude") + ggtitle("World map", subtitle = paste0("(", length(unique(world$NAME)), " countries)")) ggplot(data = world) + geom_sf() + coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ") library(ggspatial) ggplot(data = world) + geom_sf() + annotation_scale(location = "bl", width_hint = 0.5) + annotation_north_arrow(location = "bl", which_north = "true", pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"), style = north_arrow_fancy_orienteering) + coord_sf(xlim = c(-180, 180), ylim = c(-100, 100)) ###Salinity C Flux Plots library(ggplot2) library(gridExtra) library(ggpmisc) q1<-ggplot(meta, aes(Salinity, GEPmedian)) + ylim(0, 40) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) q2<-ggplot(meta, aes(Salinity, ERmedian)) + ylim(0, 10) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) q3<-ggplot(meta, aes(Salinity, CH4median)) + ylim(0, 1) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) q4<-ggplot(meta, aes(Salinity, NEPmedian)) + ylim(0, 40) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) grid.arrange(q1, q2, q3, q4, nrow=1, ncol=4) p1 <- ggplot(data = meta, aes(x = Salinity.change, y = GEPmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p1 q1 <- p1 + stat_smooth(method = "gam", formula = y ~ s(x, k = 12), size = 1) q1 p2 <- ggplot(data = meta, aes(x = Salinity.change, y = ERmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p2 q2 <- p2 + stat_smooth(method = "gam", formula = y ~ s(x, k = 12), size = 1) q2 p3 <- ggplot(data = meta, aes(x = Salinity.change, y = CH4mediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p3 q3 <- p3 + stat_smooth(method = "gam", formula = y ~ s(x, k = 13), size = 1) q3 p4 <- ggplot(data = meta, aes(x = Salinity.change, y = NEPmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p4 q4 <- p4 + stat_smooth(method = "gam", formula = y ~ s(x, k = 15), size = 1) q4 grid.arrange(q1, q2, q3, q4, nrow=1, ncol=4) ###Experimental Salinity p1 <- ggplot(data = exp.meta, aes(x = Salinity, y = GEPmedian)) + ylim(0, 40) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p1 q1 <- p1 + stat_smooth(method = "gam", formula = y ~ s(x, k = 15), size = 1) q1 p1 <- ggplot(data = exp.meta, aes(x = Salinity, y = GEPmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p1 q1 <- p1 + stat_smooth(method = "gam", formula = y ~ s(x, k = 6), size = 1) q1 p2 <- ggplot(data = exp.meta, aes(x = Salinity, y = ERmedian)) + ylim(0, 10) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p2 q2 <- p2 + stat_smooth(method = "gam", formula = y ~ s(x, k = 17), size = 1) q2 p2 <- ggplot(data = exp.meta, aes(x = Salinity, y = ERmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p2 q2 <- p2 + stat_smooth(method = "gam", formula = y ~ s(x, k = 9), size = 1) q2 p3 <- ggplot(data = exp.meta, aes(x = Salinity, y = CH4median)) + ylim(-1, 1) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p3 q3 <- p3 + stat_smooth(method = "gam", formula = y ~ s(x, k = 12), size = 1) q3 p3 <- ggplot(data = exp.meta, aes(x = Salinity, y = CH4mediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p3 q3 <- p3 + stat_smooth(method = "gam", formula = y ~ s(x, k = 5), size = 1) q3 p4 <- ggplot(data = exp.meta, aes(x = Salinity, y = NEPmedian)) + ylim(-1, 22) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p4 q4 <- p4 + stat_smooth(method = "gam", formula = y ~ s(x, k = 14), size = 1) q4 p4 <- ggplot(data = exp.meta, aes(x = Salinity, y = NEPmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p4 q4 <- p4 + stat_smooth(method = "gam", formula = y ~ s(x, k = 8), size = 1) q4 ####Salinity Gradients p5 <- ggplot(data = gradient.meta, aes(x = Salinity, y = GEPmedian)) + ylim(0, 40) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p5 q5 <- p5 + stat_smooth(method = "gam", formula = y ~ s(x, k = 8), size = 1) q5 p5 <- ggplot(data = gradient.meta, aes(x = Salinity, y = GEPmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p5 q5 <- p5 + stat_smooth(method = "gam", formula = y ~ s(x, k = 5), size = 1) q5 p6 <- ggplot(data = gradient.meta, aes(x = Salinity, y = ERmedian)) + ylim(0, 10) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p6 q6 <- p6 + stat_smooth(method = "gam", formula = y ~ s(x, k = 10), size = 1) q6 p6 <- ggplot(data = gradient.meta, aes(x = Salinity, y = ERmediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p6 q6 <- p6 + stat_smooth(method = "gam", formula = y ~ s(x, k = 5), size = 1) q6 p7 <- ggplot(data = gradient.meta, aes(x = Salinity, y = CH4median)) + ylim(-1, 1) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p7 q7 <- p7 + stat_smooth(method = "gam", formula = y ~ s(x, k = 18), size = 1) q7 p7 <- ggplot(data = gradient.meta, aes(x = Salinity, y = CH4mediandifference)) + ylim(-150, 150) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p7 q7 <- p7 + stat_smooth(method = "gam", formula = y ~ s(x, k = 12), size = 1) q7 p8 <- ggplot(data = gradient.meta, aes(x = Salinity, y = NEPmedian)) + ylim(-1, 22) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p8 q8 <- p8 + stat_smooth(method = "gam", formula = y ~ s(x, k = 9), size = 1) q8 p8 <- ggplot(data = gradient.meta, aes(x = Salinity, y = NEPmediandifference)) + ylim(-250, 250) + xlim(0, 35) + geom_point(shape=21, fill="light grey", color="black", size=4) + theme(text = element_text(size=15)) p8 q8 <- p8 + stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) q8 grid.arrange(q1, q2, q3, q4, q5, q6, q7, q8, nrow=2, ncol=4) # correlation plots: library(ggplot2) library(GGally) citation("GGally") my_fn <- function(data, mapping, method="loess", ...){ p <- ggplot(data = data, mapping = mapping) + geom_point() + geom_smooth(method=method, ...) p } png("AM_Correlations_all.daily.png", width = 900, height = 500) ggpairs(with(daily, data.frame(Salinity, GEPmedian, ERmedian, CH4median, NEPmedian)), lower = list(continuous = my_fn)) dev.off() warnings() png("AM_Correlations_all.monthly.png", width = 900, height = 500) ggpairs(with(monthly, data.frame(Salinity, GEPmedian, ERmedian, CH4median, NEPmedian)), lower = list(continuous = my_fn)) dev.off() warnings() # GAM Models https://m-clark.github.io/generalized-additive-models/application.html m.gep <- gam( GEP ~ factor(Salinity) + s(temp.mean, bs="cr", k=10) + s(par.sum, bs="cr", k=5) +s( by=factor(Salinity), k=30, bs="cr") , data=data, method = 'REML', correlation = corARMA(form=~1| Date, p=5, q=30)) gam.check(m.gep) # you want all large p-values summary(m.gep) acf( m.gep$residuals) pacf( m.gep$residuals) library(mgcv) #All Data gep <- gam(GEPmedian ~ s(Salinity), data=meta, method = "REML") summary(gep) er <- gam(ERmedian ~ s(Salinity), data=meta, method = "REML") summary(er) methane <- gam(CH4median ~ s(Salinity), data=meta, method = "REML") summary(methane) nep <- gam(NEPmedian ~ s(Salinity), data=meta, method = "REML") summary(nep) #Experimental Manipulations gep <- gam(GEPmediandifference ~ s(Salinity), data=exp.meta, method = "REML") summary(gep) er <- gam(ERmediandifference ~ s(Salinity), data=exp.meta, method = "REML") summary(er) methane <- gam(CH4mediandifference ~ s(Salinity), data=exp.meta, method = "REML") summary(methane) nep <- gam(NEPmediandifference ~ s(Salinity), data=exp.meta, method = "REML") summary(nep) #Salinity Gradients gep <- gam(GEPmediandifference ~ s(Salinity), data=gradient.meta, method = "REML") summary(gep) er <- gam(ERmediandifference ~ s(Salinity), data=gradient.meta, method = "REML") summary(er) methane <- gam(CH4mediandifference ~ s(Salinity), data=gradient.meta, method = "REML") summary(methane) nep <- gam(NEPmediandifference ~ s(Salinity), data=gradient.meta, method = "REML") summary(nep) #All Data Differences & Salinity Changes gep <- gam(GEPmediandifference ~ s(Salinity.change), data=meta, method = "REML") summary(gep) er <- gam(ERmediandifference ~ s(Salinity.change), data=meta, method = "REML") summary(er) methane <- gam(CH4mediandifference ~ s(Salinity.change), data=meta, method = "REML") summary(methane) nep <- gam(NEPmediandifference ~ s(Salinity.change), data=meta, method = "REML") summary(nep)