Quadrat Diversity
Alan Ivory
2026-03-03
pacman::p_load(tidyverse, readxl, raster, vegan, ggplot2,
tigris, sf, sp, plotly, ggrepel, rstanarm, mgcv,
brms, lme4, caret, MuMIn, VSURF, loo, patchwork)
EnvironmentalOutputs <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/03_Biodiversity/06_Processing/Environmental_Outputs.xlsx")
# 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")
# Veg Data
Veg_Cover <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Veg_Cover")
# Shrub Cover Data
shrub_data <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/Field_Data_FL_AL_MS.xlsx",
sheet = "Shrub_Cover")
# Site Data
CogonSites <- read_excel("C:/Users/DrewIvory/OneDrive - University of Florida/Desktop/School/PHD/01_Projects/05_SharedData/CogonSites_FL_AL_MS.xlsx")Number of quadrats sampled per plot
quadrat_count <- Veg_Cover %>%
group_by(Plot) %>%
summarize(total_quadrats = n_distinct(Quadrat), .groups = "drop")Filter all data to only include specified species (Per PLANTS database)
tree_data <- dplyr::filter(tree_data, Growth == "Tree")
Veg_Cover <- dplyr::filter(Veg_Cover, Growth != "Shrub" & Growth != "Tree")
shrub_data <- dplyr::filter(shrub_data, Growth == "Shrub" | Growth == "Tree")Join site info into Veg_Cover
Veg_Cover <- Veg_Cover %>%
mutate(Plot_Quadrat = paste(Plot, Quadrat, sep = "_")) %>%
left_join(
CogonSites %>% select(-Site),
by = "Plot"
)QUADRAT-LEVEL ANALYSIS
Extract cogongrass cover at the quadrat level
cogon_quadrat <- Veg_Cover %>%
filter(Species_Name %in% c("Imperata_cylindrica")) %>%
group_by(Plot, Quadrat) %>%
summarise(Cogongrass_Cover_Quadrat = sum(Cover_Per, na.rm = TRUE), .groups = "drop")Build quadrat-level species matrix (herbaceous, native only, no cogongrass)
veg_quadrat <- Veg_Cover %>%
filter(
!Species_Name %in% c("Imperata_cylindrica", "Imperata_cylindrica_Live"),
InvStatus == "Native"
) %>%
group_by(Plot, Quadrat, Species_Name) %>%
summarise(Cover_Per = sum(Cover_Per, na.rm = TRUE), .groups = "drop")
# Wide format
species_matrix_quadrat <- veg_quadrat %>%
pivot_wider(
names_from = Species_Name,
values_from = Cover_Per,
values_fill = 0
)Calculate Shannon diversity per quadrat
shannon_quadrat <- species_matrix_quadrat %>%
select(-Plot, -Quadrat) %>%
vegan::diversity(index = "shannon") %>%
as.data.frame() %>%
setNames("Shannon_Diversity") %>%
mutate(
Plot = species_matrix_quadrat$Plot,
Quadrat = species_matrix_quadrat$Quadrat
)Merge quadrat Shannon with cogongrass cover and site info
plot_quadrat_site <- Veg_Cover %>%
distinct(Plot, Quadrat, Site)
quadrat_model_data <- shannon_quadrat %>%
left_join(plot_quadrat_site, by = c("Plot", "Quadrat")) %>%
left_join(cogon_quadrat, by = c("Plot", "Quadrat")) %>%
mutate(Cogongrass_Cover_Quadrat = ifelse(
is.na(Cogongrass_Cover_Quadrat), 0, Cogongrass_Cover_Quadrat)
) %>%
left_join(
EnvironmentalOutputs %>% select(-any_of("Site")),
by = "Plot"
)
quadrat_model_data$Site <- as.factor(quadrat_model_data$Site)
summary(quadrat_model_data)## Shannon_Diversity Plot Quadrat Site
## Min. :0.0000 Length:1407 Min. : 0.00 WSF :298
## 1st Qu.:0.8574 Class :character 1st Qu.: 5.00 Jay :286
## Median :1.3863 Mode :character Median :15.00 BCNWR :217
## Mean :1.2778 Mean :14.98 DSNF :165
## 3rd Qu.:1.7380 3rd Qu.:25.00 CNF :156
## Max. :2.9444 Max. :30.00 Escambia: 77
## (Other) :208
## Cogongrass_Cover_Quadrat BurnYear Camera Latitude
## Min. : 0.00 Min. :2000 Length:1407 Min. :28.48
## 1st Qu.: 0.00 1st Qu.:2014 Class :character 1st Qu.:29.20
## Median : 0.00 Median :2022 Mode :character Median :30.77
## Mean :16.34 Mean :2018 Mean :30.38
## 3rd Qu.:15.00 3rd Qu.:2022 3rd Qu.:30.98
## Max. :97.50 Max. :2024 Max. :31.59
## NA's :287
## Longitude Muname NLCD_LandCover Region
## Min. :-89.80 Length:1407 Length:1407 Length:1407
## 1st Qu.:-88.91 Class :character Class :character Class :character
## Median :-87.14 Mode :character Mode :character Mode :character
## Mean :-86.72
## 3rd Qu.:-82.50
## Max. :-81.94
##
## SoilType Status TreeCanopy_NLCD aspect
## Length:1407 Length:1407 Min. : 0.00 Min. : 0.0
## Class :character Class :character 1st Qu.:53.00 1st Qu.: 90.0
## Mode :character Mode :character Median :76.00 Median :180.0
## Mean :66.54 Mean :171.0
## 3rd Qu.:87.00 3rd Qu.:270.0
## Max. :98.00 Max. :350.6
##
## elevation ndvi pr slope
## Min. : 17.00 Min. :0.2728 Min. :3.850 Min. : 0.000
## 1st Qu.: 38.00 1st Qu.:0.4320 1st Qu.:4.286 1st Qu.: 2.147
## Median : 59.00 Median :0.4636 Median :4.592 Median : 2.983
## Mean : 57.71 Mean :0.4507 Mean :4.566 Mean : 3.584
## 3rd Qu.: 75.00 3rd Qu.:0.4847 3rd Qu.:4.796 3rd Qu.: 4.684
## Max. :111.00 Max. :0.5357 Max. :5.294 Max. :18.429
##
## sph srad tmmn agb
## Min. :0.01051 Min. :204.7 Min. :285.7 Min. : 5.00
## 1st Qu.:0.01074 1st Qu.:206.8 1st Qu.:286.6 1st Qu.: 74.00
## Median :0.01091 Median :207.0 Median :287.2 Median : 95.00
## Mean :0.01145 Mean :210.5 Mean :287.4 Mean : 99.77
## 3rd Qu.:0.01194 3rd Qu.:218.5 3rd Qu.:289.0 3rd Qu.:133.00
## Max. :0.01325 Max. :224.2 Max. :290.3 Max. :211.00
## Quadrat-level model fitting: Shannon ~ Cogongrass
# Linear mixed model: Plot nested in Site as random effects
q_lmer_linear <- lmer(
Shannon_Diversity ~ Cogongrass_Cover_Quadrat + (1 | Site),
data = quadrat_model_data, REML = FALSE
)
# Quadratic term
q_lmer_quad <- lmer(
Shannon_Diversity ~ Cogongrass_Cover_Quadrat + I(Cogongrass_Cover_Quadrat^2) +
(1 | Site),
data = quadrat_model_data, REML = FALSE
)## Warning: Some predictor variables are on very different scales: consider
## rescaling# GAM with smoothed cogon term and random site effect
q_gam <- gam(
Shannon_Diversity ~ s(Cogongrass_Cover_Quadrat, k = 7) + s(Site, bs = "re"),
data = quadrat_model_data, method = "ML"
)
# Compare models
AIC(q_lmer_linear, q_lmer_quad)## df AIC
## q_lmer_linear 4 2317.000
## q_lmer_quad 5 2310.837summary(q_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Shannon_Diversity ~ s(Cogongrass_Cover_Quadrat, k = 7) + s(Site,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.19180 0.09517 12.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Cogongrass_Cover_Quadrat) 3.670 4.354 29.11 <2e-16 ***
## s(Site) 8.644 9.000 32.50 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.222 Deviance explained = 22.9%
## -ML = 1150.5 Scale est. = 0.29081 n = 1407Quadratic Linear Model Plot
# Prediction grid
cogon_seq <- data.frame(
Cogongrass_Cover_Quadrat = seq(
min(quadrat_model_data$Cogongrass_Cover_Quadrat, na.rm = TRUE),
max(quadrat_model_data$Cogongrass_Cover_Quadrat, na.rm = TRUE),
length.out = 300
),
Site = levels(quadrat_model_data$Site)[1]
)
lmer_preds <- predict(
q_lmer_quad,
newdata = cogon_seq,
re.form = NA
)
# 95% CI
set.seed(97)
boot_preds <- bootMer(
q_lmer_quad,
FUN = function(m) predict(m, newdata = cogon_seq, re.form = NA),
nsim = 500,
use.u = FALSE,
type = "parametric"
)
cogon_seq <- cogon_seq %>%
mutate(
fit = lmer_preds,
lwr = apply(boot_preds$t, 2, quantile, probs = 0.025),
upr = apply(boot_preds$t, 2, quantile, probs = 0.975)
)
# Thinned/jittered points
set.seed(97)
thinned_pts <- quadrat_model_data %>%
mutate(bin = cut(Cogongrass_Cover_Quadrat, breaks = 40, include.lowest = TRUE)) %>%
group_by(bin) %>%
slice_sample(n = 10, replace = FALSE) %>%
ungroup()
# Color palette
col_point <- "#8B6F57"
col_ribbon <- "#C9AC87"
col_line <- "#3E2210"
col_text <- "#1E1208"
# Figure
fig_lmer <- ggplot() +
geom_jitter(
data = thinned_pts,
aes(x = Cogongrass_Cover_Quadrat, y = Shannon_Diversity),
colour = col_point,
alpha = 0.45,
size = 1.3,
shape = 16,
width = 0.6,
height = 0
) +
geom_ribbon(
data = cogon_seq,
aes(x = Cogongrass_Cover_Quadrat, ymin = lwr, ymax = upr),
fill = col_ribbon,
alpha = 0.5
) +
geom_line(
data = cogon_seq,
aes(x = Cogongrass_Cover_Quadrat, y = fit),
colour = col_line,
linewidth = 1.0
) +
labs(
x = expression("Cogongrass" ~ "cover (% per quadrat)"),
y = expression("Shannon diversity" ~ (italic(H)*"'"))
) +
theme_classic(base_size = 12) +
theme(
axis.title = element_text(colour = col_text, size = 12),
axis.text = element_text(colour = col_text, size = 10),
axis.line = element_line(colour = col_text, linewidth = 0.4),
axis.ticks = element_line(colour = col_text, linewidth = 0.4),
axis.ticks.length = unit(2.5, "mm"),
panel.grid = element_blank(),
panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(12, 16, 10, 12)
) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
scale_x_continuous(
expand = expansion(mult = c(0.01, 0.02))
)
print(fig_lmer)
# Save
ggsave(
filename = "figure_lmer_quadrat_shannon.tiff",
plot = fig_lmer,
device = "tiff",
width = 3.5,
height = 3.2,
units = "in",
dpi = 600,
compression = "lzw"
)
ggsave(
filename = "figure_lmer_quadrat_shannon.pdf",
plot = fig_lmer,
device = "pdf",
width = 3.5,
height = 3.2,
units = "in"
)
message("Saved: figure_lmer_quadrat_shannon.tiff / .pdf")## Saved: figure_lmer_quadrat_shannon.tiff / .pdfGAM Plot
# Smooth predictions
cogon_seq <- data.frame(
Cogongrass_Cover_Quadrat = seq(
min(quadrat_model_data$Cogongrass_Cover_Quadrat, na.rm = TRUE),
max(quadrat_model_data$Cogongrass_Cover_Quadrat, na.rm = TRUE),
length.out = 300
),
Site = levels(quadrat_model_data$Site)[1]
)
preds <- predict(
q_gam,
newdata = cogon_seq,
exclude = "s(Site)",
se.fit = TRUE,
type = "response"
)
cogon_seq <- cogon_seq %>%
mutate(
fit = preds$fit,
se = preds$se.fit,
lwr = fit - 1.96 * se,
upr = fit + 1.96 * se
)
# Thin points
set.seed(42)
thinned_pts <- quadrat_model_data %>%
mutate(bin = cut(Cogongrass_Cover_Quadrat, breaks = 40, include.lowest = TRUE)) %>%
group_by(bin) %>%
slice_sample(n = 10, replace = FALSE) %>%
ungroup()
# Color palette
col_point <- "#8B6F57"
col_ribbon <- "#C9AC87"
col_line <- "#3E2210"
col_text <- "#1E1208"
fig_gam <- ggplot() +
geom_jitter(
data = thinned_pts,
aes(x = Cogongrass_Cover_Quadrat, y = Shannon_Diversity),
colour = col_point,
alpha = 0.45,
size = 1.3,
shape = 16,
width = 0.6,
height = 0
) +
# 95% CI ribbon
geom_ribbon(
data = cogon_seq,
aes(x = Cogongrass_Cover_Quadrat, ymin = lwr, ymax = upr),
fill = col_ribbon,
alpha = 0.5
) +
# GAM smooth line
geom_line(
data = cogon_seq,
aes(x = Cogongrass_Cover_Quadrat, y = fit),
colour = col_line,
linewidth = 1.0
) +
labs(
x = expression("Cogongrass" ~ "cover (% per quadrat)"),
y = expression("Shannon diversity" ~ (italic(H)*"'"))
) +
theme_classic(base_size = 12) +
theme(
axis.title = element_text(colour = col_text, size = 12),
axis.text = element_text(colour = col_text, size = 10),
axis.line = element_line(colour = col_text, linewidth = 0.4),
axis.ticks = element_line(colour = col_text, linewidth = 0.4),
axis.ticks.length = unit(2.5, "mm"),
panel.grid = element_blank(),
panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(12, 16, 10, 12)
) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
scale_x_continuous(
expand = expansion(mult = c(0.01, 0.02))
)
print(fig_gam)
# Save
ggsave(
filename = "figure_gam_quadrat_shannon.tiff",
plot = fig_gam,
device = "tiff",
width = 3.5,
height = 3.2,
units = "in",
dpi = 600,
compression = "lzw"
)
ggsave(
filename = "figure_gam_quadrat_shannon.pdf",
plot = fig_gam,
device = "pdf",
width = 3.5,
height = 3.2,
units = "in"
)PLOT-LEVEL ANALYSIS (retained from original)
Shrub cover conversion
shrub_cover <- shrub_data %>%
mutate(Cover = Line_End - Line_Start) %>%
group_by(Species_Name, Plot, InvStatus) %>%
summarise(Total_Cover = sum(Cover, na.rm = TRUE), .groups = "drop") %>%
mutate(Percent_Cover = Total_Cover / 3000 * 100)Herbaceous cover conversion (plot level)
veg_cover_summed <- Veg_Cover %>%
group_by(Plot, Species_Name, InvStatus) %>%
summarize(total_cover = sum(Cover_Per, na.rm = TRUE), .groups = "drop")
avg_species_cover <- veg_cover_summed %>%
left_join(quadrat_count, by = "Plot") %>%
mutate(avg_cover = total_cover / total_quadrats)Merge herb + shrub cover (plot level)
combined_cover <- avg_species_cover %>%
full_join(
shrub_cover %>% select(Plot, Species_Name, Percent_Cover, InvStatus),
by = c("Plot", "Species_Name")
) %>%
mutate(
overlap_flag = ifelse(!is.na(avg_cover) & !is.na(Percent_Cover), TRUE, FALSE),
final_cover = case_when(
!is.na(avg_cover) & is.na(Percent_Cover) ~ avg_cover,
is.na(avg_cover) & !is.na(Percent_Cover) ~ Percent_Cover,
TRUE ~ NA_real_
)
)Extract plot-level cogongrass cover
cogongrass_cover <- combined_cover %>%
filter(Species_Name %in% c("Imperata_cylindrica")) %>%
group_by(Plot) %>%
summarise(Cogongrass_Cover = sum(final_cover, na.rm = TRUE), .groups = "drop")Plot-level species matrix and Shannon diversity
species_matrix_plot <- combined_cover %>%
filter(
!Species_Name %in% c("Imperata_cylindrica", "Imperata_cylindrica_Live"),
InvStatus.x == "Native"
) %>%
select(Plot, Species_Name, final_cover) %>%
pivot_wider(names_from = Species_Name, values_from = final_cover, values_fill = 0)
shannon_plot <- species_matrix_plot %>%
select(-Plot) %>%
vegan::diversity(index = "shannon") %>%
as.data.frame() %>%
setNames("Shannon_Diversity") %>%
mutate(Plot = species_matrix_plot$Plot)Merge plot-level data
model_data <- shannon_plot %>%
left_join(cogongrass_cover, by = "Plot") %>%
left_join(EnvironmentalOutputs, by = "Plot") %>%
mutate(Cogongrass_Cover = ifelse(is.na(Cogongrass_Cover), 0, Cogongrass_Cover))
model_data$Site <- as.factor(model_data$Site)
summary(model_data)## Shannon_Diversity Plot Cogongrass_Cover Site
## Min. :0.000 Length:206 Min. : 0.00 WSF :44
## 1st Qu.:1.793 Class :character 1st Qu.: 0.00 Jay :42
## Median :2.201 Mode :character Median : 0.00 BCNWR :32
## Mean :2.153 Mean :17.08 DSNF :24
## 3rd Qu.:2.605 3rd Qu.:29.82 CNF :23
## Max. :3.360 Max. :92.14 Escambia:11
## (Other) :30
## BurnYear Camera Latitude Longitude
## Min. :2000 Length:206 Min. :28.48 Min. :-89.80
## 1st Qu.:2014 Class :character 1st Qu.:29.56 1st Qu.:-88.91
## Median :2022 Mode :character Median :30.77 Median :-87.14
## Mean :2018 Mean :30.38 Mean :-86.72
## 3rd Qu.:2022 3rd Qu.:30.98 3rd Qu.:-83.58
## Max. :2024 Max. :31.59 Max. :-81.94
## NA's :42
## Muname NLCD_LandCover Region SoilType
## Length:206 Length:206 Length:206 Length:206
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Status TreeCanopy_NLCD aspect elevation
## Length:206 Min. : 0.00 Min. : 0.0 Min. : 17.00
## Class :character 1st Qu.:51.50 1st Qu.: 90.0 1st Qu.: 38.00
## Mode :character Median :75.50 Median :180.0 Median : 59.00
## Mean :66.22 Mean :171.7 Mean : 57.68
## 3rd Qu.:86.75 3rd Qu.:270.0 3rd Qu.: 75.00
## Max. :98.00 Max. :350.6 Max. :111.00
##
## ndvi pr slope sph
## Min. :0.2728 Min. :3.850 Min. : 0.000 Min. :0.01051
## 1st Qu.:0.4303 1st Qu.:4.345 1st Qu.: 2.149 1st Qu.:0.01074
## Median :0.4624 Median :4.592 Median : 2.984 Median :0.01091
## Mean :0.4503 Mean :4.566 Mean : 3.601 Mean :0.01145
## 3rd Qu.:0.4846 3rd Qu.:4.796 3rd Qu.: 4.688 3rd Qu.:0.01183
## Max. :0.5357 Max. :5.294 Max. :18.429 Max. :0.01325
##
## srad tmmn agb
## Min. :204.7 Min. :285.7 Min. : 5.00
## 1st Qu.:206.8 1st Qu.:286.6 1st Qu.: 74.00
## Median :207.0 Median :287.2 Median : 96.00
## Mean :210.6 Mean :287.4 Mean : 99.92
## 3rd Qu.:216.6 3rd Qu.:288.6 3rd Qu.:133.00
## Max. :224.2 Max. :290.3 Max. :211.00
## Plot-level model fitting
This is a part of Alan’s main biodiversity paper.
shannon_lmer_ml <- lmer(Shannon_Diversity ~ Cogongrass_Cover + (1 | Site),
data = model_data, REML = FALSE)
shannon_lmer_quad_ml <- lmer(Shannon_Diversity ~ Cogongrass_Cover + I(Cogongrass_Cover^2) +
(1 | Site), data = model_data, REML = FALSE)## Warning: Some predictor variables are on very different scales: consider
## rescalingshannon_gam_ml <- gam(Shannon_Diversity ~ s(Cogongrass_Cover) + s(Site, bs = "re"),
data = model_data, method = "ML")BIOTIC RESISTANCE / IMPACT SENSITIVITY ANALYSIS
Two complementary questions:
Part 1 — Biotic Resistance: Does plot-level native diversity predict cogongrass cover?
(diversity → invasion intensity)
Part 2 — Impact Sensitivity: Does pre-invasion diversity moderate the slope of
cogongrass impact on quadrat-level Shannon diversity?
(diversity × cogon interaction → quadrat Shannon)
Part 1: Biotic Resistance
Does plot-level native diversity predict cogongrass cover?
# model_data already has plot-level Shannon and Cogongrass_Cover
# Here Shannon_Diversity is the predictor; Cogongrass_Cover is the response
biotic_resistance_lmer <- lmer(
Cogongrass_Cover ~ Shannon_Diversity + (1 | Site),
data = model_data, REML = FALSE
)
biotic_resistance_quad <- lmer(
Cogongrass_Cover ~ Shannon_Diversity + I(Shannon_Diversity^2) + (1 | Site),
data = model_data, REML = FALSE
)
biotic_resistance_gam <- gam(
Cogongrass_Cover ~ s(Shannon_Diversity) + s(Site, bs = "re"),
data = model_data, method = "ML"
)
AIC(biotic_resistance_lmer, biotic_resistance_quad)## df AIC
## biotic_resistance_lmer 4 1901.693
## biotic_resistance_quad 5 1898.681summary(biotic_resistance_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Cogongrass_Cover ~ s(Shannon_Diversity) + s(Site, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.772 2.468 6.797 1.21e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Shannon_Diversity) 1.768 2.246 3.667 0.0183 *
## s(Site) 3.819 9.000 1.056 0.0148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0771 Deviance explained = 10.2%
## -ML = 946.88 Scale est. = 552.99 n = 206