多组临床试验，样本量到底怎么算？

在进行科学研究，尤其是临床试验时，样本量的估算无疑是研究设计阶段最关键、也最令人头疼的环节之一。一个合适的样本量是确保研究结果可靠、结论具有说服力的基石。

我们大多数人都熟悉两组研究（例如，一个试验组 vs. 一个安慰剂组）的样本量计算。但如果你的研究设计更为复杂，包含了三个甚至更多的组别，比如需要比较不同剂量的药物疗效时，情况就变得不那么直观了。

这篇文章将探讨多组研究设计中样本量计算的两种核心思路，根据不同的研究目的，选择最合适的方法。

两组研究的样本量计算

在深入多组问题之前，先快速回顾一下最常见的两组比较。无论是比较两组的均数（如血压、肿瘤大小）还是率（如有效率、死亡率），样本量计算都依赖于以下几个核心参数：

• α (Alpha) 水平：即第一类错误的概率，通常设为 0.05，这是我们愿意承担的“假阳性”风险。
• 检验效能 (Power, 1-β)：即第二类错误的概率 β 的补集，代表了研究能够成功发现真实差异的能力，通常设为 0.80 或 0.90。
• 效应量 (Effect Size)：这是最关键的参数，代表了我们期望检测到的组间差异大小。例如，两组均数的差值、两组率的差值等。效应量需要基于现有文献或预实验来合理估计。
• 数据变异性 (Variability)：对于连续性变量，这通常指数据的标准差（Standard Deviation）。

对于两组设计，有大量成熟的公式和软件可以帮助我们轻松完成计算。

多组研究的样本量计算

假设一个研究需要评估一种新药的疗效，设计了四个组：

1. 安慰剂组 (Control Group)
2. 低剂量组 (Low-dose Group)
3. 中剂量组 (Medium-dose Group)
4. 高剂量组 (High-dose Group)

此时，我们的研究问题不再是简单的 “A 组 vs. B 组”。它可能包含两个层面的问题：

• 总体问题：这四个组的疗效在总体上有没有差异？
• 具体问题：低、中、高剂量组是否分别比安慰剂组更有效？

这两个不同的问题层面，直接导向了两种不同的样本量计算方法。

方法一，整体法 (Overall Approach)

这种方法主要回答上述的总体问题。

这种方法的目标是确保研究有足够的效能来检测出所有组别之间是否存在任何统计学差异。它对应的是多组数据分析时最先做的那个检验 —— 方差分析 (ANOVA) 的 F 检验（用于均数比较）或 卡方检验 (Chi-square Test)（用于率的比较）。

适用场景

• 动物实验：在动物实验中，研究者往往首先想证明某个干预措施在不同组别中“确实产生了效果”，而不急于深究具体是哪两组的差异。
• 探索性研究：在研究初期，当主要目的是筛选有潜力的剂量或干预方案时，这种方法可以提供一个初步的样本量估计。

局限性

这种方法虽然计算相对简单，但它的主要缺点是不够聚焦。即使方差分析的结果是显著的（p < 0.05），它也只能告诉你“至少有两组之间存在差异”，但无法保证你最关心的那个特定比较（例如，低剂量组 vs. 安慰剂组）也一定是显著的。这可能导致研究虽然得出了一个“有效”的笼统结论，却无法回答最核心的科学问题。

方法二，两两比较法 (Pairwise Comparison Approach)

这种方法是目前高质量临床试验中的“金标准”，它聚焦于回答上述的具体问题。

这种方法认为，一个多组研究的最终目的，通常是为了验证一个或多个特定的、有临床意义的组间差异。因此，样本量的计算也应该围绕这些核心比较来进行。

计算逻辑

确定核心比较

首先，明确研究中最关键的几个两两比较。在我们的例子中，最重要的比较通常是：

• 低剂量组 vs. 安慰剂组
• 中剂量组 vs. 安慰剂组
• 高剂量组 vs. 安慰剂组

多重比较校正 (调整 α 水平)

由于我们进行了多次比较，为了控制总体第一类错误率（即“假阳性”累积的风险），需要对 α 水平进行校正。最简单常用的方法是 Bonferroni 校正，即将原始 α 水平（如0.05）除以比较的次数（k）。

例如，在上面的例子中，有 3 个核心比较，所以校正后的 α’ = 0.05 / 3 ≈ 0.0167。

为每个核心比较计算样本量

现在，我们将这个多组问题拆解为多个两组问题。使用校正后的 α’ 水平，为每一个核心比较单独计算所需的样本量。

• 计算“低剂量 vs. 安慰剂”所需的样本量 N₁。
• 计算“中剂量 vs. 安慰剂”所需的样本量 N₂。
• 计算“高剂量 vs. 安慰剂”所需的样本量 N₃。

选取最大值

在计算出的 N₁, N₂, N₃ 中，选择最大的那个值，作为整个研究中每一组最终需要的样本量。

为什么选择最大值？

这个逻辑非常关键。选择最大值，意味着我们的研究对那个“效应量最小、最难检测出差异”的核心比较，也保证了足够的检验效能。

如果连最难的比较都有足够把握，那么那些效应量更大、更容易检测的比较，其效能自然也得到了保证。

适用场景

• 确证性临床试验，尤其是向药品监管机构（如 FDA, NMPA）提交注册申请或计划发表在高影响力期刊上的研究，这种方法是必须的。

• 目标明确的研究，当研究设计之初就明确了希望验证哪几个处理优于对照时。

我该如何选择？

那么，面对这两种方法，我们应该如何选择呢？答案取决于你的研究目的和研究阶段。

特性	整体法 Overall Approach	两两比较法 Pairwise Comparison Approach
研究问题	各组间总体有无差异？	特定组别间有无差异？
统计基础	方差分析 (ANOVA) / 卡方检验	两样本 t 检验 / 率的比较
优点	计算相对简单，所需样本量通常较小	结论更具体、严谨，统计效能有保障
缺点	结论笼统，可能无法回答核心问题	所需样本量通常更大，计算更繁琐
推荐场景	动物实验、探索性研究、初步筛选	人群临床试验、确证性研究、高质量论文

总而言之，对于严谨的人群临床试验，强烈推荐使用“两两比较法”。 它虽然要求更多的样本量，但确保了研究结果的可靠性和临床价值，让你在投入大量资源后，能够得到明确且有意义的结论。

在撰写研究计划书或论文时，务必清晰地说明你采用了哪种方法，并详细列出所有计算参数（α, Power, 效应量及其依据等），这是体现研究设计严谨性的重要一环。

代码示例

两两比较法

核心概念的转变 —— 从“检验效能”到“事件每变量 (EPV)”

对于预测模型，尤其是逻辑回归或其多分类扩展模型，最常用、最经典的样本量估算经验法则是**“事件每变量”原则 (Events Per Variable, EPV)**。

• 事件 (Events)：指的是模型要预测的结果中，发生数量最少的那个类别的样本数。在您的研究中，三个组别是“健康”、“前列腺增生 (BPH)”和“前列腺癌 (PCa)”。通常，在就诊人群中，“前列腺癌”患者的数量是最少的。因此，“事件数”就是指您样本中前列腺癌的病例数。
• 变量 (Variables)：指的是您计划纳入模型中进行筛选的所有候选预测变量的数量。这包括人口统计学特征（如年龄）、临床指标（如 PSA 总值、游离PSA 比率）、影像学特征、基因标记物等。关键点：这是一开始投入分析的变量总数，而不是最终留在模型里的变量数。

EPV 规则是什么？

这是一个经验法则，用来确保模型不会因为变量太多而样本太少导致过拟合。

• 传统标准 (EPV ≥ 10)：这是最广为人知的标准。它要求您的数据中，最少类别的样本数（事件数）至少是候选预测变量数量的 10 倍。
• 现代推荐 (EPV ≥ 20)：随着研究的深入，许多统计学家发现 EPV=10 仍然可能导致模型不稳定或系数估计不准。因此，为了构建更稳健、更可靠的模型，现在普遍推荐 EPV ≥ 20。

第一步：确定候选预测变量的数量 (p)

这是研究设计的第一步。您需要和您的团队一起，基于文献和临床经验，列出所有可能与前列腺癌诊断相关的候选变量。

• 示例
  • 年龄 (Age)
  • 总PSA (tPSA)
  • 游离/总PSA比率 (f/t PSA)
  • 前列腺体积 (Prostate Volume)
  • 直肠指检结果 (DRE result)
  • 经直肠超声特征 (TRUS features)
  • 某个新的血液生物标志物 A
  • 某个新的尿液生物标志物 B
  • …等等

假设您经过筛选，最终确定了 15个 候选预测变量。所以，p = 15。

第二步：确定最少事件组并计算其所需样本量

1. 确定最少事件组：在您的三组中，几乎可以肯定前列腺癌 (PCa) 组是样本量最少的。

2. 选择 EPV 规则：我们采用更严格的 EPV ≥ 20 标准来保证模型质量。

3. 计算最少事件组的样本量 (N_event)：
   N_event = p * EPV
N_pc = 15 (变量数) * 20 (EPV规则) = 300

   这意味着，为了构建一个稳健的模型，您的研究中至少需要 300 例前列腺癌患者。

第三步：根据各组比例，估算总样本量

现在知道了最核心的“瓶颈” —— 需要 300 例 PCa 患者。总样本量则取决于研究人群中“健康”、“BPH”和“PCa”的自然比例。这个比例需要通过回顾性数据或文献来估计。

• 场景假设：假设在您的多中心就诊人群中，经过初步筛查，这三类人群的大致比例是：

  • 健康人群: 40%
  • 前列腺增生 (BPH): 50%
  • 前列腺癌 (PCa): 10%

  这个比例说明，PCa确实是最少的组（10%）。

• 计算总样本量 (N_total)：
  既然 300 例 PCa 患者占总样本的 10%，那么总样本量可以推算出来：
  N_total = N_pc / 比例_pc
N_total = 300 / 0.10 = 3000

• 计算各组所需样本量：

  • 健康组: N_total * 比例_健康 = 3000 * 0.40 = 1200 例
  • BPH 组: N_total * 比例_BPH = 3000 * 0.50 = 1500 例
  • PCa 组: N_total * 比例_PCa = 3000 * 0.10 = 300 例

  结论：根据 EPV ≥ 20 的原则和预估的人群比例，大约需要总计 3000 名受试者，其中包括 1200 名健康人，1500 名 BPH 患者和 300 名 PCa 患者。

代码实现

import math

def calculate_prediction_model_samplesize(
    num_candidate_predictors: int,
    group_prevalences: dict,
    epv_rule: int = 20
) -> dict:
    """
    根据“事件每变量”(EPV)规则，为多分类预测模型估算样本量。

    Args:
        num_candidate_predictors (int): 计划纳入模型筛选的候选预测变量总数。
        group_prevalences (dict): 一个字典，包含各组的名称及其预估的患病率/比例。
                                  例如: {'Healthy': 0.4, 'BPH': 0.5, 'PCa': 0.1}
                                  所有比例之和应约等于1。
        epv_rule (int, optional): 每个变量需要的最少事件数。默认为20，这是推荐的稳健标准。

    Returns:
        dict: 一个包含详细计算结果的字典。
    """
    # 检查比例总和是否接近1
    if not math.isclose(sum(group_prevalences.values()), 1.0):
        print(f"警告: 各组比例之和为 {sum(group_prevalences.values()):.3f}，不等于1。请检查您的比例设置。")

    # 找到最少事件组（比例最低的组）
    min_prevalence_group = min(group_prevalences, key=group_prevalences.get)
    min_prevalence = group_prevalences[min_prevalence_group]

    if min_prevalence <= 0:
        raise ValueError("所有组的比例必须大于0。")

    # --- 核心计算步骤 ---
    # 1. 计算最少事件组所需的样本量
    required_samples_in_min_group = num_candidate_predictors * epv_rule

    # 2. 根据最少事件组的比例，反推出总样本量
    total_sample_size = required_samples_in_min_group / min_prevalence

    # 3. 计算每个组具体需要的样本量
    samples_per_group = {
        group: total_sample_size * prevalence
        for group, prevalence in group_prevalences.items()
    }
    
    # 为了实际招募，将结果取整
    required_samples_in_min_group_ceil = math.ceil(required_samples_in_min_group)
    total_sample_size_ceil = math.ceil(total_sample_size)
    samples_per_group_ceil = {
        group: math.ceil(size) for group, size in samples_per_group.items()
    }

    # --- 准备返回结果 ---
    results = {
        "parameters": {
            "num_candidate_predictors": num_candidate_predictors,
            "epv_rule": epv_rule,
            "group_prevalences": group_prevalences,
        },
        "calculation_summary": {
            "least_frequent_group": min_prevalence_group,
            "least_frequent_group_prevalence": min_prevalence,
            "required_samples_in_min_group": required_samples_in_min_group_ceil,
            "estimated_total_sample_size": total_sample_size_ceil,
        },
        "required_samples_per_group": samples_per_group_ceil
    }

    return results

def print_results(results: dict):
    """格式化并打印样本量计算结果。"""
    params = results['parameters']
    summary = results['calculation_summary']
    per_group = results['required_samples_per_group']

    print("--- 预测模型样本量估算结果 ---")
    print("\n[输入参数]")
    print(f"  - 候选预测变量数量: {params['num_candidate_predictors']}")
    print(f"  - EPV 经验法则: {params['epv_rule']}")
    print(f"  - 各组预估比例: {params['group_prevalences']}")

    print("\n[计算摘要]")
    print(f"  - 研究的瓶颈（最少事件组）: '{summary['least_frequent_group']}' (比例: {summary['least_frequent_group_prevalence']:.2%})")
    print(f"  - '{summary['least_frequent_group']}' 组至少需要样本数: {summary['required_samples_in_min_group']}")
    print(f"  - 由此估算的总样本量: {summary['estimated_total_sample_size']}")

    print("\n[各组所需样本量（向上取整）]")
    for group, count in per_group.items():
        print(f"  - {group}: {count} 例")
    print("\n--- 估算结束 ---")


# --- 主程序入口：在这里修改您的参数 ---
if __name__ == "__main__":
    
    # 1. 设置您研究的候选预测变量总数
    my_num_predictors = 15

    # 2. 设置您研究人群中各组的预估比例
    my_group_prevalences = {
        'Healthy': 0.40,
        'BPH': 0.50,
        'PCa': 0.10
    }
    
    # 3. (可选) 修改 EPV 规则，推荐使用 20
    my_epv_rule = 20

    # 执行计算
    try:
        final_results = calculate_prediction_model_samplesize(
            num_candidate_predictors=my_num_predictors,
            group_prevalences=my_group_prevalences,
            epv_rule=my_epv_rule
        )
        # 打印结果
        print_results(final_results)
    except ValueError as e:
        print(f"计算出错: {e}")

--- 预测模型样本量估算结果 ---

[输入参数]
  - 候选预测变量数量: 15
  - EPV 经验法则: 20
  - 各组预估比例: {'Healthy': 0.4, 'BPH': 0.5, 'PCa': 0.1}

[计算摘要]
  - 研究的瓶颈（最少事件组）: 'PCa' (比例: 10.00%)
  - 'PCa' 组至少需要样本数: 300
  - 由此估算的总样本量: 3000

[各组所需样本量（向上取整）]
  - Healthy: 1200 例
  - BPH: 1500 例
  - PCa: 300 例

--- 估算结束 ---

已知三分类样本数计算 (大模型)

基于 Riley 2020 框架

import numpy as np
from typing import List, Union, Dict

def calculate_expected_shrinkage_riley2020(
    n: int,
    p: int,
    assumed_r_squared_cs: Union[List[float], float]
) -> Dict[str, str]:
    """
    基于 Riley et al. (BMJ, 2020) 的公式，计算二元模型的预期收缩因子 (Expected shrinkage factor)。

    参数:
        n (int): 总样本量。
        p (int): 模型中候选预测变量的数量（包含截距除外）。
        assumed_r_squared_cs (float | list): 假设的 Cox–Snell R² 值。

    返回:
        dict: 各个 R² 假设下的收缩因子结果。
    """
    if n <= 0 or p <= 0:
        raise ValueError("n 和 p 必须为正整数。")

    if isinstance(assumed_r_squared_cs, float):
        assumed_r_squared_cs = [assumed_r_squared_cs]

    results = {}
    for r2 in assumed_r_squared_cs:
        if not (0 < r2 < 1):
            results[f"R²_cs = {r2}"] = "无效 (R² 必须在 0 与 1 之间)"
            continue

        # 计算模型的理论卡方值近似
        model_chi_sq_approx = n * (r2 / (1 - r2))

        if model_chi_sq_approx <= p:
            shrinkage = 0.0
        else:
            shrinkage = 1 - (p / model_chi_sq_approx)

        results[f"R²_cs = {r2}"] = f"{shrinkage:.3f}"

    return results


def evaluate_multinomial_sample_size_riley2020(
    group_counts: Dict[str, int],
    reference_category: str,
    p: int,
    assumed_r_squared_cs: Union[List[float], float]
):
    """
    依据 Riley 2020 样本量–收缩关系评估多分类逻辑回归开发的样本量充分性。
    每个类别与参考类别成对比较，并取最小收缩因子作为瓶颈。

    参数:
        group_counts (dict): 各类别样本量，例如 {"Healthy": 250, "BPH": 220, "Cancer": 180}
        reference_category (str): 参考类别名称。
        p (int): 候选预测变量数量。
        assumed_r_squared_cs (float | list): 假设的 Cox–Snell R²。
    """
    if reference_category not in group_counts:
        raise ValueError(f"参考类别 '{reference_category}' 不在 group_counts 中。")

    if isinstance(assumed_r_squared_cs, float):
        assumed_r_squared_cs = [assumed_r_squared_cs]

    print("=== 多分类模型样本量充分性评估 (Riley 2020 框架) ===")
    print(f"模型参数数量 p = {p}")
    print(f"参考类别：'{reference_category}' (n={group_counts[reference_category]})")
    print("-" * 55)

    n_ref = group_counts[reference_category]
    overall_min_shrinkage = {r2: 1.0 for r2 in assumed_r_squared_cs}
    bottleneck_info = {}

    for category_name, n_event in group_counts.items():
        if category_name == reference_category:
            continue

        n_total = n_ref + n_event

        print(f"\n对比: '{category_name}' (n={n_event}) vs. '{reference_category}' (n={n_ref})")
        evaluation = calculate_expected_shrinkage_riley2020(
            n=n_total,
            p=p,
            assumed_r_squared_cs=assumed_r_squared_cs
        )

        for scenario, s_str in evaluation.items():
            print(f"  - {scenario}: S = {s_str}")

            try:
                s_val = float(s_str)
                r2_val = float(scenario.split('=')[-1].strip())
                if s_val < overall_min_shrinkage[r2_val]:
                    overall_min_shrinkage[r2_val] = s_val
                    bottleneck_info[r2_val] = f"{category_name} vs {reference_category}"
            except ValueError:
                continue

    print("\n" + "-" * 55)
    print("=== 整体模型收缩因子评估结论 ===")

    for r2_val, min_s in overall_min_shrinkage.items():
        judgment = (
            "良好 (S ≥ 0.9, 过拟合风险低)" if min_s >= 0.9
            else "可接受 (0.85 ≤ S < 0.9, 有一定风险)" if min_s >= 0.85
            else "有风险 (S < 0.85, 明显过拟合风险)"
        )
        bottle = bottleneck_info.get(r2_val, "N/A")
        print(f"R²_cs = {r2_val:.2f} -> 最低收缩因子 S ≈ {min_s:.3f} [瓶颈: {bottle}]")
        print(f"  => 结论: {judgment}")


# === 实际示例 ===

your_groups = {
    "健康人群": 253,
    "前列腺增生": 215,
    "前列腺癌": 181
}

p_assumed = 12
r2_scenarios = [0.15, 0.25, 0.35]

evaluate_multinomial_sample_size_riley2020(
    group_counts=your_groups,
    reference_category="健康人群",
    p=p_assumed,
    assumed_r_squared_cs=r2_scenarios
)

已知三分类样本数计算 (大模型, 两阶段模型)

import numpy as np
from typing import List, Union, Dict, Tuple, Any, Optional


# =========================
# 基础工具函数
# =========================

def _ensure_list(x: Union[float, int, List[float], Tuple[float, ...]]) -> List[float]:
    """把标量或列表统一为列表（float）。"""
    if isinstance(x, (float, int)):
        return [float(x)]
    return [float(v) for v in x]


def _clamp01(x: float, eps: float = 1e-12) -> float:
    """把数值限制在 (0,1) 的开区间，避免 log 边界问题。"""
    return min(1 - eps, max(eps, x))


def r2_cs_max(prevalence: float) -> float:
    """
    计算在给定结局发生率（或阳性率）π 下，Cox–Snell R² 的最大值。
    R²_cs,max = 1 − exp{ 2 [ π ln π + (1 − π) ln (1 − π) ] }
    """
    pi = _clamp01(prevalence)
    return 1.0 - np.exp(2.0 * (pi * np.log(pi) + (1.0 - pi) * np.log(1.0 - pi)))


def nagelkerke_to_coxsnell(r2_nag: float, prevalence: float) -> float:
    """
    把 Nagelkerke R² 换算为 Cox–Snell R²：R²_cs = R²_nag × R²_cs,max(π)
    其中 π 为结局发生率（或阳性率）。
    """
    if not (0.0 < r2_nag < 1.0):
        raise ValueError("R²_nag 必须在 (0, 1) 之间")
    cs_max = r2_cs_max(prevalence)
    return float(r2_nag * cs_max)


def lr_chi_square_from_r2_cs(n: int, r2_cs: float) -> float:
    """
    由 Cox–Snell R² 与样本量 n 得到模型似然比 χ²：
    χ² = −n · ln(1 − R²_cs)
    """
    r2 = float(r2_cs)
    if not (0.0 < r2 < 1.0):
        raise ValueError("R²_cs 必须在 (0, 1) 之间")
    return -n * np.log1p(-r2)


def shrinkage_from_n_p_r2cs(n: int, p: int, r2_cs: float) -> float:
    """
    期望收缩因子（校准斜率的期望）：S = 1 − p/χ²，χ² 由 R²_cs 换算。
    """
    chi_sq = lr_chi_square_from_r2_cs(n, r2_cs)
    return max(0.0, 1.0 - (p / chi_sq)) if chi_sq > 0 else 0.0


def expected_delta_r2_cs(n: int, p: int, r2_cs_apparent: float) -> float:
    """
    期望的表观与“外部验证/乐观校正”后 Cox–Snell R² 的差值 ΔR²：
    利用 χ²_val ≈ χ²_app − p 的近似，有
    R²_app = 1 − exp(−χ²/n), R²_val = 1 − exp(−(χ² − p)/n)
    推得 ΔR² = (1 − R²_app) · (exp(p/n) − 1)
    """
    r2 = float(r2_cs_apparent)
    if not (0.0 < r2 < 1.0):
        raise ValueError("R²_cs 必须在 (0, 1) 之间")
    return (1.0 - r2) * (np.exp(p / n) - 1.0)


def halfwidth_prevalence_ci(n: int, prevalence: float, z: float = 1.96) -> float:
    """
    总体风险（结局/阳性率）π 的 95% 置信区间半宽度（正态近似）。
    halfwidth = z * sqrt(π(1−π)/n)
    """
    pi = _clamp01(prevalence)
    return float(z * np.sqrt(pi * (1.0 - pi) / n))


def required_n_for_target_shrinkage(p: int, r2_cs: float, S_target: float = 0.90) -> float:
    """
    给定 p、R²_cs 与目标收缩因子 S_target，反推所需总样本量 n。
    n ≥ [ p / (1 − S_target) ] / [ −ln(1 − R²_cs) ]
    """
    if not (0.0 < r2_cs < 1.0):
        raise ValueError("R²_cs 必须在 (0, 1) 之间")
    if not (0.0 < S_target < 1.0):
        raise ValueError("S_target 必须在 (0, 1) 之间")
    denom = -np.log1p(-r2_cs)  # > 0
    return (p / (1.0 - S_target)) / denom


def required_n_for_delta_r2(p: int, r2_cs: float, delta_target: float = 0.05) -> float:
    """
    给定 p、R²_cs 与目标 ΔR²（表观−外部/校正，建议 ≤0.05），反推所需 n。
    由 ΔR² = (1−R²) [exp(p/n) − 1] ≤ δ 推得：
    n ≥ p / ln( 1 + δ/(1−R²) )
    """
    if not (0.0 < r2_cs < 1.0):
        raise ValueError("R²_cs 必须在 (0, 1) 之间")
    if not (0.0 < delta_target < 1.0):
        raise ValueError("delta_target 必须在 (0, 1) 之间")
    denom = np.log(1.0 + (delta_target / (1.0 - r2_cs)))
    return p / denom


def required_n_for_prevalence_precision(prevalence: float, delta_target: float = 0.05, z: float = 1.96) -> float:
    """
    总体风险（结局/阳性率）π 的精度要求：95% CI 半宽度 ≤ δ。
    n ≥ z² · π(1−π) / δ²
    """
    pi = _clamp01(prevalence)
    return (z ** 2) * pi * (1.0 - pi) / (delta_target ** 2)


# =========================
# Riley 2020: 期望收缩因子计算（修正版）
# =========================

def calculate_expected_shrinkage_riley2020(
    n: int,
    p: int,
    assumed_r_squared_cs: Union[List[float], float]
) -> Dict[str, str]:
    """
    按 Riley 等的公式计算期望收缩因子 S（基于 Cox–Snell R²）。
    修正点：χ² = −n · ln(1 − R²_cs)
    """
    if n <= 0 or p <= 0:
        raise ValueError("n 和 p 必须为正整数。")
    r2_list = _ensure_list(assumed_r_squared_cs)

    results: Dict[str, str] = {}
    for r2 in r2_list:
        if not (0.0 < r2 < 1.0):
            results[f"R²_cs = {r2:.2f}"] = "无效 (应在 0~1 之间)"
            continue
        model_chi_sq = lr_chi_square_from_r2_cs(n, r2)
        shrinkage = max(0.0, 1.0 - (p / model_chi_sq)) if model_chi_sq > 0 else 0.0
        results[f"R²_cs = {r2:.2f}"] = f"{shrinkage:.3f}"
    return results


# =========================
# 两阶段模型评估（扩展版，含三条标准）
# =========================

def _classify_shrinkage(s: float) -> str:
    if s >= 0.90:
        return "良好 (S ≥ 0.90)"
    elif s >= 0.85:
        return "可接受 (0.85 ≤ S < 0.90)"
    else:
        return "有风险 (S < 0.85)"


def _normalize_assumed_r2_input(
    assumed_r2: Union[float, List[float], Dict[str, Union[float, List[float]]]]
) -> Tuple[List[float], List[float], bool]:
    """
    归一化 R² 输入。
    - 若传入标量或列表：两阶段共用同一组 R²，返回 pairwise=True
    - 若传入字典 {'stage1': ..., 'stage2': ...}：分别使用，返回 pairwise=False
    """
    if isinstance(assumed_r2, dict):
        r2_stage1 = _ensure_list(assumed_r2.get("stage1", []))
        r2_stage2 = _ensure_list(assumed_r2.get("stage2", []))
        return r2_stage1, r2_stage2, False
    else:
        r2_list = _ensure_list(assumed_r2)
        return r2_list, r2_list, True


def _convert_r2_to_cs_list(
    r2_list: List[float],
    r2_type: str,
    prevalence: float
) -> List[float]:
    """按类型把 R² 转换为 Cox–Snell 列表。"""
    if r2_type.lower() in ("cs", "coxsnell", "cox-snell"):
        return r2_list
    elif r2_type.lower() in ("nag", "nagelkerke"):
        return [nagelkerke_to_coxsnell(r2, prevalence) for r2 in r2_list]
    else:
        raise ValueError("r2_type 必须是 'cs' 或 'nagelkerke'。")


def evaluate_two_stage_models(
    n_healthy: int,
    n_bph: int,
    n_cancer: int,
    p_stage1: int,
    p_stage2: int,
    assumed_r2: Union[float, List[float], Dict[str, Union[float, List[float]]]],
    r2_type: str = "cs",
    target_shrinkage: float = 0.90,
    target_delta_r2: float = 0.05,
    target_delta_pi: float = 0.05,
    z_value: float = 1.96,
) -> Dict[str, Any]:
    """
    两阶段机器学习/回归预测模型样本量充分性评估（Riley 2020 框架）
    - 阶段1：健康 vs 异常（BPH + CA）
    - 阶段2：BPH vs CA（在异常者中）
    同时评估三条常用标准：
      1) 目标收缩因子 S ≥ target_shrinkage（默认 0.90）
      2) 表观与校正后 R² 的差 ΔR² ≤ target_delta_r2（默认 0.05）
      3) 总体风险（结局率）π 的 95% CI 半宽度 ≤ target_delta_pi（默认 0.05）

    参数
    - assumed_r2: 可为
        · 标量/列表：两阶段共用同一组 R²
        · 字典：{'stage1': 标量/列表, 'stage2': 标量/列表}
    - r2_type: 'cs'（Cox–Snell）或 'nagelkerke'
    """
    # 基本校验
    for x in [n_healthy, n_bph, n_cancer, p_stage1, p_stage2]:
        if not isinstance(x, (int, np.integer)) or x <= 0:
            raise ValueError("所有样本量与参数数必须为正整数。")
    if r2_type.lower() not in ("cs", "coxsnell", "cox-snell", "nag", "nagelkerke"):
        raise ValueError("r2_type 必须是 'cs' 或 'nagelkerke'。")

    # 规模与发生率
    n_stage1 = int(n_healthy + n_bph + n_cancer)
    n_stage2 = int(n_bph + n_cancer)
    n_abnormal = int(n_bph + n_cancer)

    pi_stage1 = n_abnormal / n_stage1  # “异常率”
    pi_stage2 = n_cancer / n_stage2    # “癌症率（在异常者中）”

    # R² 输入处理
    r2_s1_raw, r2_s2_raw, pairwise_same = _normalize_assumed_r2_input(assumed_r2)
    r2_s1_cs_list = _convert_r2_to_cs_list(r2_s1_raw, r2_type, pi_stage1)
    r2_s2_cs_list = _convert_r2_to_cs_list(r2_s2_raw, r2_type, pi_stage2)

    # 打印标题
    print("\n=== 两阶段机器学习预测模型样本量充分性评估 (Riley 2020 框架) ===\n")

    # 阶段1输出
    print("【第一阶段】健康 vs 异常 (BPH + CA)")
    print(f"  - 样本量: 健康={n_healthy}, 异常={n_abnormal}, 总={n_stage1}")
    print(f"  - 结局率(异常率) π₁ = {pi_stage1:.3f}")
    print(f"  - 模型参数数 p₁ = {p_stage1}")
    print(f"  - R² 类型: {'Cox–Snell' if r2_type.lower() in ('cs','coxsnell','cox-snell') else 'Nagelkerke (已换算为 Cox–Snell 用于计算)'}")

    stage1_items = []
    for r2 in r2_s1_cs_list:
        if not (0.0 < r2 < 1.0):
            print(f"  - R²_cs = {r2:.2f}: 无效 (应在 0~1 之间)")
            continue
        s = shrinkage_from_n_p_r2cs(n_stage1, p_stage1, r2)
        delta_r2 = expected_delta_r2_cs(n_stage1, p_stage1, r2)
        chi_sq = lr_chi_square_from_r2_cs(n_stage1, r2)
        p_max_for_targetS = chi_sq * (1.0 - target_shrinkage)  # 达到目标 S 时允许的最大自由度

        n_req_S = required_n_for_target_shrinkage(p_stage1, r2, target_shrinkage)
        n_req_dR2 = required_n_for_delta_r2(p_stage1, r2, target_delta_r2)

        print(f"  - R²_cs = {r2:.2f}: S = {s:.3f}  ->  {_classify_shrinkage(s)}")
        print(f"               ΔR²(表观-校正) ≈ {delta_r2:.3f}  ->  {'通过' if delta_r2 <= target_delta_r2 else '不通过'} (阈值 {target_delta_r2})")
        print(f"               建议: 若追求 S≥{target_shrinkage:.2f}，当前 R² 下 p≤{p_max_for_targetS:.1f}；或所需 n≈{np.ceil(n_req_S):.0f}")
        print(f"                     若追求 ΔR²≤{target_delta_r2:.2f}，所需 n≈{np.ceil(n_req_dR2):.0f}")

        stage1_items.append({
            "r2_cs": r2,
            "S": s,
            "delta_r2": delta_r2,
            "chi_sq": chi_sq,
            "p_max_for_targetS": p_max_for_targetS,
            "n_required_for_S": n_req_S,
            "n_required_for_deltaR2": n_req_dR2,
        })

    # 阶段1：总体风险精度
    hw1 = halfwidth_prevalence_ci(n_stage1, pi_stage1, z=z_value)
    n_req_pi1 = required_n_for_prevalence_precision(pi_stage1, target_delta_pi, z=z_value)
    print(f"  - 总体风险(异常率)精度: 95% CI 半宽度 ≈ {hw1:.3f}  ->  {'通过' if hw1 <= target_delta_pi else '不通过'} (阈值 {target_delta_pi})")
    print(f"    所需 n 以满足半宽度 ≤ {target_delta_pi:.2f}: 约 {np.ceil(n_req_pi1):.0f}")

    # 阶段2输出
    print("\n【第二阶段】前列腺增生 (BPH) vs 前列腺癌 (CA)")
    print(f"  - 样本量: BPH={n_bph}, 癌={n_cancer}, 总={n_stage2}")
    print(f"  - 结局率(癌症率) π₂ = {pi_stage2:.3f}")
    print(f"  - 模型参数数 p₂ = {p_stage2}")
    print(f"  - R² 类型: {'Cox–Snell' if r2_type.lower() in ('cs','coxsnell','cox-snell') else 'Nagelkerke (已换算为 Cox–Snell 用于计算)'}")

    stage2_items = []
    for r2 in r2_s2_cs_list:
        if not (0.0 < r2 < 1.0):
            print(f"  - R²_cs = {r2:.2f}: 无效 (应在 0~1 之间)")
            continue
        s = shrinkage_from_n_p_r2cs(n_stage2, p_stage2, r2)
        delta_r2 = expected_delta_r2_cs(n_stage2, p_stage2, r2)
        chi_sq = lr_chi_square_from_r2_cs(n_stage2, r2)
        p_max_for_targetS = chi_sq * (1.0 - target_shrinkage)

        n_req_S = required_n_for_target_shrinkage(p_stage2, r2, target_shrinkage)
        n_req_dR2 = required_n_for_delta_r2(p_stage2, r2, target_delta_r2)

        print(f"  - R²_cs = {r2:.2f}: S = {s:.3f}  ->  {_classify_shrinkage(s)}")
        print(f"               ΔR²(表观-校正) ≈ {delta_r2:.3f}  ->  {'通过' if delta_r2 <= target_delta_r2 else '不通过'} (阈值 {target_delta_r2})")
        print(f"               建议: 若追求 S≥{target_shrinkage:.2f}，当前 R² 下 p≤{p_max_for_targetS:.1f}；或所需 n≈{np.ceil(n_req_S):.0f}")
        print(f"                     若追求 ΔR²≤{target_delta_r2:.2f}，所需 n≈{np.ceil(n_req_dR2):.0f}")

        stage2_items.append({
            "r2_cs": r2,
            "S": s,
            "delta_r2": delta_r2,
            "chi_sq": chi_sq,
            "p_max_for_targetS": p_max_for_targetS,
            "n_required_for_S": n_req_S,
            "n_required_for_deltaR2": n_req_dR2,
        })

    # 阶段2：总体风险精度
    hw2 = halfwidth_prevalence_ci(n_stage2, pi_stage2, z=z_value)
    n_req_pi2 = required_n_for_prevalence_precision(pi_stage2, target_delta_pi, z=z_value)
    print(f"  - 总体风险(癌症率)精度: 95% CI 半宽度 ≈ {hw2:.3f}  ->  {'通过' if hw2 <= target_delta_pi else '不通过'} (阈值 {target_delta_pi})")
    print(f"    所需 n 以满足半宽度 ≤ {target_delta_pi:.2f}: 约 {np.ceil(n_req_pi2):.0f}")

    # 综合评估（仅当两阶段使用同一组 R² 时可逐一配对）
    print("\n=== 综合评估 ===")
    if pairwise_same and len(stage1_items) == len(stage2_items) and len(stage1_items) > 0:
        for i, (it1, it2) in enumerate(zip(stage1_items, stage2_items)):
            r2_label = f"R²_cs = {it1['r2_cs']:.2f}"
            s1, s2 = it1["S"], it2["S"]
            overall_min = min(s1, s2)
            worst_stage = "第一阶段" if s1 < s2 else "第二阶段"
            judgment = (
                "良好 (过拟合风险低)" if overall_min >= 0.90
                else "可接受 (轻度风险)" if overall_min >= 0.85
                else "有风险 (需增加样本或减少变量)"
            )
            print(f"{r2_label} -> 最低 S ≈ {overall_min:.3f} (瓶颈: {worst_stage})")
            print(f"  -> 结论: {judgment}")
    else:
        print("两阶段使用了不同的 R² 场景或数量，已分别汇报各自结果；不做配对综合判定。")

    print("\n=== 结束 ===\n")

    # 返回结构化结果，便于程序化使用
    return {
        "stage1": {
            "n": n_stage1,
            "p": p_stage1,
            "prevalence": pi_stage1,
            "items": stage1_items,
            "halfwidth_pi": hw1,
            "n_required_for_pi_precision": n_req_pi1,
        },
        "stage2": {
            "n": n_stage2,
            "p": p_stage2,
            "prevalence": pi_stage2,
            "items": stage2_items,
            "halfwidth_pi": hw2,
            "n_required_for_pi_precision": n_req_pi2,
        },
        "pairwise_same_r2": pairwise_same,
        "targets": {
            "shrinkage": target_shrinkage,
            "delta_r2": target_delta_r2,
            "delta_pi": target_delta_pi,
            "z_value": z_value,
            "r2_type_input": r2_type,
        }
    }


# =========================
# 示例运行
# =========================
# 1.
# Predictors of Bone Metastases at 68Ga-PSMA-11 PET/CT in Hormone-Sensitive Prostate Cancer (HSPC) Patients with Early Biochemical Recurrence or Persistence
# In the multivariate binary logistic regression model (accuracy: 71.2%, Nagelkerke-R2: 13%),
# 2.
# The Role of Ga-68 PSMA PET/CT Scan on Differentiating of Oligometastatic and High Risk Prostate Cancer
# The explanatory ratio of these factors to the multiple metastasis category was sufficiently high (Nagelkerke R2=0.459). The contribution of SUVmax to the model was positive [odds ratio (OR=1.42) and significant (p=0.038; Table 4).
# 3.
# Medical Help-Seeking for Sexual Concerns in Prostate Cancer Survivors
# had moderate to mild ED (B = 0.17, Exp(B) = 1.18, 95% CI = 1.13–1.24, P < .001) had higher rates of recent treatment use for ED, χ2(5) = 93.93, P < .001 (Nagelkerke R2 = 0.31).
# had higher rates of sexual help-seeking prior to the study, χ2(5) = 45.61, P < .001 (Nagelkerke R2 = 0.15).

if __name__ == "__main__":
    # 示例：与您先前给定的数字一致，且 R² 视为 Cox–Snell
    evaluate_two_stage_models(
        n_healthy=253,
        n_bph=215,
        n_cancer=181,
        p_stage1=8,
        p_stage2=8,
        assumed_r2=[0.15, 0.25, 0.35],  # 若为 Nagelkerke，请把 r2_type 设为 'nagelkerke'
        r2_type="cs",
        target_shrinkage=0.90,
        target_delta_r2=0.05,
        target_delta_pi=0.05,
        z_value=1.96,
    )