summaries <- map_dfc(
  list(TotalN = SoilDF$TotalN, TotalC = SoilDF$TotalC, 
       Nitrate = SoilDF$Nitrate, AmmoniumN = SoilDF$AmmoniumN, P = SoilDF$P),
  summary
)

pivot_table <- SoilDF %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  select(Month, Severity, TotalN, TotalC, Nitrate, Nitratelbs, P, AmmoniumN, Ammoniumlbs) %>%
  group_by(Severity, Month) %>%
  summarise(across(
    c(TotalN, TotalC, Nitrate, Nitratelbs, P, AmmoniumN, Ammoniumlbs),
    list(Avg = ~ mean(.x, na.rm = TRUE)),
    .names = "Avg_{col}"
  ), .groups = "drop")

kable(pivot_table)
Severity Month Avg_TotalN Avg_TotalC Avg_Nitrate Avg_Nitratelbs Avg_P Avg_AmmoniumN Avg_Ammoniumlbs
High May 1391.0606 1.788879 0.9939394 2.303030 6.272727 4.607273 11.090909
High July 1482.7941 2.136029 2.0941176 4.941177 7.985294 4.943529 11.941176
Low May 839.1818 1.958182 0.8727273 2.090909 3.181818 3.764546 9.000000
Low July 930.4545 1.976909 1.3818182 3.363636 4.100000 2.157273 5.181818
Medium May 1223.8214 1.811893 1.1500000 2.785714 4.114286 4.272143 10.214286
Medium July 1327.3333 1.990704 0.8555556 2.148148 4.907407 3.057778 7.296296
Unburned May 1488.9412 2.257588 0.5529412 1.352941 3.864706 4.165882 10.058824
Unburned July 1622.8235 2.295118 0.7882353 1.764706 3.117647 3.724118 9.058824
Very high May 1836.1818 2.307909 0.7909091 1.818182 7.872727 4.927273 11.818182
Very high July 1566.2727 2.273364 1.8363636 4.363636 5.472727 3.830909 9.272727 
ggplot(SoilDF, aes(x = TotalN)) +
  geom_histogram(binwidth = 100, fill = "skyblue", color = "black", alpha = 1) + facet_wrap(~Month)+
  labs(title = "Total Nitrogen",
       x = "ppm",
       y= "Frequency") +
  theme_minimal()

ggplot(SoilDF, aes(x = TotalC)) +
  geom_histogram(binwidth = .2, fill = "skyblue", color = "black", alpha = 1) +
  facet_wrap(~Month)+
  labs(title = "Total Carbon",
       x = "Percentage %",
       y= "Frequency") +
  theme_minimal()

ggplot(SoilDF, aes(x = Nitrate)) +
  geom_histogram(binwidth = .5, fill = "skyblue", color = "black", alpha = 1) + facet_wrap(~Month)+
  labs(title = "Total Nitrate",
       x = "ppm",
       y= "Frequency") +
  theme_minimal()

ggplot(SoilDF, aes(x = P)) +
  geom_histogram(binwidth = 1, fill = "skyblue", color = "black", alpha = 1) + facet_wrap(~Month)+
  labs(title = "Total Phosporus",
       x = "ppm",
       y= "Frequency") +
  theme_minimal()

ggplot(SoilDF, aes(x = AmmoniumN)) +
  geom_histogram(binwidth = 1, fill = "skyblue", color = "black", alpha = 1) + facet_wrap(~Month)+
  labs(title = "Total AmmoniumN",
       x = "ppm",
       y= "Frequency") +
  theme_minimal()

SoilDF_long <- SoilDF %>%
  pivot_longer(cols = c(TotalN), names_to = "Variable", values_to = "Value")


ggplot(SoilDF_long, aes(x = Value, y = Severity, fill = Variable)) +
  geom_boxplot() +
  labs(
    title = "Total Nitrogen",
    x = "Value",
    y = "Severity",
    fill = "Nutrient"
  ) +
  theme_minimal()

# Reshape the data for combined box plots
SoilDF_long <- SoilDF %>%
  pivot_longer(
    cols = c(Nitrate, AmmoniumN, P),
    names_to = "Variable",
    values_to = "Value"
  )

# Create a combined box plot
ggplot(SoilDF_long, aes(x = Value, y = Severity, fill = Variable)) +
  geom_boxplot() +
  labs(
    title = "Combined Box Plot for Nitrate, AmmoniumN, and P",
    x = "PPM",
    y = "Severity",
    fill = "Nutrient"
  ) +
  theme_minimal()

Nut_grid <- SoilDATA %>%
  mutate(
    Month = factor(Month, levels = month.name),
    Severity = factor(Severity, levels = c("Unburned", "Low", "Med", "High", "Very High"))
  ) %>%
  pivot_longer(
    cols = c("Nitrate", "Nitratelbs", "P", "AmmoniumN", "Ammoniumlbs", "TotalN", "TotalC"),
    names_to = "Nutrient",
    values_to = "Nutr_val",
    values_drop_na = FALSE
  ) %>%
  group_by(Month, Severity, Nutrient) %>%
  summarise(
    nut_mn = mean(Nutr_val, na.rm = TRUE),
    nut_sd = sd(Nutr_val, na.rm = TRUE),
    count = n(),
    nut_se = nut_sd / sqrt(count),
    .groups = "drop"
  ) %>%
  filter(Nutrient %in% c("AmmoniumN", "Nitrate", "P", "TotalC", "TotalN")) %>%
  ggplot(aes(x = Severity, y = nut_mn, color = Month)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin = nut_mn - nut_se, ymax = nut_mn + nut_se), width = 0.2, position = position_dodge(0.05)) +
  labs(
    title = "Nutrient vs Fire",
    x = "Severity",
    y = "Mean ± SE"
  ) +
  theme_bw() +
  facet_wrap(~ Nutrient, ncol = 3, scales = "free")

print(Nut_grid)

# Forage Quality Distributions

# Factorize the Month variable and create the histogram in a single pipeline
ForageAmended %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  ggplot(aes(x = Moisture)) +
  geom_histogram(binwidth = 0.1, fill = "skyblue", color = "black", alpha = 1) +
  facet_wrap(~ Month) +
  labs(
    title = "Moisture Content",
    x = "Percentage As Received",
    y = "Frequency"
  ) +
  theme_minimal()

Prot_grid <- ForageAmended %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  pivot_longer(
    cols = c("ProteinDW"),
    names_to = "Crude",
    values_to = "Prot_val",
    values_drop_na = TRUE
  ) %>%
  group_by(Month, EcoSite, Fire, Crude) %>%
  summarise(
    Prot_mn = mean(Prot_val, na.rm = TRUE),
    Prot_sd = sd(Prot_val, na.rm = TRUE),
    count = n(),
    Prot_se = Prot_sd / sqrt(count),
    .groups = "drop"
  ) %>%
  filter(Crude == "ProteinDW") %>%
  ggplot(aes(x = Month, y = Prot_mn, group = EcoSite, color = Fire)) +
  theme_light() +
  geom_point() +
  geom_errorbar(aes(ymin = Prot_mn - Prot_se, ymax = Prot_mn + Prot_se), width = 0.2, position = position_dodge(0.05)) +
  labs(
    title = "ProteinDW",
    x = "Month",
    y = "Mean ± SE"
  ) +
  scale_color_manual(values = c('red', 'blue', 'orange', 'green')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Add angled month labels
  facet_wrap(~EcoSite, ncol = 2, scales = "free")

print(Prot_grid)

When Ecosite is removed

Process data and generate ProteinDW grid plot

Prot_grid <- ForageAmended %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  pivot_longer(
    cols = c("ProteinDW"),
    names_to = "Crude",
    values_to = "Prot_val",
    values_drop_na = TRUE
  ) %>%
  group_by(Month, Fire, Crude) %>%
  summarise(
    Prot_mn = mean(Prot_val, na.rm = TRUE),
    Prot_sd = sd(Prot_val, na.rm = TRUE),
    count = n(),
    Prot_se = Prot_sd / sqrt(count),
    .groups = "drop"
  ) %>%
  filter(Crude == "ProteinDW") %>%
  ggplot(aes(x = Month, y = Prot_mn, group = Fire, color = Fire)) +
  theme_light() +
  geom_point() +
  geom_errorbar(aes(ymin = Prot_mn - Prot_se, ymax = Prot_mn + Prot_se), width = 0.2, position = position_dodge(0.05)) +
  labs(
    title = "ProteinDW",
    x = "Month",
    y = "Mean ± SE"
  ) +
  scale_color_manual(values = c('red', 'blue', 'orange', 'green')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Angled month labels
  facet_wrap(~ Fire, ncol = 3, scales = "free")

print(Prot_grid)

# Process data and generate Dry Matter grid plot with tilted month labels
DM_grid <- ForageAmended %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  pivot_longer(
    cols = c("DM"),
    names_to = "Dry",
    values_to = "Dry_val",
    values_drop_na = TRUE
  ) %>%
  group_by(Month, EcoSite, Fire, Dry) %>%
  summarise(
    DM_mn = mean(Dry_val, na.rm = TRUE),
    DM_sd = sd(Dry_val, na.rm = TRUE),
    count = n(),
    DM_se = DM_sd / sqrt(count),
    .groups = "drop"
  ) %>%
  filter(Dry == "DM") %>%
  ggplot(aes(x = Month, y = DM_mn, group = EcoSite, color = Fire)) +
  theme_light() +
  geom_point() +
  geom_errorbar(aes(ymin = DM_mn - DM_se, ymax = DM_mn + DM_se), width = 0.2, position = position_dodge(0.05)) +
  labs(
    title = "Dry Matter",
    x = "Month",
    y = "Mean ± SE"
  ) +
  scale_color_manual(values = c('red', 'blue', 'orange', 'green')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Tilt month labels
  facet_wrap(~ EcoSite, ncol = 2, scales = "free")

print(DM_grid)

break

# Process data and generate TDNDW grid plot
TD_grid <- ForageAmended %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  pivot_longer(
    cols = c("TDNDW"),
    names_to = "Energy",
    values_to = "TD_val",
    values_drop_na = TRUE
  ) %>%
  group_by(Month, EcoSite, Fire, Energy) %>%
  summarise(
    TD_mn = mean(TD_val, na.rm = TRUE),
    TD_sd = sd(TD_val, na.rm = TRUE),
    count = n(),
    TD_se = TD_sd / sqrt(count),
    .groups = "drop"
  ) %>%
  filter(Energy == "TDNDW") %>%
  ggplot(aes(x = Month, y = TD_mn, group = EcoSite, color = Fire)) +
  theme_light() +
  geom_point() +
  geom_errorbar(aes(ymin = TD_mn - TD_se, ymax = TD_mn + TD_se), width = 0.2, position = position_dodge(0.05)) +
  labs(
    title = "TDNDW(%)",
    x = "Month",
    y = "Mean ± SE"
  ) +
  scale_color_manual(values = c('red', 'blue', 'orange', 'green')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Tilt month labels
  facet_wrap(~ EcoSite, ncol = 3, scales = "free")

print(TD_grid)

# Process data and generate LactationDW grid plot
Lac_grid <- ForageAmended %>%
  mutate(Month = factor(Month, levels = month.name)) %>%
  pivot_longer(
    cols = c("LactationDW"),
    names_to = "Lactate",
    values_to = "Lac_val",
    values_drop_na = TRUE
  ) %>%
  group_by(Month, EcoSite, Fire, Lactate) %>%
  summarise(
    Lac_mn = mean(Lac_val, na.rm = TRUE),
    Lac_sd = sd(Lac_val, na.rm = TRUE),
    count = n(),
    Lac_se = Lac_sd / sqrt(count),
    .groups = "drop"
  ) %>%
  filter(Lactate == "LactationDW") %>%
  ggplot(aes(x = Month, y = Lac_mn, group = EcoSite, color = Fire)) +
  theme_light() +
  geom_point() +
  geom_errorbar(aes(ymin = Lac_mn - Lac_se, ymax = Lac_mn + Lac_se), width = 0.2, position = position_dodge(0.05)) +
  labs(
    title = "LactationDW(Mcal/lbs)",
    x = "Month",
    y = "Mean ± SE"
  ) +
  scale_color_manual(values = c('red', 'blue', 'orange', 'green')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Tilt month labels
  facet_wrap(~ EcoSite, ncol = 3, scales = "free")

print(Lac_grid)

ManLong <- ForageAmended %>%
  pivot_longer(
    cols = c("MaintDW"),
    names_to = "Maintenance",
    values_to = "Man-val",
    values_drop_na = TRUE
  )

Man_mean <- ManLong %>%
  group_by(Month, EcoSite, Fire, Maintenance ) %>%
  summarise(
    Man_mn = mean(`Man-val`, na.rm = TRUE),
    Man_sd = sd(`Man-val`, na.rm = TRUE),
    count = n(),
    Man_se = Man_sd / sqrt(count)
  ) %>% 
  ungroup()


Man_dat <- Man_mean %>% filter(Maintenance=="MaintDW")
Man_dat$Month <- factor(Man_dat$Month, levels = month.name)


Man<- ggplot(Man_dat, aes(x=Month, y=Man_mn, group=EcoSite, color=Fire)) + 
  theme_light() +
  geom_point()+
  geom_errorbar(aes(ymin=Man_mn-Man_se, ymax=Man_mn+Man_se), width=.2,
                position=position_dodge(0.05)) +
  labs(title="MaintenanceDW(Mcal/lbs)", x="Month", y = "MeanSE+-")+
  
  scale_color_manual(values=c('red','blue','orange','green'))

Man_grid <- Man + facet_wrap(~EcoSite, ncol = 3, scales = "free")



print(Man_grid)

GALong <- ForageAmended %>%
  pivot_longer(
    cols = c("GainDW"),
    names_to = "Gain",
    values_to = "GA-val",
    values_drop_na = TRUE
  )




GA_mean <- GALong %>%
  group_by(Month, EcoSite, Fire, Gain ) %>%
  summarise(
    GA_mn = mean(`GA-val`, na.rm = TRUE),
    GA_sd = sd(`GA-val`, na.rm = TRUE),
    count = n(),
    GA_se = GA_sd / sqrt(count)
  ) %>% 
  ungroup()


GA_dat <- GA_mean %>% filter(Gain=="GainDW")

GA_dat$Month <- factor(GA_dat$Month, levels = month.name)


GA<- ggplot(GA_dat, aes(x=Month, y=GA_mn, group=EcoSite, color=Fire)) + 
  theme_light() +
  geom_point()+
  geom_errorbar(aes(ymin=GA_mn-GA_se, ymax=GA_mn+GA_se), width=.2,
                position=position_dodge(0.05)) +
  labs(title="GainDW(Mcal/lbs)", x="Month", y = "MeanSE+-")+
  
  scale_color_manual(values=c('red','blue','orange','green'))

GA_grid <- GA + facet_wrap(~EcoSite, ncol = 3, scales = "free")



print(GA_grid)