Fast Redundancy Analysis with fastrda

Introduction

fastrda is a high-performance implementation of Redundancy Analysis (RDA) written in C++ with Armadillo and OpenMP. It is designed for large-scale ecological, genomic, and other multivariate datasets where computational efficiency is critical.

Key Features:

High-performance C++ implementation: Uses Armadillo for linear algebra and OpenMP for parallel processing.
Partial RDA: Condition out covariates using efficient QR decomposition.
Permutation Tests: Parallel permutation tests via anova_fastrda() and the S3 anova() method, supporting both overall and axis-wise tests.
Comprehensive S3 Interface: Seamless integration with R standards (print(), summary(), predict(), scores(), plot(), and biplot()).
Flexible Workspace Modes: Choose between "minimal", "compact", "full", or "none" to balance memory and functionality.
Multiple Scaling Options: Four scaling types (0–3) for sample and response scores.
Biplot Visualization: Plots via biplotrda() or S3 generics plot() / biplot().
Memory Efficient: Uses scope-based RAII and avoids unnecessary copies.
Thread Control: Set OpenMP threads directly from R.

Installation

# Install from CRAN (when available)
install.packages("fastrda", dep=TRUE)

System Requirements:

OpenMP is required for parallel processing. On Linux/macOS this is usually pre-installed. On Windows, ensure you have RTools with OpenMP support.
BLAS library (e.g., OpenBLAS, MKL) is recommended for optimal performance.

Quick Start

To demonstrate fastrda, we will use the classic mite (oribatid mite) ecological dataset from the vegan package.

Load the package:

library(fastrda)

Prepare the data:

# Load real-world ecological data: Mite dataset
if (!requireNamespace("vegan", quietly = TRUE)) {
  stop("The 'vegan' package is required to run this example.")
}
data(mite, package = "vegan")
data(mite.env, package = "vegan")

# Apply Hellinger standardization to species abundance (response) data
Y <- vegan::decostand(mite, "hellinger")

# Select environmental variables (e.g., Substrate Density and Water Content)
# The intercept is removed because RDA expects only explanatory variables.
X <- model.matrix(~ SubsDens + WatrCont, mite.env)[, -1]

Basic Redundancy Analysis

Fit a standard RDA model using fastrda:

fit <- fastrda(
  response = Y,
  predictors = X,
  axes = 2,
  scaling = 2,
  keep_workspace = "minimal",
  threads = 1,  # Increase to use multiple CPU cores
  verbose = TRUE
)

Inspect the results using S3 methods:

# Inspect the components of model fit
names(fit)
#>  [1] "raw_sample_scores"     "sample_scores"         "response_scores"      
#>  [4] "scaling_constants"     "loadings"              "coefficients"         
#>  [7] "lc_coefficients_raw"   "eigenvalues"           "prop_total"           
#> [10] "prop_constrained"      "total_inertia"         "conditioned_inertia"  
#> [13] "constrained_inertia"   "unconstrained_inertia" "R2"                   
#> [16] "adj_R2"                "rank"                  "n_axes"               
#> [19] "n_samples"             "df_model"              "pseudo_F"             
#> [22] "df_resid"              "df_cond"               "X_rcond"              
#> [25] "X_rank_deficient"      "Z_rcond"               "Z_rank_deficient"     
#> [28] "scaling"               "X_center"              "X_scale"              
#> [31] "Y_center"              "Y_scale"               "Z_center"             
#> [34] "Z_scale"               "center"                "scale_x"              
#> [37] "scale_y"               "workspace"             "call"                 
#> [40] "pred_names"            "response_names"        "sample_names"         
#> [43] "site_names"

# Print concise model output
print(fit)
#> 
#> Call:
#> fastrda(response = Y, predictors = X, axes = 2, scaling = 2, 
#>     keep_workspace = "minimal", threads = 1, verbose = TRUE)
#> 
#>               Inertia Proportion Rank Df
#> Total          35.000     1.0000   NA 69
#> Constrained     8.477     0.2422    2  2
#> Unconstrained  26.523     0.7578   NA 67
#> 
#> Eigenvalues for constrained axes:
#>        [,1]
#> [1,] 7.3817
#> [2,] 1.0953
#> attr(,"names")
#> [1] "RDA1" "RDA2"
#> 
#> Pseudo-F     : 10.7070
#> Adjusted R2  : 0.2196

# Summary of ordination model
summary(fit)
#> 
#> Call:
#> fastrda(response = Y, predictors = X, axes = 2, scaling = 2, 
#>     keep_workspace = "minimal", threads = 1, verbose = TRUE)
#> 
#>               Inertia Proportion Rank Df
#> Total          35.000     1.0000   NA 69
#> Constrained     8.477     0.2422    2  2
#> Unconstrained  26.523     0.7578   NA 67
#> 
#> Eigenvalues for constrained axes:
#>        [,1]
#> [1,] 7.3817
#> [2,] 1.0953
#> attr(,"names")
#> [1] "RDA1" "RDA2"
#> 
#> Pseudo-F     : 10.7070
#> Adjusted R2  : 0.2196

# Eigenvalues
fit$eigenvalues
#>          [,1]
#> [1,] 7.381739
#> [2,] 1.095284
#> attr(,"names")
#> [1] "RDA1" "RDA2"

# R-squared values
fit$R2
#> [1] 0.2422007
fit$adj_R2
#> [1] 0.2195798

# Pre-computed sample scores
head(fit$sample_scores)
#>        RDA1        RDA2
#> 1 0.3758482  0.06010833
#> 2 0.4178976 -1.06204778
#> 3 0.4903470 -0.41542980
#> 4 0.6376417 -0.54602140
#> 5 0.7254019  1.22438141
#> 6 1.2759345 -1.10787286

# Pre-computed response scores
head(fit$response_scores)
#>               RDA1         RDA2
#> Brachy   0.1950925  0.191384927
#> PHTH     0.7445179  0.008595452
#> HPAV    -0.1803905  0.307530506
#> RARD     0.6724093  0.092191148
#> SSTR     0.4388376 -0.369165090
#> Protopl  0.4654368  0.000886920

# Canonical loadings (correlations)
head(fit$loadings)
#>                RDA1          RDA2
#> Brachy   0.06059908  0.1543295346
#> PHTH     0.23126000  0.0069312260
#> HPAV    -0.05603237  0.2479873453
#> RARD     0.20886184  0.0743413663
#> SSTR     0.13631046 -0.2976884206
#> Protopl  0.14457262  0.0007151972

Workspace Modes

fastrda offers four workspace modes to balance memory usage and functionality:

Mode	Stores	Recommended use
`"minimal"`	Q and QtY	Permutation testing (recommended default)
`"compact"`	QtY and R	Prediction for new data
`"full"`	X, Q, QtY, Y_res, R	Prediction and permutation testing
`"none"`	Nothing	Lowest memory usage (no permutation testing)

# Minimal workspace (default) - supports permutation tests
fit_min <- fastrda(response = Y, predictors = X, keep_workspace = "minimal", verbose = FALSE)

# Compact workspace - supports prediction with newdata
fit_compact <- fastrda(response = Y, predictors = X, keep_workspace = "compact", verbose = FALSE)

# Full workspace - supports everything
fit_full <- fastrda(response = Y, predictors = X, keep_workspace = "full", verbose = FALSE)

# No workspace - fastest, no permutation tests
fit_none <- fastrda(response = Y, predictors = X, keep_workspace = "none", verbose = FALSE)

Partial RDA (Conditioning on Covariates)

Often you want to remove the effect of known covariates (e.g., spatial structure, microhabitat type like Shrub) before testing predictor variables.

# Define conditioning matrix Z (e.g., Shrub presence/absence or type)
Z <- model.matrix(~ Shrub, mite.env)[, -1]

# Partial RDA: Y ~ SubsDens + WatrCont | Shrub
fit_partial <- fastrda(
  response = Y,
  predictors = X,
  covariates = Z,
  axes = 2,
  scaling = 2,
  keep_workspace = "minimal",
  threads = 1,  
  verbose = FALSE
)

print(fit_partial)
#> 
#> Call:
#> fastrda(response = Y, predictors = X, covariates = Z, axes = 2, 
#>     scaling = 2, keep_workspace = "minimal", threads = 1, verbose = FALSE)
#> 
#>               Inertia Proportion Rank Df
#> Total         35.0000     1.0000   NA 69
#> Conditioned    5.5976     0.1599    2  2
#> Constrained    4.4394     0.1268    2  2
#> Unconstrained 24.9629     0.7132   NA 65
#> 
#> Eigenvalues for constrained axes:
#>        [,1]
#> [1,] 3.3923
#> [2,] 1.0472
#> attr(,"names")
#> [1] "RDA1" "RDA2"
#> 
#> Pseudo-F     : 5.7798
#> Adjusted R2  : 0.1080

The conditioned inertia represents the variation explained exclusively by the conditioning variables before fitting the constrained model.

Permutation Tests

Assess the statistical significance of the constrained model using anova_fastrda() or the S3 generic anova().

Overall Model Test

# 999 permutations correspond to the minimum attainable p-value of 0.001.
# Larger numbers of permutations provide finer p-value resolution at the expense of longer computation time.

res_overall <- anova(fit, permutations = 999, threads = 1)
print(res_overall)
#> 
#> Permutation Test for Redundancy Analysis
#> Model: Constrained Inertia = 8.4770 (R2 = 0.2422)
#> Test Statistic (Pseudo-F): 10.7070
#> Degrees of Freedom: Model = 2, Residuals = 67
#> Permutations: 999
#> Pr(>F) / p-value: 0.0010 ***

Axis-wise Tests

# Test each constrained axis individually
res_axis <- fastrda::anova_fastrda(
  fit,
  by = "axis",
  permutations = 999,
  threads = 1
)
print(res_axis)
#> 
#> Marginal Permutation Test for Redundancy Analysis Axes
#> Permutations: 999
#> 
#>  Axis Eigenvalue Variance (%) Marginal p Sig
#>  RDA1     7.3817       87.08%     0.0010 ***
#>  RDA2     1.0953       12.92%     0.0010 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Disclaimer:
#>  Axis-wise permutation p-values quantify the extremeness of each constrained eigenvalue relative to its permutation distribution.

Axis-wise tests are marginal, not sequential. They should be interpreted descriptively rather than as formal significance tests for axis inclusion. Each canonical axis is evaluated against its own permutation distribution.

Scaling Options

The scaling parameter controls how sample and response scores are scaled:

Scaling	Sample Scores	Response Scores	Description
0	Unscaled	Unscaled	Raw scores
1	`s * const`	`/ const`	Samples scaled by `sqrt(eigenvalue)`
2	`* const`	`s / const`	Response scaled by `sqrt(eigenvalue)` (default)
3	`sqrt(s) * const`	`sqrt(s) / const`	Symmetric scaling

Where s = sqrt(eigenvalue) and const = ((n - 1) * total_inertia)^(0.25).

# Compare different scaling options
# Scaling = 2 is the default because it emphasizes response variable relationships and matches 
# the most commonly used scaling in ecological RDA applications.

fit0 <- fastrda(response = Y, predictors = X, scaling = 0, axes = 2, keep_workspace = "none", verbose = FALSE)
fit2 <- fastrda(response = Y, predictors = X, scaling = 2, axes = 2, keep_workspace = "none", verbose = FALSE)

head(fit2$sample_scores)
#>        RDA1        RDA2
#> 1 0.3758482  0.06010833
#> 2 0.4178976 -1.06204778
#> 3 0.4903470 -0.41542980
#> 4 0.6376417 -0.54602140
#> 5 0.7254019  1.22438141
#> 6 1.2759345 -1.10787286
head(fit2$response_scores)
#>               RDA1         RDA2
#> Brachy   0.1950925  0.191384927
#> PHTH     0.7445179  0.008595452
#> HPAV    -0.1803905  0.307530506
#> RARD     0.6724093  0.092191148
#> SSTR     0.4388376 -0.369165090
#> Protopl  0.4654368  0.000886920

Score Extraction

Extract sample, response, or predictor scores using standard S3 methods or dedicated helper functions:

# Extract both sample and response scores via S3 generic scores()
sc_both <- scores(fit, display = "both", choices = 1:2)
head(sc_both$samples)
#>        RDA1        RDA2
#> 1 0.3758482  0.06010833
#> 2 0.4178976 -1.06204778
#> 3 0.4903470 -0.41542980
#> 4 0.6376417 -0.54602140
#> 5 0.7254019  1.22438141
#> 6 1.2759345 -1.10787286
head(sc_both$response)
#>               RDA1         RDA2
#> Brachy   0.1950925  0.191384927
#> PHTH     0.7445179  0.008595452
#> HPAV    -0.1803905  0.307530506
#> RARD     0.6724093  0.092191148
#> SSTR     0.4388376 -0.369165090
#> Protopl  0.4654368  0.000886920

# Get sample scores with dynamic re-scaling (e.g., scaling = 1)
sample_sc1 <- get_sample_scores(fit, scaling = 1)
head(sample_sc1)
#>       RDA1        RDA2
#> 1 1.021155  0.06290687
#> 2 1.135401 -1.11149491
#> 3 1.332241 -0.43477150
#> 4 1.732431 -0.57144322
#> 5 1.970870  1.28138651
#> 6 3.466632 -1.15945353

# Get vegan-compatible response and predictor biplot scores
bp_scores <- biplot_scores(fit, type = "both", scaling = 2)
head(bp_scores$response)
#>               RDA1         RDA2
#> Brachy   0.1950925  0.191384927
#> PHTH     0.7445179  0.008595452
#> HPAV    -0.1803905  0.307530506
#> RARD     0.6724093  0.092191148
#> SSTR     0.4388376 -0.369165090
#> Protopl  0.4654368  0.000886920
head(bp_scores$predictors)
#>                RDA1        RDA2
#> SubsDens  0.1989174 -0.14000060
#> WatrCont -0.4103032 -0.02218149

Prediction

Use predict.fastrda() to get linear combination (LC) scores or fitted response values. The columns of newdata must have the same variables, order, and encoding as the predictor matrix used to fit the model.

# In-sample LC scores
lc_scores <- predict(fit, type = "lc")
head(lc_scores)
#>        RDA1        RDA2
#> 1 0.3758482  0.06010833
#> 2 0.4178976 -1.06204778
#> 3 0.4903470 -0.41542980
#> 4 0.6376417 -0.54602140
#> 5 0.7254019  1.22438141
#> 6 1.2759345 -1.10787286

# Predict for new predictor data
new_X <- head(X, 10)
new_lc <- predict(fit, newdata = new_X, type = "lc")
head(new_lc)
#>        RDA1       RDA2
#> 1  8.482354  0.5225437
#> 2  9.431348 -9.2327702
#> 3 11.066428 -3.6114834
#> 4 14.390656 -4.7467639
#> 5 16.371279 10.6439958
#> 6 28.796007 -9.6311444

# Reconstruct response matrix (type = "response")
pred_resp <- predict(fit, type = "response", rank = 2)
head(pred_resp[, 1:5])
#>      Brachy       PHTH      HPAV       RARD        SSTR
#> 1 0.2333160 0.06887523 0.2321238 0.06351154 0.019324967
#> 2 0.2048422 0.07041546 0.1887400 0.05790424 0.035560321
#> 3 0.2238497 0.07464611 0.2115936 0.06598406 0.027796892
#> 4 0.2243656 0.08236401 0.2033849 0.07247331 0.032030942
#> 5 0.2734324 0.08809205 0.2684112 0.08904840 0.008954043
#> 6 0.2267085 0.11581156 0.1679666 0.10062199 0.050322724

Biplot Visualization

fastrda provides publication-ready ggplot2-based biplots through direct calls to biplotrda() or S3 methods plot() and biplot(). Predictor arrows indicate the direction of increasing values for each explanatory variable, whereas sample and response positions summarize their relationships in canonical space.

Basic Biplot

# Create biplot with all components using S3 plot generic
plot(fit, axes = 1:2, scaling = 2, title = "Mite RDA Biplot (S3 method)")

Direct Function Call with biplotrda()

# Create biplot directly
biplotrda(fit, axes = 1:2, scaling = 2, title = "Mite RDA Biplot (Direct Call)")

Customizing Biplots

# Samples and predictors only
biplotrda(fit,
          type = "samples_predictors",
          sample_col = "darkblue",
          pred_col = "darkred",
          title = "Mite RDA - Samples and Predictors")

# Biplot with labels enabled
biplotrda(fit,
          show_ids = TRUE,
          max_labels = 20,
          title = "Labeled Mite RDA Biplot")

Performance Tuning

OpenMP Threads

Control the number of parallel threads used by C++ routines in fastrda:

# Use all cores (default)
fit <- fastrda(response = Y, predictors = X, threads = parallel::detectCores())

# Use 4 threads
fit <- fastrda(response = Y, predictors = X, threads = 4)

Internal benchmarks on synthetic datasets containing up to 10,000 response variables demonstrated median speedups of approximately 100×.

Help and Package Manual

For detailed documentation, argument descriptions, and additional examples for specific functions, you can access the built-in help pages directly in R:

help(package = "fastrda")

?fastrda
?anova.fastrda
?predict.fastrda
?scores.fastrda
?biplotrda

Summary

fastrda provides a fast, memory-efficient implementation of redundancy analysis for large ecological, genomic, and other multivariate datasets. The package supports partial RDA, permutation testing, prediction, multiple scaling options, and flexible visualization through a consistent S3 interface.

References

Legendre, P., & Gallagher, E. D. (2001). Ecologically meaningful transformations for ordination of species data. Oecologia, 129(2), 271–280. https://doi.org/10.1007/s004420100716
Legendre, P., & Legendre, L. (2012). Numerical ecology (3rd English ed.). Elsevier.
Oksanen, J., et al. (2022). vegan: Community ecology package (R package version 2.6-4) [Computer software]. https://CRAN.R-project.org/package=vegan
R Core Team. (2026). R: A language and environment for statistical computing. R Foundation for Statistical Computing.
ter Braak, C. J. F. (1994). Canonical community ordination. Part I: Basic theory and linear methods. Ecoscience, 1(2), 127–140. https://doi.org/10.1080/11956860.1994.11682237