PKD_Subsample_Study

library(dplyr)
library(vegan)
library(F2Adonis2)
library(kableExtra)
#> 
#> Attaching package: 'kableExtra'
#> The following object is masked from 'package:dplyr':
#> 
#>     group_rows
library(ggplot2)
library(combinat)
#> 
#> Attaching package: 'combinat'
#> The following object is masked from 'package:utils':
#> 
#>     combn

# 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)
#> 
#> Permutation test for homogeneity of multivariate dispersions
#> Permutation: free
#> Number of permutations: 99
#> 
#> Response: Distances
#>           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)   
#> Groups     1  321.27  321.27 6.8192     99   0.01 **
#> Residuals 28 1319.17   47.11                        
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 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)
#> Permutation test for adonis under reduced model
#> Permutation: free
#> Number of permutations: 99
#> 
#> adonis2(formula = dist_sub ~ Group, data = sub_metadata, permutations = 99, method = "euclidean")
#>          Df SumOfSqs      R2      F Pr(>F)  
#> Model     1     2406 0.05511 1.6332   0.05 *
#> Residual 28    41250 0.94489                
#> Total    29    43656 1.00000                
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(pkd_fit_f2_sub)
#> Permutation test for F2-based PERMANOVA
#>          Df SumOfSqs     F2 Pr(>F2)
#> Model     1     2406 1.4102     0.1
#> Residual 28    41250               
#> Total    29    43656

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