## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(

  collapse = TRUE,

  comment = "#>",

  fig.width = 7,

  fig.height = 5

)

## ----setup--------------------------------------------------------------------
library(ggplot2)
library(dplyr)
library(F2Adonis2)

## -----------------------------------------------------------------------------
################################################################################
# F2Adonis2 Vignette: Univariate Validation and Welch's F Comparison
################################################################################

library(ggplot2)
library(dplyr)
library(vegan)
library(F2Adonis2)
library(onewaytests)

# ==============================================================================
# 1. Simulate three unbalanced groups, same means, different variances:
# ==============================================================================
set.seed(93939)
Xf <- factor(c(rep("Lev1", 15), rep("Lev2", 10), rep("Lev3", 5)))
y1 <- rnorm(n = 30, mean = rep(0, 30), sd = 20 * as.numeric(Xf))
df1 <- data.frame(y1, Xf)

# Visualize the heterogeneity (The Behrens-Fisher Problem)
df1 %>% ggplot(aes(x = Xf, y = y1, fill = Xf)) +
  geom_violin(alpha = 0.3) +
  geom_point() +
  theme_minimal() +
  labs(title = "Simulation: Same Means, Heterogeneous Variance")

# Standard ANOVA - Strong evidence (p=0.031) because it is biased by the variance
anova(lm(y1 ~ Xf, data = df1))

# ==============================================================================
# 2. Welch's F-test and Permutational Null
# ==============================================================================
welch_res <- oneway.test(y1 ~ Xf, data = df1)
print(welch_res)

get_welch_f <- function(y, groups) {
  oneway.test(y ~ groups)$statistic
}

# Generate Permutational Null Distribution for Welch's F
set.seed(123)
n_perm <- 99
welch_obs <- get_welch_f(y1, Xf)
welch_perms <- replicate(n_perm, {
  get_welch_f(y1, sample(Xf))
})

p_welch_perm <- (sum(welch_perms >= welch_obs) + 1) / (n_perm + 1)
cat("Permutational Welch's F p-value:", p_welch_perm, "\n")

## ----eval = FALSE-------------------------------------------------------------
# # ==============================================================================
# # 3. F2 (F2Adonis2) Univariate Implementation & Manual Null Capture
# # ==============================================================================
# # Get Observed F2 - Using permutations = 1 to avoid subscript errors
# set.seed(123)
# fit_f2 <- F2_adonis2(y1 ~ Xf, data = df1, method = "euclidean", permutations = 1, plot.dist = FALSE)
# f2_obs <- fit_f2["Model", "F2"]
# 
# # Manually capture the F2 null distribution by shuffling labels and calling
# # F2_adonis2() with permutations = 1 for each shuffle
# cat("Generating F2 null distribution...\n")
# set.seed(123)
# f2_perms <- numeric(n_perm)
# 
# for(i in 1:n_perm) {
#   shuffled_df <- df1
#   shuffled_df$Xf <- sample(df1$Xf)
# 
#   tmp <- F2_adonis2(y1 ~ Xf, data = shuffled_df, method = "euclidean", permutations = 1, plot.dist = FALSE)
#   f2_perms[i] <- tmp["Model", "F2"]
# }
# 
# # ==============================================================================
# # 4. Compare Percentiles and Distributions
# # ==============================================================================
# perc_welch <- mean(welch_perms <= welch_obs, na.rm = TRUE)
# perc_f2    <- mean(f2_perms <= f2_obs, na.rm = TRUE)
# 
# cat("\n--- Comparison of Observed Quantiles ---\n")
# cat(sprintf("Welch's F Observed Percentile: %.4f\n", perc_welch))
# cat(sprintf("F2 Statistic Observed Percentile: %.4f\n", perc_f2))
# 
# null_df <- data.frame(
#   Value = c(welch_perms, f2_perms),
#   Test = c(rep("Welch F", length(welch_perms)), rep("F2 Stat", length(f2_perms)))
# )
# 
# ggplot(null_df, aes(x = Value, fill = Test)) +
#   geom_density(alpha = 0.4) +
#   geom_vline(xintercept = welch_obs, color = "blue", linetype = "dashed", linewidth = 1) +
#   geom_vline(xintercept = f2_obs, color = "red", linetype = "dashed", linewidth = 1) +
#   theme_minimal() +
#   labs(title = "Permutational Distributions: Welch F vs F2",
#        subtitle = "Dashed lines represent observed statistics",
#        x = "Statistic Value", y = "Density")

## -----------------------------------------------------------------------------

library(F2Adonis2)

library(vegan)



# Load package datasets

data(norway_pa)

data(norway_metadata)



# Unbalanced design overview

table(norway_metadata$Area)

## -----------------------------------------------------------------------------
dist_jaccard <- vegdist(norway_pa, method = "jaccard")

mod_disp <- betadisper(dist_jaccard, norway_metadata$Area)

# Plotting dispersions
plot(mod_disp, main = "Multivariate Dispersions (Norway Data)")


library(yarrr)

# Extract data into a frame for plotting
disp_data <- data.frame(
  Distance = mod_disp$distances,
  Area = mod_disp$group
)

# Enhanced Pirate Plot
pirateplot(formula = Distance ~ Area,
           data = disp_data,
           main = "Multivariate Dispersions (Norway Data)",
           xlab = "Geographic Area",
           ylab = "Distance to Centroid",
           theme = 1,          # Clean theme
           pal = "basel",      # Professional palette
           point.o = 0.9,      # Transparency of raw points
           avg.line.o = 0.8,   # High visibility for mean
           bean.f.o = 0.2,     # Light density bean
           inf.method = "hdi") # High-density interval for inference

## -----------------------------------------------------------------------------
set.seed(123)

# Conventional perMANOVA pseudo-F:

fit_results1 <- adonis2(

  formula = norway_pa ~ Area,

  data = norway_metadata,

  method = "jaccard",

  permutations = 99

)

fit_results1


# Running the F2 test

fit_results2 <- F2_adonis2(

  formula = norway_pa ~ Area,

  data = norway_metadata,

  method = "jaccard",

  permutations = 99,

  bootstrap = TRUE,

  bias.adjust = TRUE,

  plot.dist = FALSE

)



fit_results2

## -----------------------------------------------------------------------------
# ==============================================================================
# High-Dimensional Case Study (PKD Metabolomics)
# ==============================================================================
# Demonstrating scaling to 1,554 variables where Welch's F is unavailable.
library(dplyr)
library(F2Adonis2)

set.seed(123)

cat("\n--- Loading PKD Data ---\n")

# Load pre-packaged internal data objects
data("pkd_metadata", package = "F2Adonis2")
data("pkd_metabolites", package = "F2Adonis2")

group_labels <- factor(pkd_metadata$Group)
Y_pkd_final <- as.matrix(pkd_metabolites)

cat(sprintf("Final Processed Data: %d samples x %d metabolites\n", 
            nrow(Y_pkd_final), ncol(Y_pkd_final)))

# ==============================================================================
# PERFORM F2 ANALYSIS
# ==============================================================================
cat("\n--- Running F2 PERMANOVA ---\n")

pkd_meta <- data.frame(Group = group_labels)

# Run F2 with 99 permutations and bias-adjusted bootstrap.
# permutations = 99 is used here to keep the vignette fast to build; the
# accompanying JOSS paper reports results at permutations = 999.
pkd_f2_results <- F2_adonis2(
  formula = Y_pkd_final ~ Group,
  data = pkd_meta,
  method = "euclidean",
  permutations = 99,
  bootstrap = TRUE,
  bias.adjust = TRUE,
  plot.dist = FALSE
)

print(pkd_f2_results)