Quadrat Inside Outside
Alan Ivory
2026-04-07
pacman::p_load(tidyverse, readxl, vegan, mgcv, lme4,
marginaleffects, gratia, lmerTest,
ggdist, emmeans)
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$Cogon_Cover_Class <- factor(
quadrat_model_data$Cogongrass_Cover_Quadrat,
levels = c(0, 2.5, 15, 37.5, 62.5, 85, 97.5)
)
quadrat_model_data$Cogon_Cover_Class <- relevel(
quadrat_model_data$Cogon_Cover_Class,
ref = "0"
)
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
##
## Status SoilType Soil_Group TreeCanopy_NLCD
## Length:1407 Length:1407 Min. : 18.0 Min. : 0.00
## Class :character Class :character 1st Qu.:377.0 1st Qu.:53.00
## Mode :character Mode :character Median :381.0 Median :76.00
## Mean :330.4 Mean :66.54
## 3rd Qu.:389.0 3rd Qu.:87.00
## Max. :389.0 Max. :98.00
##
## aspect elevation ndvi pr
## Min. : 0.0 Min. : 17.00 Min. :0.2728 Min. :3.850
## 1st Qu.: 90.0 1st Qu.: 38.00 1st Qu.:0.4320 1st Qu.:4.286
## Median :180.0 Median : 59.00 Median :0.4636 Median :4.592
## Mean :171.0 Mean : 57.71 Mean :0.4507 Mean :4.566
## 3rd Qu.:270.0 3rd Qu.: 75.00 3rd Qu.:0.4847 3rd Qu.:4.796
## Max. :350.6 Max. :111.00 Max. :0.5357 Max. :5.294
##
## slope sph srad tmmn
## Min. : 0.000 Min. :0.01051 Min. :204.7 Min. :285.7
## 1st Qu.: 2.147 1st Qu.:0.01074 1st Qu.:206.8 1st Qu.:286.6
## Median : 2.983 Median :0.01091 Median :207.0 Median :287.2
## Mean : 3.584 Mean :0.01145 Mean :210.5 Mean :287.4
## 3rd Qu.: 4.684 3rd Qu.:0.01194 3rd Qu.:218.5 3rd Qu.:289.0
## Max. :18.429 Max. :0.01325 Max. :224.2 Max. :290.3
##
## agb Cogon_Cover_Class
## Min. : 3.45 0 :905
## 1st Qu.:10.60 2.5 :107
## Median :17.91 15 : 98
## Mean :19.17 37.5: 78
## 3rd Qu.:25.44 62.5: 61
## Max. :56.03 85 : 71
## NA's :208 97.5: 87Presence/Absence Option
Binary presence/absence of cogongrass
quadrat_model_data <- quadrat_model_data %>%
mutate(Cogon_Present = ifelse(Cogongrass_Cover_Quadrat > 0, 1, 0))Identify Mixed Plots
# A qualifying plot must have at least one invaded AND one uninvaded quadrat
mixed_plots <- quadrat_model_data %>%
group_by(Plot) %>%
summarise(
n_invaded = sum(Cogon_Present),
n_uninvaded = sum(1 - Cogon_Present),
.groups = "drop"
) %>%
filter(n_invaded > 0, n_uninvaded > 0) %>%
pull(Plot)
cat("Number of qualifying plots:", length(mixed_plots), "\n")## Number of qualifying plots: 69subset_data <- quadrat_model_data %>%
filter(Plot %in% mixed_plots) %>%
mutate(Cogon_Present = factor(Cogon_Present, levels = c(0, 1),
labels = c("Absent", "Present")))
cat("Quadrats in analysis:", nrow(subset_data), "\n")## Quadrats in analysis: 476cat("Plots in analysis: ", n_distinct(subset_data$Plot), "\n")## Plots in analysis: 69Fit linear mixed model with random intercept for Plot
model <- lmer(Shannon_Diversity ~ Cogon_Present + (1 | Plot) + (1 | Site),
data = subset_data,
REML = TRUE)Results
summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon_Diversity ~ Cogon_Present + (1 | Plot) + (1 | Site)
## Data: subset_data
##
## REML criterion at convergence: 709.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95513 -0.65177 0.04372 0.62746 2.68291
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.06096 0.2469
## Site (Intercept) 0.08471 0.2910
## Residual 0.21127 0.4596
## Number of obs: 476, groups: Plot, 69; Site, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.43644 0.10740 8.61779 13.374 4.64e-07 ***
## Cogon_PresentPresent -0.29733 0.04585 446.48812 -6.485 2.36e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Cgn_PrsntPr -0.255# Fixed effect estimate and 95% CI
confint(model, parm = "Cogon_PresentPresent", method = "Wald")## 2.5 % 97.5 %
## Cogon_PresentPresent -0.3871922 -0.2074695# Random effect variance
VarCorr(model)## Groups Name Std.Dev.
## Plot (Intercept) 0.24689
## Site (Intercept) 0.29105
## Residual 0.45964# Proportion of variance explained by plot
icc_plot <- as.numeric(VarCorr(model)$Plot) /
(as.numeric(VarCorr(model)$Plot) + sigma(model)^2)
cat("ICC (plot-level):", round(icc_plot, 3), "\n")## ICC (plot-level): 0.224icc_site <- as.numeric(VarCorr(model)$Site) /
(as.numeric(VarCorr(model)$Site) + sigma(model)^2)
cat("ICC (site-level):", round(icc_site, 3), "\n")## ICC (site-level): 0.286Plot
plot_means <- subset_data %>%
group_by(Plot, Cogon_Present) %>%
summarise(mean_H = mean(Shannon_Diversity), .groups = "drop")
ggplot(plot_means, aes(x = Cogon_Present, y = mean_H, group = Plot)) +
geom_line(alpha = 0.35, colour = "grey50") +
geom_point(aes(colour = Cogon_Present), size = 2.5, alpha = 0.8) +
scale_colour_manual(values = c("Absent" = "#1D9E75", "Present" = "#D85A30")) +
labs(x = "Cogongrass", y = "Mean Shannon diversity (H')",
title = "Within-plot diversity contrast",
subtitle = paste0("n = ", length(mixed_plots), " mixed plots"),
colour = NULL) +
theme_bw(base_size = 13) +
theme(legend.position = "top")
Alternative Plot
# Marginal means
emm <- emmeans(model, ~ Cogon_Present)
emm_df <- as.data.frame(emm)
ggplot() +
geom_jitter(data = subset_data,
aes(x = Cogon_Present, y = Shannon_Diversity, colour = Cogon_Present),
width = 0.15, alpha = 0.25, size = 1.5) +
ggdist::stat_halfeye(data = subset_data,
aes(x = Cogon_Present, y = Shannon_Diversity, fill = Cogon_Present),
adjust = 0.8, width = 0.5, .width = 0,
point_colour = NA, alpha = 0.4, side = "left") +
geom_pointrange(data = emm_df,
aes(x = Cogon_Present, y = emmean,
ymin = lower.CL, ymax = upper.CL),
size = 0.9, linewidth = 1.2, colour = "black") +
scale_colour_manual(values = c("Absent" = "#1D9E75", "Present" = "#D85A30")) +
scale_fill_manual( values = c("Absent" = "#1D9E75", "Present" = "#D85A30")) +
labs(x = "Cogongrass", y = "Shannon diversity (H')",
title = "Effect of cogongrass presence on understory plant diversity") +
theme_bw(base_size = 13) +
theme(legend.position = "none")
Cogongrass Cover Option
Data Cleaning
subset_data <- subset_data %>%
mutate(
Cover_cat = factor(Cogongrass_Cover_Quadrat,
levels = c(0, 2.5, 15, 37.5, 62.5, 85, 97.5),
labels = c("0", "2.5", "15", "37.5", "62.5", "85", "97.5"),
ordered = TRUE)
)Cover category model
# Treatment contrasts: 0% cover is the baseline
contrasts(subset_data$Cover_cat) <- contr.treatment(7)
model_cover_cat <- lmer(Shannon_Diversity ~ Cover_cat + (1 | Plot) + (1 | Site),
data = subset_data,
REML = TRUE)
summary(model_cover_cat)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon_Diversity ~ Cover_cat + (1 | Plot) + (1 | Site)
## Data: subset_data
##
## REML criterion at convergence: 672.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.10013 -0.64888 0.04355 0.62442 2.79030
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.05539 0.2354
## Site (Intercept) 0.11061 0.3326
## Residual 0.18984 0.4357
## Number of obs: 476, groups: Plot, 69; Site, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.42242 0.11800 8.64401 12.054 1.06e-06 ***
## Cover_cat2 -0.12970 0.06242 433.72052 -2.078 0.03832 *
## Cover_cat3 -0.19634 0.07013 435.42784 -2.800 0.00534 **
## Cover_cat4 -0.16029 0.07811 439.29768 -2.052 0.04076 *
## Cover_cat5 -0.37105 0.08804 443.03431 -4.214 3.04e-05 ***
## Cover_cat6 -0.42236 0.08596 449.31657 -4.913 1.26e-06 ***
## Cover_cat7 -0.76503 0.08333 459.52515 -9.181 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Cvr_c2 Cvr_c3 Cvr_c4 Cvr_c5 Cvr_c6
## Cover_cat2 -0.145
## Cover_cat3 -0.145 0.275
## Cover_cat4 -0.125 0.230 0.198
## Cover_cat5 -0.119 0.209 0.193 0.172
## Cover_cat6 -0.111 0.190 0.178 0.184 0.164
## Cover_cat7 -0.104 0.199 0.173 0.167 0.142 0.185# Marginal means
emm_cover <- emmeans(model_cover_cat, ~ Cover_cat)
emm_cover_df <- as.data.frame(emm_cover)
# Pairwise contrasts vs. 0% cover class
contrast(emm_cover, method = "trt.vs.ctrl", ref = 1)## contrast estimate SE df t.ratio p.value
## Cover_cat2.5 - Cover_cat0 -0.130 0.0625 434 -2.074 0.1713
## Cover_cat15 - Cover_cat0 -0.196 0.0703 436 -2.794 0.0285
## Cover_cat37.5 - Cover_cat0 -0.160 0.0783 439 -2.048 0.1808
## Cover_cat62.5 - Cover_cat0 -0.371 0.0883 443 -4.204 0.0002
## Cover_cat85 - Cover_cat0 -0.422 0.0862 450 -4.899 <.0001
## Cover_cat97.5 - Cover_cat0 -0.765 0.0837 460 -9.143 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: dunnettx method for 6 tests# Model comparison: presence/absence vs. cover category
model_pres_ML <- update(model, REML = FALSE)
model_cover_cat_ML <- update(model_cover_cat, REML = FALSE)
AIC(model_pres_ML, model_cover_cat_ML)## df AIC
## model_pres_ML 5 712.0281
## model_cover_cat_ML 10 669.5590ggplot() +
geom_jitter(data = subset_data,
aes(x = Cover_cat, y = Shannon_Diversity),
width = 0.15, alpha = 0.2, size = 1.5, colour = "#888") +
ggdist::stat_halfeye(data = subset_data,
aes(x = Cover_cat, y = Shannon_Diversity),
adjust = 0.8, width = 0.5, .width = 0,
fill = "#D85A30", alpha = 0.3, side = "left",
point_colour = NA) +
geom_pointrange(data = emm_cover_df,
aes(x = Cover_cat, y = emmean,
ymin = lower.CL, ymax = upper.CL),
size = 0.9, linewidth = 1.2, colour = "black") +
geom_line(data = emm_cover_df,
aes(x = Cover_cat, y = emmean, group = 1),
colour = "black", linewidth = 0.7, linetype = "dashed") +
labs(x = "Cogongrass cover class (%)",
y = "Shannon diversity (H')",
title = "Understory diversity across cogongrass cover classes") +
theme_bw(base_size = 13) +
theme(legend.position = "none")