如何用Python实现R同款Bland-Altman统计指标计算

核心需求：已可生成Bland-Altman图，需计算偏倚（Bias）、一致性限（Limits of Agreement, LOA）及对应标准误、95%置信区间，和R的输出结果完全匹配。

R语言参考输出结果：

Standard deviation of bias:  6.52818 

Standard error of bias:  0.06593449 
Standard error for limits of agreement:  0.1126868 

Bias:  -5.855469e-13 
Bias- upper 95% CI:  0.1292452 
Bias- lower 95% CI:  -0.1292452 

Upper limit of agreement:  12.79523 
Upper LOA- upper 95% CI:  13.01612 
Upper LOA- lower 95% CI:  12.57434

计算逻辑 #

所有指标基于两种测量方法的差值序列计算，公式与R常用Bland-Altman分析包默认逻辑完全对齐：

• 差值 d = 方法A测量值 - 方法B测量值，样本量为n
• 偏倚（Bias）：d的均值
• 偏倚的标准差：d的样本标准差（自由度为n-1）
• 偏倚的标准误：偏倚标准差 / sqrt(n)
• 一致性限（LOA）：Bias ± 1.96 * 偏倚标准差
• LOA的标准误：偏倚标准差 * sqrt(3/n)
• 95%置信区间：估计值 ± t分布临界值（自由度n-1，95%置信度） * 对应标准误，大样本下临界值近似为1.96

Python实现代码 #

import numpy as np
from scipy.stats import t

def bland_altman_stats(method1, method2, ci=0.95):
    # 计算两组测量的差值
    d = method1 - method2
    n = len(d)
    # 核心基础指标计算
    bias = np.mean(d)
    sd_d = np.std(d, ddof=1) # 采用样本标准差（自由度n-1）和R默认计算逻辑一致
    se_bias = sd_d / np.sqrt(n)
    loa_upper = bias + 1.96 * sd_d
    loa_lower = bias - 1.96 * sd_d
    se_loa = sd_d * np.sqrt(3/n)
    # 计算t分布临界值
    alpha = 1 - ci
    t_crit = t.ppf(1 - alpha/2, df=n-1)
    # 偏倚的95%置信区间
    bias_ci_upper = bias + t_crit * se_bias
    bias_ci_lower = bias - t_crit * se_bias
    # 上一致性限的95%置信区间
    loa_upper_ci_upper = loa_upper + t_crit * se_loa
    loa_upper_ci_lower = loa_upper - t_crit * se_loa
    # 下一致性限的95%置信区间，如需输出可自行补充打印
    loa_lower_ci_upper = loa_lower + t_crit * se_loa
    loa_lower_ci_lower = loa_lower - t_crit * se_loa

    # 打印与R输出格式完全一致的结果
    print(f"Standard deviation of bias:  {sd_d:.5f}\n")
    print(f"Standard error of bias:  {se_bias:.8f}")
    print(f"Standard error for limits of agreement:  {se_loa:.8f}\n")
    print(f"Bias:  {bias:.7e}")
    print(f"Bias- upper 95% CI:  {bias_ci_upper:.7f}")
    print(f"Bias- lower 95% CI:  {bias_ci_lower:.7f}\n")
    print(f"Upper limit of agreement:  {loa_upper:.5f}")
    print(f"Upper LOA- upper 95% CI:  {loa_upper_ci_upper:.5f}")
    print(f"Upper LOA- lower 95% CI:  {loa_upper_ci_lower:.5f}")

    # 返回字典格式结果方便后续调用
    return {
        "bias": bias,
        "sd_bias": sd_d,
        "se_bias": se_bias,
        "bias_95ci": (bias_ci_lower, bias_ci_upper),
        "loa_upper": loa_upper,
        "loa_lower": loa_lower,
        "se_loa": se_loa,
        "loa_upper_95ci": (loa_upper_ci_lower, loa_upper_ci_upper),
        "loa_lower_95ci": (loa_lower_ci_lower, loa_lower_ci_upper)
    }

调用示例 #

# 替换为实际的两组测量数据序列
method1 = np.array([1.2, 3.5, 5.1, ...])
method2 = np.array([1.3, 3.4, 5.0, ...])
stats_result = bland_altman_stats(method1, method2)

内容的提问来源于stack exchange，提问作者Carla Martins