674 微生物组效应量与统计功效¶

# 微生物组效应量计算
library(vegan)        # 生态学统计
library(phyloseq)     # 微生物组数据
library(pwr)          # 统计功效分析

# ============ 1. PERMANOVA效应量（R²） ============
# R²是PERMANOVA最直接的效应量指标
# 准备数据
otu <- as(otu_table(ps), "matrix")       # OTU表
if (!taxa_are_rows(ps)) otu <- t(otu)    # 确保物种在行
otu_t <- t(otu)                          # 转置（样本在行）
meta <- as(sample_data(ps), "data.frame")  # 元数据

# 计算Bray-Curtis距离
bc_dist <- vegdist(otu_t, method = "bray")  # Bray-Curtis距离

# PERMANOVA检验
perm_result <- adonis2(
  bc_dist ~ Group + Age + BMI,           # 公式：距离 ~ 分组 + 协变量
  data = meta,                           # 元数据
  permutations = 999,                    # 置换次数
  by = "margin"                          # 边际检验（控制其他变量）
)
print(perm_result)  # 查看结果

# 提取R²（效应量）
r2_group <- perm_result$R2[1]            # Group的R²
cat("Group的效应量(R²):", r2_group, "\n")
# 微生物组中R²通常2-5%是正常的！
# R² > 10% 在微生物组研究中已经很大了

# ============ 2. ω²（Omega-squared，校正后的效应量） ============
# R²会高估效应量，ω²是校正后的无偏估计
omega_squared <- function(adonis_result, term_idx = 1) {
  ss_term <- adonis_result$SumOfSqs[term_idx]  # 组间平方和
  ss_total <- sum(adonis_result$SumOfSqs)       # 总平方和
  ms_residual <- adonis_result$SumOfSqs[nrow(adonis_result)] /
                 adonis_result$Df[nrow(adonis_result)]  # 残差均方
  n <- sum(adonis_result$Df) + 1                # 总样本数
  df_term <- adonis_result$Df[term_idx]         # 组间自由度

  omega2 <- (ss_term - df_term * ms_residual) /
            (ss_total + ms_residual)             # ω²公式
  return(omega2)
}

w2 <- omega_squared(perm_result)
cat("校正效应量(ω²):", w2, "\n")

# ============ 3. 单变量效应量（Cohen's d） ============
# 对单个物种计算效应量
cohens_d <- function(x1, x2) {
  n1 <- length(x1)                       # 组1样本数
  n2 <- length(x2)                       # 组2样本数
  s_pooled <- sqrt(((n1-1)*var(x1) + (n2-1)*var(x2)) /
                    (n1+n2-2))           # 合并标准差
  d <- (mean(x1) - mean(x2)) / s_pooled  # Cohen's d
  return(d)
}

# 对每个物种计算Cohen's d
group1 <- meta$Group == "Disease"        # 疾病组
group2 <- meta$Group == "Control"        # 对照组
effect_sizes <- apply(otu_t, 2, function(x) {
  cohens_d(x[group1], x[group2])         # 计算每个物种的d值
})
top_effects <- sort(abs(effect_sizes), decreasing = TRUE)[1:10]
cat("\n效应量最大的10个物种:\n")
print(round(top_effects, 3))

# 统计功效分析
library(pwr)          # 功效分析包
library(micropower)   # 微生物组专用功效分析（如可用）

# ============ 1. 基于Cohen's d的功效分析 ============
# 场景：两组比较（t检验），估算所需样本量

# 已知效应量，求样本量
power_result <- pwr.t.test(
  d = 0.5,            # 效应量（中等效应）
  sig.level = 0.05,   # 显著性水平
  power = 0.80,       # 期望功效（80%）
  type = "two.sample", # 双样本t检验
  alternative = "two.sided"  # 双侧检验
)
cat("中等效应量(d=0.5)需要每组", ceiling(power_result$n), "个样本\n")

# 不同效应量的样本量需求
effect_sizes <- c(0.2, 0.3, 0.5, 0.8, 1.0)  # 效应量范围
sample_sizes <- sapply(effect_sizes, function(d) {
  res <- pwr.t.test(d = d, sig.level = 0.05, power = 0.80,
                     type = "two.sample")
  ceiling(res$n)                         # 每组所需样本数
})
power_table <- data.frame(
  effect_size = effect_sizes,
  samples_per_group = sample_sizes,
  total_samples = sample_sizes * 2
)
print(power_table)
# 输出：
# effect_size samples_per_group total_samples
#         0.2               394           788  ← 小效应量需要很多样本
#         0.3               176           352
#         0.5                64           128  ← 中效应量
#         0.8                26            52
#         1.0                17            34  ← 大效应量

# ============ 2. PERMANOVA功效分析 ============
# 基于模拟的PERMANOVA功效估算
permanova_power <- function(n_per_group, n_groups = 2,
                             effect_size_r2 = 0.05,
                             n_vars = 100, n_sims = 100) {
  """通过模拟估算PERMANOVA的统计功效"""
  significant <- 0                       # 显著次数计数

  for (i in 1:n_sims) {
    n_total <- n_per_group * n_groups    # 总样本数
    groups <- rep(1:n_groups, each = n_per_group)  # 分组

    # 模拟数据：组间有一定差异
    data_matrix <- matrix(rnorm(n_total * n_vars), nrow = n_total)
    # 给不同组添加效应
    shift <- sqrt(effect_size_r2 / (1 - effect_size_r2)) * sqrt(n_vars)
    data_matrix[groups == 2, 1:5] <- data_matrix[groups == 2, 1:5] + shift

    # PERMANOVA检验
    dist_mat <- vegdist(abs(data_matrix), method = "bray")
    result <- adonis2(dist_mat ~ factor(groups), permutations = 199)
    if (result$`Pr(>F)`[1] < 0.05) {
      significant <- significant + 1     # 统计显著次数
    }
  }
  power <- significant / n_sims          # 功效=显著比例
  return(power)
}

# 不同样本量的功效
cat("\nPERMANOVA功效分析 (R²=0.05):\n")
for (n in c(10, 20, 30, 50, 100)) {
  pwr <- permanova_power(n_per_group = n, effect_size_r2 = 0.05, n_sims = 50)
  cat(sprintf("  每组%d个样本: 功效=%.2f\n", n, pwr))
}

# ============ 3. 功效曲线可视化 ============
library(ggplot2)      # 绑图

# 生成功效曲线数据
power_data <- expand.grid(
  n = seq(10, 200, by = 10),             # 样本量范围
  d = c(0.2, 0.3, 0.5, 0.8)             # 效应量
)
power_data$power <- mapply(function(n, d) {
  pwr.t.test(n = n, d = d, sig.level = 0.05,
             type = "two.sample")$power
}, power_data$n, power_data$d)

# 绘制功效曲线
ggplot(power_data, aes(x = n, y = power, color = factor(d))) +
  geom_line(linewidth = 1.2) +           # 线图
  geom_hline(yintercept = 0.80, linetype = "dashed", color = "red") +
  annotate("text", x = 180, y = 0.82, label = "80% 功效线",
           color = "red") +               # 标注80%线
  labs(x = "每组样本量", y = "统计功效",
       color = "效应量 (Cohen's d)",
       title = "统计功效 vs 样本量") +
  theme_minimal() +                       # 简洁主题
  scale_y_continuous(limits = c(0, 1))    # y轴范围
ggsave("power_curve.png", width = 8, height = 5, dpi = 150)

# 微生物组效应量的特殊性分析
import numpy as np   # 数值计算
import pandas as pd  # 数据处理

# 1. 微生物组研究中的典型效应量
typical_effects = pd.DataFrame({
    "研究类型": [
        "疾病 vs 健康（IBD）",
        "疾病 vs 健康（T2D）",
        "饮食干预",
        "抗生素干预",
        "地理/种族差异",
        "季节变化",
        "个体内时间变异"
    ],
    "典型R²": [0.05, 0.02, 0.03, 0.15, 0.10, 0.01, 0.30],
    "效应大小": ["中", "小", "小-中", "大", "中-大", "小", "大"],
    "建议每组样本量": ["50-100", "100-200", "80-150", "20-40",
                       "30-60", "200+", "20-30"]
})
print("微生物组研究典型效应量:")
print(typical_effects.to_string(index=False))

# 2. 多重比较校正对功效的影响
def adjusted_power(n, d, n_tests, alpha=0.05):
    """考虑多重比较校正后的功效"""
    from scipy.stats import norm  # 正态分布
    # Bonferroni校正后的alpha
    alpha_adj = alpha / n_tests          # 校正后的显著性水平
    # 非中心参数
    ncp = d * np.sqrt(n / 2)             # 非中心参数
    # 临界值
    z_crit = norm.ppf(1 - alpha_adj / 2)  # 临界z值
    # 功效
    power = 1 - norm.cdf(z_crit - ncp) + norm.cdf(-z_crit - ncp)
    return power

# 不同物种数的功效影响
print("\n多重比较对功效的影响（n=50/组, d=0.5）:")
for n_taxa in [10, 50, 100, 500, 1000]:
    pwr = adjusted_power(n=50, d=0.5, n_tests=n_taxa)
    print(f"  {n_taxa}个物种: 校正后功效={pwr:.3f}")

# 3. 效应量的置信区间
def cohens_d_ci(d, n1, n2, alpha=0.05):
    """Cohen's d的置信区间（近似）"""
    from scipy.stats import norm
    se = np.sqrt((n1+n2)/(n1*n2) + d**2/(2*(n1+n2)))  # 标准误
    z = norm.ppf(1 - alpha/2)            # z值
    ci_lower = d - z * se                # 下界
    ci_upper = d + z * se                # 上界
    return ci_lower, ci_upper

# 示例
d = 0.5                                  # 点估计
ci = cohens_d_ci(d, n1=30, n2=30)        # 30/组
print(f"\nCohen's d = {d}, 95% CI: [{ci[0]:.3f}, {ci[1]:.3f}]")
print("→ 效应量的不确定性很大，需要足够的样本量")

# 效应量解读
Cohen's d: 0.2(小) 0.5(中) 0.8(大)
R²: 0.01(小) 0.06(中) 0.14(大)
微生物组R²: 0.02-0.05 是很常见的

# 样本量经验法则
小效应(d=0.2): 每组~400个样本
中效应(d=0.5): 每组~65个样本
大效应(d=0.8): 每组~26个样本

# 微生物组样本量建议
横断面研究: 50-200/组（取决于效应量）
纵向研究: 20-50/组（重复测量增加功效）
干预研究: 30-80/组

# 功效分析步骤
1. 估算效应量（文献/预实验）
2. 设定α=0.05, power=0.80
3. 用pwr包或模拟计算样本量
4. 考虑多重比较的功效损失
5. 加20%余量应对样本丢失

# 报告效应量的最佳实践
- 同时报告p值和效应量
- 给出效应量的置信区间
- 对比文献中同类研究的效应量