input to be able to knit

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 errors
library(lubridate)
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.5.2
library(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':
## 
##     stamp
library(rlang)
## 
## Attaching package: 'rlang'
## 
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
library(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_rows
library(nlme)
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(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'
library(gt)
## Warning: package 'gt' was built under R version 4.5.3
## 
## Attaching package: 'gt'
## 
## The following object is masked from 'package:cowplot':
## 
##     as_gtable

This first chunk (chamber maturation stage figure) - most of this is not needed - just needed it to be able to knit

# Add temperature zone to raw data
temp_raw <- Cleaned_Data %>%
  filter(!is.na(pile_temp)) %>%
  mutate(
    Pile_label = ifelse(Pile == "C", "Control (2×/month)", "Experimental (1×/month)"),
    temp_zone = case_when(
      pile_temp < 40                      ~ "Mesophilic",
      pile_temp >= 40 & pile_temp < 45    ~ "Transition",
      pile_temp >= 45 & pile_temp <= 70   ~ "Thermophilic",
      pile_temp > 70                      ~ "Hyper-thermophilic"
    ),
    temp_zone = factor(temp_zone, levels = c("Mesophilic", "Transition",
                                             "Thermophilic", "Hyper-thermophilic"))
  )

# Daily pile means with zone
temp_pile_means <- daily_pile %>%
  filter(valid) %>%
  mutate(
    Pile_label = ifelse(Pile == "Control", "Control (2×/month)", "Experimental (1×/month)")
  )

zone_cols <- c(
  "Mesophilic"         = "#2166ac",
  "Transition"         = "#74add1",
  "Thermophilic"       = "#f46d43",
  "Hyper-thermophilic" = "#a50026"
)

turns_both <- c(159, 186, 208, 234, 262, 290)
turns_control_only <- c(173, 191, 222, 249, 276)

# Separate turning lines for each panel
turns_control <- data.frame(
  xint = c(turns_both, turns_control_only),
  ltype = c(rep("solid", length(turns_both)), rep("solid", length(turns_control_only))),
  Pile_label = "Control (2×/month)"
)

turns_experimental <- data.frame(
  xint = turns_both,
  ltype = rep("solid", length(turns_both)),
  Pile_label = "Experimental (1×/month)"
)

turns_df <- bind_rows(turns_control, turns_experimental)

# Zone labels (right side of each panel)
zone_labels <- expand.grid(
  Pile_label = c("Control (2×/month)", "Experimental (1×/month)"),
  stringsAsFactors = FALSE
) %>%
  cross_join(data.frame(
    y = c(30, 42.5, 57.5, 75),
    label = c("Mesophilic", "Transition", "Thermophilic", "Hyper-\nthermophilic"),
    colour = c("#2166ac", "#74add1", "#f46d43", "#a50026")
  ))


make_pile_plot <- function(pile_filter, turn_doys, title_txt) {
  ggplot() +
    geom_vline(xintercept = turn_doys, colour = "grey70", linewidth = 0.2) +
# dont need raw for this one    geom_point(data = temp_raw %>% filter(Pile_label == pile_filter),
#               aes(x = DOY.initial_value, y = pile_temp, colour = temp_zone),
#               alpha = 0.08, size = 0.4) +
    geom_smooth(data = temp_pile_means %>% filter(Pile_label == pile_filter),
                aes(x = DOY_day, y = pile_daily_mean),
                method = "loess", se = FALSE, span = 0.4,
                linewidth = 0.8, colour = "black") +
    geom_point(data = temp_pile_means %>% filter(Pile_label == pile_filter),
               aes(x = DOY_day, y = pile_daily_mean, colour = temp_zone),
               shape = 17, size = 2.5) +
    geom_hline(yintercept = c(40, 45, 70),
               linetype = "dotted", colour = "grey60", linewidth = 0.3) +
    annotate("text", x = 302, y = 75,   label = "Hyper-therm.", size = 2.5, hjust = 0, colour = "#a50026", fontface = "bold") +
    annotate("text", x = 302, y = 57.5, label = "Thermophilic", size = 2.5, hjust = 0, colour = "#f46d43", fontface = "bold") +
    annotate("text", x = 302, y = 42.5, label = "Transition",   size = 2.5, hjust = 0, colour = "#74add1", fontface = "bold") +
    annotate("text", x = 302, y = 30,   label = "Mesophilic",   size = 2.5, hjust = 0, colour = "#2166ac", fontface = "bold") +
    scale_colour_manual(values = zone_cols, drop = FALSE, guide = "none") +
    coord_cartesian(xlim = c(150, 315), clip = "off") +
    labs(title = title_txt, x = "Day of Year",
         y = expression("Temperature ("*degree*"C)")) +
    theme_classic() +
    theme(legend.position = "none",
          plot.title = element_text(size = 11, face = "bold"),
          plot.margin = margin(5, 55, 5, 5))
}

fig_pile_ctrl <- make_pile_plot("Control (2×/month)",      c(turns_both, turns_control_only), "Control (2×/month)")
fig_pile_exp  <- make_pile_plot("Experimental (1×/month)", turns_both,                         "Experimental (1×/month)")

left_col  <- chamber_plots[["C1"]] / chamber_plots[["C2"]] / chamber_plots[["C3"]] / fig_pile_ctrl
right_col <- chamber_plots[["C4"]] / chamber_plots[["C5"]] / chamber_plots[["C6"]] / fig_pile_exp

final_fig <- (left_col | right_col)
plot(final_fig)
## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 51 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'

# PDF — THIS IS KEY TO SAVE AS PDF OR OTHER VECTOR FORMAT 
ggsave("temp_figure.pdf", plot = final_fig,
       width = 12, height = 12, units = "in", device = cairo_pdf)
## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 51 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'

The figure above highlights the variability across chambers within each pile while also revealing the overall pile-level trends.A parallel figure for the ratio data would follow the same logic. We could overlay zone thresholds — for example, ranges suggestive of aerobic vs. anaerobic conditions, or another biologically meaningful cutoff — though I’m open to alternatives depending on what best communicates the data. Whether or not we include color-coded thresholds, the layout would mirror the temperature figure: a 2 × 4 grid, with each column representing a pile and each row showing a chamber, with raw observations overlaid on DOY averages. The bottom panel of each column would show only the DOY averages with a fitted trendline, summarizing pile-level behavior.

Below would be the format if you want to show the area under the curve - as we are looking at ratio the area under the curve does not matter. Still this format could be used to show how the piles differ ? Personally I like the other one above but this is your section and you can choose how to present following the theme of this one of the other or if you have a better one go for it. This is just an example I would save it using ggsave to add it to the paper.

turns_both <- c(159, 186, 208, 234, 262, 290)
turns_control_only <- c(173, 191, 222, 249, 276)


# Layer 1: All individual observations (use raw DOY for precise timing)
co2_all_obs <- co2_clean %>%
  filter(!is.na(FCO2_DRY.LIN)) %>%
  mutate(flux_g = FCO2_DRY.LIN * conv_co2)

# Layer 2: Daily treatment means
co2_plot_means <- co2_treatment_daily %>%
  mutate(flux_g = daily_treatment_mean * conv_co2)

ggplot() +
  geom_vline(xintercept = turns_both, 
             linetype = "solid", colour = "grey50", linewidth = 0.3) +
  geom_vline(xintercept = turns_control_only, 
             linetype = "dashed", colour = "grey50", linewidth = 0.3) +
  geom_area(data = co2_plot_means,
            aes(x = DOY_int, y = flux_g, fill = Pile),
            alpha = 0.15, position = "identity") +
  geom_point(data = co2_all_obs,
             aes(x = DOY, y = flux_g, colour = Pile),
             alpha = 0.15, size = 0.8) +
  geom_line(data = co2_plot_means,
            aes(x = DOY_int, y = flux_g, colour = Pile),
            linewidth = 0.5) +
  geom_point(data = co2_plot_means,
             aes(x = DOY_int, y = flux_g, colour = Pile),
             shape = 17, size = 2) +
  scale_colour_manual(values = c("C" = "#4477AA", "E" = "#EE7733"),
                      labels = c("Control (2×/month)", "Experimental (1×/month)")) +
  scale_fill_manual(values = c("C" = "#4477AA", "E" = "#EE7733"),
                    labels = c("Control (2×/month)", "Experimental (1×/month)")) +
  coord_cartesian(ylim = c(0, 3500)) +
  labs(
    x = "Day of Year",
    y = expression(CO[2]~Flux~(g~m^{-2}~d^{-1})),
    colour = "Treatment",
    fill = "Treatment"
  ) +
  theme_classic() +
  theme(
    legend.position = "bottom",
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    legend.text = element_text(size = 10)
  )