## ----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
)


## ----loadlib------------------------------------------------------------------
library(fastrda)


## ----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]


## -----------------------------------------------------------------------------
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 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)


## -----------------------------------------------------------------------------
# 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)


## -----------------------------------------------------------------------------
# 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)


## -----------------------------------------------------------------------------
# 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)


## -----------------------------------------------------------------------------
# Test each constrained axis individually
res_axis <- anova_fastrda(
  fit,
  by = "axis",
  permutations = 999,
  threads = 1
)
print(res_axis)


## -----------------------------------------------------------------------------
# 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)


## -----------------------------------------------------------------------------
# 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)


## -----------------------------------------------------------------------------
# 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])


## ----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)")


## ----fig.width=7, fig.height=5------------------------------------------------
# Create biplot directly
biplotrda(fit, axes = 1:2, scaling = 2, title = "Mite RDA Biplot (Direct Call)")


## ----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")


## ----fig.width=7, fig.height=5------------------------------------------------
# Biplot with labels enabled
biplotrda(fit,
          show_ids = TRUE,
          max_labels = 20,
          title = "Labeled Mite RDA Biplot")


## ----eval=FALSE---------------------------------------------------------------
# # Use all cores (default)
# fit <- fastrda(Y, X, threads = parallel::detectCores())
# 
# # Use 4 threads
# fit <- fastrda(Y, X, threads = 4)
# 

## ----eval=FALSE---------------------------------------------------------------
# help(package = "fastrda")
# 
# ?fastrda
# ?anova.fastrda
# ?predict.fastrda
# ?scores.fastrda
# ?biplotrda
#