Zenn
💊

江島法に基づくサンプルサイズの厳密解をPythonで導出【BE試験】

に公開
Python
データ分析
統計学
臨床試験
サンプルサイズ設計
tech

江島らのBE試験サンプルサイズ厳密解をPythonで再現してみた

生物学的同等性試験（Bioequivalence試験、以下BE試験）では、被験薬と基準薬が同等であることを統計的に検証する必要があります。その際、サンプルサイズの設計は非常に重要な工程の一つです。

多くの場合、近似的な計算式（正規分布近似など）で済まされることが多いのですが、江島らの論文では 非心t分布を用いて検出力から“厳密に”必要なサンプルサイズを求める方法 が紹介されています。

本記事では、その方法を Python で再現 し、近似ではなく「数理的に正しい最小n」を導出するプロセスを紹介します。

方法の概要

サンプルサイズnを求めるにあたって、次のようなパラメータを使います：

  • α：有意水準（今回は 0.10）
  • β：検出力の補数（1 - β = 0.80）
  • MSW：SASで算出された誤差分散（MS(ERROR)）
  • 群数：2群（被験薬群、基準薬群）

このとき、次の流れでnを導出します：

  1. 初期nをMSWから仮決め
  2. 自由度dofを計算
  3. t分布の臨界値（両側）を求める
  4. 検出力が条件を満たす非心度δを逆算
  5. 効果量λとの大小比較で収束判定（λ ≥ δ）

このロジックを、以下のPythonコードで実装します。

実装コード

# coding: UTF-8
""" 江島ら論文によるサンプルサイズ算出(BE) """
""" インポート """
import math
import numpy as np
from scipy.stats import t, nct
from scipy.optimize import brentq
from loguru import logger
# ===========================================================================================
""" 実行関数 """
# 暫定的なnを求める（収束がなかなかしない場合は適宜調整してください）
def numbers(msw: float) -> int:
    return max(3, int(msw * 100))


# 自由度
def degree_of_freedom(n: int, group_num: int) -> int:
    return (n - 1) * group_num


# 有効サンプルサイズ
def effect_size(n: int, group_num: int) -> float:
    n_total_1 = 1
    n_total_2 = 0
    for _ in range(group_num):
        n_total_1 = n_total_1 * n
        n_total_2 += n
    return n_total_1 / n_total_2


# 両側検定臨界値(t)
def t_critical(alpha: float, dof: int) -> float:
    return t.ppf(1 - alpha / 2, dof)


# 検出力誤差関数
def power_diff(delta: float):
    # 両側： P(|t| > t_crit) = P(t > t_crit) + P(t < -t_crit)
    power = 1 - nct.cdf(t_crit, dof, delta) + nct.cdf(-t_crit, dof, delta)
    return power - power_target


# λcalc
def lambda_calc(n: int, msw: float) -> float:
    """
    効果量=math.log10(1.2)
    これはBE試験における標準的な等価範囲（±20%）であることを示す
    """
    return np.sqrt(n) * math.log10(1.2) / np.sqrt(msw)

# ===========================================================================================
""" Main """
if __name__ == "__main__":
    # パラメータ設定
    alpha: float = 0.10
    beta: float = 0.2
    power_target: float = 1 - beta
    group_num: int = 2
    msw: float = 0.0125

    # 実行処理
    # 基準サンプルサイズ
    n = numbers(msw=msw)
    # 表のδを超えるまで繰り返す
    while True:
        # 分布表の値
        # 自由度
        dof = degree_of_freedom(n=n, group_num=group_num)
        # 有効サンプルサイズ
        n_eff = effect_size(n=n, group_num=group_num)
        # 両側検定臨界値(t)
        t_crit = t_critical(alpha=alpha, dof=dof)
        # δを求める（非心度）
        delta = brentq(power_diff, 0, 30)
        # 効果量 d = δ / √n_eff
        d = delta / np.sqrt(n_eff)

        # 実際の値（λcalc）
        lam_c = lambda_calc(n=n, msw=msw)
        # 収束判定
        if lam_c >= delta:
            logger.debug("収束しました。")
            break
        else:
            n += 1

    # デバッグ
    logger.debug(f"=============================================")
    logger.debug(f"n={n}")
    logger.debug(f"自由度={dof}")
    logger.debug(f"effect size={n_eff}")
    logger.debug(f"t critical={t_crit: 4f}")
    logger.debug(f"分布表に値するδ={delta: 4f}、効果量d={d: 4f}")
    logger.debug(f"λcalc={lam_c: 4f}、n={n}")

出力ログ

2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:81 - 収束しました。
2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:87 - =============================================
2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:88 - n=14
2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:89 - 自由度=26
2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:90 - effect size=7.0
2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:91 - t critical= 1.705618
2025-07-02 21:56:14.259 | DEBUG    | __main__:<module>:92 - 分布表に値するδ= 2.553768、効果量d= 0.965234

🔍 解説

  • n = 14：最終的に求まった最小のサンプルサイズ（各群で14例、合計28例）
  • 自由度 = 26
  • t critical = 1.7056：α = 0.10 に対応する両側のt臨界値
  • δ = 2.5538：この自由度で目的の検出力を得るための非心度（非心t分布表の値になります）
  • 効果量 d ≒ 0.965：δを有効サンプルサイズで割ったもの
  • λcalc ≒ 2.65：指定のMSWにおける理論上の効果量

収束判定条件である λcalc >= δ を満たしており、このnが厳密なサンプルサイズの最小値と判断されます。

この結果は、近似法では得られない「非心度に基づいた厳密解」として、江島らの手法に沿った実装が正しく機能していることを示しています。

おわりに

非心t分布を用いた厳密なサンプルサイズ設計は、近似計算では見落とされがちな「境界値」を丁寧に拾い上げることができます。

今回は BE 試験を想定し、効果量を log10(1.2)（±20%の等価範囲）として計算しましたが、これはあくまでBE特有の等価範囲設定に基づいたものです。

一般的な分散分析型の検定（例：群間差や交互作用の検出）では、効果量の定義も異なり、MSW（群内分散）が大きい場合は、所定の効果量では収束しないこともあります。

このような場合、効果量の再定義や検出力の見直しが必要になるため、適切な前提設定が不可欠です。

BE試験以外の応用（一般的な群間比較など）については、別記事にて紹介予定です。

この記事が、より信頼できる試験設計の一助となれば幸いです。

Discussion