tooluniverse-statistical-modeling

生物医学数据分析的统计建模

规则零 — 首先检查预计算结果

在遵循以下任何指示之前，请先扫描数据文件夹，查找：

• *_executed.ipynb → 使用 tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' 读取，并将其单元格输出作为权威答案引用
• 预计算结果文件（名称类似 *results*、*deseq*、*enrich*、*stats*、*_simplified.csv 的 CSV/TSV 文件）→ 直接读取并报告所需数值
• 规范分析脚本（analysis.R、run_*.py、find_*.R、*.Rmd）→ 原样执行并读取输出

仅当上述文件均不存在时，才按照下方此技能的重新分析流程操作。从原始数据重新运行得到的数值会与已发布答案不同，且速度慢得多（所需轮数通常是 5–10 倍）。

---

主要脚本 — 优先使用

这些脚本在 scripts/ 中编码了针对具体问题的注意事项，并输出带标签、可解析的结果。优先使用它们，而不是临时编写 statsmodels / scipy 代码。

脚本	适用场景
r_natural_spline_regression.py	任何提及 R 语法 lm(y ~ ns(x, df = K))、“自然样条”或要求提供样条 R²/F/峰值预测置信区间的问题。始终调用 Rscript，以确保与 splines::ns() 的结果一致。
spline_model_compare.py	“二次、三次和自然样条中拟合效果最佳的模型” / “最优 x 处的最大菌落面积”。在 R 中拟合全部三种模型，按调整后 R²/AIC/BIC 排序，并报告最佳模型的峰值 (x*、y*) 及其 95% 置信区间。
logistic_regression_or.py	二元或有序 Logistic 回归，且答案为 OR（或 OR + 95% 置信区间）时。处理标签编码、显式的 Placebo=0/BCG=1 映射以及交互项（--interaction A:B → 创建 A_B = A*B）。为每个系数输出 OR + 置信区间，并为指定的 --coef-name 输出 SCALARS 区块。
power_analysis.py	“每组最小样本量”、“TTestIndPower”或“给定 Cohen’s d，power=0.8 时的 N”。从 CSV 计算合并标准差的 Cohen’s d（或接受 --effect-size），然后调用 TTestIndPower.solve_power。
expression_anova.py	跨细胞类型或样本组进行逐基因 ANOVA / 中位数 LFC（不要跨基因汇总 — 见下方警告）。
prepare_ae_cohort.py	SDTM DM/AE 文件上的临床试验 AE 严重程度检验（卡方 / 有序）；使用 encoding='latin1'，并计算每位受试者所有 AE 的 max(AESEV)（不筛选 AEPT）。
stat_tests.py	仅使用标准库的卡方拟合优度检验、Fisher 精确检验和简单 OLS。适用于 scipy/statsmodels 不可用的情况。

具体调用示例

自然样条回归（R^2、整体 F 检验 p 值、峰值 Y + 95% 置信区间）：

python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
  --csv data.csv --y-col Area \
  --ratio-col Ratio --new-x-col Frequency_strain \
  --filter "StrainNumber not in ['1', '98']" \
  --df 4 --workdir /tmp/spline_run

二次 vs 三次 vs 自然样条比较 + 最佳模型峰值：

python skills/tooluniverse-statistical-modeling/scripts/spline_model_compare.py \
  --csv data.csv --y-col Area \
  --ratio-col Ratio --new-x-col Frequency_strain \
  --filter "StrainNumber not in ['1', '98']" \
  --ns-df 4 --workdir /tmp/spline_cmp

**报告峰值位置（x*）时，使用拟合 x 变量的单位，而不是派生标签。**如果模型使用频率/比例列（例如 Frequency_strain，取值为 0–1）进行拟合，那么对于“最大值对应的比例/频率是多少”这一问题，答案应是该分数（例如 0.909），而不是由其推导出的冒号比值（例如 10:1）。如果问题要求 0–1 的值，或拟合的 x 列是分数，则将冒号比值 a:b 转换为分数 a/(a+b)。

含交互项的有序 Logistic 回归（例如试验 AE 严重程度）：

python skills/tooluniverse-statistical-modeling/scripts/logistic_regression_or.py \
  --csv merged.csv --outcome AESEV --outcome-type ordinal --outcome-order "1,2,3,4" \
  --predictors TRTGRP,expect_interact,patients_seen,MHONGO \
  --encode TRTGRP,expect_interact,patients_seen \
  --encode-map "TRTGRP:Placebo=0,BCG=1" \
  --interaction MHONGO:TRTGRP_cat \
  --coef-name TRTGRP_cat

基于试点 CSV 的双样本功效分析：

python skills/tooluniverse-statistical-modeling/scripts/power_analysis.py \
  --csv pilot.csv --value-col MeasuredValue --group-col Group \
  --group-a Treatment --group-b Control \
  --power 0.8 --alpha 0.05

---

工作区隔离（至关重要）

任何分析的输入数据文件夹都必须保持不变，以确保能够复现分析。会写入中间文件（R 驱动脚本、预处理后的 CSV、比较表）的脚本必须将文件写入 /tmp/，或写入通过 --workdir 指定的目录。此技能中的两个基于 R 的脚本都会拒绝运行：如果 --workdir 解析后指向输入 CSV 的父目录（或其任何祖先目录）。

# OK
--workdir /tmp/spline_run

# Refused:
--workdir <path-equal-to-or-containing-the-input-csv>/...

---

重要 — 编写任何代码前请阅读

1. 临床试验 AE 分析（回归、卡方检验、任何严重程度检验）：使用捆绑脚本（或封装该脚本的 clinical_trial_ae_severity_test ToolUniverse 工具）：

   tu run clinical_trial_ae_severity_test '{"dm_file":"DM.csv","ae_file":"AE.csv","test":"chi-square","group_col":"TRTGRP"}'
   
   # Or directly:
   python skills/tooluniverse-statistical-modeling/scripts/prepare_ae_cohort.py \
     --dm DM.csv --ae AE.csv --test chi-square --group TRTGRP \
     --subgroup "expect_interact=Yes"   # optional
   

   该脚本/工具会处理：SDTM CSV 使用 encoding='latin1'；对每位受试者取所有 AE 中的 max(AESEV)（不按 AEPT 筛选）；与 DM 执行内连接；可选亚组筛选；可选使用协变量进行有序 Logistic 回归。

   为什么不按 AEPT 筛选 — AESEV 是 AE 表中由方案定义的严重程度量表。按 AEPT 筛选 AE（例如只保留 AEPT == "COVID-19"）会排除严重程度最高的事件记录在其他 AEPT 标签下的受试者，大幅改变列联表，甚至可能使检验结果反转。短语“COVID-19 严重程度”描述的是结局，而不是筛选条件。

   • ❌ 错误：ae[ae['AEPT'].str.contains('COVID-19')].groupby('USUBJID')['AESEV'].max() — 仅筛选 COVID-19 事件
   • ✅ 正确：ae.groupby('USUBJID')['AESEV'].max() — 使用所有 AE 记录

2. 多特征数据的表达量 ANOVA／倍数变化（基因 × 样本矩阵）： 如果要得到单个“F 统计量”或“倍数变化”，应按基因分别运行后再汇总 — 切勿将所有基因的 expr.values.ravel() 合并。

   • 对于 F 统计量：推导每个样本的量（例如 DESeq2 计算两种细胞类型之间每个基因的 LFC，然后对各组的这些 LFC 做 ANOVA），或者仅针对单个目标基因运行。
   • 对于两组之间的 中位数／平均 log2 倍数变化：使用 design=~<group> 运行 DESeq2，提取每个基因的 log2FoldChange（如果流程使用收缩，则应用收缩），然后计算所有基因的中位数／平均值。

   ❌ 错误（汇总）：按样本计算 log2(sum_counts_groupA / sum_counts_groupB) 后再汇总 — 这得到的是总量之比，会被高表达基因主导。 ✅ 正确（按基因）：DESeq2 → results_df['log2FoldChange'].median()。

   合理性检查经验法则：对于少数几组的生物学 ANOVA，F > 50 意味着你做了汇总（典型的生物学 F 值为 0.5–10）。相似组之间 |中位数 LFC| > 2 意味着你做了汇总（典型的 |中位数| < 1）。

   使用捆绑脚本：python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py（或 expression_anova_per_gene ToolUniverse 工具）。

3. 样条模型 — R splines::ns(x, df=K) ≠ Python patsy.dmatrix("cr(x, df=K)")。由于内部结点位置、边界结点位置和基函数正交化方式不同，它们生成的设计矩阵也不同。任何提及 R 语法（如 lm(y ~ ns(x, df = 4))）的问题，都应通过 Rscript 运行 R。使用捆绑封装器：

   python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
     --csv data.csv --y-col Y --x-col X --df 4 --workdir /tmp/spline_run
   

   对于“菌株 X 频率”的共培养模型，纳入纯目标菌株（freq=1），但排除非目标纯菌株（freq=0）。

4. CSV 编码：临床试验 CSV 通常需要 encoding='latin1'。

5. 计数类变量与长度类变量之间的 Pearson 相关性：当一个变量的取值跨越多个数量级（原始读取计数、TPM、基因长度、转录本丰度）时，即使对数转换后的 r 为中等水平，原始 Pearson r 也往往接近 0。最终答案正文中务必计算并明确报告全部四种变体： r(x, y)、r(log10(x+1), y)、r(x, log10(y+1))、r(log10(x+1), log10(y+1))。 以表格列出，并标明一个作为首选。已发表的答案可能是其中任意一个，问题文本很少会明确使用了哪种变换组合。

   ## Primary answer: r = X.XXX (transform: <name>)
   
   ## Sensitivity (all 4 transform combinations):
   - r(x, y)                   = ...
   - r(log10(x+1), y)          = ...
   - r(x, log10(y+1))          = ...
   - r(log10(x+1), log10(y+1)) = ...
   

   背景 — 对于任何单一变换版本：

   import numpy as np
   from scipy.stats import pearsonr
   r_raw, _ = pearsonr(x, y)
   r_log, _ = pearsonr(np.log10(x + 1), y)
   print(f"r_raw={r_raw:.4f}  r_log10={r_log:.4f}")
   

   默认规则：

   • 问题提到“对数转换”／“对数表达量” → 报告 r_log10
   • 问题未指定，但变量是基因表达量／RNA 计数 → 同样报告 r_log10 作为标准答案（大多数已发表的基因长度与表达量相关性都采用对数尺度）
   • 当 |r_raw| < 0.1 且 |r_log10| > 0.2 时，优先选择 r_log10

   ❌ 错误：只报告 r_raw ≈ 0.05，而对数变换后的值为 0.35 ✅ 正确：“r_raw = 0.05；r_log10 = 0.35（对数分布表达量的标准结果）”

---

计算，而不是描述

每次进行统计分析，都要通过 Bash 编写并运行 Python 代码。不要只描述你「会怎么做」——直接执行。使用 pandas 处理数据，使用 statsmodels 进行回归，使用 scipy 进行检验，使用 matplotlib 绘图。执行代码并报告实际数值（β、p 值、置信区间、N）。

查证，而不是猜测

对任何不确定的科学事实，先搜索数据库，而不是凭记忆推理。

功能

• 线性回归 - 对连续结局进行 OLS 回归，并执行诊断检验
• 逻辑回归 - 二元、有序和多项模型，并计算优势比
• 生存分析 - Cox 比例风险模型和 Kaplan-Meier 曲线
• 混合效应模型 - 用于层级数据或重复测量数据的 LMM/GLMM
• ANOVA - 单因素/双因素 ANOVA，以及用于组学数据的逐特征 ANOVA
• 模型诊断 - 假设检验、拟合统计量、残差分析
• 统计检验 - t 检验、卡方检验、Mann-Whitney 检验、Kruskal-Wallis 检验等

适用场景

在用户提出以下请求时应用此技能：

• 「X 与 Y 相关的优势比是多少？」
• 「治疗的风险比是多少？」
• 「对 Y 关于 X1、X2、X3 拟合线性回归」
• 「对严重程度结局进行有序逻辑回归」
• 「时间 T 时的 Kaplan-Meier 生存估计是多少？」
• 「调整混杂因素后，优势比降低了百分之多少？」
• 「运行带随机截距的混合效应模型」
• 「计算 A 和 B 之间的交互项」
• 「比较各组的 ANOVA F 统计量是多少？」
• 「检验基因/miRNA 表达是否在不同细胞类型间存在差异」

模型选择决策树

START: What type of outcome variable?
|
+-- CONTINUOUS (height, blood pressure, score)
|   +-- Independent observations -> Linear Regression (OLS)
|   +-- Repeated measures -> Mixed-Effects Model (LMM)
|   +-- Count data -> Poisson/Negative Binomial
|
+-- BINARY (yes/no, disease/healthy)
|   +-- Independent observations -> Logistic Regression
|   +-- Repeated measures -> Logistic Mixed-Effects (GLMM/GEE)
|   +-- Rare events -> Firth logistic regression
|
+-- ORDINAL (mild/moderate/severe, stages I/II/III/IV)
|   +-- Ordinal Logistic Regression (Proportional Odds)
|
+-- MULTINOMIAL (>2 unordered categories)
|   +-- Multinomial Logistic Regression
|
+-- TIME-TO-EVENT (survival time + censoring)
    +-- Regression -> Cox Proportional Hazards
    +-- Survival curves -> Kaplan-Meier

工作流

阶段 0：数据验证

目标：加载数据，识别变量类型，并检查缺失值。

关键：先识别结局变量

在进行任何分析之前，先确认你实际要预测的是什么：

1. 通读完整问题 - 查找「预测 [结局]」、「建模 [结局]」或「因变量」等表述
2. 检查可用列 - 列出数据集中的所有列
3. 将问题与数据匹配 - 找出与所述结局相对应的列
4. 确认结局存在 - 不要根据预测变量创建结局变量

常见错误：问题提到「肥胖」→ 假定结局 = BMI >= 30（这会导致 BMI 预测变量与结局之间产生循环逻辑）。务必先检查数据列：print(df.columns.tolist())

import pandas as pd
import numpy as np

df = pd.read_csv('data.csv')
print(f"Observations: {len(df)}, Variables: {len(df.columns)}, Missing: {df.isnull().sum().sum()}")

for col in df.columns:
    n_unique = df[col].nunique()
    if n_unique == 2:
        print(f"{col}: binary")
    elif n_unique <= 10 and df[col].dtype == 'object':
        print(f"{col}: categorical ({n_unique} levels)")
    elif df[col].dtype in ['float64', 'int64']:
        print(f"{col}: continuous (mean={df[col].mean():.2f})")

阶段 1：模型拟合

目标：根据结局类型拟合适当的模型。

使用上方的决策树选择模型类型，然后参考相应的参考文件查看详细代码：

• 线性回归：references/linear_models.md
• 逻辑回归（二元）：references/logistic_regression.md
• 有序逻辑回归：references/ordinal_logistic.md
• Cox 比例风险模型：references/cox_regression.md
• ANOVA / 统计检验：anova_and_tests.md

常用模型速查：

import statsmodels.formula.api as smf
import numpy as np

# Linear regression
model = smf.ols('outcome ~ predictor1 + predictor2', data=df).fit()

# Logistic regression (odds ratios)
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)
ors = np.exp(model.params)
ci = np.exp(model.conf_int())

# Cox proportional hazards
from lifelines import CoxPHFitter
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age']], duration_col='time', event_col='event')
hr = cph.hazard_ratios_['treatment']

确定性工具替代方案——对于不需要自定义 lifelines 诊断的 Kaplan-Meier / log-rank / Cox 问题，Survival_kaplan_meier、Survival_log_rank_test 和 Survival_cox_regression（ToolUniverse 工具，纯计算，无需导入 lifelines）可通过一次调用给出相同的统计量；对于一次性数值计算，也比手动编写代码更不容易出错：

tu run Survival_kaplan_meier '{"durations":[5,10,15,20,25,30,12,8,22,18],"event_observed":[1,0,1,1,0,1,1,0,0,1]}'
# -> median_survival_time, per-time survival_table (survival_probability,
#    at_risk, ci_lower_95/ci_upper_95 via Greenwood log-log), n_events, n_censored

tu run Survival_log_rank_test '{"durations_a":[...],"events_a":[...],"durations_b":[...],"events_b":[...]}'
# -> chi2_statistic, p_value, observed/expected events per group, interpretation string

tu run Survival_cox_regression '{"durations":[...],"event_observed":[...],"covariates":{"age":[...],"stage":[...]}}'
# -> per-covariate hazard_ratio + 95% CI + p_value, log_likelihood, convergence flag

典型工作流：先使用 Survival_kaplan_meier 获取描述性曲线/中位生存期，再使用 Survival_log_rank_test 检验两组是否存在差异，最后使用 Survival_cox_regression 调整协变量并计算风险比（已实时验证：stage 的 HR=1.71，95% CI 为 0.46–6.28；当每个协变量对应的事件数低于 10 时，同时出现了 low_epv_warning——务必指出这一警告，不要单独将 HR 报告为可靠结果）。对于这些工具未涵盖的手动编写 lifelines 诊断（比例风险检验、Schoenfeld 残差），请使用 references/cox_regression.md。

阶段 1b：多特征数据的 ANOVA

当数据包含多个特征（基因、miRNA、代谢物）时，使用逐特征 ANOVA（而非聚合分析）。这是基因组学中最常见的模式。

完整的决策树、两种方法和示例，请参阅 anova_and_tests.md。

基因表达数据的默认方法：逐特征 ANOVA（方法 B）。

阶段 2：模型诊断

目标：检查模型假设和拟合质量。

不同模型类型的关键诊断：

• OLS：Shapiro-Wilk（正态性）、Breusch-Pagan（异方差性）、VIF（多重共线性）
• Cox：通过 cph.check_assumptions() 进行比例风险检验
• Logistic：Hosmer-Lemeshow、ROC/AUC

诊断代码和常见问题请参阅 references/troubleshooting.md。

阶段 3：结果解读

目标：生成适合发表的摘要。

对于每项结果，都要报告：效应量（OR/HR/系数）、95% CI、p 值和模型拟合统计量。常见问答模式请参阅 common_patterns_summary.md。

常见模式

模式	问题类型	关键步骤
1	有序回归的比值比	拟合 OrderedModel，计算 exp(coef)
2	OR 的百分比降幅	比较粗略模型与调整后模型
3	交互效应	拟合 A * B，提取 A:B 系数
4	风险比	Cox PH 模型，计算 exp(coef)
5	多特征 ANOVA	逐特征 F 统计量（而非聚合分析）

每种模式的解决方案代码请参阅 common_patterns_summary.md。 15 种以上详细的问题模式请参阅 references/common_patterns.md。

Statsmodels 与 Scikit-learn

使用场景	库	原因
推断（p 值、CI、OR）	statsmodels	提供完整的统计输出
预测（准确率、AUC）	scikit-learn	预测工具更完善
混合效应模型	statsmodels	唯一选择
正则化（LASSO、Ridge）	scikit-learn	优化效果更好
生存分析	lifelines	专用库

一般规则：对于统计推断问题（p 值、OR、HR），使用 statsmodels。

Python 包要求

statsmodels>=0.14.0
scikit-learn>=1.3.0
lifelines>=0.27.0
pandas>=2.0.0
numpy>=1.24.0
scipy>=1.10.0

核心原则

1. 数据优先——建模前始终检查并验证数据
2. 根据结局类型选择模型——使用上述决策树
3. 检查假设——验证模型假设（线性关系、比例风险等）
4. 完整报告——始终报告效应量、CI、p 值和模型拟合统计量
5. 关注混杂因素——在指定或具有临床意义时调整混杂因素
6. 可复现分析——所有代码都必须是确定性的且可复现
7. 稳健的错误处理——妥善处理收敛失败、分离和共线性
8. 正确舍入——遵循要求的精度（通常为小数点后 2–4 位）

结果解读的推理框架

证据分级

等级	标准	示例
强	p < 0.001，效应量具有临床意义，满足模型假设	OR = 3.5（95% CI：2.1-5.8），p < 0.001，Hosmer-Lemeshow p > 0.05
中等	p < 0.05，效应量合理，假设存在轻微问题	HR = 1.8（95% CI：1.1-2.9），p = 0.02，PH 检验结果临界
弱	p < 0.05，但 CI 较宽、效应较小或违反假设	OR = 1.2（95% CI：1.01-1.43），p = 0.04，某个协变量的 VIF > 5
证据不足	p >= 0.05，或模型未能收敛/通过诊断	系数不显著，且模型出现分离警告

解读指南

• 模型诊断（R-squared）：对于 OLS，在生物医学数据中，R-squared > 0.7 表示拟合良好；0.3–0.7 表示拟合中等。调整后 R-squared 会对新增预测变量进行惩罚。对于逻辑模型，请使用伪 R-squared（McFadden > 0.2 表示可接受）和 AUC（> 0.7 = 可接受，> 0.8 = 区分能力良好）。
• 使用 AIC/BIC 比较模型：数值越低越好。模型间的 AIC 差异 > 10，强烈支持选择 AIC 较低的模型。BIC 对复杂度的惩罚比 AIC 更重，因此更倾向于简单模型。以预测为重点进行选择时使用 AIC，以推断为重点时使用 BIC。
• 系数显著性阈值：请报告精确的 p 值，而非仅使用显著性星号。对于多个预测变量，请进行 Bonferroni 或 FDR 校正。在包含 20 个预测变量的模型中，若未经校正的系数 p = 0.049，很可能是假阳性。
• 生存分析中 HR 的解读：HR > 1 表示暴露组的风险增加（结局更差）。HR = 2.0 表示事件发生的瞬时风险是原来的两倍。务必检验比例风险假设——如果该假设不成立，HR 就是随时间变化的平均值，可能会产生误导。请同时报告中位生存时间和 HR，以便临床解读。
• 比值比的解读：OR = 1.0 表示无关联。OR > 1 表示比值增加。95% CI 不包含 1.0 即可确认统计学显著性。对于罕见结局，OR 近似于相对风险；对于常见结局（患病率 > 10%），OR 会高估相对风险。
• 混杂评估：比较粗略 OR/HR 与调整后 OR/HR。对协变量进行调整后，效应估计值变化 > 10% 表示该变量可能存在混杂。

综合问题

1. 模型诊断（残差图、Hosmer-Lemeshow 检验、PH 检验）是否支持所选模型的有效性，还是假设违背要求采用其他方法（例如稳健标准误、分层模型）？
2. 对于调整模型，纳入混杂因素后，主要效应估计值的变化是否超过 10%，表明存在有意义的混杂？
3. 考虑到预测变量和结局的量级，所报告的效应量（OR、HR、系数）除了具有统计学显著性之外，在临床上是否也有意义？
4. 使用 AIC/BIC 比较嵌套模型时，更复杂的模型是否显著改善了拟合，还是更简单的模型更符合简约原则？
5. 对于生存分析，整个随访期间是否均满足比例风险假设，还是 Schoenfeld 残差提示效应随时间变化？

---

完整性检查清单

完成任何统计分析前：

• 已确定结局变量：已确认哪一列是真正的结局变量
• 已验证数据：已确认样本量 N、缺失值和变量类型
• 已识别多特征数据：如果存在多个特征，请按特征分别处理
• 模型适用：结局类型与模型族相匹配
• 已检查假设：已执行相关诊断
• 已报告效应量：OR/HR/Cohen's d 及其 CI
• 已报告 P 值：必要时进行了适当校正
• 已评估模型拟合度：R-squared、AIC/BIC、一致性
• 已解读结果：以通俗语言进行解读
• 精度正确：数字已按适当精度取整

捆绑脚本

这些可直接运行的脚本位于 skills/tooluniverse-statistical-modeling/scripts/。 请通过 Bash 工具使用它们——对于上述文档中反复出现的问题模式，这些脚本可提供确定性答案。

r_natural_spline_regression.py — R 中的自然样条回归

调用 Rscript，使用 splines::ns() 拟合 lm(y ~ ns(x, df=K))。 输出 R²、调整后 R²、含 df1/df2 的 F 统计量、整体 F 检验 p 值、残差标准误、 系数表（估计值、标准误、t 值、p 值），以及预测网格上的峰值和 来自 predict.lm(..., interval='confidence') 的 95% 置信区间。支持 使用 --ratio-col 快捷选项，将“a:b”字符串比例转换为频率分数 a/(a+b)。 拒绝将文件写入输入 CSV 的父目录。

spline_model_compare.py — 二次模型、三次模型与自然样条的比较

在 R 中对相同的 x、y 拟合全部三个模型，按调整后 R²、AIC 和 BIC 排序，并输出最佳模型的峰值 (x*，y*) 及其 95% 置信区间。适用于 “拟合效果最佳的模型”以及“最优 x 处预测 y 的最大值”等问题。

logistic_regression_or.py — 二元或有序 Logistic 回归及 OR

拟合 sm.Logit（二元）或 OrderedModel（有序比例优势模型），并为每个系数 输出 OR（exp(coef)）、95% 置信区间和 p 值。 支持标签编码（--encode A,B,C）、显式值映射 （--encode-map TRTGRP:Placebo=0,BCG=1）以及交互列 （--interaction A:B -> 创建 A_B = A*B）。使用 --coef-name <NAME> 还会输出一个标记为所请求系数的 SCALARS 块。

power_analysis.py — t 检验的双样本所需样本量

根据 CSV 中的 --value-col、--group-col、--group-a、--group-b 计算 Cohen's d（合并标准差）， 然后调用 TTestIndPower.solve_power，并使用 --alpha、 --power、--alternative 参数。适用于“每组最小样本量”这类效能分析 问题。返回原始 N 和向上取整后的 N。

stat_tests.py — 基础统计检验（纯标准库，无需 scipy）

实现卡方拟合优度检验、Fisher 精确检验和简单线性回归， 无需任何外部依赖。所有 p 值均从原理出发计算，使用 伽马函数（卡方检验）或超几何枚举（Fisher 检验）。

# Chi-square goodness-of-fit
python stat_tests.py --type chi_square --observed 100,50,25 --expected 87.5,50,37.5

# Fisher's exact test (2×2 table)
python stat_tests.py --type fisher_exact --a 10 --b 5 --c 3 --d 20
python stat_tests.py --type fisher_exact --a 10 --b 5 --c 3 --d 20 --alternative greater

# Simple linear regression (OLS)
python stat_tests.py --type regression --x "1,2,3,4,5" --y "2.1,4.0,5.9,8.1,10.0"

关键公式：

• chi_square：χ² = Σ (O−E)²/E；通过上正则化不完全伽马函数 Q(df/2, χ²/2) 计算 p 值
• fisher_exact：超几何 PMF；p 值 = 概率 ≤ P(观测值) 的所有概率之和
• regression：b1 = Sxy/Sxx；b0 = ȳ − b1x̄；R² = 1 − SSR/SST；包含 SE 和 t 统计量

输出包括：完整的列联表/数据表、逐步算式、显著性说明， 以及每项检验的往返验证。

何时使用 stat_tests.py，何时使用 statsmodels：

• 需要快速健全性检查且无需导入模块，或环境中没有 scipy/statsmodels 时，使用 stat_tests.py。
• 需要多变量回归、Logistic 模型或生存分析时，使用 statsmodels。

format_statistical_output.py — 格式化报告结果

用于将拟合后的 statsmodels 结果格式化为可用于发表的表格的实用函数。 从分析脚本中导入并调用；不是独立的 CLI 工具。

model_diagnostics.py — 自动化模型诊断

对拟合模型执行假设检验（正态性、异方差性、多重共线性）。 从分析脚本中导入并调用；不是独立的 CLI 工具。

---

文件结构

tooluniverse-statistical-modeling/
+-- SKILL.md                          # This file (workflow guide)
+-- QUICK_START.md                    # 8 quick examples
+-- EXAMPLES.md                       # Legacy examples
+-- TOOLS_REFERENCE.md                # ToolUniverse tool catalog
+-- anova_and_tests.md                # ANOVA decision tree and code
+-- common_patterns_summary.md         # Common solution patterns
+-- test_skill.py                     # Test suite
+-- references/
|   +-- logistic_regression.md        # Detailed logistic examples
|   +-- ordinal_logistic.md           # Ordinal logit guide
|   +-- cox_regression.md             # Survival analysis guide
|   +-- linear_models.md              # OLS and mixed-effects
|   +-- common_patterns.md            # 15+ question patterns
|   +-- troubleshooting.md            # Diagnostic issues
+-- scripts/
    +-- r_natural_spline_regression.py # lm(y ~ ns(x, df=K)) via Rscript
    +-- spline_model_compare.py        # quadratic vs cubic vs natural-spline (Rscript)
    +-- logistic_regression_or.py      # binary / ordinal logistic + ORs + interactions
    +-- power_analysis.py              # TTestIndPower required-N from CSV
    +-- expression_anova.py            # per-gene ANOVA / log2FC summary
    +-- prepare_ae_cohort.py           # SDTM AE/DM cohort prep
    +-- stat_tests.py                  # Chi-square, Fisher's exact, OLS (stdlib)
    +-- format_statistical_output.py   # Format results for reporting
    +-- model_diagnostics.py           # Automated diagnostics

ToolUniverse 集成

此技能主要用于计算，但 ToolUniverse 工具也可提供数据：

使用场景	工具
临床试验数据	search_clinical_trials
药物安全性结局	FAERS_calculate_disproportionality
基因-疾病关联	OpenTargets_target_disease_evidence
生物标志物数据	fda_pharmacogenomic_biomarkers

完整工具目录请参阅 TOOLS_REFERENCE.md。

参考资料

• statsmodels：https://www.statsmodels.org/
• lifelines：https://lifelines.readthedocs.io/
• scikit-learn：https://scikit-learn.org/
• 有序模型：statsmodels.miscmodels.ordinal_model.OrderedModel

分析惯例

这些惯例是经过验证的最佳实践。适用于数据集/问题与之匹配的情况。

强制要求：表达量 ANOVA / 倍数变化分析请使用捆绑脚本

针对按基因进行 ANOVA 或中位数 log2FC 的问题，请使用捆绑脚本：

python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --exclude-groups PBMC --mode anova

python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --group-a CD14 --group-b CD19 --mode fold_change

不要自行编写 pandas ANOVA——聚合层级（按基因，而非按样本）至关重要，而且很容易弄错。

CSV 编码

临床试验导出文件（SDTM）通常采用 latin1 编码。如果 pd.read_csv() 因 UnicodeDecodeError 失败，请使用 encoding='latin1' 重试。

临床试验 AE 分析——适用于回归、卡方检验及任何严重程度检验

针对临床试验 AE 严重程度结局与试验协变量之间的任何统计检验（卡方检验、有序 / 逻辑回归、Mann-Whitney 等），请在报告了任何 AE 的队列上拟合 / 检验（并与人口统计数据进行内连接）。不要按 AEPT 预先筛选 AE 记录——即使问题提到“COVID-19 严重程度”或任何特定病症，也不要这样做。对于每位受试者，请使用其所有 AE 记录中的 max(AESEV)，不论 AEPT 为何。

原因：AE 表中的 AESEV 反映了研究方案定义的严重程度量表。按特定 AEPT 值预先筛选（例如只保留某些病症标签），会剔除那些最严重程度记录在其他 AEPT 标签下的受试者，从而大幅改变列联表，并可能使结果从显著变为不显著。

不要将没有 AE 的受试者填充为 AESEV=0——这会稀释信号。

dm = pd.read_csv("DM.csv", encoding='latin1')
ae = pd.read_csv("AE.csv", encoding='latin1')
sev = ae.groupby('USUBJID')['AESEV'].max().reset_index()
df = dm.merge(sev, on='USUBJID', how='inner').dropna(subset=['AESEV'])
df['AESEV'] = df['AESEV'].astype(int)
# Now ordinal regression / chi-square on df

优势比偏差

“优势比降低百分比”指 (1 − OR) × 100%——即相对于 OR=1（无效应）的偏差。OR=0.75 → 相对于参照值降低 25%。不要将其解释为 (unadjusted_OR − adjusted_OR) / unadjusted_OR；由于调整几乎不会使设定良好的 OR 发生变化，该值几乎总是 ≈0%。

F 统计量与 p 值

scipy.stats.f_oneway(g1, g2, ...) 返回 (F, p)。如果 GT 看起来像 (0.76, 0.78)，而你计算出的 F=91.6，那么你返回的是 F 统计量，但问题问的是 p 值（或反之）。务必重新阅读问题，确认答案需要哪一个。

对不同组的表达量进行 ANOVA——聚合方式很重要

当问题询问“比较不同细胞类型 / 组之间表达量的 F 统计量”时：

• 每个基因 / miRNA 是一个观测值。对于 K 个组中的 N 个基因，每组有 N 个值（该基因在该组样本中的平均或中位表达量）
• 运行 f_oneway(group1_values, group2_values, ...)，其中每组包含 N 个基因级别的值
• 不要先按样本对所有基因求和——这会得到总 RNA 含量，是完全不同的量
• 当大多数基因在组间没有差异时，F 统计量通常较低（0.5–2.0），而不是较高（50+）
• 如果 F 统计量 >10，但生物学背景表明“没有显著差异”，你可能按样本而不是按基因进行了聚合

组间 log2 倍数变化——先按基因计算，再汇总

当问题询问“组 A 与组 B 之间的中位数 log2 倍数变化”时：

• 按基因计算 log2FC：对每个基因，计算 log2(mean_expr_A / mean_expr_B)
• 然后对所有基因取中位数
• 不要先按样本对所有基因求和——这会得到单个总表达量比值，而不是按基因计算的中位数倍数变化
• 中位数 log2FC 接近 0，表示大多数基因在两组中的表达量相近（当两组属于相似细胞类型时，这是预期结果）

ANOVA 和倍数变化分析均使用捆绑脚本：

# Per-gene ANOVA across cell types
python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --exclude-groups PBMC --mode anova

# Per-gene log2FC between two groups
python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --group-a CD14 --group-b CD19 --mode fold_change

基于应变共培养数据的自然样条回归

在拟合应变共培养频率数据的模型时：

1. 比例 → 频率转换：如果 Ratio 列包含 "10:1" 这样的字符串，请将其转换为频率分数：first / (first + second) = 0.909。以 [0, 1] 范围内的分数报告，而不是比例表示法。

2. 纯菌株端点：是否纳入纯菌株数据（freq=0 和 freq=1）取决于模型类型：

   • 三次／多项式模型：仅使用共培养行拟合（排除纯菌株）。三次模型描述混合群体的响应曲线；纯菌株属于根本不同的生物学状态，纳入它们通常会降低 R²。
   • 自然样条模型（ns(freq, df=4)）：纳入焦点菌株的纯菌株端点（即模型要预测其频率的菌株），但排除非焦点纯菌株。例如，在建模“ΔrhlI 占总群体的频率”时，纳入纯 ΔrhlI（freq=1.0），但排除纯 ΔlasI（freq=0.0）。这会在焦点菌株占主导的高端锚定样条。

3. R 与 Python 样条：R 的 ns()（来自 splines 包）和 Python 的 patsy.cr() 或 scipy BSpline 产生的节点位置和边界条件不同。如果问题涉及 R 的 lm(y ~ ns(x, df=4))，请使用通过 Rscript 运行 R 的内置包装器：

   python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
     --csv data.csv --y-col Area --x-col Frequency \
     --df 4 --workdir /tmp/spline_run
   

   该包装器会输出 R²、调整后 R²、整体 F 检验 p 值、系数表，以及包含峰值（x*、y*）和峰值处 95% 置信区间的预测网格。

   不要用 Python 的 patsy.dmatrix("cr(x, df=4)") 替代——它会得出不同的 R²、F 统计量和预测值。

4. 峰值处的预测：内置 R 包装器已在精细的 1000 点网格上进行预测，并报告 which.max(pred[,"fit"]) 及其通过 interval="confidence" 得到的置信区间。如果必须手动计算，请遵循相同的步骤；不要使用粗略网格。

5. 在{二次、三次、自然样条}中选最佳模型：如果问题询问“二次、三次和自然样条中拟合效果最佳的模型”或“最佳频率下的最大菌落面积”，请使用在 R 中拟合全部三种模型并进行排名的比较包装器：

   python skills/tooluniverse-statistical-modeling/scripts/spline_model_compare.py \
     --csv data.csv --y-col Area --x-col Frequency \
     --ns-df 4 --workdir /tmp/spline_cmp
   

   输出中的 BEST_BY_ADJ_R2 行会指出胜出的模型，BEST_PEAK_Y 是需要报告的峰值 y。

支持

有关详细示例和故障排除，请参阅：

• 逻辑回归：references/logistic_regression.md
• 有序模型：references/ordinal_logistic.md
• 生存分析：references/cox_regression.md
• 线性／混合模型：references/linear_models.md
• 常见模式：references/common_patterns.md
• ANOVA 和检验：anova_and_tests.md
• 诊断：references/troubleshooting.md