CTD_26_07_2025
2025-07-26
This R markdown document is a draft of the analysis of 2023 emissions data from Churchtown Dairy’s composting windrows. A control pile and experimental pile were turned at different rates, and emissions/physical properties measured. The control pile was turned every 2 weeks and the experimental pile every 4 weeks during the summer and fall of 2023.
Research Question: How does turning frequency influence properties of compost windrows over time, and how do these properties explain trace gas emissions (e.g., CO₂, CH₄, N₂O) in both control and treatment piles?
To improve understanding of the underlying data here are a series of visual representations of physical and chemical properties of the composting windrows throughout the experimental cycle:
Figure 1: Rate of emissions of carbon dioxide throughout the experimental cycle - lines represent turning dates: Color of lines represent treatment:
Co2_plot_by_chamber <- ggplot(FCO2_Cleaned, aes(x= DOY, y = FCO2_DRY.LIN, group = Chamber_Corrected, color = Pile)) +
geom_point() +
# Vertical lines for sampling events
geom_vline(xintercept = c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines for turning events
geom_vline(xintercept = c(159, 186, 208, 234, 262, 290),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = " Carbon Dioxide Flux",
subtitle = "Chamber Fluxes by Treatment Group",
x = "Day of Year (DOY)",
y = expression(CO[2]~Flux~(umol~m^{-2}~s^{-1})),
color = "Turning Regime"
) +
# Theme customization
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Display the plot
print(Co2_plot_by_chamber)
Note: Read dashed lines are turning of both the experimental and control
piles. Black dashed lines represent turning only the control pile.
During turning events spikes in emissions are apparent.
Let’s look at this for each individual Chamber to see if we see similar trends across Chamber:
library(ggplot2)
Co2_plot_by_chamber <- ggplot(FCO2_Cleaned, aes(x = DOY, y = FCO2_DRY.LIN, color = Pile)) +
# Points per chamber
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
# Sampling events
geom_vline(xintercept = c(173, 191, 222, 249, 276),
linetype = "dashed", color = "black", size = 0.3) +
# Turning events
geom_vline(xintercept = c(159, 186, 208, 234, 262, 290),
linetype = "dashed", color = "red", size = 0.4) +
# Facet by chamber
facet_wrap(~ Chamber_Corrected, ncol = 2, scales = "free_y") +
# Labels
labs(
title = "Carbon Dioxide Flux",
subtitle = "Chamber Fluxes by Treatment Group",
x = "Day of Year (DOY)",
y = expression(CO[2]~Flux~(μmol~m^{-2}~s^{-1})), # using μmol instead of 'umol'
color = "Turning Regime"
) +
# Theme
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Co2_plot_by_chamber)
To improve interpretability of emissions data - weekly averages have also been plotted.
# creating weekly average:
FCO2_weekly <- FCO2_Cleaned %>%
mutate(Week = floor(DOY/7)) %>%
group_by(Pile, Chamber_Corrected, Week) %>%
summarise(
FCO2_weekly_mean = mean(FCO2_DRY.LIN, na.rm = TRUE),
.groups = "drop"
)library(ggplot2)
Co2_plot_weekly <- ggplot(FCO2_weekly, aes(x = Week, y = FCO2_weekly_mean, color = Pile)) +
# Points for weekly average per chamber
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
# Optional: Add smoother if desired later, e.g., geom_smooth(method = "loess", se = FALSE)
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Weekly-Averaged Carbon Dioxide Flux by Chamber",
subtitle = "Grouped by Treatment Regime",
x = "Week of Year",
y = expression(CO[2]~Flux~(μmol~m^{-2}~s^{-1})),
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Co2_plot_weekly)
Separating by chamber shows the spatial heterogeneity of the pile. To better understand larger trends between treatments not chambers, here is the average by treatment group, not chamber.
# creating weekly average:
FCO2_weekly_pile <- FCO2_Cleaned %>%
mutate(Week = floor(DOY/7)) %>%
group_by(Pile, Week) %>%
summarise(
FCO2_weekly_mean = mean(FCO2_DRY.LIN, na.rm = TRUE),
.groups = "drop"
)library(ggplot2)
Co2_plot_weekly_pile <- ggplot(FCO2_weekly_pile, aes(x = Week, y = FCO2_weekly_mean, color = Pile)) +
# Points for weekly average per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
# Optional: Add smoother if desired later, e.g., geom_smooth(method = "loess", se = FALSE)
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Weekly-Averaged Carbon Dioxide Flux by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "Week of Year",
y = expression(CO[2]~Flux~(μmol~m^{-2}~s^{-1})),
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Co2_plot_weekly_pile)
A line can also be added to show trends:
library(ggplot2)
Co2_plot_weekly_pile_line <- ggplot(FCO2_weekly_pile, aes(x = Week, y = FCO2_weekly_mean, color = Pile)) +
# Points for weekly average per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
geom_line(aes(group = Pile), linewidth = 1) +
# Optional: Add smoother if desired later, e.g., geom_smooth(method = "loess", se = FALSE)
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Weekly-Averaged Carbon Dioxide Flux by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "Week of Year",
y = expression(CO[2]~Flux~(μmol~m^{-2}~s^{-1})),
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Co2_plot_weekly_pile_line)
Moving on from Carbon Dioxide - the same graphs are created for each
trace gas (still need to clean these gases)
Here we see towards then end of the cycle a marked difference in emissions of co2 from the control and the experimental with the hypothetically more areated control pile emitting more co2 and if you look at the temp graph also warmer.
Lets move on to methane:
# this is a quick methane cleaning: not the actual cleaning - that will be done separately
# note these are only points that also have co2 - Im going to clean the other gases as there may be values that have poor co2 values but good ch4 values
Methane_data_for_plots <- FCO2_Cleaned %>%
filter(FCH4_DRY.LIN_R2 > 0.9 & FCH4_DRY.LIN > 0)
Methane_plot_points <- ggplot(Methane_data_for_plots, aes(x= DOY, y = FCH4_DRY.LIN, group = Chamber_Corrected, color = Pile)) +
geom_point() +
# Vertical lines for sampling events
geom_vline(xintercept = c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines for turning events
geom_vline(xintercept = c(159, 186, 208, 234, 262, 290),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = " Methane Flux",
subtitle = "Chamber Fluxes by Treatment Group",
x = "Day of Year (DOY)",
y = expression(Ch[4]~Flux~(nmol~m^{-2}~s^{-1})),
color = "Turning Regime"
) +
# Theme customization
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Methane_plot_points)
Similiar to co2 lets look at chambers individually
Ch4_plot_by_chamber <- ggplot(Methane_data_for_plots, aes(x = DOY, y = FCH4_DRY.LIN, color = Pile)) +
# Points per chamber
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
# Sampling events
geom_vline(xintercept = c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290),
linetype = "dashed", color = "black", size = 0.3) +
# Turning events
geom_vline(xintercept = c(159, 186, 208, 234, 262, 290),
linetype = "dashed", color = "red", size = 0.4) +
# Facet by chamber
facet_wrap(~ Chamber_Corrected, ncol = 2, scales = "free_y") +
# Labels
labs(
title = "Methane Flux",
subtitle = "Chamber Fluxes by Treatment Group",
x = "Day of Year (DOY)",
y = expression(Ch[4]~Flux~(nmol~m^{-2}~s^{-1})), # using μmol instead of 'umol'
color = "Turning Regime"
) +
# Theme
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Ch4_plot_by_chamber)
# creating weekly average by pile:
Weekly_ch4_Chambers_pile <- Methane_data_for_plots %>%
mutate(Week = floor(DOY/7)) %>%
group_by(Pile, Week) %>%
summarise(
weekly_ch4_by_chamber = mean(FCH4_DRY.LIN, na.rm = TRUE),
.groups = "drop"
)
Ch4_plot_weekly_pile_line <- ggplot(Weekly_ch4_Chambers_pile, aes(x = Week, y = weekly_ch4_by_chamber, color = Pile)) +
# Points for weekly average per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
geom_line(aes(group = Pile), linewidth = 1) +
# Optional: Add smoother if desired later, e.g., geom_smooth(method = "loess", se = FALSE)
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Weekly-Averaged Methane Flux by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "Week of Year",
y = expression(CH[4]~Flux~(nmol~m^{-2}~s^{-1})),
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Ch4_plot_weekly_pile_line)
Visually this is interesting, while the higher co2 values fit with the
hypothesis of more aeration with higher turning frequency - the high ch4
emissions do not as we would expect a lower ch4 emissions, this does not
appear to be the case. Overall the more turns the more emissions
visually, though that might not be significant.
Let’s move on to other variables, starting with temperature:
Temp_by_chamber_plot <- ggplot(FCO2_Cleaned, aes(x = DOY, y = TS_1.initial_value,group = Chamber_Corrected, color = Pile)) +
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Temperature under Chamber by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "DOY",
y = "Temp C",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Temp_by_chamber_plot)## Warning: Removed 1213 rows containing missing values or values outside the scale range
## (`geom_point()`).
Another version with seperate chambers:
library(ggplot2)
Temp_by_chamber_plot <- ggplot(FCO2_Cleaned, aes(x = DOY, y = TS_1.initial_value, color = Pile)) +
# Points by chamber
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
# Sampling events
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Turning events
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Facet by chamber
facet_wrap(~ Chamber_Corrected, ncol = 2, scales = "free_y") +
# Labels and formatting
labs(
title = "Chamber Temperature by Treatment",
subtitle = "Faceted by Chamber, Colored by Treatment",
x = "Day of Year (DOY)",
y = "Temperature (°C)",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Temp_by_chamber_plot)## Warning: Removed 1213 rows containing missing values or values outside the scale range
## (`geom_point()`).
As with other graphs - it is useful to simplify and average for ease of interpretability: Here we average the temperature by week for each chamber separately
# creating weekly average:
Weekly_Temp_Chambers <- FCO2_Cleaned %>%
mutate(Week = floor(DOY/7)) %>%
group_by(Pile,Chamber_Corrected, Week) %>%
summarise(
weekly_temp_by_chamber = mean(TS_1.initial_value, na.rm = TRUE),
.groups = "drop"
)library(ggplot2)
Temp_by_chamber_weekly <- ggplot(Weekly_Temp_Chambers, aes(x = Week, y = weekly_temp_by_chamber, group = Chamber_Corrected,color = Pile)) +
# Points for weekly average per chamber
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
# Optional: Add smoother if desired later, e.g., geom_smooth(method = "loess", se = FALSE)
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Weekly Average Temperature by Treatment (individual Chamber values)",
subtitle = "Grouped by Treatment Regime",
x = "Week of Year",
y = "Temp C",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Temp_by_chamber_weekly)## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
If we average by pile not chamber we get the following result:
# creating weekly average by pile:
Weekly_Temp_Chambers_pile <- FCO2_Cleaned %>%
mutate(Week = floor(DOY/7)) %>%
group_by(Pile, Week) %>%
summarise(
weekly_temp_by_chamber = mean(TS_1.initial_value, na.rm = TRUE),
.groups = "drop"
)library(ggplot2)
Temp_by_chamber_weekly_pile <- ggplot(Weekly_Temp_Chambers_pile, aes(x = Week, y = weekly_temp_by_chamber, group = Pile ,color = Pile)) +
# Points for weekly average per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
geom_line()+
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "Weekly Average Temperature by Treatment (Pile)",
subtitle = "Grouped by Treatment Regime",
x = "Week of Year",
y = "Temp C",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(Temp_by_chamber_weekly_pile)
Here a lower temp for the experimental pile fits with the lower
emissions seen both the ch4 and co2 data towards the end of the
experimental cycle.
Moving on to Soil Water Content measured in each Pile SWC
SWC_by_chamber_plot <- ggplot(FCO2_Cleaned, aes(x = DOY, y = SWC_1.initial_value ,group = Chamber_Corrected, color = Pile)) +
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "SWC under Chamber by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "DOY",
y = "[m+3m-3]",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(SWC_by_chamber_plot)## Warning: Removed 1213 rows containing missing values or values outside the scale range
## (`geom_point()`).
Something is up here - lets look at these individually to see if we can figure out what is going on. I took a look at the values and there are many that are .95 - I find it unlikely that so many values are exactly the same - I’m going to remove values of .95 as I believe they indicate an error with the sensor:
FCO2_Cleaned <- FCO2_Cleaned %>%
mutate(SWC_1.initial_value = ifelse(SWC_1.initial_value > 0.9, NA, SWC_1.initial_value))library(ggplot2)
SWC_by_chamber_plot <- ggplot(FCO2_Cleaned, aes(x = DOY, y = SWC_1.initial_value, color = Pile)) +
# Points per chamber
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
# Vertical lines: black for sampling
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines: red for turning
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Labels and formatting
labs(
title = "Soil Water Content (SWC) Under Each Chamber",
subtitle = "Faceted by Chamber, Colored by Treatment",
x = "Day of Year (DOY)",
y = "SWC (m³ m⁻³)",
color = "Turning Regime"
) +
# Facet by chamber
facet_wrap(~ Chamber_Corrected, ncol = 2) + # Try ncol = 2 or 3 for layout
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12), # chamber labels
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(SWC_by_chamber_plot)## Warning: Removed 3973 rows containing missing values or values outside the scale range
## (`geom_point()`).
There seems to be more variability in the control pile with the experimental pile having much higher SWC over time - with values close to 1 for large swathes of time. When values > .95 removed there is not much continous data for the experimental pile.
What if we look at FH20
FCO2_Cleaned_H2O <- FCO2_Cleaned %>%
mutate(FH2O = ifelse(FH2O < 0, NA, FH2O))FH2O_plot <- ggplot(FCO2_Cleaned_H2O, aes(x = DOY, y = FH2O ,group = Chamber_Corrected, color = Pile)) +
geom_point(aes(group = Chamber_Corrected), size = 2, alpha = 0.7) +
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "FH2O by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "DOY",
y = "[Moles h20 _ not sure units will check]",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(FH2O_plot)## Warning: Removed 5641 rows containing missing values or values outside the scale range
## (`geom_point()`).
Averaging by week/pile
# creating weekly average by pile:
FCO2_Cleaned_H2O_weekly <- FCO2_Cleaned_H2O %>%
mutate(Week = floor(DOY/7)) %>%
group_by(Pile, Week) %>%
summarise(
weekly_fh2o_by_chamber = mean(FH2O, na.rm = TRUE),
.groups = "drop"
)
FH2O_plot_weekly <- ggplot(FCO2_Cleaned_H2O_weekly, aes(x = Week, y = weekly_fh2o_by_chamber ,group = Pile, color = Pile)) +
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
# Vertical lines for events (convert DOY to weeks for alignment)
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290) / 7),
linetype = "dashed", color = "black", size = 0.3) +
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290) / 7),
linetype = "dashed", color = "red", size = 0.4) +
labs(
title = "FH20 weekly by Treatment",
subtitle = "Grouped by Treatment Regime",
x = "DOY",
y = "[Moles h20 _ not sure units will check]",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
print(FH2O_plot_weekly)## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
These values are hard to interpret - I have concerns over them.Lets look at the more discrete lab measured moisture content:
Moisture_data <- read.csv("C:/Users/KAUFMADA/Documents/CTD/Moisture_Content_2023.csv")Oven_moisture_data_plot <- ggplot(Moisture_data, aes(x = DOY, y = Moisture, color = Pile)) +
# Points per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
# Vertical lines: black for sampling
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines: red for turning
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Labels and formatting
labs(
title = "Soil Water Content (SWC) by Pile",
subtitle = " Colored by Treatment",
x = "Day of Year (DOY)",
y = "Moisture %",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12), # chamber labels
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Oven_moisture_data_plot)
It is hard to compare the two given the units are different - but we do NOT see the same trend of the Experimental pile having a higher moisture content as in the SWC sensor data. Further discussion of SWC is warranted - and cleaning needed this is something we can discuss as a group. The control pile which is turned more frequently having higher moisture content is not expected. The only explanation I can think of is turning redistributes the moisture, and the outer layer of the experimental pile is more dried out, yet as these measurements occured after turning, that would not explain this data.
Moving on to Bulk Density,Chamber elevation, 2d,3d surface area and volume.
Bulk_Density_Data <- read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_raw/BD_2023.csv")Bulk_Density_data_plot <- ggplot(Bulk_Density_Data, aes(x = DOY, y = grams.cm.3, color = Pile)) +
# Points per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
# Vertical lines: black for sampling
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines: red for turning
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Labels and formatting
labs(
title = "Bulk Density Pile",
subtitle = " Colored by Treatment",
x = "Day of Year (DOY)",
y = "grams per cm3",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12), # chamber labels
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Bulk_Density_data_plot)
Chamber Elevation Data:
chamber_elevation <- read.csv("C:/Users/KAUFMADA/Documents/R_Files/windrow-fluxes-co2-model-improvements/windrow-fluxes-co2-model-improvements/data_raw/Chamber_Elevation.csv")
chamber_elevation_after <- chamber_elevation %>%
filter(Turning_Status == "After")
Chamber_elevation_plot <- ggplot(chamber_elevation_after, aes(x = DOY, y = Chamber_Elevation, group = Chamber, color = Pile)) +
# Points per chamber
geom_point(aes(group = Chamber), size = 2, alpha = 0.7) +
# Vertical lines: black for sampling
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines: red for turning
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Labels and formatting
labs(
title = "Chamber Elevation",
subtitle = " Colored by Treatment",
x = "Day of Year (DOY)",
y = "Meters",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12), # chamber labels
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(Chamber_elevation_plot)
Pile Volume
surface_area_volume <- read.csv("C:/Users/KAUFMADA/Documents/CTD/2023 All Surface Area and Volume.csv")
volume_plot <- ggplot(surface_area_volume, aes(x = Date, y = Volume.m.3, color = Pile)) +
# Points per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
# Vertical lines: black for sampling
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines: red for turning
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Labels and formatting
labs(
title = "Pile Volume",
subtitle = " Colored by Treatment",
x = "Day of Year (DOY)",
y = "M^3",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12), # chamber labels
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(volume_plot)
Pile Surface Area - 2d
surface_area_plot <- ggplot(surface_area_volume, aes(x = Date, y = Area_2D.m.2, color = Pile)) +
# Points per chamber
geom_point(aes(group = Pile), size = 2, alpha = 0.7) +
# Vertical lines: black for sampling
geom_vline(xintercept = floor(c(159, 173, 186, 191, 208, 222, 234, 249, 262, 276, 290)),
linetype = "dashed", color = "black", size = 0.3) +
# Vertical lines: red for turning
geom_vline(xintercept = floor(c(159, 186, 208, 234, 262, 290)),
linetype = "dashed", color = "red", size = 0.4) +
# Labels and formatting
labs(
title = "Pile Area",
subtitle = " Colored by Treatment",
x = "Day of Year (DOY)",
y = "M^2",
color = "Turning Regime"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold"),
strip.text = element_text(size = 12), # chamber labels
legend.position = "top",
legend.title = element_text(face = "bold"),
legend.key.width = unit(1.5, "cm")
)
# Display the plot
print(surface_area_plot)
Research Question: How does turning frequency influence the physical and chemical properties of compost windrows over time, and how do these properties explain trace gas emissions (e.g., CO₂, CH₄, N₂O) in both control and treatment piles?
Lets begin by modeling emissions based on our known environmental drivers - which drivers explain the most variability in our measured values, what interactions between variables are most important. Does this differ by pile (treatment)?
# Plot the first histogram
hist(FCO2_Cleaned$FCO2_DRY.LIN, breaks = "FD",
main = "FCO2 Hist", xlab = "FCO2 umols/m^2/s")
Log transforming data
FCO2_Cleaned <- FCO2_Cleaned %>%
filter(FCO2_DRY.LIN > 0) %>% # Remove rows with 0 or negative values
mutate(log_FCO2_DRY.LIN = log(FCO2_DRY.LIN))# Plot the first histogram
hist(FCO2_Cleaned$log_FCO2_DRY.LIN, breaks = "FD",
main = "log FCO2 Hist", xlab = "FCO2 umols/m^2/s")
that looks way better
Scaling the variables
FCO2_Cleaned_scaled <- FCO2_Cleaned %>%
mutate(un_scaled_DOY = DOY) %>% # Save unscaled DOY for later use (e.g., autocorrelation)
mutate(across(where(is.numeric) & !matches("DOY"), scale))
# Step 2: Convert character variables to factors
FCO2_Cleaned_scaled_categorized <- FCO2_Cleaned_scaled %>%
mutate(across(where(is.character), as.factor))
# Step 3: Fill NA in numeric columns with median value
FCO2_Cleaned_scaled_categorized <- FCO2_Cleaned_scaled_categorized %>%
mutate(across(
where(is.numeric),
~ ifelse(is.na(.), median(., na.rm = TRUE), .)
))Variables scaled - median value used to fill NA:
Initial Full Model will include the following variables
Pile - C or E for control or experimental Chamber_Corrected -> chamber value (corrected for a switch in port between chamber 4-6) Chamber_Elevation -> average elevation of chamber used to correct for amount of material contributing to flux under the chamber air_temp -> air temp measured at eddy covar station half a mile away from pile - we also can use tmpa which is the air temp from local weather station TMPA -> air temp atlocal watchdog station VWCB -> soil moisture from sensor in watchdog, can tell us about soil moisture conditions. HMD -> humiditiy air_pressure -> bar from ec station BAR -> pressure from watchdog station TA.range -> the range in temp of the air measured in the long term chambers SWC_1.initial_value -> volumemetric water content measured by stevens probe TS_1.initial_value -> temp in pile from stevens probe at start of measurement FN2O_DRY.LIN -> lin flux of n2o FCH4_DRY.LIN -> lin flux of ch4 FH2O -> flux of h2o BulkDensity -> g/cm^3 DaysSinceLastTurn -> days since pile was turned TotalTurns -> total turns of the pile at each point in time
and of course Fco2_Dry.lin
# vars with interaction
interaction_terms <- c("TMPA", "SWC_1.initial_value", "TS_1.initial_value")
# All predictors to include
predictors <- c(
"TMPA", "VWCB", "HMD", "BAR", "TA.range",
"SWC_1.initial_value", "TS_1.initial_value",
"BulkDensity", "DaysSinceLastTurn", "TotalTurns", "Chamber_Elevation"
)# Variables to use in additive form only
additive_only <- setdiff(predictors, interaction_terms)
# Build interaction part
interaction_formula <- paste("Pile *", interaction_terms, collapse = " + ")
# Build additive part
additive_formula <- paste(additive_only, collapse = " + ")
# Combine into full formula
fixed_formula <- as.formula(
paste("log_FCO2_DRY.LIN ~", interaction_formula, "+", additive_formula)
)library(dplyr)
FCO2_Cleaned_scaled_categorized <- FCO2_Cleaned_scaled_categorized %>%
arrange(Chamber_Corrected, un_scaled_DOY) %>% # make sure it's ordered
group_by(Chamber_Corrected) %>%
mutate(time_index = row_number()) %>% # time index unique within each chamber
ungroup()
library(nlme)
full_model <- lme(
fixed = fixed_formula,
random = ~1 | Chamber_Corrected,
correlation = corAR1(form = ~ time_index | Chamber_Corrected),
data = FCO2_Cleaned_scaled_categorized,
na.action = na.exclude
)Full model has been run, lets look at the summary
summary(full_model)## Linear mixed-effects model fit by REML
## Data: FCO2_Cleaned_scaled_categorized
## AIC BIC logLik
## 7783.387 7918.929 -3872.694
##
## Random effects:
## Formula: ~1 | Chamber_Corrected
## (Intercept) Residual
## StdDev: 0.4958339 0.7916648
##
## Correlation Structure: AR(1)
## Formula: ~time_index | Chamber_Corrected
## Parameter estimate(s):
## Phi
## 0.8872216
## Fixed effects: list(fixed_formula)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00426348 0.2909507 9259 -0.014654 0.9883
## PileE -0.01342878 0.4140089 4 -0.032436 0.9757
## TMPA -0.03024323 0.0106116 9259 -2.850019 0.0044
## SWC_1.initial_value -0.04118356 0.0150080 9259 -2.744112 0.0061
## TS_1.initial_value 0.25882738 0.0370615 9259 6.983736 0.0000
## VWCB 0.01015890 0.0156065 9259 0.650939 0.5151
## HMD -0.05453176 0.0077048 9259 -7.077620 0.0000
## BAR -0.00611257 0.0031476 9259 -1.941957 0.0522
## TA.range 0.01771255 0.0087221 9259 2.030775 0.0423
## BulkDensity 0.02968563 0.0432064 9259 0.687066 0.4921
## DaysSinceLastTurn -0.15848449 0.0213511 9259 -7.422771 0.0000
## TotalTurns -0.20823218 0.0542780 9259 -3.836401 0.0001
## Chamber_Elevation 0.11701382 0.0334029 9259 3.503100 0.0005
## PileE:TMPA 0.03361953 0.0137148 9259 2.451332 0.0143
## PileE:SWC_1.initial_value 0.10323006 0.0285819 9259 3.611729 0.0003
## PileE:TS_1.initial_value 0.18416277 0.0510444 9259 3.607897 0.0003
## Correlation:
## (Intr) PileE TMPA SWC_1. TS_1._ VWCB HMD
## PileE -0.709
## TMPA -0.005 0.008
## SWC_1.initial_value 0.011 -0.007 -0.014
## TS_1.initial_value -0.022 0.029 -0.050 -0.012
## VWCB -0.001 0.002 -0.002 0.004 0.003
## HMD 0.003 -0.003 0.431 0.002 -0.082 -0.140
## BAR 0.002 -0.002 0.049 0.002 0.005 -0.031 0.143
## TA.range 0.003 -0.005 -0.015 0.000 -0.041 0.033 -0.088
## BulkDensity 0.072 -0.106 -0.013 0.009 0.002 -0.132 0.066
## DaysSinceLastTurn -0.007 0.012 0.030 0.043 0.072 -0.107 0.007
## TotalTurns -0.079 0.120 0.088 0.029 0.246 -0.140 0.005
## Chamber_Elevation 0.007 -0.007 -0.017 0.016 0.051 0.052 -0.042
## PileE:TMPA 0.006 -0.010 -0.622 0.012 0.015 -0.023 0.017
## PileE:SWC_1.initial_value -0.004 0.003 0.013 -0.521 0.025 0.004 0.002
## PileE:TS_1.initial_value 0.023 -0.025 0.040 0.016 -0.681 0.024 0.059
## BAR TA.rng BlkDns DysSLT TtlTrn Chmb_E PE:TMP
## PileE
## TMPA
## SWC_1.initial_value
## TS_1.initial_value
## VWCB
## HMD
## BAR
## TA.range -0.021
## BulkDensity 0.035 -0.015
## DaysSinceLastTurn 0.037 0.061 -0.103
## TotalTurns 0.008 0.016 -0.618 0.369
## Chamber_Elevation 0.004 0.027 0.068 0.139 0.002
## PileE:TMPA 0.126 -0.002 0.058 -0.029 -0.076 -0.004
## PileE:SWC_1.initial_value -0.030 -0.004 0.044 0.016 -0.003 0.020 -0.030
## PileE:TS_1.initial_value 0.039 -0.036 0.167 -0.049 -0.209 0.046 -0.057
## PE:SWC
## PileE
## TMPA
## SWC_1.initial_value
## TS_1.initial_value
## VWCB
## HMD
## BAR
## TA.range
## BulkDensity
## DaysSinceLastTurn
## TotalTurns
## Chamber_Elevation
## PileE:TMPA
## PileE:SWC_1.initial_value
## PileE:TS_1.initial_value 0.025
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.73409155 -0.54473886 0.09249508 0.62959221 3.77562209
##
## Number of Observations: 9279
## Number of Groups: 6Looking at the significant P>0,05 predictors + interactions and un loging them we can see there effect size
coefs <- summary(full_model)$tTable
coef_df <- as.data.frame(coefs)
coef_df$Variable <- rownames(coef_df)
coef_df$exp_Estimate <- exp(coef_df$Value)
coef_df$percent_change <- (coef_df$exp_Estimate - 1) * 100
coef_df %>%
select(Variable, Value, percent_change, `p-value`)## Variable Value percent_change
## (Intercept) (Intercept) -0.004263481 -0.4254405
## PileE PileE -0.013428778 -1.3339015
## TMPA TMPA -0.030243233 -2.9790482
## SWC_1.initial_value SWC_1.initial_value -0.041183563 -4.0347043
## TS_1.initial_value TS_1.initial_value 0.258827384 29.5410177
## VWCB VWCB 0.010158900 1.0210676
## HMD HMD -0.054531761 -5.3071567
## BAR BAR -0.006112571 -0.6093927
## TA.range TA.range 0.017712551 1.7870349
## BulkDensity BulkDensity 0.029685635 3.0130646
## DaysSinceLastTurn DaysSinceLastTurn -0.158484494 -14.6563803
## TotalTurns TotalTurns -0.208232177 -18.7981517
## Chamber_Elevation Chamber_Elevation 0.117013819 12.4134963
## PileE:TMPA PileE:TMPA 0.033619532 3.4191056
## PileE:SWC_1.initial_value PileE:SWC_1.initial_value 0.103230056 10.8746453
## PileE:TS_1.initial_value PileE:TS_1.initial_value 0.184162773 20.2211495
## p-value
## (Intercept) 9.883088e-01
## PileE 9.756784e-01
## TMPA 4.381312e-03
## SWC_1.initial_value 6.079179e-03
## TS_1.initial_value 3.072295e-12
## VWCB 5.151020e-01
## HMD 1.573129e-12
## BAR 5.217265e-02
## TA.range 4.230633e-02
## BulkDensity 4.920586e-01
## DaysSinceLastTurn 1.248192e-13
## TotalTurns 1.256796e-04
## Chamber_Elevation 4.620481e-04
## PileE:TMPA 1.425120e-02
## PileE:SWC_1.initial_value 3.057724e-04
## PileE:TS_1.initial_value 3.103156e-04Here we can see the ” effect size of different predictors and interactions”
The largest predictor by some margin was temperature of the pile followed by Days since last turning, and total turns (though this is going to be autocorrelated with DOY and time), Then Chamber elevation
Interacion Terms
PileE TMPA - 3 percent PileE SWC_1 10 - the effect of swc was 10 percent higher on co2 flux in the e pile PileE TS_1 20 - the effect of temp was 20 percent higher on co2 flux in the e pile
library(performance)
print(r2(full_model))## # R2 for Mixed Models
##
## Conditional R2: 0.415
## Marginal R2: 0.186Most of the variance is just between the chambers with a conditional r^2 of 0.186 Our environmental drivers do not explain most of the variance in emissions.
library(ggeffects)
ggpredict(full_model, terms = c("TMPA [all]", "Pile")) ## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## TMPA | Predicted | 95% CI
## -------------------------------
## -5.46 | 0.17 | -0.41, 0.75
## -1.10 | 0.04 | -0.53, 0.61
## -0.41 | 0.02 | -0.55, 0.59
## 0.26 | 0.00 | -0.57, 0.57
## 0.94 | -0.02 | -0.59, 0.55
## 3.86 | -0.11 | -0.69, 0.46
##
## Pile: E
##
## TMPA | Predicted | 95% CI
## -------------------------------
## -5.46 | -0.04 | -0.62, 0.55
## -1.10 | -0.02 | -0.60, 0.55
## -0.41 | -0.02 | -0.59, 0.55
## 0.26 | -0.02 | -0.59, 0.55
## 0.94 | -0.02 | -0.59, 0.56
## 3.86 | -0.01 | -0.59, 0.57
##
## Adjusted for:
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("VWCB [all]", "Pile")) # Soil moisture## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## VWCB | Predicted | 95% CI
## -------------------------------
## -0.92 | -0.01 | -0.58, 0.57
## -0.45 | 0.00 | -0.57, 0.57
## 0.06 | 0.00 | -0.57, 0.57
## 0.59 | 0.01 | -0.56, 0.58
## 1.27 | 0.02 | -0.55, 0.59
## 2.88 | 0.03 | -0.54, 0.61
##
## Pile: E
##
## VWCB | Predicted | 95% CI
## -------------------------------
## -0.92 | -0.03 | -0.60, 0.54
## -0.45 | -0.02 | -0.60, 0.55
## 0.06 | -0.02 | -0.59, 0.55
## 0.59 | -0.01 | -0.59, 0.56
## 1.27 | -0.01 | -0.58, 0.57
## 2.88 | 0.01 | -0.57, 0.59
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("HMD [all]", "Pile")) # Humidity## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## HMD | Predicted | 95% CI
## -----------------------------------
## -3.04 | 0.17 | -0.40, 0.74
## -2.15 | 0.12 | -0.45, 0.69
## -1.43 | 0.08 | -0.49, 0.65
## -0.71 | 0.04 | -0.53, 0.61
## -1.00e-03 | 0.00 | -0.57, 0.57
## 1.43 | -0.07 | -0.64, 0.50
##
## Pile: E
##
## HMD | Predicted | 95% CI
## -----------------------------------
## -3.04 | 0.15 | -0.43, 0.72
## -2.15 | 0.10 | -0.48, 0.67
## -1.43 | 0.06 | -0.51, 0.63
## -0.71 | 0.02 | -0.55, 0.59
## -1.00e-03 | -0.02 | -0.59, 0.55
## 1.43 | -0.10 | -0.67, 0.47
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("BAR [all]", "Pile")) # Pressure## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## BAR | Predicted | 95% CI
## -------------------------------
## -1.48 | 0.01 | -0.56, 0.58
## -1.07 | 0.01 | -0.56, 0.58
## -0.68 | 0.01 | -0.56, 0.58
## -0.27 | 0.01 | -0.56, 0.58
## 0.14 | 0.00 | -0.57, 0.57
## 87.70 | -0.53 | -1.32, 0.25
##
## Pile: E
##
## BAR | Predicted | 95% CI
## -------------------------------
## -1.48 | -0.01 | -0.58, 0.56
## -1.07 | -0.01 | -0.59, 0.56
## -0.68 | -0.02 | -0.59, 0.56
## -0.27 | -0.02 | -0.59, 0.55
## 0.14 | -0.02 | -0.59, 0.55
## 87.70 | -0.56 | -1.34, 0.23
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("TA.range [all]", "Pile")) # Temp range## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## TA.range | Predicted | 95% CI
## ----------------------------------
## -1.16 | -0.02 | -0.59, 0.55
## -0.36 | 0.00 | -0.57, 0.57
## 0.43 | 0.01 | -0.56, 0.58
## 1.23 | 0.03 | -0.54, 0.60
## 2.12 | 0.04 | -0.53, 0.61
## 7.54 | 0.14 | -0.45, 0.72
##
## Pile: E
##
## TA.range | Predicted | 95% CI
## ----------------------------------
## -1.16 | -0.04 | -0.61, 0.53
## -0.36 | -0.03 | -0.60, 0.55
## 0.43 | -0.01 | -0.58, 0.56
## 1.23 | 0.00 | -0.57, 0.57
## 2.12 | 0.02 | -0.56, 0.59
## 7.54 | 0.11 | -0.47, 0.70
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("SWC_1.initial_value [all]", "Pile")) # Water content## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## SWC_1.initial_value | Predicted | 95% CI
## ---------------------------------------------
## -1.47 | 0.06 | -0.51, 0.63
## -0.83 | 0.03 | -0.54, 0.60
## -0.41 | 0.02 | -0.55, 0.59
## 0.11 | -0.01 | -0.58, 0.56
## 0.93 | -0.04 | -0.61, 0.53
## 2.13 | -0.09 | -0.66, 0.49
##
## Pile: E
##
## SWC_1.initial_value | Predicted | 95% CI
## ---------------------------------------------
## -1.47 | -0.10 | -0.68, 0.47
## -0.83 | -0.06 | -0.64, 0.51
## -0.41 | -0.04 | -0.61, 0.54
## 0.11 | 0.00 | -0.58, 0.57
## 0.93 | 0.05 | -0.53, 0.62
## 2.13 | 0.12 | -0.46, 0.70
##
## Adjusted for:
## * TMPA = -0.00
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("TS_1.initial_value [all]", "Pile")) # Internal pile temp## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## TS_1.initial_value | Predicted | 95% CI
## ---------------------------------------------
## -2.94 | -0.76 | -1.37, -0.15
## -2.38 | -0.62 | -1.21, -0.02
## -1.78 | -0.46 | -1.05, 0.13
## -0.87 | -0.22 | -0.80, 0.35
## -0.03 | -0.01 | -0.58, 0.56
## 2.29 | 0.59 | 0.00, 1.18
##
## Pile: E
##
## TS_1.initial_value | Predicted | 95% CI
## ---------------------------------------------
## -2.94 | -1.33 | -1.94, -0.72
## -2.38 | -1.08 | -1.68, -0.48
## -1.78 | -0.81 | -1.40, -0.23
## -0.87 | -0.41 | -0.99, 0.17
## -0.03 | -0.04 | -0.61, 0.53
## 2.29 | 0.99 | 0.39, 1.58
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("BulkDensity [all]", "Pile")) # Density## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## BulkDensity | Predicted | 95% CI
## -------------------------------------
## -2.32 | -0.06 | -0.65, 0.52
## -1.56 | -0.04 | -0.62, 0.53
## -0.61 | -0.01 | -0.58, 0.55
## -0.13 | 0.00 | -0.57, 0.57
## 0.20 | 0.01 | -0.56, 0.58
## 2.45 | 0.08 | -0.54, 0.70
##
## Pile: E
##
## BulkDensity | Predicted | 95% CI
## -------------------------------------
## -2.32 | -0.09 | -0.71, 0.53
## -1.56 | -0.07 | -0.66, 0.53
## -0.61 | -0.04 | -0.62, 0.54
## -0.13 | -0.02 | -0.60, 0.55
## 0.20 | -0.01 | -0.59, 0.56
## 2.45 | 0.05 | -0.54, 0.65
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("DaysSinceLastTurn [all]", "Pile"))## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## DaysSinceLastTurn | Predicted | 95% CI
## -------------------------------------------
## -1.36 | 0.22 | -0.35, 0.79
## -0.82 | 0.13 | -0.44, 0.71
## -0.28 | 0.05 | -0.52, 0.62
## 0.34 | -0.05 | -0.62, 0.52
## 0.89 | -0.14 | -0.71, 0.43
## 3.22 | -0.51 | -1.09, 0.08
##
## Pile: E
##
## DaysSinceLastTurn | Predicted | 95% CI
## -------------------------------------------
## -1.36 | 0.20 | -0.38, 0.77
## -0.82 | 0.11 | -0.46, 0.68
## -0.28 | 0.02 | -0.55, 0.60
## 0.34 | -0.07 | -0.65, 0.50
## 0.89 | -0.16 | -0.73, 0.41
## 3.22 | -0.53 | -1.12, 0.06
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * TotalTurns = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("TotalTurns [all]", "Pile"))## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## TotalTurns | Predicted | 95% CI
## ------------------------------------
## -1.38 | 0.29 | -0.31, 0.89
## -0.99 | 0.21 | -0.38, 0.80
## -0.21 | 0.05 | -0.52, 0.62
## 0.58 | -0.12 | -0.69, 0.45
## 0.97 | -0.20 | -0.77, 0.37
## 2.15 | -0.44 | -1.04, 0.15
##
## Pile: E
##
## TotalTurns | Predicted | 95% CI
## ------------------------------------
## -1.38 | 0.27 | -0.31, 0.85
## -0.99 | 0.19 | -0.39, 0.76
## -0.21 | 0.02 | -0.55, 0.59
## 0.58 | -0.14 | -0.72, 0.44
## 0.97 | -0.22 | -0.81, 0.37
## 2.15 | -0.47 | -1.10, 0.17
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * Chamber_Elevation = -0.01
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggpredict(full_model, terms = c("Chamber_Elevation [all]", "Pile"))## # Predicted values of log_FCO2_DRY.LIN
##
## Pile: C
##
## Chamber_Elevation | Predicted | 95% CI
## -------------------------------------------
## -2.45 | -0.28 | -0.87, 0.31
## -0.85 | -0.09 | -0.67, 0.48
## -0.18 | -0.02 | -0.59, 0.55
## 0.23 | 0.03 | -0.54, 0.60
## 1.13 | 0.14 | -0.44, 0.71
## 2.72 | 0.32 | -0.27, 0.92
##
## Pile: E
##
## Chamber_Elevation | Predicted | 95% CI
## -------------------------------------------
## -2.45 | -0.31 | -0.90, 0.29
## -0.85 | -0.12 | -0.69, 0.46
## -0.18 | -0.04 | -0.61, 0.53
## 0.23 | 0.01 | -0.56, 0.58
## 1.13 | 0.11 | -0.46, 0.69
## 2.72 | 0.30 | -0.30, 0.90
##
## Adjusted for:
## * TMPA = -0.00
## * SWC_1.initial_value = -0.13
## * TS_1.initial_value = 0.02
## * VWCB = -0.00
## * HMD = -0.00
## * BAR = 0.00
## * TA.range = -0.00
## * BulkDensity = 0.00
## * DaysSinceLastTurn = 0.00
## * TotalTurns = 0.00
## * Chamber_Corrected = C1
## * time_index = 774.00##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.library(MuMIn)##
## Attaching package: 'MuMIn'## The following object is masked from 'package:randomForest':
##
## importancer.squaredGLMM(full_model)## R2m R2c
## [1,] 0.185679 0.4151145Very similar to the above values: