glasshouse_dom.R
################################################################################
#### Project: Lowland plant migrations alpine soil C loss
#### Title: Glasshouse pot data
#### Author: Tom Walker (thomas.walker@usys.ethz.ch)
#### Date: 29 November 2021
#### ---------------------------------------------------------------------------
#### PROLOGUE ------------------------------------------------------------------
## Options ----
# remove objects from global environment
rm(list = ls())
# R session options (no factors, bias against scientific #s)
options(
stringsAsFactors = F,
scipen = 6
)
## Libraries ----
# standard library set
library(nlme)
library(emmeans)
library(tidyverse)
# custom functions
source("./r_code/mm2in.R")
#### DATA ----------------------------------------------------------------------
## Load ----
soil <- drake::readd(gh_soil)
## Format ----
soilPools <- soil$pot_soil
#### FORMAT FOR PLOTS ----------------------------------------------------------
## Soil pools ----
sumSoil <- soilPools %>%
# recalculate Cmic in mg g
mutate(k = k * 24) %>%
mutate(p = p/1000) %>%
# group by treatment, calculate means/se, ungroup
group_by(treatment) %>%
summarise(a350_mean = mean(DOM.a350, na.rm = T),
a350_se = sd(DOM.a350, na.rm = T)/sqrt(n()),
ftot_mean = mean(DOM.Ftot, na.rm = T),
ftot_se = sd(DOM.Ftot, na.rm = T)/sqrt(n()),
fi_mean = mean(DOM.FI, na.rm = T),
fi_se = sd(DOM.FI, na.rm = T)/sqrt(n()),
dom1_mean = mean(DOM.C1, na.rm = T),
dom1_se = sd(DOM.C1, na.rm = T)/sqrt(n()),
dom2_mean = mean(DOM.C2, na.rm = T),
dom2_se = sd(DOM.C2, na.rm = T)/sqrt(n()),
dom3_mean = mean(DOM.C3, na.rm = T),
dom3_se = sd(DOM.C3, na.rm = T)/sqrt(n()),
dom4_mean = mean(DOM.C4, na.rm = T),
dom4_se = sd(DOM.C4, na.rm = T)/sqrt(n()),
dom5_mean = mean(DOM.C5, na.rm = T),
dom5_se = sd(DOM.C5, na.rm = T)/sqrt(n()),
dom6_mean = mean(DOM.C6, na.rm = T),
dom6_se = sd(DOM.C6, na.rm = T)/sqrt(n())) %>%
ungroup
## get bare lines ----
allBare <- filter(sumSoil, treatment == "B")
sumSoil <- filter(sumSoil, treatment != "B")
#### SOIL BARE TREATMENT BOUNDS ------------------------------------------------
## Calculate ribbon for bare treatment on plots ----
soilPools %>%
filter(treatment == "B") %>%
dplyr::select(Cmic.ugC.g:b) %>%
# calculate upper and lower bounds for all variables
summarise(
across(
everything(),
function(x){
# summary statistics
mu <- mean(x, na.rm = T)
n <- sum(!is.na(x))
sdev <- sd(x, na.rm = T)
se <- sdev/sqrt(n)
# calculate upper and lower confidence intervals
high <- mu + se
low <- mu - se
out <- c(high, low)
return(out)
}
)
)
## Remove bare treatment (only relevant as reference) ----
soilPools <- soilPools %>%
filter(treatment != "B")
#### ANALYSE -------------------------------------------------------------------
## DOM a350 ----
# build model
m1 <- lme(
DOM.a350 ~ treatment,
random = ~ 1 | block,
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r1 <- residuals(m1, type = "normalized")
par(mfrow = c(1, 3))
plot(r1 ~ fitted(m1))
boxplot(r1 ~ soilPools$treatment)
hist(r1)
# test main effects
m1a <- update(m1, ~.- treatment)
anova(m1, m1a)
## DOM total fluorescence ----
# build model
m2 <- lme(
DOM.Ftot ~ treatment,
random = ~ 1 | block,
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r2 <- residuals(m2, type = "normalized")
par(mfrow = c(1, 3))
plot(r2 ~ fitted(m2))
boxplot(r2 ~ soilPools$treatment)
hist(r2)
# test main effects
m2a <- update(m2, ~.- treatment)
anova(m2, m2a)
## DOM fluorescence index ----
# build model
m3 <- lme(
DOM.FI ~ treatment,
random = ~ 1 | block,
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r3 <- residuals(m3, type = "normalized")
par(mfrow = c(1, 3))
plot(r3 ~ fitted(m3))
boxplot(r3 ~ soilPools$treatment)
hist(r3)
# test main effects
m3a <- update(m3, ~.- treatment)
anova(m3, m3a)
## DOM component 1 (named component #3 in MS) ----
# build model
m7 <- lme(
DOM.C1 ~ treatment,
random = ~ 1 | block,
weights = varIdent(form = ~ 1 | treatment),
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r7 <- residuals(m7, type = "normalized")
par(mfrow = c(1, 3))
plot(r7 ~ fitted(m7))
boxplot(r7 ~ soilPools$treatment)
hist(r7)
# test main effects
m7a <- update(m7, ~.- treatment)
anova(m7, m7a)
## DOM component 2 (named component #4 in MS) ----
# build model
m8 <- lme(
DOM.C2 ~ treatment,
random = ~ 1 | block,
weights = varIdent(form = ~ 1 | treatment),
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r8 <- residuals(m8, type = "normalized")
par(mfrow = c(1, 3))
plot(r8 ~ fitted(m8))
boxplot(r8 ~ soilPools$treatment)
hist(r8)
# test main effects
m8a <- update(m8, ~.- treatment)
anova(m8, m8a)
## DOM component 3 (named component #6 in MS) ----
# build model
m9 <- lme(
DOM.C3 ~ treatment,
random = ~ 1 | block,
weights = varIdent(form = ~ 1 | treatment),
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r9 <- residuals(m9, type = "normalized")
par(mfrow = c(1, 3))
plot(r9 ~ fitted(m9))
boxplot(r9 ~ soilPools$treatment)
hist(r9)
# test main effects
m9a <- update(m9, ~.- treatment)
anova(m9, m9a)
## DOM component 4 (named component #1 in MS) ----
# build model
m10 <- lme(
DOM.C4 ~ treatment,
random = ~ 1 | block,
weights = varIdent(form = ~ 1 | treatment),
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r10 <- residuals(m10, type = "normalized")
par(mfrow = c(1, 3))
plot(r10 ~ fitted(m10))
boxplot(r10 ~ soilPools$treatment)
hist(r10)
# test main effects
m10a <- update(m10, ~.- treatment)
anova(m10, m10a)
## DOM component 5 (named component #5 in MS) ----
# build model
m11 <- lme(
DOM.C5 ~ treatment,
random = ~ 1 | block,
weights = varIdent(form = ~ 1 | treatment),
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r11 <- residuals(m11, type = "normalized")
par(mfrow = c(1, 3))
plot(r11 ~ fitted(m11))
boxplot(r11 ~ soilPools$treatment)
hist(r11)
# test main effects
m11a <- update(m11, ~.- treatment)
anova(m11, m11a)
## DOM component 6 (named component #2 in MS) ----
# build model
m12 <- lme(
DOM.C6 ~ treatment,
random = ~ 1 | block,
weights = varIdent(form = ~ 1 | treatment),
data = soilPools,
na.action = "na.exclude",
method = "ML"
)
# diagnose model
r12 <- residuals(m12, type = "normalized")
par(mfrow = c(1, 3))
plot(r12 ~ fitted(m12))
boxplot(r12 ~ soilPools$treatment)
hist(r12)
# test main effects
m12a <- update(m12, ~.- treatment)
anova(m12, m12a)
#### PLOT ----------------------------------------------------------------------
a350 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(20, 140)) +
aes(x = treatment,
y = a350_mean,
ymax = a350_mean + a350_se,
ymin = a350_mean - a350_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$a350_mean - allBare$a350_se) +
geom_hline(yintercept = allBare$a350_mean + allBare$a350_se) +
xlab("") +
ylab(expression(paste(a[350], " (", m^{-1}, ")", sep = "")))
ftot <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(5, 32)) +
aes(x = treatment,
y = ftot_mean,
ymax = ftot_mean + ftot_se,
ymin = ftot_mean - ftot_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$ftot_mean - allBare$ftot_se) +
geom_hline(yintercept = allBare$ftot_mean + allBare$ftot_se) +
xlab("") +
ylab(expression(paste(F[tot], " (R. U.)", sep = "")))
fi <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(1.55, 1.85)) +
aes(x = treatment,
y = fi_mean,
ymax = fi_mean + fi_se,
ymin = fi_mean - fi_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$fi_mean - allBare$fi_se) +
geom_hline(yintercept = allBare$fi_mean + allBare$fi_se) +
xlab("") +
ylab("Fluorescence Index")
dom1to3 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(2, 45)) +
aes(x = treatment,
y = dom1_mean,
ymax = dom1_mean + dom1_se,
ymin = dom1_mean - dom1_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$dom1_mean - allBare$dom1_se) +
geom_hline(yintercept = allBare$dom1_mean + allBare$dom1_se) +
xlab("") +
ylab("C3 (humic-like, %)")
dom2to4 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(2, 45)) +
aes(x = treatment,
y = dom2_mean,
ymax = dom2_mean + dom2_se,
ymin = dom2_mean - dom2_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$dom2_mean - allBare$dom2_se) +
geom_hline(yintercept = allBare$dom2_mean + allBare$dom2_se) +
xlab("") +
ylab("C4 (humic-like, %)")
dom3to6 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(2, 45)) +
aes(x = treatment,
y = dom3_mean,
ymax = dom3_mean + dom3_se,
ymin = dom3_mean - dom3_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$dom3_mean - allBare$dom3_se) +
geom_hline(yintercept = allBare$dom3_mean + allBare$dom3_se) +
xlab("") +
ylab("C6 (fumic acid-like, %)")
dom4to1 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(2, 45)) +
aes(x = treatment,
y = dom4_mean,
ymax = dom4_mean + dom4_se,
ymin = dom4_mean - dom4_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$dom4_mean - allBare$dom4_se) +
geom_hline(yintercept = allBare$dom4_mean + allBare$dom4_se) +
xlab("") +
ylab("C1 (protein-like, %)")
dom5to5 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(2, 45)) +
aes(x = treatment,
y = dom5_mean,
ymax = dom5_mean + dom5_se,
ymin = dom5_mean - dom5_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$dom5_mean - allBare$dom5_se) +
geom_hline(yintercept = allBare$dom5_mean + allBare$dom5_se) +
xlab("") +
ylab("C5 (humic-like, %)")
dom6to2 <- ggplot(sumSoil) +
theme_bw() +
theme(panel.grid = element_blank()) +
coord_cartesian(ylim = c(2, 45)) +
aes(x = treatment,
y = dom6_mean,
ymax = dom6_mean + dom6_se,
ymin = dom6_mean - dom6_se) +
geom_errorbar(width = 0.2) +
geom_bar(stat = "identity", col = "black", fill = "white") +
geom_hline(yintercept = allBare$dom6_mean - allBare$dom6_se) +
geom_hline(yintercept = allBare$dom6_mean + allBare$dom6_se) +
xlab("") +
ylab("C2 (protein-like, %)")
#### BUILD PANEL ---------------------------------------------------------------
postscript(
file = "./figure_builds/glasshouse_dom.eps",
width = mm2in(150),
height = mm2in(120)
)
cowplot::plot_grid(
a350, ftot, fi,
dom4to1, dom6to2, dom1to3,
dom2to4, dom5to5, dom3to6, dom3to6,
nrow = 2,
align = "hv"
)
dev.off()
