---
title: "Comprehensive Ecological and Fisheries Data Analysis with aanova"
author: "Ataher Ali"
date: "2026-08-30"
output:
rmarkdown::html_vignette:
toc: true
toc_depth: 3
vignette: >
%\VignetteIndexEntry{Comprehensive Ecological and Fisheries Data Analysis with aanova}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
# Introduction
The `aanova` 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:
``` r
library(ggplot2)
library(dplyr)
library(emmeans)
library(multcomp)
library(MASS)
library(aanova)
1. Univariate and Factorial ANOVA Modules
One-Way ANOVA
When evaluating a single categorical driver against a continuous
response variable—such as comparing fish body weight across ecological
habitats—use one_way_anova().
res_one <- one_way_anova(
data = hilsa_weight,
factor_var = "Habitat",
numeric_var = "Weight_g",
factor_levels = c("River", "Estuary", "Marine"),
plot_type = "boxplot",
sig_display = "letters"
)
print(res_one$ANOVA_Summary)
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Habitat 2 965266 482633 43.96 6.39e-14 ***
#> Residuals 87 955142 10979
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(res_one$Plot)
Two-Way ANOVA
To examine how two independent factors and their interaction
influence stock weight parameters (e.g., Habitat combined with Fishing
Season), use two_way_anova().
res_two <- two_way_anova(
data = hilsa_two_way,
factor1_var = "Habitat",
factor2_var = "Season",
numeric_var = "Weight_g",
factor1_levels = c("River", "Estuary", "Marine"),
factor2_levels = c("Dry", "Monsoon"),
plot_type = "boxplot"
)
print(res_two$ANOVA_Summary)
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Habitat 2 1040321 520161 67.683 < 2e-16 ***
#> Season 1 335004 335004 43.590 3.45e-09 ***
#> Habitat:Season 2 16153 8076 1.051 0.354
#> Residuals 84 645564 7685
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(res_two$Plot)
Three-Way ANOVA with Faceting
For complex monitoring datasets spanning multiple life history
stages, three_way_anova() incorporates three categorical
factors and automatically generates faceted interaction plots with
compact letter displays.
res_three <- three_way_anova(
data = hilsa_three_way,
factor1_var = "Habitat",
factor2_var = "Season",
factor3_var = "Size_Class",
numeric_var = "Weight_g",
factor1_levels = c("River", "Estuary", "Marine"),
factor2_levels = c("Dry", "Monsoon"),
factor3_levels = c("Juvenile", "Adult"),
y_limits = c(0, 1600),
plot_type = "boxplot"
)
print(res_three$ANOVA_Summary)
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Habitat 2 1206114 603057 87.219 < 2e-16 ***
#> Season 1 361772 361772 52.322 7.16e-11 ***
#> Size_Class 1 5105022 5105022 738.329 < 2e-16 ***
#> Habitat:Season 2 12779 6389 0.924 0.400
#> Habitat:Size_Class 2 29225 14612 2.113 0.126
#> Season:Size_Class 1 1607 1607 0.232 0.631
#> Habitat:Season:Size_Class 2 18392 9196 1.330 0.269
#> Residuals 108 746744 6914
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(res_three$Plot)
2. Multivariate Analysis of Variance (MANOVA)
When analyzing multiple correlated morphometric traits simultaneously
(e.g., body depth, head length, and fin length across habitats),
standard univariate ANOVA is insufficient. The
manova_analysis() function handles multivariate testing and
supports both faceted trait boxplots and Canonical Discriminant Analysis
(LDA) ordination space.
res_lda <- manova_analysis(
data = hilsa_morphology,
response_vars = c("Body_Depth_cm", "Head_Length_cm", "Fin_Length_cm"),
factor_var = "Habitat",
factor_levels = c("River", "Estuary", "Marine"),
plot_type = "lda",
color_palette = "Set1"
)
print(res_lda$MANOVA_Summary)
#> Df Pillai approx F num Df den Df Pr(>F)
#> Habitat 2 0.93326 25.08 6 172 < 2.2e-16 ***
#> Residuals 87
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(res_lda$Plot)
3. Analysis of Covariance (ANCOVA)
To test group differences in a continuous response variable while
controlling for a continuous covariate (such as comparing fish weight
across habitats while adjusting for total length), use
ancova_analysis().
res_ancova <- ancova_analysis(
data = hilsa_ancova,
response_var = "Weight_g",
factor_var = "Habitat",
covariate_var = "Total_Length_cm",
factor_levels = c("River", "Estuary", "Marine"),
color_palette = "Set1"
)
print(res_ancova$ANCOVA_Table)
#> Analysis of Variance Table
#>
#> Response: Weight_g
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Total_Length_cm 1 7018460 7018460 12317.84 < 2.2e-16 ***
#> Habitat 2 233082 116541 204.54 < 2.2e-16 ***
#> Residuals 86 49001 570
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(res_ancova$Adjusted_Means)
#> Habitat emmean SE df lower.CL upper.CL
#> River 1164 4.45 86 1155 1173
#> Estuary 1223 4.36 86 1214 1232
#> Marine 1293 4.43 86 1284 1302
#>
#> Confidence level used: 0.95
print(res_ancova$Plot)
4. Multi-Model Regression Analysis
The regression_analysis() function fits linear,
logarithmic (log-log power curves for length-weight relationships), or
polynomial curves. It automatically computes \(R^2\) and \(p\)-values, generating publication-ready
plots with shaded confidence intervals.
res_reg <- regression_analysis(
data = hilsa_regression,
x_var = "Total_Length_cm",
y_var = "Weight_g",
fit_type = "logarithmic",
group_var = "Habitat",
color_palette = "Set1"
)
print(res_reg$Model_Summary)
#>
#> Call:
#> stats::lm(formula = mod_formula, data = data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.230993 -0.043236 0.004661 0.050214 0.203609
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -4.37808 0.12674 -34.543 <2e-16 ***
#> log_x 3.02537 0.03599 84.072 <2e-16 ***
#> HabitatMarine -0.33709 0.19267 -1.750 0.0835 .
#> HabitatRiver -0.19019 0.17345 -1.097 0.2757
#> log_x:HabitatMarine 0.10258 0.05610 1.829 0.0706 .
#> log_x:HabitatRiver 0.05638 0.04950 1.139 0.2576
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.074 on 94 degrees of freedom
#> Multiple R-squared: 0.9956, Adjusted R-squared: 0.9953
#> F-statistic: 4237 on 5 and 94 DF, p-value: < 2.2e-16
print(res_reg$Plot)
5. Generalized Linear Models (GLMs) for Count Data
Fisheries count data (such as catch-per-unit-effort or abundance
counts) often violate normal distribution assumptions. The
glm_analysis() function fits Poisson, Quasipoisson, or
Binomial models and outputs a forest plot of Incident Rate Ratios (IRR)
or Odds Ratios.
res_glm <- glm_analysis(
data = hilsa_catch,
response_var = "Catch_Count",
predictor_vars = c("Habitat", "Season", "Fishing_Hours"),
family_type = "quasipoisson",
color_palette = "Set1"
)
print(res_glm$Model_Summary)
#>
#> Call:
#> stats::glm(formula = formula_str, family = fam_obj, data = data)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.861524 0.061177 46.774 < 2e-16 ***
#> HabitatMarine 0.482481 0.034343 14.049 < 2e-16 ***
#> HabitatRiver -0.335685 0.041633 -8.063 7.87e-13 ***
#> SeasonMonsoon 0.545593 0.030374 17.962 < 2e-16 ***
#> Fishing_Hours 0.041905 0.006377 6.571 1.52e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for quasipoisson family taken to be 0.966984)
#>
#> Null deviance: 1051.39 on 119 degrees of freedom
#> Residual deviance: 110.45 on 115 degrees of freedom
#> AIC: NA
#>
#> Number of Fisher Scoring iterations: 4
print(res_glm$Plot)
6. Correlation Matrices and Custom Heatmaps
To explore collinearity among environmental parameters,
correlation_heatmap() computes Pearson, Spearman, or
Kendall correlation matrices and renders customizable circle or square
heatmaps with triangle layout options.
res_cor <- correlation_heatmap(
data = hilsa_env,
method = "pearson",
shape = "circle",
view = "lower",
color_palette = "RdBu"
)
print(res_cor$Correlation_Matrix)
#> SST_C Salinity_ppt Depth_m DO_mgL Chlorophyll_a
#> SST_C 1.00000000 NA NA NA NA
#> Salinity_ppt 0.03127984 1.000000000 NA NA NA
#> Depth_m -0.14477734 0.071226992 1.00000000 NA NA
#> DO_mgL 0.07423731 0.008644817 -0.04597967 1.00000000 NA
#> Chlorophyll_a 0.08427619 -0.120722592 0.05933732 0.18501354 1.0000000
#> Catch_kg -0.04228272 -0.198932998 0.10053673 0.06940787 0.1524039
#> Catch_kg
#> SST_C NA
#> Salinity_ppt NA
#> Depth_m NA
#> DO_mgL NA
#> Chlorophyll_a NA
#> Catch_kg 1
print(res_cor$Plot)
7. Integrated Mantel Test & Correlation Heatmap
The mantel_heatmap_analysis() function combines internal
environmental correlation matrices with community Mantel test linkages
into a publication-grade network-heatmap layout.
set.seed(42)
env_test <- data.frame(
pH = rnorm(40, 7.5, 0.4),
DO = rnorm(40, 6.2, 0.7),
Temp = rnorm(40, 27.5, 1.5),
Salinity = rnorm(40, 14.0, 2.5)
)
comm_test <- data.frame(
Taxon_A = rpois(40, 12),
Taxon_B = rpois(40, 18),
Taxon_C = rpois(40, 7)
)
res_mantel <- mantel_heatmap_analysis(
comm_data = comm_test,
env_data = env_test,
method = "pearson",
color_palette = "RdBu"
)
print(res_mantel$Plot)