---
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.
<img src="../man/figures/fastrdalogo.png" alt="fastrda" width="200" style="float: right; margin-left: 15px;"/>

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