临床试验中的样本量估算---实践篇

# -*- coding: utf-8 -*-
class SampleSizeCalculator:
    def __init__(self, alpha=0.05, beta=0.1, sigma=None,
                 delta=None, is_one_sample=False, is_one_side=False):

        self.alpha = alpha # 显著水平,假设检验中犯第一类错误的概率
        self.beta = beta # 假设检验中犯第二类错误的概率
        self.power = 1 - beta # 检验效能
        self.sigma = sigma # 标准差
        self.delta = delta # 可容忍误差
        self.is_one_sample = is_one_sample # 单样本 or 双样本
        self.is_one_side = is_one_side # 单侧检验 or 双侧检验

    # Method1: 通过公式计算样本量
    def get_sample_size1(self):
        import math
        from scipy.stats import norm, t

        # sigma为总体方差时
        if self.is_one_side:
            z_percent_alpha = norm.ppf(self.alpha, 0, 1)
        else:
            z_percent_alpha = norm.ppf(self.alpha / 2, 0, 1)
        z_percent_beta = norm.ppf(self.beta, 0, 1)

        if self.is_one_sample:
            z_test_num = (z_percent_alpha + z_percent_beta) ** 2 * (self.sigma / self.delta) ** 2
            z_test_num = math.ceil(z_test_num)
            if self.is_one_side:
                t_percent_alpha = t.ppf(self.alpha, z_test_num - 1)  # 单样本T检验的自由度为n-1
            else:
                t_percent_alpha = t.ppf(self.alpha / 2, z_test_num - 1)  # 单样本T检验的自由度为n-1
            t_percent_beta = t.ppf(self.beta, z_test_num - 1)
            t_test_num = (t_percent_alpha + t_percent_beta) ** 2 * (self.sigma / self.delta) ** 2
        else:
            z_test_num = 2 * (z_percent_alpha + z_percent_beta) ** 2 * (self.sigma / self.delta) ** 2
            z_test_num = math.ceil(z_test_num)
            if self.is_one_side:
                t_percent_alpha = t.ppf(self.alpha, 2*z_test_num - 2)  # 双样本T检验的自由度为n1+n2-2
            else:
                t_percent_alpha = t.ppf(self.alpha/2, 2*z_test_num - 2)  # 双样本T检验的自由度为n1+n2-2

            t_percent_beta = t.ppf(self.beta, 2*z_test_num - 2)
            t_test_num = 2 * (t_percent_alpha + t_percent_beta) ** 2 * (self.sigma / self.delta) ** 2

        return t_test_num



    # Method2: 通过三方库计算样本量
    def get_sample_size2(self, k):
        if self.is_one_sample:
            from statsmodels.stats.power import TTestPower
            ttest = TTestPower()
            if self.is_one_side:
                alternative = "larger"
            else:
                alternative = "two-sided"
            t_test_num = ttest.solve_power(effect_size=self.delta / self.sigma,
                                           alpha=self.alpha,
                                           power=self.power,
                                           alternative=alternative)
        else:
            from statsmodels.stats.power import TTestIndPower
            ttest = TTestIndPower()
            if self.is_one_side:
                alternatives = "larger"
            else:
                alternative = "two-sided"
            t_test_num = ttest.solve_power(effect_size=self.delta / self.sigma,
                                           nobs1=None,
                                           alpha=self.alpha,
                                           power=self.power,
                                           ratio=k,
                                           alternative=alternative)
        return t_test_num


if __name__ == "__main__":
    alpha = 0.05
    beta = 0.02
    sigma = 8
    delta = 5
    is_one_sample = True
    is_one_side = True

    sample_size_calculator = SampleSizeCalculator(alpha=alpha,
                                                  beta=beta,
                                                  sigma=sigma,
                                                  delta=delta,
                                                  is_one_sample=is_one_sample,
                                                  is_one_side=is_one_side)

    sample_size1 = sample_size_calculator.get_sample_size1()
    print(sample_size1)
    sample_size2 = sample_size_calculator.get_sample_size2(k=1)
    print(sample_size2)