Bayesian Framework Example
Alan Ivory
2025-08-20
pacman::p_load(dplyr, readxl, tidyverse, raster, vegan, tigris, sf, sp, plotly, ggrepel, kableExtra, brms, parameters, gt)
## Set seed
set.seed(97)
# 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"))
# 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"))
Seasonal_Fuel_Sampling <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/01_FuelDynamics/02_Data/02_Fuel_Data/Seasonal_Fuel_Sampling.xlsx",
sheet = "Fuel_Data")
Seasonal_Fuel_Sampling <- Seasonal_Fuel_Sampling %>%
filter(Authority %in% c("BRSF", "WSF", "Jay"))
# Seasonal Sampling Locations
Seasonal_Sampling_Locations <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/01_FuelDynamics/02_Data/02_Fuel_Data/Seasonal_Fuel_Sampling.xlsx",
sheet = "Sites")
Seasonal_Sampling_Locations <- Seasonal_Sampling_Locations %>%
filter(Authority %in% c("BRSF", "WSF", "Jay"))
# Bag Weights
bag_weights <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/01_FuelDynamics/02_Data/02_Fuel_Data/Seasonal_Fuel_Sampling.xlsx",
sheet = "Bag_Avg")## New names:
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`# 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"))
# Only include Florida/Alabama Sites
CogonSites <- CogonSites[CogonSites$Authority != "CNF" & CogonSites$Authority != "DSNF", ]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
Dead
Litter
Average Bag Weight Method
Live
Dead
Litter
Fill in gaps within Net method with Avg method.
Combine summer sampling (CogonSites) locations with winter sampling locations
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 Table
Objective 1- Cogongrass relationship with fuel dynamics
Data Cleaning
GLMM Live Fuel Biomass
Model Comparison
## Warning: Found 1 observations with a pareto_k > 0.7 in model
## 'LV_Bio_model_student'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.## Warning: There were 15 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: avg_live_biomass ~ Status * Season + Status * Region + (1 | Plot)
## Data: merged_sites (Number of observations: 183)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 183)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.19 0.14 0.01 0.52 1.03 125 34
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.64 0.12 0.41 0.89 1.00 1690
## StatusInvaded 1.36 0.20 0.98 1.76 1.00 1048
## SeasonSpring -0.04 0.18 -0.38 0.32 1.01 673
## SeasonWinter -0.19 0.17 -0.52 0.15 1.01 456
## RegionCF -0.05 0.14 -0.32 0.23 1.01 1317
## StatusInvaded:SeasonSpring -0.67 0.31 -1.24 -0.04 1.00 1450
## StatusInvaded:SeasonWinter -0.91 0.29 -1.46 -0.33 1.00 1364
## StatusInvaded:RegionCF 1.21 0.26 0.70 1.72 1.01 892
## Tail_ESS
## Intercept 1773
## StatusInvaded 975
## SeasonSpring 776
## SeasonWinter 523
## RegionCF 1611
## StatusInvaded:SeasonSpring 2365
## StatusInvaded:SeasonWinter 2394
## StatusInvaded:RegionCF 1616
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.58 0.09 0.36 0.76 1.03 125 34
## nu 2.48 0.74 1.44 4.28 1.01 278 176
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
Live Biomass Summary Table
Conditional Effects Live Biomass
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
Conditional Effects Live Biomass
Live Biomass Posterior Probabilities
GLMM Dead Fuel Biomass
Model Comparison
## Warning: There were 7 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: avg_dead_biomass ~ Status * Season + Status * Region + (1 | Plot)
## Data: merged_sites (Number of observations: 181)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 181)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.19 0.15 0.01 0.57 1.00 939 886
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.39 0.18 0.04 0.75 1.00 3000
## StatusInvaded 1.70 0.31 1.12 2.31 1.00 2229
## SeasonSpring 0.40 0.29 -0.16 0.98 1.00 3237
## SeasonWinter 0.39 0.28 -0.13 0.97 1.00 3074
## RegionCF 0.14 0.22 -0.29 0.58 1.00 3622
## StatusInvaded:SeasonSpring -0.30 0.47 -1.25 0.56 1.00 2693
## StatusInvaded:SeasonWinter 0.70 0.63 -0.59 1.95 1.00 2883
## StatusInvaded:RegionCF 0.91 0.52 -0.11 1.93 1.00 3177
## Tail_ESS
## Intercept 2880
## StatusInvaded 1860
## SeasonSpring 3117
## SeasonWinter 2462
## RegionCF 2642
## StatusInvaded:SeasonSpring 2269
## StatusInvaded:SeasonWinter 2825
## StatusInvaded:RegionCF 2763
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.93 0.15 0.66 1.24 1.00 1620 1459
## nu 1.89 0.49 1.19 3.03 1.00 2426 1717
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
Dead Biomass Summary Table
Conditional Effects Dead Biomass
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?