Introduction
The aquats package provides a robust, streamlined
statistical and visualization toolkit tailored specifically for
fisheries science, stock assessment workflows, and aquatic ecology. It
bridges standard parametric models (ANOVA, ANCOVA, MANOVA) with advanced
ecological tools, including non-linear morphometric regressions,
generalized linear models (GLMs) for count data, customizable
correlation heatmaps, and integrated Mantel test network linkages.
Load the package and required dependencies to begin:
0.1 Two-Sample Comparison (t-Test and F-Test)
0.1.1 Conceptual Framework
Comparing continuous ecological, biometric, or water-quality
parameters between two distinct sample groups (e.g., Reference
vs. Impacted sites, Upstream vs. Downstream) is a foundational
statistical task. While Student’s classical \(t\)-test assumes equal population
variances, ecological field data frequently violate this condition. The
compare_two_groups() function automates variance screening
via Fisher’s \(F\)-test, executes
either Student’s or Welch’s \(t\)-test
with Satterthwaite degrees of freedom, and packages diagnostics, effect
sizes (Cohen’s \(d\)), and graphics
into a single call.
0.1.2 Hypotheses and Variance Assumptions
The procedure simultaneously tests: 1. Variance Homogeneity
(Fisher’s \(F\)-test): \[H_0: \sigma_1^2 = \sigma_2^2 \quad \text{vs.}
\quad H_1: \sigma_1^2 \neq \sigma_2^2\] 2. Mean Equality
(\(t\)-test): \[H_0: \mu_1 = \mu_2 \quad \text{vs.} \quad H_1:
\mu_1 \neq \mu_2\] When var_equal = "auto", Welch’s
correction is adopted automatically if the \(F\)-test rejects homoscedasticity (\(p < 0.05\)).
0.1.3 In-Vignette Data Generation
We simulate dissolved oxygen (DO, \(\text{mg}\cdot\text{L}^{-1}\)) monitoring records across two river stations exhibiting unequal variances:
0.1.4 Visualization Styles
0.1.4.1 Style A: Publication Boxplot with Jitter
Default publication aesthetic displaying distribution medians, interquartile ranges, raw observations as jittered points, and an annotated significance bracket:
res_box <- compare_two_groups(
data = river_survey,
response_var = "DO_mgL",
group_var = "Station",
var_equal = "auto",
plot_type = "box",
palette = "Dark2",
title = "Style A: Two-Sample Boxplot with Significance Bracket"
)
print(res_box$Plot)
0.1.4.2 Style B: Faceted / Density-Informed Violin Plot
Displays underlying data density distributions combined with an internal box-whisker geometry:
res_violin <- compare_two_groups(
data = river_survey,
response_var = "DO_mgL",
group_var = "Station",
var_equal = "auto",
plot_type = "violin",
palette = "Set2",
title = "Style B: Density Violin with Internal Quartiles"
)
print(res_violin$Plot)
0.1.4.3 Style C: Minimalist Point and Mean-CI Error Bars
A clean point-and-interval layout emphasizing group means and 95% confidence intervals against jittered points:
res_jitter <- compare_two_groups(
data = river_survey,
response_var = "DO_mgL",
group_var = "Station",
var_equal = "auto",
plot_type = "jitter",
palette = "Set1",
title = "Style C: Group Means and 95% CI Error Bars"
)
print(res_jitter$Plot)
0.1.5 Model Inspection and Normality Diagnostics
Review the automated assumption checks, formal test statistics, and effect sizes:
# 1. Distribution Normality per group (Shapiro-Wilk)
knitr::kable(res_box$Normality_Test, caption = "Shapiro-Wilk Normality Evaluation")| Group | Statistic_W | p_value | Normal_05 |
|---|---|---|---|
| Downstream | 0.9716 | 0.6866 | Yes |
| Upstream | 0.9499 | 0.2495 | Yes |
# 2. Fisher's F-test for Variance Equality
knitr::kable(res_box$F_Test, caption = "Fisher's F-Test for Homogeneity of Variances")| Ratio_Var | F_Statistic | df_num | df_denom | p_value | Equal_Var_05 |
|---|---|---|---|---|---|
| 2.297 | 2.297 | 24 | 24 | 0.0468 | No |
| Test_Type | Alternative | t_Statistic | df | p_value | CI_Lower | CI_Upper |
|---|---|---|---|---|---|---|
| Welch Two Sample t-test | two.sided | -7.6683 | 41.57 | 0 | -2.5276 | -1.4741 |
# 4. Group Descriptives and Cohen's d
knitr::kable(res_box$Summary_Table, caption = "Group Statistics and Standardized Effect Sizes")| Group | N | Mean | SD | SE | CI_Lower | CI_Upper | Mean_Diff | Cohen_d |
|---|---|---|---|---|---|---|---|---|
| Downstream | 25 | 5.902 | 1.089 | 0.218 | 5.453 | 6.352 | -2.001 | 2.169 |
| Upstream | 25 | 7.903 | 0.719 | 0.144 | 7.607 | 8.200 | 2.001 | 2.169 |
0.2 Contingency Analysis and Goodness-of-Fit (Chi-Square Test)
0.2.1 Conceptual Framework
While analysis of variance and \(t\)-tests evaluate continuous measurements
across groups, ecological and aquatic assessments frequently involve
categorical observations (e.g., habitat substrate types, feeding guilds,
species maturity status, or pollution compliance tiers). The
test_chisq() function tests for significant associations
between categorical factors (Test of Independence) or evaluates whether
an observed frequency distribution conforms to theoretical expectations
(Goodness-of-Fit).
0.2.2 Hypotheses and Effect Sizes
For a two-way contingency table (\(r \times c\)): 1. Hypothesis of Independence: \[H_0: \text{Variable } X \text{ and Variable } Y \text{ are independent}\] \[H_1: \text{Variable } X \text{ and Variable } Y \text{ are associated}\] 2. Standardized Residuals: To determine which specific combinations drive a significant test result, standardized Pearson residuals (\(r_{ij}\)) are computed: \[r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij} (1 - p_{i\cdot}) (1 - p_{\cdot j})}}\] Cells with \(\vert{}r_{ij}\vert{} > 2.0\) indicate deviations significant at \(\alpha = 0.05\). 3. Strength of Association: Quantified using Cramér’s \(V\): \[V = \sqrt{\frac{\chi^2}{N \cdot \min(r - 1, c - 1)}}\]
0.2.3 In-Vignette Data Generation
We simulate fish maturity status across three sampling river reaches:
set.seed(42)
fish_maturity <- data.frame(
Reach = factor(rep(c("Upstream", "Midstream", "Estuary"), times = c(60, 70, 70))),
Maturity = factor(c(
sample(c("Immature", "Maturing", "Spent"), 60, replace = TRUE, prob = c(0.60, 0.30, 0.10)),
sample(c("Immature", "Maturing", "Spent"), 70, replace = TRUE, prob = c(0.35, 0.45, 0.20)),
sample(c("Immature", "Maturing", "Spent"), 70, replace = TRUE, prob = c(0.15, 0.35, 0.50))
))
)0.2.4 Visualization Styles
0.2.4.1 Style A: 100% Stacked Proportion Bar Chart
Displays the relative conditional percentage breakdown of categories across groups:
res_stacked <- test_chisq(
data = fish_maturity,
x = "Reach",
y = "Maturity",
plot_type = "stacked",
palette = "Set2",
title = "Style A: Maturity Proportions Across River Reaches"
)
print(res_stacked$Plot)
0.2.5 Model Inspection and Residual Diagnostics
Review test statistics, expected cell frequencies, and standardized residuals:
# 1. Statistical Test Summary and Effect Size (Cramér's V)
knitr::kable(res_stacked$Chi_Square_Test, caption = "Chi-Square Test Summary and Association Strength")| Test_Type | Chi_Square | df | p_value | Cramers_V | Phi |
|---|---|---|---|---|---|
| Pearson’s Chi-squared test | 30.9278 | 4 | 0 | 0.2781 | NA |
# 2. Observed Contingency Table
knitr::kable(res_stacked$Observed_Counts, caption = "Observed Cell Frequencies")| Reach | Immature | Maturing | Spent | Total |
|---|---|---|---|---|
| Estuary | 8 | 29 | 33 | 70 |
| Midstream | 30 | 29 | 11 | 70 |
| Upstream | 29 | 18 | 13 | 60 |
| Total | 67 | 76 | 57 | 200 |
# 3. Expected Cell Frequencies
knitr::kable(res_stacked$Expected_Counts, caption = "Expected Frequencies Under Independence")| Reach | Immature | Maturing | Spent |
|---|---|---|---|
| Estuary | 23.45 | 26.6 | 19.95 |
| Midstream | 23.45 | 26.6 | 19.95 |
| Upstream | 20.10 | 22.8 | 17.10 |
# 4. Standardized Pearson Residuals (|r| > 2 denotes significant deviation)
knitr::kable(res_stacked$Standardized_Residuals, caption = "Standardized Pearson Residuals")| Reach | Var1 | Var2 | Freq |
|---|---|---|---|
| 1 | Estuary | Immature | -4.853 |
| 2 | Midstream | Immature | 2.057 |
| 3 | Upstream | Immature | 2.910 |
| 4 | Estuary | Maturing | 0.733 |
| 5 | Midstream | Maturing | 0.733 |
| 6 | Upstream | Maturing | -1.526 |
| 7 | Estuary | Spent | 4.286 |
| 8 | Midstream | Spent | -2.939 |
| 9 | Upstream | Spent | -1.401 |
While the Chi-Square test evaluates associations among discrete frequency classes, comparing continuous metric responses across three or more treatment groups requires partitioning continuous variance, addressed in the next section.
0.3 One-Way Analysis of Variance (ANOVA)
The one_way_anova() function provides an end-to-end
analytical and graphical pipeline for evaluating differences in a
continuous response variable across levels of a single categorical
factor. In aquatic and fisheries research, this is commonly employed to
compare water quality parameters across sampling stations, fish
morphometric indices across river reaches, or nutrient concentrations
across seasons.
The function handles defensive data cleaning, model fitting, computation of global F-tests, extraction of Estimated Marginal Means (EMMs), pairwise post-hoc contrasts with compact letter displays (CLD), and automated generation of publication-ready ggplot2 figures.
0.3.1 Conceptual Framework and Statistical Formulation
A classical fixed-effects One-Way ANOVA partitions the total variance in the response variable into between-group (explained) variance and within-group (unexplained error) variance:
\[y_{ij} = \mu + \alpha_i + \varepsilon_{ij}\]
where: - \(y_{ij}\) represents the \(j\)-th observation from the \(i\)-th treatment or habitat group (\(i = 1, 2, \dots, k\); \(j = 1, 2, \dots, n_i\)), - \(\mu\) represents the grand population mean, - \(\alpha_i = \mu_i - \mu\) is the differential effect associated with group \(i\), subject to the constraint \(\sum n_i \alpha_i = 0\), - \(\varepsilon_{ij} \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)\) represents random residual error.
0.3.1.1 Hypotheses and Test Statistics
The omnibus null hypothesis states that all group means are identical:
\[H_0: \mu_1 = \mu_2 = \dots = \mu_k\]
\[H_1: \text{At least one group mean differs from the others.}\]
The hypothesis is evaluated using the F-statistic:
\[F = \frac{\text{MS}_{\text{between}}}{\text{MS}_{\text{within}}} = \frac{\text{SS}_{\text{between}} / (k - 1)}{\text{SS}_{\text{within}} / (N - k)}\]
When \(H_0\) is rejected (\(p < 0.05\)), pairwise differences are resolved using Tukey’s Honest Significant Difference (HSD) test:
\[q = \frac{\bar{y}_i - \bar{y}_j}{\sqrt{\frac{\text{MS}_{\text{within}}}{2} \left(\frac{1}{n_i} + \frac{1}{n_j}\right)}}\]
Groups sharing the same letter in the Compact Letter Display (CLD) are not statistically distinguishable at the specified significance threshold (\(\alpha = 0.05\)).
0.3.2 In-Vignette Data Generation
To illustrate the functionality, we simulate a standard riverine
biomonitoring dataset evaluating Dissolved Oxygen (DO, mg/L) across
three river zones: Upper Basin (forested reference),
Mid Reach (agricultural run-off), and
Lower Estuary (downstream industrial/tidal influence).
set.seed(101)
n_obs <- 30
do_survey <- data.frame(
Zone = factor(
rep(c("Upper Basin", "Mid Reach", "Lower Estuary"), each = n_obs),
levels = c("Upper Basin", "Mid Reach", "Lower Estuary")
),
DO_mgL = c(
rnorm(n_obs, mean = 8.50, sd = 0.45),
rnorm(n_obs, mean = 6.80, sd = 0.55),
rnorm(n_obs, mean = 4.90, sd = 0.60)
)
)
head(do_survey, 6)
#> Zone DO_mgL
#> 1 Upper Basin 8.353284
#> 2 Upper Basin 8.748608
#> 3 Upper Basin 8.196275
#> 4 Upper Basin 8.596462
#> 5 Upper Basin 8.639846
#> 6 Upper Basin 9.0282850.3.3 Visualization Styles
The one_way_anova() function supports multiple layout
types depending on reporting requirements.
0.3.3.1 Style A: Publication Boxplot with Jitter and Mean Marker (Default)
This style shows median lines, interquartile ranges, raw sample
distributions via jittered points, and group means displayed as diamonds
(show_mean = TRUE). Tukey post-hoc letters are
automatically placed above each group distribution:
res_box <- one_way_anova(
data = do_survey,
factor_var = "Zone",
numeric_var = "DO_mgL",
plot_type = "boxplot",
add_jitter = TRUE,
show_mean = TRUE,
mean_color = "darkred",
sig_display = "letters",
color_palette = "Blues"
)
res_box$Plot
0.3.3.2 Style B: Grouped Barplot with Standard Error Bars (Mean ± SE)
For aquaculture trials or fisheries management reports,
plot_type = "barplot" aggregates observations to their mean
and attaches directional error bars (error_type = "se" or
error_type = "sd"):
res_bar <- one_way_anova(
data = do_survey,
factor_var = "Zone",
numeric_var = "DO_mgL",
plot_type = "barplot",
error_type = "se",
sig_display = "letters",
color_palette = "Blues"
)
res_bar$Plot
0.3.3.3 Style C: Minimalist Point-Range Geometry
plot_type = "pointrange" provides a clean, modern
alternative for comparing group means without column bars:
res_point <- one_way_anova(
data = do_survey,
factor_var = "Zone",
numeric_var = "DO_mgL",
plot_type = "pointrange",
error_type = "se",
sig_display = "letters",
color_palette = "Dark2"
)
res_point$Plot
0.3.4 Model Inspection and Numerical Summaries
The function returns a comprehensive list containing the underlying ANOVA model, post-hoc pairwise comparisons, Estimated Marginal Means, and a summary statistics table:
# Global ANOVA Table (F-statistic and p-value)
res_box$ANOVA_Summary
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Zone 2 184.71 92.36 339.9 <2e-16 ***
#> Residuals 87 23.64 0.27
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Pairwise Tukey HSD Differences
res_box$TukeyHSD
#> Tukey multiple comparisons of means
#> 95% family-wise confidence level
#>
#> Fit: stats::aov(formula = formula_str, data = clean_data)
#>
#> $Zone
#> diff lwr upr p adj
#> Mid Reach-Upper Basin -1.740017 -2.060938 -1.419097 0
#> Lower Estuary-Upper Basin -3.509114 -3.830034 -3.188193 0
#> Lower Estuary-Mid Reach -1.769096 -2.090017 -1.448176 0
# Estimated Marginal Means (EMMs)
res_box$EMMs
#> Zone emmean SE df lower.CL upper.CL
#> Upper Basin 8.46 0.0952 87 8.27 8.65
#> Mid Reach 6.72 0.0952 87 6.53 6.91
#> Lower Estuary 4.95 0.0952 87 4.76 5.14
#>
#> Confidence level used: 0.95
# Computed Summary Statistics Table
res_box$Summary_Table
#> # A tibble: 3 × 6
#> Zone mean_val sd_val se_val error_margin upper_limit
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Upper Basin 8.46 0.376 0.0686 0.0686 8.53
#> 2 Mid Reach 6.72 0.564 0.103 0.103 6.83
#> 3 Lower Estuary 4.95 0.596 0.109 0.109 5.060.4 Two-Way Analysis of Variance (ANOVA)
The two_way_anova() function performs factorial analysis
of variance to simultaneously test the main effects of two categorical
factors as well as their multiplicative interaction on a continuous
aquatic response variable. In fisheries biology and aquatic ecology,
factorial designs are standard practice for assessing seasonal variation
across river basins, depth profiles across salinity gradients, or growth
performance across diet and temperature treatments.
The function handles automated missing value filtering, model parameter estimation, computation of interaction sum-of-squares, post-hoc contrast grouping (via Estimated Marginal Means and Tukey’s HSD), and renders both grouped interaction plots and faceted sub-panel displays.
0.4.1 Conceptual Framework and Statistical Formulation
A classical fixed-effects balanced or unbalanced Two-Way ANOVA models the continuous response \(y_{ijk}\) as:
\[y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk}\]
where: - \(y_{ijk}\) is the \(k\)-th replicate observation within level \(i\) of Factor A and level \(j\) of Factor B (\(i = 1, \dots, a\); \(j = 1, \dots, b\); \(k = 1, \dots, n_{ij}\)), - \(\mu\) is the overall population grand mean, - \(\alpha_i\) is the main effect of Factor A, subject to the identifiability constraint \(\sum n_i \alpha_i = 0\), - \(\beta_j\) is the main effect of Factor B, subject to the identifiability constraint \(\sum n_j \beta_j = 0\), - \((\alpha\beta)_{ij}\) represents the interaction term capturing non-additive synergies or antagonisms between factors, - \(\varepsilon_{ijk} \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)\) is the random error term.
0.4.1.1 Hypotheses and Variance Partitioning
The test partitions the total sum of squares (\(\text{SS}_{\text{Total}}\)) into orthogonal components:
\[\text{SS}_{\text{Total}} = \text{SS}_{A} + \text{SS}_{B} + \text{SS}_{AB} + \text{SS}_{\text{Residual}}\]
The three core hypotheses evaluated by the model are:
Main Effect of Factor A: \[H_0^{(A)}: \alpha_1 = \dots = \alpha_a = 0 \quad \text{vs.} \quad F_A = \frac{\text{MS}_A}{\text{MS}_{\text{Residual}}}\]
Main Effect of Factor B: \[H_0^{(B)}: \beta_1 = \dots = \beta_b = 0 \quad \text{vs.} \quad F_B = \frac{\text{MS}_B}{\text{MS}_{\text{Residual}}}\]
Interaction Effect (\(A \times B\)): \[H_0^{(AB)}: (\alpha\beta)_{11} = \dots = (\alpha\beta)_{ab} = 0 \quad \text{vs.} \quad F_{AB} = \frac{\text{MS}_{AB}}{\text{MS}_{\text{Residual}}}\]
If the interaction effect is statistically significant (\(p < 0.05\)), main effects cannot be interpreted in isolation, and post-hoc pairwise contrasts are computed across the joint interaction levels using Estimated Marginal Means (EMMs).
0.4.2 In-Vignette Data Generation
We simulate a biomonitoring dataset evaluating fish condition factor
(\(K\)) across three river basins
(Upper, Central, Lower) sampled
over two distinct hydrological periods (Wet Season vs
Dry Season):
set.seed(42)
n_reps <- 15
basins <- c("Upper Reach", "Central Basin", "Lower Estuary")
seasons <- c("Wet Season", "Dry Season")
fish_condition <- data.frame(
Basin = factor(rep(rep(basins, each = n_reps), times = 2), levels = basins),
Season = factor(rep(seasons, each = n_reps * length(basins)), levels = seasons),
Condition_Factor = c(
# Wet Season: Upper, Central, Lower
rnorm(n_reps, mean = 1.35, sd = 0.08),
rnorm(n_reps, mean = 1.18, sd = 0.09),
rnorm(n_reps, mean = 1.05, sd = 0.10),
# Dry Season: Upper, Central, Lower
rnorm(n_reps, mean = 1.20, sd = 0.07),
rnorm(n_reps, mean = 0.98, sd = 0.08),
rnorm(n_reps, mean = 0.82, sd = 0.09)
)
)
head(fish_condition, 6)
#> Basin Season Condition_Factor
#> 1 Upper Reach Wet Season 1.459677
#> 2 Upper Reach Wet Season 1.304824
#> 3 Upper Reach Wet Season 1.379050
#> 4 Upper Reach Wet Season 1.400629
#> 5 Upper Reach Wet Season 1.382341
#> 6 Upper Reach Wet Season 1.3415100.4.3 Visualization Styles
The two_way_anova() function renders both grouped
interaction plots (dodged) and faceted sub-panel displays across all
standard geometries.
0.4.3.1 Style A: Dodged Boxplot with Interaction Letters (Default)
This style groups the primary factor along the x-axis, using color fill to contrast the secondary factor side-by-side with raw sample jitters and mean indicators:
res_two_box <- two_way_anova(
data = fish_condition,
factor1_var = "Basin",
factor2_var = "Season",
numeric_var = "Condition_Factor",
plot_type = "boxplot",
add_jitter = TRUE,
show_mean = TRUE,
mean_color = "darkred",
color_palette = "Set2"
)
res_two_box$Plot
0.4.3.2 Style B: Faceted Barplot with Standard Error Bars (Mean ± SE)
By specifying facet_var = "Season", the visualization
cleanly segments the interaction into independent panels, preventing
crowded error bars and overlapping categories:
res_two_bar <- two_way_anova(
data = fish_condition,
factor1_var = "Basin",
factor2_var = "Season",
numeric_var = "Condition_Factor",
facet_var = "Season",
plot_type = "barplot",
error_type = "se",
color_palette = "Set2"
)
res_two_bar$Plot
0.4.3.3 Style C: Minimalist Point-Range Interaction
The point-range display plots mean coordinates and directional standard error bars, making it straightforward to inspect interaction divergence between levels:
res_two_point <- two_way_anova(
data = fish_condition,
factor1_var = "Basin",
factor2_var = "Season",
numeric_var = "Condition_Factor",
plot_type = "pointrange",
error_type = "se",
color_palette = "Dark2"
)
res_two_point$Plot
0.4.4 Model Inspection and Numerical Summaries
The function returns the comprehensive ANOVA table showing degrees of freedom, sum-of-squares, mean-squares, \(F\)-values, and \(p\)-values for both main factors and their interaction:
# Factorial ANOVA Table (Main effects and Basin:Season interaction)
res_two_box$ANOVA_Summary
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Basin 2 2.1146 1.0573 136.277 < 2e-16 ***
#> Season 1 0.6614 0.6614 85.244 2.02e-14 ***
#> Basin:Season 2 0.0048 0.0024 0.312 0.733
#> Residuals 84 0.6517 0.0078
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Post-hoc pairwise differences (Tukey HSD)
res_two_box$TukeyHSD
#> Tukey multiple comparisons of means
#> 95% family-wise confidence level
#>
#> Fit: stats::aov(formula = formula_str, data = clean_data)
#>
#> $Basin
#> diff lwr upr p adj
#> Central Basin-Upper Reach -0.2244000 -0.2786628 -0.1701371 0
#> Lower Estuary-Upper Reach -0.3728956 -0.4271585 -0.3186328 0
#> Lower Estuary-Central Basin -0.1484957 -0.2027586 -0.0942328 0
#>
#> $Season
#> diff lwr upr p adj
#> Dry Season-Wet Season -0.171446 -0.208373 -0.134519 0
#>
#> $`Basin:Season`
#> diff lwr
#> Central Basin:Wet Season-Upper Reach:Wet Season -0.23997417 -0.33377867
#> Lower Estuary:Wet Season-Upper Reach:Wet Season -0.37294488 -0.46674938
#> Upper Reach:Dry Season-Upper Reach:Wet Season -0.18186162 -0.27566612
#> Central Basin:Dry Season-Upper Reach:Wet Season -0.39068738 -0.48449188
#> Lower Estuary:Dry Season-Upper Reach:Wet Season -0.55470804 -0.64851254
#> Lower Estuary:Wet Season-Central Basin:Wet Season -0.13297070 -0.22677521
#> Upper Reach:Dry Season-Central Basin:Wet Season 0.05811255 -0.03569195
#> Central Basin:Dry Season-Central Basin:Wet Season -0.15071321 -0.24451771
#> Lower Estuary:Dry Season-Central Basin:Wet Season -0.31473387 -0.40853837
#> Upper Reach:Dry Season-Lower Estuary:Wet Season 0.19108326 0.09727876
#> Central Basin:Dry Season-Lower Estuary:Wet Season -0.01774250 -0.11154700
#> Lower Estuary:Dry Season-Lower Estuary:Wet Season -0.18176316 -0.27556766
#> Central Basin:Dry Season-Upper Reach:Dry Season -0.20882576 -0.30263026
#> Lower Estuary:Dry Season-Upper Reach:Dry Season -0.37284642 -0.46665092
#> Lower Estuary:Dry Season-Central Basin:Dry Season -0.16402066 -0.25782516
#> upr p adj
#> Central Basin:Wet Season-Upper Reach:Wet Season -0.14616967 0.0000000
#> Lower Estuary:Wet Season-Upper Reach:Wet Season -0.27914038 0.0000000
#> Upper Reach:Dry Season-Upper Reach:Wet Season -0.08805712 0.0000031
#> Central Basin:Dry Season-Upper Reach:Wet Season -0.29688288 0.0000000
#> Lower Estuary:Dry Season-Upper Reach:Wet Season -0.46090354 0.0000000
#> Lower Estuary:Wet Season-Central Basin:Wet Season -0.03916620 0.0011495
#> Upper Reach:Dry Season-Central Basin:Wet Season 0.15191705 0.4670331
#> Central Basin:Dry Season-Central Basin:Wet Season -0.05690871 0.0001515
#> Lower Estuary:Dry Season-Central Basin:Wet Season -0.22092937 0.0000000
#> Upper Reach:Dry Season-Lower Estuary:Wet Season 0.28488776 0.0000009
#> Central Basin:Dry Season-Lower Estuary:Wet Season 0.07606200 0.9937363
#> Lower Estuary:Dry Season-Lower Estuary:Wet Season -0.08795866 0.0000031
#> Central Basin:Dry Season-Upper Reach:Dry Season -0.11502126 0.0000001
#> Lower Estuary:Dry Season-Upper Reach:Dry Season -0.27904192 0.0000000
#> Lower Estuary:Dry Season-Central Basin:Dry Season -0.07021616 0.0000301
# Estimated Marginal Means with compact letter displays
res_two_box$EMMs
#> Basin Season emmean SE df lower.CL upper.CL
#> Upper Reach Wet Season 1.389 0.0227 84 1.344 1.434
#> Central Basin Wet Season 1.149 0.0227 84 1.104 1.194
#> Lower Estuary Wet Season 1.016 0.0227 84 0.971 1.061
#> Upper Reach Dry Season 1.207 0.0227 84 1.162 1.252
#> Central Basin Dry Season 0.998 0.0227 84 0.953 1.043
#> Lower Estuary Dry Season 0.834 0.0227 84 0.789 0.879
#>
#> Confidence level used: 0.95
# Summary statistics table
res_two_box$Summary_Table
#> # A tibble: 6 × 7
#> Basin Season mean_val sd_val se_val error_margin upper_limit
#> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Upper Reach Wet Season 1.39 0.0820 0.0212 0.0212 1.41
#> 2 Upper Reach Dry Season 1.21 0.0762 0.0197 0.0197 1.23
#> 3 Central Basin Wet Season 1.15 0.122 0.0315 0.0315 1.18
#> 4 Central Basin Dry Season 0.998 0.0648 0.0167 0.0167 1.01
#> 5 Lower Estuary Wet Season 1.02 0.0997 0.0257 0.0257 1.04
#> 6 Lower Estuary Dry Season 0.834 0.0704 0.0182 0.0182 0.8520.5 Three-Way Analysis of Variance (ANOVA)
The three_way_anova() function evaluates the
simultaneous influence of three categorical factors and their complete
interaction structure on a continuous aquatic response variable. In
limnology, marine biology, and ecotoxicology, three-way factorial
designs are commonly implemented to examine community or physiological
responses across spatial zones, hydrological seasons, and water column
depths simultaneously (e.g., Basin × Season × Depth).
The function automates factorial model fitting, calculation of main and interaction \(F\)-tests, derivation of Estimated Marginal Means across the complete three-way interaction lattice, post-hoc Tukey HSD letter grouping, and publication-ready faceted plotting.
0.5.1 Conceptual Framework and Statistical Formulation
The full three-way fixed-effects factorial model partitions the response variable \(y_{ijkl}\) into three main factor effects, three two-way interaction terms, and one three-way interaction term:
\[y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha\beta)_{ij} + (\alpha\gamma)_{ik} + (\beta\gamma)_{jk} + (\alpha\beta\gamma)_{ijk} + \varepsilon_{ijkl}\]
where: - \(y_{ijkl}\) represents the \(l\)-th observation under treatment combination \((i, j, k)\), - \(\mu\) represents the grand population mean, - \(\alpha_i\), \(\beta_j\), and \(\gamma_k\) represent the main effects of Factor 1 (e.g., Basin), Factor 2 (e.g., Season), and Factor 3 (e.g., Depth Layer), - \((\alpha\beta)_{ij}\), \((\alpha\gamma)_{ik}\), and \((\beta\gamma)_{jk}\) represent the two-way interaction effects, - \((\alpha\beta\gamma)_{ijk}\) represents the three-way interaction capturing higher-order non-linear dependencies across all factor levels, - \(\varepsilon_{ijkl} \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)\) is the random error term.
0.5.1.1 Variance Partitioning and Hypotheses
Total variation is partitioned into seven distinct explained components alongside the residual error:
\[\text{SS}_{\text{Total}} = \text{SS}_A + \text{SS}_B + \text{SS}_C + \text{SS}_{AB} + \text{SS}_{AC} + \text{SS}_{BC} + \text{SS}_{ABC} + \text{SS}_{\text{Residual}}\]
Each source of variation is tested against the residual mean square (\(\text{MS}_{\text{Residual}}\)). If the three-way interaction term \((\alpha\beta\gamma)_{ijk}\) is significant (\(p < 0.05\)), the response to Factor A depends simultaneously on the levels of both Factor B and Factor C. In this scenario, lower-order terms cannot be interpreted independently, and post-hoc contrasts are conducted across all joint treatment combinations.
0.5.2 In-Vignette Data Generation
We simulate a limnological biomonitoring dataset evaluating
Phytoplankton Chlorophyll-a concentrations (\(\mu\text{g/L}\)) across two Lake Basins
(North Basin, South Basin), two Seasons
(Pre-Monsoon, Post-Monsoon), and two Depth
Zones (Epilimnion, Hypolimnion):
set.seed(303)
n_reps <- 10
basins <- c("North Basin", "South Basin")
seasons <- c("Pre-Monsoon", "Post-Monsoon")
depths <- c("Epilimnion", "Hypolimnion")
design_grid <- expand.grid(
Basin = basins,
Season = seasons,
Depth = depths
)
chla_data <- do.call(rbind, lapply(seq_len(nrow(design_grid)), function(idx) {
row_info <- design_grid[idx, ]
# Baseline Chlorophyll-a dynamic
base_val <- 14.0
if (row_info$Basin == "South Basin") base_val <- base_val + 4.5
if (row_info$Season == "Post-Monsoon") base_val <- base_val + 6.0
if (row_info$Depth == "Hypolimnion") base_val <- base_val - 5.5
# Three-way synergistic interaction in productive shallow waters
if (row_info$Basin == "South Basin" && row_info$Season == "Post-Monsoon" && row_info$Depth == "Epilimnion") {
base_val <- base_val + 8.0
}
data.frame(
Basin = factor(rep(row_info$Basin, n_reps), levels = basins),
Season = factor(rep(row_info$Season, n_reps), levels = seasons),
Depth = factor(rep(row_info$Depth, n_reps), levels = depths),
Chla_ugL = rnorm(n_reps, mean = base_val, sd = 1.2)
)
}))
head(chla_data, 6)
#> Basin Season Depth Chla_ugL
#> 1 North Basin Pre-Monsoon Epilimnion 13.40367
#> 2 North Basin Pre-Monsoon Epilimnion 14.33461
#> 3 North Basin Pre-Monsoon Epilimnion 14.46136
#> 4 North Basin Pre-Monsoon Epilimnion 14.08230
#> 5 North Basin Pre-Monsoon Epilimnion 13.30615
#> 6 North Basin Pre-Monsoon Epilimnion 15.576670.5.3 Visualization Styles
The three_way_anova() function maps Factor 1 to the
x-axis, Factor 2 to fill/color groupings, and uses
facet_wrap to partition panels by Factor 3.
0.5.3.1 Style A: Faceted Boxplot with Jitter and Means (Default)
Displays medians, IQRs, observation jitters, diamond mean markers, and three-way interaction letters:
res_three_box <- three_way_anova(
data = chla_data,
factor1_var = "Basin",
factor2_var = "Season",
factor3_var = "Depth",
numeric_var = "Chla_ugL",
plot_type = "boxplot",
add_jitter = TRUE,
show_mean = TRUE,
mean_color = "darkred",
color_palette = "Set1"
)
res_three_box$Plot
0.5.3.2 Style B: Faceted Barplot with Standard Error Bars (Mean ± SE)
Displays aggregated group means with directional standard error bars across faceted sub-panels:
res_three_bar <- three_way_anova(
data = chla_data,
factor1_var = "Basin",
factor2_var = "Season",
factor3_var = "Depth",
numeric_var = "Chla_ugL",
plot_type = "barplot",
error_type = "se",
color_palette = "Set1"
)
res_three_bar$Plot
0.5.3.3 Style C: Minimalist Faceted Point-Range Geometry
Presents mean values and error bounds across faceted panels:
res_three_point <- three_way_anova(
data = chla_data,
factor1_var = "Basin",
factor2_var = "Season",
factor3_var = "Depth",
numeric_var = "Chla_ugL",
plot_type = "pointrange",
error_type = "se",
color_palette = "Dark2"
)
res_three_point$Plot
0.5.4 Model Inspection and Numerical Summaries
Extract the global three-way ANOVA summary, pairwise Tukey HSD differences, and Estimated Marginal Means:
# Global Three-Way ANOVA Table (Main effects, 2-way, and 3-way interactions)
res_three_box$ANOVA_Summary
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Basin 1 842.8 842.8 588.35 < 2e-16 ***
#> Season 1 1355.0 1355.0 945.99 < 2e-16 ***
#> Depth 1 1127.1 1127.1 786.88 < 2e-16 ***
#> Basin:Season 1 88.7 88.7 61.95 2.67e-11 ***
#> Basin:Depth 1 77.7 77.7 54.27 2.32e-10 ***
#> Season:Depth 1 55.3 55.3 38.58 3.03e-08 ***
#> Basin:Season:Depth 1 44.3 44.3 30.90 4.33e-07 ***
#> Residuals 72 103.1 1.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Post-hoc pairwise differences (Tukey HSD)
res_three_box$TukeyHSD
#> Tukey multiple comparisons of means
#> 95% family-wise confidence level
#>
#> Fit: stats::aov(formula = formula_str, data = clean_data)
#>
#> $Basin
#> diff lwr upr p adj
#> South Basin-North Basin 6.49138 5.95789 7.02487 0
#>
#> $Season
#> diff lwr upr p adj
#> Post-Monsoon-Pre-Monsoon 8.231143 7.697653 8.764633 0
#>
#> $Depth
#> diff lwr upr p adj
#> Hypolimnion-Epilimnion -7.507095 -8.040585 -6.973605 0
#>
#> $`Basin:Season`
#> diff lwr
#> South Basin:Pre-Monsoon-North Basin:Pre-Monsoon 4.384920 3.3895175
#> North Basin:Post-Monsoon-North Basin:Pre-Monsoon 6.124683 5.1292805
#> South Basin:Post-Monsoon-North Basin:Pre-Monsoon 14.722523 13.7271207
#> North Basin:Post-Monsoon-South Basin:Pre-Monsoon 1.739763 0.7443605
#> South Basin:Post-Monsoon-South Basin:Pre-Monsoon 10.337603 9.3422007
#> South Basin:Post-Monsoon-North Basin:Post-Monsoon 8.597840 7.6024377
#> upr p adj
#> South Basin:Pre-Monsoon-North Basin:Pre-Monsoon 5.380323 0.0000000
#> North Basin:Post-Monsoon-North Basin:Pre-Monsoon 7.120086 0.0000000
#> South Basin:Post-Monsoon-North Basin:Pre-Monsoon 15.717926 0.0000000
#> North Basin:Post-Monsoon-South Basin:Pre-Monsoon 2.735166 0.0001037
#> South Basin:Post-Monsoon-South Basin:Pre-Monsoon 11.333006 0.0000000
#> South Basin:Post-Monsoon-North Basin:Post-Monsoon 9.593243 0.0000000
#>
#> $`Basin:Depth`
#> diff lwr
#> South Basin:Epilimnion-North Basin:Epilimnion 8.462844 7.467441
#> North Basin:Hypolimnion-North Basin:Epilimnion -5.535631 -6.531034
#> South Basin:Hypolimnion-North Basin:Epilimnion -1.015715 -2.011117
#> North Basin:Hypolimnion-South Basin:Epilimnion -13.998475 -14.993877
#> South Basin:Hypolimnion-South Basin:Epilimnion -9.478558 -10.473961
#> South Basin:Hypolimnion-North Basin:Hypolimnion 4.519916 3.524514
#> upr p adj
#> South Basin:Epilimnion-North Basin:Epilimnion 9.45824641 0.0000000
#> North Basin:Hypolimnion-North Basin:Epilimnion -4.54022843 0.0000000
#> South Basin:Hypolimnion-North Basin:Epilimnion -0.02031205 0.0437494
#> North Basin:Hypolimnion-South Basin:Epilimnion -13.00307228 0.0000000
#> South Basin:Hypolimnion-South Basin:Epilimnion -8.48315590 0.0000000
#> South Basin:Hypolimnion-North Basin:Hypolimnion 5.51531894 0.0000000
#>
#> $`Season:Depth`
#> diff lwr
#> Post-Monsoon:Epilimnion-Pre-Monsoon:Epilimnion 9.8934143 8.8980118
#> Pre-Monsoon:Hypolimnion-Pre-Monsoon:Epilimnion -5.8448235 -6.8402261
#> Post-Monsoon:Hypolimnion-Pre-Monsoon:Epilimnion 0.7240484 -0.2713541
#> Pre-Monsoon:Hypolimnion-Post-Monsoon:Epilimnion -15.7382379 -16.7336404
#> Post-Monsoon:Hypolimnion-Post-Monsoon:Epilimnion -9.1693659 -10.1647685
#> Post-Monsoon:Hypolimnion-Pre-Monsoon:Hypolimnion 6.5688720 5.5734694
#> upr p adj
#> Post-Monsoon:Epilimnion-Pre-Monsoon:Epilimnion 10.888817 0.0000000
#> Pre-Monsoon:Hypolimnion-Pre-Monsoon:Epilimnion -4.849421 0.0000000
#> Post-Monsoon:Hypolimnion-Pre-Monsoon:Epilimnion 1.719451 0.2317211
#> Pre-Monsoon:Hypolimnion-Post-Monsoon:Epilimnion -14.742835 0.0000000
#> Post-Monsoon:Hypolimnion-Post-Monsoon:Epilimnion -8.173963 0.0000000
#> Post-Monsoon:Hypolimnion-Pre-Monsoon:Hypolimnion 7.564275 0.0000000
#>
#> $`Basin:Season:Depth`
#> diff
#> South Basin:Pre-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 4.8686400
#> North Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 6.2992104
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 18.3562582
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion -5.3611036
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion -1.4599035
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 0.5890520
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 5.7276847
#> North Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 1.4305705
#> South Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 13.4876182
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -10.2297436
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -6.3285435
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -4.2795879
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion 0.8590448
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Post-Monsoon:Epilimnion 12.0570477
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -11.6603140
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -7.7591139
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -5.7101584
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -0.5715257
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -23.7173618
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -19.8161617
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -17.7672061
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -12.6285734
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 3.9012001
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 5.9501557
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 11.0887883
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 2.0489556
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 7.1875883
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Hypolimnion 5.1386327
#> lwr
#> South Basin:Pre-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 3.1977257
#> North Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 4.6282962
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 16.6853440
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion -7.0320178
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion -3.1308177
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion -1.0818622
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 4.0567705
#> North Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion -0.2403437
#> South Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 11.8167040
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -11.9006578
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -7.9994577
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -5.9505021
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -0.8118694
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Post-Monsoon:Epilimnion 10.3861335
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -13.3312282
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -9.4300281
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -7.3810726
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -2.2424399
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -25.3882760
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -21.4870759
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -19.4381203
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -14.2994876
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 2.2302859
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 4.2792415
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 9.4178741
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 0.3780414
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 5.5166741
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Hypolimnion 3.4677185
#> upr
#> South Basin:Pre-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 6.5395542
#> North Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 7.9701246
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 20.0271724
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion -3.6901894
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 0.2110107
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 2.2599663
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 7.3985989
#> North Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 3.1014847
#> South Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 15.1585324
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -8.5588294
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -4.6576293
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion -2.6086737
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion 2.5299590
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Post-Monsoon:Epilimnion 13.7279620
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -9.9893998
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -6.0881997
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion -4.0392442
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion 1.0993885
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -22.0464476
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -18.1452475
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -16.0962919
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion -10.9576592
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 5.5721143
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 7.6210699
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 12.7597026
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 3.7198698
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 8.8585025
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Hypolimnion 6.8095469
#> p adj
#> South Basin:Pre-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 0.0000000
#> North Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 0.0000000
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Pre-Monsoon:Epilimnion 0.0000000
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 0.0000000
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 0.1311252
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 0.9546834
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Epilimnion 0.0000000
#> North Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 0.1480327
#> South Basin:Post-Monsoon:Epilimnion-South Basin:Pre-Monsoon:Epilimnion 0.0000000
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion 0.0000000
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion 0.0000000
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion 0.0000000
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Epilimnion 0.7457880
#> South Basin:Post-Monsoon:Epilimnion-North Basin:Post-Monsoon:Epilimnion 0.0000000
#> North Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion 0.0000000
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion 0.0000000
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion 0.0000000
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Epilimnion 0.9614041
#> North Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion 0.0000000
#> South Basin:Pre-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion 0.0000000
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion 0.0000000
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Post-Monsoon:Epilimnion 0.0000000
#> South Basin:Pre-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 0.0000000
#> North Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 0.0000000
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Pre-Monsoon:Hypolimnion 0.0000000
#> North Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 0.0063183
#> South Basin:Post-Monsoon:Hypolimnion-South Basin:Pre-Monsoon:Hypolimnion 0.0000000
#> South Basin:Post-Monsoon:Hypolimnion-North Basin:Post-Monsoon:Hypolimnion 0.0000000
# Estimated Marginal Means across the 3-way interaction grid
res_three_box$EMMs
#> Basin Season Depth emmean SE df lower.CL upper.CL
#> North Basin Pre-Monsoon Epilimnion 13.92 0.378 72 13.17 14.68
#> South Basin Pre-Monsoon Epilimnion 18.79 0.378 72 18.04 19.55
#> North Basin Post-Monsoon Epilimnion 20.22 0.378 72 19.47 20.98
#> South Basin Post-Monsoon Epilimnion 32.28 0.378 72 31.52 33.03
#> North Basin Pre-Monsoon Hypolimnion 8.56 0.378 72 7.81 9.32
#> South Basin Pre-Monsoon Hypolimnion 12.46 0.378 72 11.71 13.22
#> North Basin Post-Monsoon Hypolimnion 14.51 0.378 72 13.76 15.27
#> South Basin Post-Monsoon Hypolimnion 19.65 0.378 72 18.90 20.41
#>
#> Confidence level used: 0.95
# Aggregated summary statistics
res_three_box$Summary_Table
#> # A tibble: 8 × 8
#> Basin Season Depth mean_val sd_val se_val error_margin upper_limit
#> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 North Basin Pre-Monsoon Epil… 13.9 0.898 0.284 0.284 14.2
#> 2 North Basin Pre-Monsoon Hypo… 8.56 1.47 0.466 0.466 9.03
#> 3 North Basin Post-Monsoon Epil… 20.2 1.21 0.384 0.384 20.6
#> 4 North Basin Post-Monsoon Hypo… 14.5 1.11 0.352 0.352 14.9
#> 5 South Basin Pre-Monsoon Epil… 18.8 0.964 0.305 0.305 19.1
#> 6 South Basin Pre-Monsoon Hypo… 12.5 1.52 0.481 0.481 12.9
#> 7 South Basin Post-Monsoon Epil… 32.3 0.910 0.288 0.288 32.6
#> 8 South Basin Post-Monsoon Hypo… 19.7 1.31 0.413 0.413 20.10.6 Analysis of Covariance (ANCOVA)
The ancova_analysis() function combines linear
regression and Analysis of Variance to test for differences in a
continuous response variable across categorical treatment groups while
statistically controlling for a continuous nuisance variable
(covariate). In fisheries stock assessment and aquatic ecology, ANCOVA
is routinely used to compare fish weight across sampling locations while
controlling for body length, or to compare metabolic excretion rates
across diet formulations while controlling for individual biomass.
The function handles missing value removal, fits an additive general linear model, extracts the omnibus ANCOVA table, computes Estimated Marginal Means (EMMs) adjusted to the mean value of the covariate, and renders publication-quality scatter plots with parallel regression trend lines.
0.6.1 Conceptual Framework and Statistical Formulation
The additive fixed-effects ANCOVA model represents the response variable \(y_{ij}\) as:
\[y_{ij} = \mu + \alpha_i + \beta (x_{ij} - \bar{x}) + \varepsilon_{ij}\]
where: - \(y_{ij}\) is the response observed for the \(j\)-th individual in group \(i\) (\(i = 1, \dots, k\); \(j = 1, \dots, n_i\)), - \(\mu\) is the overall adjusted grand mean, - \(\alpha_i\) is the categorical treatment effect of group \(i\), subject to the constraint \(\sum n_i \alpha_i = 0\), - \(\beta\) is the common pooled regression slope relating the continuous covariate \(x\) to response \(y\), - \(x_{ij}\) is the covariate value for individual \(j\) in group \(i\), centered around the grand covariate mean \(\bar{x}\), - \(\varepsilon_{ij} \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)\) is the normally distributed random error term.
0.6.1.1 Homogeneity of Slopes and Adjusted Means
A core assumption of additive ANCOVA is the homogeneity of regression slopes, meaning the relationship between the covariate (\(x\)) and the response (\(y\)) does not vary across factor levels (\(\beta_1 = \beta_2 = \dots = \beta_k = \beta\)).
When this assumption holds, the categorical factor effect evaluates whether the elevations (y-intercepts) of the parallel regression lines differ significantly:
\[H_0: \alpha_1 = \alpha_2 = \dots = \alpha_k = 0\] \[H_1: \text{At least one group elevation differs.}\]
Because raw group means are confounded by variations in the
covariate, ancova_analysis() computes Estimated Marginal
Means (EMMs). The adjusted mean for group \(i\) (\(\hat{\mu}_i\)) is evaluated at the grand
mean of the covariate (\(\bar{x}_{\bullet\bullet}\)):
\[\hat{\mu}_i = \bar{y}_{i\bullet} - \hat{\beta}(\bar{x}_{i\bullet} - \bar{x}_{\bullet\bullet})\]
0.6.2 In-Vignette Data Generation
We simulate a fisheries length-weight dataset evaluating somatic body
weight (\(\text{Weight}\), \(\text{g}\)) across three distinct aquatic
habitats (Estuary, River, Marine)
while controlling for continuous total body length (\(\text{Length}\), \(\text{cm}\)):
set.seed(2026)
n_per_habitat <- 35
habitats <- c("Estuary", "River", "Marine")
fish_growth <- do.call(rbind, lapply(habitats, function(hab) {
# Generate continuous covariate: fish total length (cm)
len <- runif(n_per_habitat, min = 12.0, max = 38.0)
# Group elevation shifts (intercept differences)
elevation_shift <- switch(hab,
"Estuary" = -12.0,
"River" = -4.0,
"Marine" = 8.0
)
# Additive response: common slope (beta = 3.25) + habitat elevation + noise
wt <- 3.25 * len + elevation_shift + rnorm(n_per_habitat, mean = 0, sd = 4.2)
data.frame(
Habitat = factor(rep(hab, n_per_habitat), levels = habitats),
Length = len,
Weight = wt
)
}))
head(fish_growth, 6)
#> Habitat Length Weight
#> 1 Estuary 30.16551 85.25238
#> 2 Estuary 26.46979 73.78666
#> 3 Estuary 15.64364 38.28039
#> 4 Estuary 19.42881 55.97057
#> 5 Estuary 26.43959 73.29607
#> 6 Estuary 12.65341 21.492770.6.3 Visualization Styles
ancova_analysis() plots the empirical observations
alongside model-fitted regression lines. This preserves the parallel
slopes of the additive model, avoiding the misleading non-parallel lines
produced by unconstrained smoothers.
0.6.3.1 Style A: Single-Panel Parallel Slopes Comparison (Default)
Presents all groups on a shared coordinate grid
(facet = FALSE), allowing direct visual comparison of
vertical line separations (elevation differences):
res_ancova <- ancova_analysis(
data = fish_growth,
response_var = "Weight",
factor_var = "Habitat",
covariate_var = "Length",
facet = FALSE,
color_palette = "Dark2"
)
# Render publication plot
res_ancova$Plot
0.6.3.2 Style B: Faceted Panel View
When working with dense datasets or multiple overlapping factor
levels, setting facet = TRUE separates each group into its
own panel using facet_wrap, while keeping the coordinate
scales fixed for consistent comparison:
res_ancova_facet <- ancova_analysis(
data = fish_growth,
response_var = "Weight",
factor_var = "Habitat",
covariate_var = "Length",
facet = TRUE,
color_palette = "Dark2"
)
res_ancova_facet$Plot
0.6.4 Model Inspection and Numerical Summaries
The returned list object provides the ANOVA table testing the covariate and categorical factor effects, as well as the covariate-adjusted Estimated Marginal Means:
# 1. ANCOVA Table (Testing covariate and factor significance)
res_ancova$ANCOVA_Table
#> Analysis of Variance Table
#>
#> Response: Weight
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Length 1 69814 69814 4111.69 < 2.2e-16 ***
#> Habitat 2 6473 3236 190.61 < 2.2e-16 ***
#> Residuals 101 1715 17
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2. Estimated Marginal Means (EMMs adjusted to mean Length)
res_ancova$Adjusted_Means
#> Habitat emmean SE df lower.CL upper.CL
#> Estuary 61.6 0.699 101 60.2 63.0
#> River 69.2 0.698 101 67.8 70.6
#> Marine 81.0 0.703 101 79.6 82.3
#>
#> Confidence level used: 0.950.7 Multivariate Analysis of Variance (MANOVA)
The manova_analysis() function extends analysis of
variance to multi-response hypotheses, testing whether mean vectors of
several correlated continuous variables differ across levels of a
categorical grouping variable. In fisheries biology and aquatic ecology,
MANOVA is commonly applied in morphometric stock discrimination (e.g.,
comparing head length, body depth, and caudal peduncle length
simultaneously across river populations) or multi-metric water quality
assessments (e.g., evaluating simultaneous variations in nitrates,
phosphates, and ammonium across wetland sites).
The function conducts defensive data verification, fits the multivariate general linear model, evaluates omnibus significance using Pillai’s trace criterion, derives univariate ANOVA decompositions, and provides two distinct visualization geometries: faceted univariate trait boxplots and Canonical Discriminant Analysis (LDA) scatter plots with 95% confidence ellipses.
0.7.1 Conceptual Framework and Statistical Formulation
Let \(\mathbf{y}_{ij} = [y_{ij1}, y_{ij2}, \dots, y_{ijp}]^T\) represent a \(p\)-dimensional vector of continuous response measurements recorded for the \(j\)-th specimen in treatment group \(i\) (\(i = 1, \dots, k\); \(j = 1, \dots, n_i\)). The one-way multivariate linear model is defined as:
\[\mathbf{y}_{ij} = \boldsymbol{\mu} + \boldsymbol{\alpha}_i + \boldsymbol{\varepsilon}_{ij}\]
where: - \(\boldsymbol{\mu} = [\mu_1, \dots, \mu_p]^T\) represents the grand centroid vector across all \(p\) traits, - \(\boldsymbol{\alpha}_i = [\alpha_{i1}, \dots, \alpha_{ip}]^T\) represents the differential group effect vector for group \(i\), - \(\boldsymbol{\varepsilon}_{ij} \overset{\text{i.i.d.}}{\sim} \mathcal{N}_p(\mathbf{0}, \boldsymbol{\Sigma})\) represents the multivariate normal random error vector with variance-covariance matrix \(\boldsymbol{\Sigma}\).
0.7.1.1 Hypothesis Testing and Test Statistics
The omnibus multivariate null hypothesis states that group centroids are identical across all \(p\) dimensions:
\[H_0: \boldsymbol{\mu}_1 = \boldsymbol{\mu}_2 = \dots = \boldsymbol{\mu}_k\] \[H_1: \text{At least two group centroid vectors differ.}\]
The total variance-covariance matrix is partitioned into the Between-Group Sum of Squares and Cross-Products matrix (\(\mathbf{H}\)) and the Within-Group (Error) Sum of Squares and Cross-Products matrix (\(\mathbf{E}\)):
\[\mathbf{H} = \sum_{i=1}^k n_i (\bar{\mathbf{y}}_{i\bullet} - \bar{\mathbf{y}}_{\bullet\bullet})(\bar{\mathbf{y}}_{i\bullet} - \bar{\mathbf{y}}_{\bullet\bullet})^T\]
\[\mathbf{E} = \sum_{i=1}^k \sum_{j=1}^{n_i} (\mathbf{y}_{ij} - \bar{\mathbf{y}}_{i\bullet})(\mathbf{y}_{ij} - \bar{\mathbf{y}}_{i\bullet})^T\]
While several test criteria exist (Wilks’ Lambda, Hotelling-Lawley
trace, Roy’s largest root), manova_analysis() defaults to
Pillai’s Trace (\(V\)) due to its
robustness against departures from multivariate normality and moderate
covariance heterogeneity:
\[V = \text{tr}\left[\mathbf{H}(\mathbf{H} + \mathbf{E})^{-1}\right] = \sum_{m=1}^s \frac{\lambda_m}{1 + \lambda_m}\]
where \(\lambda_m\) represents the eigenvalues of \(\mathbf{E}^{-1}\mathbf{H}\), and \(s = \min(p, k - 1)\). The statistic is converted to an approximate \(F\)-distribution to establish statistical significance.
0.7.2 In-Vignette Data Generation
We simulate a morphometric stock identification survey measuring four
anatomical truss traits (Total Length TL, Head Length
HL, Body Depth BD, and Caudal Width
CW in cm) across three wild fish populations
(River_A, River_B, and
River_C):
set.seed(505)
n_fish <- 25
populations <- c("River_A", "River_B", "River_C")
# Simulate multivariate normal morphometric data
morpho_data <- do.call(rbind, lapply(populations, function(pop) {
# Baseline centroid shifts per population
shift <- switch(pop,
"River_A" = c(TL = 24.0, HL = 5.2, BD = 6.5, CW = 2.1),
"River_B" = c(TL = 28.5, HL = 6.4, BD = 7.8, CW = 2.6),
"River_C" = c(TL = 21.0, HL = 4.8, BD = 5.9, CW = 1.9)
)
data.frame(
Population = factor(rep(pop, n_fish), levels = populations),
TL = rnorm(n_fish, mean = shift["TL"], sd = 1.5),
HL = rnorm(n_fish, mean = shift["HL"], sd = 0.4),
BD = rnorm(n_fish, mean = shift["BD"], sd = 0.5),
CW = rnorm(n_fish, mean = shift["CW"], sd = 0.2)
)
}))
head(morpho_data, 6)
#> Population TL HL BD CW
#> 1 River_A 22.31822 5.435375 5.605484 1.812157
#> 2 River_A 22.07691 5.041592 6.145151 1.972692
#> 3 River_A 20.94099 5.282860 5.595230 1.829734
#> 4 River_A 22.59340 5.254322 6.134524 2.257909
#> 5 River_A 23.23476 5.786759 6.319575 2.392531
#> 6 River_A 23.38073 4.793155 5.212262 1.9675720.7.3 Visualization Styles
The manova_analysis() function provides two
complementary visual perspectives: faceted univariate trait
distributions and multivariate canonical ordination.
0.7.3.1 Style A: Faceted Univariate Trait Boxplots (Default)
Decomposes the multivariate space into separate panels for each continuous trait, showing median values, interquartile ranges, and individual observation distributions:
res_manova_box <- manova_analysis(
data = morpho_data,
response_vars = c("TL", "HL", "BD", "CW"),
factor_var = "Population",
plot_type = "boxplot",
color_palette = "Set2"
)
# Render faceted traits
res_manova_box$Plot
0.7.3.2 Style B: Canonical Discriminant Analysis (LDA) Ordination
Projects the multi-trait measurements into canonical discriminant
space (Linear Discriminant axes LD1 and LD2),
displaying multivariate centroid separation and 95% group confidence
ellipses:
res_manova_lda <- manova_analysis(
data = morpho_data,
response_vars = c("TL", "HL", "BD", "CW"),
factor_var = "Population",
plot_type = "lda",
color_palette = "Set2"
)
# Render canonical discriminant ordination
res_manova_lda$Plot
0.7.4 Model Inspection and Numerical Summaries
The returned list delivers both the global multivariate Pillai’s trace summary and the individual univariate ANOVA breakdowns:
# 1. Omnibus MANOVA Summary (Pillai's Trace, F-approximation, p-value)
res_manova_box$MANOVA_Summary
#> Df Pillai approx F num Df den Df Pr(>F)
#> Population 2 0.94978 15.826 8 140 < 2.2e-16 ***
#> Residuals 72
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2. Univariate ANOVA Breakdowns for Each Morphometric Trait
res_manova_box$Univariate_ANOVAs
#> Response TL :
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Population 2 755.7 377.85 162.71 < 2.2e-16 ***
#> Residuals 72 167.2 2.32
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Response HL :
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Population 2 26.439 13.2195 93.691 < 2.2e-16 ***
#> Residuals 72 10.159 0.1411
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Response BD :
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Population 2 47.020 23.5101 69.044 < 2.2e-16 ***
#> Residuals 72 24.517 0.3405
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Response CW :
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Population 2 7.6740 3.837 78.335 < 2.2e-16 ***
#> Residuals 72 3.5267 0.049
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 10.8 Parametric & Non-Linear Regression Modeling
The regression_analysis() function models functional
dependencies between continuous aquatic parameters—including
morphometric fish growth, allometric length-weight relationships (\(W = a \cdot L^b\)), nutrient saturation
kinetics, chemical calibration lines, and thermal performance
curves.
The function provides: 1. Dynamic Equation Generation: Automatically extracts model coefficients, evaluates statistical significance (\(R^2\) and \(p\)-value), and prints the parsed algebraic formula directly on the plot canvas. 2. Dual Uncertainty Bands: Toggles between Confidence Intervals (measuring precision of the expected mean response \(\mu_{y\mid x}\)) and Prediction Intervals (capturing individual biological variation of future single specimens). 3. Comprehensive Functional Forms: Fits linear (\(y = \beta_0 + \beta_1 x\)), allometric power (\(y = a \cdot x^b\)), exponential (\(y = a \cdot e^{kx}\)), logarithmic (\(y = \beta_0 + \beta_1 \ln(x)\)), quadratic polynomial (\(y = \beta_0 + \beta_1 x + \beta_2 x^2\)), or cubic polynomial (\(y = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3\)) models across single-panel and multi-panel faceted designs.
0.8.1 Conceptual Framework & Confidence vs. Prediction Bounds
For a linear or linearized model \(Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i\) with residual variance \(s^2\):
0.8.1.1 Confidence Interval for the Mean Response
Measures uncertainty in estimating the true conditional population mean \(\mathbb{E}(Y \mid x_0)\) at a specified predictor value \(x_0\):
\[\text{CI}_{95\%}(x_0) = \hat{y}_0 \pm t_{n-2, \, 0.975} \times s \sqrt{\frac{1}{n} + \frac{(x_0 - \bar{x})^2}{\sum_{i=1}^n (x_i - \bar{x})^2}}\]
As sample size increases (\(n \to \infty\)), the confidence band contracts toward zero width around the fitted trajectory.
0.8.1.2 Prediction Interval for Individual Observations
Predicts the boundary within which a single newly sampled organism \(Y_{\text{new}}\) is expected to fall:
\[\text{PI}_{95\%}(x_0) = \hat{y}_0 \pm t_{n-2, \, 0.975} \times s \sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar{x})^2}{\sum_{i=1}^n (x_i - \bar{x})^2}}\]
The additional unit term (\(1\)) inside the radical accounts for irreducible biological variance (\(\sigma^2\)), making prediction bands substantially wider than confidence intervals.
0.8.2 In-Vignette Data Generation
We simulate a riverine fisheries biometric dataset tracking Total
Length (cm) and Body Weight (g) across 75 specimens collected from three
river basins (River A, River B,
River C):
set.seed(606)
n_fish <- 25
rivers <- c("River A", "River B", "River C")
fish_allometry <- do.call(rbind, lapply(rivers, function(rv) {
len <- runif(n_fish, min = 12, max = 38)
# Allometric power parameters: W = a * L^b
a_param <- switch(rv, "River A" = 0.012, "River B" = 0.009, "River C" = 0.015)
b_param <- switch(rv, "River A" = 3.05, "River B" = 3.18, "River C" = 2.92)
wt <- a_param * (len^b_param) * exp(rnorm(n_fish, mean = 0, sd = 0.12))
data.frame(
River = rv,
Length = len,
Weight = wt,
stringsAsFactors = FALSE
)
}))
head(fish_allometry, 6)
#> River Length Weight
#> 1 River A 26.15128 276.34037
#> 2 River A 33.53511 573.18190
#> 3 River A 16.57117 56.49807
#> 4 River A 32.41410 516.93445
#> 5 River A 14.74483 40.48794
#> 6 River A 37.76574 750.740980.8.3 Visualization Styles Across Functional Models
regression_analysis() handles linear, power-law,
exponential, logarithmic, and higher-order polynomial relationships.
0.8.3.1 Style A: Allometric Power-Law Model (\(W = a \cdot L^b\)) with 95% Prediction Interval
Fits an allometric power curve, calculates wide 95% prediction intervals to accommodate individual specimen growth scatter, and plots the power formula on the canvas:
res_power_reg <- regression_analysis(
data = fish_allometry,
x_var = "Length",
y_var = "Weight",
fit_type = "power",
interval = "prediction",
show_equation = TRUE,
title = "Allometric Fish Length-Weight Relationship (95% PI)"
)
# Render allometric power plot
res_power_reg$Plot
0.8.3.2 Style B: Exponential Growth/Decay (\(y = a \cdot e^{kx}\)) with 95% Confidence Interval
Models exponential trajectories (e.g., mortality catch curves, light attenuation) with confidence intervals around the expected conditional mean:
res_exp_reg <- regression_analysis(
data = fish_allometry,
x_var = "Length",
y_var = "Weight",
fit_type = "exponential",
interval = "confidence",
show_equation = TRUE,
title = "Exponential Length-Weight Curve (95% CI)"
)
# Render exponential curve
res_exp_reg$Plot
0.8.3.3 Style C: Faceted Quadratic Polynomial (\(y = \beta_0 + \beta_1 x + \beta_2 x^2\))
Fits 2nd-order polynomial curves across river systems with
facet = TRUE, placing panel-specific equations and \(R^2\) values in each sub-panel:
res_poly_reg <- regression_analysis(
data = fish_allometry,
x_var = "Length",
y_var = "Weight",
fit_type = "polynomial",
group_var = "River",
facet = TRUE,
interval = "confidence",
show_equation = TRUE,
color_palette = "Dark2",
title = "Faceted Quadratic Growth Trajectories by Basin"
)
# Render faceted quadratic curves
res_poly_reg$Plot
0.8.3.4 Style D: Faceted Cubic Polynomial (\(y = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3\))
Accommodates complex, multi-inflection trajectories (e.g., seasonal thermal tolerance limits or multi-stage cohort growth) using a 3rd-order polynomial:
res_cubic_reg <- regression_analysis(
data = fish_allometry,
x_var = "Length",
y_var = "Weight",
fit_type = "cubic",
group_var = "River",
facet = TRUE,
interval = "confidence",
show_equation = TRUE,
color_palette = "Set1",
title = "Faceted Cubic Growth Polynomials Across Basins"
)
# Render faceted cubic curves
res_cubic_reg$Plot
0.8.4 Model Inspection & Diagnostics
The returned list object provides access to algebraic formulas,
goodness-of-fit statistics (\(R^2\),
\(F\)-test \(p\)-values), fitted lm model
objects, and evaluation grids:
# 1. Tidy summary of model equations, R2 values, and p-values
res_poly_reg$Equation_Table
#> Group Fit_Type Equation R_Squared P_Value
#> 1 River A polynomial y == '117.91' - '19.05'*x + '0.948'*x^2 0.9766 0
#> 2 River B polynomial y == '163.88' - '24.71'*x + '1.134'*x^2 0.9846 0
#> 3 River C polynomial y == '66.86' - '11.83'*x + '0.671'*x^2 0.9795 0
#> Label
#> 1 atop(y == '117.91' - '19.05'*x + '0.948'*x^2, paste(R^2 == '0.98', ', ', p < 0.001))
#> 2 atop(y == '163.88' - '24.71'*x + '1.134'*x^2, paste(R^2 == '0.98', ', ', p < 0.001))
#> 3 atop(y == '66.86' - '11.83'*x + '0.671'*x^2, paste(R^2 == '0.98', ', ', p < 0.001))
# 2. Inspect underlying linear model coefficients for River A
summary(res_poly_reg$Models[["River A"]])
#>
#> Call:
#> stats::lm(formula = y ~ x + x2, data = df_fit)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -89.230 -13.681 1.514 8.028 60.074
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 117.9122 83.6118 1.410 0.1725
#> x -19.0516 7.3204 -2.603 0.0163 *
#> x2 0.9484 0.1476 6.424 1.83e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 36.18 on 22 degrees of freedom
#> Multiple R-squared: 0.9766, Adjusted R-squared: 0.9745
#> F-statistic: 459 on 2 and 22 DF, p-value: < 2.2e-16
# 3. Preview calculated prediction grid with 95% confidence intervals
head(res_poly_reg$Predictions, 6)
#> x y_fit y_lwr y_upr Group
#> 1 12.70514 28.95552 -7.0132498 64.92429 River A
#> 2 12.91573 30.06073 -4.6373997 64.75886 River A
#> 3 13.12632 31.25007 -2.2206115 64.72076 River A
#> 4 13.33692 32.52354 0.2359541 64.81112 River A
#> 5 13.54751 33.88113 2.7310749 65.03118 River A
#> 6 13.75810 35.32284 5.2634742 65.38221 River A0.9 Generalized Linear Models (GLMs)
The glm_analysis() function handles non-normal aquatic
survey data where standard ordinary least squares (OLS) assumptions are
violated. In fisheries surveys, biodiversity assessments, and habitat
modeling, empirical data frequently feature discrete integer counts with
positive skew (e.g., fish catch per unit effort, abundance of invasive
species) or binary presence/absence states (e.g., detection of
endangered fauna or disease occurrence).
The function fits models across Poisson, Quasipoisson (for overdispersed counts), or Binomial error distributions using link functions (\(g(\cdot)\)). It extracts Wald parameter estimates, derives multiplicative effect metrics—Incident Rate Ratios (IRR) for count families or Odds Ratios (OR) for binary families—computes 95% Wald confidence bounds, and generates a publication-ready forest plot.
0.9.1 Conceptual Framework and Statistical Formulation
A Generalized Linear Model consists of three components: 1. Random Component: An exponential dispersion family distribution for the response variable \(Y_i\). 2. Systematic Component: A linear predictor combination of \(p\) explanatory variables: \[\eta_i = \beta_0 + \sum_{j=1}^p \beta_j x_{ij}\] 3. Link Function: A monotonic function \(g(\cdot)\) mapping the expected mean \(\mu_i = \mathbb{E}(Y_i)\) to the linear predictor: \[\eta_i = g(\mu_i)\]
0.9.1.1 Poisson and Quasipoisson Regression (Count Dynamics)
For count data (e.g., CPUE or species richness), the log link ensures strictly positive predicted yields:
\[\ln(\mu_i) = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}\]
Exponentiating coefficients yields the Incident Rate Ratio (IRR):
\[\text{IRR}_j = \exp(\beta_j)\]
An \(\text{IRR} > 1\) denotes a
multiplicative increase in the expected count per unit increase in \(x_j\), while an \(\text{IRR} < 1\) represents a
proportional reduction. When the variance exceeds the mean (\(\text{Var}(Y) = \phi \mu\) with dispersion
parameter \(\phi > 1\)), the
quasipoisson family inflates parameter standard errors by
\(\sqrt{\phi}\) to account for
overdispersion without altering point estimates.
0.9.1.2 Binomial Logistic Regression (Presence / Absence)
For binary occurrence states \(Y_i \in \{0, 1\}\), the logit link maps probability \(p_i = P(Y_i = 1)\) to continuous space:
\[\text{logit}(p_i) = \ln\left(\frac{p_i}{1 - p_i}\right) = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}\]
Exponentiating estimates yields the Odds Ratio (OR):
\[\text{OR}_j = \exp(\beta_j)\]
The 95% Wald confidence intervals for effect sizes are calculated as:
\[\text{CI}_{95\%} = \exp(\hat{\beta}_j \pm 1.96 \times \text{SE}(\hat{\beta}_j))\]
If the confidence interval spans \(1.0\), the predictor does not exert a statistically detectable multiplicative effect at \(\alpha = 0.05\).
0.9.2 In-Vignette Data Generation
We simulate an estuarine fish survey monitoring the abundance (count
per haul) of a juvenile target species across 120 sampling trawls. The
response count is modeled as a function of Water Temperature (°C),
Salinity (PSU), Dissolved Oxygen (mg/L), and Substrate Type
(Mud vs. Sand):
set.seed(707)
n_trawls <- 120
temp_c <- runif(n_trawls, min = 18, max = 30)
sal_psu <- runif(n_trawls, min = 5, max = 32)
do_mgl <- runif(n_trawls, min = 3.5, max = 9.0)
substr <- sample(c("Mud", "Sand"), size = n_trawls, replace = TRUE)
# Multiplicative linear predictor
log_mu <- -1.8 +
0.08 * temp_c -
0.04 * sal_psu +
0.22 * do_mgl +
ifelse(substr == "Sand", 0.35, 0.0)
# Simulate Poisson counts
catch_counts <- rpois(n_trawls, lambda = exp(log_mu))
trawl_survey <- data.frame(
Abundance = catch_counts,
Temperature = temp_c,
Salinity = sal_psu,
DO = do_mgl,
Substrate = factor(substr, levels = c("Mud", "Sand"))
)
head(trawl_survey, 6)
#> Abundance Temperature Salinity DO Substrate
#> 1 1 19.19934 31.872330 4.420170 Mud
#> 2 3 23.25641 7.830847 5.040356 Sand
#> 3 1 21.22099 16.444491 5.220697 Sand
#> 4 2 24.38437 5.188938 4.339394 Sand
#> 5 4 27.09255 16.775980 5.855344 Mud
#> 6 0 26.38856 12.958779 5.937050 Mud0.9.3 Visualization Styles
The forest plot produced by glm_analysis() centers
predictor effects relative to the reference null line at \(\text{Effect Size} = 1.0\).
0.9.3.1 Style A: Count Regression Forest Plot (Poisson Model)
Evaluates the Incident Rate Ratio (IRR) across environmental drivers, with horizontal error bars showing 95% confidence intervals:
res_glm_count <- glm_analysis(
data = trawl_survey,
response_var = "Abundance",
predictor_vars = c("Temperature", "Salinity", "DO", "Substrate"),
family_type = "poisson",
color_palette = "Dark2"
)
# Render forest plot
res_glm_count$Plot
0.9.3.2 Style B: Binary Occurrence Forest Plot (Binomial Model)
When the response variable is converted into occurrence
(presence/absence), family_type = "binomial" generates Odds
Ratios (OR):
# Create binary presence/absence column
trawl_survey$Presence <- ifelse(trawl_survey$Abundance > 0, 1, 0)
res_glm_binary <- glm_analysis(
data = trawl_survey,
response_var = "Presence",
predictor_vars = c("Temperature", "Salinity", "DO", "Substrate"),
family_type = "binomial",
color_palette = "Set1"
)
res_glm_binary$Plot
0.9.4 Model Inspection and Numerical Summaries
The returned list provides the GLM model summary and a coefficient table with calculated effect sizes:
# 1. Global GLM Summary (Deviance, Dispersion, Estimates, z-values, p-values)
res_glm_count$Model_Summary
#>
#> Call:
#> stats::glm(formula = formula_str, family = fam_obj, data = clean_data)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -1.896996 0.476145 -3.984 6.77e-05 ***
#> Temperature 0.073725 0.015557 4.739 2.15e-06 ***
#> Salinity -0.033417 0.007253 -4.607 4.08e-06 ***
#> DO 0.230767 0.035161 6.563 5.27e-11 ***
#> SubstrateSand 0.374581 0.112705 3.324 0.000889 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for poisson family taken to be 1)
#>
#> Null deviance: 225.54 on 119 degrees of freedom
#> Residual deviance: 116.91 on 115 degrees of freedom
#> AIC: 432.46
#>
#> Number of Fisher Scoring iterations: 5
# 2. Tidy Coefficient Table with Exponentiated Effect Sizes (IRR) and 95% CIs
res_glm_count$Coefficient_Table
#> Estimate Std_Error z_value p_value Term
#> Temperature 0.07372515 0.015556870 4.739074 2.146976e-06 Temperature
#> Salinity -0.03341746 0.007253455 -4.607110 4.083048e-06 Salinity
#> DO 0.23076685 0.035161128 6.563124 5.269193e-11 DO
#> SubstrateSand 0.37458124 0.112705395 3.323543 8.888167e-04 SubstrateSand
#> Effect_Size Lower_CI Upper_CI
#> Temperature 1.0765109 1.0441819 1.1098408
#> Salinity 0.9671347 0.9534825 0.9809825
#> DO 1.2595655 1.1756851 1.3494305
#> SubstrateSand 1.4543823 1.1661161 1.81390830.10 Non-metric Multidimensional Scaling (NMDS)
The nmds_analysis() function performs unconstrained
multivariate ordination for community composition datasets (e.g.,
phytoplankton, zooplankton, macroinvertebrates, and fish assemblages).
Unlike eigenvalue-based ordination techniques (e.g., PCA), NMDS avoids
assuming linear relationships between species abundances and underlying
environmental gradients. Instead, it relies on rank-order
dissimilarities, making it resilient to zero-inflated ecological count
matrices.
The function supports input ingestion (data frames, matrices, lists),
implements standard ecological transformations (e.g., Hellinger,
Wisconsin double standardization), runs iterative ordination using
vegan::metaMDS, monitors convergence stress, overlays 95%
normal confidence ellipses, and projects significant environmental
vectors using vegan::envfit.
0.10.1 Conceptual Framework and Mathematical Formulation
Let \(\mathbf{Y}_{n \times p}\) denote a community matrix of \(n\) sampling sites and \(p\) taxa. Pairwise ecological dissimilarity between site \(j\) and site \(k\) is calculated using an appropriate index (typically Bray-Curtis dissimilarity for abundance data):
\[d_{jk} = \frac{\sum_{i=1}^p |y_{ij} - y_{ik}|}{\sum_{i=1}^p (y_{ij} + y_{ik})}\]
NMDS seeks a low-dimensional configuration (\(k = 2\) or \(3\) dimensions) of site coordinates \(\mathbf{X}_{n \times k}\) such that the Euclidean distances between points in ordination space (\(\hat{d}_{jk} = \Vert{}\mathbf{x}_j - \mathbf{x}_k\Vert{}\)) preserve the monotonic rank order of the original dissimilarities (\(d_{jk}\)):
\[d_{jk} \approx f(\hat{d}_{jk})\]
where \(f\) is a monotonic step or regression function fitted using isotonic regression.
0.10.1.1 Stress Evaluation and Quality Criteria
The goodness-of-fit of the configuration is measured by Kruskal’s Stress (\(S\)):
\[S = \sqrt{\frac{\sum_{j < k} (d_{jk} - \hat{d}_{jk})^2}{\sum_{j < k} d_{jk}^2}}\]
Stress measures the rank-order discrepancy between original ecological dissimilarities and low-dimensional mapped distances: - \(S < 0.05\): Excellent representation with minimal distortion. - \(0.05 \le S < 0.10\): Good ordination with reliable configuration. - \(0.10 \le S \le 0.20\): Usable ordination; 2D configuration provides an adequate representation, but higher-order details should be interpreted with care. - \(S > 0.20\): High stress; results should be interpreted cautiously, and increasing ordination dimensionality (\(k = 3\)) may be warranted.
0.10.1.2 Pre-Ordination Hellinger Transformation
To alleviate the “horseshoe” or “arch” effect common in raw Euclidean
distance models, transform = "hellinger" applies the
Hellinger metric:
\[y'_{ij} = \sqrt{\frac{y_{ij}}{\sum_{m=1}^p y_{mj}}}\]
This downweights dominant taxa and eliminates the zero-inflation paradox without dropping absent species.
0.10.2 In-Vignette Data Generation
We simulate a river basin macroinvertebrate community survey across
45 sampling stations distributed among three Habitat types
(Riffle, Run, Pool). Seven
diagnostic taxa counts are generated, along with three continuous water
quality parameters (Flow Velocity, Dissolved Oxygen, and Water
Temperature):
set.seed(808)
n_sites <- 15
habitats <- c("Riffle", "Run", "Pool")
# Macroinvertebrate community matrix (7 diagnostic taxa)
comm_list <- lapply(habitats, function(hab) {
# Baseline taxon abundance profiles
means <- switch(hab,
"Riffle" = c(Ephem = 28, Plecop = 22, Trichop = 18, Chirono = 4, Oligoch = 1, Odonat = 2, Gastrop = 1),
"Run" = c(Ephem = 12, Plecop = 8, Trichop = 14, Chirono = 16, Oligoch = 6, Odonat = 5, Gastrop = 3),
"Pool" = c(Ephem = 2, Plecop = 0, Trichop = 4, Chirono = 32, Oligoch = 24, Odonat = 12, Gastrop = 18)
)
# Simulate Poisson counts
mat <- do.call(cbind, lapply(means, function(m) rpois(n_sites, lambda = m)))
df <- as.data.frame(mat)
df$Habitat <- hab
df
})
macro_data <- do.call(rbind, comm_list)
# Place metadata as first column
macro_data <- macro_data[, c("Habitat", setdiff(names(macro_data), "Habitat"))]
# Corresponding environmental matrix
env_data <- data.frame(
Velocity_ms = c(rnorm(n_sites, 0.85, 0.1), rnorm(n_sites, 0.45, 0.08), rnorm(n_sites, 0.12, 0.04)),
DO_mgl = c(rnorm(n_sites, 9.2, 0.4), rnorm(n_sites, 7.8, 0.5), rnorm(n_sites, 5.4, 0.6)),
Temp_C = c(rnorm(n_sites, 19.5, 0.8), rnorm(n_sites, 22.0, 0.7), rnorm(n_sites, 25.8, 1.1))
)
head(macro_data[, 1:6], 6)
#> Habitat Ephem Plecop Trichop Chirono Oligoch
#> 1 Riffle 35 15 23 3 2
#> 2 Riffle 26 21 19 3 2
#> 3 Riffle 39 22 19 2 0
#> 4 Riffle 19 18 24 4 1
#> 5 Riffle 25 31 19 0 0
#> 6 Riffle 26 29 25 1 00.10.3 Visualization Styles
nmds_analysis() displays sample site scores color-coded
by habitat with 95% confidence ellipses, and can overlay significant
environmental gradient vectors.
0.10.3.1 Style A: Community Ordination with 95% Confidence Ellipses
Computes Hellinger-transformed Bray-Curtis NMDS ordination and projects cluster confidence regions:
res_nmds <- nmds_analysis(
data = macro_data,
group_col = "Habitat",
transform = "hellinger",
dist_method = "bray",
palette = "Dark2"
)
# Render ordination plot
res_nmds$Plot
0.10.3.2 Style B:
Community Ordination with Environmental Vector Fitting
(envfit)
Supplying env_data automatically runs
vegan::envfit (999 permutations) and projects significant
(\(p \le 0.05\)) environmental
directional arrows into the ordination space:
res_nmds_env <- nmds_analysis(
data = macro_data,
group_col = "Habitat",
env_data = env_data,
transform = "hellinger",
dist_method = "bray",
palette = "Set2",
title = "Macroinvertebrate Community NMDS with Environmental Drivers"
)
res_nmds_env$Plot
0.10.4 Model Inspection and Numerical Summaries
The returned list provides the stress value, tidy site coordinate
scores, underlying metaMDS object, and permutation vector
fits:
# 1. Final Ordination Stress
res_nmds$Stress
#> [1] 0.03217201
# 2. Tidy NMDS Axis Coordinates
head(res_nmds$Plot_Data, 6)
#> NMDS1 NMDS2 Habitat
#> 1 -0.5252891 0.037631731 Riffle
#> 2 -0.4970030 -0.007769052 Riffle
#> 3 -0.5858803 -0.092789932 Riffle
#> 4 -0.4764525 -0.030577205 Riffle
#> 5 -0.6635751 -0.245218978 Riffle
#> 6 -0.7915610 -0.054026622 Riffle
# 3. Environmental Vector Fits (r2, p-values)
res_nmds_env$Envfit_Results
#>
#> ***VECTORS
#>
#> NMDS1 NMDS2 r2 Pr(>r)
#> Velocity_ms -0.92515 -0.37960 0.9032 0.001 ***
#> DO_mgl -0.99241 0.12297 0.9357 0.001 ***
#> Temp_C 0.97695 0.21348 0.8872 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Permutation: free
#> Number of permutations: 9990.11 Permutational Multivariate Analysis of Variance (PERMANOVA)
The permanova_analysis() function performs
non-parametric multivariate hypothesis testing on ecological distance
matrices using vegan::adonis2. While ordination techniques
such as NMDS graphically portray patterns of community resemblance,
PERMANOVA tests whether spatial, seasonal, or experimental groups occupy
statistically distinct regions in multivariate ecological space.
Because community data are non-normal, zero-inflated, and contain complex species co-dependencies, standard parametric MANOVA assumptions (multivariate normality and homogeneous covariance matrices) are regularly violated. PERMANOVA resolves this by partitioning dissimilarity matrices and calculating significance via permutation testing.
0.11.1 Conceptual Framework and Statistical Formulation
Let \(\mathbf{D} = [d_{ij}]\) represent an \(N \times N\) symmetric ecological distance matrix (e.g., Bray-Curtis or Euclidean distances) among \(N\) samples across \(g\) distinct groups. The total sum of squared dissimilarities (\(\text{SS}_{\text{Total}}\)) is partitioned into within-group (\(\text{SS}_W\)) and between-group (\(\text{SS}_A\)) components:
\[\text{SS}_{\text{Total}} = \frac{1}{N} \sum_{i=1}^{N-1} \sum_{j=i+1}^N d_{ij}^2\]
\[\text{SS}_W = \sum_{k=1}^g \frac{1}{n_k} \sum_{i=1}^{n_k-1} \sum_{j=i+1}^{n_k} d_{ij(k)}^2\]
\[\text{SS}_A = \text{SS}_{\text{Total}} - \text{SS}_W\]
where \(n_k\) is the sample size in group \(k\).
0.11.1.1 The Pseudo-\(F\) Test Statistic
The non-parametric pseudo-\(F\) statistic compares the between-group variance to within-group variance:
\[F = \frac{\text{SS}_A / (g - 1)}{\text{SS}_W / (N - g)}\]
Because the null distribution of pseudo-\(F\) does not follow an ordinary parametric \(F\)-distribution under arbitrary distance metrics, significance (\(p\)-value) is determined through Monte Carlo permutation of sample labels:
\[p = \frac{\sum_{m=1}^M \mathbb{I}(F^{\pi}_m \ge F_{\text{obs}}) + 1}{M + 1}\]
where \(M\) is the number of permutations (typically 999), \(F_{\text{obs}}\) is the pseudo-\(F\) calculated on the intact dataset, and \(F^{\pi}_m\) is the pseudo-\(F\) obtained from the \(m\)-th randomly permuted configuration.
0.11.2 In-Vignette Data Generation
We simulate a marine benthic community monitoring study assessing
macrobenthic invertebrate abundances across three seabed sediment
classes (Coarse Gravel, Fine Sand, and
Silty Mud):
set.seed(909)
n_reps <- 12
substrates <- c("Coarse Gravel", "Fine Sand", "Silty Mud")
benthic_survey <- do.call(rbind, lapply(substrates, function(sub) {
# Diagnostic species abundance baselines
means <- switch(sub,
"Coarse Gravel" = c(Polychaeta = 14, Bivalvia = 6, Amphipoda = 28, Decapoda = 12, Echinodermata = 16),
"Fine Sand" = c(Polychaeta = 22, Bivalvia = 18, Amphipoda = 15, Decapoda = 8, Echinodermata = 4),
"Silty Mud" = c(Polychaeta = 35, Bivalvia = 26, Amphipoda = 4, Decapoda = 2, Echinodermata = 1)
)
mat <- do.call(cbind, lapply(means, function(m) rpois(n_reps, lambda = m)))
df <- as.data.frame(mat)
df$Substrate <- sub
df
}))
# Organize Substrate as grouping metadata column
benthic_survey <- benthic_survey[, c("Substrate", setdiff(names(benthic_survey), "Substrate"))]
head(benthic_survey, 6)
#> Substrate Polychaeta Bivalvia Amphipoda Decapoda Echinodermata
#> 1 Coarse Gravel 17 10 23 11 15
#> 2 Coarse Gravel 14 15 33 9 21
#> 3 Coarse Gravel 14 12 40 15 23
#> 4 Coarse Gravel 11 6 31 13 14
#> 5 Coarse Gravel 19 5 26 13 15
#> 6 Coarse Gravel 9 7 23 10 140.11.3 Statistical Execution and Permutation Testing
We test for differences in community structure across substrate types using the Bray-Curtis distance metric with Hellinger pre-transformation:
0.11.4 Model Inspection and Numerical Summaries
The returned list provides the ANOVA-style summary table detailing degrees of freedom, sum-of-squares, \(R^2\) (variation explained by the grouping factor), pseudo-\(F\) values, and permutation-based \(p\)-values:
# 1. PERMANOVA Table (adonis2 results)
res_permanova$PERMANOVA_Table
#> Permutation test for adonis under reduced model
#> Permutation: free
#> Number of permutations: 999
#>
#> vegan::adonis2(formula = formula_obj, data = design_df, permutations = permutations, method = dist_method)
#> Df SumOfSqs R2 F Pr(>F)
#> Model 2 0.86869 0.89866 146.32 0.001 ***
#> Residual 33 0.09796 0.10134
#> Total 35 0.96665 1.00000
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2. Distance Metric Applied
res_permanova$Distance_Method
#> [1] "bray"
# 3. Cleaned and Transformed Data Preview
head(res_permanova$Cleaned_Data, 4)
#> Group Polychaeta Bivalvia Amphipoda Decapoda Echinodermata
#> 1 Coarse Gravel 0.4729527 0.3627381 0.5501196 0.3804430 0.4442617
#> 2 Coarse Gravel 0.3900947 0.4037864 0.5989121 0.3127716 0.4777665
#> 3 Coarse Gravel 0.3668997 0.3396831 0.6201737 0.3797773 0.4702700
#> 4 Coarse Gravel 0.3829708 0.2828427 0.6429101 0.4163332 0.43204940.12 Similarity Percentages (SIMPER) Analysis
The simper_analysis() function identifies which
individual taxa or species drive community-level differences between
ecological groups. While NMDS visualizes spatial or experimental
separation and PERMANOVA evaluates whether those differences are
statistically significant, SIMPER (Similarity Percentages) decomposes
the overall Bray-Curtis dissimilarity into percentage contributions
attributable to each individual taxon.
In fisheries ecology, benthic surveys, and aquatic biomonitoring, SIMPER serves as an indicator species detection tool—revealing which particular taxa characterize undisturbed reference sites versus human-impacted habitats, or which species drive seasonal community turnover.
0.12.1 Conceptual Framework and Mathematical Formulation
Consider two groups of samples, Group \(1\) and Group \(2\). The Bray-Curtis dissimilarity between sample \(j\) in Group \(1\) and sample \(k\) in Group \(2\) is given by:
\[d_{jk} = \frac{\sum_{i=1}^p |y_{ij} - y_{ik}|}{\sum_{i=1}^p (y_{ij} + y_{ik})}\]
SIMPER decomposes this total dissimilarity into the specific contribution of the \(i\)-th taxon:
\[\delta_{jk}^{(i)} = \frac{|y_{ij} - y_{ik}|}{\sum_{m=1}^p (y_{mj} + y_{mk})}\]
0.12.1.1 Average Contribution and Diagnostic Ratios
The average contribution of species \(i\) (\(\bar{\delta}_i\)) across all \(n_1 \times n_2\) pairwise comparisons between Group \(1\) and Group \(2\) is:
\[\bar{\delta}_i = \frac{1}{n_1 n_2} \sum_{j=1}^{n_1} \sum_{k=1}^{n_2} \delta_{jk}^{(i)}\]
Taxa are ranked in descending order of \(\bar{\delta}_i\). The cumulative percentage contribution of the top \(m\) species is:
\[\text{Cumulative Contribution (\%)} = \frac{\sum_{i=1}^m \bar{\delta}_i}{\sum_{i=1}^p \bar{\delta}_i} \times 100\]
A high average contribution combined with a consistent ratio of contribution to standard deviation indicates that a species is a reliable discriminator between the two environmental or biological groups.
0.12.2 In-Vignette Data Generation
We simulate a lake fish assemblage survey monitoring abundances
across three habitat zones (Littoral, Pelagic,
and Profundal). Six representative fish taxa are surveyed
across 36 sampling stations:
set.seed(1010)
n_stations <- 12
zones <- c("Littoral", "Pelagic", "Profundal")
lake_survey <- do.call(rbind, lapply(zones, function(zn) {
# Distinct abundance profiles across lake zones
means <- switch(zn,
"Littoral" = c(Perch = 32, Pike = 8, Roach = 45, Bream = 22, Bleak = 12, Zander = 2),
"Pelagic" = c(Perch = 10, Pike = 1, Roach = 18, Bream = 6, Bleak = 55, Zander = 14),
"Profundal" = c(Perch = 2, Pike = 0, Roach = 4, Bream = 1, Bleak = 2, Zander = 28)
)
mat <- do.call(cbind, lapply(means, function(m) rpois(n_stations, lambda = m)))
df <- as.data.frame(mat)
df$Zone <- zn
df
}))
# Organize Zone as grouping metadata column
lake_survey <- lake_survey[, c("Zone", setdiff(names(lake_survey), "Zone"))]
head(lake_survey, 6)
#> Zone Perch Pike Roach Bream Bleak Zander
#> 1 Littoral 32 10 45 34 8 0
#> 2 Littoral 26 8 49 17 12 3
#> 3 Littoral 37 6 38 25 11 3
#> 4 Littoral 43 7 47 24 15 0
#> 5 Littoral 36 6 30 20 17 1
#> 6 Littoral 37 11 42 21 11 10.12.3 Visualization Styles
The simper_analysis() function extracts the top \(N\) taxa contributing to dissimilarity for
each pairwise comparison and displays them using faceted horizontal bar
charts.
0.12.3.1 Style A: Faceted Top Contributing Taxa (Default)
Decomposes all pairwise group comparisons
(Littoral_Pelagic, Littoral_Profundal,
Pelagic_Profundal) and displays the top 5 discriminating
species per comparison:
res_simper <- simper_analysis(
data = lake_survey,
group_col = "Zone",
top_n = 5,
transform = "hellinger",
permutations = 999,
seed = 42
)
# Render faceted contribution plot
res_simper$Plot
0.12.4 Model Inspection and Numerical Summaries
The returned list contains the raw vegan::simper object,
along with a tidy data frame listing each species’ average contribution,
cumulative percentage, and permutation \(p\)-value:
# 1. Tidy summary table of top discriminating taxa
head(res_simper$Summary_Table, 10)
#> Comparison Species Contribution Cumulative p_value
#> 1 Littoral_Pelagic Bleak 0.08738423 0.2736940 0.001
#> 2 Littoral_Pelagic Zander 0.06289710 0.4706925 1.000
#> 3 Littoral_Pelagic Roach 0.04713252 0.6183151 0.010
#> 4 Littoral_Pelagic Pike 0.04134981 0.7478258 0.001
#> 5 Littoral_Pelagic Bream 0.04036191 0.8742423 0.976
#> 6 Littoral_Profundal Zander 0.19514771 0.3706790 0.001
#> 7 Littoral_Profundal Bream 0.08759999 0.5370734 0.001
#> 8 Littoral_Profundal Perch 0.08316557 0.6950446 0.001
#> 9 Littoral_Profundal Roach 0.07154487 0.8309426 0.001
#> 10 Littoral_Profundal Pike 0.06068598 0.9462144 0.001
# 2. Inspect raw vegan SIMPER summary for a specific comparison
summary(res_simper$SIMPER_Object)$Littoral_Pelagic[1:5, ]
#> average sd ratio ava avb cumsum p
#> Bleak 0.08738423 0.01342109 6.510964 0.3272274 0.70611825 0.2736940 0.001
#> Zander 0.06289710 0.01653698 3.803422 0.1018942 0.37505851 0.4706925 1.000
#> Roach 0.04713252 0.01538010 3.064513 0.6116993 0.40663320 0.6183151 0.010
#> Pike 0.04134981 0.01655756 2.497337 0.2377041 0.05898552 0.7478258 0.001
#> Bream 0.04036191 0.01918986 2.103294 0.4205971 0.24640228 0.8742423 0.9760.13 Principal Component Analysis (PCA) Community Biplot
The pca_community_biplot() function performs
unconstrained linear ordination using Principal Component Analysis
(vegan::rda()) to project aquatic communities into a
reduced orthogonal subspace. While non-metric multidimensional scaling
(NMDS) is non-parametric and rank-based, PCA constructs synthetic axes
that maximize total explained variance.
A key advantage of PCA for ecological communities is the simultaneous display of both sampling sites (points) and taxa (loadings) as directional vectors—forming a classical biplot. In this representation, the cosine of the angle between two species vectors approximates their linear correlation, while the projection of site points onto a species arrow indicates the relative abundance of that taxon at those sampling stations.
0.13.1 Conceptual Framework and Mathematical Formulation
Let \(\mathbf{Y}_{n \times p}\) denote an abundance or biomass matrix of \(n\) sampling sites and \(p\) aquatic species. PCA decomposes the centered (and optionally standardized) variance-covariance matrix \(\mathbf{S} = \frac{1}{n-1} \mathbf{Y}^T \mathbf{Y}\) via singular value decomposition (SVD) or spectral decomposition:
\[\mathbf{S} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^T\]
where: - \(\boldsymbol{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_p)\) is the diagonal matrix of eigenvalues sorted in descending order (\(\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0\)), - \(\mathbf{V} = [\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_p]\) is the matrix of orthonormal eigenvectors containing the species loadings.
The percentage of total community variance explained by the \(k\)-th principal component (\(\text{PC}_k\)) is:
\[\text{Variance Explained (\%)} = \frac{\lambda_k}{\sum_{m=1}^p \lambda_m} \times 100\]
Site coordinates in the reduced 2D ordination space are computed as:
\[\mathbf{F} = \mathbf{Y} \mathbf{V}_2\]
where \(\mathbf{V}_2\) denotes the first two eigenvector columns corresponding to \(\lambda_1\) and \(\lambda_2\).
0.13.1.1 Transformation-Based PCA (tb-PCA)
Raw Euclidean distance applied to unstandardized species counts can
suffer from the “horseshoe” or “arch” artifact because double zeros
(joint absences) are treated as evidence of ecological similarity. To
address this, transform = "hellinger" applies Legendre
& Gallagher’s Hellinger transformation prior to ordination:
\[y'_{ij} = \sqrt{\frac{y_{ij}}{\sum_{m=1}^p y_{mj}}}\]
This tb-PCA preserves Hellinger distances among samples, making PCA well suited for compositional abundance data.
0.13.2 In-Vignette Data Generation
We simulate an estuarine zooplankton community survey evaluating
counts of eight diagnostic taxa across three distinct salinity regimes
(Oligohaline, Mesohaline, and
Polyhaline):
set.seed(1111)
n_samples <- 14
salinity_regimes <- c("Oligohaline", "Mesohaline", "Polyhaline")
zooplankton_data <- do.call(rbind, lapply(salinity_regimes, function(reg) {
# Diagnostic zooplankton abundance profiles across salinity regimes
means <- switch(reg,
"Oligohaline" = c(Rotifer_A = 48, Rotifer_B = 35, Cladocera_A = 26, Cladocera_B = 18, Cyclopoid = 12, Calanoid = 4, Mysid = 1, Cirripede = 0),
"Mesohaline" = c(Rotifer_A = 12, Rotifer_B = 8, Cladocera_A = 14, Cladocera_B = 22, Cyclopoid = 32, Calanoid = 28, Mysid = 9, Cirripede = 5),
"Polyhaline" = c(Rotifer_A = 2, Rotifer_B = 0, Cladocera_A = 3, Cladocera_B = 4, Cyclopoid = 15, Calanoid = 42, Mysid = 24, Cirripede = 38)
)
mat <- do.call(cbind, lapply(means, function(m) rpois(n_samples, lambda = m)))
df <- as.data.frame(mat)
df$Salinity_Zone <- reg
df
}))
# Place metadata as first column
zooplankton_data <- zooplankton_data[, c("Salinity_Zone", setdiff(names(zooplankton_data), "Salinity_Zone"))]
head(zooplankton_data[, 1:6], 6)
#> Salinity_Zone Rotifer_A Rotifer_B Cladocera_A Cladocera_B Cyclopoid
#> 1 Oligohaline 47 34 26 15 9
#> 2 Oligohaline 57 29 28 12 12
#> 3 Oligohaline 52 48 19 15 16
#> 4 Oligohaline 56 30 26 15 10
#> 5 Oligohaline 48 42 29 16 13
#> 6 Oligohaline 27 36 27 23 140.13.3 Visualization Styles
pca_community_biplot() plots site coordinates, adds 95%
normal confidence ellipses around environmental groups, and displays
directional species loading vectors.
0.13.3.1 Style A: Hellinger-Transformed PCA Community Biplot (Default)
Applies Hellinger transformation to avoid double-zero distortions while visualizing all species loading vectors alongside salinity cluster groupings:
res_pca_biplot <- pca_community_biplot(
data = zooplankton_data,
group_col = "Salinity_Zone",
transform = "hellinger",
scale = FALSE,
color_palette = "Dark2",
arrow_multiplier = 0.8
)
# Render community biplot
res_pca_biplot$Plot
0.13.3.2 Style B: Biplot with Top-Contributing Species Filtered
In species-rich datasets, plotting all vectors can lead to visual
clutter. Setting top_n_species = 5 selectively displays
only the taxa with the longest loadings in the ordination plane:
res_pca_top5 <- pca_community_biplot(
data = zooplankton_data,
group_col = "Salinity_Zone",
transform = "hellinger",
top_n_species = 5,
color_palette = "Set1",
arrow_multiplier = 0.85,
title = "Zooplankton PCA Biplot (Top 5 Indicator Species)"
)
res_pca_top5$Plot
0.13.4 Model Inspection and Numerical Summaries
The returned list provides the underlying vegan::rda
model, coordinates of site scores with group labels, and species loading
vector magnitudes:
# 1. Inspect species vector coordinates and loading lengths
head(res_pca_biplot$Species_Scores, 6)
#> PC1 PC2 Species Vector_Length
#> Rotifer_A -0.55197481 0.078411347 Rotifer_A 0.5575164
#> Rotifer_B -0.59577347 0.007136423 Rotifer_B 0.5958162
#> Cladocera_A -0.34527594 -0.064900831 Cladocera_A 0.3513226
#> Cladocera_B -0.21734139 -0.190612221 Cladocera_B 0.2890853
#> Cyclopoid 0.05250865 -0.259792134 Cyclopoid 0.2650455
#> Calanoid 0.52674486 -0.185815423 Calanoid 0.5585584
# 2. Site ordination coordinates (PC1 and PC2 scores)
head(res_pca_biplot$Plot_Data, 6)
#> PC1 PC2 Salinity_Zone
#> sit1 -0.3314896 0.22274565 Oligohaline
#> sit2 -0.3336086 0.25922694 Oligohaline
#> sit3 -0.3051546 0.17602933 Oligohaline
#> sit4 -0.3545453 0.24684867 Oligohaline
#> sit5 -0.3411016 0.24157699 Oligohaline
#> sit6 -0.3295956 0.05113759 Oligohaline
# 3. Eigenvalues and total variance summary
summary(res_pca_biplot$PCA_Object)$cont$importance[, 1:2]
#> PC1 PC2
#> Eigenvalue 0.2023484 0.01946336
#> Proportion Explained 0.8698830 0.08367175
#> Cumulative Proportion 0.8698830 0.953554720.14 PCA Correlation Circle for Environmental Parameters
The pca_correlation_circle() function visualizes the
relationship among multidimensional environmental, limnological, and
hydrochemical variables using a PCA variables factor map (correlation
circle). When analyzing water quality matrices (e.g., pH, dissolved
oxygen, nutrients, heavy metal concentrations, conductivity), variables
are measured on different physical scales and often exhibit
collinearities.
The function standardizes environmental variables to unit variance (\(\mu = 0, \sigma = 1\)), extracts variable correlations with the leading principal components, evaluates the quality of representation using the squared cosine (\(\cos^2\)) metric, and projects the variables within a unit correlation circle.
0.14.1 Conceptual Framework and Mathematical Formulation
Let \(\mathbf{X}_{n \times p}\) denote an environmental dataset containing \(n\) water sampling stations and \(p\) physico-chemical parameters. Standardizing columns to unit variance yields:
\[z_{ij} = \frac{x_{ij} - \bar{x}_j}{s_j}\]
PCA on the standardized matrix decomposes the Pearson correlation matrix \(\mathbf{R} = \frac{1}{n-1} \mathbf{Z}^T \mathbf{Z}\):
\[\mathbf{R} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^T\]
The coordinate of variable \(j\) along principal component \(k\) represents the Pearson correlation coefficient between that environmental variable and component axis \(k\):
\[r(x_j, \text{PC}_k) = v_{jk} \sqrt{\lambda_k}\]
where \(v_{jk}\) is the eigenvector loading and \(\lambda_k\) is the eigenvalue of the \(k\)-th principal component.
0.14.1.1 The Unit Circle and \(\cos^2\) Metric
Because the variables are normalized, the maximum theoretical length of any vector in the full \(p\)-dimensional space is \(1.0\). When projected onto the 2D plane defined by \(\text{PC}_1\) and \(\text{PC}_2\), the squared distance from the origin represents the squared cosine (\(\cos^2\)):
\[\cos^2_j = r(x_j, \text{PC}_1)^2 + r(x_j, \text{PC}_2)^2\]
Interpretation criteria: - \(\cos^2 \to 1.0\) (Vector touches the unit circle): The variable is well represented in the selected 2D ordination plane; its position and angles accurately reflect relationships with other variables. - \(\cos^2 \to 0\) (Vector remains near the center): The variable is poorly represented in this 2D projection; inferences about its dynamics should be confirmed using higher-order principal components. - Angle Between Vectors: - \(\theta \approx 0^\circ\) (\(\cos \theta \to +1\)): Strong positive correlation. - \(\theta \approx 90^\circ\) (\(\cos \theta \to 0\)): Orthogonal (uncorrelated) variables. - \(\theta \approx 180^\circ\) (\(\cos \theta \to -1\)): Strong inverse correlation.
0.14.2 In-Vignette Data Generation
We simulate a limnological survey measuring eight physico-chemical water quality parameters across 40 sampling stations along a river-to-reservoir transition:
set.seed(1212)
n_stations <- 40
# Simulate correlated water quality gradients
temperature <- rnorm(n_stations, mean = 24.5, sd = 3.0)
do_mgl <- 14.5 - (0.35 * temperature) + rnorm(n_stations, 0, 0.5)
ph_val <- 7.2 + (0.04 * temperature) + rnorm(n_stations, 0, 0.3)
turbidity <- runif(n_stations, 5, 45)
secchi_cm <- 120 - (2.1 * turbidity) + rnorm(n_stations, 0, 4)
bod_mgl <- 1.2 + (0.08 * turbidity) + rnorm(n_stations, 0, 0.4)
nitrate_mgl <- rnorm(n_stations, mean = 2.4, sd = 0.6)
phosphate <- 0.15 + (0.05 * nitrate_mgl) + rnorm(n_stations, 0, 0.03)
water_quality <- data.frame(
StationID = paste0("ST_", sprintf("%02d", 1:n_stations)),
Temperature = temperature,
DO = do_mgl,
pH = ph_val,
Turbidity = turbidity,
Secchi = secchi_cm,
BOD = bod_mgl,
Nitrate = nitrate_mgl,
Phosphate = phosphate
)
head(water_quality, 6)
#> StationID Temperature DO pH Turbidity Secchi BOD Nitrate
#> 1 ST_01 22.61274 6.874055 7.905515 10.44149 104.99044 2.405120 2.033309
#> 2 ST_02 30.01227 3.943311 7.917997 16.77372 88.49811 2.137321 2.090744
#> 3 ST_03 25.53104 5.307665 8.238213 11.67667 97.99624 2.446044 2.602819
#> 4 ST_04 20.41971 6.597412 8.412772 18.24393 85.96159 2.811601 2.239365
#> 5 ST_05 26.10744 5.302458 7.879704 32.13295 45.15147 3.525980 2.781280
#> 6 ST_06 24.61160 5.835160 8.248662 19.75954 74.96844 3.245458 1.809519
#> Phosphate
#> 1 0.2599896
#> 2 0.2041083
#> 3 0.2793542
#> 4 0.2462373
#> 5 0.2738273
#> 6 0.24174230.14.3 Visualization Styles
The pca_correlation_circle() function generates a unit
circle with variables colored according to their representation quality
(\(\cos^2\)).
0.14.3.1 Style A: Standard Correlation Circle with \(\cos^2\) Gradient (PC1 vs PC2)
Excludes the identifier column
(exclude_cols = "StationID"), projects variables onto \(\text{PC}_1\) and \(\text{PC}_2\), and applies a color scale
based on representation quality:
res_pca_circle <- pca_correlation_circle(
data = water_quality,
exclude_cols = "StationID",
dim1 = 1,
dim2 = 2,
palette = c("#2A9D8F", "#E9C46A", "#F4A261", "#E76F51"),
title = "Water Quality Parameters Correlation Circle"
)
# Render correlation circle
res_pca_circle$Plot
0.14.4 Model Inspection and Numerical Summaries
The returned list provides the variable coordinate correlations, \(\cos^2\) representation values, and the variance breakdown for each principal component:
# 1. Variable coordinates along PC1 and PC2 with cos2 representation quality
res_pca_circle$Variable_Stats
#> x y Variable cos2
#> Temperature 0.6248845 0.6942111 Temperature 0.87240966
#> DO -0.5404295 -0.7708735 DO 0.88631008
#> pH 0.2395007 0.5537403 pH 0.36398892
#> Turbidity 0.9181992 -0.3522868 Turbidity 0.96719586
#> Secchi -0.9160075 0.3283661 Secchi 0.94689408
#> BOD 0.8790032 -0.3569626 BOD 0.90006888
#> Nitrate -0.1358684 0.2706143 Nitrate 0.09169232
#> Phosphate -0.0388324 0.2119688 Phosphate 0.04643871
# 2. Percentage and cumulative variance explained per principal component
res_pca_circle$Variance_Summary
#> Dimension Eigenvalue Variance_Percent Cumulative_Variance
#> 1 PC1 3.21467972 40.1834965 40.18350
#> 2 PC2 1.86031879 23.2539849 63.43748
#> 3 PC3 1.61039716 20.1299644 83.56745
#> 4 PC4 0.79743347 9.9679184 93.53536
#> 5 PC5 0.35996347 4.4995434 98.03491
#> 6 PC6 0.10354818 1.2943522 99.32926
#> 7 PC7 0.04227383 0.5284228 99.85768
#> 8 PC8 0.01138539 0.1423173 100.00000
# 3. Individual station ordination scores (first 6 samples)
head(res_pca_circle$Sample_Scores[, 1:4], 6)
#> PC1 PC2 PC3 PC4
#> 1 -1.7931263 -0.34086385 0.70640904 0.5390157
#> 2 0.1295898 1.55184446 2.16313645 1.4574737
#> 3 -0.8845816 1.32301245 -0.06146701 0.2188447
#> 4 -0.8460891 -0.38276489 0.51581956 -1.1101072
#> 5 1.5072193 -0.22394477 -0.70224218 0.9365320
#> 6 0.2438725 -0.05787865 1.22990777 -0.15597140.15 Redundancy Analysis (RDA)
The rda_analysis() function performs direct gradient
analysis (constrained linear ordination) to model the variation in
aquatic community composition directly explained by a set of continuous
environmental variables. While unconstrained methods (PCA, NMDS)
ordinate community matrices independently of external drivers,
Redundancy Analysis combines multivariate multiple linear regression
with principal component analysis.
The resulting ordination axes represent linear combinations of explanatory environmental predictors, ordered by the amount of community variance they explain. The function coordinates data cleaning across both matrices, evaluates global and marginal variable significance through non-parametric Monte Carlo permutations, and renders a publication-ready ordination Triplot (simultaneously presenting sampling sites, species loadings, and environmental vectors).
0.15.1 Conceptual Framework and Mathematical Formulation
Let \(\mathbf{Y}_{n \times p}\) denote the community abundance matrix (Hellinger-transformed) across \(n\) sampling stations and \(p\) aquatic species, and let \(\mathbf{X}_{n \times m}\) denote the standardized matrix of \(m\) environmental predictors.
RDA executes in two mathematical phases:
0.15.1.1 Phase 1: Multivariate Multiple Linear Regression
Each species column \(\mathbf{y}_j\) is regressed onto the complete environmental matrix \(\mathbf{X}\):
\[\hat{\mathbf{Y}} = \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y} = \mathbf{H} \mathbf{Y}\]
where \(\hat{\mathbf{Y}}\) contains the fitted values explained by the environmental variables, and \(\mathbf{H}\) represents the projection “hat” matrix. The residual matrix \(\mathbf{Y}_{\text{res}} = \mathbf{Y} - \hat{\mathbf{Y}}\) captures unconstrained community variation.
0.15.1.2 Phase 2: Principal Component Analysis on Fitted Values
PCA is conducted on the matrix of fitted values \(\hat{\mathbf{Y}}\) via singular value decomposition:
\[\frac{1}{n-1} \hat{\mathbf{Y}}^T \hat{\mathbf{Y}} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^T\]
The canonical eigenvalues \(\lambda_k \in \boldsymbol{\Lambda}\) measure the variance explained by each constrained canonical axis (\(\text{RDA}_k\)). The constrained proportion of total community inertia is given by:
\[R^2_{\text{biplot}} = \frac{\sum_{k=1}^m \lambda_k}{\text{Total Community Inertia}}\]
0.15.1.3 Permutation Hypothesis Testing
Significance is determined via permutation testing
(permutations = 999) using pseudo-\(F\) statistics:
\[F = \frac{\text{SS}_{\text{constrained}} / m}{\text{SS}_{\text{residual}} / (n - m - 1)}\]
Marginal permutation testing (by = "margin") tests the
unique partial contribution of each environmental variable while
controlling for all other predictors in the model.
0.15.2 In-Vignette Data Generation
We simulate a riverine fish community monitoring study assessing the
abundances of seven freshwater fish species across 36 sampling stations
spanning three distinct River Zones (Headwater,
Mid-River, Lowland Floodplain). The community
matrix is matched with four continuous water quality parameters (Water
Temperature, Dissolved Oxygen, Flow Velocity, and Conductivity):
set.seed(1313)
n_sites <- 12
zones <- c("Headwater", "Mid-River", "Lowland")
# Simulated environmental matrix
env_list <- lapply(zones, function(zn) {
switch(zn,
"Headwater" = data.frame(
Temp_C = rnorm(n_sites, 14.5, 0.8),
DO_mgL = rnorm(n_sites, 10.2, 0.4),
Flow_ms = rnorm(n_sites, 0.95, 0.1),
Cond_uS = rnorm(n_sites, 120, 15)
),
"Mid-River" = data.frame(
Temp_C = rnorm(n_sites, 19.0, 0.7),
DO_mgL = rnorm(n_sites, 8.1, 0.5),
Flow_ms = rnorm(n_sites, 0.55, 0.08),
Cond_uS = rnorm(n_sites, 280, 25)
),
"Lowland" = data.frame(
Temp_C = rnorm(n_sites, 24.5, 1.0),
DO_mgL = rnorm(n_sites, 5.8, 0.6),
Flow_ms = rnorm(n_sites, 0.18, 0.05),
Cond_uS = rnorm(n_sites, 490, 35)
)
)
})
env_survey <- do.call(rbind, env_list)
# Matched community counts (7 fish species)
comm_list <- lapply(zones, function(zn) {
means <- switch(zn,
"Headwater" = c(Trout = 30, Sculpin = 22, Dace = 12, Minnow = 4, Barb = 1, Catfish = 0, Carp = 0),
"Mid-River" = c(Trout = 6, Sculpin = 8, Dace = 24, Minnow = 28, Barb = 14, Catfish = 4, Carp = 2),
"Lowland" = c(Trout = 0, Sculpin = 1, Dace = 5, Minnow = 12, Barb = 32, Catfish = 26, Carp = 38)
)
mat <- do.call(cbind, lapply(means, function(m) rpois(n_sites, lambda = m)))
df <- as.data.frame(mat)
df$Zone <- zn
df
})
comm_survey <- do.call(rbind, comm_list)
# Metadata as first column
comm_survey <- comm_survey[, c("Zone", setdiff(names(comm_survey), "Zone"))]
head(comm_survey[, 1:6], 4)
#> Zone Trout Sculpin Dace Minnow Barb
#> 1 Headwater 25 22 11 7 0
#> 2 Headwater 26 23 19 6 0
#> 3 Headwater 32 15 12 3 2
#> 4 Headwater 29 17 7 6 1
head(env_survey, 4)
#> Temp_C DO_mgL Flow_ms Cond_uS
#> 1 15.45148 10.011911 0.9047721 125.05351
#> 2 14.94111 9.966365 0.9534446 127.99818
#> 3 14.78291 10.219464 0.9291478 135.87260
#> 4 13.93682 10.100078 1.1271965 95.687840.15.3 Visualization Styles
The rda_analysis() function creates a triplot
integrating site distributions, taxon vectors, and environmental
directional arrows.
0.15.3.1 Style A: Constrained Hellinger RDA Triplot (Default)
Applies Hellinger pre-transformation to the community matrix, standardizes environmental variables, projects 95% group confidence ellipses, and displays repelled labels for species and environmental drivers:
res_rda <- rda_analysis(
comm_data = comm_survey,
env_data = env_survey,
group_col = "Zone",
transform = "hellinger",
color_palette = "Dark2",
species_arrow_mult = 1.0,
env_arrow_mult = 1.2
)
# Render publication triplot
res_rda$Plot
0.15.4 Model Inspection and Numerical Summaries
The returned list provides the underlying vegan::rda
object, global model permutation significance, and marginal variable
significance:
# 1. Global RDA Model Permutation Test (adonis/anova)
res_rda$Model_Significance
#> Permutation test for rda under reduced model
#> Permutation: free
#> Number of permutations: 999
#>
#> Model: rda(formula = comm_mat ~ Temp_C + DO_mgL + Flow_ms + Cond_uS, data = env_mat_scaled)
#> Df Variance F Pr(>F)
#> Model 4 0.268865 40.447 0.001 ***
#> Residual 31 0.051517
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2. Marginal Environmental Term Significance (unique partial effects)
res_rda$Variable_Significance
#> Permutation test for rda under reduced model
#> Marginal effects of terms
#> Permutation: free
#> Number of permutations: 999
#>
#> Model: rda(formula = comm_mat ~ Temp_C + DO_mgL + Flow_ms + Cond_uS, data = env_mat_scaled)
#> Df Variance F Pr(>F)
#> Temp_C 1 0.010027 6.0334 0.014 *
#> DO_mgL 1 0.005364 3.2277 0.066 .
#> Flow_ms 1 0.002476 1.4900 0.232
#> Cond_uS 1 0.002079 1.2511 0.271
#> Residual 31 0.051517
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3. Constrained vs Unconstrained Inertia Breakdown
summary(res_rda$RDA_Object)$tot.chi
#> [1] 0.3203818
summary(res_rda$RDA_Object)$concont$importance[, 1:2]
#> RDA1 RDA2
#> Eigenvalue 0.2582838 0.009324137
#> Proportion Explained 0.9606456 0.034679644
#> Cumulative Proportion 0.9606456 0.9953252810.16 Distance-based Redundancy Analysis (db-RDA)
The dbrda_analysis() function performs constrained
ordination using non-Euclidean ecological distance metrics via
Distance-based Redundancy Analysis (vegan::dbrda).
Classical Redundancy Analysis (RDA) relies on Euclidean distances, which
can lead to distortion when analyzing zero-skewed ecological abundance
matrices.
Distance-based RDA resolves this limitation by conducting Principal Coordinate Analysis (PCoA) on an arbitrary dissimilarity matrix (e.g., Bray-Curtis, Jaccard, or Kulczynski) and using the resulting metric coordinate axes as the response matrix in a standard RDA. This enables direct multivariate regression of semi-metric ecological dissimilarities against continuous environmental predictors while avoiding negative eigenvalue artifacts.
0.16.1 Conceptual Framework and Mathematical Formulation
Let \(\mathbf{D} = [d_{ij}]_{n \times n}\) denote an ecological distance matrix among \(n\) sampling sites. The matrix \(\mathbf{D}\) is centered via Gower’s centering transformation:
\[\mathbf{G} = -\frac{1}{2} \left( \mathbf{I} - \frac{1}{n} \mathbf{11}^T \right) (\mathbf{D}^{(2)}) \left( \mathbf{I} - \frac{1}{n} \mathbf{11}^T \right)\]
where \(\mathbf{D}^{(2)} = [d_{ij}^2]\) is the matrix of squared dissimilarities.
Spectral decomposition of \(\mathbf{G}\) produces eigenvalues \(\boldsymbol{\Lambda}\) and eigenvectors \(\mathbf{V}\):
\[\mathbf{G} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^T\]
Principal coordinate vectors with positive eigenvalues \(\mathbf{Z} = \mathbf{V}_+ \boldsymbol{\Lambda}_+^{1/2}\) represent the Euclidean embedding of the original dissimilarity matrix.
0.16.1.1 Constrained Regression and Model Permutation
The coordinate matrix \(\mathbf{Z}\) is regressed onto the standardized environmental predictor matrix \(\mathbf{X}_{n \times m}\):
\[\hat{\mathbf{Z}} = \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Z}\]
PCA is subsequently applied to \(\hat{\mathbf{Z}}\) to construct the canonical axes (\(\text{dbRDA}_1, \text{dbRDA}_2\)). Overall model significance and the marginal contributions of individual environmental variables are evaluated through permutation testing (pseudo-\(F\) tests with 999 Monte Carlo permutations):
\[F = \frac{\text{tr}(\hat{\mathbf{Z}}^T \hat{\mathbf{Z}}) / m}{\text{tr}(\mathbf{E}^T \mathbf{E}) / (n - m - 1)}\]
where \(\mathbf{E} = \mathbf{Z} - \hat{\mathbf{Z}}\) represents the residual matrix of unconstrained variation.
0.16.2 In-Vignette Data Generation
We simulate an estuarine benthic macroinvertebrate community across
36 stations classified into three Hydrological Zones
(Upper Estuary, Middle Basin,
Outer Marine). The community dataset is paired with four
continuous sediment and water column metrics (Salinity, Total Organic
Carbon, Silt Fraction, and Water Depth):
set.seed(1414)
n_stations <- 12
zones <- c("Upper Estuary", "Middle Basin", "Outer Marine")
# Simulated environmental predictors
env_list <- lapply(zones, function(zn) {
switch(zn,
"Upper Estuary" = data.frame(
Salinity_PSU = rnorm(n_stations, 8.5, 1.2),
TOC_percent = rnorm(n_stations, 4.2, 0.5),
Silt_percent = rnorm(n_stations, 68, 6),
Depth_m = rnorm(n_stations, 3.5, 0.4)
),
"Middle Basin" = data.frame(
Salinity_PSU = rnorm(n_stations, 18.0, 1.5),
TOC_percent = rnorm(n_stations, 2.6, 0.4),
Silt_percent = rnorm(n_stations, 42, 5),
Depth_m = rnorm(n_stations, 8.2, 0.8)
),
"Outer Marine" = data.frame(
Salinity_PSU = rnorm(n_stations, 32.5, 0.9),
TOC_percent = rnorm(n_stations, 1.1, 0.2),
Silt_percent = rnorm(n_stations, 15, 3),
Depth_m = rnorm(n_stations, 18.5, 1.5)
)
)
})
env_estuary <- do.call(rbind, env_list)
# Simulated benthic community matrix (6 diagnostic taxa)
comm_list <- lapply(zones, function(zn) {
means <- switch(zn,
"Upper Estuary" = c(Oligochaeta = 42, Chironomidae = 35, Polychaeta_A = 6, Bivalvia_A = 2, Amphipoda_A = 4, Decapoda_A = 0),
"Middle Basin" = c(Oligochaeta = 14, Chironomidae = 18, Polychaeta_A = 28, Bivalvia_A = 22, Amphipoda_A = 16, Decapoda_A = 5),
"Outer Marine" = c(Oligochaeta = 1, Chironomidae = 0, Polychaeta_A = 48, Bivalvia_A = 34, Amphipoda_A = 30, Decapoda_A = 18)
)
mat <- do.call(cbind, lapply(means, function(m) rpois(n_stations, lambda = m)))
df <- as.data.frame(mat)
df$Estuary_Zone <- zn
df
})
comm_estuary <- do.call(rbind, comm_list)
# Set grouping metadata as first column
comm_estuary <- comm_estuary[, c("Estuary_Zone", setdiff(names(comm_estuary), "Estuary_Zone"))]
head(comm_estuary[, 1:6], 4)
#> Estuary_Zone Oligochaeta Chironomidae Polychaeta_A Bivalvia_A Amphipoda_A
#> 1 Upper Estuary 45 32 11 0 3
#> 2 Upper Estuary 35 26 6 3 11
#> 3 Upper Estuary 38 29 4 1 1
#> 4 Upper Estuary 47 37 3 2 3
head(env_estuary, 4)
#> Salinity_PSU TOC_percent Silt_percent Depth_m
#> 1 7.599730 3.919065 56.09989 3.218138
#> 2 8.395406 3.957087 79.05169 3.487593
#> 3 9.941565 3.902267 64.37852 3.872967
#> 4 8.293515 4.436544 50.99208 3.2626030.16.3 Visualization Styles
The dbrda_analysis() function displays the sample
coordinates, 95% group confidence ellipses, species vectors, and
environmental directional arrows on a single coordinate grid.
0.16.3.1 Style A: Bray-Curtis db-RDA Triplot (Default)
Computes Bray-Curtis dissimilarity, applies environmental scaling, and produces the publication triplot:
res_dbrda <- dbrda_analysis(
comm_data = comm_estuary,
env_data = env_estuary,
group_col = "Estuary_Zone",
dist_method = "bray",
transform = "none",
color_palette = "Dark2",
species_arrow_mult = 1.0,
env_arrow_mult = 1.1
)
# Render publication triplot
res_dbrda$Plot
0.16.4 Model Inspection and Numerical Summaries
The returned list provides the underlying vegan::dbrda
model object, the global permutation significance test, and marginal
term significance:
# 1. Global db-RDA Permutation Test (adonis/anova)
res_dbrda$Model_Significance
#> Permutation test for dbrda under reduced model
#> Permutation: free
#> Number of permutations: 999
#>
#> Model: vegan::dbrda(formula = comm_mat ~ Salinity_PSU + TOC_percent + Silt_percent + Depth_m, data = env_mat_scaled, distance = dist_method)
#> Df SumOfSqs F Pr(>F)
#> Model 4 4.2675 60.5 0.001 ***
#> Residual 31 0.5467
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2. Marginal Environmental Predictor Significance (unique partial variance)
res_dbrda$Variable_Significance
#> Permutation test for dbrda under reduced model
#> Marginal effects of terms
#> Permutation: free
#> Number of permutations: 999
#>
#> Model: vegan::dbrda(formula = comm_mat ~ Salinity_PSU + TOC_percent + Silt_percent + Depth_m, data = env_mat_scaled, distance = dist_method)
#> Df SumOfSqs F Pr(>F)
#> Salinity_PSU 1 0.04813 2.7295 0.108
#> TOC_percent 1 0.12682 7.1919 0.012 *
#> Silt_percent 1 -0.00053 -0.0303 0.998
#> Depth_m 1 0.02884 1.6355 0.225
#> Residual 31 0.54666
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3. Constrained vs Unconstrained Variance Breakdown
summary(res_dbrda$dbRDA_Object)$tot.chi
#> [1] 4.814109
summary(res_dbrda$dbRDA_Object)$concont$importance[, 1:2]
#> dbRDA1 dbRDA2
#> Eigenvalue 4.2176612 0.05124736
#> Proportion Explained 0.9883327 0.01200889
#> Cumulative Proportion 0.9883327 1.000341610.17 Hierarchical Clustering Dendrogram
The cluster_dendrogram() function performs agglomerative
hierarchical cluster analysis (HCA) on ecological community abundance
matrices or multi-parameter environmental datasets. While ordination
approaches (e.g., NMDS, PCA, RDA) display continuous gradients across
reduced orthogonal dimensions, hierarchical clustering groups sampling
units into discrete nested clusters based on pairwise resemblance.
The function computes dissimilarity matrices using standard
ecological distance metrics (e.g., Bray-Curtis, Jaccard, Euclidean),
applies optional community transformations, conducts agglomerative
clustering across user-selected linkage criteria (UPGMA, Ward’s minimum
variance, Complete), and generates a native ggplot2 tree
diagram complete with color-coded group leaves, an optional cut-off
threshold line, and sample cluster assignments.
0.17.1 Conceptual Framework and Mathematical Formulation
Let \(\mathbf{Y}_{n \times p}\) denote a community matrix of \(n\) sampling stations and \(p\) aquatic taxa. Pairwise ecological dissimilarity \(d_{ij}\) is computed using the Bray-Curtis metric:
\[d_{ij} = \frac{\sum_{k=1}^p |y_{ik} - y_{jk}|}{\sum_{k=1}^p (y_{ik} + y_{jk})}\]
Agglomerative hierarchical clustering begins with each sampling unit treated as a single-element cluster (\(C_i = \{i\}\)). At each sequential step, the two clusters with the minimum inter-cluster distance are merged until all samples are unified into a single root node.
0.17.1.1 Agglomerative Linkage Criteria
The definition of distance between two composite clusters \(A\) and \(B\) depends on the chosen linkage method:
Unweighted Pair-Group Method with Arithmetic Mean (UPGMA /
cluster_method = "average"): Calculates the mean pairwise dissimilarity between all members of cluster \(A\) and cluster \(B\): \[D(A, B) = \frac{1}{|A| \cdot |B|} \sum_{i \in A} \sum_{j \in B} d_{ij}\] Preserves cophenetic correlation with minimal spatial distortion.Ward’s Minimum Variance (
cluster_method = "ward.D2"): Merges the cluster pair that minimizes the increase in total within-cluster error sum of squares (ESS): \[\Delta \text{ESS}_{AB} = \frac{|A| \cdot |B|}{|A| + |B|} \|\bar{\mathbf{y}}_A - \bar{\mathbf{y}}_B\|^2\] Tends to produce compact, spherical, and balanced cluster topologies.Complete Linkage (
cluster_method = "complete"): Defines distance as the maximum pairwise dissimilarity between clusters: \[D(A, B) = \max \{d_{ij} : i \in A, \, j \in B\}\]
0.17.2 In-Vignette Data Generation
We simulate a limnological biomonitoring survey assessing
phytoplankton community counts across 24 sampling stations situated in
three lake zones (Littoral, Pelagic,
Hypolimnetic Inlet):
set.seed(1515)
n_per_zone <- 8
lake_zones <- c("Littoral", "Pelagic", "Hypolimnetic")
phyto_list <- lapply(lake_zones, function(zn) {
means <- switch(zn,
"Littoral" = c(Diatom_A = 35, Diatom_B = 24, Cyanob_A = 4, Cyanob_B = 1, Green_A = 18, Green_B = 12),
"Pelagic" = c(Diatom_A = 12, Diatom_B = 8, Cyanob_A = 42, Cyanob_B = 30, Green_A = 25, Green_B = 18),
"Hypolimnetic" = c(Diatom_A = 4, Diatom_B = 2, Cyanob_A = 2, Cyanob_B = 0, Green_A = 6, Green_B = 4)
)
mat <- do.call(cbind, lapply(means, function(m) rpois(n_per_zone, lambda = m)))
df <- as.data.frame(mat)
df$Zone <- zn
df
})
phyto_survey <- do.call(rbind, phyto_list)
rownames(phyto_survey) <- paste0("ST_", sprintf("%02d", seq_len(nrow(phyto_survey))))
# Organize grouping metadata column
phyto_survey <- phyto_survey[, c("Zone", setdiff(names(phyto_survey), "Zone"))]
head(phyto_survey[, 1:5], 5)
#> Zone Diatom_A Diatom_B Cyanob_A Cyanob_B
#> ST_01 Littoral 29 24 7 3
#> ST_02 Littoral 33 16 3 0
#> ST_03 Littoral 40 14 1 1
#> ST_04 Littoral 40 30 1 3
#> ST_05 Littoral 34 29 4 50.17.3 Visualization Styles
The cluster_dendrogram() function generates vertical and
horizontal dendrogram views with group-coded leaf tips.
0.17.3.1 Style A: Vertical Tree Layout with Cluster Partitioning Line (Default)
Applies Hellinger pre-transformation with Bray-Curtis dissimilarity
and UPGMA (average) clustering, displaying a threshold line
at \(k = 3\):
res_dendro <- cluster_dendrogram(
data = phyto_survey,
group_col = "Zone",
dist_method = "bray",
transform = "hellinger",
cluster_method = "average",
k = 3,
palette = "Dark2",
orientation = "vertical"
)
# Render publication tree
res_dendro$Plot
0.17.3.2 Style B: Horizontal Tree Layout
Using orientation = "horizontal" flips the coordinate
axes, making station labels easier to read when sample sizes are
larger:
res_dendro_horiz <- cluster_dendrogram(
data = phyto_survey,
group_col = "Zone",
dist_method = "bray",
transform = "hellinger",
cluster_method = "ward.D2",
k = 3,
palette = "Dark2",
orientation = "horizontal",
title = "Phytoplankton Assemblage Clusters (Ward.D2 Linkage)"
)
res_dendro_horiz$Plot
0.17.4 Model Inspection and Numerical Summaries
The returned list provides the sample-to-cluster assignment table,
dissimilarity matrix, and underlying hclust object:
# 1. Sample Cluster Membership Assignments
head(res_dendro$Cluster_Assignments, 10)
#> Sample Zone Cluster
#> 1 ST_01 Littoral Cluster_1
#> 2 ST_02 Littoral Cluster_1
#> 3 ST_03 Littoral Cluster_1
#> 4 ST_04 Littoral Cluster_1
#> 5 ST_05 Littoral Cluster_1
#> 6 ST_06 Littoral Cluster_1
#> 7 ST_07 Littoral Cluster_1
#> 8 ST_08 Littoral Cluster_1
#> 9 ST_09 Pelagic Cluster_2
#> 10 ST_10 Pelagic Cluster_2
# 2. Inspect Agglomerative Cluster Merges
res_dendro$HClust_Object$merge[1:6, ]
#> [,1] [,2]
#> [1,] -13 -16
#> [2,] -14 1
#> [3,] -4 -7
#> [4,] -5 -6
#> [5,] -15 2
#> [6,] -9 5
# 3. Distance Matrix Preview (First 4 Samples)
as.matrix(res_dendro$Distance_Matrix)[1:4, 1:4]
#> ST_01 ST_02 ST_03 ST_04
#> ST_01 0.00000000 0.1173526 0.11920702 0.06805391
#> ST_02 0.11735263 0.0000000 0.10490610 0.10833905
#> ST_03 0.11920702 0.1049061 0.00000000 0.08964695
#> ST_04 0.06805391 0.1083391 0.08964695 0.000000000.18 Correlation Heatmaps for Environmental & Morphometric Matrices
The correlation_heatmap() function generates
publication-ready correlation heatmaps for multidimensional
environmental, limnological, and fisheries datasets. In aquatic
ecosystem monitoring, evaluating pairwise associations among water
quality parameters (e.g., Temperature, Dissolved Oxygen, Nutrients,
Turbidity) is an essential exploratory step before conducting
constrained ordinations or regression modeling.
The function computes correlation coefficients across Pearson,
Spearman, or Kendall formulations, extracts pairwise significance \(p\)-values, supports directional matrix
views (lower, upper, or full),
provides alternative tile geometries (square or
circle), and dynamically adjusts annotation font contrast
across diverging color palettes.
0.18.1 Conceptual Framework and Statistical Formulation
For two continuous variables \(X\) and \(Y\) observed across \(n\) sampling stations:
0.18.1.1 Pearson Product-Moment Correlation
Measures linear dependency between normally distributed continuous variables:
\[r_{xy} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^n (x_i - \bar{x})^2 \sum_{i=1}^n (y_i - \bar{y})^2}}\]
The statistical significance of \(r_{xy}\) is evaluated using a two-tailed Student’s \(t\)-distribution with \(\nu = n - 2\) degrees of freedom:
\[t = r_{xy} \sqrt{\frac{n - 2}{1 - r_{xy}^2}}\]
0.18.1.2 Spearman Rank-Order Correlation
Measures monotonic relationships between non-normally distributed or ordinal parameters by computing Pearson’s coefficient on ranked values:
\[\rho = 1 - \frac{6 \sum_{i=1}^n d_i^2}{n(n^2 - 1)}\]
where \(d_i = \text{rank}(x_i) - \text{rank}(y_i)\).
Pairwise significance thresholds are mapped to standard compact asterisks: - \(\text{ns}: p \ge 0.05\) - \(*: p < 0.05\) - \(**: p < 0.01\) - \(***: p < 0.001\)
0.18.2 In-Vignette Data Generation
We simulate a multi-parameter aquatic environmental survey evaluating seven correlated physico-chemical variables recorded across 50 limnological stations:
set.seed(1616)
n_stations <- 50
temp_c <- rnorm(n_stations, mean = 23.5, sd = 2.8)
do_mgl <- 13.8 - (0.32 * temp_c) + rnorm(n_stations, 0, 0.4)
ph_val <- 7.1 + (0.05 * temp_c) + rnorm(n_stations, 0, 0.25)
turb_ntu <- runif(n_stations, 4, 38)
secchi_m <- 2.8 - (0.055 * turb_ntu) + rnorm(n_stations, 0, 0.2)
tp_mgl <- 0.04 + (0.006 * turb_ntu) + rnorm(n_stations, 0, 0.02)
tn_mgl <- 0.85 + (12.5 * tp_mgl) + rnorm(n_stations, 0, 0.15)
limno_env <- data.frame(
Temperature = temp_c,
DO = do_mgl,
pH = ph_val,
Turbidity = turb_ntu,
Secchi = secchi_m,
Total_P = tp_mgl,
Total_N = tn_mgl
)
head(limno_env, 6)
#> Temperature DO pH Turbidity Secchi Total_P Total_N
#> 1 22.19114 5.982808 8.407656 8.327564 2.3928749 0.09957881 1.912191
#> 2 26.59665 6.100091 8.240560 24.557620 1.2303458 0.15707745 2.792743
#> 3 25.11247 6.161079 8.381893 34.516946 0.6536182 0.24425880 3.752906
#> 4 22.64647 6.449037 8.274918 9.993436 2.0100879 0.11489375 2.498569
#> 5 23.06031 7.278349 8.215123 30.431676 1.0822492 0.19719694 3.357934
#> 6 30.20779 2.871502 8.630534 9.455666 2.2882730 0.09409752 2.1901120.18.3 Visualization Styles
The correlation_heatmap() function provides tile-based
and bubble-based visual styles.
0.18.3.1 Style A: Lower-Triangle Square Heatmap with Coefficients & Asterisks (Default)
Displays the lower triangle with correlation values and significance stars, automatically adjusting text color for high-contrast visibility:
res_cor_square <- correlation_heatmap(
data = limno_env,
method = "pearson",
shape = "square",
view = "lower",
diag = FALSE,
label_type = "both",
color_palette = "RdBu"
)
# Render lower-triangle heatmap
res_cor_square$Plot
0.18.3.2 Style B: Circle Bubble Heatmap with Significance Asterisks
Uses circle markers where bubble diameter corresponds to the absolute correlation magnitude (\(\vert{}r\vert{}\)) and color corresponds to directional correlation sign:
res_cor_circle <- correlation_heatmap(
data = limno_env,
method = "spearman",
shape = "circle",
view = "lower",
diag = FALSE,
label_type = "stars",
color_palette = "RdYlBu",
title = "Environmental Matrix Spearman Correlations"
)
res_cor_circle$Plot
0.18.4 Model Inspection and Numerical Summaries
The returned list provides the symmetric correlation coefficient matrix, exact \(p\)-values, and a tidy summary data frame:
# 1. Pearson Correlation Matrix Preview
round(res_cor_square$Correlation_Matrix[1:4, 1:4], 3)
#> Temperature DO pH Turbidity
#> Temperature 1.000 -0.915 0.652 0.218
#> DO -0.915 1.000 -0.618 -0.110
#> pH 0.652 -0.618 1.000 0.036
#> Turbidity 0.218 -0.110 0.036 1.000
# 2. Exact P-Value Matrix Preview
signif(res_cor_square$P_Value_Matrix[1:4, 1:4], 3)
#> Temperature DO pH Turbidity
#> Temperature 0.00e+00 1.67e-20 2.99e-07 0.129
#> DO 1.67e-20 0.00e+00 1.71e-06 0.446
#> pH 2.99e-07 1.71e-06 0.00e+00 0.802
#> Turbidity 1.29e-01 4.46e-01 8.02e-01 0.000
# 3. Tidy Pairwise Summary Table
head(res_cor_square$Summary_Table, 8)
#> Var1 Var2 Correlation P_Value Significance
#> 3 DO Temperature -0.91457645 1.674026e-20 ***
#> 8 Secchi DO 0.16852240 2.420479e-01 ns
#> 10 Secchi Temperature -0.27450426 5.370824e-02 ns
#> 13 Secchi Turbidity -0.93266093 6.826376e-23 ***
#> 14 Secchi pH -0.05415836 7.087424e-01 ns
#> 22 Total_N DO -0.08583431 5.533923e-01 ns
#> 23 Total_N Secchi -0.86826669 3.197710e-16 ***
#> 24 Total_N Temperature 0.13998938 3.322344e-01 ns0.19 Integrated Environmental Correlation and Mantel Network Analysis
The mantel_heatmap_analysis() function models and
visualizes the coupling between abiotic environmental gradients and
biological community structures. In aquatic ecology, researchers often
analyze pairwise physico-chemical parameters alongside multivariate
dissimilarity matrices across multiple trophic levels (e.g.,
Phytoplankton, Zooplankton, Benthos, Fish).
This function integrates these two workflows into a single display: 1. A pairwise environmental correlation matrix (Pearson or Spearman) rendered as a lower-triangular heatmap. 2. Community-wide Mantel tests connecting biological community nodes to corresponding environmental parameters using curved or straight linkage vectors.
Linkage attributes convey statistical properties: line thickness indicates the magnitude of the Mantel correlation coefficient (\(r\)), while color and line style denote permutation-derived significance thresholds (\(p < 0.01\), \(p < 0.05\), \(p \ge 0.05\)).
0.19.1 Conceptual Framework and Statistical Formulation
The analytical framework combines pairwise variable correlation with matrix permutation testing.
0.19.1.1 Mantel Matrix Correlation
Let \(\mathbf{D}_M\) represent the ecological dissimilarity matrix (e.g., Bray-Curtis) among \(n\) sampling stations for a biological community:
\[d_{ij}^{(M)} = \frac{\sum_{k=1}^p |y_{ik} - y_{jk}|}{\sum_{k=1}^p (y_{ik} + y_{jk})}\]
Let \(\mathbf{D}_E\) denote the standardized Euclidean distance matrix computed across one or more environmental variables between the same sampling stations:
\[d_{ij}^{(E)} = \sqrt{\sum_{m=1}^q (z_{im} - z_{jm})^2}\]
The Mantel test computes the Pearson correlation coefficient between the lower triangular elements (\(j < i\)) of the two unfolded distance matrices:
\[r_M = \frac{\sum_{i > j} (d_{ij}^{(M)} - \bar{d}_M)(d_{ij}^{(E)} - \bar{d}_E)}{\sqrt{\sum_{i > j} (d_{ij}^{(M)} - \bar{d}_M)^2 \sum_{i > j} (d_{ij}^{(E)} - \bar{d}_E)^2}}\]
where \(r_M \in [-1, 1]\).
0.19.1.2 Permutation Significance Testing
Because elements within a distance matrix are not mutually independent, standard parametric tests for correlation significance are invalid. Instead, significance is assessed through Monte Carlo permutation of the rows and columns of one distance matrix (e.g., \(\mathbf{D}_M\)) while holding the second (\(\mathbf{D}_E\)) fixed:
\[p = \frac{\sum_{\pi=1}^B \mathbb{I}(r_M^{(\pi)} \ge r_M^{\text{obs}}) + 1}{B + 1}\]
where \(B\) is the number of random permutations (default: 999).
0.19.2 In-Vignette Data Generation
We simulate a limnological survey across 35 sampling stations containing six continuous environmental variables (Water Temperature, Dissolved Oxygen, pH, Turbidity, Total Phosphorus, Total Nitrogen) paired with two biological trophic communities: Phytoplankton (6 taxa) and Zooplankton (5 taxa):
set.seed(1717)
n_stations <- 35
# 1. Environmental Gradients
temp_c <- rnorm(n_stations, mean = 24.0, sd = 2.5)
do_mgl <- 13.5 - (0.30 * temp_c) + rnorm(n_stations, 0, 0.4)
ph_val <- 7.2 + (0.04 * temp_c) + rnorm(n_stations, 0, 0.2)
turb_ntu <- runif(n_stations, 5, 40)
tp_mgl <- 0.03 + (0.007 * turb_ntu) + rnorm(n_stations, 0, 0.015)
tn_mgl <- 0.70 + (14.0 * tp_mgl) + rnorm(n_stations, 0, 0.12)
env_matrix <- data.frame(
Temp = temp_c,
DO = do_mgl,
pH = ph_val,
Turb = turb_ntu,
Total_P = tp_mgl,
Total_N = tn_mgl
)
# 2. Phytoplankton Community Counts (sensitive to nutrients and turbidity)
phyto_mat <- data.frame(
Cyanob_1 = rpois(n_stations, lambda = pmax(1, round(5 + 80 * tp_mgl))),
Cyanob_2 = rpois(n_stations, lambda = pmax(1, round(3 + 60 * tp_mgl))),
Diatom_1 = rpois(n_stations, lambda = pmax(1, round(28 - 0.4 * turb_ntu))),
Diatom_2 = rpois(n_stations, lambda = pmax(1, round(20 - 0.3 * turb_ntu))),
Chloroph_1 = rpois(n_stations, lambda = pmax(1, round(12 + 2.5 * ph_val))),
Chloroph_2 = rpois(n_stations, lambda = pmax(1, round(10 + 1.8 * ph_val)))
)
# 3. Zooplankton Community Counts (sensitive to temperature and dissolved oxygen)
zoo_mat <- data.frame(
Cladocera_1 = rpois(n_stations, lambda = pmax(1, round(2 + 2.2 * do_mgl))),
Cladocera_2 = rpois(n_stations, lambda = pmax(1, round(4 + 1.8 * do_mgl))),
Copepod_1 = rpois(n_stations, lambda = pmax(1, round(8 + 0.8 * temp_c))),
Copepod_2 = rpois(n_stations, lambda = pmax(1, round(6 + 0.6 * temp_c))),
Rotifer_1 = rpois(n_stations, lambda = pmax(1, round(15 + 0.3 * turb_ntu)))
)
# Bundle biological communities into a named list
aquatic_trophic_comms <- list(
"Phytoplankton" = phyto_mat,
"Zooplankton" = zoo_mat
)
head(env_matrix, 4)
#> Temp DO pH Turb Total_P Total_N
#> 1 20.36277 7.331757 7.902365 7.094070 0.07979408 2.036248
#> 2 19.72894 7.526729 8.056922 27.895527 0.20880550 3.590681
#> 3 29.25114 3.755476 8.352574 10.632814 0.09551539 1.887045
#> 4 21.42924 7.085186 8.155138 6.198184 0.08728681 2.0025460.19.3 Visualization Styles
The mantel_heatmap_analysis() function provides both an
integrated network linkage diagram and a side-by-side barplot
comparison.
0.19.3.1 Style A: Curved Linkage Network Diagram (Default)
Projects environmental correlations as lower-triangular tiles while connecting community nodes to environmental anchors with curved vectors:
res_mantel_net <- mantel_heatmap_analysis(
comm_data = aquatic_trophic_comms,
env_data = env_matrix,
plot_type = "network",
line_style = "curve",
method = "pearson",
spec_dist = "bray",
transform = "hellinger",
only_significant = FALSE,
line_color_by = "significance",
permutations = 999,
seed = 42,
show_cor_text = TRUE,
color_palette = "RdBu"
)
# Render publication network plot
res_mantel_net$Plot
0.19.3.2 Style B: Side-by-Side Correlation Heatmap and Mantel Bars
Displays the full environmental correlation matrix alongside horizontal bars showing Mantel \(r\) effect sizes faceted across trophic communities:
res_mantel_bars <- mantel_heatmap_analysis(
comm_data = aquatic_trophic_comms,
env_data = env_matrix,
plot_type = "bars",
method = "pearson",
spec_dist = "bray",
transform = "hellinger",
permutations = 999,
seed = 42,
color_palette = "RdBu"
)
res_mantel_bars$Plot
0.19.4 Model Inspection and Numerical Summaries
The returned list provides the pairwise environmental correlation matrix and a tidy data frame of community-specific Mantel statistics:
# 1. Environmental Correlation Matrix (First 4 Variables)
round(res_mantel_net$Correlation_Matrix[1:4, 1:4], 3)
#> Temp DO pH Turb
#> Temp 1.000 -0.926 0.623 0.151
#> DO -0.926 1.000 -0.593 -0.145
#> pH 0.623 -0.593 1.000 0.004
#> Turb 0.151 -0.145 0.004 1.000
# 2. Tidy Mantel Test Statistics and Permutation P-Values
head(res_mantel_net$Mantel_Results, 8)
#> Community Variable Mantel_r p_value Significance Mantel_R_Tier Line_Type
#> 1 Phytoplankton Temp 0.032 0.303 p ≥ 0.05 < 0.20 dashed
#> 2 Phytoplankton DO 0.074 0.177 p ≥ 0.05 < 0.20 dashed
#> 3 Phytoplankton pH 0.043 0.282 p ≥ 0.05 < 0.20 dashed
#> 4 Phytoplankton Turb 0.718 0.001 p < 0.01 ≥ 0.40 solid
#> 5 Phytoplankton Total_P 0.709 0.001 p < 0.01 ≥ 0.40 solid
#> 6 Phytoplankton Total_N 0.686 0.001 p < 0.01 ≥ 0.40 solid
#> 7 Zooplankton Temp 0.148 0.045 p < 0.05 < 0.20 solid
#> 8 Zooplankton DO 0.120 0.120 p ≥ 0.05 < 0.20 dashed0.20 Ecological Alpha Diversity Indices & Comparative Testing
The calc_diversity() and plot_diversity()
functions provide an integrated quantitative suite for estimating,
testing, and visualizing community alpha diversity patterns in fisheries
and aquatic ecology. Understanding local species richness, compositional
dominance, and structural evenness is foundational for assessing habitat
degradation, monitoring marine protected areas (MPAs), and comparing
stocking strategies across culture environments.
The package calculates seven primary ecological metrics simultaneously: 1. Total Abundance (\(N\)): Overall individual organism counts. 2. Species Richness (\(S\)): Total count of unique taxa observed. 3. Shannon-Wiener Index (\(H'\)): Information-theoretic entropy capturing both richness and proportional abundance. 4. Gini-Simpson Index (\(1 - D\)): The probability that two randomly sampled individuals belong to different species. 5. Margalef’s Richness (\(d\)): Species richness normalized relative to the natural logarithm of sample size. 6. Menhinick’s Index (\(D_{mn}\)): Species richness adjusted relative to the square root of total abundance. 7. Pielou’s Evenness (\(J'\)): The ratio of observed Shannon diversity relative to maximum possible diversity.
The visualization function plot_diversity() pairs these
indices with parametric (One-Way / Two-Way Factorial ANOVA) or
non-parametric (Kruskal-Wallis) omnibus tests, calculates
1,000-iteration bootstrap confidence intervals, and formats
publication-ready boxplots with raw observation overlays.
0.20.1 Conceptual Framework and Mathematical Formulation
Let a community sample contain \(S\) observed species and total organism abundance \(N = \sum_{i=1}^S n_i\), where \(n_i\) represents the count of the \(i\)-th taxon. The relative abundance (proportion) of species \(i\) is \(p_i = n_i / N\).
0.20.1.1 Shannon-Wiener Diversity Index (\(H'\))
Derived from communication theory, \(H'\) measures the uncertainty in predicting the species identity of a randomly sampled individual:
\[H' = -\sum_{i=1}^S p_i \ln(p_i)\]
Values typically range between \(1.5\) and \(3.5\) in natural aquatic systems, with higher values reflecting richer and more equitable communities.
0.20.1.2 Gini-Simpson Index (\(1 - D\))
Simpson’s index of dominance measures the probability that two randomly selected individuals belong to the same species: \(D = \sum_{i=1}^S p_i^2\). The complementary Gini-Simpson diversity metric is:
\[1 - D = 1 - \sum_{i=1}^S p_i^2\]
This metric ranges from \(0\) (monospecific culture) to \(1.0\) (infinite community of equally abundant species).
0.20.1.3 Species Richness Normalizations
Because raw observed richness (\(S\)) increases monotonically with sampling effort (\(N\)), Margalef (\(d\)) and Menhinick (\(D_{mn}\)) normalize richness across sample size discrepancies:
\[d = \frac{S - 1}{\ln(N)} \quad (N > 1)\]
\[D_{mn} = \frac{S}{\sqrt{N}} \quad (N > 0)\]
0.20.1.4 Pielou’s Equitability Index (\(J'\))
Measures how close the observed distribution of abundances is to maximum theoretical evenness (\(H'_{\max} = \ln(S)\)):
\[J' = \frac{H'}{\ln(S)} \quad (S > 1)\]
where \(J' \in [0, 1]\). Values approaching \(1.0\) indicate equal abundances across taxa, whereas lower values denote heavy dominance by one or two species.
0.20.2 In-Vignette Data Generation
We simulate a riverine fish biomonitoring survey assessing the
abundances of 10 freshwater species across 36 sampling reaches
categorized by River Zone (Headwater,
Mid-River, Lowland) and Sampling Season
(Dry Season, Wet Season):
set.seed(1818)
n_reps <- 6
zones <- c("Headwater", "Mid-River", "Lowland")
seasons <- c("Dry Season", "Wet Season")
design <- expand.grid(Zone = zones, Season = seasons)
comm_runs <- lapply(seq_len(nrow(design)), function(idx) {
zn <- as.character(design$Zone[idx])
sn <- as.character(design$Season[idx])
# Base abundance profiles
means <- switch(zn,
"Headwater" = c(Sp1 = 45, Sp2 = 38, Sp3 = 12, Sp4 = 2, Sp5 = 0, Sp6 = 0, Sp7 = 0, Sp8 = 0, Sp9 = 0, Sp10 = 0),
"Mid-River" = c(Sp1 = 18, Sp2 = 22, Sp3 = 25, Sp4 = 16, Sp5 = 14, Sp6 = 8, Sp7 = 4, Sp8 = 1, Sp9 = 0, Sp10 = 0),
"Lowland" = c(Sp1 = 6, Sp2 = 8, Sp3 = 12, Sp4 = 24, Sp5 = 28, Sp6 = 32, Sp7 = 18, Sp8 = 15, Sp9 = 10, Sp10 = 6)
)
if (sn == "Wet Season") {
means <- round(means * 1.35 + 2) # Seasonal flood pulse enhancement
}
mat <- do.call(cbind, lapply(means, function(m) rpois(n_reps, lambda = m)))
df <- as.data.frame(mat)
df$Zone <- zn
df$Season <- sn
df
})
fish_comm_data <- do.call(rbind, comm_runs)
# Organize grouping metadata columns at the front
fish_comm_data <- fish_comm_data[, c("Zone", "Season", paste0("Sp", 1:10))]
head(fish_comm_data[, 1:7], 6)
#> Zone Season Sp1 Sp2 Sp3 Sp4 Sp5
#> 1 Headwater Dry Season 38 39 16 0 0
#> 2 Headwater Dry Season 49 38 7 3 0
#> 3 Headwater Dry Season 49 47 14 2 0
#> 4 Headwater Dry Season 46 41 11 1 0
#> 5 Headwater Dry Season 56 47 17 4 0
#> 6 Headwater Dry Season 33 23 9 3 00.20.3 Index Calculation and Data Summary
First, calculate the complete alpha diversity table across all
sampling stations using calc_diversity():
alpha_table <- calc_diversity(
data = fish_comm_data,
group_col = c("Zone", "Season")
)
# Display sample output with computed diversity indices
head(alpha_table, 6)
#> Sample Zone Season Abundance Richness Shannon Simpson Margalef
#> 1 1 Headwater Dry Season 93 3 1.033 0.628 0.441
#> 2 2 Headwater Dry Season 97 4 1.009 0.585 0.656
#> 3 3 Headwater Dry Season 112 4 1.058 0.617 0.636
#> 4 4 Headwater Dry Season 99 4 1.012 0.600 0.653
#> 5 5 Headwater Dry Season 124 4 1.110 0.633 0.622
#> 6 6 Headwater Dry Season 68 4 1.123 0.631 0.711
#> Menhinick Pielou
#> 1 0.311 0.940
#> 2 0.406 0.728
#> 3 0.378 0.763
#> 4 0.402 0.730
#> 5 0.359 0.801
#> 6 0.485 0.8100.20.4 Visualization Styles
plot_diversity() generates single-factor and two-factor
(factorial) comparative displays.
0.20.4.1 Style A: Single-Factor Diversity Comparison with Kruskal-Wallis Test (Default)
Compares Shannon diversity across river zones, displaying raw sample jitter overlays and non-parametric test statistics:
res_div_single <- plot_diversity(
data = fish_comm_data,
group_col = "Zone",
index = "shannon",
test = "kruskal",
color_palette = "Dark2",
permutations = 1000
)
# Render single-factor diversity plot
res_div_single$Plot
0.20.4.2 Style B: Two-Factor Factorial Interaction Boxplot with Two-Way ANOVA
Evaluates how Pielou’s Evenness (\(J'\)) responds across the interaction of Zone and Season using dodged boxplots and parametric Two-Way ANOVA:
res_div_factorial <- plot_diversity(
data = fish_comm_data,
group_col = c("Zone", "Season"),
index = "pielou",
test = "anova",
color_palette = "Set2",
permutations = 1000
)
# Render factorial diversity plot
res_div_factorial$Plot
0.20.5 Model Inspection and Numerical Summaries
The returned list provides the underlying tidy diversity dataframe, bootstrap 95% confidence intervals, and hypothesis test output:
# 1. Omnibus Hypothesis Test Result (Kruskal-Wallis / ANOVA)
res_div_single$Test_Result
#>
#> Kruskal-Wallis rank sum test
#>
#> data: Metric_Value by Zone
#> Kruskal-Wallis chi-squared = 30.298, df = 2, p-value = 2.635e-07
# 2. Bootstrap Summary Statistics (Means, 95% Lower and Upper Confidence Limits)
res_div_single$Group_Summary
#> Group Index N Mean LCI UCI
#> 1 Headwater Shannon 12 1.242 1.139 1.348
#> 2 Lowland Shannon 12 2.138 2.108 2.160
#> 3 Mid-River Shannon 12 1.911 1.862 1.955
# 3. Two-Way ANOVA Factorial Summary Table
res_div_factorial$Test_Result
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Zone 2 0.29240 0.14620 97.19 7.80e-14 ***
#> Season 1 0.02987 0.02987 19.86 0.000108 ***
#> Zone:Season 2 0.04179 0.02090 13.89 5.37e-05 ***
#> Residuals 30 0.04513 0.00150
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 10.21 Species Accumulation Curves & Sampling Effort Evaluation
The specaccum_plot() function constructs and visualizes
species accumulation curves (collector’s curves) to evaluate inventory
completeness, sampling sufficiency, and species richness scaling in
aquatic biomonitoring and fisheries assessments. When conducting
community surveys (e.g., river electrofishing, benthic grab surveys,
plankton trawls), raw observed species counts are constrained by
sampling intensity.
Species accumulation curves model the expected rate at which new species are recorded as additional sampling units or hauls are aggregated. Asymptotic leveling of the curve signifies sampling saturation—demonstrating that the sampling protocol has captured the majority of common and rare taxa within the target habitat.
The function computes curves using vegan::specaccum,
accommodates multiple analytical accumulation methods (e.g.,
exact analytical solution, coleman
rarefaction, or random permutations), calculates 95% confidence
intervals, and generates publication-grade visualizations for single
pooled inventories or comparative multi-habitat groupings
(by_group = TRUE).
0.21.1 Conceptual Framework and Mathematical Formulation
Let a community dataset comprise \(N\) total sampling units and \(S_{\text{obs}}\) total observed species across the pooled inventory.
0.21.1.1 The Exact Analytical Expectation (Mao Tau Formulation)
Rather than relying solely on empirical randomizations,
method = "exact" calculates the expected species richness
\(\bar{S}_n\) for any sub-sample of
\(n\) sampling units (\(1 \le n \le N\)) using combinatorial
probability:
\[\bar{S}_n = S_{\text{obs}} - \sum_{i=1}^{S_{\text{obs}}} \frac{\binom{N - f_i}{n}}{\binom{N}{n}}\]
where \(f_i\) denotes the incidence frequency (number of sampling units containing species \(i\)). The combinatorial term represents the hyper-geometric probability that species \(i\) is absent from a random selection of \(n\) samples.
0.21.1.2 Coleman Rarefaction Formulation
Under method = "coleman", the non-interactive
expectation models individuals randomly distributed across samples
independently of spatial aggregation:
\[S_n = S_{\text{obs}} - \sum_{i=1}^{S_{\text{obs}}} \left(1 - \frac{n}{N}\right)^{n_i}\]
where \(n_i\) represents the total organism abundance of the \(i\)-th species.
0.21.1.3 Analytical Variance and Confidence Bounds
For methods with analytical variance derivations, the 95% confidence boundaries around expected richness \(\bar{S}_n\) are computed as:
\[\text{CI}_{95\%} = \bar{S}_n \pm 1.96 \times \sigma(\bar{S}_n)\]
where \(\sigma(\bar{S}_n)\) represents the standard deviation of species richness at sample size \(n\).
0.21.2 In-Vignette Data Generation
We simulate a riverine fish biodiversity survey comprising 45 electrofishing reaches divided across three distinct aquatic habitat zones: Highland Headwater (15 reaches), Midland Reach (15 reaches), and Lowland Floodplain (15 reaches). A total of 14 fish species with varying rarity and spatial distribution are surveyed:
set.seed(1919)
n_samples_per_hab <- 15
habitats <- c("Headwater", "Midland", "Lowland")
sim_list <- lapply(habitats, function(hab) {
# Diagnostic baseline probabilities of encounter
lambdas <- switch(hab,
"Headwater" = c(Sp01 = 25, Sp02 = 18, Sp03 = 8, Sp04 = 2, Sp05 = 0, Sp06 = 0, Sp07 = 0, Sp08 = 0, Sp09 = 0, Sp10 = 0, Sp11 = 0, Sp12 = 0, Sp13 = 0, Sp14 = 0),
"Midland" = c(Sp01 = 8, Sp02 = 12, Sp03 = 20, Sp04 = 15, Sp05 = 10, Sp06 = 6, Sp07 = 3, Sp08 = 1, Sp09 = 0, Sp10 = 0, Sp11 = 0, Sp12 = 0, Sp13 = 0, Sp14 = 0),
"Lowland" = c(Sp01 = 1, Sp02 = 2, Sp03 = 6, Sp04 = 14, Sp05 = 22, Sp06 = 28, Sp07 = 18, Sp08 = 12, Sp09 = 8, Sp10 = 6, Sp11 = 4, Sp12 = 3, Sp13 = 2, Sp14 = 1)
)
mat <- do.call(cbind, lapply(lambdas, function(l) rpois(n_samples_per_hab, lambda = l)))
df <- as.data.frame(mat)
df$Habitat <- hab
df
})
fish_inventory <- do.call(rbind, sim_list)
# Place metadata at the front
fish_inventory <- fish_inventory[, c("Habitat", setdiff(names(fish_inventory), "Habitat"))]
head(fish_inventory[, 1:7], 6)
#> Habitat Sp01 Sp02 Sp03 Sp04 Sp05 Sp06
#> 1 Headwater 24 16 11 5 0 0
#> 2 Headwater 23 12 8 1 0 0
#> 3 Headwater 22 10 12 2 0 0
#> 4 Headwater 29 11 8 1 0 0
#> 5 Headwater 31 21 6 1 0 0
#> 6 Headwater 25 14 5 2 0 00.21.3 Visualization Styles
The specaccum_plot() function supports both
whole-ecosystem pooled curves and multi-habitat comparative curves.
0.21.3.1 Style A: Pooled Whole-Ecosystem Curve with 95% Confidence Band (Default)
Evaluates overall survey completeness across all pooled sampling
units using the analytical exact method:
res_acc_pooled <- specaccum_plot(
data = fish_inventory,
group_col = "Habitat",
by_group = FALSE,
method = "exact",
line_color = "#2A9D8F",
ci = 1.96,
show_ci = TRUE,
title = "River Basin Fish Species Accumulation Curve"
)
# Render pooled curve
res_acc_pooled$Plot
0.21.3.2 Style B: Comparative Multi-Habitat Accumulation Curves
By specifying by_group = TRUE, the function estimates
accumulation curves separately for each habitat level, allowing direct
visual comparison of species richness saturation rates across
environmental strata:
res_acc_grouped <- specaccum_plot(
data = fish_inventory,
group_col = "Habitat",
by_group = TRUE,
method = "exact",
palette = "Dark2",
show_ci = TRUE,
title = "Comparative Species Accumulation Across River Zones"
)
res_acc_grouped$Plot
0.21.4 Model Inspection and Numerical Summaries
The returned list provides the underlying
vegan::specaccum object and the tidy data frame containing
exact site counts, expected richness, standard deviations, and
confidence intervals:
# 1. Inspect accumulation curve trajectory table (First 8 sample intervals)
head(res_acc_pooled$Curve_Data, 8)
#> Sites Richness SD Lower Upper
#> 1 1 8.088889 3.656484 0.9221799 15.25560
#> 2 2 10.459596 3.597171 3.4091399 17.51005
#> 3 3 11.723679 3.137314 5.5745433 17.87281
#> 4 4 12.487298 2.638540 7.3157596 17.65884
#> 5 5 12.978136 2.187067 8.6914857 17.26479
#> 6 6 13.304413 1.798105 9.7801273 16.82870
#> 7 7 13.525083 1.469828 10.6442192 16.40595
#> 8 8 13.675544 1.195670 11.3320315 16.01906
# 2. Inspect grouped curve summary table preview
head(res_acc_grouped$Curve_Data, 6)
#> Group Sites Richness SD Lower Upper
#> 1 Headwater 1 4 0 4 4
#> 2 Headwater 2 4 0 4 4
#> 3 Headwater 3 4 0 4 4
#> 4 Headwater 4 4 0 4 4
#> 5 Headwater 5 4 0 4 4
#> 6 Headwater 6 4 0 4 40.22 Weighted Arithmetic Water Quality Index (WAWQI)
The calc_wqi() and plot_wqi() functions
calculate and visualize the Weighted Arithmetic Water Quality Index
(WAWQI) across spatial monitoring stations and temporal sampling
campaigns. In fisheries management and aquatic environmental impact
assessments, raw multi-parameter water chemistry matrices (e.g., pH,
Dissolved Oxygen, Biological Oxygen Demand, Total Dissolved Solids,
Turbidity) can be difficult to interpret individually.
The WAWQI methodology aggregates multiple physico-chemical metrics into a single score that quantifies overall ecosystem water quality. The functions compute parameter-specific unit weights, derive quality ratings, handle station replicates with Student’s \(t\)-distribution 95% confidence intervals, assign ecological suitability tiers, and generate publication-ready bar and lollipop plots overlaid with regulatory benchmark bands.
0.22.1 Conceptual Framework and Mathematical Formulation
Let \(n\) denote the number of evaluated water quality parameters. Each parameter is evaluated against recommended regulatory thresholds (e.g., WHO, BIS, or national environmental standards).
0.22.1.1 Unit Weighting Factor (\(W_i\))
The unit weight \(W_i\) for the \(i\)-th parameter is inversely proportional to its standard permissible limit (\(S_i\)), ensuring that toxic or sensitive parameters (such as low dissolved oxygen thresholds or high pesticide levels) carry greater weight in the final index:
\[W_i = \frac{K}{S_i} = \frac{1}{S_i \sum_{j=1}^n \frac{1}{S_j}}\]
where \(K = \frac{1}{\sum_{j=1}^n (1 / S_j)}\) is the proportionality normalization constant, ensuring that \(\sum_{i=1}^n W_i = 1.0\).
0.22.1.2 Quality Rating Scale (\(Q_i\))
The sub-index rating \(Q_i\) for parameter \(i\) reflects the observed concentration (\(V_i\)) relative to its ideal standard value (\(V_0\)) and standard limit (\(S_i\)):
\[Q_i = \frac{V_i - V_0}{S_i - V_0} \times 100\]
For pollutants and nutrient parameters (e.g., BOD, TDS, Turbidity), the ideal zero-effect concentration is \(V_0 = 0\text{ mg/L}\). For physiological parameters with specific natural equilibria: - pH: \(V_0 = 7.0\) (neutral baseline) - Dissolved Oxygen (DO): \(V_0 = 14.6\text{ mg/L}\) (theoretical saturation baseline at standard temperature and pressure), formulated as: \[Q_{\text{DO}} = \frac{14.6 - V_{\text{DO}}}{14.6 - S_{\text{DO}}} \times 100\] Lower dissolved oxygen concentrations yield higher \(Q_i\) penalty scores.
0.22.1.3 Overall Index Aggregation and Classification
The aggregate WAWQI is computed as the weighted arithmetic mean of all sub-indices:
\[\text{WAWQI} = \frac{\sum_{i=1}^n Q_i W_i}{\sum_{i=1}^n W_i} = \sum_{i=1}^n Q_i W_i\]
Ecological water quality ratings are classified according to established benchmark boundaries: - \(0 \le \text{WAWQI} \le 25\): Excellent Water Quality (Pristine aquatic habitat) - \(26 \le \text{WAWQI} \le 50\): Good Water Quality (Suitable for aquaculture and fisheries) - \(51 \le \text{WAWQI} \le 75\): Poor Water Quality (Threatened aquatic biodiversity) - \(76 \le \text{WAWQI} \le 100\): Very Poor Water Quality (Severe eutrophication and habitat degradation) - \(\text{WAWQI} > 100\): Unsuitable for Aquatic Life (Acute toxicity / critical hypoxia)
0.22.2 In-Vignette Data Generation
We simulate a seasonal river biomonitoring dataset evaluating five
standard water quality parameters (pH, DO, BOD, TDS, and Turbidity)
across four monitoring stations (Station_1 to
Station_4), with 6 replicate water samples collected at
each station:
set.seed(2026)
n_reps <- 6
stations <- c("Station_1 (Upstream)", "Station_2 (Agricultural)", "Station_3 (Urban Discharge)", "Station_4 (Estuarine Mouth)")
wqi_survey <- do.call(rbind, lapply(stations, function(stn) {
# Parameter baseline values reflecting an anthropogenic pollution gradient
param_means <- switch(stn,
"Station_1 (Upstream)" = c(pH = 7.4, DO = 8.8, BOD = 1.8, TDS = 140, Turbidity = 2.5),
"Station_2 (Agricultural)" = c(pH = 7.8, DO = 6.8, BOD = 3.6, TDS = 320, Turbidity = 4.8),
"Station_3 (Urban Discharge)" = c(pH = 8.3, DO = 4.1, BOD = 8.5, TDS = 680, Turbidity = 9.4),
"Station_4 (Estuarine Mouth)" = c(pH = 7.9, DO = 5.4, BOD = 5.2, TDS = 510, Turbidity = 6.2)
)
data.frame(
Station = stn,
pH = rnorm(n_reps, mean = param_means["pH"], sd = 0.2),
DO = rnorm(n_reps, mean = param_means["DO"], sd = 0.35),
BOD = rnorm(n_reps, mean = param_means["BOD"], sd = 0.4),
TDS = rnorm(n_reps, mean = param_means["TDS"], sd = 25),
Turbidity = rnorm(n_reps, mean = param_means["Turbidity"], sd = 0.5)
)
}))
head(wqi_survey, 6)
#> Station pH DO BOD TDS Turbidity
#> 1 Station_1 (Upstream) 7.504118 8.542699 1.7114254 119.8820 2.524103
#> 2 Station_1 (Upstream) 7.184062 8.442957 1.7096734 157.2437 3.454042
#> 3 Station_1 (Upstream) 7.427848 8.839744 0.7812474 131.7832 3.365472
#> 4 Station_1 (Upstream) 7.383050 8.634173 2.3388006 135.8830 2.529073
#> 5 Station_1 (Upstream) 7.266672 8.657125 2.0465633 105.1993 2.822664
#> 6 Station_1 (Upstream) 6.896782 8.544348 1.8870257 176.6456 3.3628140.22.3 Calculating Index Scores and Station Summaries
Run calc_wqi() to derive individual replicate
sub-indices and station-level summary statistics:
res_wqi <- calc_wqi(
data = wqi_survey,
id_col = "Station",
standards = c(
pH = 8.5,
DO = 5.0,
BOD = 5.0,
TDS = 500,
Turbidity = 5.0
),
ci = 0.95
)
# Station-level summary table with 95% Confidence Intervals
res_wqi$Summary_Data
#> Station N WQI_Mean WQI_SD WQI_SE WQI_LCI WQI_UCI
#> Mean Station_1 (Upstream) 6 47.06 1.53 0.63 45.45 48.67
#> Mean1 Station_2 (Agricultural) 6 76.88 4.29 1.75 72.37 81.39
#> Mean2 Station_3 (Urban Discharge) 6 146.40 4.36 1.78 141.83 150.98
#> Mean3 Station_4 (Estuarine Mouth) 6 100.18 4.02 1.64 95.96 104.39
#> Status
#> Mean Good
#> Mean1 Very Poor
#> Mean2 Unsuitable
#> Mean3 Unsuitable0.22.4 Visualization Styles
plot_wqi() provides both columnar bar charts and
lollipop geometries with optional background benchmark classification
shading.
0.22.4.1 Style A: Status-Filled Bar Chart with Benchmark Bands and 95% CI Whiskers (Default)
Displays station scores against regulatory classification tiers, attaching 95% confidence intervals to convey measurement variability:
p_wqi_bar <- plot_wqi(
data = res_wqi,
id_col = "Station",
type = "bar",
show_ci = TRUE,
show_bands = TRUE,
title = "River Basin Water Quality Index (WAWQI) with 95% CI"
)
# Render publication plot
p_wqi_bar
0.22.4.2 Style B: Lollipop Display
Provides an alternative lollipop format without background shading, focusing on discrete station comparisons:
p_wqi_lollipop <- plot_wqi(
data = res_wqi,
id_col = "Station",
type = "lollipop",
show_ci = TRUE,
show_bands = FALSE,
title = "Station-Level Water Quality Ratings"
)
p_wqi_lollipop
0.22.5 Model Inspection and Replicate Data
Inspect individual replicate sub-indices (\(Q_i\)) and corresponding final sample scores:
# 1. Preview individual replicate evaluations (first 6 samples)
head(res_wqi$Replicate_Data, 6)
#> Station pH DO BOD TDS Turbidity Q_pH
#> 1 Station_1 (Upstream) 7.504118 8.542699 1.7114254 119.8820 2.524103 33.61
#> 2 Station_1 (Upstream) 7.184062 8.442957 1.7096734 157.2437 3.454042 12.27
#> 3 Station_1 (Upstream) 7.427848 8.839744 0.7812474 131.7832 3.365472 28.52
#> 4 Station_1 (Upstream) 7.383050 8.634173 2.3388006 135.8830 2.529073 25.54
#> 5 Station_1 (Upstream) 7.266672 8.657125 2.0465633 105.1993 2.822664 17.78
#> 6 Station_1 (Upstream) 6.896782 8.544348 1.8870257 176.6456 3.362814 0.00
#> Q_DO Q_BOD Q_TDS Q_Turbidity WQI
#> 1 63.10 34.23 23.98 50.48 46.64
#> 2 64.14 34.19 31.45 69.08 48.62
#> 3 60.00 15.62 26.36 67.31 44.46
#> 4 62.14 46.78 27.18 50.58 48.58
#> 5 61.90 40.93 21.04 56.45 47.23
#> 6 63.08 37.74 35.33 67.26 46.81
# 2. Inspect station-level summary statistics
res_wqi$Summary_Data[, c("Station", "N", "WQI_Mean", "WQI_SE", "Status")]
#> Station N WQI_Mean WQI_SE Status
#> Mean Station_1 (Upstream) 6 47.06 0.63 Good
#> Mean1 Station_2 (Agricultural) 6 76.88 1.75 Very Poor
#> Mean2 Station_3 (Urban Discharge) 6 146.40 1.78 Unsuitable
#> Mean3 Station_4 (Estuarine Mouth) 6 100.18 1.64 Unsuitable0.23 Heavy Metal Pollution & Ecotoxicological Indices
The calc_metal_indices() and
plot_metal_indices() functions evaluate heavy metal
pollution in aquatic ecosystems (water columns, riverbed sediments, and
fish tissue matrices). Aquatic toxicology and environmental impact
assessments require multi-element risk indices to convert raw trace
element concentrations into standardized toxicological risk scores.
The package computes five core international heavy metal contamination indices simultaneously: 1. Contamination Factor (\(CF_i\)): Ratio of observed metal concentration to its regulatory permissible limit. 2. Heavy Metal Pollution Index (HPI): Weighted arithmetic rating model balancing individual toxicities and guideline limits. 3. Heavy Metal Evaluation Index (HEI): Additive sum of contamination factors evaluating collective toxicity. 4. Degree of Contamination (\(C_d\)): Cumulative contamination degradation reflecting net excessive pollution above unity. 5. Pollution Load Index (PLI): Geometric mean of contamination factors quantifying multi-metal bioaccumulation. 6. Comprehensive Nemerow Pollution Index (\(P_N\)): Weighted combination highlighting peak extreme pollutant concentrations.
The visualization function plot_metal_indices() renders
column and lollipop graphics with 95% confidence intervals and dashed
benchmark thresholds.
0.23.1 Conceptual Framework and Mathematical Formulation
Let \(C_i\) denote the observed concentration of the \(i\)-th heavy metal, and let \(S_i\) denote its corresponding national or international drinking/aquatic life permissible standard limit.
0.23.1.1 Contamination Factor (\(CF_i\))
The base metric for single-metal contamination is the ratio:
\[CF_i = \frac{C_i}{S_i}\]
0.23.1.2 Heavy Metal Evaluation Index (HEI) & Contamination Degree (\(C_d\))
HEI expresses the cumulative additive load across all \(n\) analyzed trace metals:
\[\text{HEI} = \sum_{i=1}^n CF_i = \sum_{i=1}^n \frac{C_i}{S_i}\]
Benchmark classification: \(\text{HEI} < 10\) (Low); \(10 \le \text{HEI} \le 20\) (Medium); \(\text{HEI} > 20\) (High pollution).
The Degree of Contamination (\(C_d\)) computes the net excess above the permissible threshold (\(CF_i = 1\)):
\[C_d = \sum_{i=1}^n (CF_i - 1)\]
Benchmark classification: \(C_d < 6\) (Low); \(6 \le C_d < 12\) (Moderate); \(C_d \ge 12\) (High contamination).
0.23.1.3 Pollution Load Index (PLI)
Originally formulated by Tomlinson et al., PLI calculates the geometric mean of contamination factors:
\[\text{PLI} = \left( \prod_{i=1}^n CF_i \right)^{\frac{1}{n}} = \sqrt[n]{CF_1 \times CF_2 \times \dots \times CF_n}\]
Interpretation threshold: \(\text{PLI} \le 1.0\) indicates baseline unpolluted conditions, whereas \(\text{PLI} > 1.0\) indicates progressive heavy metal deterioration.
0.23.1.4 Comprehensive Nemerow Pollution Index (\(P_N\))
Evaluates total quality while giving special weight to the worst single pollutant:
\[P_N = \sqrt{\frac{\overline{CF}^2 + CF_{\max}^2}{2}}\]
where \(\overline{CF} = \frac{1}{n} \sum CF_i\) is the average contamination factor and \(CF_{\max}\) is the maximum observed single contamination factor.
0.23.1.5 Heavy Metal Pollution Index (HPI)
Computes a weighted arithmetic aggregation:
\[\text{HPI} = \frac{\sum_{i=1}^n W_i Q_i}{\sum_{i=1}^n W_i}\]
where unit weights are inversely proportional to guidelines (\(W_i = 1 / S_i\)), and sub-index ratings are \(Q_i = (C_i / S_i) \times 100\). The critical threshold is \(\text{HPI} = 100\); values exceeding \(100\) indicate critical metal pollution unsuitable for aquatic habitats.
0.23.2 In-Vignette Data Generation
We simulate a river basin heavy metal survey measuring six trace
elements (\(\mu\text{g/L}\)): Arsenic
(As), Cadmium (Cd), Chromium
(Cr), Copper (Cu), Lead (Pb), and
Zinc (Zn) across five monitoring reaches with 5 sample
replicates per reach:
set.seed(2121)
n_reps <- 5
sites <- c("Site_1 (Forest)", "Site_2 (Mining)", "Site_3 (Smelting)", "Site_4 (Urban)", "Site_5 (Estuary)")
metal_survey <- do.call(rbind, lapply(sites, function(stn) {
# Baseline metal concentrations (ug/L)
means <- switch(stn,
"Site_1 (Forest)" = c(As = 1.2, Cd = 0.2, Cr = 4.5, Cu = 45, Pb = 1.1, Zn = 120),
"Site_2 (Mining)" = c(As = 18.5, Cd = 5.8, Cr = 68.0, Cu = 1850, Pb = 24.5, Zn = 4200),
"Site_3 (Smelting)" = c(As = 26.0, Cd = 8.2, Cr = 95.0, Cu = 2600, Pb = 38.0, Zn = 5800),
"Site_4 (Urban)" = c(As = 6.5, Cd = 1.4, Cr = 28.0, Cu = 420, Pb = 8.5, Zn = 1450),
"Site_5 (Estuary)" = c(As = 4.2, Cd = 0.9, Cr = 16.5, Cu = 210, Pb = 4.2, Zn = 850)
)
data.frame(
Station = stn,
As = pmax(0.1, rnorm(n_reps, mean = means["As"], sd = means["As"] * 0.12)),
Cd = pmax(0.05, rnorm(n_reps, mean = means["Cd"], sd = means["Cd"] * 0.15)),
Cr = pmax(0.5, rnorm(n_reps, mean = means["Cr"], sd = means["Cr"] * 0.10)),
Cu = pmax(2.0, rnorm(n_reps, mean = means["Cu"], sd = means["Cu"] * 0.12)),
Pb = pmax(0.1, rnorm(n_reps, mean = means["Pb"], sd = means["Pb"] * 0.14)),
Zn = pmax(5.0, rnorm(n_reps, mean = means["Zn"], sd = means["Zn"] * 0.10))
)
}))
head(metal_survey, 6)
#> Station As Cd Cr Cu Pb Zn
#> 1 Site_1 (Forest) 1.085621 0.1844308 4.358769 45.40370 1.167778 132.3617
#> 2 Site_1 (Forest) 1.081720 0.1949841 4.678270 49.65574 1.101673 134.3707
#> 3 Site_1 (Forest) 1.209242 0.1466376 3.586198 46.37561 1.233846 124.5558
#> 4 Site_1 (Forest) 1.160399 0.1911586 4.680002 48.77459 1.350155 120.0001
#> 5 Site_1 (Forest) 1.220738 0.1577921 4.466422 41.96124 1.161292 118.0802
#> 6 Site_2 (Mining) 18.100377 4.5200358 71.436150 1746.83159 20.417828 4330.67340.23.3 Calculating Heavy Metal Indices
Run calc_metal_indices() using standard international
drinking and aquatic life guideline limits:
res_metals <- calc_metal_indices(
data = metal_survey,
id_col = "Station",
standards = c(
As = 10,
Cd = 3,
Cr = 50,
Cu = 2000,
Pb = 10,
Zn = 3000
),
ci = 0.95
)
# Station summary statistics table (Means, 95% CIs, and assigned status)
head(res_metals$Summary_Data[, c("Station", "HPI_Mean", "HPI_Status", "PLI_Mean", "PLI_Status")], 5)
#> Station HPI_Mean HPI_Status PLI_Mean PLI_Status
#> Mean Site_1 (Forest) 8.08 Suitable 0.06 Unpolluted
#> Mean1 Site_2 (Mining) 205.36 Polluted/Critical 1.62 Polluted
#> Mean2 Site_3 (Smelting) 295.52 Polluted/Critical 2.30 Polluted
#> Mean3 Site_4 (Urban) 56.92 Suitable 0.49 Unpolluted
#> Mean4 Site_5 (Estuary) 34.26 Suitable 0.29 Unpolluted0.23.4 Visualization Styles
plot_metal_indices() creates status-filled charts with
95% confidence intervals and dashed benchmark thresholds.
0.23.4.1 Style A: Heavy Metal Pollution Index (HPI) with Critical Threshold (Default)
Displays mean HPI values with 95% confidence intervals against the critical threshold line (\(\text{HPI} = 100\)):
p_hpi <- plot_metal_indices(
data = res_metals,
metric = "HPI",
id_col = "Station",
type = "bar",
show_ci = TRUE,
show_threshold = TRUE,
palette = "auto",
title = "Heavy Metal Pollution Index (HPI) Across River Stations"
)
# Render HPI plot
p_hpi
0.23.4.2 Style B: Pollution Load Index (PLI) Lollipop Plot
Visualizes the multi-metal Pollution Load Index (PLI) relative to the baseline threshold (\(\text{PLI} = 1.0\)):
p_pli <- plot_metal_indices(
data = res_metals,
metric = "PLI",
id_col = "Station",
type = "lollipop",
show_ci = TRUE,
show_threshold = TRUE,
palette = "auto",
title = "Pollution Load Index (PLI) Across Sampling Reaches"
)
p_pli
0.23.4.3 Style C: Degree of Contamination (Cd)
Evaluates cumulative excess contamination (\(C_d\)) with contamination limits (\(C_d = 6, 12\)):
p_cd <- plot_metal_indices(
data = res_metals,
metric = "Cd",
id_col = "Station",
type = "bar",
show_ci = TRUE,
show_threshold = TRUE,
palette = "auto"
)
p_cd
0.23.5 Model Inspection and Replicate Contamination Factors
Inspect sample-level contamination factors (\(CF_i\)) and overall indices:
# 1. Preview replicate contamination factors and computed indices (first 4 samples)
head(res_metals$Replicate_Data[, c("Station", "CF_As", "CF_Cd", "CF_Pb", "HPI", "PLI")], 4)
#> Station CF_As CF_Cd CF_Pb HPI PLI
#> 1 Site_1 (Forest) 0.11 0.06 0.12 8.08 0.06
#> 2 Site_1 (Forest) 0.11 0.06 0.11 8.19 0.07
#> 3 Site_1 (Forest) 0.12 0.05 0.12 7.61 0.06
#> 4 Site_1 (Forest) 0.12 0.06 0.14 8.71 0.07
# 2. Comprehensive Station Summary
res_metals$Summary_Data[, c("Station", "HEI_Mean", "HEI_Status", "Nemerow_PN_Mean", "Nemerow_Status")]
#> Station HEI_Mean HEI_Status Nemerow_PN_Mean Nemerow_Status
#> Mean Site_1 (Forest) 0.45 Low 0.10 Clean
#> Mean1 Site_2 (Mining) 10.14 Medium 2.14 Moderately Polluted
#> Mean2 Site_3 (Smelting) 14.49 Medium 3.11 Heavily Polluted
#> Mean3 Site_4 (Urban) 3.20 Low 0.73 Warning Limit
#> Mean4 Site_5 (Estuary) 1.85 Low 0.37 Clean0.24 Descriptive Statistics & Automated Table Summaries
The summarize_obs() function generates comprehensive
descriptive summaries for multi-parameter limnological, water quality,
and fisheries biometric matrices. Before applying formal parametric or
multivariate models, empirical researchers require rigorous tabular
summaries detailing central tendency, variance, standard errors, and
sample dispersion across spatial or experimental factor levels.
The function computes both parametric metrics (Sample Size \(N\), Arithmetic Mean \(\bar{x}\), Standard Deviation \(s\), Standard Error of the Mean \(\text{SE}\)) and non-parametric metrics
(Median \(\tilde{x}\), Minimum,
Maximum), constructs ready-to-publish composite notation (\(\bar{x} \pm s\)), accommodates multi-factor
grouping strata, and supports direct export to formatted Excel workbooks
(.xlsx) or comma-separated tables (.csv).
0.24.1 Conceptual Framework and Statistical Formulation
Let \(\{x_1, x_2, \dots, x_n\}\) denote a set of \(n\) complete observations recorded for a continuous aquatic variable \(X\) within a given sampling stratum.
0.24.1.1 Parametric Estimators of Central Tendency and Dispersion
Arithmetic Mean (\(\bar{x}\)): \[\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\]
Sample Standard Deviation (\(s\)): Measures sample dispersion around the mean using Bessel’s correction for unbiased estimation: \[s = \sqrt{\frac{1}{n - 1} \sum_{i=1}^n (x_i - \bar{x})^2} \quad (n > 1)\]
Standard Error of the Mean (\(\text{SE}\)): Quantifies the theoretical precision of the sample mean as an estimator of the true population parameter \(\mu\): \[\text{SE} = \frac{s}{\sqrt{n}}\]
0.24.1.2 Non-Parametric Order Statistics
For non-normal, skewed, or zero-inflated aquatic measurements (e.g., bacterial coliform densities, chemical trace contaminants):
Sample Median (\(\tilde{x}\)): The 50th percentile order statistic, robust against extreme outlier observations: \[\tilde{x} = \begin{cases} x_{\left(\frac{n+1}{2}\right)} & \text{if } n \text{ is odd} \\ \frac{1}{2}\left(x_{\left(\frac{n}{2}\right)} + x_{\left(\frac{n}{2} + 1\right)}\right) & \text{if } n \text{ is even} \end{cases}\]
Sample Extrema: \[\text{Min} = x_{(1)}, \quad \text{Max} = x_{(n)}\]
0.24.2 In-Vignette Data Generation
We simulate a river basin fisheries survey evaluating fish Total
Length (cm), Total Weight (g), and Condition Factor (\(K\)) across 60 specimens sampled across
three distinct River Systems (Karnafuli,
Halda, Sangu) during two hydrologic periods
(Dry, Wet):
set.seed(2222)
n_per_combo <- 10
systems <- c("Karnafuli", "Halda", "Sangu")
seasons <- c("Dry", "Wet")
survey_grid <- expand.grid(River = systems, Season = seasons)
survey_rows <- lapply(seq_len(nrow(survey_grid)), function(i) {
rv <- as.character(survey_grid$River[i])
sn <- as.character(survey_grid$Season[i])
mean_len <- switch(rv, "Karnafuli" = 28.5, "Halda" = 34.0, "Sangu" = 22.0)
if (sn == "Wet") mean_len <- mean_len * 1.15
len <- rnorm(n_per_combo, mean = mean_len, sd = 3.2)
# Allometric length-weight relationship: W = 0.012 * L^3.05
wt <- 0.012 * (len^3.05) + rnorm(n_per_combo, mean = 0, sd = 15)
k <- (wt / (len^3)) * 100
data.frame(
River = rv,
Season = sn,
Total_Length_cm = pmax(10, len),
Body_Weight_g = pmax(20, wt),
Condition_Factor = k,
stringsAsFactors = FALSE
)
})
fish_morph_survey <- do.call(rbind, survey_rows)
head(fish_morph_survey, 6)
#> River Season Total_Length_cm Body_Weight_g Condition_Factor
#> 1 Karnafuli Dry 27.41820 324.7617 1.575607
#> 2 Karnafuli Dry 31.50533 432.1249 1.381839
#> 3 Karnafuli Dry 34.06070 552.2470 1.397567
#> 4 Karnafuli Dry 30.72824 369.5970 1.273841
#> 5 Karnafuli Dry 29.97935 401.8828 1.491533
#> 6 Karnafuli Dry 27.49172 304.6220 1.4660730.24.3 Analytical Execution
summarize_obs() supports global overall summaries as
well as single-factor and multi-factor stratified analyses.
0.24.3.1 Style A: Single-Factor Grouped Summary (Across River Systems)
Summarizes length, weight, and condition factor across river systems:
res_summary_river <- summarize_obs(
data = fish_morph_survey,
vars = c("Total_Length_cm", "Body_Weight_g", "Condition_Factor"),
group = "River",
digits = 2
)
# Display tabular summary
res_summary_river[, c("River", "Variable", "N", "Mean_SD", "SE", "Median", "Min", "Max")]
#> River Variable N Mean_SD SE Median Min Max
#> 1 Halda Total_Length_cm 20 35.47 ± 3.74 0.84 34.54 30.10 43.62
#> 2 Karnafuli Total_Length_cm 20 31.87 ± 3.09 0.69 32.24 25.00 37.11
#> 3 Sangu Total_Length_cm 20 23.56 ± 4.16 0.93 23.28 17.32 31.51
#> 4 Halda Body_Weight_g 20 663.07 ± 227.50 50.87 600.84 396.03 1211.61
#> 5 Karnafuli Body_Weight_g 20 474.78 ± 129.99 29.07 480.26 228.20 730.16
#> 6 Sangu Body_Weight_g 20 195.32 ± 112.50 25.16 177.42 59.98 452.93
#> 7 Halda Condition_Factor 20 1.44 ± 0.04 0.01 1.44 1.34 1.49
#> 8 Karnafuli Condition_Factor 20 1.43 ± 0.06 0.01 1.44 1.27 1.58
#> 9 Sangu Condition_Factor 20 1.34 ± 0.13 0.03 1.38 0.99 1.500.24.3.2 Style B: Two-Factor Stratified Summary (River x Season Interaction)
Evaluates biometric traits across both River System and Season simultaneously:
res_summary_factorial <- summarize_obs(
data = fish_morph_survey,
vars = c("Total_Length_cm", "Body_Weight_g"),
group = c("River", "Season"),
digits = 2
)
head(res_summary_factorial[, c("River", "Season", "Variable", "Mean_SD", "SE", "Median")], 6)
#> River Season Variable Mean_SD SE Median
#> 1 Halda Dry Total_Length_cm 32.87 ± 1.58 0.50 33.18
#> 2 Karnafuli Dry Total_Length_cm 30.96 ± 2.52 0.80 31.12
#> 3 Sangu Dry Total_Length_cm 20.95 ± 2.50 0.79 20.82
#> 4 Halda Wet Total_Length_cm 38.07 ± 3.47 1.10 37.24
#> 5 Karnafuli Wet Total_Length_cm 32.78 ± 3.45 1.09 32.82
#> 6 Sangu Wet Total_Length_cm 26.17 ± 3.91 1.24 26.170.24.4 Inspection and Formatted Output
The returned data frame provides publication-ready string representations alongside raw numeric columns suitable for downstream reporting or plotting:
# 1. Structure of the exported summary table
str(res_summary_river)
#> 'data.frame': 9 obs. of 10 variables:
#> $ River : chr "Halda" "Karnafuli" "Sangu" "Halda" ...
#> $ Variable: chr "Total_Length_cm" "Total_Length_cm" "Total_Length_cm" "Body_Weight_g" ...
#> $ N : int 20 20 20 20 20 20 20 20 20
#> $ Mean : num 35.5 31.9 23.6 663.1 474.8 ...
#> $ SD : num 3.74 3.09 4.16 227.5 129.99 ...
#> $ SE : num 0.84 0.69 0.93 50.87 29.07 ...
#> $ Median : num 34.5 32.2 23.3 600.8 480.3 ...
#> $ Min : num 30.1 25 17.3 396 228.2 ...
#> $ Max : num 43.6 37.1 31.5 1211.6 730.2 ...
#> $ Mean_SD : chr "35.47 ± 3.74" "31.87 ± 3.09" "23.56 ± 4.16" "663.07 ± 227.50" ...
# 2. Inspect Mean +/- SD column for direct manuscript inclusion
res_summary_river[, c("River", "Variable", "Mean_SD")]
#> River Variable Mean_SD
#> 1 Halda Total_Length_cm 35.47 ± 3.74
#> 2 Karnafuli Total_Length_cm 31.87 ± 3.09
#> 3 Sangu Total_Length_cm 23.56 ± 4.16
#> 4 Halda Body_Weight_g 663.07 ± 227.50
#> 5 Karnafuli Body_Weight_g 474.78 ± 129.99
#> 6 Sangu Body_Weight_g 195.32 ± 112.50
#> 7 Halda Condition_Factor 1.44 ± 0.04
#> 8 Karnafuli Condition_Factor 1.43 ± 0.06
#> 9 Sangu Condition_Factor 1.34 ± 0.130.25 Distributional Histograms & Density Curves
The plot_histogram() function generates
publication-grade frequency distributions, empirical kernel density
estimates (KDE), and theoretical normal curves for fisheries
morphometrics and environmental water quality metrics. Inspecting
distributional skewness, modality (e.g., polymodal fish cohorts), and
kurtosis is a prerequisite for evaluating parametric model assumptions
(ANOVA, linear regression) and detecting recruitment pulses in
length-frequency data.
The function supports automatic dynamic scaling between raw counts and probability density space, overlays theoretical Gaussian profiles with empirical kernel fits, displays parametric mean (solid red line) and non-parametric median (dashed blue line) references, and supports multi-panel faceting across environmental strata.
0.25.1 Conceptual Framework and Statistical Formulation
Let \(\{x_1, x_2, \dots, x_n\}\) denote \(n\) independent observations of a continuous fisheries variable \(X\) (e.g., fork length or body weight).
0.25.1.1 Empirical Kernel Density Estimation (KDE)
To estimate the probability density function \(f(x)\) continuously without discretization artifacts, kernel smoothing is defined as:
\[\hat{f}_h(x) = \frac{1}{n h} \sum_{i=1}^n K\left(\frac{x - x_i}{h}\right)\]
where \(K(\cdot)\) is a standard Gaussian kernel function \(K(u) = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}u^2}\) and \(h > 0\) is the bandwidth smoothing parameter.
0.25.1.2 Theoretical Gaussian Overlay
When add_normal = TRUE is enabled, the function projects
the parameterized normal distribution curve evaluated at sample
estimates \((\bar{x}, s)\):
\[f(x; \bar{x}, s) = \frac{1}{s \sqrt{2\pi}} \exp\left( -\frac{(x - \bar{x})^2}{2 s^2} \right)\]
Discrepancies between the empirical kernel curve (red) and the theoretical Gaussian curve (green dashed) highlight skewness and kurtosis deviations from normality.
0.25.1.3 Central Tendency Reference Lines
The function evaluates asymmetry using the divergence between: - Sample Mean: \(\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\) (solid red line) - Sample Median: \(\tilde{x} = x_{(0.5)}\) (dashed blue line)
In right-skewed data (common in juvenile fish counts or pollutant spikes), \(\bar{x} > \tilde{x}\), whereas left-skewed distributions exhibit \(\bar{x} < \tilde{x}\).
0.25.2 In-Vignette Data Generation
We simulate a length-frequency fisheries survey measuring Total
Length (cm) across 250 specimens of an estuarine fish population sampled
across three estuarine habitats (Mangrove Creek,
Mudflat, Open Channel):
set.seed(2323)
n_fish <- 250
habitats <- sample(c("Mangrove Creek", "Mudflat", "Open Channel"), size = n_fish, replace = TRUE, prob = c(0.4, 0.35, 0.25))
# Simulate bimodal cohort growth
lengths <- vapply(habitats, function(hab) {
cohort <- rbinom(1, size = 1, prob = 0.65) # Age-0 vs Age-1 cohort
if (hab == "Mangrove Creek") {
if (cohort == 1) rnorm(1, 14.5, 1.8) else rnorm(1, 22.0, 2.1)
} else if (hab == "Mudflat") {
if (cohort == 1) rnorm(1, 18.2, 2.0) else rnorm(1, 26.5, 2.4)
} else {
if (cohort == 1) rnorm(1, 22.8, 2.2) else rnorm(1, 31.0, 2.8)
}
}, numeric(1))
estuary_fish <- data.frame(
Habitat = factor(habitats),
Total_Length = pmax(8.0, lengths)
)
head(estuary_fish, 6)
#> Habitat Total_Length
#> 1 Open Channel 29.48856
#> 2 Open Channel 20.06740
#> 3 Mangrove Creek 14.41802
#> 4 Mangrove Creek 26.79226
#> 5 Mudflat 19.81244
#> 6 Open Channel 27.198290.25.3 Visualization Styles
plot_histogram() provides single-distribution
diagnostics with distribution curves and comparative multi-panel faceted
views.
0.25.3.1 Style A: Overall Distribution with Density, Normal Fit & Reference Baselines
Overlays the empirical kernel density curve, theoretical Gaussian curve, and mean/median reference lines:
p_hist_single <- plot_histogram(
data = estuary_fish,
x = "Total_Length",
bins = 25,
add_density = TRUE,
add_normal = TRUE,
show_stats = TRUE,
xlab = "Total Length (cm)",
ylab = "Density",
title = "Fish Length-Frequency Distribution Diagnostics"
)
# Render distribution plot
p_hist_single
0.25.3.2 Style B: Faceted Distribution Across Habitats
Using facet = TRUE separates the length-frequency
distribution across habitats into separate panels:
p_hist_faceted <- plot_histogram(
data = estuary_fish,
x = "Total_Length",
group = "Habitat",
facet = TRUE,
bins = 20,
add_density = TRUE,
palette = "Dark2",
xlab = "Total Length (cm)",
title = "Comparative Length Distributions Across Habitats"
)
p_hist_faceted
0.25.4 Inspection of Distributional Symmetry
Examine the divergence between mean and median to evaluate skewness:
# 1. Compare central tendency benchmarks
mean_len <- mean(estuary_fish$Total_Length)
median_len <- median(estuary_fish$Total_Length)
sd_len <- sd(estuary_fish$Total_Length)
data.frame(
Metric = c("Mean (cm)", "Median (cm)", "Std. Deviation (cm)", "Skewness Offset"),
Value = round(c(mean_len, median_len, sd_len, mean_len - median_len), 3)
)
#> Metric Value
#> 1 Mean (cm) 21.033
#> 2 Median (cm) 20.625
#> 3 Std. Deviation (cm) 5.279
#> 4 Skewness Offset 0.4080.26 Publication Bar & Lollipop Plots with Grouped Comparisons
The plot_bar() and plot_lollipop()
functions render standardized bar and lollipop charts for categorical
inventories, survey frequency counts, and metric aggregations across
fisheries and aquatic ecology studies. Whether reporting catch
frequencies across fishing gear types, species occurrence counts across
sampling zones, or aggregated landings per river basin, these functions
provide direct comparative visualization.
The functions provide three operating modes: 1. Frequency
Tabulation (y = NULL): Computes categorical
occurrence counts automatically. 2. Metric Summary
Aggregation (stat = "mean" or "sum"):
Aggregates continuous numeric responses across single or multi-factor
grouping levels. 3. Identity Plotting
(stat = "identity"): Plots pre-computed values
directly.
In addition to standard bars (type = "bar"), lollipop
geometry (type = "lollipop" or via
plot_lollipop()) replaces solid bars with a slender stem
and circular head, reducing ink-to-data ratio while enhancing clarity
when comparing numerous categories. Both modes support side-by-side
grouped arrangements (position = "dodge"), cumulative
stacking (position = "stack"), 100% proportional fills
(position = "fill"), automated numeric value labelling,
coordinate flipping, and multi-column batch plotting.
0.26.1 Conceptual Framework and Aggregation Mechanics
Let a survey dataset contain observations indexed by \(i = 1, \dots, n\) across categorical predictor \(X\) with \(K\) unique levels.
0.26.1.1 Frequency Tabulation
When no response metric \(Y\) is supplied, the plotted height for factor level \(k\) is the sample frequency:
\[f_k = \sum_{i=1}^n \mathbb{I}(x_i = k)\]
0.26.1.2 Continuous Metric Aggregation
When numeric variable \(Y\) is
specified: - Mean Aggregation (stat = "mean"): \[\bar{y}_k = \frac{1}{n_k} \sum_{i: x_i = k}
y_i\] - Sum Aggregation (stat = "sum"): \[Y_k = \sum_{i: x_i = k} y_i\]
0.26.1.3 Proportional
Stacking (position = "fill")
When examining proportional composition across a secondary grouping factor \(G\) with \(M\) categories:
\[p_{km} = \frac{\sum_{i: x_i = k, g_i = m} y_i}{\sum_{m'=1}^M \sum_{i: x_i = k, g_i = m'} y_i}\]
The height of each stack normalizes to \(1.0\) (or 100%), allowing comparison of relative composition independently of differences in total sample size.
0.26.2 In-Vignette Data Generation
We simulate a marine small-scale fisheries monitoring survey
evaluating fish landings (kg per trip) and catch values across 90
fishing trips distributed among three Gear Types (Gillnet,
Longline, Trap) across three Coastal Sectors
(North Sector, Central Sector,
South Sector):
set.seed(2424)
n_trips <- 90
gears <- sample(c("Gillnet", "Longline", "Trap"), size = n_trips, replace = TRUE, prob = c(0.45, 0.35, 0.20))
sectors <- sample(c("North Sector", "Central Sector", "South Sector"), size = n_trips, replace = TRUE)
landings_kg <- vapply(seq_len(n_trips), function(i) {
g <- gears[i]
base <- switch(g, "Gillnet" = 120, "Longline" = 210, "Trap" = 65)
round(pmax(15, rnorm(1, mean = base, sd = base * 0.22)), 1)
}, numeric(1))
fisheries_catch_survey <- data.frame(
Gear_Type = factor(gears),
Sector = factor(sectors),
Landings_kg = landings_kg
)
head(fisheries_catch_survey, 6)
#> Gear_Type Sector Landings_kg
#> 1 Longline Central Sector 227.8
#> 2 Longline South Sector 225.8
#> 3 Gillnet North Sector 146.2
#> 4 Longline Central Sector 178.5
#> 5 Longline Central Sector 244.7
#> 6 Trap North Sector 75.60.26.3 Visualization Styles
0.26.3.1 Style A: Mean Aggregation with Grouped Dodge Layout and Value Labels (Default)
Calculates mean landings per trip across Gear Types grouped by Coastal Sector, placing numeric value labels above each bar:
p_bar_dodge <- plot_bar(
data = fisheries_catch_survey,
x = "Gear_Type",
y = "Landings_kg",
group = "Sector",
stat = "mean",
position = "dodge",
show_labels = TRUE,
palette = "Dark2",
title_prefix = "Small-Scale Fisheries"
)
p_bar_dodge
0.26.3.2 Style B:
Proportional 100% Stacked Bar Chart
(position = "fill")
Visualizes the relative proportion of fishing effort (trip count) per Gear Type across coastal sectors:
p_bar_fill <- plot_bar(
data = fisheries_catch_survey,
x = "Sector",
group = "Gear_Type",
position = "fill",
palette = "Set2",
ylab = "Proportion of Fishing Trips",
title_prefix = "Gear Composition"
)
p_bar_fill
0.26.3.3 Style C: Horizontal Bar Chart for Station Inventories
When category labels are long, horizontal = TRUE flips
coordinate axes to preserve legibility:
p_bar_horiz <- plot_bar(
data = fisheries_catch_survey,
x = "Gear_Type",
y = "Landings_kg",
stat = "sum",
horizontal = TRUE,
show_labels = TRUE,
palette = "Blues",
ylab = "Total Catch Landings (kg)"
)
p_bar_horiz
0.26.3.4 Style D: Clean Lollipop Plot for Total Catch Aggregations
Lollipop plots provide high contrast for ranked or discrete totals while eliminating visual clutter:
p_lollipop <- plot_lollipop(
data = fisheries_catch_survey,
x = "Gear_Type",
y = "Landings_kg",
stat = "sum",
point_size = 4.0,
line_size = 1.0,
show_labels = TRUE,
palette = "Set1",
ylab = "Total Catch Landings (kg)",
title_prefix = "Lollipop Summary"
)
p_lollipop
0.26.3.5 Style E: Multi-Factor Dodged Lollipop Plot
Compares mean landings across gears dodged by coastal sector using lollipop geometry:
p_lollipop_grouped <- plot_bar(
data = fisheries_catch_survey,
x = "Gear_Type",
y = "Landings_kg",
group = "Sector",
type = "lollipop",
stat = "mean",
point_size = 3.8,
line_size = 0.9,
palette = "Dark2",
title_prefix = "Dodged Lollipop Comparison"
)
p_lollipop_grouped
0.26.4 Inspection of Aggregated Values
# Summary table of Mean Landings per Gear and Sector
aggregate(Landings_kg ~ Gear_Type + Sector, data = fisheries_catch_survey, FUN = function(x) round(mean(x), 2))
#> Gear_Type Sector Landings_kg
#> 1 Gillnet Central Sector 126.12
#> 2 Longline Central Sector 227.89
#> 3 Trap Central Sector 70.24
#> 4 Gillnet North Sector 130.36
#> 5 Longline North Sector 262.84
#> 6 Trap North Sector 62.66
#> 7 Gillnet South Sector 104.42
#> 8 Longline South Sector 205.10
#> 9 Trap South Sector 62.730.27 Compositional Graphics & Survey Perception Analysis
In fisheries management, aquaculture surveys, and aquatic biomonitoring, reporting categorical proportions is standard practice. Whether portraying the percentage contribution of distinct species to total commercial landings, evaluating multi-tiered habitat classifications, or quantifying fisher community perceptions regarding marine protected areas (MPAs), standard plots often underperform in clarity or visual engagement.
The aquats package provides four specialized
compositional and perception graphics: 1. plot_pie3d():
Native isometric 3D pie charts with perspective compression, extruded
depth rims, and shadow rendering. 2. plot_circular_bar():
Radial bar charts with concentric reference gridlines and dynamic label
angles for many-category comparisons. 3. plot_donut():
Single-, double-, or triple-tiered nested sunburst/donut charts
representing hierarchical ecological structures. 4.
plot_likert(): Diverging and 100% stacked horizontal bar
charts designed for socio-economic survey evaluation.
0.27.1 3D Depth Pie
Charts (plot_pie3d)
The plot_pie3d() function renders an isometric 3D pie
chart directly in ggplot2 by calculating geometric polygon
coordinates for both the top circular face (compressed by perspective
factor \(\lambda \in [0.2, 0.8]\)) and
the extruded bottom rim:
\[x(\theta) = r \cos(\theta), \quad y(\theta) = r \sin(\theta) \times \lambda\]
The front rim is extruded over the interval \(\theta \in [\pi, 2\pi]\) (where \(y \le 0\)) with vertical depth \(\delta\):
\[y_{\text{rim\_bottom}}(\theta) = y(\theta) - \delta\]
We simulate commercial fisheries catch composition data across five key target species:
catch_composition <- data.frame(
Species = c("Hilsa Shad", "Rohu Carp", "Catla", "Prawn", "Catfish"),
Landings_MT = c(420, 210, 150, 95, 45)
)
p_pie <- plot_pie3d(
data = catch_composition,
x = "Species",
y = "Landings_MT",
stat = "identity",
tilt = 0.52,
depth = 0.22,
show_labels = TRUE,
min_percent = 5,
palette = "Set2",
title = "Annual Landings Composition by Target Taxa (MT)"
)
# Render 3D pie chart
p_pie
0.27.2 Circular Radial
Bar Plots (plot_circular_bar)
When a categorical variable contains eight or more levels (e.g.,
river stations, fishing gears, sub-districts), standard vertical bar
charts require horizontal label rotation that impairs readability.
plot_circular_bar() maps values along a polar coordinate
axis (coord_polar), reserving a central void fraction
(\(\rho\)) and rendering dashed
concentric reference rings with scale labels:
\[\theta_i = \frac{2\pi (i - 0.5)}{K}, \quad r_i = y_i\]
We simulate gear-specific catch per unit effort (CPUE in kg/haul) across 10 artisanal gear types:
gear_cpue <- data.frame(
Gear = c("Current Net", "Set Bag Net", "Longline", "Cast Net", "Gillnet",
"Push Net", "Drag Net", "Stake Net", "Traps", "Hook & Line"),
CPUE = c(42.5, 38.2, 29.0, 14.8, 24.5, 11.2, 18.6, 21.4, 8.5, 16.0)
)
p_circular <- plot_circular_bar(
data = gear_cpue,
x = "Gear",
y = "CPUE",
stat = "identity",
inner_radius = 0.35,
show_grid = TRUE,
palette = "Spectral",
title = "Artisanal Gear Catch per Unit Effort (kg/haul)"
)
# Render circular bar plot
p_circular
0.27.3 Hierarchical
Nested Donut Charts (plot_donut)
The plot_donut() function displays hierarchical
groupings (e.g., Habitat Zone \(\to\)
Feeding Guild \(\to\) Target Species)
across up to three concentric rings. Each layer aligns angular
boundaries (\(\theta_{\min},
\theta_{\max}\)) with parent categories while varying radial
boundaries (\(r_{\text{in}},
r_{\text{out}}\)):
\[\theta_{\text{start}}^{(k)} = \sum_{j=1}^{k-1} p_j \times 2\pi, \quad \theta_{\text{end}}^{(k)} = \sum_{j=1}^k p_j \times 2\pi\]
We simulate an estuarine survey categorizing fish biomass across Zone (Inner vs Outer Estuary) and Trophic Guild (Carnivore, Herbivore, Planktivore):
hierarchical_survey <- data.frame(
Zone = c("Inner", "Inner", "Inner", "Outer", "Outer", "Outer"),
Guild = c("Carnivore", "Herbivore", "Planktivore", "Carnivore", "Herbivore", "Planktivore"),
Biomass = c(120, 85, 145, 210, 60, 95)
)
p_donut <- plot_donut(
data = hierarchical_survey,
levels = c("Zone", "Guild"),
value = "Biomass",
show_labels = TRUE,
min_percent = 5,
palette = "Dark2",
title = "Hierarchical Estuarine Biomass (Zone > Trophic Guild)"
)
# Render nested donut chart
p_donut
0.27.4 Survey Perception
Likert Visualizations (plot_likert)
Socio-economic evaluations in fisheries (e.g., fisher compliance,
willingness to support seasonal bans, perception of catch declines) are
recorded on ordinal Likert scales. plot_likert() aligns
responses using either a diverging layout (centered on the neutral
mid-point) or a 100% stacked layout.
For an odd number of categories \(K\) with mid-point \(M = \lceil K/2 \rceil\), the diverging coordinate offset is:
\[x_{\text{start}} = -\left( \sum_{k=1}^{M-1} p_k + \frac{1}{2} p_M \right)\]
This splits the neutral category symmetrically across the vertical \(x = 0\) axis while anchoring negative responses (Disagree) to the left and positive responses (Agree) to the right.
We simulate a community perception survey evaluating fisher responses across four socio-economic items:
set.seed(2525)
n_fishers <- 100
likert_levels <- c("Strongly Disagree", "Disagree", "Neutral", "Agree", "Strongly Agree")
survey_perceptions <- data.frame(
Sanctuary_Support = sample(likert_levels, n_fishers, replace = TRUE, prob = c(0.10, 0.15, 0.10, 0.35, 0.30)),
Seasonal_Ban_Compliance = sample(likert_levels, n_fishers, replace = TRUE, prob = c(0.08, 0.12, 0.15, 0.40, 0.25)),
Catch_Decline_Observed = sample(likert_levels, n_fishers, replace = TRUE, prob = c(0.05, 0.08, 0.12, 0.35, 0.40)),
Subsidy_Sufficiency = sample(likert_levels, n_fishers, replace = TRUE, prob = c(0.35, 0.30, 0.15, 0.12, 0.08))
)
p_likert_diverging <- plot_likert(
data = survey_perceptions,
items = c("Sanctuary_Support", "Seasonal_Ban_Compliance", "Catch_Decline_Observed", "Subsidy_Sufficiency"),
levels = likert_levels,
type = "diverging",
palette = "RdYlBu",
show_labels = TRUE,
clean_labels = TRUE,
title = "Fisher Community Conservation Perception Survey (Diverging %)"
)
# Render diverging Likert plot
p_likert_diverging
0.28 Ecological Flow and Trophic Visualizations (Sankey Diagram)
0.28.1 Conceptual Framework
Visualizing interconnected pathways—such as energy transfer across trophic webs, ontogenetic dietary shifts across life stages, contaminant routing through environmental compartments, or fish migration across river catchments—requires representing both categorical transitions and volumetric flow magnitudes.
Alluvial and Sankey diagrams represent discrete categories as
vertically stacked blocks (strata) distributed across
consecutive parallel stages (axes), interconnected by
curvilinear ribbons (alluvia) whose heights are
proportional to flow frequency, biomass, or energy flux. The
plot_sankey() function provides a publication-ready
interface built on top of ggalluvial and
ggplot2, offering automated frequency aggregation, weighted
flows, missing-case filtering, and adaptive color palette
interpolation.
0.28.2 Flow Allocation and Stratum Architecture
The diagram decomposes relational data across \(k\) ordered categorical axes: 1.
Strata (Nodes): Represent discrete factor levels \(A_{i, j}\) at stage \(i\). 2. Alluvia (Flow
Ribbons): Trace the path taken by observations across all
stages \(i \in \{1, \ldots, k\}\). 3.
Weighting Metrics: * Unweighted Frequency
(\(N\)): When
value = NULL, ribbon thickness represents sample counts:
\[W = \sum 1\] * Weighted Flux
(\(V\)): When
value is assigned (e.g., ingested biomass \(\text{g}\), energy flow \(\text{t}\cdot\text{km}^{-2}\cdot\text{year}^{-1}\),
or contaminant concentration), ribbon thickness scales to cumulative
quantity: \[W = \sum V\]
0.28.3 In-Vignette Data Generation
We simulate an aquatic dietary study evaluating stomach content analysis across life stages, predatory feeding guilds, and preferred prey taxa:
set.seed(42)
diet_survey <- data.frame(
Life_Stage = factor(
sample(c("Juvenile", "Sub-Adult", "Adult"), 200, replace = TRUE, prob = c(0.40, 0.35, 0.25)),
levels = c("Juvenile", "Sub-Adult", "Adult")
),
Feeding_Guild = factor(
sample(c("Planktivore", "Benthivore", "Piscivore"), 200, replace = TRUE, prob = c(0.45, 0.35, 0.20)),
levels = c("Planktivore", "Benthivore", "Piscivore")
),
Prey_Taxa = factor(
sample(c("Cladocera", "Chironomidae", "Small Teleost", "Detritus"), 200, replace = TRUE),
levels = c("Cladocera", "Chironomidae", "Small Teleost", "Detritus")
),
Biomass_g = round(runif(200, min = 0.5, max = 18.5), 2)
)
head(diet_survey)
#> Life_Stage Feeding_Guild Prey_Taxa Biomass_g
#> 1 Adult Piscivore Detritus 9.47
#> 2 Adult Benthivore Detritus 5.59
#> 3 Juvenile Piscivore Detritus 14.48
#> 4 Adult Planktivore Cladocera 5.97
#> 5 Sub-Adult Planktivore Detritus 9.78
#> 6 Sub-Adult Planktivore Chironomidae 9.100.28.4 Visualization Styles
0.28.4.1 Style A: Biomass-Weighted Multi-Tier Alluvial
Flow ribbons are scaled to total prey biomass (\(\text{g}\)) and colored according to the
primary origin axis (Life_Stage):
p_weighted <- plot_sankey(
data = diet_survey,
axes = c("Life_Stage", "Feeding_Guild", "Prey_Taxa"),
value = "Biomass_g",
fill_axis = "Life_Stage",
alpha = 0.70,
width = 0.28,
palette = "Dark2",
title = "Style A: Biomass-Weighted Dietary Flux Across Life History Stages"
)
print(p_weighted)
0.28.4.2 Style B: Unweighted Observation Frequency (Count Flow)
When individual sample records or capture occurrences are tracked
without biomass weighting (value = NULL), ribbons reflect
observation frequencies and are colored by the intermediate axis
(Feeding_Guild):
p_counts <- plot_sankey(
data = diet_survey,
axes = c("Life_Stage", "Feeding_Guild", "Prey_Taxa"),
value = NULL, # Automatically computes record counts
fill_axis = "Feeding_Guild", # Colors ribbons by guild category
alpha = 0.60,
palette = "Set2",
title = "Style B: Specimen Frequency Pathways Across Trophic Strata"
)
print(p_counts)
0.28.4.3 Style C: Two-Tier Catchment-to-Estuary Transition
Sankey diagrams can also model spatial routing or migration between two distinct environmental zones:
# Simulate upstream-to-downstream fish passage counts
set.seed(101)
passage_data <- data.frame(
Release_Reach = sample(c("Headwater", "Mainstem Reservoir", "Tributary"), 150, replace = TRUE),
Recapture_Zone = sample(c("Estuary Lower", "Brackish Bay", "Nearshore Marine"), 150, replace = TRUE)
)
p_passage <- plot_sankey(
data = passage_data,
axes = c("Release_Reach", "Recapture_Zone"),
fill_axis = 1,
stratum_fill = "grey96",
stratum_color = "black",
label_size = 4.0,
palette = "Set1",
title = "Style C: Two-Tier Longitudinal Catchment Migration"
)
print(p_passage)
0.28.5 Customizing and Exporting Diagrams
Because plot_sankey() returns a native
ggplot2 object, it integrates with standard
ggplot2 layers and themes:
# Add custom axis adjustments and styling layers
p_custom <- p_weighted +
ggplot2::theme(
legend.position = "bottom",
plot.title = ggplot2::element_text(color = "navy", face = "bold", size = 14),
axis.title.y = ggplot2::element_text(size = 11, face = "bold")
) +
ggplot2::guides(fill = ggplot2::guide_legend(nrow = 1))
print(p_custom)
To export the diagram at publication resolution:
ggplot2::ggsave(
filename = "trophic_sankey_diagram.png",
plot = p_custom,
width = 9,
height = 5.5,
dpi = 300
)0.29 Temporal Tracking & Multi-Factor Time-Series Analysis
The plot_timeseries() function provides
publication-grade line tracking and longitudinal time-series
visualization for limnological water chemistry, estuarine hydrology, and
fisheries stock surveys. In long-term aquatic studies, monitoring
parameters across temporal gradients (e.g., daily tidal cycles, monthly
flood pulses, annual recruitment oscillations) requires clear
visualization of trends while accounting for within-season sampling
error.
The function handles: 1. Diverse Time-Step Ingestion: Parses date
formats (Date, POSIXct), continuous numeric
steps (e.g., Year: 2018–2026), or ordered calendar factors
("Jan", "Feb", "Mar" or seasonal
pulses: "Pre-Monsoon", "Monsoon",
"Post-Monsoon"). 2. Multi-Factor Hierarchical Partitioning:
- Single-Factor: Unstratified ecosystem temporal tracking. - Two-Factor:
Multi-station or inter-species trajectories differentiated by color
(color_var). - Three-Factor (Factorial Interaction):
High-dimensional tracking combining color, line styles
(linetype_var), and multi-panel faceting
(facet_var). 3. Automated Replicate Aggregation:
Automatically averages replicate samples per time-step and constructs
shaded uncertainty bands: 95% Confidence Intervals ("ci"),
Standard Error ("se"), or Standard Deviation
("sd").
0.29.1 Conceptual Framework & Uncertainty Band Mechanics
Let \(y_{t, i}\) denote the \(i\)-th sample replicate (\(i = 1, \dots, n_t\)) observed at temporal interval \(t\).
0.29.1.1 Replicate Arithmetic Mean
\[\bar{y}_t = \frac{1}{n_t} \sum_{i=1}^{n_t} y_{t, i}\]
0.29.1.2 Standard Error of the Mean (SE Ribbon)
Quantifies the precision of the estimated mean trajectory across repeated sampling visits: \[\text{SE}_t = \frac{s_t}{\sqrt{n_t}}, \quad \text{Ribbon}_{\text{SE}} = \bar{y}_t \pm \text{SE}_t\]
0.29.1.3 95% Student’s \(t\) Confidence Band (CI Ribbon)
Constructs the confidence boundary containing the true population parameter \(\mu_t\) at level \(\alpha = 0.05\): \[\text{CI}_{95\%, \, t} = \bar{y}_t \pm t_{n_t - 1, \, 0.975} \times \frac{s_t}{\sqrt{n_t}}\]
0.29.2 In-Vignette Data Generation
We simulate a 12-month estuarine monitoring survey tracking Dissolved
Oxygen (mg/L) across three River Stations (Upper River,
Middle Basin, Estuary Mouth) at two distinct
Water Depths (Surface vs. Bottom), collecting
4 replicate water column grabs per visit:
set.seed(2929)
months_seq <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
stations <- c("Upper River", "Middle Basin", "Estuary Mouth")
depths <- c("Surface", "Bottom")
n_reps <- 4
ts_grid <- expand.grid(Month = months_seq, Station = stations, Depth = depths)
ts_rows <- lapply(seq_len(nrow(ts_grid)), function(i) {
m <- as.character(ts_grid$Month[i])
st <- as.character(ts_grid$Station[i])
dp <- as.character(ts_grid$Depth[i])
# Base seasonal oxygen curve (higher in winter, lower during summer stratification)
m_idx <- match(m, months_seq)
base_do <- 7.5 + 2.2 * cos((m_idx - 1) * (2 * pi / 12))
# Spatial and depth adjustments
st_offset <- switch(st, "Upper River" = 0.8, "Middle Basin" = 0.0, "Estuary Mouth" = -0.6)
dp_offset <- if (dp == "Bottom") -1.4 else 0.0
sim_do <- pmax(2.0, rnorm(n_reps, mean = base_do + st_offset + dp_offset, sd = 0.35))
data.frame(
Month = factor(m, levels = months_seq),
Station = st,
Depth = dp,
DO_mgL = round(sim_do, 2),
stringsAsFactors = FALSE
)
})
aquatic_timeseries <- do.call(rbind, ts_rows)
head(aquatic_timeseries, 8)
#> Month Station Depth DO_mgL
#> 1 Jan Upper River Surface 10.63
#> 2 Jan Upper River Surface 10.49
#> 3 Jan Upper River Surface 10.59
#> 4 Jan Upper River Surface 10.39
#> 5 Feb Upper River Surface 9.88
#> 6 Feb Upper River Surface 9.65
#> 7 Feb Upper River Surface 9.64
#> 8 Feb Upper River Surface 10.590.29.3 Visualization Styles
plot_timeseries() provides single-series trajectories,
multi-line comparisons, and full three-factor factorial faceted
layouts.
0.29.3.1 Style A: Single-Series Temporal Trajectory with 95% CI Ribbon (Default)
Tracks overall monthly mean dissolved oxygen pooled across the watershed with 95% Confidence Interval whiskers:
res_ts_single <- plot_timeseries(
data = aquatic_timeseries,
time_var = "Month",
y_var = "DO_mgL",
ribbon = "ci",
point_size = 2.8,
line_width = 1.1,
ylab = "Dissolved Oxygen (mg/L)",
title = "Seasonal Dissolved Oxygen Dynamics (Mean \u00b1 95% CI)"
)
# Render single-series plot
res_ts_single$Plot
0.29.3.2 Style B: Two-Factor Comparison with Color-Coded Stations & SE Ribbons
Compares monthly trajectories across three river monitoring stations using color-coded lines and Standard Error (\(\pm 1\text{ SE}\)) ribbons:
res_ts_stations <- plot_timeseries(
data = aquatic_timeseries,
time_var = "Month",
y_var = "DO_mgL",
color_var = "Station",
ribbon = "se",
color_palette = "Dark2",
point_size = 2.4,
line_width = 1.0,
ylab = "Dissolved Oxygen (mg/L)",
title = "Comparative Station Dissolved Oxygen Trends"
)
res_ts_stations$Plot
0.29.3.3 Style C: Three-Factor Factorial Layout (Color + Linetype + Facet)
Demonstrates a high-dimensional layout tracking Month \(\times\) Station \(\times\) Depth: Stations are color-coded,
Water Depths are mapped to distinct line types (Surface
solid vs. Bottom dashed), and stations are faceted into
individual panels:
res_ts_factorial <- plot_timeseries(
data = aquatic_timeseries,
time_var = "Month",
y_var = "DO_mgL",
color_var = "Depth",
linetype_var = "Depth",
facet_var = "Station",
ribbon = "ci",
color_palette = "Set1",
point_size = 2.2,
ylab = "Dissolved Oxygen (mg/L)",
title = "Spatial and Bathymetric Dissolved Oxygen Profiles Across Seasons"
)
res_ts_factorial$Plot
0.29.4 Summary Table Inspection
The returned list contains the aggregated statistical table, detailing replicate counts (\(N\)), Means, Standard Deviations, Standard Errors, and Confidence Limits for each discrete temporal coordinate:
# 1. Preview aggregated summary statistics table (First 6 records)
head(res_ts_factorial$Summary_Data, 6)
#> Month Depth Depth.1 Station N Mean SD SE LCI
#> 1 Jan Bottom Bottom Estuary Mouth 4 7.7525 0.1992277 0.09961384 7.435484
#> 2 Feb Bottom Bottom Estuary Mouth 4 7.6200 0.2141650 0.10708252 7.279216
#> 3 Mar Bottom Bottom Estuary Mouth 4 6.4825 0.3854327 0.19271633 5.869191
#> 4 Apr Bottom Bottom Estuary Mouth 4 5.4100 0.6069596 0.30347982 4.444192
#> 5 May Bottom Bottom Estuary Mouth 4 4.3450 0.2272297 0.11361485 3.983427
#> 6 Jun Bottom Bottom Estuary Mouth 4 3.6150 0.3055596 0.15277980 3.128786
#> UCI
#> 1 8.069516
#> 2 7.960784
#> 3 7.095809
#> 4 6.375808
#> 5 4.706573
#> 6 4.101214