GreenFeed Final Data Report
In greenfeedr: Process and Report 'GreenFeed' Data

library(dplyr)
library(ggplot2)

Study: `r exp`

GreenFeed unit(s): `r unit`

Experimental days: `r floor(as.numeric(difftime(max(df$StartTime), min(as.Date(df$StartTime)), units = "days") +1))` days

Number of animals with gas records: r if (!is.null(rfid_file) && is.data.frame(rfid_file) && nrow(rfid_file) > 0) { paste0(nrow(dplyr::semi_join(rfid_file, df, by = "RFID")), " out of ", nrow(rfid_file)) } else { length(unique(df$RFID)) }

\newpage

Gas records per Day

cols_to_convert <- c("CH4GramsPerDay", "CO2GramsPerDay", "O2GramsPerDay", "H2GramsPerDay")
df[cols_to_convert] <- lapply(df[cols_to_convert], as.numeric)

# Plot 1: Total number of production records per day
plot1 <- ggplot(as.data.frame(table(as.Date(df$StartTime))), aes(x = Var1, y = Freq)) +
  geom_col(color = "black") +
  labs(
    title = "Total Records Per Day",
    x = "",
    y = "Total Records"
  ) +
  geom_text(aes(label = Freq), vjust = -0.5, color = "black", size = 2.2, position = position_dodge(width = 0.9)) +
  theme_classic() +
  theme(
    plot.title = element_text(size = 11, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1.05, size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.position = "none"
  )

plot1


generate_combined_plot <- function(df, plot_opt) {
  # Convert to lowercase to avoid case sensitivity issues
  plot_opt <- tolower(plot_opt)

  if ("all" %in% plot_opt) {
    options_selected <- c("ch4", "o2", "co2", "h2")
  } else {
    options_selected <- plot_opt
  }

  # Normalize the data
  df_normalized <- df %>%
    dplyr::mutate(
      Normalized_CH4 = scale(CH4GramsPerDay),
      Normalized_CO2 = scale(CO2GramsPerDay),
      Normalized_O2 = scale(O2GramsPerDay),
      Normalized_H2 = scale(H2GramsPerDay)
    )

  # Create a base plot
  combined_plot <- ggplot(df_normalized[df_normalized$HourOfDay <= 23, ], aes(x = HourOfDay)) +
    theme_classic() +
    theme(
      plot.title = element_text(size = 11, face = "bold"),
      axis.text.x = element_text(angle = 0, size = 9),
      axis.text.y = element_text(angle = 0, size = 9),
      axis.title.y = element_text(size = 10, face = "bold"),
      axis.title.x = element_text(size = 10, face = "bold"),
      legend.position = "bottom"
    ) +
    labs(
      title = "Gas Production Across The Day",
      x = "",
      y = "Normalized Gas Value",
      color = "Gas type"
    ) +
    scale_x_continuous(
      breaks = seq(0, 23),
      labels = c(
        "12 AM", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",
        "12 PM", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11"
      )
    )

  # Initialize empty lists for points and smooth lines
  points_list <- list()
  smooth_list <- list()

  # Loop through selected options to add points and smooth lines
  for (gas in options_selected) {
    if (gas == "ch4") {
      points_list <- c(points_list, list(geom_point(aes(y = Normalized_CH4), color = "#E0E0E0")))
      smooth_list <- c(smooth_list, list(geom_smooth(aes(y = Normalized_CH4, color = "CH4"), se = FALSE)))
    }
    if (gas == "co2") {
      points_list <- c(points_list, list(geom_point(aes(y = Normalized_CO2), color = "#E0E0E0")))
      smooth_list <- c(smooth_list, list(geom_smooth(aes(y = Normalized_CO2, color = "CO2"), se = FALSE)))
    }
    if (gas == "o2") {
      points_list <- c(points_list, list(geom_point(aes(y = Normalized_O2), color = "#E0E0E0")))
      smooth_list <- c(smooth_list, list(geom_smooth(aes(y = Normalized_O2, color = "O2"), se = FALSE)))
    }
    if (gas == "h2") {
      points_list <- c(points_list, list(geom_point(aes(y = Normalized_H2), color = "#E0E0E0")))
      smooth_list <- c(smooth_list, list(geom_smooth(aes(y = Normalized_H2, color = "H2"), se = FALSE)))
    }
  }

  # Add points and smooth lines to the combined plot
  combined_plot <- combined_plot + do.call("list", points_list)
  combined_plot <- combined_plot + do.call("list", smooth_list)

  # Check for valid options
  valid_options <- c("all", "ch4", "co2", "o2", "h2")
  if (!all(options_selected %in% valid_options)) {
    stop("Invalid plot option selected.")
  }

  return(combined_plot)
}

plot_combined <- generate_combined_plot(df, plot_opt)
message(plot_combined)

\newpage

Gas records per animal

# Assuming RFIDfile is provided or not, set the grouping variable
group_var <- if (!is.null(rfid_file) && is.data.frame(rfid_file) && nrow(rfid_file) > 0) "FarmName" else "RFID"

# Determine the order of FarmName or RFID based on daily_CH4
farmname_order <- df %>%
  dplyr::mutate(day = as.Date(EndTime)) %>%
  dplyr::group_by(!!sym(group_var), day) %>%
  dplyr::summarise(
    n = n(),
    daily_CH4 = weighted.mean(CH4GramsPerDay, GoodDataDuration, na.rm = TRUE)
  ) %>%
  dplyr::group_by(!!sym(group_var)) %>%
  dplyr::summarise(
    n = sum(n),
    daily_CH4 = mean(daily_CH4, na.rm = TRUE)
  ) %>%
  dplyr::arrange(desc(daily_CH4)) %>%
  dplyr::pull(!!sym(group_var))


# Plot 1: Total number of records per animal
plot1 <- df %>%
  dplyr::mutate(day = as.Date(EndTime)) %>%
  dplyr::group_by(!!sym(group_var), day) %>%
  dplyr::summarise(
    n = n(),
    daily_CH4 = weighted.mean(CH4GramsPerDay, GoodDataDuration, na.rm = TRUE)
  ) %>%
  dplyr::group_by(!!sym(group_var)) %>%
  dplyr::summarise(
    n = sum(n),
    daily_CH4 = mean(daily_CH4, na.rm = TRUE)
  ) %>%
  ggplot(aes(x = factor(!!sym(group_var), levels = farmname_order), y = n)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  labs(
    title = "Total Records Per Animal",
    x = "",
    y = "Total Records"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(size = 11, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1.05, size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.position = "none"
  ) +
  geom_text(aes(label = n), vjust = -1, color = "black", position = position_dodge(width = 0.9), size = 2.2)


# Plot distribution of records throughout the day
plot2 <- df %>%
  dplyr::mutate(AMPM = case_when(
    HourOfDay >= 22 ~ "10PM-4AM",
    HourOfDay < 4 ~ "10PM-4AM",
    HourOfDay >= 4 & HourOfDay < 10 ~ "4AM-10AM",
    HourOfDay >= 10 & HourOfDay < 16 ~ "10AM-4PM",
    HourOfDay >= 16 & HourOfDay < 22 ~ "4PM-10PM",
    TRUE ~ NA_character_
  )) %>%
  dplyr::group_by(!!sym(group_var), AMPM) %>%
  dplyr::summarise(n = n()) %>%
  ggplot(aes(
    x = factor(!!sym(group_var), levels = farmname_order), y = n,
    fill = factor(AMPM, levels = c("10PM-4AM", "4AM-10AM", "10AM-4PM", "4PM-10PM"))
  )) +
  geom_bar(stat = "identity", position = "fill") +
  labs(
    title = "Daily Records Distribution",
    x = "",
    y = "Percentage of total records",
    fill = "Time-Windows (24h)"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(size = 11, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1.05, size = 8),
    legend.position = "bottom",
    axis.title.y = element_text(size = 10, face = "bold")
  ) +
  scale_fill_brewer(palette = "BrBG") +
  scale_y_continuous(breaks = c(0, 0.25, 0.50, 0.75, 1), labels = c("0%", "25%", "50%", "75%", "100%"), expand = c(0, 0))


# Create the plots
plot1
plot2

\newpage

Daily gas production per Animal

## The daily averages for different gases is calculated as the weighted mean using visit time to the GreenFeed (or 'GoodDataDuration')

# Function to generate plots
generate_plots <- function(data, group_var, plot_opt = c("All", "CH4", "O2", "CO2", "H2")) {
  # Convert to lowercase to avoid case sensitivity issues
  plot_opt <- tolower(plot_opt)

  if ("all" %in% plot_opt) {
    options_selected <- c("ch4", "o2", "co2", "h2")
  } else {
    options_selected <- plot_opt
  }

  plots <- list()

  if ("ch4" %in% options_selected) {
    p1 <- df %>%
      dplyr::mutate(day = as.Date(EndTime)) %>%
      dplyr::group_by(!!sym(group_var), day) %>%
      dplyr::summarise(daily_CH4 = weighted.mean(CH4GramsPerDay, GoodDataDuration, na.rm = TRUE)) %>%
      {
        ggplot(., aes(x = reorder(!!sym(group_var), -daily_CH4), y = daily_CH4, color = daily_CH4)) +
          geom_boxplot(fatten = NULL, outlier.shape = NA) +
          stat_summary(
            fun = mean, geom = "errorbar",
            aes(ymax = ..y.., ymin = ..y..), width = 0.75, size = 0.7,
            color = "black", linetype = "solid"
          ) +
          labs(
            title = "Methane (CH4) Production Per Animal",
            x = "",
            y = "CH4 (g/d)"
          ) +
          theme_classic() +
          theme(
            plot.title = element_text(size = 11, face = "bold"),
            axis.text.x = element_text(angle = 45, hjust = 1.05, size = 5),
            axis.title.y = element_text(size = 10, face = "bold"),
            legend.position = "none"
          ) +
          geom_hline(yintercept = mean(.$daily_CH4), linetype = "dashed", color = "blue", linewidth = 0.6) +
          scale_y_continuous(breaks = seq(0, max(.$daily_CH4), by = 100))
      }
    plots <- c(plots, list(p1))
  }

  if ("co2" %in% options_selected) {
    p3 <- df %>%
      dplyr::mutate(day = as.Date(EndTime)) %>%
      dplyr::group_by(!!sym(group_var), day) %>%
      dplyr::summarise(daily_CO2 = weighted.mean(CO2GramsPerDay, GoodDataDuration, na.rm = TRUE)) %>%
      {
        ggplot(., aes(x = reorder(!!sym(group_var), -daily_CO2), y = daily_CO2, color = daily_CO2)) +
          geom_boxplot(fatten = NULL, outlier.shape = NA) +
          stat_summary(
            fun = mean, geom = "errorbar",
            aes(ymax = ..y.., ymin = ..y..), width = 0.75, size = 0.7,
            color = "black", linetype = "solid"
          ) +
          labs(
            title = "Carbon Dioxide (CO2) Production Per Animal",
            x = "",
            y = "CO2 (g/d)"
          ) +
          theme_classic() +
          theme(
            plot.title = element_text(size = 11, face = "bold"),
            axis.text.x = element_text(angle = 45, hjust = 1.05, size = 5),
            axis.title.y = element_text(size = 10, face = "bold"),
            legend.position = "none"
          ) +
          geom_hline(yintercept = mean(.$daily_CO2), linetype = "dashed", color = "red", linewidth = 0.6) +
          scale_y_continuous(breaks = seq(0, max(.$daily_CO2), by = 2000))
      }
    plots <- c(plots, list(p3))
  }

  if ("o2" %in% options_selected) {
    p2 <- df %>%
      dplyr::mutate(day = as.Date(EndTime)) %>%
      dplyr::group_by(!!sym(group_var), day) %>%
      dplyr::summarise(daily_O2 = weighted.mean(O2GramsPerDay, GoodDataDuration, na.rm = TRUE)) %>%
      {
        ggplot(., aes(x = reorder(!!sym(group_var), -daily_O2), y = daily_O2, color = daily_O2)) +
          geom_boxplot(fatten = NULL, outlier.shape = NA) +
          stat_summary(
            fun = mean, geom = "errorbar",
            aes(ymax = ..y.., ymin = ..y..), width = 0.75, size = 0.7,
            color = "black", linetype = "solid"
          ) +
          labs(
            title = "Oxygen (O2) Production Per Animal",
            x = "",
            y = "O2 (g/d)"
          ) +
          theme_classic() +
          theme(
            plot.title = element_text(size = 11, face = "bold"),
            axis.text.x = element_text(angle = 45, hjust = 1.05, size = 5),
            axis.title.y = element_text(size = 10, face = "bold"),
            legend.position = "none"
          ) +
          geom_hline(yintercept = mean(.$daily_O2), linetype = "dashed", color = "orange", linewidth = 0.6) +
          scale_y_continuous(breaks = seq(0, max(.$daily_O2), by = 2000))
      }
    plots <- c(plots, list(p2))
  }

  if ("h2" %in% options_selected) {
    if (any(df$H2GramsPerDay > 0)) {
      p4 <- df %>%
        dplyr::mutate(day = as.Date(EndTime)) %>%
        dplyr::group_by(!!sym(group_var), day) %>%
        dplyr::summarise(daily_H2 = weighted.mean(H2GramsPerDay, GoodDataDuration, na.rm = TRUE)) %>%
        {
          ggplot(., aes(x = reorder(!!sym(group_var), -daily_H2), y = daily_H2, color = daily_H2)) +
            geom_boxplot(fatten = NULL, outlier.shape = NA) +
            stat_summary(
              fun = mean, geom = "errorbar",
              aes(ymax = ..y.., ymin = ..y..), width = 0.75, size = 0.7,
              color = "black", linetype = "solid"
            ) +
            labs(
              title = "Hydrogen (H2) Production Per Animal",
              x = "",
              y = "H2 (g/d)"
            ) +
            theme_classic() +
            theme(
              plot.title = element_text(size = 11, face = "bold"),
              axis.text.x = element_text(angle = 45, hjust = 1.05, size = 5),
              axis.title.y = element_text(size = 10, face = "bold"),
              legend.position = "none"
            ) +
            geom_hline(yintercept = mean(.$daily_H2), linetype = "dashed", color = "purple", linewidth = 0.6) +
            scale_y_continuous(breaks = seq(0, max(.$daily_H2), by = 0.5))
        }
      plots <- c(plots, list(p4))
    } else {
      message("No hydrogen (H2) data from GreenFeed")
    }
  }

  return(plots)
}

# Call the function and display the plots
plots <- generate_plots(df, group_var, plot_opt)
for (p in plots) {
  print(p)
}