Christmas_CTD_Final
Daniel Kaufman
2025-12-20
Downloading libraries for project:
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.5.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lubridate)
library(patchwork)## Warning: package 'patchwork' was built under R version 4.5.2library(ggplot2)
library(readr)
library(cowplot)## Warning: package 'cowplot' was built under R version 4.5.2##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:patchwork':
##
## align_plots
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(rlang)##
## Attaching package: 'rlang'
##
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splicelibrary(knitr)
library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.5.2##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(nlme)##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(purrr)
library(emmeans)## Warning: package 'emmeans' was built under R version 4.5.2## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'Introduction - Need to do more Literature Review for the Narrative: Revisit when other sections are done: Methods - Re- do diagrams, explain flux theory, and theory needed for stats, locations, data cleaning, and more revisit when Results is done.
Results Section
Functions for creating figures:
# importing data sets
Flux_Data_With_Covars <- read_csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/Flux_Data_With_Covars.csv")## Rows: 21090 Columns: 82
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): Pile, Chamber_Corrected, Turning_Status, Date, DATE_TIME
## dbl (77): Chamber_Elevation, DOY, air_temperature, RH, VPD, Tdew, wind_speed...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Functions to clean data
# function used to remove NA values which LICOR uses -9999 for:
clean_numeric <- function(x,
min_val = -Inf,
max_val = Inf,
bad_values = c(-9999, -9999, 9999, 99999)) {
x <- as.numeric(x)
x[x %in% bad_values] <- NA_real_
x[x < min_val | x > max_val] <- NA_real_
x
}
# function that calculates the SE of a function:
se <- function(x) sd(x, na.rm = TRUE) / sqrt(sum(!is.na(x)))
Cleaned_Data <- Flux_Data_With_Covars %>%
mutate(
# DOY has decimals, so collapse to whole day
DOY = floor(DOY),
pile_temp = clean_numeric(TS_1.initial_value, min_val = -10, max_val = 90),
pile_moist = clean_numeric(SWC_1.initial_value, min_val = 0.01, max_val = 0.98),
BulkDensity = clean_numeric(BulkDensity, min_val = 0.01, max_val = 1)
)
daily_pile_doy <- Cleaned_Data %>%
group_by(Pile, DOY) %>%
summarise(
n_obs = sum(!is.na(pile_temp) | !is.na(pile_moist)),
pile_temp_mean = mean(pile_temp, na.rm = TRUE),
pile_temp_se = se(pile_temp),
pile_moist_mean = mean(pile_moist, na.rm = TRUE),
pile_moist_se = se(pile_moist),
.groups = "drop"
)
# making weekly figures to reduce noise for final paper figures:
summarise_weekly_mean_se <- function(df,
doy,
value,
group,
min_n = 2) {
doy <- enquo(doy)
value <- enquo(value)
group <- enquo(group)
df %>%
filter(!is.na(!!value), !is.na(!!doy), !is.na(!!group)) %>%
mutate(
week = floor((!!doy - 1) / 7) + 1
) %>%
group_by(!!group, week) %>%
summarise(
DOY = mean(!!doy, na.rm = TRUE), # x-position for plotting
mean = mean(!!value, na.rm = TRUE),
sd = sd(!!value, na.rm = TRUE),
n = sum(!is.na(!!value)),
se = sd / sqrt(n),
.groups = "drop"
) %>%
filter(n >= min_n)
}
weekly_moist <- summarise_weekly_mean_se(
df = Cleaned_Data,
doy = DOY,
value = pile_moist,
group = Pile,
min_n = 2
)
weekly_temp <- summarise_weekly_mean_se(
df = Cleaned_Data,
doy = DOY,
value = pile_temp,
group = Pile,
min_n = 2
)Plotting Moisture and Temperature for the experimental cycle.
plot_points_err <- function(df,
x, y_mean, y_se,
group,
cols = NULL,
shapes = NULL,
point_size = 2.0,
point_alpha = 0.85,
err_alpha = 0.6,
err_lwd = 0.4,
err_width = 0,
base_size = 12,
legend_position = "top",
drop_na = TRUE,
x_lab = NULL,
y_lab = NULL,
legend_title = NULL,
vline = NULL,
vlines_red = NULL,
vlines_black = NULL,
vline_lwd = 0.4,
vline_alpha = 0.8,
vline_linetype = "dashed",
vline_red_linetype = "solid",
vline_black_linetype = "dashed",
vline_red_colour = "#9b1c31",
vline_black_colour = "#000000"){
# capture column names
x <- enquo(x)
y_mean <- enquo(y_mean)
y_se <- enquo(y_se)
group <- enquo(group)
# optionally drop NA rows for required columns
if (drop_na) {
df <- df[complete.cases(
df[, c(as_name(x), as_name(y_mean), as_name(y_se), as_name(group))]
), ]
}
p <- ggplot(df, aes(x = !!x, y = !!y_mean)) +
# error bars first (so points draw on top)
geom_errorbar(
aes(ymin = (!!y_mean) - (!!y_se),
ymax = (!!y_mean) + (!!y_se),
colour = !!group),
width = err_width,
alpha = err_alpha,
linewidth = err_lwd
)
if (!is.null(vlines_red)) {
p <- p + geom_vline(
xintercept = vlines_red,
colour = vline_red_colour,
linetype = vline_red_linetype,
linewidth = vline_lwd,
alpha = vline_alpha
)
}
if (!is.null(vlines_black)) {
p <- p + geom_vline(
xintercept = vlines_black,
colour = vline_black_colour,
linetype = vline_black_linetype,
linewidth = vline_lwd,
alpha = vline_alpha
)
}
p <- p +
# points (no lines)
geom_point(
aes(colour = !!group, shape = !!group),
size = point_size,
alpha = point_alpha
) +
theme_bw(base_size = base_size) +
theme(
legend.position = legend_position,
panel.grid.minor = element_blank()
)
# manual scales if provided (keeps consistent styling across plots)
if (!is.null(cols)) {
p <- p + scale_colour_manual(values = cols)
}
if (!is.null(shapes)) {
p <- p + scale_shape_manual(values = shapes)
}
# ── LABELS SECTION ─────────────────────────────
if (!is.null(x_lab) || !is.null(y_lab) || !is.null(legend_title)) {
p <- p + labs(
x = x_lab,
y = y_lab,
colour = legend_title,
shape = legend_title
)
}
return(p)
}Pile Temperature Plot using Daily Average Values Pile Moisture Using Daily Average Values
#turning data to add to plots if needed
turns_both_red <- c(159, 186, 208, 234, 262, 290)
turns_control_black <- c(173, 191, 222, 249, 276)
turns_both_red_week <- c(159, 186, 208, 234, 262, 290)/7
turns_control_black_week <- c(173, 191, 222, 249, 276)/7
pile_cols <- c(C = "#000000", E = "#9b1c31")
pile_shapes <- c(C = 16, E = 17)
fig_moist_doy <- plot_points_err(
daily_pile_doy,
x = DOY,
y_mean = pile_moist_mean,
y_se = pile_moist_se,
group = Pile,
cols = pile_cols,
shapes = pile_shapes,
point_size = 1.8,
x_lab = "Day of year",
y_lab = expression("Volumetric water content (m"^3*" m"^-3*")"),
legend_title = "Pile",
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_moist_doy
fig_moist_weekly <- plot_points_err(
weekly_moist,
x = week,
y_mean = mean,
y_se = se,
group = Pile,
cols = pile_cols,
shapes = pile_shapes,
point_size = 1.8,
x_lab = "Week of year",
y_lab = expression("Volumetric water content (m"^3*" m"^-3*")"),
legend_title = "Pile",
vlines_red = turns_both_red_week,
vlines_black = turns_control_black_week,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_moist_weekly
fig_temp_doy <- plot_points_err(
daily_pile_doy,
x = DOY,
y_mean =pile_temp_mean,
y_se = pile_temp_se,
group = Pile,
cols = pile_cols,
shapes = pile_shapes,
point_size = 1.8,
x_lab = "Day of year",
y_lab = expression("Pile Temperature Degrees Celsius"),
legend_title = "Pile",
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_temp_doy
fig_temp_weekly <- plot_points_err(
weekly_temp,
x = week,
y_mean =mean,
y_se = se,
group = Pile,
cols = pile_cols,
shapes = pile_shapes,
point_size = 1.8,
x_lab = "Week of year",
y_lab = expression("Pile Temperature Degrees Celsius"),
legend_title = "Pile",
vlines_red = turns_both_red_week,
vlines_black = turns_control_black_week,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_temp_weekly
Creating Tables with Pile Variables:
pile_vars_temp <- Cleaned_Data %>%
filter(!is.na(pile_temp)) %>%
group_by(Pile) %>%
summarise(
n_obs = n(),
mean_temp = mean(pile_temp),
sd_temp = sd(pile_temp),
se_temp = sd_temp/sqrt(n_obs)
)
pile_vars_water <- Cleaned_Data %>%
filter(!is.na(pile_moist)) %>%
group_by(Pile) %>%
summarise(
n_obs = n(),
mean_moist = mean(pile_moist),
sd_moist = sd(pile_moist),
se_moist = sd_moist/sqrt(n_obs)
)Displaying Tables of Summary Stats for Moisture and Temperature
| Pile | n_obs | mean_moist | sd_moist | se_moist |
|---|---|---|---|---|
| C | 9177 | 0.490 | 0.327 | 0.003 |
| E | 9403 | 0.787 | 0.244 | 0.003 |
| Pile | n_obs | mean_temp | sd_temp | se_temp |
|---|---|---|---|---|
| C | 9192 | 54.069 | 8.992 | 0.094 |
| E | 9420 | 51.778 | 12.812 | 0.132 |
Examining Measured Fluxes by Pile
library(readxl)
co2_raw <- read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/CO2_R_2_filtered.csv")
ch4_raw <- read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/CH4_R_2_filtered.csv")
n2o_raw <- read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_clean/N2O_R_2_filtered.csv")Making Chart with Average values by Pile and SD,SE
pile_vars_co2 <- co2_raw %>%
filter(!is.na(FCO2_DRY.LIN)) %>%
group_by(Pile) %>%
summarise(
n_obs = n(),
mean_co2 = mean(FCO2_DRY.LIN),
sd_co2 = sd(FCO2_DRY.LIN),
se_co2 = sd_co2/sqrt(n_obs)
)
pile_vars_ch4 <- ch4_raw %>%
filter(!is.na(FCH4_DRY.LIN)) %>%
group_by(Pile) %>%
summarise(
n_obs = n(),
mean_ch4 = mean(FCH4_DRY.LIN),
sd_ch4 = sd(FCH4_DRY.LIN),
se_ch4 = sd_ch4/sqrt(n_obs)
)
pile_vars_n2o <- n2o_raw %>%
filter(!is.na(FN2O_DRY.LIN)) %>%
group_by(Pile) %>%
summarise(
n_obs = n(),
mean_n2o = mean(FN2O_DRY.LIN),
sd_n2o = sd(FN2O_DRY.LIN),
se_n2o = sd_n2o/sqrt(n_obs)
)Displaying Average Values for Each Trace Gas:
kable(
pile_vars_co2,
digits = 3,
caption = "Summary statistics for co2 flux by pile."
) %>%
kable_styling(
full_width = FALSE,
position = "center",
bootstrap_options = c("bordered", "striped")
)| Pile | n_obs | mean_co2 | sd_co2 | se_co2 |
|---|---|---|---|---|
| C | 5064 | 190.076 | 202.288 | 2.843 |
| E | 4215 | 191.757 | 193.987 | 2.988 |
kable(
pile_vars_ch4,
digits = 3,
caption = "Summary statistics for ch4 flux by pile."
) %>%
kable_styling(
full_width = FALSE,
position = "center",
bootstrap_options = c("bordered", "striped")
)| Pile | n_obs | mean_ch4 | sd_ch4 | se_ch4 |
|---|---|---|---|---|
| C | 5355 | 7055.016 | 11741.80 | 160.456 |
| E | 4937 | 6594.771 | 10324.23 | 146.935 |
kable(
pile_vars_n2o,
digits = 3,
caption = "Summary statistics for n2o flux by pile."
) %>%
kable_styling(
full_width = FALSE,
position = "center",
bootstrap_options = c("bordered", "striped")
)| Pile | n_obs | mean_n2o | sd_n2o | se_n2o |
|---|---|---|---|---|
| C | 6906 | 157.089 | 295.036 | 3.550 |
| E | 6911 | 98.481 | 245.689 | 2.955 |
Creating Daily and Weekly Plots for all Fluxes
# N2O
daily_n2o <- n2o_raw %>%
mutate(DOY = floor(DOY)) %>% # collapse fractional DOY
group_by(Pile, DOY) %>%
summarise(
n_obs = sum(!is.na(FN2O_DRY.LIN)),
n2o_mean = mean(FN2O_DRY.LIN, na.rm = TRUE),
n2o_se = se(FN2O_DRY.LIN),
.groups = "drop"
) %>%
filter(n_obs >= 3) # tweak threshold if needed
# CO2
daily_co2 <- co2_raw %>%
mutate(DOY = floor(DOY)) %>%
group_by(Pile, DOY) %>%
summarise(
n_obs = sum(!is.na(FCO2_DRY.LIN)),
co2_mean = mean(FCO2_DRY.LIN, na.rm = TRUE),
co2_se = se(FCO2_DRY.LIN),
.groups = "drop"
) %>%
filter(n_obs >= 3)
# CH4
daily_ch4 <- ch4_raw %>%
mutate(DOY = floor(DOY)) %>%
group_by(Pile, DOY) %>%
summarise(
n_obs = sum(!is.na(FCH4_DRY.LIN)),
ch4_mean = mean(FCH4_DRY.LIN, na.rm = TRUE),
ch4_se = se(FCH4_DRY.LIN),
.groups = "drop"
) %>%
filter(n_obs >= 3)
weekly_n2o <- n2o_raw %>%
mutate(
week = floor((DOY - 1) / 7) + 1
) %>%
group_by(Pile, week) %>%
summarise(
n_obs = sum(!is.na(FN2O_DRY.LIN)),
n2o_mean = mean(FN2O_DRY.LIN, na.rm = TRUE),
n2o_se = sd(FN2O_DRY.LIN, na.rm = TRUE) / sqrt(n_obs),
DOY = mean(DOY, na.rm = TRUE), # x-position for plotting (optional but useful)
.groups = "drop"
) %>%
filter(n_obs >= 3)
# CO2
weekly_co2 <- co2_raw %>%
mutate(
week = floor((DOY - 1) / 7) + 1
) %>%
group_by(Pile, week) %>%
summarise(
n_obs = sum(!is.na(FCO2_DRY.LIN)),
co2_mean = mean(FCO2_DRY.LIN, na.rm = TRUE),
co2_se = sd(FCO2_DRY.LIN, na.rm = TRUE) / sqrt(n_obs),
DOY = mean(DOY, na.rm = TRUE),
.groups = "drop"
) %>%
filter(n_obs >= 3)
# CH4
weekly_ch4 <- ch4_raw %>%
mutate(
week = floor((DOY - 1) / 7) + 1
) %>%
group_by(Pile, week) %>%
summarise(
n_obs = sum(!is.na(FCH4_DRY.LIN)),
ch4_mean = mean(FCH4_DRY.LIN, na.rm = TRUE),
ch4_se = sd(FCH4_DRY.LIN, na.rm = TRUE) / sqrt(n_obs),
DOY = mean(DOY, na.rm = TRUE),
.groups = "drop"
) %>%
filter(n_obs >= 3)Plotting Values
fig_co2_daily <- plot_points_err(
daily_co2,
x = DOY, y_mean = co2_mean, y_se = co2_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_co2_daily
fig_co2_weekly <- plot_points_err(
weekly_co2,
x = week, y_mean = co2_mean, y_se = co2_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red_week,
vlines_black = turns_control_black_week,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_co2_weekly
fig_ch4_daily <- plot_points_err(
daily_ch4,
x = DOY, y_mean = ch4_mean, y_se = ch4_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_ch4_daily
fig_ch4_weekly <- plot_points_err(
weekly_ch4,
x = week, y_mean = ch4_mean, y_se = ch4_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red_week,
vlines_black = turns_control_black_week,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_ch4_weekly
fig_n2o_daily <- plot_points_err(
daily_n2o,
x = DOY, y_mean = n2o_mean, y_se = n2o_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_n2o_daily
fig_n2o_weekly <- plot_points_err(
weekly_n2o,
x = week, y_mean = n2o_mean, y_se = n2o_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red_week,
vlines_black = turns_control_black_week,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_n2o_weekly
Looking at A Combination of all three gases using established 100 year conversion, these can easily be changed:
Co2 1 Ch4 28 N2o 265
To combine Trace Gases as cleaning was done on a per measurement basis:
Method 1: Looking at trace gases where all there trace gases where taken at once Method 2: Averaging Per Day and Comparing Daily Average GHG Potential
Method to keep only rows with all three trace gases
Flux_All_Gases_QC <- Flux_Data_With_Covars %>%
filter(
# CO2
FCO2_DRY.LIN > 0,
FCO2_DRY.LIN_R2 > 0.9,
# CH4
FCH4_DRY.LIN > 0,
FCH4_DRY.LIN_R2 > 0.9,
# N2O
FN2O_DRY.LIN > 0,
FN2O_DRY.LIN_R2 > 0.9
)
Flux_Data_With_Covars %>%
summarise(
total_obs = n(),
pass_all = sum(
FCO2_DRY.LIN > 0 & FCO2_DRY.LIN_R2 > 0.9 &
FCH4_DRY.LIN > 0 & FCH4_DRY.LIN_R2 > 0.9 &
FN2O_DRY.LIN > 0 & FN2O_DRY.LIN_R2 > 0.9,
na.rm = TRUE
),
percent_retained = 100 * pass_all / total_obs
)## # A tibble: 1 × 3
## total_obs pass_all percent_retained
## <int> <int> <dbl>
## 1 21090 5558 26.4Flux_All_Gases_QC <- Flux_All_Gases_QC %>%
mutate(
# keep units explicit
CO2_umol = FCO2_DRY.LIN, # µmol m-2 s-1
N2O_umol = FN2O_DRY.LIN, # µmol m-2 s-1
CH4_umol = FCH4_DRY.LIN * 0.001, # nmol -> µmol
# GWP100 CO2-equivalent (µmol CO2-eq m-2 s-1)
GHG_CO2e_umol = CO2_umol + 28 * CH4_umol + 265 * N2O_umol
)Summary Stats of Co2e by pile when looking at measurements where all three gases passed qa/qc
pile_vars_co2e <- Flux_All_Gases_QC %>%
filter(!is.na(GHG_CO2e_umol)) %>%
group_by(Pile) %>%
summarise(
n_obs = n(),
mean_CO2e = mean(GHG_CO2e_umol),
sd_CO2e = sd(GHG_CO2e_umol),
se_CO2e = sd_CO2e/sqrt(n_obs)
)
kable(
pile_vars_co2e,
digits = 3,
caption = "Summary statistics for co2e flux by pile."
) %>%
kable_styling(
full_width = FALSE,
position = "center",
bootstrap_options = c("bordered", "striped")
)| Pile | n_obs | mean_CO2e | sd_CO2e | se_CO2e |
|---|---|---|---|---|
| C | 3269 | 59570.66 | 101672.24 | 1778.259 |
| E | 2289 | 40046.13 | 96185.23 | 2010.414 |
We are going to write the code to plot this but we will plot all three together - as this is only the first method of plotting the values.
- Lets first create plots for this method then we will work on the other methods.
# daily co2e
daily_co2e <- Flux_All_Gases_QC %>%
mutate(DOY = floor(DOY)) %>%
group_by(Pile, DOY) %>%
summarise(
n_obs = sum(!is.na(GHG_CO2e_umol)),
co2e_mean = mean(GHG_CO2e_umol, na.rm = TRUE),
co2e_se = sd(GHG_CO2e_umol, na.rm = TRUE)/sqrt(n_obs),
.groups = "drop"
) %>%
filter(n_obs >= 3)
# weekly co2e
weekly_co2e <- Flux_All_Gases_QC %>%
mutate(
week = floor((DOY - 1) / 7) + 1
) %>%
group_by(Pile, week) %>%
summarise(
n_obs = sum(!is.na(GHG_CO2e_umol)),
co2e_mean = mean(GHG_CO2e_umol, na.rm = TRUE),
co2e_se = sd(GHG_CO2e_umol, na.rm = TRUE) / sqrt(n_obs),
DOY = mean(DOY, na.rm = TRUE),
.groups = "drop"
) %>%
filter(n_obs >= 3)
# CO2
fig__co2e_daily <- plot_points_err(
daily_co2e,
x = DOY, y_mean = co2e_mean, y_se = co2e_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig__co2e_daily
fig_co2e_weekly <- plot_points_err(
weekly_co2e,
x = week, y_mean =co2e_mean, y_se = co2e_se,
group = Pile,
cols = pile_cols, shapes = pile_shapes,
point_size = 1.4,
vlines_red = turns_both_red_week,
vlines_black = turns_control_black_week,
vline_red_linetype = "solid",
vline_black_linetype = "dashed"
)
fig_co2e_weekly
Method 2: Averaging Per Day and Comparing Daily Average GHG Potential
daily_all3 <- daily_co2 %>%
select(Pile, DOY, n_obs_co2 = n_obs, co2_mean, co2_se) %>%
inner_join(
daily_ch4 %>% select(Pile, DOY, n_obs_ch4 = n_obs, ch4_mean, ch4_se),
by = c("Pile", "DOY")
) %>%
inner_join(
daily_n2o %>% select(Pile, DOY, n_obs_n2o = n_obs, n2o_mean, n2o_se),
by = c("Pile", "DOY")
)
daily_all3 <- daily_all3 %>%
mutate(
ch4_mean_umol = ch4_mean * 0.001,
ch4_se_umol = ch4_se * 0.001,
# CO2-equivalent in µmol CO2-eq m-2 s-1 (using your factors)
co2e_mean = co2_mean + 28 * ch4_mean_umol + 273 * n2o_mean
)
# not going to error bars for this as it is messy but for now we assume they are independent so each one is one measurement and not correlated which is not true:
daily_all3 <- daily_all3 %>%
mutate(
co2e_se = sqrt(
co2_se^2 +
(28 * ch4_se_umol)^2 +
(273 * n2o_se)^2
)
)
fig_co2e_daily_all3 <- plot_points_err(
daily_all3,
x = DOY,
y_mean = co2e_mean,
y_se = co2e_se,
group = Pile,
cols = pile_cols,
shapes = pile_shapes,
point_size = 1.6,
vlines_red = turns_both_red,
vlines_black = turns_control_black,
vline_red_linetype = "solid",
vline_black_linetype = "dashed",
x_lab = "Day of year",
y_lab = expression("CO"[2]*"-eq flux ("*mu*"mol m"^-2*" s"^-1*")"),
legend_title = "Pile"
)
fig_co2e_daily_all3
This is good, I don’t think I need to do a weekly version. Depending on
what we want to do we can revisit.
Now we have looked at the differences in each trace gas and all of them together as well as the baseline differences between the piles (temp, moisture). We did not plot other variables like bulk density but the limited number of samples makes this a little difficult. - Not sure how we should include this, in other parts I have this as just a mean , sd , se.
Variable Unit Value SD SE n Bulk Density C g/cm^3 0.17 3.58 1.03 12 Bulk Density E g/cm^3 0.20 5.17 1.49 12
Basic Model Building:
Numerous models have been employed throughout the process of writting this paper: Here I will present all the models, explain them, and then we can choose the model which best explains the data.
Model 1: Linear Mixed Effects Model:
Here we are looking at the magnitude of fluxes between piles taking into account random effects from chambers and an auto correlation structure.
# log transforming flux data for each flux, grouping chambers together and arranging by doy
dat_co2 <- co2_raw %>%
mutate(log_flux = log(FCO2_DRY.LIN))%>%
arrange(Chamber_Corrected, DOY) %>%
group_by(Chamber_Corrected) %>%
mutate(t_index = row_number()) %>%
ungroup()
# applying model
m_co2 <- lme(
log_flux ~ Pile,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2,
method = "REML"
)
# log transforming flux data for each flux, grouping chambers together and arranging by doy
dat_ch4 <- ch4_raw %>%
mutate(log_flux = log(FCH4_DRY.LIN))%>%
arrange(Chamber_Corrected, DOY) %>%
group_by(Chamber_Corrected) %>%
mutate(t_index = row_number()) %>%
ungroup()
# applying model
m_ch4 <- lme(
log_flux ~ Pile,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4,
method = "REML"
)
# log transforming flux data for each flux, grouping chambers together and arranging by doy
dat_n2o <- n2o_raw %>%
mutate(log_flux = log(FN2O_DRY.LIN))%>%
arrange(Chamber_Corrected, DOY) %>%
group_by(Chamber_Corrected) %>%
mutate(t_index = row_number()) %>%
ungroup()
# applying model
m_n2o <- lme(
log_flux ~ Pile,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o,
method = "REML"
)
dat_co2e <- Flux_All_Gases_QC %>%
mutate(log_flux = log(GHG_CO2e_umol))%>%
arrange(Chamber_Corrected, DOY) %>%
group_by(Chamber_Corrected) %>%
mutate(t_index = row_number()) %>%
ungroup()
# applying model
m_co2e <- lme(
log_flux ~ Pile,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e,
method = "REML"
)
summary(m_co2)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 4.85541274 0.2502765 9273 19.40019048 3.231981e-82
## PileE 0.01324912 0.3567358 4 0.03713986 9.721531e-01summary(m_ch4)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 7.2578757 0.6867876 10286 10.5678613 5.704952e-26
## PileE -0.4425029 1.0106864 4 -0.4378241 6.841178e-01summary(m_n2o)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 4.0477206 0.2135944 13811 18.950502 4.396067e-79
## PileE -0.4668237 0.3152567 4 -1.480773 2.127801e-01summary(m_co2e)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 9.9295460 0.3194308 5552 31.0851212 1.003075e-195
## PileE -0.3951411 0.4905084 4 -0.8055747 4.656491e-01extract_pile_effect <- function(model, gas_name) {
tt <- summary(model)$tTable
beta <- tt["PileE", "Value"]
se <- tt["PileE", "Std.Error"]
pval <- tt["PileE", "p-value"]
data.frame(
Gas = gas_name,
Ratio_E_over_C = exp(beta),
CI_low = exp(beta - 1.96 * se),
CI_high = exp(beta + 1.96 * se),
p_value = pval
)
}
results_lme <- bind_rows(
extract_pile_effect(m_co2, "CO2"),
extract_pile_effect(m_ch4, "CH4"),
extract_pile_effect(m_n2o, "N2O"),
extract_pile_effect(m_co2e, "CO2e")
)
results_lme## Gas Ratio_E_over_C CI_low CI_high p_value
## 1 CO2 1.0133373 0.50361004 2.038983 0.9721531
## 2 CH4 0.6424265 0.08861552 4.657331 0.6841178
## 3 N2O 0.6269906 0.33799410 1.163089 0.2127801
## 4 CO2e 0.6735850 0.25755097 1.761658 0.4656491knitr::kable(
results_lme,
digits = 3,
caption = "Pile effect estimates from log-scale linear mixed-effects models."
)| Gas | Ratio_E_over_C | CI_low | CI_high | p_value |
|---|---|---|---|---|
| CO2 | 1.013 | 0.504 | 2.039 | 0.972 |
| CH4 | 0.642 | 0.089 | 4.657 | 0.684 |
| N2O | 0.627 | 0.338 | 1.163 | 0.213 |
| CO2e | 0.674 | 0.258 | 1.762 | 0.466 |
adding chamber elevation to the model as this impact flux magnitude and does not have to do with mechanism or treatment:
dat_co2 <- dat_co2 %>%
group_by(Chamber_Corrected) %>%
mutate(
Chamber_Elevation_filled = ifelse(
is.na(Chamber_Elevation),
median(Chamber_Elevation, na.rm = TRUE), # or mean(...)
Chamber_Elevation
)
) %>%
ungroup()
dat_co2 <- dat_co2 %>%
mutate(Chamber_Elevation_sc = as.numeric(scale(Chamber_Elevation_filled)))
m_co2_elev <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2,
method = "REML",
na.action = na.omit
)
dat_ch4 <- dat_ch4 %>%
group_by(Chamber_Corrected) %>%
mutate(
Chamber_Elevation_filled = ifelse(
is.na(Chamber_Elevation),
median(Chamber_Elevation, na.rm = TRUE), # or mean(...)
Chamber_Elevation
)
) %>%
ungroup()
dat_ch4 <- dat_ch4 %>%
mutate(Chamber_Elevation_sc = as.numeric(scale(Chamber_Elevation_filled)))
m_ch4_elev <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4,
method = "REML",
na.action = na.omit
)
dat_n2o <- dat_n2o %>%
group_by(Chamber_Corrected) %>%
mutate(
Chamber_Elevation_filled = ifelse(
is.na(Chamber_Elevation),
median(Chamber_Elevation, na.rm = TRUE), # or mean(...)
Chamber_Elevation
)
) %>%
ungroup()
dat_n2o <- dat_n2o %>%
mutate(Chamber_Elevation_sc = as.numeric(scale(Chamber_Elevation_filled)))
m_n2o_elev <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o,
method = "REML",
na.action = na.omit
)
dat_co2e <- dat_co2e %>%
group_by(Chamber_Corrected) %>%
mutate(
Chamber_Elevation_filled = ifelse(
is.na(Chamber_Elevation),
median(Chamber_Elevation, na.rm = TRUE), # or mean(...)
Chamber_Elevation
)
) %>%
ungroup()
dat_co2e <- dat_co2e %>%
mutate(Chamber_Elevation_sc = as.numeric(scale(Chamber_Elevation_filled)))
m_co2e_elev <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e,
method = "REML",
na.action = na.omit
)
results_lme_elv <- bind_rows(
extract_pile_effect(m_co2e_elev, "CO2"),
extract_pile_effect(m_ch4_elev, "CH4"),
extract_pile_effect(m_n2o_elev, "N2O"),
extract_pile_effect(m_co2e_elev, "CO2e")
)
results_lme_elv## Gas Ratio_E_over_C CI_low CI_high p_value
## 1 CO2 0.6089703 0.20520084 1.807229 0.4220016
## 2 CH4 0.6297099 0.06327159 6.267183 0.7133416
## 3 N2O 0.6643068 0.32098921 1.374824 0.3322370
## 4 CO2e 0.6089703 0.20520084 1.807229 0.4220016dat_co2_int <- dat_co2 %>%
mutate(Pile = factor(Pile, levels = c("C", "E"))) %>%
# remove invalid temperature values
mutate(
TS_1.initial_value = ifelse(
TS_1.initial_value <= 0,
NA,
TS_1.initial_value
)
) %>%
# fill missing values with the mean
mutate(
TS_1_filled = ifelse(
is.na(TS_1.initial_value),
mean(TS_1.initial_value, na.rm = TRUE),
TS_1.initial_value
)
) %>%
# scale AFTER filling
mutate(
pile_temp_sc = as.numeric(scale(TS_1_filled))
)
m_co2_temp_int <- lme(
log_flux ~ Pile * pile_temp_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2_int,
method = "REML",
na.action = na.omit
)
summary(m_co2_temp_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 4.86755246 0.26679958 9270 18.24422825 4.331843e-73
## PileE 0.02531670 0.37894125 4 0.06680903 9.499398e-01
## pile_temp_sc 0.22594883 0.02694158 9270 8.38662226 5.734823e-17
## Chamber_Elevation_sc 0.09157048 0.02282059 9270 4.01262472 6.051956e-05
## PileE:pile_temp_sc 0.14611501 0.04231313 9270 3.45318321 5.564882e-04m_base_ml <- update(m_co2_elev, method = "ML")
m_int_ml <- update(m_co2_temp_int, method = "ML")
anova(m_base_ml , m_int_ml)## Model df AIC BIC logLik Test L.Ratio p-value
## m_base_ml 1 7 5631.331 5681.279 -2808.665
## m_int_ml 2 9 5465.668 5529.888 -2723.834 1 vs 2 169.6624 <.0001# temperature range (scaled)
temp_seq <- seq(
min(dat_co2_int$pile_temp_sc),
max(dat_co2_int$pile_temp_sc),
length.out = 100
)
newdat <- expand.grid(
pile_temp_sc = temp_seq,
Pile = levels(dat_co2_int$Pile)
)
# hold elevation at its mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict population-level (fixed effects only)
newdat$pred_log <- predict(m_co2_temp_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_co2_int,
aes(x = pile_temp_sc,
y = exp(log_flux),
colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat,
aes(x = pile_temp_sc,
y = pred_flux,
colour = Pile),
linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile temperature (scaled)",
y = expression(CO[2]~flux~(back~transformed)),
colour = "Pile"
)
This was an interesting plot, lets do it for the other trace gases and
then move forward. At lower temps the piles converge but there is more
flux per increase in temp for the E pile
dat_ch4_int <- dat_ch4 %>%
mutate(Pile = factor(Pile, levels = c("C", "E"))) %>%
# remove invalid temperature values
mutate(
TS_1.initial_value = ifelse(
TS_1.initial_value <= 0,
NA,
TS_1.initial_value
)
) %>%
# fill missing values with the mean
mutate(
TS_1_filled = ifelse(
is.na(TS_1.initial_value),
mean(TS_1.initial_value, na.rm = TRUE),
TS_1.initial_value
)
) %>%
# scale AFTER filling
mutate(
pile_temp_sc = as.numeric(scale(TS_1_filled))
)
m_ch4_temp_int <- lme(
log_flux ~ Pile * pile_temp_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4_int,
method = "REML",
na.action = na.omit
)
summary(m_ch4_temp_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 7.23817774 0.78662856 10283 9.2015191 4.214351e-20
## PileE -0.38927899 1.14064600 4 -0.3412794 7.500671e-01
## pile_temp_sc 0.03203206 0.03755505 10283 0.8529361 3.937146e-01
## Chamber_Elevation_sc 0.26676320 0.03865930 10283 6.9003628 5.492551e-12
## PileE:pile_temp_sc 0.54664546 0.06184593 10283 8.8388265 1.125410e-18m_base_ml <- update(m_ch4_elev, method = "ML")
m_int_ml <- update(m_ch4_temp_int, method = "ML")
anova(m_base_ml , m_int_ml)## Model df AIC BIC logLik Test L.Ratio p-value
## m_base_ml 1 7 11026.02 11076.70 -5506.012
## m_int_ml 2 9 10894.13 10959.28 -5438.066 1 vs 2 135.8926 <.0001# temperature range (scaled)
temp_seq <- seq(
min(dat_ch4_int$pile_temp_sc),
max(dat_ch4_int$pile_temp_sc),
length.out = 100
)
newdat <- expand.grid(
pile_temp_sc = temp_seq,
Pile = levels(dat_ch4_int$Pile)
)
# hold elevation at its mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict population-level (fixed effects only)
newdat$pred_log <- predict(m_ch4_temp_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_ch4_int,
aes(x = pile_temp_sc,
y = exp(log_flux),
colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat,
aes(x = pile_temp_sc,
y = pred_flux,
colour = Pile),
linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile temperature (scaled)",
y = expression(Ch[4]~flux~(back~transformed)),
colour = "Pile"
)
For both ch4 and co2 the E pile has a stronger response to an increase in temp.
dat_n2o_int <- dat_n2o %>%
mutate(Pile = factor(Pile, levels = c("C", "E"))) %>%
# remove invalid temperature values
mutate(
TS_1.initial_value = ifelse(
TS_1.initial_value <= 0,
NA,
TS_1.initial_value
)
) %>%
# fill missing values with the mean
mutate(
TS_1_filled = ifelse(
is.na(TS_1.initial_value),
mean(TS_1.initial_value, na.rm = TRUE),
TS_1.initial_value
)
) %>%
# scale AFTER filling
mutate(
pile_temp_sc = as.numeric(scale(TS_1_filled))
)
m_n2o_temp_int <- lme(
log_flux ~ Pile * pile_temp_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o_int,
method = "REML",
na.action = na.omit
)
summary(m_n2o_temp_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 4.03288450 0.25385310 13808 15.8866862 2.471057e-56
## PileE -0.41692329 0.37023339 4 -1.1261094 3.231086e-01
## pile_temp_sc -0.01431110 0.03251010 13808 -0.4402047 6.597958e-01
## Chamber_Elevation_sc 0.25640388 0.03787607 13808 6.7695482 1.344010e-11
## PileE:pile_temp_sc -0.02935484 0.05247626 13808 -0.5593929 5.759027e-01m_base_ml <- update(m_n2o_elev, method = "ML")
m_int_ml <- update(m_n2o_temp_int, method = "ML")
anova(m_base_ml , m_int_ml)## Model df AIC BIC logLik Test L.Ratio p-value
## m_base_ml 1 7 10630.19 10682.92 -5308.093
## m_int_ml 2 9 10632.88 10700.68 -5307.438 1 vs 2 1.310477 0.5193# temperature range (scaled)
temp_seq <- seq(
min(dat_n2o_int$pile_temp_sc),
max(dat_n2o_int$pile_temp_sc),
length.out = 100
)
newdat <- expand.grid(
pile_temp_sc = temp_seq,
Pile = levels(dat_n2o_int$Pile)
)
# hold elevation at its mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict population-level (fixed effects only)
newdat$pred_log <- predict(m_n2o_temp_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_n2o_int,
aes(x = pile_temp_sc,
y = exp(log_flux),
colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat,
aes(x = pile_temp_sc,
y = pred_flux,
colour = Pile),
linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile temperature (scaled)",
y = expression(n2o~flux~(back~transformed)),
colour = "Pile"
)
dat_co2e_int <- dat_co2e %>%
mutate(Pile = factor(Pile, levels = c("C", "E"))) %>%
# remove invalid temperature values
mutate(
TS_1.initial_value = ifelse(
TS_1.initial_value <= 0,
NA,
TS_1.initial_value
)
) %>%
# fill missing values with the mean
mutate(
TS_1_filled = ifelse(
is.na(TS_1.initial_value),
mean(TS_1.initial_value, na.rm = TRUE),
TS_1.initial_value
)
) %>%
# scale AFTER filling
mutate(
pile_temp_sc = as.numeric(scale(TS_1_filled))
)
m_co2e_temp_int <- lme(
log_flux ~ Pile * pile_temp_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e_int,
method = "REML",
na.action = na.omit
)
summary(m_co2e_temp_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 9.9584957 0.41351778 5549 24.0823880 5.753037e-122
## PileE -0.5311861 0.61399085 4 -0.8651369 4.357648e-01
## pile_temp_sc 0.1202260 0.04005672 5549 3.0013927 2.699437e-03
## Chamber_Elevation_sc 0.2199742 0.03335744 5549 6.5944558 4.663712e-11
## PileE:pile_temp_sc -0.3739876 0.06384617 5549 -5.8576365 4.964783e-09m_base_ml <- update(m_co2e_elev, method = "ML")
m_int_ml <- update(m_co2e_temp_int, method = "ML")
anova(m_base_ml , m_int_ml)## Model df AIC BIC logLik Test L.Ratio p-value
## m_base_ml 1 7 4827.596 4873.957 -2406.798
## m_int_ml 2 9 4797.290 4856.897 -2389.645 1 vs 2 34.30643 <.0001# temperature range (scaled)
temp_seq <- seq(
min(dat_co2e_int$pile_temp_sc),
max(dat_co2e_int$pile_temp_sc),
length.out = 100
)
newdat <- expand.grid(
pile_temp_sc = temp_seq,
Pile = levels(dat_co2e_int$Pile)
)
# hold elevation at its mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict population-level (fixed effects only)
newdat$pred_log <- predict(m_co2e_temp_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_co2e_int,
aes(x = pile_temp_sc,
y = exp(log_flux),
colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat,
aes(x = pile_temp_sc,
y = pred_flux,
colour = Pile),
linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile temperature (scaled)",
y = expression(Co2e~flux~(back~transformed)),
colour = "Pile"
)
# ---- 1) Build modeling dataframe (CO2 + moisture) ----
dat_co2_moist_int <- dat_co2 %>%
mutate(
Pile = factor(Pile, levels = c("C", "E")),
# moisture raw (your chosen column)
moist_raw = SWC_1.initial_value,
# OPTIONAL: screen impossible moisture values (adjust if needed)
moist_raw = ifelse(moist_raw < 0.01 | moist_raw > 0.99, NA, moist_raw),
# mean-fill, then scale
moist_filled = ifelse(
is.na(moist_raw),
mean(moist_raw, na.rm = TRUE),
moist_raw
),
pile_moist_sc = as.numeric(scale(moist_filled))
)
# sanity checks
stopifnot(nrow(dat_co2_moist_int) > 0)
stopifnot(all(!is.na(dat_co2_moist_int$pile_moist_sc)))
stopifnot(!is.null(levels(dat_co2_moist_int$Pile)))
# ---- 2) Fit moisture interaction model (REML) ----
m_co2_moist_int <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2_moist_int,
method = "REML",
na.action = na.omit
)
summary(m_co2_moist_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 4.84217146 0.29294469 9270 16.52930249 1.649691e-60
## PileE -0.01294063 0.41710942 4 -0.03102454 9.767363e-01
## pile_moist_sc -0.09582767 0.02889467 9270 -3.31644826 9.152002e-04
## Chamber_Elevation_sc 0.07924892 0.02398781 9270 3.30371615 9.577405e-04
## PileE:pile_moist_sc 0.13491124 0.04839345 9270 2.78779942 5.317547e-03# ---- 3) ML refits for valid model comparison ----
m_co2_elev_ml <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2_moist_int, # IMPORTANT: same data object
method = "ML",
na.action = na.omit
)
m_co2_moist_int_ml <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2_moist_int, # same data object
method = "ML",
na.action = na.omit
)
anova(m_co2_elev_ml, m_co2_moist_int_ml)## Model df AIC BIC logLik Test L.Ratio p-value
## m_co2_elev_ml 1 7 5631.331 5681.279 -2808.665
## m_co2_moist_int_ml 2 9 5623.971 5688.190 -2802.985 1 vs 2 11.36004 0.0034# ---- 4) Prediction grid + plot ----
moist_seq <- seq(
min(dat_co2_moist_int$pile_moist_sc, na.rm = TRUE),
max(dat_co2_moist_int$pile_moist_sc, na.rm = TRUE),
length.out = 100
)
newdat <- expand.grid(
pile_moist_sc = moist_seq,
Pile = levels(dat_co2_moist_int$Pile)
)
# hold elevation at mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict fixed effects only
newdat$pred_log <- predict(m_co2_moist_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_co2_moist_int,
aes(x = pile_moist_sc, y = exp(log_flux), colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat, aes(y = pred_flux), linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile moisture (scaled)",
y = expression(CO[2]~flux~(back~transformed)),
colour = "Pile"
)
dat_ch4_moist_int <- dat_ch4 %>%
mutate(
Pile = factor(Pile, levels = c("C", "E")),
# moisture raw (your chosen column)
moist_raw = SWC_1.initial_value,
# OPTIONAL: screen impossible moisture values (adjust if needed)
moist_raw = ifelse(moist_raw < 0.01 | moist_raw > 0.99, NA, moist_raw),
# mean-fill, then scale
moist_filled = ifelse(
is.na(moist_raw),
mean(moist_raw, na.rm = TRUE),
moist_raw
),
pile_moist_sc = as.numeric(scale(moist_filled))
)
# sanity checks
stopifnot(nrow(dat_ch4_moist_int) > 0)
stopifnot(all(!is.na(dat_ch4_moist_int$pile_moist_sc)))
stopifnot(!is.null(levels(dat_ch4_moist_int$Pile)))
# ---- 2) Fit moisture interaction model (REML) ----
m_ch4_moist_int <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4_moist_int,
method = "REML",
na.action = na.omit
)
summary(m_ch4_moist_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 7.4061972 0.72914951 10283 10.1573094 3.995447e-24
## PileE -0.6500348 1.06818848 4 -0.6085394 5.756931e-01
## pile_moist_sc 0.3063050 0.03759737 10283 8.1469782 4.164514e-16
## Chamber_Elevation_sc 0.2132810 0.03879621 10283 5.4974713 3.944808e-08
## PileE:pile_moist_sc -0.2441262 0.07226154 10283 -3.3783704 7.318748e-04# ---- 3) ML refits for valid model comparison ----
m_ch4_elev_ml <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4_moist_int, # IMPORTANT: same data object
method = "ML",
na.action = na.omit
)
m_ch4_moist_int_ml <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4_moist_int, # same data object
method = "ML",
na.action = na.omit
)
anova(m_ch4_elev_ml, m_ch4_moist_int_ml)## Model df AIC BIC logLik Test L.Ratio p-value
## m_ch4_elev_ml 1 7 11026.02 11076.70 -5506.012
## m_ch4_moist_int_ml 2 9 10962.77 11027.93 -5472.387 1 vs 2 67.2504 <.0001# ---- 4) Prediction grid + plot ----
moist_seq <- seq(
min(dat_ch4_moist_int$pile_moist_sc, na.rm = TRUE),
max(dat_ch4_moist_int$pile_moist_sc, na.rm = TRUE),
length.out = 100
)
newdat <- expand.grid(
pile_moist_sc = moist_seq,
Pile = levels(dat_ch4_moist_int$Pile)
)
# hold elevation at mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict fixed effects only
newdat$pred_log <- predict(m_ch4_moist_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_ch4_moist_int,
aes(x = pile_moist_sc, y = exp(log_flux), colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat, aes(y = pred_flux), linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile moisture (scaled)",
y = expression(CH[4]~flux~(back~transformed)),
colour = "Pile"
)
dat_n2o_moist_int <- dat_n2o %>%
mutate(
Pile = factor(Pile, levels = c("C", "E")),
# moisture raw (your chosen column)
moist_raw = SWC_1.initial_value,
# OPTIONAL: screen impossible moisture values (adjust if needed)
moist_raw = ifelse(moist_raw < 0.01 | moist_raw > 0.99, NA, moist_raw),
# mean-fill, then scale
moist_filled = ifelse(
is.na(moist_raw),
mean(moist_raw, na.rm = TRUE),
moist_raw
),
pile_moist_sc = as.numeric(scale(moist_filled))
)
# sanity checks
stopifnot(nrow(dat_n2o_moist_int) > 0)
stopifnot(all(!is.na(dat_n2o_moist_int$pile_moist_sc)))
stopifnot(!is.null(levels(dat_n2o_moist_int$Pile)))
# ---- 2) Fit moisture interaction model (REML) ----
m_n2o_moist_int <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o_moist_int,
method = "REML",
na.action = na.omit
)
summary(m_n2o_moist_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 4.02060940 0.26031732 13808 15.4450322 2.275968e-53
## PileE -0.39629575 0.37916930 4 -1.0451683 3.549443e-01
## pile_moist_sc -0.02186211 0.03052206 13808 -0.7162724 4.738353e-01
## Chamber_Elevation_sc 0.26227703 0.03812119 13808 6.8800849 6.239340e-12
## PileE:pile_moist_sc 0.01957934 0.04830110 13808 0.4053601 6.852191e-01# ---- 3) ML refits for valid model comparison ----
m_n2o_elev_ml <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o_moist_int, # IMPORTANT: same data object
method = "ML",
na.action = na.omit
)
m_n2o_moist_int_ml <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o_moist_int, # same data object
method = "ML",
na.action = na.omit
)
anova(m_n2o_elev_ml, m_n2o_moist_int_ml)## Model df AIC BIC logLik Test L.Ratio
## m_n2o_elev_ml 1 7 10630.19 10682.92 -5308.093
## m_n2o_moist_int_ml 2 9 10633.75 10701.56 -5307.876 1 vs 2 0.4346655
## p-value
## m_n2o_elev_ml
## m_n2o_moist_int_ml 0.8047# ---- 4) Prediction grid + plot ----
moist_seq <- seq(
min(dat_n2o_moist_int$pile_moist_sc, na.rm = TRUE),
max(dat_n2o_moist_int$pile_moist_sc, na.rm = TRUE),
length.out = 100
)
newdat <- expand.grid(
pile_moist_sc = moist_seq,
Pile = levels(dat_n2o_moist_int$Pile)
)
# hold elevation at mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict fixed effects only
newdat$pred_log <- predict(m_n2o_moist_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_n2o_moist_int,
aes(x = pile_moist_sc, y = exp(log_flux), colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat, aes(y = pred_flux), linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile moisture (scaled)",
y = expression(n2o~flux~(back~transformed)),
colour = "Pile"
)
dat_co2e_moist_int <- dat_co2e %>%
mutate(
Pile = factor(Pile, levels = c("C", "E")),
# moisture raw (your chosen column)
moist_raw = SWC_1.initial_value,
# OPTIONAL: screen impossible moisture values (adjust if needed)
moist_raw = ifelse(moist_raw < 0.01 | moist_raw > 0.99, NA, moist_raw),
# mean-fill, then scale
moist_filled = ifelse(
is.na(moist_raw),
mean(moist_raw, na.rm = TRUE),
moist_raw
),
pile_moist_sc = as.numeric(scale(moist_filled))
)
# sanity checks
stopifnot(nrow(dat_co2e_moist_int) > 0)
stopifnot(all(!is.na(dat_co2e_moist_int$pile_moist_sc)))
stopifnot(!is.null(levels(dat_co2e_moist_int$Pile)))
# ---- 2) Fit moisture interaction model (REML) ----
m_co2e_moist_int <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e_moist_int,
method = "REML",
na.action = na.omit
)
summary(m_co2e_moist_int)$tTable## Value Std.Error DF t-value p-value
## (Intercept) 10.0183537 0.33351685 5549 30.0385233 8.157644e-184
## PileE -0.4570616 0.50581505 4 -0.9036142 4.172990e-01
## pile_moist_sc 0.1221855 0.04464400 5549 2.7368848 6.222153e-03
## Chamber_Elevation_sc 0.1935003 0.03342559 5549 5.7889862 7.469444e-09
## PileE:pile_moist_sc -0.2787880 0.07960032 5549 -3.5023475 4.648106e-04# ---- 3) ML refits for valid model comparison ----
m_co2e_elev_ml <- lme(
log_flux ~ Pile + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e_moist_int, # IMPORTANT: same data object
method = "ML",
na.action = na.omit
)
m_co2e_moist_int_ml <- lme(
log_flux ~ Pile * pile_moist_sc + Chamber_Elevation_sc,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e_moist_int, # same data object
method = "ML",
na.action = na.omit
)
anova(m_co2e_elev_ml, m_co2e_moist_int_ml)## Model df AIC BIC logLik Test L.Ratio
## m_co2e_elev_ml 1 7 4827.596 4873.957 -2406.798
## m_co2e_moist_int_ml 2 9 4818.227 4877.834 -2400.114 1 vs 2 13.36917
## p-value
## m_co2e_elev_ml
## m_co2e_moist_int_ml 0.0013# ---- 4) Prediction grid + plot ----
moist_seq <- seq(
min(dat_co2e_moist_int$pile_moist_sc, na.rm = TRUE),
max(dat_co2e_moist_int$pile_moist_sc, na.rm = TRUE),
length.out = 100
)
newdat <- expand.grid(
pile_moist_sc = moist_seq,
Pile = levels(dat_co2e_moist_int$Pile)
)
# hold elevation at mean (0 after scaling)
newdat$Chamber_Elevation_sc <- 0
# predict fixed effects only
newdat$pred_log <- predict(m_co2e_moist_int, newdata = newdat, level = 0)
newdat$pred_flux <- exp(newdat$pred_log)
ggplot(dat_co2e_moist_int,
aes(x = pile_moist_sc, y = exp(log_flux), colour = Pile)) +
geom_point(alpha = 0.15, size = 1) +
geom_line(data = newdat, aes(y = pred_flux), linewidth = 1.2) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank()) +
labs(
x = "Pile moisture (scaled)",
y = expression(Co2e~flux~(back~transformed)),
colour = "Pile"
)
Looking at impacts post turning -
dat_co2_pt <- dat_co2_int %>%
mutate(
PostTurn_0_1d = if_else(!is.na(DaysSinceLastTurn) &
DaysSinceLastTurn > 0 &
DaysSinceLastTurn <= 1,
1L, 0L)
)
count(dat_co2_pt, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 4723
## 2 C 1 341
## 3 E 0 4089
## 4 E 1 126dat_co2_pt_use <- dat_co2_pt %>%
filter(!is.na(FCO2_DRY.LIN), is.finite(FCO2_DRY.LIN))
count(dat_co2_pt_use, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 4723
## 2 C 1 341
## 3 E 0 4089
## 4 E 1 126summary_co2 <- dat_co2_pt_use %>%
group_by(Pile, PostTurn_0_1d) %>%
summarise(
n = n(),
mean = mean(FCO2_DRY.LIN),
sd = sd(FCO2_DRY.LIN),
se = sd / sqrt(n),
median = median(FCO2_DRY.LIN),
.groups = "drop"
)
summary_co2## # A tibble: 4 × 7
## Pile PostTurn_0_1d n mean sd se median
## <fct> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 C 0 4723 188. 179. 2.61 132.
## 2 C 1 341 217. 403. 21.8 94.0
## 3 E 0 4089 193. 195. 3.04 146.
## 4 E 1 126 164. 174. 15.5 130.co2_effect <- summary_co2 %>%
select(Pile, PostTurn_0_1d, mean, median) %>%
tidyr::pivot_wider(names_from = PostTurn_0_1d,
values_from = c(mean, median),
names_prefix = "PostTurn=") %>%
mutate(
mean_diff = `mean_PostTurn=1` - `mean_PostTurn=0`,
mean_ratio = `mean_PostTurn=1` / `mean_PostTurn=0`,
median_diff = `median_PostTurn=1` - `median_PostTurn=0`,
median_ratio= `median_PostTurn=1` / `median_PostTurn=0`
)
co2_effect## # A tibble: 2 × 9
## Pile `mean_PostTurn=0` `mean_PostTurn=1` `median_PostTurn=0`
## <fct> <dbl> <dbl> <dbl>
## 1 C 188. 217. 132.
## 2 E 193. 164. 146.
## # ℹ 5 more variables: `median_PostTurn=1` <dbl>, mean_diff <dbl>,
## # mean_ratio <dbl>, median_diff <dbl>, median_ratio <dbl>ggplot(dat_co2_pt_use,
aes(x = factor(PostTurn_0_1d),
y = FCO2_DRY.LIN,
fill = Pile)) +
geom_boxplot(outlier.alpha = 0.3) +
scale_x_discrete(labels = c("0" = "Not post-turn",
"1" = "0–1 d post-turn")) +
labs(x = "", y = "CO₂ flux") +
theme_bw()
m_co2_postturn <- lme(
fixed = log(FCO2_DRY.LIN) ~ Pile * PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2_pt_use,
method = "REML",
na.action = na.omit
)
m_co2_no_int <- lme(
fixed = log(FCO2_DRY.LIN) ~ Pile + PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2_pt_use,
method = "REML",
na.action = na.omit
)
anova(m_co2_no_int, m_co2_postturn)## Warning in anova.lme(m_co2_no_int, m_co2_postturn): fitted objects with
## different fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## m_co2_no_int 1 7 5634.897 5684.843 -2810.448
## m_co2_postturn 2 8 5610.898 5667.978 -2797.449 1 vs 2 25.99915 <.0001summary(m_co2_postturn)## Linear mixed-effects model fit by REML
## Data: dat_co2_pt_use
## AIC BIC logLik
## 5610.898 5667.978 -2797.449
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 0.4255526 0.7555908
##
## Correlation Structure: AR(1)
## Formula: ~t_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.917621
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Pile
## Parameter estimates:
## C E
## 1.000000 1.202184
## Fixed effects: log(FCO2_DRY.LIN) ~ Pile * PostTurn_0_1d
## Value Std.Error DF t-value p-value
## (Intercept) 4.853508 0.2509395 9271 19.341350 0.0000
## PileE -0.005522 0.3576634 4 -0.015440 0.9884
## PostTurn_0_1d 0.024753 0.0475013 9271 0.521105 0.6023
## PileE:PostTurn_0_1d 0.518734 0.0964963 9271 5.375683 0.0000
## Correlation:
## (Intr) PileE PT_0_1
## PileE -0.702
## PostTurn_0_1d -0.014 0.010
## PileE:PostTurn_0_1d 0.007 -0.013 -0.492
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.19017184 -0.59941032 0.07751574 0.70676900 3.28926220
##
## Number of Observations: 9279
## Number of Groups: 6So if we look at day after turning between the pile there is a huge difference that is statistically sig, if you change this to 3 days there is too much noise but the mean is still higher. This is something that can be used to say that there are differences in how the piles respond to turning based on the treatment. Not sure what to say about this, but we can decide what this is.
dat_ch4_pt <- dat_ch4_int %>%
mutate(
PostTurn_0_1d = if_else(!is.na(DaysSinceLastTurn) &
DaysSinceLastTurn > 0 &
DaysSinceLastTurn <= 1,
1L, 0L)
)
count(dat_ch4_pt, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 4957
## 2 C 1 398
## 3 E 0 4806
## 4 E 1 131dat_ch4_pt_use <- dat_ch4_pt %>%
filter(!is.na(FCH4_DRY.LIN), is.finite(FCH4_DRY.LIN))
count(dat_ch4_pt_use, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 4957
## 2 C 1 398
## 3 E 0 4806
## 4 E 1 131summary_ch4 <- dat_ch4_pt_use %>%
group_by(Pile, PostTurn_0_1d) %>%
summarise(
n = n(),
mean = mean(FCH4_DRY.LIN),
sd = sd(FCH4_DRY.LIN),
se = sd / sqrt(n),
median = median(FCH4_DRY.LIN),
.groups = "drop"
)
summary_ch4## # A tibble: 4 × 7
## Pile PostTurn_0_1d n mean sd se median
## <fct> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 C 0 4957 7101. 11525. 164. 1914.
## 2 C 1 398 6488. 14168. 710. 1377.
## 3 E 0 4806 6710. 10434. 151. 2527.
## 4 E 1 131 2362. 2198. 192. 1665.ch4_effect <- summary_ch4 %>%
select(Pile, PostTurn_0_1d, mean, median) %>%
tidyr::pivot_wider(names_from = PostTurn_0_1d,
values_from = c(mean, median),
names_prefix = "PostTurn=") %>%
mutate(
mean_diff = `mean_PostTurn=1` - `mean_PostTurn=0`,
mean_ratio = `mean_PostTurn=1` / `mean_PostTurn=0`,
median_diff = `median_PostTurn=1` - `median_PostTurn=0`,
median_ratio= `median_PostTurn=1` / `median_PostTurn=0`
)
ch4_effect## # A tibble: 2 × 9
## Pile `mean_PostTurn=0` `mean_PostTurn=1` `median_PostTurn=0`
## <fct> <dbl> <dbl> <dbl>
## 1 C 7101. 6488. 1914.
## 2 E 6710. 2362. 2527.
## # ℹ 5 more variables: `median_PostTurn=1` <dbl>, mean_diff <dbl>,
## # mean_ratio <dbl>, median_diff <dbl>, median_ratio <dbl>ggplot(dat_ch4_pt_use,
aes(x = factor(PostTurn_0_1d),
y = FCH4_DRY.LIN,
fill = Pile)) +
geom_boxplot(outlier.alpha = 0.3) +
scale_x_discrete(labels = c("0" = "Not post-turn",
"1" = "0–1 d post-turn")) +
labs(x = "", y = "CO₂ flux") +
theme_bw()
m_ch4_postturn <- lme(
fixed = log(FCH4_DRY.LIN) ~ Pile * PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4_pt_use,
method = "REML",
na.action = na.omit
)
m_ch4_no_int <- lme(
fixed = log(FCH4_DRY.LIN) ~ Pile + PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ch4_pt_use,
method = "REML",
na.action = na.omit
)
anova(m_ch4_no_int, m_ch4_postturn)## Warning in anova.lme(m_ch4_no_int, m_ch4_postturn): fitted objects with
## different fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## m_ch4_no_int 1 7 11069.28 11119.95 -5527.641
## m_ch4_postturn 2 8 11061.47 11119.38 -5522.733 1 vs 2 9.816498 0.0017summary(m_ch4_postturn)## Linear mixed-effects model fit by REML
## Data: dat_ch4_pt_use
## AIC BIC logLik
## 11061.47 11119.38 -5522.733
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 1.136322 1.575652
##
## Correlation Structure: AR(1)
## Formula: ~t_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.9783855
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Pile
## Parameter estimates:
## C E
## 1.000000 1.640851
## Fixed effects: log(FCH4_DRY.LIN) ~ Pile * PostTurn_0_1d
## Value Std.Error DF t-value p-value
## (Intercept) 7.262763 0.6869372 10284 10.572675 0.0000
## PileE -0.462141 1.0108066 4 -0.457200 0.6713
## PostTurn_0_1d -0.057435 0.0507006 10284 -1.132831 0.2573
## PileE:PostTurn_0_1d 0.484199 0.1402659 10284 3.452011 0.0006
## Correlation:
## (Intr) PileE PT_0_1
## PileE -0.680
## PostTurn_0_1d -0.006 0.004
## PileE:PostTurn_0_1d 0.002 -0.006 -0.361
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.0399901 -0.5922974 0.2570463 0.7412248 2.5679548
##
## Number of Observations: 10292
## Number of Groups: 6These values don’t make that much sense - we would expect post turn = 0 for the values to be similar with E having more ch4 then we have a much high mean for C after turning which I dont know what to make of that but we have a lower C pile median than E pile which makes me think there are huge spikes in the control and then it is lower where the E is just higher overall but not a huge spike ?
dat_n2o_pt <- dat_n2o_int %>%
mutate(
PostTurn_0_1d = if_else(!is.na(DaysSinceLastTurn) &
DaysSinceLastTurn > 0 &
DaysSinceLastTurn <= 1,
1L, 0L)
)
count(dat_n2o_pt, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 6397
## 2 C 1 509
## 3 E 0 6642
## 4 E 1 269dat_n2o_pt_use <- dat_n2o_pt %>%
filter(!is.na(FN2O_DRY.LIN), is.finite(FN2O_DRY.LIN))
count(dat_n2o_pt_use, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 6397
## 2 C 1 509
## 3 E 0 6642
## 4 E 1 269summary_n2o <- dat_n2o_pt_use %>%
group_by(Pile, PostTurn_0_1d) %>%
summarise(
n = n(),
mean = mean(FN2O_DRY.LIN),
sd = sd(FN2O_DRY.LIN),
se = sd / sqrt(n),
median = median(FN2O_DRY.LIN),
.groups = "drop"
)
summary_n2o## # A tibble: 4 × 7
## Pile PostTurn_0_1d n mean sd se median
## <fct> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 C 0 6397 154. 298. 3.73 61.5
## 2 C 1 509 200. 247. 10.9 111.
## 3 E 0 6642 96.2 247. 3.03 34.7
## 4 E 1 269 154. 195. 11.9 89.9n2o_effect <- summary_n2o %>%
select(Pile, PostTurn_0_1d, mean, median) %>%
tidyr::pivot_wider(names_from = PostTurn_0_1d,
values_from = c(mean, median),
names_prefix = "PostTurn=") %>%
mutate(
mean_diff = `mean_PostTurn=1` - `mean_PostTurn=0`,
mean_ratio = `mean_PostTurn=1` / `mean_PostTurn=0`,
median_diff = `median_PostTurn=1` - `median_PostTurn=0`,
median_ratio= `median_PostTurn=1` / `median_PostTurn=0`
)
n2o_effect## # A tibble: 2 × 9
## Pile `mean_PostTurn=0` `mean_PostTurn=1` `median_PostTurn=0`
## <fct> <dbl> <dbl> <dbl>
## 1 C 154. 200. 61.5
## 2 E 96.2 154. 34.7
## # ℹ 5 more variables: `median_PostTurn=1` <dbl>, mean_diff <dbl>,
## # mean_ratio <dbl>, median_diff <dbl>, median_ratio <dbl>ggplot(dat_n2o_pt_use,
aes(x = factor(PostTurn_0_1d),
y = FN2O_DRY.LIN,
fill = Pile)) +
geom_boxplot(outlier.alpha = 0.3) +
scale_x_discrete(labels = c("0" = "Not post-turn",
"1" = "0–1 d post-turn")) +
labs(x = "", y = "CO₂ flux") +
theme_bw()
m_n2o_postturn <- lme(
fixed = log(FN2O_DRY.LIN) ~ Pile * PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o_pt_use,
method = "REML",
na.action = na.omit
)
m_n2o_no_int <- lme(
fixed = log(FN2O_DRY.LIN) ~ Pile + PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_n2o_pt_use,
method = "REML",
na.action = na.omit
)
anova(m_n2o_no_int, m_n2o_postturn)## Warning in anova.lme(m_n2o_no_int, m_n2o_postturn): fitted objects with
## different fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## m_n2o_no_int 1 7 10525.12 10577.85 -5255.561
## m_n2o_postturn 2 8 10529.83 10590.10 -5256.914 1 vs 2 2.707689 0.0999summary(m_n2o_postturn)## Linear mixed-effects model fit by REML
## Data: dat_n2o_pt_use
## AIC BIC logLik
## 10529.83 10590.1 -5256.914
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 0.3133339 1.222624
##
## Correlation Structure: AR(1)
## Formula: ~t_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.9670161
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Pile
## Parameter estimates:
## C E
## 1.000000 1.289316
## Fixed effects: log(FN2O_DRY.LIN) ~ Pile * PostTurn_0_1d
## Value Std.Error DF t-value p-value
## (Intercept) 4.002188 0.21311569 13809 18.779415 0.0000
## PileE -0.447245 0.31504610 4 -1.419618 0.2287
## PostTurn_0_1d 0.537164 0.04795271 13809 11.201953 0.0000
## PileE:PostTurn_0_1d -0.022472 0.10048281 13809 -0.223644 0.8230
## Correlation:
## (Intr) PileE PT_0_1
## PileE -0.676
## PostTurn_0_1d -0.019 0.013
## PileE:PostTurn_0_1d 0.009 -0.019 -0.477
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.96298572 -0.67720660 0.05171918 0.69990626 3.22264388
##
## Number of Observations: 13817
## Number of Groups: 6Hmmmmm this one is weird so I’m not sure how to interpret it, but there is a stat difference p value of .09 at 1 day and .08 at 3 days and then at 7 days we see 0.00 something so unlike the other gases there is less noise and more a pronounced difference a week after turning unlike the other gases. So it seems n2o does not change immediately after turning, there is a delay vs co2 and ch4 where there is an immediate change, this might be mechanistic as it might be more microbe mediated vs built up release and aeration related ?
So there is a period after turning where n2o is different that is from like 2-3 days to about a week and then they are more similar, this pulse is very interesting.
dat_co2e_pt <- dat_co2e_int %>%
mutate(
PostTurn_0_1d = if_else(!is.na(DaysSinceLastTurn) &
DaysSinceLastTurn > 0 &
DaysSinceLastTurn <= 1,
1L, 0L)
)
count(dat_co2e_pt, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 3012
## 2 C 1 257
## 3 E 0 2201
## 4 E 1 88dat_co2e_pt_use <- dat_co2e_pt %>%
filter(!is.na(GHG_CO2e_umol), is.finite(GHG_CO2e_umol))
count(dat_co2e_pt_use, Pile, PostTurn_0_1d)## # A tibble: 4 × 3
## Pile PostTurn_0_1d n
## <fct> <int> <int>
## 1 C 0 3012
## 2 C 1 257
## 3 E 0 2201
## 4 E 1 88summary_co2e <- dat_co2e_pt_use %>%
group_by(Pile, PostTurn_0_1d) %>%
summarise(
n = n(),
mean = mean(GHG_CO2e_umol),
sd = sd(GHG_CO2e_umol),
se = sd / sqrt(n),
median = median(GHG_CO2e_umol),
.groups = "drop"
)
summary_co2e## # A tibble: 4 × 7
## Pile PostTurn_0_1d n mean sd se median
## <fct> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 C 0 3012 59676. 104318. 1901. 22688.
## 2 C 1 257 58331. 62977. 3928. 33204.
## 3 E 0 2201 39549. 97303. 2074. 16789.
## 4 E 1 88 52479. 61025. 6505. 34833.co2e_effect <- summary_co2e %>%
select(Pile, PostTurn_0_1d, mean, median) %>%
tidyr::pivot_wider(names_from = PostTurn_0_1d,
values_from = c(mean, median),
names_prefix = "PostTurn=") %>%
mutate(
mean_diff = `mean_PostTurn=1` - `mean_PostTurn=0`,
mean_ratio = `mean_PostTurn=1` / `mean_PostTurn=0`,
median_diff = `median_PostTurn=1` - `median_PostTurn=0`,
median_ratio= `median_PostTurn=1` / `median_PostTurn=0`
)
co2e_effect## # A tibble: 2 × 9
## Pile `mean_PostTurn=0` `mean_PostTurn=1` `median_PostTurn=0`
## <fct> <dbl> <dbl> <dbl>
## 1 C 59676. 58331. 22688.
## 2 E 39549. 52479. 16789.
## # ℹ 5 more variables: `median_PostTurn=1` <dbl>, mean_diff <dbl>,
## # mean_ratio <dbl>, median_diff <dbl>, median_ratio <dbl>ggplot(dat_co2e_pt_use,
aes(x = factor(PostTurn_0_1d),
y = GHG_CO2e_umol,
fill = Pile)) +
geom_boxplot(outlier.alpha = 0.3) +
scale_x_discrete(labels = c("0" = "Not post-turn",
"1" = "0–1 d post-turn")) +
labs(x = "", y = "CO₂e flux") +
theme_bw()
m_co2e_postturn <- lme(
fixed = log(GHG_CO2e_umol) ~ Pile * PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e_pt_use,
method = "REML",
na.action = na.omit
)
m_co2e_no_int <- lme(
fixed = log(GHG_CO2e_umol) ~ Pile + PostTurn_0_1d,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_co2e_pt_use,
method = "REML",
na.action = na.omit
)
anova(m_co2e_no_int, m_co2e_postturn)## Warning in anova.lme(m_co2e_no_int, m_co2e_postturn): fitted objects with
## different fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## m_co2e_no_int 1 7 4835.454 4881.811 -2410.727
## m_co2e_postturn 2 8 4822.706 4875.684 -2403.353 1 vs 2 14.7479 1e-04summary(m_co2e_postturn)## Linear mixed-effects model fit by REML
## Data: dat_co2e_pt_use
## AIC BIC logLik
## 4822.706 4875.684 -2403.353
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 0.4855702 1.182501
##
## Correlation Structure: AR(1)
## Formula: ~t_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.9616903
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Pile
## Parameter estimates:
## C E
## 1.000000 1.396841
## Fixed effects: log(GHG_CO2e_umol) ~ Pile * PostTurn_0_1d
## Value Std.Error DF t-value p-value
## (Intercept) 9.913977 0.3164929 5550 31.324482 0.0000
## PileE -0.414718 0.4862342 4 -0.852919 0.4418
## PostTurn_0_1d 0.210790 0.0533706 5550 3.949553 0.0001
## PileE:PostTurn_0_1d 0.536631 0.1300331 5550 4.126880 0.0000
## Correlation:
## (Intr) PileE PT_0_1
## PileE -0.651
## PostTurn_0_1d -0.013 0.008
## PileE:PostTurn_0_1d 0.005 -0.014 -0.410
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2912859 -0.6082548 0.1009931 0.6584155 3.2024992
##
## Number of Observations: 5558
## Number of Groups: 6It is significant for all three - which makes sense - sig goes down as I increase the time period we are looking at after turning.
Looking at ratios of gases:
While the means are not significantly different the ratio of trace gases might be, we would expect the exp pile to have more a higher ch4/co2 ratio compared to the more aerated c pile.
# -----------------------------
# 1) User-defined QC thresholds
# -----------------------------
r2_thresh <- 0.90 # example; set to whatever you use in your workflow
# ---------------------------------------------------
# 2) Build the ratio dataset with strict QC + units
# ---------------------------------------------------
# Key choices:
# - Only keep rows where BOTH CO2 and CH4 pass the linear-fit R² filter
# - Convert CH4 from nmol -> µmol so CO2 and CH4 are in comparable units
# - Keep only strictly positive fluxes for BOTH gases (required for log-ratios)
# - Create a per-chamber time index so AR(1) is defined correctly
dat_ratio <- Flux_All_Gases_QC %>%
# ---- QC filter: both gases must pass the linear model R² threshold ----
filter(
FCO2_DRY.LIN_R2 >= r2_thresh, # <-- replace with your actual CO2 R² column name
FCH4_DRY.LIN_R2 >= r2_thresh # <-- replace with your actual CH4 R² column name
) %>%
# ---- Units: put both gases in µmol m-2 s-1 ----
mutate(
CO2_umol = FCO2_DRY.LIN, # CO2 already in µmol m-2 s-1
CH4_umol = FCH4_DRY.LIN * 0.001 # CH4 in nmol -> µmol
) %>%
# ---- Ratio requires positive fluxes (log undefined at <= 0) ----
filter(
CO2_umol > 0,
CH4_umol > 0
) %>%
# ---- Define the ratio and its log-transform ----
# log_ratio is what we model because ratios are skewed and multiplicative.
mutate(
ratio_CO2_CH4 = CO2_umol / CH4_umol,
log_ratio = log(ratio_CO2_CH4)
) %>%
# ---- Time ordering + AR(1) index within chamber ----
arrange(Chamber_Corrected, DOY) %>%
group_by(Chamber_Corrected) %>%
mutate(t_index = row_number()) %>%
ungroup()
# Optional quick sanity checks
dat_ratio %>% summarise(
n_obs = n(),
n_chambers = n_distinct(Chamber_Corrected),
min_ratio = min(ratio_CO2_CH4),
max_ratio = max(ratio_CO2_CH4)
)## # A tibble: 1 × 4
## n_obs n_chambers min_ratio max_ratio
## <int> <int> <dbl> <dbl>
## 1 5558 6 2.23 115095.# -----------------------------------------
# 3) Fit the mixed model on the log-ratio
# -----------------------------------------
# Model interpretation:
# - Fixed effect (Pile): does the mean log(CO2/CH4) differ by pile?
# - Random intercept (Chamber_Corrected): each chamber has its own baseline ratio
# - corAR1: accounts for strong autocorrelation within chamber time-series
# - varIdent: allows different residual variance in each pile
m_ratio <- lme(
log_ratio ~ Pile,
random = ~ 1 | Chamber_Corrected,
correlation = corAR1(form = ~ t_index | Chamber_Corrected),
weights = varIdent(form = ~ 1 | Pile),
data = dat_ratio,
method = "REML"
)
summary(m_ratio)## Linear mixed-effects model fit by REML
## Data: dat_ratio
## AIC BIC logLik
## 5382.935 5422.671 -2685.467
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 1.052017 1.252931
##
## Correlation Structure: AR(1)
## Formula: ~t_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.9668886
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Pile
## Parameter estimates:
## C E
## 1.000000 1.633941
## Fixed effects: log_ratio ~ Pile
## Value Std.Error DF t-value p-value
## (Intercept) 4.728213 0.6300319 5552 7.504720 0.0000
## PileE 0.899644 0.9315278 4 0.965772 0.3888
## Correlation:
## (Intr)
## PileE -0.676
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.3740851 -0.6680709 -0.1878222 0.5465996 3.0176256
##
## Number of Observations: 5558
## Number of Groups: 6# ---------------------------------------------------------
# 4) Report pile-wise means on the RESPONSE SCALE (ratio)
# ---------------------------------------------------------
# emmeans avoids confusion about back-transformations and gives:
# - estimated mean ratio per pile
# - confidence intervals
# - correct degrees-of-freedom for the experimental design
# 4) Compare piles on the RESPONSE SCALE (ratio E/C)
# ============================================================
# Your model is fit on the log scale:
# log_ratio = log(CO2/CH4)
#
# The pile contrast from emmeans (E - C) is therefore on the log scale.
# If we exponentiate that log difference, we get a multiplicative ratio:
# exp( (log_ratio_E) - (log_ratio_C) ) = (CO2/CH4)_E / (CO2/CH4)_C
#
# This is the cleanest way to report: "Pile E has X times the CO2/CH4 ratio of pile C",
# along with a 95% confidence interval for that multiplicative factor.
# ============================================================
# 1) Estimated marginal means on the LINK (log) scale
emm_link <- emmeans(m_ratio, ~ Pile) # returns E[log(CO2/CH4)] by pile
# 2) Pairwise contrast on the LINK scale: (E - C) in log units
ctr_link <- pairs(emm_link)
# 3) Get the contrast estimate + 95% CI on the LINK scale (this guarantees CI columns exist)
ctr_link_ci <- confint(ctr_link) # returns estimate + lower.CL + upper.CL (and p.value)
ctr_link_ci## contrast estimate SE df lower.CL upper.CL
## C - E -0.9 0.932 4 -3.49 1.69
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95# This table is:
# estimate = (log_ratio_E - log_ratio_C)
# lower.CL, upper.CL = confidence interval for that difference in logs
# 4) Convert to RESPONSE SCALE: multiplicative ratio (E/C) with CI
ratio_E_over_C <- as.data.frame(ctr_link_ci) %>%
mutate(
ratio = exp(estimate),
ratio_low = exp(lower.CL),
ratio_high = exp(upper.CL)
)
ratio_E_over_C## contrast estimate SE df lower.CL upper.CL ratio ratio_low
## 1 C - E -0.8996436 0.9315278 4 -3.48598 1.686692 0.4067146 0.03062375
## ratio_high
## 1 5.401585# Interpretation:
# ratio > 1 => E has higher CO2/CH4 than C (less CH4 per unit CO2)
# ratio < 1 => E has lower CO2/CH4 than C (more CH4 per unit CO2)
# ------------------------------------------------------------
# Optional: also report pile means on the RESPONSE scale
# ------------------------------------------------------------
# These are the estimated mean CO2/CH4 ratios (not contrasts)
emm_resp <- emmeans(m_ratio, ~ Pile, type = "response")
emm_resp## Pile emmean SE df lower.CL upper.CL
## C 4.73 0.630 5 3.11 6.35
## E 5.63 0.686 4 3.72 7.53
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95confint(emm_resp)## Pile emmean SE df lower.CL upper.CL
## C 4.73 0.630 5 3.11 6.35
## E 5.63 0.686 4 3.72 7.53
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95# -----------------------------------------
# 5) Optional: quick diagnostic plots
# -----------------------------------------
plot(m_ratio) # residuals vs fitted
qqnorm(resid(m_ratio))
qqline(resid(m_ratio))
We do not have to use specific figures we can edit this, what Zaie has
done on the ratios is similiar. There is not a statistical difference in
the ratio of co2 to ch4 throughout the whole experimental cycle as there
just are not enough df and there is too much variance, what we do see if
a large effect size for the whole experiment.
Still A negative coefficient for PileE in a model of log(CO₂ / CH₄) means that the experimental pile has a lower CO₂/CH₄ ratio than the control pile,
Pile E has ~59% lower CO₂/CH₄ than pile C
Therefore:
More CH₄ per unit CO₂
Greater relative anaerobic (methanogenic) contribution
which in turn means more CH₄ emitted per unit CO₂ in the experimental pile. If you look at the end where the piles seem to stablize there are stat significant results but that is cherry picking times ?