Make this into a function:

harvest_stats = describeBy(harvest_data %>%
                             select(Year, DEC_BASIS, PEO_BASIS, ESL_BASIS, FUT_1, FUT_2, Crop_Progression), harvest_data$Year, digits=2, mat=TRUE)


harvest_means <- harvest_stats %>%
  filter(vars %in% c(2, 3, 4, 5, 6, 7)) %>%  # Adjust the 'vars' numbers to match your actual variables
  select(item, group1, mean, min, max) %>% 
  mutate(Year = as.numeric(group1)) %>% 
  select(-group1) %>% 
  mutate(
    variable = gsub("[0-9]+", "", rownames(.)),  # Remove numeric suffixes
    variable = case_when(
      grepl("^FUT_1", rownames(.)) ~ "FUT_1",
      grepl("^FUT_2", rownames(.)) ~ "FUT_2",
      TRUE ~ variable  # Retain the previously assigned value for other variables
    )
  )

# Clean the variable names in harvest_means to match those in Harvest_min_max
harvest_means_clean <- harvest_means %>%
  mutate(variable = case_when(
    grepl("^DEC_BASIS", variable) ~ "DEC_BASIS",
    grepl("^PEO_BASIS", variable) ~ "PEO_BASIS",
    grepl("^ESL_BASIS", variable) ~ "ESL_BASIS",
    grepl("^FUT_1", variable) ~ "FUT_1",
    grepl("^FUT_2", variable) ~ "FUT_2",
    grepl("^Crop_Progression", variable) ~ "Crop_Progression",
    TRUE ~ variable
  ))

# Join harvest_means with Harvest_min_max based on Year and variable
harvest_means_with_dates <- harvest_means_clean %>%
  left_join(Harvest_min_max, by = "Year") %>%
  mutate(min_date = case_when(
    variable == "DEC_BASIS" ~ DEC_BASIS_date_min,
    variable == "PEO_BASIS" ~ PEO_BASIS_date_min,
    variable == "ESL_BASIS" ~ ESL_BASIS_date_min,
    variable == "FUT_1" ~ FUT_1_date_min,
    variable == "FUT_2" ~ FUT_2_date_min,
    variable == "Crop_Progression" ~ Crop_Progression_date_min,
    TRUE ~ NA_POSIXct_
  ),
  max_date = case_when(
    variable == "DEC_BASIS" ~ DEC_BASIS_date_max,
    variable == "PEO_BASIS" ~ PEO_BASIS_date_max,
    variable == "ESL_BASIS" ~ ESL_BASIS_date_max,
    variable == "FUT_1" ~ FUT_1_date_max,
    variable == "FUT_2" ~ FUT_2_date_max,
    variable == "Crop_Progression" ~ Crop_Progression_date_max,
    TRUE ~ NA_POSIXct_
  )) %>%
  select(Year, variable, mean, min, min_date, max, max_date)

# Print the result
print(harvest_means_with_dates)

# -------------------------------------------------------------

non_harvest_stats = describeBy(non_harvest_data %>%
                                 select(Year, DEC_BASIS, PEO_BASIS, ESL_BASIS, FUT_1, FUT_2, Crop_Progression), non_harvest_data$Year, digits=2, mat=TRUE)


non_harvest_means <- non_harvest_stats %>%
  filter(vars %in% c(2, 3, 4, 5, 6, 7)) %>%  # Adjust the 'vars' numbers to match your actual variables
  select(item, group1, mean, min, max) %>% 
  mutate(Year = as.numeric(group1)) %>% 
  select(-group1) %>% 
  mutate(
    variable = gsub("[0-9]+", "", rownames(.)),  # Remove numeric suffixes
    variable = case_when(
      grepl("^FUT_1", rownames(.)) ~ "FUT_1",
      grepl("^FUT_2", rownames(.)) ~ "FUT_2",
      TRUE ~ variable  # Retain the previously assigned value for other variables
    )
  )

# Clean the variable names in non_harvest_means to match those in Harvest_min_max
non_harvest_means_clean <- non_harvest_means %>%
  mutate(variable = case_when(
    grepl("^DEC_BASIS", variable) ~ "DEC_BASIS",
    grepl("^PEO_BASIS", variable) ~ "PEO_BASIS",
    grepl("^ESL_BASIS", variable) ~ "ESL_BASIS",
    grepl("^FUT_1", variable) ~ "FUT_1",
    grepl("^FUT_2", variable) ~ "FUT_2",
    grepl("^Crop_Progression", variable) ~ "Crop_Progression",
    TRUE ~ variable
  ))

# Join non_harvest_means with Harvest_min_max based on Year and variable
non_harvest_means_with_dates <- non_harvest_means_clean %>%
  left_join(non_Harvest_min_max, by = "Year") %>%
  mutate(min_date = case_when(
    variable == "DEC_BASIS" ~ DEC_BASIS_date_min,
    variable == "PEO_BASIS" ~ PEO_BASIS_date_min,
    variable == "ESL_BASIS" ~ ESL_BASIS_date_min,
    variable == "FUT_1" ~ FUT_1_date_min,
    variable == "FUT_2" ~ FUT_2_date_min,
    variable == "Crop_Progression" ~ Crop_Progression_date_min,
    TRUE ~ NA_POSIXct_
  ),
  max_date = case_when(
    variable == "DEC_BASIS" ~ DEC_BASIS_date_max,
    variable == "PEO_BASIS" ~ PEO_BASIS_date_max,
    variable == "ESL_BASIS" ~ ESL_BASIS_date_max,
    variable == "FUT_1" ~ FUT_1_date_max,
    variable == "FUT_2" ~ FUT_2_date_max,
    variable == "Crop_Progression" ~ Crop_Progression_date_max,
    TRUE ~ NA_POSIXct_
  )) %>%
  select(Year, variable, mean, min, min_date, max, max_date)

# Print the result
print(non_harvest_means_with_dates)