library(readr)
library(dplyr)
library(tidyr)
library(purrr)
library(tibble)
library(ggplot2)
library(showtext)
library(car)
library(emmeans)
library(multcomp)
library(broom)
library(ggcorrplot)
font_add("Times New Roman",
regular = "/Library/Fonts/Times New Roman.ttf",
bold = "/Library/Fonts/Times New Roman Bold.ttf",
italic = "/Library/Fonts/Times New Roman Italic.ttf")
showtext_auto()
showtext_opts(dpi = 600)
data_path <- "/Users/mrswinstead/Downloads/mnScript_Data.csv"
output_dir <- "/Users/mrswinstead/Desktop/Manuscript Images"
stats_dir <- file.path(output_dir, "Stats")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(stats_dir, recursive = TRUE, showWarnings = FALSE)
dat <- read_csv(data_path, show_col_types = FALSE)
df_base <- dat %>%
filter(!is.na(Graz), !is.na(UpLow), !is.na(BurnInterval)) %>%
mutate(
Graz2 = case_when(
tolower(trimws(as.character(Graz))) %in% c("no","ungrazed","u","0") ~ "Ungrazed",
tolower(trimws(as.character(Graz))) %in% c("yes","grazed","bison","b","1") ~ "Bison grazed",
TRUE ~ NA_character_
),
BurnInterval = factor(BurnInterval, levels = c("01Y","02Y","04Y","20Y")),
UpLow = case_when(
UpLow %in% c("Up","UP","Upland","UPLAND","UPLANDS") ~ "UPLANDS",
UpLow %in% c("Low","LOW","Lowland","LOWLAND","LOWLANDS") ~ "LOWLANDS",
TRUE ~ as.character(UpLow)
),
UpLow = factor(UpLow, levels = c("UPLANDS","LOWLANDS")),
Graz2 = factor(Graz2, levels = c("Ungrazed","Bison grazed"))
)
plot_specs <- tibble(
source_col = c("GWC","pH","TotalC","TotalN","total_CN",
"ExtNH4N","ExtNO3N","ExtTIN","NH4NO3",
"PMN_7d","PMN_28d","ACEprotein",
"BG_drysoil","NAG_drysoil","LAP_drysoil","BGNAGLAP"),
file_name = c("GWC","pH","TotalC","TotalN","total_CN",
"ExtNH4N","ExtNO3N","ExtTIN","NH4NO3",
"PMN_7d","PMN_28d","ACEprotein",
"BG_drysoil","NAG_drysoil","LAP_drysoil","BGNAGLAP"),
log_transform = c(TRUE, FALSE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE)
)
y_labels <- list(
GWC = expression(bold(paste("Soil water content (g g", ""^-1, ")"))),
pH = expression(bold("Soil pH")),
TotalC = expression(bold("Total soil C (%)")),
TotalN = expression(bold("Total soil N (%)")),
total_CN = expression(bold("Soil molar C:N")),
ExtNH4N = expression(bold(paste("Extractable ", NH[4]^"+", "-N (", mu, "g g", ""^-1, ")"))),
ExtNO3N = expression(bold(paste("Extractable ", NO[3]^"-", "-N (", mu, "g g", ""^-1, ")"))),
ExtTIN = expression(bold(paste("Total inorganic N (", mu, "g g", ""^-1, ")"))),
NH4NO3 = expression(bold(paste(NH[4]^"+", "-N : ", NO[3]^"-", "-N"))),
PMN_7d = expression(bold(paste("PMN-7d (", mu, "g N kg", ""^-1, ")"))),
PMN_28d = expression(bold(paste("PMN-28d (", mu, "g N kg", ""^-1, ")"))),
ACEprotein = expression(bold(paste("ACE protein (mg g", ""^-1, ")"))),
BG_drysoil = expression(bold(paste(beta, "G activity (nmol g", ""^-1, " h", ""^-1, ")"))),
NAG_drysoil = expression(bold(paste("NAG activity (nmol g", ""^-1, " h", ""^-1, ")"))),
LAP_drysoil = expression(bold(paste("LAP activity (nmol g", ""^-1, " h", ""^-1, ")"))),
BGNAGLAP = expression(bold(paste("ln(", beta, "G) / ln(NAG+LAP)")))
)
fig_theme <- theme_bw(base_size = 10, base_family = "Times New Roman") +
theme(
text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
axis.title.x = element_text(family = "Times New Roman", size = 16, color = "black", face = "bold"),
axis.title.y = element_text(family = "Times New Roman", size = 16, color = "black", face = "bold"),
axis.text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
axis.text.x = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
axis.text.y = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
strip.text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
legend.text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
legend.title = element_blank(),
strip.background = element_blank(),
legend.position = "bottom",
panel.grid.minor = element_blank()
)
fill_palette <- c("Ungrazed" = "#E69F00", "Bison grazed" = "#56B4E9")
for (i in seq_len(nrow(plot_specs))) {
source_col <- plot_specs$source_col[i]
file_name <- plot_specs$file_name[i]
y_label <- y_labels[[file_name]]
if (!source_col %in% names(df_base)) next
df_plot <- df_base %>%
mutate(raw_value = as.numeric(.data[[source_col]])) %>%
filter(!is.na(raw_value), !is.na(BurnInterval), !is.na(Graz2), !is.na(UpLow))
p <- ggplot(df_plot, aes(x = BurnInterval, y = raw_value, fill = Graz2)) +
geom_boxplot(width = 0.75,
position = position_dodge(width = 0.8),
outlier.size = 0.8) +
facet_wrap(~ UpLow, ncol = 1) +
scale_fill_manual(values = fill_palette) +
labs(x = "Burn interval", y = y_label, fill = NULL) +
fig_theme
ggsave(file.path(output_dir, paste0(file_name, "_boxplot.pdf")),
p, width = 5, height = 7, device = cairo_pdf)
ggsave(file.path(output_dir, paste0(file_name, "_boxplot.png")),
p, width = 5, height = 7, dpi = 600)
ggsave(file.path(output_dir, paste0(file_name, "_boxplot.tiff")),
p, width = 5, height = 7, dpi = 600, compression = "lzw")
}
master_anova <- list()
master_letters <- list()
master_summary <- list()
master_diag <- list()
for (i in seq_len(nrow(plot_specs))) {
source_col <- plot_specs$source_col[i]
file_name <- plot_specs$file_name[i]
log_it <- plot_specs$log_transform[i]
if (!source_col %in% names(df_base)) next
df_var <- df_base %>%
mutate(raw_value = as.numeric(.data[[source_col]])) %>%
filter(!is.na(raw_value), !is.na(BurnInterval), !is.na(Graz2), !is.na(UpLow))
if (log_it) {
df_var <- df_var %>%
filter(raw_value > 0) %>%
mutate(stats_value = log(raw_value))
} else {
df_var <- df_var %>%
mutate(stats_value = raw_value)
}
options(contrasts = c("contr.sum", "contr.poly"))
mod_t3 <- lm(stats_value ~ BurnInterval * Graz2 * UpLow, data = df_var)
anova_table <- car::Anova(mod_t3, type = 3) %>%
broom::tidy() %>%
mutate(Variable = file_name,
Transformation = ifelse(log_it, "log", "raw"),
across(where(is.numeric), ~ signif(.x, 4))) %>%
relocate(Variable, Transformation)
write_csv(anova_table, file.path(stats_dir, paste0(file_name, "_ANOVA.csv")))
master_anova[[file_name]] <- anova_table
options(contrasts = c("contr.treatment", "contr.poly"))
mod_em <- lm(stats_value ~ BurnInterval * Graz2 * UpLow, data = df_var)
shapiro_p <- shapiro.test(residuals(mod_em))$p.value
levene_p <- car::leveneTest(stats_value ~ BurnInterval * Graz2 * UpLow,
data = df_var)$`Pr(>F)`[1]
diag_tbl <- tibble(
Variable = file_name,
Transformation = ifelse(log_it, "log", "raw"),
Shapiro_p = signif(shapiro_p, 4),
Levene_p = signif(levene_p, 4)
)
master_diag[[file_name]] <- diag_tbl
emm_letters <- emmeans::emmeans(mod_em, ~ BurnInterval * Graz2 | UpLow)
cld_tbl <- multcomp::cld(emm_letters, Letters = letters, adjust = "tukey") %>%
as.data.frame() %>%
rename(LetterGroup = .group) %>%
mutate(Variable = file_name,
LetterGroup = trimws(LetterGroup),
across(where(is.numeric), ~ signif(.x, 4))) %>%
relocate(Variable) %>%
arrange(UpLow, Graz2, BurnInterval)
write_csv(cld_tbl, file.path(stats_dir, paste0(file_name, "_Tukey_letters.csv")))
master_letters[[file_name]] <- cld_tbl
pairs_graz <- contrast(
emmeans::emmeans(mod_em, ~ Graz2 | BurnInterval * UpLow),
method = "revpairwise", adjust = "tukey"
) %>%
broom::tidy() %>%
mutate(Variable = file_name,
across(where(is.numeric), ~ signif(.x, 4))) %>%
relocate(Variable)
write_csv(pairs_graz, file.path(stats_dir, paste0(file_name, "_grazing_contrasts.csv")))
pairs_burn <- contrast(
emmeans::emmeans(mod_em, ~ BurnInterval | Graz2 * UpLow),
method = "pairwise", adjust = "tukey"
) %>%
broom::tidy() %>%
mutate(Variable = file_name,
across(where(is.numeric), ~ signif(.x, 4))) %>%
relocate(Variable)
write_csv(pairs_burn, file.path(stats_dir, paste0(file_name, "_burn_contrasts.csv")))
summary_tbl <- df_var %>%
group_by(UpLow, Graz2, BurnInterval) %>%
summarise(
n = n(),
mean = mean(raw_value, na.rm = TRUE),
sd = sd(raw_value, na.rm = TRUE),
median = median(raw_value, na.rm = TRUE),
min = min(raw_value, na.rm = TRUE),
max = max(raw_value, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(Variable = file_name,
across(c(mean, sd, median, min, max), ~ signif(.x, 4))) %>%
relocate(Variable)
write_csv(summary_tbl, file.path(stats_dir, paste0(file_name, "_summary_stats.csv")))
master_summary[[file_name]] <- summary_tbl
}
options(contrasts = c("contr.treatment", "contr.poly"))
write_csv(bind_rows(master_anova), file.path(stats_dir, "MASTER_ANOVA.csv"))
write_csv(bind_rows(master_letters), file.path(stats_dir, "MASTER_Tukey_letters.csv"))
write_csv(bind_rows(master_summary), file.path(stats_dir, "MASTER_summary_stats.csv"))
write_csv(bind_rows(master_diag), file.path(stats_dir, "MASTER_diagnostics.csv"))
heatmap_vars <- c(
"NH4NO3", "NAG_drysoil", "ExtTIN", "GWC", "PMN_28d",
"ACEprotein", "LAP_drysoil", "BGNAGLAP",
"TotalN", "PMN_7d", "total_CN", "pH"
)
df_corr <- df_base %>%
dplyr::select(all_of(heatmap_vars)) %>%
mutate(across(everything(), ~ as.numeric(.x)))
df_pearson <- df_corr %>%
mutate(across(-pH, ~ ifelse(.x > 0, log(.x), NA_real_)))
pearson_corr <- cor(df_pearson, use = "pairwise.complete.obs", method = "pearson")
pearson_p <- cor_pmat(df_pearson, method = "pearson")
spearman_corr <- cor(df_corr, use = "pairwise.complete.obs", method = "spearman")
spearman_p <- cor_pmat(df_corr, method = "spearman")
fdr_pmat <- function(pmat) {
upper_idx <- upper.tri(pmat)
p_upper <- pmat[upper_idx]
p_adj <- p.adjust(p_upper, method = "fdr")
out <- pmat
out[upper_idx] <- p_adj
out[lower.tri(out)] <- t(out)[lower.tri(out)]
diag(out) <- 0
out
}
pearson_p_fdr <- fdr_pmat(pearson_p)
spearman_p_fdr <- fdr_pmat(spearman_p)
heatmap_theme <- theme(
text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
axis.text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
axis.text.x = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold",
angle = 45, hjust = 1),
axis.text.y = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
legend.text = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
legend.title = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold"),
plot.title = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold",
hjust = 0.5),
plot.subtitle = element_text(family = "Times New Roman", size = 10, color = "black", face = "bold",
hjust = 0.5),
panel.grid = element_blank()
)
heatmap_colors <- c("#56B4E9", "white", "#E69F00")
p_pearson <- ggcorrplot(
pearson_corr,
p.mat = pearson_p_fdr,
type = "lower",
lab = TRUE,
lab_size = 3,
insig = "blank",
sig.level = 0.05,
hc.order = FALSE,
colors = heatmap_colors,
outline.color = "white",
legend.title = "Corr"
) +
labs(
title = "Pearson correlations on ln-transformed N-focused variables",
subtitle = "Blank cells: not significant after FDR correction (p > 0.05)"
) +
heatmap_theme
p_spearman <- ggcorrplot(
spearman_corr,
p.mat = spearman_p_fdr,
type = "lower",
lab = TRUE,
lab_size = 3,
insig = "blank",
sig.level = 0.05,
hc.order = FALSE,
colors = heatmap_colors,
outline.color = "white",
legend.title = "Corr"
) +
labs(
title = "Spearman correlations among N-focused soil variables",
subtitle = "Blank cells: not significant after FDR correction (p > 0.05)"
) +
heatmap_theme
ggsave(file.path(output_dir, "Heatmap_Pearson.pdf"), p_pearson,
width = 8, height = 7, device = cairo_pdf)
ggsave(file.path(output_dir, "Heatmap_Pearson.png"), p_pearson,
width = 8, height = 7, dpi = 600)
ggsave(file.path(output_dir, "Heatmap_Pearson.tiff"), p_pearson,
width = 8, height = 7, dpi = 600, compression = "lzw")
ggsave(file.path(output_dir, "Heatmap_Spearman.pdf"), p_spearman,
width = 8, height = 7, device = cairo_pdf)
ggsave(file.path(output_dir, "Heatmap_Spearman.png"), p_spearman,
width = 8, height = 7, dpi = 600)
ggsave(file.path(output_dir, "Heatmap_Spearman.tiff"), p_spearman,
width = 8, height = 7, dpi = 600, compression = "lzw")
corr_long <- function(corr_mat, p_mat, method_label) {
vars <- rownames(corr_mat)
out <- list()
for (i in 2:length(vars)) {
for (j in 1:(i - 1)) {
out[[length(out) + 1]] <- tibble(
Method = method_label,
Variable_1 = vars[i],
Variable_2 = vars[j],
Coefficient = signif(corr_mat[i, j], 4),
P_FDR = signif(p_mat[i, j], 4)
)
}
}
bind_rows(out)
}
correlations_table <- bind_rows(
corr_long(pearson_corr, pearson_p_fdr, "Pearson (ln-transformed)"),
corr_long(spearman_corr, spearman_p_fdr, "Spearman (raw)")
) %>%
mutate(Significant = ifelse(P_FDR < 0.05, "yes", "no"))
write_csv(correlations_table, file.path(stats_dir, "MASTER_correlations.csv"))