Chapter4 Objective1
Alan Ivory
2025-06-13
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"))
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", "BCNWR", "Escambia", "Hooper", "Howell", "ONF"))
# 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
Comb_Live_Data_Net <- Combined_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 coercionComb_Live_Data_Net <- Comb_Live_Data_Net %>%
mutate(biomass = Net_Live)
Comb_Live_Data_Net <- Comb_Live_Data_Net %>%
filter(biomass >= 0)
Comb_Live_Data_Net <- Comb_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 <- Comb_Live_Data_Net %>%
group_by(Plot, Season, 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
Comb_Dead_Data_Net <- Combined_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
)
Comb_Dead_Data_Net <- Comb_Dead_Data_Net %>%
mutate(biomass = Net_Dead)
Comb_Dead_Data_Net <- Comb_Dead_Data_Net %>%
filter(biomass >= 0)
Comb_Dead_Data_Net <- Comb_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 <- Comb_Dead_Data_Net %>%
group_by(Plot, Season, 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
Comb_Litter_Data_Net <- Combined_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
)
Comb_Litter_Data_Net <- Comb_Litter_Data_Net %>%
mutate(biomass = Net_Litter)
Comb_Litter_Data_Net <- Comb_Litter_Data_Net %>%
filter(biomass >= 0)
Comb_Litter_Data_Net <- Comb_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 <- Comb_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
Comb_Live_Data_Avg <- Combined_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.Comb_Live_Data_Avg <- Comb_Live_Data_Avg %>%
mutate(biomass = Live_Weight_Post - Dry_LiveBag)
Comb_Live_Data_Avg <- Comb_Live_Data_Avg %>%
filter(biomass >= 0)
Comb_Live_Data_Avg <- Comb_Live_Data_Avg %>%
mutate(relative_moisture_content = ifelse(biomass > bioT, (Live_Weight_Initial - Live_Weight_Post) / biomass * 100, NA))
avg_live_values_Avg <- Comb_Live_Data_Avg %>%
group_by(Plot, Season, 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
Comb_Dead_Data_Avg <- Combined_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 coercionComb_Dead_Data_Avg <- Comb_Dead_Data_Avg %>%
mutate(biomass = Dead_Weight_Post - Dry_DeadBag)
Comb_Dead_Data_Avg <- Comb_Dead_Data_Avg %>%
filter(biomass >= 0)
Comb_Dead_Data_Avg <- Comb_Dead_Data_Avg %>%
mutate(relative_moisture_content = ifelse(biomass > bioT, (Dead_Weight_Initial - Dead_Weight_Post) / biomass * 100, NA))
avg_dead_values_Avg <- Comb_Dead_Data_Avg %>%
group_by(Plot, Season, 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
Comb_Litter_Data_Avg <- Combined_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 coercionComb_Litter_Data_Avg <- Comb_Litter_Data_Avg %>%
mutate(biomass = Litter_Weight_Post - Dry_LitterBag)
Comb_Litter_Data_Avg <- Comb_Litter_Data_Avg %>%
filter(biomass >= 0)
Comb_Litter_Data_Avg <- Comb_Litter_Data_Avg %>%
mutate(relative_moisture_content = ifelse(biomass > bioT, (Litter_Weight_Initial - Litter_Weight_Post) / biomass * 100, NA))
avg_litter_values_Avg <- Comb_Litter_Data_Avg %>%
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")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),
Season = coalesce(Season_Net, Season_Avg),
Status = coalesce(Status_Net, Status_Avg)
) %>%
select(Plot, Season, 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),
Season = coalesce(Season_Net, Season_Avg),
Status = coalesce(Status_Net, Status_Avg)
) %>%
select(Plot, Season, 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),
Season = coalesce(Season_Net, Season_Avg),
Status = coalesce(Status_Net, Status_Avg)
) %>%
select(Plot, Season, Status, avg_litter_biomass, avg_litter_moisture_content, avg_soil_moisture)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_model_quantiles <- avg_fuel_values %>%
group_by(Status, Season) %>%
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 | Season | 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 | Green_Up | 1.1562667 | 2.161600 | 3.4720000 | 1.8922667 | 2.4192000 | 3.1082667 | 3.486400 | 4.974400 | 9.215467 | 54.66557 | 117.06158 | 133.0470 | 9.195465 | 10.39843 | 16.57708 | 9.487179 | 11.44906 | 14.92138 | 4.833333 | 5.766667 | 7.100000 | 0.5600000 | 0.7533333 | 0.9500000 |
| Invaded | Summer | 1.9560000 | 2.609067 | 3.7072000 | 1.3474667 | 2.4234667 | 4.0458667 | 2.652533 | 4.274133 | 6.774667 | 133.52828 | 148.91918 | 167.7775 | 13.037626 | 17.99824 | 27.12177 | 12.411588 | 19.05224 | 32.37791 | 5.833333 | 10.366667 | 13.016667 | 0.6850000 | 0.8400000 | 1.0300000 |
| Invaded | Winter | 0.9458667 | 1.497600 | 2.4898667 | 1.9476000 | 3.4101333 | 6.1982667 | 3.023867 | 5.620267 | 7.916800 | 95.76815 | 119.49396 | 151.1411 | 13.788625 | 19.72430 | 36.73561 | 15.979345 | 27.86209 | 47.96388 | 5.933333 | 10.916667 | 13.891667 | 0.7116667 | 0.8500000 | 0.9733333 |
| Non_Invaded | Green_Up | 0.3685333 | 0.508800 | 0.7546667 | 0.3370667 | 0.7482667 | 1.3616000 | 3.924267 | 6.187200 | 8.342933 | 51.67528 | 80.33146 | 146.1973 | 9.750584 | 15.46700 | 17.59368 | 7.290356 | 12.30107 | 17.27920 | 3.066667 | 4.033333 | 9.300000 | 0.2333333 | 0.3333333 | 0.4366667 |
| Non_Invaded | Summer | 0.3402667 | 0.636800 | 1.0261333 | 0.2090667 | 0.4392000 | 0.8254667 | 3.882600 | 6.036267 | 7.348533 | 104.15527 | 140.30346 | 176.0684 | 7.881773 | 16.09442 | 24.74012 | 10.345724 | 17.65342 | 30.36058 | 3.991667 | 6.333333 | 9.108333 | 0.1533333 | 0.2366667 | 0.3366667 |
| Non_Invaded | Winter | 0.2112000 | 0.272800 | 0.5125333 | 0.3168000 | 0.7242667 | 1.8066667 | 4.612133 | 6.603733 | 9.257733 | 30.64949 | 85.18529 | 184.4349 | 11.679660 | 16.25337 | 49.32975 | 15.910721 | 42.65335 | 56.80946 | 6.233333 | 9.566667 | 12.375000 | 0.1116667 | 0.2366667 | 0.3583333 |
Objective 1- Cogongrass relationship with fuel dynamics
Data Cleaning
GLMM Live Fuel Biomass
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_LV_Bio <- loo(LV_Bio_model_gaussian)## Warning: Found 108 observations with a pareto_k > 0.7 in model
## 'LV_Bio_model_gaussian'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.loo_student_LV_Bio <- loo(LV_Bio_model_student)## Warning: Found 9 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.# Model comparison
loo_compare(loo_gaussian_LV_Bio, loo_student_LV_Bio)## elpd_diff se_diff
## LV_Bio_model_student 0.0 0.0
## LV_Bio_model_gaussian -24.4 16.2summary(LV_Bio_model_student)## Warning: There were 123 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 + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 246)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 246)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.20 0.14 0.01 0.50 1.04 108 50
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.25 0.15 -0.05 0.54 1.00 1992
## StatusInvaded 1.02 0.23 0.58 1.50 1.00 1918
## SeasonGreen_Up 0.17 0.21 -0.23 0.58 1.01 299
## SeasonSummer 0.25 0.16 -0.06 0.57 1.01 1835
## RegionCF 0.34 0.12 0.11 0.59 1.00 1167
## RegionSM 0.44 0.16 0.13 0.75 1.00 2039
## StatusInvaded:SeasonGreen_Up 0.19 0.37 -0.52 0.95 1.00 1913
## StatusInvaded:SeasonSummer 0.86 0.27 0.31 1.37 1.01 1455
## Tail_ESS
## Intercept 1937
## StatusInvaded 2329
## SeasonGreen_Up 642
## SeasonSummer 2121
## RegionCF 1300
## RegionSM 1486
## StatusInvaded:SeasonGreen_Up 1744
## StatusInvaded:SeasonSummer 1816
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.59 0.08 0.43 0.75 1.01 230 311
## nu 2.68 0.69 1.64 4.41 1.01 921 872
##
## 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).# Plot the posterior predictive checks
pp_check(LV_Bio_model_student) +
labs(title = "Posterior Predictive Check for Live Biomass Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Rearrange "High", "Mean", "Low" to "Low", "Mean", "High"
# Plot the conditional effects
#conditional_effects(LV_Bio_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Live Biomass Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_live_biomass, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Live Biomass by Invasion Status and Season",
x = "Invasion Status",
y = "Live Biomass (tonne/ha)",
color = "Season") +
theme_minimal()## 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()`).
GLMM Dead Fuel Biomass
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_D_Bio <- loo(D_Bio_model_gaussian)## Warning: Found 86 observations with a pareto_k > 0.7 in model
## 'D_Bio_model_gaussian'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.loo_student_D_Bio <- loo(D_Bio_model_student)
# Model comparison
loo_compare(loo_gaussian_D_Bio, loo_student_D_Bio)## elpd_diff se_diff
## D_Bio_model_student 0.0 0.0
## D_Bio_model_gaussian -34.4 13.2summary(D_Bio_model_student)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: avg_dead_biomass ~ Status * Season + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 241)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 241)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.13 0.10 0.00 0.36 1.01 1211 1524
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.59 0.21 0.20 1.03 1.00 4284
## StatusInvaded 2.40 0.60 1.34 3.74 1.00 1924
## SeasonGreen_Up 0.08 0.31 -0.53 0.70 1.00 3913
## SeasonSummer -0.30 0.22 -0.76 0.12 1.00 3885
## RegionCF 0.31 0.17 -0.01 0.65 1.00 4009
## RegionSM 0.44 0.18 0.10 0.80 1.00 5165
## StatusInvaded:SeasonGreen_Up -0.72 0.69 -2.17 0.55 1.00 1995
## StatusInvaded:SeasonSummer -0.61 0.62 -1.98 0.48 1.00 2040
## Tail_ESS
## Intercept 3317
## StatusInvaded 1984
## SeasonGreen_Up 2859
## SeasonSummer 2988
## RegionCF 3151
## RegionSM 2968
## StatusInvaded:SeasonGreen_Up 2210
## StatusInvaded:SeasonSummer 1767
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.72 0.10 0.53 0.93 1.00 1922 1544
## nu 1.55 0.29 1.08 2.25 1.00 2061 1317
##
## 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).# Plot the posterior predictive checks
pp_check(D_Bio_model_student) +
labs(title = "Posterior Predictive Check for Dead Biomass Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Rearrange "High", "Mean", "Low" to "Low", "Mean", "High"
# Plot the conditional effects
#conditional_effects(D_Bio_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Dead Biomass Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_dead_biomass, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Dead Biomass by Invasion Status and Season",
x = "Invasion Status",
y = "Dead Biomass (tonne/ha)",
color = "Season") +
theme_minimal()## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
GLMM Litter Fuel Biomass
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_LT_Bio <- loo(LT_Bio_model_gaussian)## Warning: Found 57 observations with a pareto_k > 0.7 in model
## 'LT_Bio_model_gaussian'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.loo_student_LT_Bio <- loo(LT_Bio_model_student)## Warning: Found 5 observations with a pareto_k > 0.7 in model
## 'LT_Bio_model_student'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.# Model comparison
loo_compare(loo_gaussian_LT_Bio, loo_student_LT_Bio)## elpd_diff se_diff
## LT_Bio_model_student 0.0 0.0
## LT_Bio_model_gaussian -3.4 5.1summary(LT_Bio_model_student)## Warning: There were 59 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_litter_biomass ~ Status * Season + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 249)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 249)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.11 0.65 0.07 2.30 1.04 83 221
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.04 0.72 5.64 8.45 1.00 1063
## StatusInvaded -0.83 0.97 -2.75 1.05 1.01 847
## SeasonGreen_Up -0.30 1.05 -2.30 1.83 1.01 651
## SeasonSummer -0.98 0.76 -2.46 0.49 1.01 1141
## RegionCF -0.60 0.44 -1.47 0.27 1.00 2200
## RegionSM -0.33 0.65 -1.54 0.95 1.00 2069
## StatusInvaded:SeasonGreen_Up 0.29 1.47 -2.59 3.10 1.01 713
## StatusInvaded:SeasonSummer -0.21 1.09 -2.29 1.96 1.01 898
## Tail_ESS
## Intercept 1663
## StatusInvaded 1643
## SeasonGreen_Up 738
## SeasonSummer 1952
## RegionCF 1964
## RegionSM 1652
## StatusInvaded:SeasonGreen_Up 931
## StatusInvaded:SeasonSummer 1563
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 2.52 0.41 1.60 3.20 1.04 93 110
## nu 5.97 3.64 2.46 15.52 1.02 273 439
##
## 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).# Plot the posterior predictive checks
pp_check(LT_Bio_model_student) +
labs(title = "Posterior Predictive Check for Live Biomass Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# Plot the conditional effects
#conditional_effects(LT_Bio_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of litter Biomass Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_litter_biomass, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "litter Biomass by Invasion Status and Season",
x = "Invasion Status",
y = "Dead Biomass (tonne/ha)",
color = "Season") +
theme_minimal()
GLMM Live Moisture
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_LV_Moisture <- loo(LV_Moisture_model_gaussian)## Warning: Found 79 observations with a pareto_k > 0.7 in model
## 'LV_Moisture_model_gaussian'. We recommend to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.loo_student_LV_Moisture <- loo(LV_Moisture_model_student)## Warning: Found 12 observations with a pareto_k > 0.7 in model
## 'LV_Moisture_model_student'. We recommend to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.# Model comparison
loo_compare(loo_gaussian_LV_Moisture, loo_student_LV_Moisture)## elpd_diff se_diff
## LV_Moisture_model_gaussian 0.0 0.0
## LV_Moisture_model_student -9.3 3.1summary(LV_Moisture_model_student)## Warning: There were 1 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_moisture_content ~ Status * Season + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 187)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 14.43 9.46 0.67 33.98 1.04 101 237
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 103.95 28.42 49.05 162.66 1.00 1248
## StatusInvaded 19.87 29.62 -41.90 77.26 1.00 1239
## SeasonGreen_Up -9.44 32.29 -76.11 49.91 1.00 1356
## SeasonSummer 46.60 28.70 -12.31 102.28 1.00 1245
## RegionCF -3.40 7.51 -17.99 11.26 1.00 3458
## RegionSM -26.66 9.52 -45.10 -8.33 1.00 3439
## StatusInvaded:SeasonGreen_Up -9.20 35.87 -75.55 64.64 1.00 1367
## StatusInvaded:SeasonSummer -13.11 30.74 -72.02 49.81 1.00 1246
## Tail_ESS
## Intercept 1492
## StatusInvaded 1415
## SeasonGreen_Up 1539
## SeasonSummer 1671
## RegionCF 2645
## RegionSM 2954
## StatusInvaded:SeasonGreen_Up 1395
## StatusInvaded:SeasonSummer 1369
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 35.87 5.80 22.79 45.35 1.03 153 201
## nu 9.72 8.12 2.75 31.56 1.00 868 1141
##
## 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).# Plot the posterior predictive checks
pp_check(LV_Moisture_model_student) +
labs(title = "Posterior Predictive Check for Live Moisture Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# Plot the conditional effects
#conditional_effects(LV_Moisture_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Live Moisture Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_live_moisture_content, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Live Moisture by Invasion Status and Season",
x = "Invasion Status",
y = "Live Moisture (%)",
color = "Season") +
theme_minimal()## Warning: Removed 62 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 62 rows containing missing values or values outside the scale range
## (`geom_point()`).
GLMM Dead Moisture
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_D_Moisture <- loo(D_Moisture_model_gaussian)## Warning: Found 67 observations with a pareto_k > 0.7 in model
## 'D_Moisture_model_gaussian'. We recommend to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.loo_student_D_Moisture <- loo(D_Moisture_model_student)
# Model comparison
loo_compare(loo_gaussian_D_Moisture, loo_student_D_Moisture)## elpd_diff se_diff
## D_Moisture_model_student 0.0 0.0
## D_Moisture_model_gaussian -76.8 15.6summary(D_Moisture_model_student)## Warning: There were 8 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_moisture_content ~ Status * Season + 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) 1.94 1.26 0.09 4.64 1.01 374 288
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 15.60 2.32 11.26 20.42 1.00 1383
## StatusInvaded 3.65 2.84 -1.87 9.32 1.00 1189
## SeasonGreen_Up -1.48 2.93 -7.41 4.22 1.00 1467
## SeasonSummer -2.86 2.67 -8.06 2.41 1.00 1534
## RegionCF -1.71 1.43 -4.51 1.04 1.00 2343
## RegionSM 3.97 2.08 -0.05 8.08 1.00 1553
## StatusInvaded:SeasonGreen_Up -5.16 3.79 -12.27 2.42 1.00 1234
## StatusInvaded:SeasonSummer 0.20 3.43 -6.57 6.69 1.00 1335
## Tail_ESS
## Intercept 2152
## StatusInvaded 1852
## SeasonGreen_Up 1890
## SeasonSummer 2205
## RegionCF 2312
## RegionSM 807
## StatusInvaded:SeasonGreen_Up 1926
## StatusInvaded:SeasonSummer 1746
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.67 0.86 4.06 7.49 1.01 483 214
## nu 1.30 0.19 1.02 1.74 1.01 681 436
##
## 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).# Plot the posterior predictive checks
pp_check(D_Moisture_model_student) +
labs(title = "Posterior Predictive Check for Dead Moisture Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# Plot the conditional effects
#conditional_effects(D_Moisture_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Dead Moisture Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_dead_moisture_content, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Dead Moisture by Invasion Status and Season",
x = "Invasion Status",
y = "Dead Moisture (%)",
color = "Season") +
theme_minimal()## Warning: Removed 66 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_point()`).
GLMM Litter Moisture
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_LT_Moisture <- loo(LT_Moisture_model_gaussian)## Warning: Found 107 observations with a pareto_k > 0.7 in model
## 'LT_Moisture_model_gaussian'. We recommend to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.loo_student_LT_Moisture <- loo(LT_Moisture_model_student)
# Model comparison
loo_compare(loo_gaussian_LT_Moisture, loo_student_LT_Moisture)## elpd_diff se_diff
## LT_Moisture_model_student 0.0 0.0
## LT_Moisture_model_gaussian -30.4 13.8summary(LT_Moisture_model_student)## Warning: There were 2 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_litter_moisture_content ~ Status * Season + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 248)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 248)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.50 1.79 0.09 6.82 1.01 438 310
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 40.37 5.67 28.03 50.40 1.00 1487
## StatusInvaded -13.40 6.59 -25.30 0.22 1.00 1290
## SeasonGreen_Up -25.65 6.30 -36.80 -12.24 1.00 1496
## SeasonSummer -22.77 6.18 -33.57 -9.81 1.00 1264
## RegionCF -4.32 1.98 -8.14 -0.51 1.00 4958
## RegionSM -0.48 2.45 -5.35 4.24 1.00 2757
## StatusInvaded:SeasonGreen_Up 12.81 7.45 -2.64 26.52 1.00 1313
## StatusInvaded:SeasonSummer 15.24 6.93 0.84 28.14 1.00 1318
## Tail_ESS
## Intercept 2029
## StatusInvaded 1990
## SeasonGreen_Up 2076
## SeasonSummer 1043
## RegionCF 3098
## RegionSM 2450
## StatusInvaded:SeasonGreen_Up 2003
## StatusInvaded:SeasonSummer 1791
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 9.26 1.27 6.81 11.78 1.00 1023 988
## nu 1.52 0.26 1.11 2.12 1.00 1808 1726
##
## 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).# Plot the posterior predictive checks
pp_check(LT_Moisture_model_student) +
labs(title = "Posterior Predictive Check for Litter Moisture Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# Plot the conditional effects
#conditional_effects(LT_Moisture_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Litter Moisture Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_litter_moisture_content, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Litter Moisture by Invasion Status and Season",
x = "Invasion Status",
y = "Litter Moisture (%)",
color = "Season") +
theme_minimal()## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
GLMM Soil Moisture
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_Soil_Moisture <- loo(Soil_Moisture_model_gaussian)## Warning: Found 82 observations with a pareto_k > 0.7 in model
## 'Soil_Moisture_model_gaussian'. We recommend to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.loo_student_Soil_Moisture <- loo(Soil_Moisture_model_student)## Warning: Found 2 observations with a pareto_k > 0.7 in model
## 'Soil_Moisture_model_student'. We recommend to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.# Model comparison
loo_compare(loo_gaussian_Soil_Moisture, loo_student_Soil_Moisture)## elpd_diff se_diff
## Soil_Moisture_model_student 0.0 0.0
## Soil_Moisture_model_gaussian -81.6 20.0summary(Soil_Moisture_model_student)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: avg_soil_moisture ~ Status * Season + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 249)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 249)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.13 0.88 0.19 3.35 1.05 70 278
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 11.07 0.90 9.33 12.81 1.00 943
## StatusInvaded 0.25 1.24 -2.18 2.57 1.00 812
## SeasonGreen_Up -2.91 1.32 -5.50 -0.34 1.00 1120
## SeasonSummer -3.20 0.96 -4.99 -1.31 1.00 1032
## RegionCF -4.76 0.57 -5.85 -3.65 1.00 1558
## RegionSM 0.30 0.80 -1.29 1.86 1.00 1719
## StatusInvaded:SeasonGreen_Up -0.57 1.79 -4.10 2.90 1.00 916
## StatusInvaded:SeasonSummer 2.33 1.38 -0.31 5.04 1.00 926
## Tail_ESS
## Intercept 1380
## StatusInvaded 1503
## SeasonGreen_Up 1857
## SeasonSummer 1908
## RegionCF 1894
## RegionSM 1260
## StatusInvaded:SeasonGreen_Up 1602
## StatusInvaded:SeasonSummer 1758
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 2.09 0.69 0.94 3.36 1.05 66 183
## nu 1.68 0.42 1.04 2.64 1.04 84 175
##
## 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).# Plot the posterior predictive checks
pp_check(Soil_Moisture_model_student) +
labs(title = "Posterior Predictive Check for Soil Moisture Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# Plot the conditional effects
#conditional_effects(Soil_Moisture_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Soil Moisture Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_soil_moisture, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Soil Moisture by Invasion Status and Season",
x = "Invasion Status",
y = "Soil Moisture (%)",
color = "Season") +
theme_minimal()
GLMM Veg Height
Model Comparison
# Compare models with Leave-One-Out Cross-Validation
loo_gaussian_Veg_Height <- loo(Veg_Height_model_gaussian)## Warning: Found 68 observations with a pareto_k > 0.7 in model
## 'Veg_Height_model_gaussian'. We recommend to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.loo_student_Veg_Height <- loo(Veg_Height_model_student)## Warning: Found 19 observations with a pareto_k > 0.7 in model
## 'Veg_Height_model_student'. We recommend to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.# Model comparison
loo_compare(loo_gaussian_Veg_Height, loo_student_Veg_Height)## elpd_diff se_diff
## Veg_Height_model_gaussian 0.0 0.0
## Veg_Height_model_student -1.1 3.5summary(Veg_Height_model_student)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 333 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_height ~ Status * Season + Region + (1 | Plot)
## Data: merged_sites (Number of observations: 246)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Plot (Number of levels: 246)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.04 0.00 0.12 1.06 48 54
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.17 0.04 0.09 0.26 1.06 45
## StatusInvaded 0.59 0.06 0.48 0.70 1.05 49
## SeasonGreen_Up 0.11 0.06 0.00 0.22 1.06 45
## SeasonSummer 0.02 0.05 -0.07 0.11 1.07 35
## RegionCF 0.15 0.02 0.10 0.19 1.02 1094
## RegionSM 0.07 0.04 0.01 0.15 1.01 587
## StatusInvaded:SeasonGreen_Up -0.18 0.08 -0.34 -0.03 1.04 68
## StatusInvaded:SeasonSummer 0.01 0.06 -0.11 0.15 1.06 45
## Tail_ESS
## Intercept 208
## StatusInvaded 173
## SeasonGreen_Up 386
## SeasonSummer 223
## RegionCF 981
## RegionSM 2237
## StatusInvaded:SeasonGreen_Up 951
## StatusInvaded:SeasonSummer 61
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.15 0.02 0.11 0.18 1.07 35 110
## nu 8.94 6.01 3.21 26.29 1.04 69 133
##
## 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).# Plot the posterior predictive checks
pp_check(Veg_Height_model_student) +
labs(title = "Posterior Predictive Check for Veg Height Model") +
theme_minimal()## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# Plot the conditional effects
#conditional_effects(Veg_Height_model_student) %>%
# plot() +
# labs(title = "Conditional Effects of Veg Height Model") +
# theme_minimal()
# Plot the relationship between invasion status and surface rate of spread
merged_sites %>%
ggplot(aes(x = Status, y = avg_height, color = Season)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.2) +
labs(title = "Veg Height by Invasion Status and Season",
x = "Invasion Status",
y = "Veg Height (m)",
color = "Season") +
theme_minimal()## 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()`).