Chap1_Obj2

pacman::p_load(dplyr, readxl, tidyr, raster, vegan, tigris, sf, sp, plotly, ggrepel, kableExtra, brms)

# Tree PCQ Data
tree_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Tree_PCQ")
tree_data <- tree_data %>%
  filter(Authority %in% c("BRSF", "WSF", "Jay", "BCNWR", "Escambia", "Hooper", "Howell", "ONF"))

# Soil Data
fuel_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
                        sheet = "Fuel_Sampling")
fuel_data <- fuel_data %>%
  filter(Authority %in% c("BRSF", "WSF", "Jay", "BCNWR", "Escambia", "Hooper", "Howell", "ONF"))

# Site Data
CogonSites <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/CogonSites_FL_AL_MS.xlsx")
CogonSites <- CogonSites %>%
  filter(Authority %in% c("BRSF", "WSF", "Jay", "BCNWR", "Escambia", "Hooper", "Howell", "ONF"))

# Only include Florida/Alabama Sites
CogonSites <- CogonSites[CogonSites$Authority != "CNF" & CogonSites$Authority != "DSNF", ]

Merge

Filter All data to only include specified species (Per PLANTS database)

Filter all data to only include species found at 3% of all sites

#Fuel Dynamics ## Combine seasonal fuel data

Net Weight Method

Live

fuel_Live_Data_Net <- fuel_data_Net %>%
  mutate(
    Live_Bag = as.numeric(Live_Bag),
    Live_Weight_Post = as.numeric(Live_Weight_Post),
    Live_Weight_Initial = as.numeric(Live_Weight_Initial),
    Live_Height = as.numeric(Height),
    Net_Live = as.numeric(Net_Live),
    Status = as.character(Status)  # Status as a character
  )

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Live_Height = as.numeric(Height)`.
## Caused by warning:
## ! NAs introduced by coercion

fuel_Live_Data_Net <- fuel_Live_Data_Net %>%
  mutate(biomass = Net_Live)

fuel_Live_Data_Net <- fuel_Live_Data_Net %>%
  filter(biomass >= 0)

fuel_Live_Data_Net <- fuel_Live_Data_Net %>%
  mutate(relative_moisture_content = ifelse(biomass > bioT, ((Live_Weight_Initial - Live_Bag) - Net_Live) / Net_Live * 100, NA))

avg_live_values_Net <- fuel_Live_Data_Net %>%
  group_by(Plot, Status) %>%
  summarize(avg_live_biomass = mean(biomass, na.rm = TRUE),
            avg_live_moisture_content = mean(relative_moisture_content, na.rm = TRUE),
            avg_soil_moisture = mean(Soil_Moisture, na.rm = TRUE),
            avg_height = mean(Live_Height, na.rm = TRUE) / 100,
            .groups = "drop")

Dead

fuel_Dead_Data_Net <- fuel_data_Net %>%
  mutate(
    Dead_Bag = as.numeric(Dead_Bag),
    Dead_Weight_Post = as.numeric(Dead_Weight_Post),
    Dead_Weight_Initial = as.numeric(Dead_Weight_Initial),
    Net_Dead = as.numeric(Net_Dead),
    Status = as.character(Status)  # Status as a character
  )

fuel_Dead_Data_Net <- fuel_Dead_Data_Net %>%
  mutate(biomass = Net_Dead)

fuel_Dead_Data_Net <- fuel_Dead_Data_Net %>%
  filter(biomass >= 0)

fuel_Dead_Data_Net <- fuel_Dead_Data_Net %>%
  mutate(relative_moisture_content = ifelse(biomass > bioT, ((Dead_Weight_Initial - Dead_Bag) - Net_Dead) / Net_Dead * 100, NA))

avg_dead_values_Net <- fuel_Dead_Data_Net %>%
  group_by(Plot, Status) %>%
  summarize(avg_dead_biomass = mean(biomass, na.rm = TRUE),
            avg_dead_moisture_content = mean(relative_moisture_content, na.rm = TRUE),
            avg_soil_moisture = mean(Soil_Moisture, na.rm = TRUE),
            .groups = "drop")

Litter

fuel_Litter_Data_Net <- fuel_data_Net %>%
  mutate(
    Litter_Bag = as.numeric(Litter_Bag),
    Litter_Weight_Post = as.numeric(Litter_Weight_Post),
    Litter_Weight_Initial = as.numeric(Litter_Weight_Initial),
    Net_Litter = as.numeric(Net_Litter),
    Status = as.character(Status)  # Status as a character
  )

fuel_Litter_Data_Net <- fuel_Litter_Data_Net %>%
  mutate(biomass = Net_Litter)

fuel_Litter_Data_Net <- fuel_Litter_Data_Net %>%
  filter(biomass >= 0)

fuel_Litter_Data_Net <- fuel_Litter_Data_Net %>%
  mutate(relative_moisture_content = ifelse(biomass > bioT, ((Litter_Weight_Initial - Litter_Bag) - Net_Litter) / Net_Litter * 100, NA))

avg_litter_values_Net <- fuel_Litter_Data_Net %>%
  group_by(Plot, Season, Status) %>%
  summarize(avg_litter_biomass = mean(biomass, na.rm = TRUE),
            avg_litter_moisture_content = mean(relative_moisture_content, na.rm = TRUE),
            avg_soil_moisture = mean(Soil_Moisture, na.rm = TRUE),
            .groups = "drop")

Average Bag Weight Method

Live

fuel_Live_Data_Avg <- fuel_data_Avg %>%
  mutate(
    Live_Bag = as.numeric(Live_Bag),
    Live_Weight_Post = as.numeric(Live_Weight_Post),
    Live_Weight_Initial = as.numeric(Live_Weight_Initial),
    Live_Height = as.numeric(Height),
    Dry_LiveBag = as.numeric(Dry_LiveBag),
    Status = as.character(Status)  # Status as a character
  )

## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Live_Height = as.numeric(Height)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

fuel_Live_Data_Avg <- fuel_Live_Data_Avg %>%
  mutate(biomass = Live_Weight_Post - Dry_LiveBag)

fuel_Live_Data_Avg <- fuel_Live_Data_Avg %>%
  filter(biomass >= 0)

fuel_Live_Data_Avg <- fuel_Live_Data_Avg %>%
  mutate(relative_moisture_content = ifelse(biomass > bioT, (Live_Weight_Initial - Live_Weight_Post) / biomass * 100, NA))

avg_live_values_Avg <- fuel_Live_Data_Avg %>%
  group_by(Plot, Status) %>%
  summarize(avg_live_biomass = mean(biomass, na.rm = TRUE),
            avg_live_moisture_content = mean(relative_moisture_content, na.rm = TRUE),
            avg_soil_moisture = mean(Soil_Moisture, na.rm = TRUE),
            avg_height = mean(Live_Height, na.rm = TRUE) / 100,
            .groups = "drop")

Dead

fuel_Dead_Data_Avg <- fuel_data_Avg %>%
  mutate(
    Dead_Bag = as.numeric(Dead_Bag),
    Dead_Weight_Post = as.numeric(Dead_Weight_Post),
    Dead_Weight_Initial = as.numeric(Dead_Weight_Initial),
    Dry_DeadBag = as.numeric(Dry_DeadBag),
    Status = as.character(Status)  # Status as a character
  )

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Dry_DeadBag = as.numeric(Dry_DeadBag)`.
## Caused by warning:
## ! NAs introduced by coercion

fuel_Dead_Data_Avg <- fuel_Dead_Data_Avg %>%
  mutate(biomass = Dead_Weight_Post - Dry_DeadBag)

fuel_Dead_Data_Avg <- fuel_Dead_Data_Avg %>%
  filter(biomass >= 0)

fuel_Dead_Data_Avg <- fuel_Dead_Data_Avg %>%
  mutate(relative_moisture_content = ifelse(biomass > bioT, (Dead_Weight_Initial - Dead_Weight_Post) / biomass * 100, NA))

avg_dead_values_Avg <- fuel_Dead_Data_Avg %>%
  group_by(Plot, Status) %>%
  summarize(avg_dead_biomass = mean(biomass, na.rm = TRUE),
            avg_dead_moisture_content = mean(relative_moisture_content, na.rm = TRUE),
            avg_soil_moisture = mean(Soil_Moisture, na.rm = TRUE),
            .groups = "drop")

Litter

fuel_Litter_Data_Avg <- fuel_data_Avg %>%
  mutate(
    Litter_Bag = as.numeric(Litter_Bag),
    Litter_Weight_Post = as.numeric(Litter_Weight_Post),
    Litter_Weight_Initial = as.numeric(Litter_Weight_Initial),
    Dry_LitterBag = as.numeric(Dry_LitterBag),
    Status = as.character(Status)  # Status as a character
  )

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Dry_LitterBag = as.numeric(Dry_LitterBag)`.
## Caused by warning:
## ! NAs introduced by coercion

fuel_Litter_Data_Avg <- fuel_Litter_Data_Avg %>%
  mutate(biomass = Litter_Weight_Post - Dry_LitterBag)

fuel_Litter_Data_Avg <- fuel_Litter_Data_Avg %>%
  filter(biomass >= 0)

fuel_Litter_Data_Avg <- fuel_Litter_Data_Avg %>%
  mutate(relative_moisture_content = ifelse(biomass > bioT, (Litter_Weight_Initial - Litter_Weight_Post) / biomass * 100, NA))

avg_litter_values_Avg <- fuel_Litter_Data_Avg %>%
  group_by(Plot, Status) %>%
  summarize(avg_litter_biomass = mean(biomass, na.rm = TRUE),
            avg_litter_moisture_content = mean(relative_moisture_content, na.rm = TRUE),
            avg_soil_moisture = mean(Soil_Moisture, na.rm = TRUE),
            .groups = "drop")

Fill in gaps within Net method with Avg method.

# Live
avg_live_values_Combined <- avg_live_values_Net %>%
  full_join(avg_live_values_Avg, by = "Plot", suffix = c("_Net", "_Avg")) %>%
  mutate(
    avg_live_biomass = coalesce(avg_live_biomass_Net, avg_live_biomass_Avg),
    avg_live_moisture_content = coalesce(avg_live_moisture_content_Net, avg_live_moisture_content_Avg),
    avg_soil_moisture = coalesce(avg_soil_moisture_Net, avg_soil_moisture_Avg),
    avg_height = coalesce(avg_height_Net, avg_height_Avg),
    Status = coalesce(Status_Net, Status_Avg)
  ) %>%
  select(Plot, Status, avg_live_biomass, avg_live_moisture_content, avg_soil_moisture, avg_height)

# Dead
avg_dead_values_Combined <- avg_dead_values_Net %>%
  full_join(avg_dead_values_Avg, by = "Plot", suffix = c("_Net", "_Avg")) %>%
  mutate(
    avg_dead_biomass = coalesce(avg_dead_biomass_Net, avg_dead_biomass_Avg),
    avg_dead_moisture_content = coalesce(avg_dead_moisture_content_Net, avg_dead_moisture_content_Avg),
    avg_soil_moisture = coalesce(avg_soil_moisture_Net, avg_soil_moisture_Avg),
    Status = coalesce(Status_Net, Status_Avg)
  ) %>%
  select(Plot, Status, avg_dead_biomass, avg_dead_moisture_content, avg_soil_moisture)

# Litter
avg_litter_values_Combined <- avg_litter_values_Net %>%
  full_join(avg_litter_values_Avg, by = "Plot", suffix = c("_Net", "_Avg")) %>%
  mutate(
    avg_litter_biomass = coalesce(avg_litter_biomass_Net, avg_litter_biomass_Avg),
    avg_litter_moisture_content = coalesce(avg_litter_moisture_content_Net, avg_litter_moisture_content_Avg),
    avg_soil_moisture = coalesce(avg_soil_moisture_Net, avg_soil_moisture_Avg),
    Status = coalesce(Status_Net, Status_Avg)
  ) %>%
  select(Plot, Status, avg_litter_biomass, avg_litter_moisture_content, avg_soil_moisture)

Merge avg_live_values, avg_dead_values, and avg_litter_values

Average Fuel Values by Invasion Status and Season

There are 10,000 square meters in a hectare. Biomass is from 25 cm by 25 cm quadrats, so we have 0.0625 square meters. Therefore, 10,000/0.0625 = 160,000. So biomass gets multiplied by 160,000 and divided by 1,000,000 to convert from grams to tonnes.

25th, 50th and 75 quantiles of fuel values

Fuel_model_quantiles <- avg_fuel_values %>%
  group_by(Status) %>%
  summarize(avg_live_biomass_25 = quantile(avg_live_biomass, 0.25, na.rm = TRUE) * 0.16,
            avg_live_biomass_50 = quantile(avg_live_biomass, 0.50, na.rm = TRUE) * 0.16,
            avg_live_biomass_75 = quantile(avg_live_biomass, 0.75, na.rm = TRUE) * 0.16,
            avg_dead_biomass_25 = quantile(avg_dead_biomass, 0.25, na.rm = TRUE) * 0.16,
            avg_dead_biomass_50 = quantile(avg_dead_biomass, 0.50, na.rm = TRUE) * 0.16,
            avg_dead_biomass_75 = quantile(avg_dead_biomass, 0.75, na.rm = TRUE) * 0.16,
            avg_litter_biomass_25 = quantile(avg_litter_biomass, 0.25, na.rm = TRUE) * 0.16,
            avg_litter_biomass_50 = quantile(avg_litter_biomass, 0.50, na.rm = TRUE) * 0.16,
            avg_litter_biomass_75 = quantile(avg_litter_biomass, 0.75, na.rm = TRUE) * 0.16,
            avg_live_moisture_content_25 = quantile(avg_live_moisture_content, 0.25, na.rm = TRUE),
            avg_live_moisture_content_50 = quantile(avg_live_moisture_content, 0.50, na.rm = TRUE),
            avg_live_moisture_content_75 = quantile(avg_live_moisture_content, 0.75, na.rm = TRUE),
            avg_dead_moisture_content_25 = quantile(avg_dead_moisture_content, 0.25, na.rm = TRUE),
            avg_dead_moisture_content_50 = quantile(avg_dead_moisture_content, 0.50, na.rm = TRUE),
            avg_dead_moisture_content_75 = quantile(avg_dead_moisture_content, 0.75, na.rm = TRUE),
            avg_litter_moisture_content_25 = quantile(avg_litter_moisture_content, 0.25, na.rm = TRUE),
            avg_litter_moisture_content_50 = quantile(avg_litter_moisture_content, 0.50, na.rm = TRUE),
            avg_litter_moisture_content_75 = quantile(avg_litter_moisture_content, 0.75, na.rm = TRUE),
            avg_soil_moisture_25 = quantile(avg_soil_moisture, 0.25, na.rm = TRUE),
            avg_soil_moisture_50 = quantile(avg_soil_moisture, 0.50, na.rm = TRUE),
            avg_soil_moisture_75 = quantile(avg_soil_moisture, 0.75, na.rm = TRUE),
            avg_height_25 = quantile(avg_height, 0.25, na.rm = TRUE),
            avg_height_50 = quantile(avg_height, 0.50, na.rm = TRUE),
            avg_height_75 = quantile(avg_height, 0.75, na.rm = TRUE),
            .groups = "drop")

# Kable table of quantiles
kable(Fuel_model_quantiles)

Status	avg_live_biomass_25	avg_live_biomass_50	avg_live_biomass_75	avg_dead_biomass_25	avg_dead_biomass_50	avg_dead_biomass_75	avg_litter_biomass_25	avg_litter_biomass_50	avg_litter_biomass_75	avg_live_moisture_content_25	avg_live_moisture_content_50	avg_live_moisture_content_75	avg_dead_moisture_content_25	avg_dead_moisture_content_50	avg_dead_moisture_content_75	avg_litter_moisture_content_25	avg_litter_moisture_content_50	avg_litter_moisture_content_75	avg_soil_moisture_25	avg_soil_moisture_50	avg_soil_moisture_75	avg_height_25	avg_height_50	avg_height_75
Invaded	1.9420000	2.5928000	3.707200	1.4458667	2.441333	4.0536000	2.5756	4.253333	6.818267	133.3854	148.8094	167.7036	12.906657	17.98529	26.59409	12.39905	18.90885	32.23530	5.616667	10.366667	12.866667	0.6825000	0.8383333	1.0291667
Non_Invaded	0.3309333	0.6501333	1.038133	0.2074667	0.457600	0.8298667	3.8672	6.037333	7.357867	104.1553	140.3035	176.0684	7.881773	16.09442	24.74012	10.34319	17.65008	30.25146	3.866667	6.266667	8.933333	0.1533333	0.2366667	0.3366667