## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## -----------------------------------------------------------------------------
library(dplyr)
library(vegan)
library(F2Adonis2)
library(kableExtra)
library(ggplot2)
library(combinat)

# Load package data objects explicitly
data("pkd_metadata", package = "F2Adonis2")
data("pkd_metabolites", package = "F2Adonis2")

set.seed(527)

# ==============================================================================
# 1. Sub-sample the PKD dataset (10 Ctrl / 20 PKD)
# ==============================================================================
sub_metadata <- bind_rows(
  pkd_metadata %>% filter(Group == "Ctrl") %>% slice_sample(n = 10),
  pkd_metadata %>% filter(Group == "PKD") %>% slice_sample(n = 20)
)
sub_metabolites <- pkd_metabolites[rownames(sub_metadata), ]
dist_sub <- vegdist(sub_metabolites, method = "euclidean")

# ==============================================================================
# 2. Dispersion diagnostics (betadisper)
# ==============================================================================
mod_sub <- betadisper(dist_sub, sub_metadata$Group)
anova_disp <- permutest(mod_sub, permutations = 99)
print(anova_disp)

# Spatial hull plot
plot(mod_sub, hull = TRUE, ellipse = FALSE,
     main = "Spatial Dispersion: Sub-sample (10 Ctrl / 20 PKD)")

# Violin + boxplot + jittered points of distance-to-centroid
distances_df <- data.frame(
  distance = mod_sub$distances,
  group = mod_sub$group
)

ggplot(distances_df, aes(x = group, y = distance, fill = group)) +
  geom_violin(alpha = 0.2, trim = TRUE, color = NA) +
  geom_boxplot(width = 0.2, alpha = 0.6, outlier.shape = NA) +
  geom_jitter(width = 0.1, size = 1.5, alpha = 0.5) +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "none",
    panel.grid.major.x = element_blank()
  ) +
  labs(
    title = "Dispersion Check: Sub-sampled Data",
    x = "Group",
    y = "Distance to Centroid"
  )

# ==============================================================================
# 3. F1 vs. F2 inference on the subsample
# ==============================================================================
pkd_fit_f1_sub <- adonis2(dist_sub ~ Group, data = sub_metadata,
                          method = "euclidean", permutations = 99)

pkd_fit_f2_sub <- F2_adonis2(dist_sub ~ Group, data = sub_metadata,
                             method = "euclidean", permutations = 99)

f1_p_sub <- pkd_fit_f1_sub["Model", "Pr(>F)"]
f2_p_sub <- pkd_fit_f2_sub["Model", "Pr(>F2)"]

sub_results <- data.frame(
  Sample_Size = "10 Ctrl / 20 PKD",
  F1_Pval = f1_p_sub,
  F2_Pval = f2_p_sub,
  Verdict = ifelse(f2_p_sub > 0.05 & f1_p_sub < 0.05, "Conflicted", "Consistent")
)

sub_results %>%
  mutate(
    F1_Pval = sprintf("%.4f%s", F1_Pval, ifelse(F1_Pval < 0.05, "*", "")),
    F2_Pval = sprintf("%.4f", F2_Pval)
  ) %>%
  kable(
    col.names = c(
      "Subsample Configuration",
      "Standard F1 ($p$)",
      "F2 ($p$)",
      "Verdict"
    ),
    escape = FALSE,
    caption = "Check: Subsampled PKD Analysis"
  ) %>%
  kable_styling(full_width = FALSE)

print(pkd_fit_f1_sub)
print(pkd_fit_f2_sub)

# F1 permutation null distribution
psub_f1 <- permustats(pkd_fit_f1_sub)
f1_plot_data <- data.frame(F_vals = as.numeric(psub_f1$permutations))
f1_obs_stat <- pkd_fit_f1_sub["Model", "F"]

ggplot(f1_plot_data, aes(x = F_vals)) +
  geom_histogram(bins = 30, fill = "darkgoldenrod", color = "white", alpha = 0.7) +
  geom_vline(xintercept = f1_obs_stat, color = "red", linewidth = 1.2, linetype = "dashed") +
  annotate("text", x = f1_obs_stat, y = Inf, label = paste("Observed F =", round(f1_obs_stat, 3)),
           vjust = 2, hjust = -0.1, color = "red", fontface = "bold") +
  theme_minimal() +
  labs(
    title = "F1 Null Distribution: PKD Subsample",
    x = "Permuted Pseudo-F Statistic",
    y = "Frequency"
  )

# ==============================================================================
# 4. PCA ordination and component-level decomposition
# ==============================================================================
pca_fit <- prcomp(sub_metabolites, scale. = TRUE)
pca_scores <- pca_fit$x[, 1:5]
comp_data <- data.frame(pca_scores, Group = sub_metadata$Group)
var_exp <- (pca_fit$sdev^2 / sum(pca_fit$sdev^2)) * 100

# PCA scores plot (PC1 vs PC2)
ggplot(comp_data, aes(x = PC1, y = PC2, color = Group, fill = Group)) +
  geom_point(size = 3, alpha = 0.6) +
  stat_ellipse(geom = "polygon", alpha = 0.2) +
  theme_minimal() +
  labs(title = "PCA Scores: Subsampled Data",
       subtitle = "Visualizing the overlap of group centroids",
       x = paste0("PC1 (", round(var_exp[1], 1), "%)"),
       y = paste0("PC2 (", round(var_exp[2], 1), "%)"))

# Univariate t-tests on the top 5 PCs
univariate_tests <- lapply(1:5, function(i) {
  test <- t.test(comp_data[[i]] ~ comp_data$Group)
  data.frame(Component = paste0("PC", i), P_Value = test$p.value)
})

do.call(rbind, univariate_tests) %>%
  mutate(P_Value = sprintf("%.4f%s", P_Value, ifelse(P_Value < 0.05, "*", ""))) %>%
  kable(caption = "Decomposition: Univariate Significance of Components",
        booktabs = TRUE) %>%
  kable_styling(latex_options = "hold_position")

# Pairwise F2 sensitivity analysis across all 2D PC planes
pc_pairs <- combn(1:5, 2)

pairwise_results <- lapply(1:ncol(pc_pairs), function(i) {
  pair_idx <- pc_pairs[, i]
  dist_pair <- dist(pca_scores[, pair_idx], method = "euclidean")
  
  fit <- F2_adonis2(dist_pair ~ Group, data = sub_metadata,
                    permutations = 99, bootstrap = FALSE, plot.dist = FALSE)
  
  data.frame(
    Pair = paste0("PC", pair_idx[1], " vs PC", pair_idx[2]),
    F2_Stat = fit["Model", "F2"],
    P_Value = fit["Model", "Pr(>F2)"]
  )
})

do.call(rbind, pairwise_results) %>%
  mutate(P_Value = sprintf("%.4f%s", P_Value, ifelse(P_Value < 0.05, "*", ""))) %>%
  kable(caption = "Pairwise Sensitivity Analysis: F2 on 2D PC Planes",
        booktabs = TRUE) %>%
  kable_styling(latex_options = "hold_position")

## ----pkd-subsample-full, eval = FALSE-----------------------------------------
# pkd_fit_f1_sub <- adonis2(dist_sub ~ Group, data = sub_metadata,
#                           method = "euclidean", permutations = 999)
# 
# pkd_fit_f2_sub <- F2_adonis2(dist_sub ~ Group, data = sub_metadata,
#                              method = "euclidean", permutations = 999)
# 
# pairwise_results_full <- lapply(1:ncol(pc_pairs), function(i) {
#   pair_idx <- pc_pairs[, i]
#   dist_pair <- dist(pca_scores[, pair_idx], method = "euclidean")
# 
#   fit <- F2_adonis2(dist_pair ~ Group, data = sub_metadata,
#                     permutations = 999, bootstrap = FALSE, plot.dist = FALSE)
# 
#   data.frame(
#     Pair = paste0("PC", pair_idx[1], " vs PC", pair_idx[2]),
#     F2_Stat = fit["Model", "F2"],
#     P_Value = fit["Model", "Pr(>F2)"]
#   )
# })
# 
# do.call(rbind, pairwise_results_full)