Introduction_to

Statement of Need

In multivariate ecological and biological studies, researchers often rely on distance-based permutational methods like PERMANOVA (the \(F_1\) statistic) to detect differences among groups. However, standard PERMANOVA is strictly valid only when group dispersions are homogeneous. As demonstrated in Anderson et al. (2017), the \(F_1\) test is not reliable under heterogeneous dispersions in the space of the chosen dissimilarity measure, particularly for unbalanced designs:

Inflated Type I Error: When smaller groups have larger spreads, the \(F_1\) test becomes too liberal, falsely identifying differences in centroids.

Loss of Power: When larger groups have larger spreads, the \(F_1\) test becomes too conservative.

F2Adonis2 provides a mathematically grounded solution by modifying the PERMANOVA test statistic, ensuring valid inference even when groups exhibit heteroskedasticity in non-Euclidean space.

Methodology: The \(F_2\) Statistic

The \(F_2\) statistic addresses the BFP by replacing the pooled residual error term in the denominator with a weighted sum of group-specific mean-squared dissimilarities (\(V_i\)):

\[F_2 = \frac{SS_A}{\sum_{i=1}^g (1 - n_i/N)V_i}\]

Where \(SS_A\) is the among-group sum of squares, \(g\) is the number of groups, \(n_i\) is the sample size of group \(i\), \(N\) is the total sample size, and \(V_i\) is the group-specific variance calculated in dissimilarity space.

Resampling Methods

By default, \(p\)-values are computed by permuting group labels and recalculating the group-specific dispersions \(V_i\) for every permutation. The package also provides separate-sample residual bootstrapping (bootstrap = TRUE), which resamples residuals within each group in centered dissimilarity space so the groups share a common centroid, along with an empirical bias adjustment for the resulting bootstrap variances (bias.adjust = TRUE). Following Anderson et al. (2017), permutation is recommended for general applications, since it gave better Type I error control and equal or greater power than the bootstrap in their simulations.

Univariate Comparison

The test statistic provides a generalization of Welch’s F-test. The following code compares Welch’s F-test with the conventional ANOVA in a situation with no difference in the means but different variances.

################################################################################
# F2Adonis2 Vignette: Univariate Validation and Welch's F Comparison
################################################################################

library(ggplot2)
library(dplyr)
library(vegan)
#> Loading required package: permute
library(F2Adonis2)
library(onewaytests)

# ==============================================================================
# 1. Simulate three unbalanced groups, same means, different variances:
# ==============================================================================
set.seed(93939)
Xf <- factor(c(rep("Lev1", 15), rep("Lev2", 10), rep("Lev3", 5)))
y1 <- rnorm(n = 30, mean = rep(0, 30), sd = 20 * as.numeric(Xf))
df1 <- data.frame(y1, Xf)

# Visualize the heterogeneity (The Behrens-Fisher Problem)
df1 %>% ggplot(aes(x = Xf, y = y1, fill = Xf)) +
  geom_violin(alpha = 0.3) +
  geom_point() +
  theme_minimal() +
  labs(title = "Simulation: Same Means, Heterogeneous Variance")


# Standard ANOVA - Strong evidence (p=0.031) because it is biased by the variance
anova(lm(y1 ~ Xf, data = df1))
#> Analysis of Variance Table
#> 
#> Response: y1
#>           Df  Sum Sq Mean Sq F value  Pr(>F)  
#> Xf         2  8418.6  4209.3  3.9403 0.03152 *
#> Residuals 27 28843.5  1068.3                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ==============================================================================
# 2. Welch's F-test and Permutational Null
# ==============================================================================
welch_res <- oneway.test(y1 ~ Xf, data = df1)
print(welch_res)
#> 
#>  One-way analysis of means (not assuming equal variances)
#> 
#> data:  y1 and Xf
#> F = 4.0305, num df = 2.0000, denom df = 8.8085, p-value = 0.05716

get_welch_f <- function(y, groups) {
  oneway.test(y ~ groups)$statistic
}

# Generate Permutational Null Distribution for Welch's F
set.seed(123)
n_perm <- 99
welch_obs <- get_welch_f(y1, Xf)
welch_perms <- replicate(n_perm, {
  get_welch_f(y1, sample(Xf))
})

p_welch_perm <- (sum(welch_perms >= welch_obs) + 1) / (n_perm + 1)
cat("Permutational Welch's F p-value:", p_welch_perm, "\n")
#> Permutational Welch's F p-value: 0.04

The chunk below builds the F2 null distribution by calling F2_adonis2() once per permutation. This is not executed automatically, since 99 repeated calls (each re-parsing the formula and rebuilding the distance matrix) may increase the running time.

# ==============================================================================
# 3. F2 (F2Adonis2) Univariate Implementation & Manual Null Capture
# ==============================================================================
# Get Observed F2 - Using permutations = 1 to avoid subscript errors
set.seed(123)
fit_f2 <- F2_adonis2(y1 ~ Xf, data = df1, method = "euclidean", permutations = 1, plot.dist = FALSE)
f2_obs <- fit_f2["Model", "F2"]

# Manually capture the F2 null distribution by shuffling labels and calling
# F2_adonis2() with permutations = 1 for each shuffle
cat("Generating F2 null distribution...\n")
set.seed(123)
f2_perms <- numeric(n_perm)

for(i in 1:n_perm) {
  shuffled_df <- df1
  shuffled_df$Xf <- sample(df1$Xf)
  
  tmp <- F2_adonis2(y1 ~ Xf, data = shuffled_df, method = "euclidean", permutations = 1, plot.dist = FALSE)
  f2_perms[i] <- tmp["Model", "F2"]
}

# ==============================================================================
# 4. Compare Percentiles and Distributions
# ==============================================================================
perc_welch <- mean(welch_perms <= welch_obs, na.rm = TRUE)
perc_f2    <- mean(f2_perms <= f2_obs, na.rm = TRUE)

cat("\n--- Comparison of Observed Quantiles ---\n")
cat(sprintf("Welch's F Observed Percentile: %.4f\n", perc_welch))
cat(sprintf("F2 Statistic Observed Percentile: %.4f\n", perc_f2))

null_df <- data.frame(
  Value = c(welch_perms, f2_perms),
  Test = c(rep("Welch F", length(welch_perms)), rep("F2 Stat", length(f2_perms)))
)

ggplot(null_df, aes(x = Value, fill = Test)) +
  geom_density(alpha = 0.4) +
  geom_vline(xintercept = welch_obs, color = "blue", linetype = "dashed", linewidth = 1) +
  geom_vline(xintercept = f2_obs, color = "red", linetype = "dashed", linewidth = 1) +
  theme_minimal() +
  labs(title = "Permutational Distributions: Welch F vs F2",
       subtitle = "Dashed lines represent observed statistics",
       x = "Statistic Value", y = "Density")

Case Study 1: Norway Benthic Macrofauna

To demonstrate the effectiveness of the \(F_2\) statistic, we replicate the analysis of the Norway benthic macrofauna dataset (Ekofisk oil field) (Ellingsen and Gray 2002). This dataset consists of presence-absence data for 809 species across 101 sites, categorized into 5 Areas with unbalanced sample sizes.

Load Data and Libraries


library(F2Adonis2)

library(vegan)



# Load package datasets

data(norway_pa)

data(norway_metadata)



# Unbalanced design overview

table(norway_metadata$Area)
#> 
#>  1  2  3  4  5 
#> 16 21 25 19 20

Visualizing Dispersion Heterogeneity

Before running the test, we can confirm the presence of the Behrens-Fisher problem by visualizing the unequal group spreads using vegan::betadisper.

dist_jaccard <- vegdist(norway_pa, method = "jaccard")

mod_disp <- betadisper(dist_jaccard, norway_metadata$Area)

# Plotting dispersions
plot(mod_disp, main = "Multivariate Dispersions (Norway Data)")



library(yarrr)
#> Loading required package: jpeg
#> Loading required package: BayesFactor
#> Warning: package 'BayesFactor' was built under R version 4.5.3
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.8. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
#> 
#> Type BFManual() to open the manual.
#> ************
#> Loading required package: circlize
#> ========================================
#> circlize version 0.4.18
#> CRAN page: https://cran.r-project.org/package=circlize
#> Github page: https://github.com/jokergoo/circlize
#> Documentation: https://jokergoo.github.io/circlize_book/book/
#> 
#> If you use it in published research, please cite:
#> Gu, Z. circlize implements and enhances circular visualization
#>   in R. Bioinformatics 2014.
#> 
#> This message can be suppressed by:
#>   suppressPackageStartupMessages(library(circlize))
#> ========================================
#> 
#> Attaching package: 'yarrr'
#> The following object is masked from 'package:ggplot2':
#> 
#>     diamonds

# Extract data into a frame for plotting
disp_data <- data.frame(
  Distance = mod_disp$distances,
  Area = mod_disp$group
)

# Enhanced Pirate Plot
pirateplot(formula = Distance ~ Area,
           data = disp_data,
           main = "Multivariate Dispersions (Norway Data)",
           xlab = "Geographic Area",
           ylab = "Distance to Centroid",
           theme = 1,          # Clean theme
           pal = "basel",      # Professional palette
           point.o = 0.9,      # Transparency of raw points
           avg.line.o = 0.8,   # High visibility for mean
           bean.f.o = 0.2,     # Light density bean
           inf.method = "hdi") # High-density interval for inference

Global F2 PERMANOVA

We apply the F2_adonis2 function. Setting bootstrap = TRUE and bias.adjust = TRUE also reports the bootstrap and bias-adjusted bootstrap p-values alongside the permutation p-value. permutations = 99 is used here to reduce the running time; the accompanying JOSS paper reports results at permutations = 999.

set.seed(123)

# Conventional perMANOVA pseudo-F:

fit_results1 <- adonis2(

  formula = norway_pa ~ Area,

  data = norway_metadata,

  method = "jaccard",

  permutations = 99

)

fit_results1
#> Permutation test for adonis under reduced model
#> Permutation: free
#> Number of permutations: 99
#> 
#> adonis2(formula = norway_pa ~ Area, data = norway_metadata, permutations = 99, method = "jaccard")
#>           Df SumOfSqs      R2      F Pr(>F)   
#> Model      4   12.106 0.34934 12.886   0.01 **
#> Residual  96   22.548 0.65066                 
#> Total    100   34.654 1.00000                 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


# Running the F2 test

fit_results2 <- F2_adonis2(

  formula = norway_pa ~ Area,

  data = norway_metadata,

  method = "jaccard",

  permutations = 99,

  bootstrap = TRUE,

  bias.adjust = TRUE,

  plot.dist = FALSE

)



fit_results2
#> Permutation test for F2-based PERMANOVA
#>           Df SumOfSqs    F2 Pr(>F2) Pr(>F2.boot) Pr(>F2.ba)   
#> Model      4   12.106 13.51    0.01         0.01       0.01 **
#> Residual  96   22.548                                         
#> Total    100   34.654                                         
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Case Study 2: PKD Metabolomics

This example uses the urinary metabolomics matrix included in the package (95 samples, 1,554 metabolites), derived from the cohort in Houske et al. (2023).

# ==============================================================================
# High-Dimensional Case Study (PKD Metabolomics)
# ==============================================================================
# Demonstrating scaling to 1,554 variables where Welch's F is unavailable.
library(dplyr)
library(F2Adonis2)

set.seed(123)

cat("\n--- Loading PKD Data ---\n")
#> 
#> --- Loading PKD Data ---

# Load pre-packaged internal data objects
data("pkd_metadata", package = "F2Adonis2")
data("pkd_metabolites", package = "F2Adonis2")

group_labels <- factor(pkd_metadata$Group)
Y_pkd_final <- as.matrix(pkd_metabolites)

cat(sprintf("Final Processed Data: %d samples x %d metabolites\n", 
            nrow(Y_pkd_final), ncol(Y_pkd_final)))
#> Final Processed Data: 95 samples x 1554 metabolites

# ==============================================================================
# PERFORM F2 ANALYSIS
# ==============================================================================
cat("\n--- Running F2 PERMANOVA ---\n")
#> 
#> --- Running F2 PERMANOVA ---

pkd_meta <- data.frame(Group = group_labels)

# Run F2 with 99 permutations and bias-adjusted bootstrap.
# permutations = 99 is used here to keep the vignette fast to build; the
# accompanying JOSS paper reports results at permutations = 999.
pkd_f2_results <- F2_adonis2(
  formula = Y_pkd_final ~ Group,
  data = pkd_meta,
  method = "euclidean",
  permutations = 99,
  bootstrap = TRUE,
  bias.adjust = TRUE,
  plot.dist = FALSE
)

print(pkd_f2_results)
#> Permutation test for F2-based PERMANOVA
#>          Df SumOfSqs     F2 Pr(>F2) Pr(>F2.boot) Pr(>F2.ba)
#> Model     1     2211 1.4273    0.08         0.87       0.83
#> Residual 93   143865                                       
#> Total    94   146076

Implementation Details

The F2_adonis2 function is designed to be a “drop-in” replacement for vegan::adonis2, supporting standard R formulas and distance calculations via vegdist. Key features include:

Internal Consistency: Implements the \(V_i\) estimator directly from the distance matrix to handle high-dimensional compositional data.

Flexible Bootstrapping: Offers both standard permutation and separate-sample residual bootstrapping.

Bias Correction: Includes an empirical bias correction for bootstrap variances as derived in Anderson et al. (2017).

References

Anderson, M. J., Walsh, D. C., Clarke, K. R., Gorley, R. N., & Guerra-Castro, E. (2017). Some solutions to the multivariate Behrens–Fisher problem for dissimilarity-based analyses. Australian & New Zealand Journal of Statistics, 59(1), 57-79. doi:10.1111/anzs.12176

Dag, O., Dolgun, A., Konar, N. M. (2018). onewaytests: An R Package for One-Way Tests in Independent Groups Designs. The R Journal, 10(1), 175-199.

Ellingsen, K. E., & Gray, J. S. (2002). Spatial patterns of benthic diversity: is there a latitudinal gradient along the Norwegian continental shelf? Journal of Animal Ecology, 71(3), 373-389. doi:10.1046/j.1365-2656.2002.00606.x

Houske, E. A., Glimm, M. G., Bergstrom, A. R., Slipher, S. K., Welhaven, H. D., Greenwood, M. C., Linse, G. M., June, R. K., Yu, A. S. L., Wallace, D. P., & Hahn, A. K. (2023). Metabolomic profiling to identify early urinary biomarkers and metabolic pathway alterations in autosomal dominant polycystic kidney disease. American Journal of Physiology-Renal Physiology, 324(6), F590-F602. doi:10.1152/ajprenal.00301.2022

Welch, B. L. (1951). On the Comparison of Several Mean Values: An Alternative Approach. Biometrika, 38, 330-336.

Introduction_to_F2Adonis2

Lakviru Perera

Summary