---
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.

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 site and species 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")

# Or install development version from GitHub
# devtools::install_github("zcebeci/fastrda")

```

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 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(
  genotype = Y,
  environment = 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 site scores (samples)
head(fit$site_scores)

# Pre-computed species scores (response variables)
head(fit$species_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(Y, X, keep_workspace = "minimal", verbose = FALSE)

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

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

# No workspace - fastest, no permutation tests
fit_none <- fastrda(Y, 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 environmental 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(
  genotype = Y,
  environment = 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 <- 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 site and species scores are scaled:

| Scaling | Site Scores | Species Scores | Description |
| --- | --- | --- | --- |
| 0 | Unscaled | Unscaled | Raw scores |
| 1 | `s * const` | `/ const` | Sites scaled by `sqrt(eigenvalue)` |
| 2 | `* const` | `s / const` | Species 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 species relationships and matches the most 
# commonly used scaling in ecological RDA applications.

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

head(fit2$site_scores)
head(fit2$species_scores)

```

---

## Score Extraction

Extract site, species, or environmental scores using the standard S3 `scores()` method, `get_site_scores()`, or `biplot_scores()`:

```{r}
# Extract both site and species scores via S3 generic scores()
sc_both <- scores(fit, display = "both", choices = 1:2)
head(sc_both$sites)
head(sc_both$species)

# Get site scores with dynamic re-scaling (e.g., scaling = 1)
site_sc1 <- get_site_scores(fit, scaling = 1)
head(site_sc1)

# Get vegan-compatible species and environmental biplot scores
bp_scores <- biplot_scores(fit, type = "both", scaling = 2)
head(bp_scores$species)
head(bp_scores$environment)

```

---

## 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 environmental matrix used to fit the model.

```{r}
# In-sample LC scores
lc_scores <- predict(fit, type = "lc")
head(lc_scores)

# Predict for new environmental 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()`.
Environmental arrows indicate the direction of increasing values for each explanatory variable, whereas site and species 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}
# Sites and environment only
biplotrda(fit,
          type = "sites_environment",
          site_col = "darkblue",
          env_col = "darkred",
          title = "Mite RDA - Sites and Environment")

```

```{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(Y, X, threads = parallel::detectCores())

# Use 4 threads
fit <- fastrda(Y, 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

```