Plot-Level Analysis: Invasive Species Effects

Evaluating Impacts on Native Flora Biodiversity and Composition

Introduction

This document performs a plot-level analysis examining the effects of invasive species cover on native flora. The analysis aggregates data collected across 12 quadrats per plot by averaging cover values, resulting in a dataset representing 100 distinct plots.

1. Environment Setup & Libraries

We begin by loading all necessary packages required for data manipulation, diversity metrics, regression modeling, and ecological community analysis.

library(dplyr)
library(tidyr)
library(readxl)
library(vegan)
library(hillR)
library(segmented)
library(indicspecies)
library(ggplot2)
library(corrplot)
library(betapart)
library(lme4)
library(car)
library(emmeans)
library(scatterplot3d)
library(cowplot)
library(knitr)
library(viridis)
library(gt)

2. Data Loading & Preprocessing

Load Dataset

First, we load the raw quadrat-level dataset.

Note

Ensure your invasive_cover.xlsx file is placed directly inside your active R project directory to support reproducible parsing paths.

setwd("~/Google Drive/nduah_analyses")

invasive_cover <- read_excel("invasive_cover.xlsx")

Format Transformation & Mid-Point Conversion

We convert the ordinal cover rankings to numeric mid-point percentages using a lookup vector. Then, we transform the community data into a long format for easier aggregation.

ordinal_to_midpoint <- c(
  "1" = 0.5,    # <1%
  "2" = 2.5,    # 1‑5%
  "3" = 7.5,    # 5‑10%
  "4" = 17.5,   # 10‑25%
  "5" = 32.5,   # 25‑50%
  "6" = 62.5,   # 50‑75%
  "7" = 87.5    # 75‑100%
)

quadrat_cols <- colnames(invasive_cover)[5:ncol(invasive_cover)]

df_long <- invasive_cover %>%
  pivot_longer(cols = all_of(quadrat_cols),
               names_to = "quadrat_id",
               values_to = "cover_ord") %>%
  mutate(
    cover_pct = ifelse(cover_ord == 0, 0, ordinal_to_midpoint[as.character(cover_ord)]),
    plot_id = sub("Q.*$", "", quadrat_id)
  )

head(df_long)

# A tibble: 6 × 8
  species_name         type  form  family quadrat_id cover_ord cover_pct plot_id
  <chr>                <chr> <chr> <chr>  <chr>          <dbl>     <dbl> <chr>  
1 Abutilon mauritianum Nati… Shrub Malva… A10Q1              0         0 A10    
2 Abutilon mauritianum Nati… Shrub Malva… A10Q10             0         0 A10    
3 Abutilon mauritianum Nati… Shrub Malva… A10Q11             0         0 A10    
4 Abutilon mauritianum Nati… Shrub Malva… A10Q12             0         0 A10    
5 Abutilon mauritianum Nati… Shrub Malva… A10Q2              0         0 A10    
6 Abutilon mauritianum Nati… Shrub Malva… A10Q3              0         0 A10    

Data Aggregation

We average quadrat cover levels up to the Plot \(\times\) Species level.

plot_species <- df_long %>%
  group_by(plot_id, species_name, type, form, family) %>%
  summarise(cover_pct = mean(cover_pct, na.rm = TRUE), .groups = "drop")

head(plot_species)

# A tibble: 6 × 6
  plot_id species_name            type     form      family        cover_pct
  <chr>   <chr>                   <chr>    <chr>     <chr>             <dbl>
1 A1      Abutilon mauritianum    Native   Shrub     Malvaceae         0    
2 A1      Acanthospermum hispidum Invasive Forb/herb Asteraceae        0    
3 A1      Achyranthes aspera      Invasive Forb/herb Amaranthaceae     0    
4 A1      Ageratum conyzoides     Invasive Forb/herb Asteraceae        0.417
5 A1      Ajuga integrifolia      Native   Forb/herb Lamiaceae         0    
6 A1      Alysicarpus glumaceus   Native   Forb/herb Fabaceae          0    

3. Species per life form

# ------------------------------------------------------------
# UNIQUE SPECIES PER LIFE FORM × ORIGIN
# (count distinct species names, not plot occurrences)
# ------------------------------------------------------------
lifeform_richness <- plot_species %>%
  mutate(origin = ifelse(type == "Native", "Native", "Exotic")) %>%
  group_by(form, origin) %>%
  summarise(n_species = n_distinct(species_name), .groups = "drop") %>%
  pivot_wider(names_from = origin, values_from = n_species, values_fill = 0) %>%
  dplyr::select(form, Native, Exotic) %>%          # <-- guarantees column order
  mutate(Total = Native + Exotic,
         Exotic_pct = round(100 * Exotic / Total, 1)) %>%
  arrange(desc(Total))

lifeform_richness %>%
  gt() %>%
  tab_caption("Species richness of native and exotic flora by life form")

Species richness of native and exotic flora by life form

form	Native	Exotic	Total	Exotic_pct
Forb/herb	17	16	33	48.5
Grass	26	2	28	7.1
Shrub	11	3	14	21.4
Tree	2	1	3	33.3

# ------------------------------------------------------------
# MEAN COVER PER LIFE FORM × ORIGIN (averaged across plots)
# ------------------------------------------------------------
lifeform_cover <- plot_species %>%
  mutate(origin = ifelse(type == "Native", "Native", "Exotic")) %>%
  group_by(form, origin) %>%
  summarise(
    mean_cover = round(mean(cover_pct, na.rm = TRUE), 2),
    sd_cover   = round(sd(cover_pct,   na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(cover_label = paste0(mean_cover, " ± ", sd_cover)) %>%
  dplyr::select(form, origin, cover_label) %>%
  pivot_wider(names_from = origin, values_from = cover_label,
              values_fill = "—") %>%
  # Rename AND reorder before kable touches it — no col.names argument needed
  dplyr::transmute(
    `Life Form`              = form,
    `Native (mean ± SD %)`   = Native,
    `Exotic (mean ± SD %)`   = Exotic
  ) %>%
  arrange(`Life Form`)
lifeform_cover %>%
  gt() %>%
  tab_caption("Mean cover (%) of native and exotic flora by life form (across all plot × species records)")

Mean cover (%) of native and exotic flora by life form (across all plot × species records)

Life Form	Native (mean ± SD %)	Exotic (mean ± SD %)
Forb/herb	0.81 ± 1.88	1.33 ± 2.43
Grass	2.19 ± 9.34	0.15 ± 0.39
Shrub	0.75 ± 2.84	2.06 ± 3.85
Tree	0.04 ± 0.32	0.01 ± 0.05

Full table of species per plot

This table tells us how frequent each species of native and exotic species were and help to show who is who in this grassland community.

species_frequency <- plot_species %>%
  filter(cover_pct > 0) %>%          # presence only for frequency count
  group_by(species_name, type, form, family) %>%
  summarise(
    n_plots      = n_distinct(plot_id),
    mean_cover   = round(mean(cover_pct, na.rm = TRUE), 2),
    sd_cover     = round(sd(cover_pct,   na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(
    cover_label  = paste0(mean_cover, " ± ", sd_cover),
    freq_pct     = round(100 * n_plots / n_distinct(plot_species$plot_id), 1),
    origin       = ifelse(type == "Native", "Native", "Exotic")
  ) %>%
  dplyr::transmute(
    `Species`         = species_name,
    `Family`          = family,
    `Life Form`       = form,
    `Origin`          = origin,
    `N Plots`         = n_plots,
    `Frequency (%)`   = freq_pct,
    `Mean Cover ± SD` = cover_label
  ) %>%
  arrange(desc(`N Plots`), `Species`)

species_frequency %>%
  gt() %>%
  tab_caption("All species ordered by plot frequency (n = 100 plots)") %>%
  tab_style(
    style = cell_text(style = "italic"),
    locations = cells_body(columns = `Species`)
  ) %>%
  tab_style(
    style = cell_fill(color = "#fff3e0"),
    locations = cells_body(rows = `Origin` == "Exotic")
  ) %>%
  data_color(
    columns = `Frequency (%)`,
    method  = "numeric",
    palette = "Blues"
  )

All species ordered by plot frequency (n = 100 plots)

Species	Family	Life Form	Origin	N Plots	Frequency (%)	Mean Cover ± SD
Cynodon dactylon	Poaceae	Grass	Native	100	100	45.85 ± 14.15
Xanthium orientale	Asteraceae	Forb/herb	Exotic	100	100	4.74 ± 2.7
Cyperus odoratus	Cyperaceae	Grass	Native	99	99	2.89 ± 2.53
Dyschoriste radicans	Acanthaceae	Shrub	Native	99	99	7.45 ± 6.24
Solanum incanum	Solanaceae	Shrub	Exotic	99	99	6.21 ± 4.39
Bidens pilosa	Asteraceae	Forb/herb	Exotic	98	98	2.11 ± 1.62
Trifolium semipilosum	Fabaceae	Forb/herb	Native	98	98	4.9 ± 4
Eragrostis tenuifolia	Poaceae	Grass	Native	97	97	6.38 ± 5.79
Dichondra repens	Convolvulaceae	Forb/herb	Exotic	95	95	6.81 ± 4.72
Indigofera colutea	Fabaceae	Forb/herb	Native	93	93	2.58 ± 2.16
Euphorbia heterophylla	Euphorbiaceae	Forb/herb	Exotic	89	89	1.49 ± 1
Schkuhria pinnata	Asteraceae	Forb/herb	Exotic	83	83	1.12 ± 1.04
Cirsium vulgare	Asteraceae	Forb/herb	Exotic	78	78	1.53 ± 1.46
Tagetes minuta	Asteraceae	Forb/herb	Exotic	77	77	1.41 ± 1.07
Ageratum conyzoides	Asteraceae	Forb/herb	Exotic	75	75	0.72 ± 0.68
Cycnium tubulosum	Orobanchaceae	Forb/herb	Native	75	75	1.44 ± 1.48
Glycine wightii	Fabaceae	Forb/herb	Native	70	70	2.45 ± 2.54
Conyza stricta	Asteraceae	Forb/herb	Native	67	67	0.9 ± 0.87
Acanthospermum hispidum	Asteraceae	Forb/herb	Exotic	66	66	1.28 ± 1.14
Achyranthes aspera	Amaranthaceae	Forb/herb	Exotic	64	64	2.3 ± 2.37
Sida cuneifolia	Malvaceae	Forb/herb	Native	63	63	1.13 ± 1.18
Panicum massaiense	Poaceae	Grass	Native	54	54	1.18 ± 2.2
Teramnus uncinatus	Fabaceae	Forb/herb	Native	51	51	1.93 ± 2.55
Digitaria ciliaris	Poaceae	Grass	Exotic	50	50	0.57 ± 0.61
Sida rhombifolia	Malvaceae	Shrub	Native	50	50	1.33 ± 1.03
Commelina benghalensis	Commelinaceae	Forb/herb	Native	47	47	0.67 ± 0.81
Cucumis aculeatus	Cucurbitaceae	Forb/herb	Native	38	38	0.72 ± 0.64
Emilia discifolia	Asteraceae	Forb/herb	Native	38	38	0.72 ± 0.63
Sporobolus pyramidalis	Poaceae	Grass	Native	37	37	0.97 ± 1.04
Leonotis nepetifolia	Lamiaceae	Forb/herb	Native	36	36	0.44 ± 0.43
Brachiaria eruciformis	Poaceae	Grass	Native	35	35	1.25 ± 1.43
Convolvulus sagittatus	Convolvulaceae	Forb/herb	Exotic	33	33	0.73 ± 0.63
Hypoestes verticillaris	Acanthaceae	Forb/herb	Native	32	32	0.56 ± 0.47
Eragrostis exasperata	Poaceae	Grass	Native	22	22	1.15 ± 2.53
Erigeron canadensis	Asteraceae	Forb/herb	Exotic	21	21	0.41 ± 0.47
Gutenbergia boranensis	Asteraceae	Forb/herb	Native	21	21	0.39 ± 0.48
Oxalis stricta	Oxalidaceae	Forb/herb	Exotic	21	21	0.41 ± 0.44
Gutenbergia cordifolia	Asteraceae	Forb/herb	Native	18	18	0.66 ± 0.72
Ocimum gratissimum	Lamiaceae	Shrub	Native	17	17	0.43 ± 0.48
Panicum maximum	Poaceae	Grass	Native	15	15	0.65 ± 0.61
Silybum marianum	Asteraceae	Forb/herb	Exotic	14	14	0.44 ± 0.36
Sorghum versicolor	Poaceae	Grass	Native	14	14	0.21 ± 0.14
Cynoglossum amabile	Boraginaceae	Forb/herb	Exotic	11	11	0.53 ± 0.38
Asclepias fruticosa	Apocynaceae	Shrub	Native	9	9	0.24 ± 0.18
Lippia kituiensis	Verbenaceae	Shrub	Native	9	9	1.26 ± 1.12
Vachellia xanthophloea	Fabaceae	Tree	Native	9	9	0.8 ± 1.37
Chloris pycnothrix	Poaceae	Grass	Native	8	8	0.64 ± 0.57
Pennisetum mezianum	Poaceae	Grass	Native	8	8	0.73 ± 0.75
Bothriochloa insculpta	Poaceae	Grass	Native	7	7	0.39 ± 0.22
Ajuga integrifolia	Lamiaceae	Forb/herb	Native	6	6	0.59 ± 0.5
Nicoteba betonica	Acanthaceae	Shrub	Native	6	6	0.49 ± 0.22
Digitaria velutina	Poaceae	Grass	Exotic	4	4	0.14 ± 0.19
Brachiaria brizantha	Poaceae	Grass	Native	3	3	0.99 ± 1.17
Panicum coloratum	Poaceae	Grass	Native	3	3	0.83 ± 0.62
Abutilon mauritianum	Malvaceae	Shrub	Native	2	2	0.12 ± 0.12
Aristida kenyensis	Poaceae	Grass	Native	2	2	0.52 ± 0.15
Echinochloa colona	Poaceae	Grass	Native	2	2	0.33 ± 0.41
Indigofera spicata	Fabaceae	Forb/herb	Native	2	2	0.42 ± 0.29
Leucaena leucocephala	Fabaceae	Tree	Exotic	2	2	0.31 ± 0.15
Sida paniculata	Malvaceae	Shrub	Exotic	2	2	0.42 ± 0.29
Alysicarpus glumaceus	Fabaceae	Forb/herb	Native	1	1	1.46 ± NA
Aristida adoensis	Poaceae	Grass	Native	1	1	0.04 ± NA
Capparis tomentosa	Capparaceae	Shrub	Native	1	1	0.62 ± NA
Chloris gayana	Poaceae	Grass	Native	1	1	0.83 ± NA
Cynodon aethiopicus	Poaceae	Grass	Native	1	1	0.21 ± NA
Eragrostis namaquensis	Poaceae	Grass	Native	1	1	0.04 ± NA
Eragrostis rigidior	Poaceae	Grass	Native	1	1	0.04 ± NA
Fuerstia africana	Lamiaceae	Shrub	Native	1	1	0.04 ± NA
Harpachne schimperi	Poaceae	Grass	Native	1	1	0.62 ± NA
Leonotis mollissima	Lamiaceae	Shrub	Native	1	1	0.21 ± NA
Pogonarthria squarrosa	Poaceae	Grass	Native	1	1	0.21 ± NA
Psiadia punctulata	Asteraceae	Shrub	Native	1	1	0.21 ± NA
Senna didymobotrya	Fabaceae	Shrub	Exotic	1	1	1.67 ± NA
Setaria pumila	Poaceae	Grass	Native	1	1	0.21 ± NA
Solanum nigrum	Solanaceae	Forb/herb	Exotic	1	1	0.21 ± NA
Sporobolus filipes	Poaceae	Grass	Native	1	1	1.25 ± NA
Sporobolus fimbriatus	Poaceae	Grass	Native	1	1	0.62 ± NA
Vernonia amygdalina	Asteraceae	Tree	Native	1	1	0.04 ± NA

Interpretation: Many of the most widespread species are exotic ones, specially herbs/forbs… tye must have a really important role in the ecosystem fucntioning

3. Metric Calculations

Native and Exotic Diversity (Hill Numbers)

We restructure native and exotic species records into parallel wide format community matrices to calculate Hill numbers (\(q = 0, 1, 2\)).

native_wide <- plot_species %>%
  filter(type == "Native") %>%
  dplyr::select(plot_id, species_name, cover_pct) %>%
  pivot_wider(names_from = species_name, values_from = cover_pct, values_fill = 0) %>%
  as.data.frame()
rownames(native_wide) <- native_wide$plot_id
native_wide_matrix <- native_wide %>% dplyr::select(-plot_id)

exotic_wide <- plot_species %>%
  filter(type != "Native") %>%
  dplyr::select(plot_id, species_name, cover_pct) %>%
  pivot_wider(names_from = species_name, values_from = cover_pct, values_fill = 0) %>%
  as.data.frame()
rownames(exotic_wide) <- exotic_wide$plot_id
exotic_wide_matrix <- exotic_wide %>% dplyr::select(-plot_id)

native_div <- data.frame(
  plot_id = rownames(native_wide_matrix),
  q0 = hill_taxa(native_wide_matrix, q = 0),
  q1 = hill_taxa(native_wide_matrix, q = 1),
  q2 = hill_taxa(native_wide_matrix, q = 2)
)

exotic_div <- data.frame(
  plot_id = rownames(exotic_wide_matrix),
  q0 = hill_taxa(exotic_wide_matrix, q = 0),
  q1 = hill_taxa(exotic_wide_matrix, q = 1),
  q2 = hill_taxa(exotic_wide_matrix, q = 2)
)

4. Visualizations & Distributions

Histograms of True Diversity

combined_wide <- bind_rows(
  native_div %>% mutate(origin = "Native"),
  exotic_div %>% mutate(origin = "Exotic")
)

div_long <- combined_wide %>%
  pivot_longer(
    cols = c(q0, q1, q2),
    names_to = "q_level",
    values_to = "diversity_value"
  ) %>%
  mutate(
    origin = factor(origin, levels = c("Native", "Exotic")),
    q_level = case_when(
      q_level == "q0" ~ "q = 0 (Rare species)",
      q_level == "q1" ~ "q = 1 (Common species)",
      q_level == "q2" ~ "q = 2 (Dominant species)"
    )
  )

ggplot(div_long, aes(x = diversity_value, fill = origin)) +
  geom_histogram(bins = 10, color = "white", alpha = 0.85) +
  facet_grid(origin ~ q_level, scales = "free") +
  scale_fill_manual(values = c("Native" = "#2e7d32", "Exotic" = "#c62828")) +
  labs(
    title = "True Diversity Value Distributions",
    x = "Diversity Metric Value",
    y = "Number of Plots (Frequency)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    strip.text = element_text(face = "bold", size = 11),
    strip.background = element_rect(fill = "gray95", color = "gray80"),
    panel.spacing = unit(1.2, "lines"),
    plot.title = element_text(face = "bold", size = 13)
  ) -> hist.all

hist.all

Interpretation: Most plots with few species (left-skewed) means that highly diverse patches are uncommon and probably patchy.

Diversity Profiles

By summarizing mean tendencies alongside 95% Confidence Intervals

profile_stats <- div_long %>%
  group_by(origin, q_level) %>%
  summarise(
    n = n(),
    mean_val = mean(diversity_value, na.rm = TRUE),
    sd_val = sd(diversity_value, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    se_val = sd_val / sqrt(n),
    ci_lower = mean_val - (1.96 * se_val),
    ci_upper = mean_val + (1.96 * se_val)
  )

pd <- position_dodge(width = 0.25)
ggplot(profile_stats, aes(x = q_level, y = mean_val, color = origin, group = origin)) +
  geom_line(size = 1.2, position = pd) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.15, size = 0.8, position = pd) +
  geom_point(size = 4, stroke = 1.5, fill = "white", shape = 21, position = pd) +
  scale_color_manual(values = c("Native" = "#2e7d32", "Exotic" = "#c62828")) +
  labs(
    title = "Taxonomic Diversity Scaling Profiles",
    x = "Order of Diversity (Hill Numbers)",
    y = "Effective Number of Species (Mean)",
    color = "Origin of Species"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold", size = 13)
  ) -> profile.all

profile.all

Intrpretation: Rare species (q=0) comprise the most of the diversiy of native species but not that exotic repond fo up to 30-40% of all species. However, common and dominant species are exotic, suggesting that ecosystem fucntioning might be heavily driven by exotic species.

5. Regression Modeling & Intercept Strengths

We map native versus exotic variables horizontally using facet_grid() to observe shifting interaction trends directly.

regression_data <- combined_wide %>%
  pivot_longer(
    cols = c(q0, q1, q2),
    names_to = "q_level",
    values_to = "diversity_value"
  ) %>%
  pivot_wider(
    names_from = origin,
    values_from = diversity_value
  ) %>%
  mutate(
    q_level = factor(q_level, levels = c("q0", "q1", "q2"),
                     labels = c("q0 (Rare species)", "q1 (Common species)", "q2 (Dominant species)"))
  )

ggplot(regression_data, aes(x = Exotic, y = Native)) +
  geom_jitter(alpha = 0.6, color = "#2b5c8f", size = 2.5) +
  geom_smooth(method = "lm", formula = y ~ x, color = "#c62828", fill = "gray85", size = 1.2) +
  facet_grid(. ~ q_level, scales = "free") +
  labs(
    x = "Exotic Diversity",
    y = "Native Diversity"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    strip.text = element_text(face = "bold", size = 11),
    strip.background = element_rect(fill = "gray95", color = "gray80"),
    panel.spacing = unit(1.5, "lines"),
    plot.title = element_text(face = "bold", size = 13)
  )

Mixed-Effects Slope Comparison

To determine if regression slopes differ significantly across Hill values while controlling for repeated measures within plots, we implement a linear mixed-effects framework.

regression_data_model <- combined_wide %>%
  pivot_longer(cols = c(q0, q1, q2), names_to = "q_level", values_to = "diversity_value") %>%
  pivot_wider(names_from = origin, values_from = diversity_value) %>%
  mutate(q_level = factor(q_level, levels = c("q0", "q1", "q2")))

mixed_model <- lmer(Native ~ Exotic * q_level + (1 | plot_id), data = regression_data_model)

cat(" Type III ANOVA Results \n")

 Type III ANOVA Results 

Anova(mixed_model, type = "III")

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: Native
                Chisq Df Pr(>Chisq)    
(Intercept)    78.645  1  < 2.2e-16 ***
Exotic         51.263  1  8.077e-13 ***
q_level        41.799  2  8.386e-10 ***
Exotic:q_level 15.688  2   0.000392 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n Estimated Slopes (Beta Coefficients) per q-level \n")

 Estimated Slopes (Beta Coefficients) per q-level 

slopes_output <- emtrends(mixed_model, ~ q_level, var = "Exotic")
print(slopes_output)

 q_level Exotic.trend    SE  df lower.CL upper.CL
 q0            0.5993 0.084 272    0.434    0.765
 q1            0.3846 0.119 289    0.150    0.619
 q2            0.0815 0.126 286   -0.167    0.330

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

cat("\n Pairwise Comparisons of Slopes (Beta Strength) \n")

 Pairwise Comparisons of Slopes (Beta Strength) 

pairs(slopes_output)

 contrast estimate    SE  df t.ratio p.value
 q0 - q1     0.215 0.122 223   1.764  0.1843
 q0 - q2     0.518 0.131 227   3.954  0.0003
 q1 - q2     0.303 0.133 197   2.283  0.0606

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 3 estimates 

Contrast	Estimate (\(\beta\) Difference)	Standard Error (SE)	Degrees of Freedom (df)	t-ratio	p-value	Signif.
q0 - q1	0.0842	0.0899	202	0.936	0.6179	ns
q0 - q2	0.2350	0.0928	201	2.533	0.0322	*
q1 - q2	0.1509	0.0928	196	1.627	0.2368	ns

Conclusion: Across rare and common species metrics (\(q_0\) and \(q_1\)), exotic richness follows native baseline capacities without clear signifiers of exclusion. However, the strength of this relationship changes significantly when comparing species richness (\(q_0\)) against dominant profiles (\(q_2\)), where exotic dominance tracks distinctly higher.

6. Multivariate Cross-Boundary Species Correlations

Community Wide Matrix Preparation

We build a flat wide table tracking species across all unique plots for ecological validation.

matrix_data <- plot_species %>%
  mutate(origin_status = ifelse(type == "Native", "Native", "Exotic")) %>%
  dplyr::select(species_name, origin_status, plot_id, cover_pct)

master_wide_matrix <- matrix_data %>%
  pivot_wider(names_from = plot_id, values_from = cover_pct, values_fill = 0) %>%
  arrange(origin_status, species_name)

print(master_wide_matrix[1:5, 1:5])

# A tibble: 5 × 5
  species_name            origin_status    A1   A10    A2
  <chr>                   <chr>         <dbl> <dbl> <dbl>
1 Acanthospermum hispidum Exotic        0      0    1.04 
2 Achyranthes aspera      Exotic        0      0    3.54 
3 Ageratum conyzoides     Exotic        0.417  0    0.208
4 Bidens pilosa           Exotic        4.17   2.71 0.833
5 Cirsium vulgare         Exotic        0.667  1.67 1.92 

Multiple Cross-Boundary Significance Mapping

We test every individual Native-to-Exotic species coverage intersection using Spearman Rank calculations due to zero-inflated tracking dimensions.

community_wide <- plot_species %>%
  mutate(origin_status = ifelse(type == "Native", "Native", "Exotic")) %>%
  dplyr::select(plot_id, species_name, cover_pct) %>%
  pivot_wider(names_from = species_name, values_from = cover_pct, values_fill = 0) %>%
  as.data.frame()

rownames(community_wide) <- community_wide$plot_id
community_wide$plot_id <- NULL

native_spp <- intersect(unique(plot_species$species_name[plot_species$type == "Native"]), colnames(community_wide))
exotic_spp <- intersect(unique(plot_species$species_name[plot_species$type != "Native"]), colnames(community_wide))

test_res <- cor.mtest(community_wide, method = "spearman", conf.level = 0.95)
full_r <- cor(community_wide, method = "spearman")
full_p <- test_res$p

rownames(full_p) <- colnames(community_wide)
colnames(full_p) <- colnames(community_wide)

sub_r <- full_r[native_spp, exotic_spp, drop = FALSE]
sub_p <- full_p[native_spp, exotic_spp, drop = FALSE]

valid_rows <- apply(sub_r, 1, function(x) !all(is.na(x)))
valid_cols <- apply(sub_r, 2, function(x) !all(is.na(x)))
sub_r <- sub_r[valid_rows, valid_cols, drop = FALSE]
sub_p <- sub_p[valid_rows, valid_cols, drop = FALSE]

corrplot(
  sub_r,
  method = "color",
  p.mat = sub_p,
  sig.level = 0.05,
  insig = "label_sig",
  pch.col = "black",
  pch.cex = 0.8,
  tl.col = "black",
  tl.srt = 45,
  tl.cex = 0.7,
  mar = c(0, 0, 2, 0),
  title = "Significant Interactions (Alpha = 0.05)"
)

Quantifying Significant Relationships

We summarize the total volume of positive, negative, and neutral species-to-species interactions.

summary_table <- data.frame(
  r_value = as.vector(sub_r),
  p_value = as.vector(sub_p)
) %>%
  filter(!is.na(r_value) & !is.na(p_value)) %>%
  mutate(
    interaction_type = case_when(
      p_value <= 0.05 & r_value < 0 ~ "Significant Negative (Competitive)",
      p_value <= 0.05 & r_value > 0 ~ "Significant Positive (Facilitative)",
      TRUE ~ "Non-Significant"
    )
  )

interaction_summary <- summary_table %>%
  group_by(interaction_type) %>%
  summarise(Count = n(), .groups = "drop") %>%
  mutate(Percentage = round((Count / sum(Count)) * 100, 2))

kable(interaction_summary, col.names = c("Interaction Category", "Count", "Percentage (%)"), 
      caption = "Summary Metrics of Evaluated Inter-Species Relationships")

Summary Metrics of Evaluated Inter-Species Relationships

Interaction Category	Count	Percentage (%)
Non-Significant	1069	86.77
Significant Negative (Competitive)	33	2.68
Significant Positive (Facilitative)	130	10.55

7. Individual Exotic Species Interaction Profiles

We evaluate which particular exotic entities hold the widest tracking footprint of significant correlations across native plants.

interaction_df <- as.data.frame(as.table(sub_r)) %>% rename(Native_Species = Var1, Exotic_Species = Var2, r_value = Freq)
p_values_df    <- as.data.frame(as.table(sub_p)) %>% rename(Native_Species = Var1, Exotic_Species = Var2, p_value = Freq)

sig_interactions <- interaction_df %>%
  left_join(p_values_df, by = c("Native_Species", "Exotic_Species")) %>%
  filter(p_value <= 0.05) %>%
  mutate(correlation_direction = ifelse(r_value < 0, "Negative (Competitive)", "Positive (Facilitative)"))

exotic_order <- sig_interactions %>%
  group_by(Exotic_Species) %>%
  summarise(total_sig = n(), .groups = "drop")

sig_interactions <- sig_interactions %>%
  left_join(exotic_order, by = "Exotic_Species") %>%
  mutate(Exotic_Species = reorder(Exotic_Species, total_sig))

ggplot(sig_interactions, aes(x = r_value, y = Exotic_Species, color = correlation_direction)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50", size = 0.8) +
  geom_jitter(width = 0, height = 0.2, size = 3, alpha = 0.75) +
  scale_color_manual(values = c("Negative (Competitive)" = "#d32f2f", "Positive (Facilitative)" = "#1976d2")) +
  xlim(-0.5, 0.5) +
  labs(
    x = "Spearman Correlation Coefficient (r)",
    y = "Exotic Species",
    color = "Correlation Direction"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.minor = element_blank(),
    axis.text.y = element_text(face = "italic"),
    plot.title = element_text(face = "bold", size = 13),
    legend.position = "bottom"
  )

Ecological Insight: The positive interaction frequencies distinctly outweigh competitive exclusions. This indicates that during early natural regeneration, shared environmental niches and facilitation play a larger role in shaping community structure than direct competition.

8. Multiplicative Regional Beta Diversity Partitioning

To determine if exotic species are driving community homogenization (reducing variation among plots) or differentiation (increasing variation among plots), we partition diversity metrics into Alpha, Beta, and Gamma components.

Partitioning Matrix Executions

native_matrix <- plot_species %>%
  filter(type == "Native") %>%
  dplyr::select(plot_id, species_name, cover_pct) %>%
  pivot_wider(names_from = species_name, values_from = cover_pct, values_fill = 0) %>%
  as.data.frame()
rownames(native_matrix) <- native_matrix$plot_id
native_matrix$plot_id <- NULL

combined_matrix <- plot_species %>%
  dplyr::select(plot_id, species_name, cover_pct) %>%
  pivot_wider(names_from = species_name, values_from = cover_pct, values_fill = 0) %>%
  as.data.frame()
rownames(combined_matrix) <- combined_matrix$plot_id
combined_matrix$plot_id <- NULL

get_hill_partition <- function(matrix_input) {
  q0_part <- hill_taxa_parti(matrix_input, q = 0)
  q1_part <- hill_taxa_parti(matrix_input, q = 1)
  q2_part <- hill_taxa_parti(matrix_input, q = 2)
  
  data.frame(
    Order  = c("q = 0", "q = 1", "q = 2"),
    Alpha  = c(q0_part$TD_alpha, q1_part$TD_alpha, q2_part$TD_alpha),
    Gamma  = c(q0_part$TD_gamma, q1_part$TD_gamma, q2_part$TD_gamma),
    Beta   = c(q0_part$TD_beta,  q1_part$TD_beta,  q2_part$TD_beta)
  )
}

native_partition   <- get_hill_partition(native_matrix) %>% mutate(Community = "Native Only")
combined_partition <- get_hill_partition(combined_matrix) %>% mutate(Community = "Combined (Nat + Exo)")

comparison_table <- bind_rows(native_partition, combined_partition) %>%
  dplyr::select(Community, Order, Alpha, Gamma, Beta) %>%
  arrange(Order, Community)

comparison_table %>%
  rename(
    `Community Profile` = Community,
    `Diversity Order` = Order,
    `Mean Alpha (α)` = Alpha,
    `Regional Gamma (γ)` = Gamma,
    `Beta Diversity (β)` = Beta
  ) %>%
  kable(format = "markdown", digits = 3, caption = "Multiplicative Hill Diversity Partitioning Summary")

Multiplicative Hill Diversity Partitioning Summary

Community Profile	Diversity Order	Mean Alpha (α)	Regional Gamma (γ)	Beta Diversity (β)
Combined (Nat + Exo)	q = 0	25.620	78.000	3.044
Native Only	q = 0	14.780	56.000	3.789
Combined (Nat + Exo)	q = 1	8.378	11.016	1.315
Native Only	q = 1	4.367	5.606	1.284
Combined (Nat + Exo)	q = 2	4.069	4.632	1.138
Native Only	q = 2	2.478	2.729	1.101

Here, it is clear that exotic adds to alpha and gamma-diversity at all levels (0,1 and 2) and reduces beta-diversity at all levels, too… more prominently to rare species (q=0)

Beta-diversity

or presence/absence data, betapart uses the Baselga decomposition:

βsor=βsim+βsne

where:

βsim = turnover component βsne = nestedness-resultant component βsor = total Sørensen beta diversity

## Presence-absence matrices
native_pa <- native_matrix
native_pa[native_pa > 0] <- 1

combined_pa <- combined_matrix
combined_pa[combined_pa > 0] <- 1

## Beta partitioning
native_beta <- beta.multi(native_pa, index.family = "sorensen")

combined_beta <- beta.multi(combined_pa, index.family = "sorensen")

## Summary table
beta_partition_summary <- data.frame(
  Community = c("Native only",
                "Native + Exotic"),
  Total_Beta = c(native_beta$beta.SOR,
                 combined_beta$beta.SOR),
  Turnover = c(native_beta$beta.SIM,
               combined_beta$beta.SIM),
  Nestedness = c(native_beta$beta.SNE,
                 combined_beta$beta.SNE)
)

knitr::kable(
  beta_partition_summary,
  digits = 3,
  caption = "Beta diversity partitioning into turnover and nestedness components"
)

Beta diversity partitioning into turnover and nestedness components

Community	Total_Beta	Turnover	Nestedness
Native only	0.947	0.917	0.030
Native + Exotic	0.938	0.906	0.032

Species turnover is high, even accounting for exotic species

NMDS

library(dplyr)
library(tidyr)
library(vegan)
library(ggplot2)
library(viridis)

# ------------------------------------------------------------
# COMMUNITY MATRIX (COVER ABUNDANCE), ALL SPECIES
# combined_matrix was already built in Section 8 (plot x species cover_pct, native + exotic)
# ------------------------------------------------------------

# ------------------------------------------------------------
# NMDS (BRAY-CURTIS, ABUNDANCE-WEIGHTED)
# ------------------------------------------------------------
set.seed(123)

nmds_combined <- metaMDS(
  combined_matrix,
  distance = "bray",
  k = 2,
  trymax = 200
)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.2442817 
Run 1 stress 0.2432797 
... New best solution
... Procrustes: rmse 0.0107429  max resid 0.07206244 
Run 2 stress 0.2724609 
Run 3 stress 0.2459239 
Run 4 stress 0.2536431 
Run 5 stress 0.2482237 
Run 6 stress 0.24897 
Run 7 stress 0.2460263 
Run 8 stress 0.2432805 
... Procrustes: rmse 0.0003377584  max resid 0.002183741 
... Similar to previous best
Run 9 stress 0.2482441 
Run 10 stress 0.2435766 
... Procrustes: rmse 0.02797078  max resid 0.257386 
Run 11 stress 0.2668087 
Run 12 stress 0.2468991 
Run 13 stress 0.2435099 
... Procrustes: rmse 0.006712688  max resid 0.0570437 
Run 14 stress 0.2432679 
... New best solution
... Procrustes: rmse 0.003132957  max resid 0.01892986 
Run 15 stress 0.2515268 
Run 16 stress 0.2509864 
Run 17 stress 0.2548257 
Run 18 stress 0.2546691 
Run 19 stress 0.2549783 
Run 20 stress 0.2475056 
Run 21 stress 0.2520907 
Run 22 stress 0.2524709 
Run 23 stress 0.249381 
Run 24 stress 0.2664936 
Run 25 stress 0.2433672 
... Procrustes: rmse 0.02775623  max resid 0.2585561 
Run 26 stress 0.2454455 
Run 27 stress 0.2505969 
Run 28 stress 0.2482504 
Run 29 stress 0.2466012 
Run 30 stress 0.2482232 
Run 31 stress 0.2486501 
Run 32 stress 0.2439886 
Run 33 stress 0.243286 
... Procrustes: rmse 0.02744327  max resid 0.2581131 
Run 34 stress 0.253868 
Run 35 stress 0.2440312 
Run 36 stress 0.2582477 
Run 37 stress 0.2482449 
Run 38 stress 0.2452475 
Run 39 stress 0.2470245 
Run 40 stress 0.246785 
Run 41 stress 0.2456615 
Run 42 stress 0.2475998 
Run 43 stress 0.2452673 
Run 44 stress 0.2482503 
Run 45 stress 0.2457945 
Run 46 stress 0.2743839 
Run 47 stress 0.2756156 
Run 48 stress 0.2488786 
Run 49 stress 0.2435061 
... Procrustes: rmse 0.006913872  max resid 0.05740086 
Run 50 stress 0.2485016 
Run 51 stress 0.2461108 
Run 52 stress 0.2563147 
Run 53 stress 0.2456616 
Run 54 stress 0.2527189 
Run 55 stress 0.2449955 
Run 56 stress 0.2485548 
Run 57 stress 0.2442993 
Run 58 stress 0.2624948 
Run 59 stress 0.2679726 
Run 60 stress 0.2467027 
Run 61 stress 0.2463945 
Run 62 stress 0.2459409 
Run 63 stress 0.2482232 
Run 64 stress 0.2445321 
Run 65 stress 0.2481346 
Run 66 stress 0.2471875 
Run 67 stress 0.2445321 
Run 68 stress 0.2492658 
Run 69 stress 0.2435007 
... Procrustes: rmse 0.007408962  max resid 0.05773538 
Run 70 stress 0.2443088 
Run 71 stress 0.2481395 
Run 72 stress 0.2434965 
... Procrustes: rmse 0.006472729  max resid 0.05685346 
Run 73 stress 0.2494021 
Run 74 stress 0.243502 
... Procrustes: rmse 0.007358048  max resid 0.05773801 
Run 75 stress 0.2739188 
Run 76 stress 0.2442712 
Run 77 stress 0.2460175 
Run 78 stress 0.2468608 
Run 79 stress 0.2476013 
Run 80 stress 0.2824677 
Run 81 stress 0.2435101 
... Procrustes: rmse 0.00736949  max resid 0.05790038 
Run 82 stress 0.2457795 
Run 83 stress 0.2435767 
... Procrustes: rmse 0.02808765  max resid 0.2572498 
Run 84 stress 0.2454827 
Run 85 stress 0.2446497 
Run 86 stress 0.2479641 
Run 87 stress 0.2714292 
Run 88 stress 0.2434965 
... Procrustes: rmse 0.006615339  max resid 0.05736531 
Run 89 stress 0.2472816 
Run 90 stress 0.2438004 
Run 91 stress 0.2485542 
Run 92 stress 0.2476216 
Run 93 stress 0.2505975 
Run 94 stress 0.2482681 
Run 95 stress 0.2437848 
Run 96 stress 0.246782 
Run 97 stress 0.2457265 
Run 98 stress 0.2446688 
Run 99 stress 0.2605548 
Run 100 stress 0.2475862 
Run 101 stress 0.2489696 
Run 102 stress 0.2441657 
Run 103 stress 0.2525667 
Run 104 stress 0.2489767 
Run 105 stress 0.2435126 
... Procrustes: rmse 0.006976336  max resid 0.05833751 
Run 106 stress 0.2482466 
Run 107 stress 0.2482231 
Run 108 stress 0.2435009 
... Procrustes: rmse 0.007443223  max resid 0.05805689 
Run 109 stress 0.4114804 
Run 110 stress 0.2442649 
Run 111 stress 0.2481378 
Run 112 stress 0.2435018 
... Procrustes: rmse 0.007299228  max resid 0.05733767 
Run 113 stress 0.2448624 
Run 114 stress 0.2468987 
Run 115 stress 0.267445 
Run 116 stress 0.2434961 
... Procrustes: rmse 0.006492487  max resid 0.05677175 
Run 117 stress 0.2596071 
Run 118 stress 0.2550836 
Run 119 stress 0.2454963 
Run 120 stress 0.247598 
Run 121 stress 0.2459239 
Run 122 stress 0.2488829 
Run 123 stress 0.2465219 
Run 124 stress 0.2583626 
Run 125 stress 0.2466422 
Run 126 stress 0.2457182 
Run 127 stress 0.262545 
Run 128 stress 0.245587 
Run 129 stress 0.2442653 
Run 130 stress 0.2554766 
Run 131 stress 0.251456 
Run 132 stress 0.2432673 
... New best solution
... Procrustes: rmse 0.0002016811  max resid 0.00168517 
... Similar to previous best
*** Best solution repeated 1 times

cat("\nNMDS Stress =", round(nmds_combined$stress, 3), "\n")

NMDS Stress = 0.243 

# ------------------------------------------------------------
# EXOTIC CONTRIBUTION TO TOTAL COVER, PER PLOT
# ------------------------------------------------------------
plot_contribution <- plot_species %>%
  mutate(origin_status = ifelse(type == "Native", "Native", "Exotic")) %>%
  group_by(plot_id) %>%
  summarise(
    # NOTE: sum(), not mean() — this is RELATIVE COVER (share of cumulative
    # plot cover attributable to exotics), the standard approach for cover-class
    # vegetation data where summed cover across species can exceed 100% due to
    # overlapping layers/patchy co-occurring species. mean() would not be
    # bounded 0-100% and would shift with species count, not actual dominance.
    total_cover  = sum(cover_pct, na.rm = TRUE),
    exotic_cover = sum(cover_pct[origin_status == "Exotic"], na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(exotic_contrib_pct = 100 * exotic_cover / total_cover)

# ------------------------------------------------------------
# MERGE WITH NMDS SITE SCORES
# ------------------------------------------------------------
nmds_scores_combined <- as.data.frame(scores(nmds_combined, display = "sites"))
nmds_scores_combined$plot_id <- rownames(nmds_scores_combined)

nmds_scores_combined <- nmds_scores_combined %>%
  left_join(plot_contribution, by = "plot_id")

# ------------------------------------------------------------
# ENVFIT VECTOR OVERLAY (OPTIONAL BUT KEEPS CONSISTENCY WITH SECTION 8)
# ------------------------------------------------------------
fit_contrib <- envfit(
  nmds_combined,
  plot_contribution["exotic_contrib_pct"],
  permutations = 999
)
print(fit_contrib)

***VECTORS

                      NMDS1    NMDS2     r2 Pr(>r)   
exotic_contrib_pct -0.77295  0.63446 0.1452  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Permutation: free
Number of permutations: 999

vec_contrib <- as.data.frame(scores(fit_contrib, display = "vectors"))
vec_contrib$label <- paste0(
  "Exotic contribution\nR\u00b2=", round(fit_contrib$vectors$r, 2),
  "\np=", signif(fit_contrib$vectors$pvals, 2)
)

# ------------------------------------------------------------
# PLOT: SIZE + COLOR BOTH ENCODE EXOTIC COVER CONTRIBUTION
# ------------------------------------------------------------
ggplot(nmds_scores_combined,
       aes(NMDS1, NMDS2,
           size = exotic_contrib_pct,
           color = exotic_contrib_pct)) +

  geom_point(alpha = 0.85) +

  geom_segment(
    data = vec_contrib,
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
    inherit.aes = FALSE,
    arrow = arrow(length = unit(0.25, "cm")),
    linewidth = 1.1,
    colour = "black"
  ) +

  geom_text(
    data = vec_contrib,
    aes(x = NMDS1, y = NMDS2, label = label),
    inherit.aes = FALSE,
    hjust = -0.1, vjust = -0.2, size = 4
  ) +

  scale_color_viridis_c(option = "plasma", name = "Exotic cover\ncontribution (%)") +
  scale_size_continuous(range = c(2, 9), name = "Exotic cover\ncontribution (%)") +

  labs(
    title = "NMDS of Plant Community Composition (All Species)",
    subtitle = paste0(
      "Bray-Curtis on cover abundance | Stress = ",
      round(nmds_combined$stress, 3)
    ),
    x = "NMDS1",
    y = "NMDS2"
  ) +

  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold"),
    legend.position = "right"
  )

Interpretation: NMDS stress is too high, so it should be interpreted cautiously but it suggests that exotic species drives plant community to somewhere the upper-left corner of the graph.