--- title: "Fast Redundancy Analysis with fastrda" author: "Zeynel Cebeci" date: "`r Sys.Date()`" output: rmarkdown::html_document: theme: flatly highlight: tango toc: true toc_float: true toc_depth: 2 number_sections: false code_folding: show df_print: paged fig_width: 7 fig_height: 5 vignette: > %\VignetteIndexEntry{Fast Redundancy Analysis with fastrda} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = TRUE, cache = FALSE, fig.width = 7, fig.height = 5, fig.align = "center", message = FALSE, warning = FALSE ) ``` ## 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. fastrda

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 ```r # 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: ```{r loadlib} library(fastrda) ``` Prepare the data: ```{r, message=FALSE} # 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`: ```{r} 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: ```{r} # Inspect the components of model fit names(fit) # Print concise model output print(fit) # Summary of ordination model summary(fit) # Eigenvalues fit$eigenvalues # R-squared values fit$R2 fit$adj_R2 # Pre-computed sample scores head(fit$sample_scores) # Pre-computed response scores head(fit$response_scores) # Canonical loadings (correlations) head(fit$loadings) ``` --- ## 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) | ```{r} # 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. ```{r} # 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) ``` 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 ```{r} # 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) ``` ### Axis-wise Tests ```{r} # Test each constrained axis individually res_axis <- fastrda::anova_fastrda( fit, by = "axis", permutations = 999, threads = 1 ) print(res_axis) ``` > 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)`. ```{r} # 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) head(fit2$response_scores) ``` --- ## Score Extraction Extract sample, response, or predictor scores using standard S3 methods or dedicated helper functions: ```{r} # Extract both sample and response scores via S3 generic scores() sc_both <- scores(fit, display = "both", choices = 1:2) head(sc_both$samples) head(sc_both$response) # Get sample scores with dynamic re-scaling (e.g., scaling = 1) sample_sc1 <- get_sample_scores(fit, scaling = 1) head(sample_sc1) # Get vegan-compatible response and predictor biplot scores bp_scores <- biplot_scores(fit, type = "both", scaling = 2) head(bp_scores$response) head(bp_scores$predictors) ``` --- ## 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. ```{r} # In-sample LC scores lc_scores <- predict(fit, type = "lc") head(lc_scores) # Predict for new predictor data new_X <- head(X, 10) new_lc <- predict(fit, newdata = new_X, type = "lc") head(new_lc) # Reconstruct response matrix (type = "response") pred_resp <- predict(fit, type = "response", rank = 2) head(pred_resp[, 1:5]) ``` --- ## 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 ```{r, fig.width=7, fig.height=5} # 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() ```{r, fig.width=7, fig.height=5} # Create biplot directly biplotrda(fit, axes = 1:2, scaling = 2, title = "Mite RDA Biplot (Direct Call)") ``` ### Customizing Biplots ```{r, fig.width=7, fig.height=5} # Samples and predictors only biplotrda(fit, type = "samples_predictors", sample_col = "darkblue", pred_col = "darkred", title = "Mite RDA - Samples and Predictors") ``` ```{r, fig.width=7, fig.height=5} # 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`: ```{r, eval=FALSE} # 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: ```{r, eval=FALSE} 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