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 で3種類すべてを当てはめ、調整済み R²/AIC/BIC で順位付けし、最良モデルのピーク（x*、y*）と95％信頼区間を報告します。
logistic_regression_or.py	回答が OR（または OR + 95％信頼区間）となる二値または順序ロジスティック回帰。ラベルエンコーディング、明示的な Placebo=0/BCG=1 の対応付け、交互作用項（--interaction A:B → A_B = A*B を作成）に対応しています。すべての係数について OR + 信頼区間を出力し、指定された --coef-name の SCALARS ブロックも出力します。
power_analysis.py	「グループごとの最小サンプルサイズ」、「TTestIndPower」、「Cohen の d を指定した場合、検出力=0.8に必要な N」。CSV からプールした標準偏差に基づく Cohen の d を計算するか（または --effect-size を受け付け）、TTestIndPower.solve_power を実行します。
expression_anova.py	細胞タイプまたはサンプルグループごとの遺伝子単位の ANOVA / 中央 LFC（遺伝子をまたいだプールではありません。以下の警告を参照）。
prepare_ae_cohort.py	SDTM DM/AE ファイルを使用した臨床試験の AE 重症度検定（カイ二乗／順序）。encoding='latin1' を使用し、AEPT フィルターを使わず、被験者ごとにすべての AE を対象として max(AESEV) を使用します。
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) の割合に変換してください。

交互作用項を含む順序ロジスティック回帰（例：臨床試験の 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 に基づく2標本の検出力分析：

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 との内部結合、任意のサブグループフィルター、共変量を含む任意の順序ロジスティック回帰に対応しています。

   AEPT でフィルターしない理由 — AESEV は AE テーブル上でプロトコルにより定義された重症度尺度です。AE を AEPT（例：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 統計量の場合：サンプルごとの量（例：2 種類の細胞型間における各遺伝子の DESeq2 LFC）を導出し、それらの LFC に対してグループ間で ANOVA を行うか、単一の対象遺伝子に対して実行してください。
   • 2 グループ間の対数 2 フォールドチェンジの中央値／平均値の場合：design=~<group> で DESeq2 を実行し、遺伝子ごとの log2FoldChange を抽出して（パイプラインで使用している場合はシュリンケージを適用）、遺伝子間の中央値／平均値を求めてください。

   ❌ 誤り（集計）：サンプルごとに log2(sum_counts_groupA / sum_counts_groupB) を計算してから要約 — 高発現遺伝子に左右される合計値の比になります。 ✅ 正しい（遺伝子ごと）：DESeq2 → results_df['log2FoldChange'].median()。

   妥当性チェックの目安：いくつかのグループ間の生物学的 ANOVA で F > 50 の場合は、集計してしまった可能性があります（生物学的な F の典型値は 0.5～10 です）。類似したグループ間で |median LFC| > 2 の場合も、集計してしまった可能性があります（|median| の典型値は 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)") は同じではありません。内部ノットの配置、境界ノットの配置、基底の直交化が異なるため、異なるデザイン行列が生成されます。lm(y ~ ns(x, df = 4)) のような R 構文に言及する質問には、必ず Rscript を使ってください。付属のラッパーを使用してください：

   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
   

   「strain X の頻度」を扱う共培養モデルでは、焦点となる株単独（freq=1）は含め、焦点ではない株単独（freq=0）は除外してください。

4. CSV のエンコーディング：臨床試験の CSV では encoding='latin1' が必要になることがよくあります。

5. カウント系変数と長さ系変数の間のピアソン相関：一方の変数が複数桁にわたる場合（未加工リードカウント、TPM、遺伝子長、転写産物量など）、対数変換後の相関が中程度であっても、未加工値のピアソン r は 0 に近くなることがよくあります。最終回答の本文では、必ず次の 4 つの値をすべて計算し、明示してください： r(x, y)、r(log10(x+1), y)、r(x, log10(y+1))、r(log10(x+1), log10(y+1))。 表にして、主要な値として選んだものに印を付けてください。公開済みの回答は 4 つのいずれでもかまいません。質問文だけでは、どの変換の組み合わせが使われたか明確でないことがほとんどです。

   ## 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)) = ...
   

   背景 — 変換を 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 を優先します

   ❌ 誤り：対数変換後が 0.35 なのに、r_raw ≈ 0.05 だけを報告する ✅ 正しい：「r_raw = 0.05；r_log10 = 0.35（対数分布する発現量の標準値）」

---

説明するのではなく、計算する

すべての統計解析で、Python コードを（Bash 経由で）記述して実行してください。「こうするつもりです」と説明するのではなく、実際に実行してください。データ加工には 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 のインポートは不要）を使うと、1 回の 呼び出しで同じ統計量が得られ、単発の数値を手作業のコードで求めるよりもエラーが起きにくくなります：

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 で 2 群間の差を検定し、その後 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()による比例ハザード検定
• ロジスティック：Hosmer-Lemeshow、ROC/AUC

診断コードと一般的な問題については、references/troubleshooting.mdを参照してください。

フェーズ 3：解釈

目標：論文掲載に適した要約を作成します。

すべての結果について、効果量（OR/HR/係数）、95% CI、p 値、モデル適合度統計量を報告します。一般的な質問と回答のパターンについては、common_patterns_summary.mdを参照してください。

一般的なパターン

パターン	質問の種類	主な手順
1	順序回帰からのオッズ比	OrderedModelを適合させ、係数を指数変換
2	ORの減少率	未調整モデルと調整済みモデルを比較
3	交互作用効果	A * Bを適合させ、A:Bの係数を抽出
4	ハザード比	Cox PHモデル、係数を指数変換
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は中程度の適合を示します。Adjusted 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は、イベントの瞬間リスクが2倍であることを意味します。比例ハザードの仮定を必ず検証してください。この仮定に違反している場合、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、欠損値、変数の型を確認済み
• 複数の特徴量を含むデータを特定：特徴量が複数ある場合は、特徴量ごとのアプローチを使用
• 適切なモデルを選択：転帰の型とモデルファミリーが一致
• 仮定を確認：関連する診断を実施済み
• 効果量を報告：CI付きのOR/HR/Cohen's d
• 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 値、残差 SE、 係数テーブル（推定値、SE、t、p）、および predict.lm(..., interval='confidence') による 95% CI 付き予測グリッドのピークを出力します。"a:b" 形式の文字列比率を 頻度分率 a/(a+b) に変換する --ratio-col ショートカットをサポートします。 入力 CSV の親ディレクトリへの書き込みは拒否します。

spline_model_compare.py — 二次 vs 三次 vs 自然スプライン

同じ x,y に対して R で 3 つのモデルすべてをフィットし、調整済み R²、AIC、BIC で 順位付けし、最良モデルのピーク (x*, y*) を 95% CI 付きで出力します。 「最良フィットモデル」の質問や「最適な x における最大予測 y」に使用します。

logistic_regression_or.py — OR 付き二値または順序ロジスティック回帰

sm.Logit（二値）または OrderedModel（順序比例オッズ）をフィットし、 すべての係数について OR（exp(coef)）と 95% CI、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 検定の二標本必要 N

--value-col、--group-col、--group-a、--group-b を指定して CSV から Cohen's d（プール標準偏差）を計算し、続いて --alpha、--power、--alternative を 用いて TTestIndPower.solve_power を実行します。「群あたりの最小サンプルサイズ」の 検出力の質問に使用します。生の N と切り上げ後の N の両方を返します。

stat_tests.py — 基本的な統計検定（純粋な stdlib、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、p 値は上部正則化不完全ガンマ Q(df/2, χ²/2) による
• 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 を使用します。
• 多変量回帰、ロジスティックモデル、または生存時間分析が必要な場合は 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または中央値ログ2FCに関する質問には、同梱スクリプトを使用してください。

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重症度の結果に対するあらゆる統計検定（カイ二乗検定、順序／ロジスティック回帰、マン・ホイットニー検定など）を行う場合は、何らかのAEを報告したコホート（人口統計データとの内部結合後）に対して適合／検定してください。AEPTでAE記録を事前に絞り込まないでください。「COVID-19の重症度」と質問に書かれている場合や、特定の症状が指定されている場合も同様です。AEPTに関係なく、被験者ごとにすべてのAE記録を対象としてAESEVの最大値を使用してください。

理由：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と解釈しないでください。適切に指定されたオッズ比は調整でほとんど変わらないため、この方法ではほぼ常に約0%になります。

F統計量とp値

scipy.stats.f_oneway(g1, g2, ...)は(F, p)を返します。GTが(0.76, 0.78)のように見えるのに、F=91.6を計算した場合、質問がp値を求めているのにF統計量を返した（またはその逆の）可能性があります。回答で何を求められているか、必ず質問を読み直してください。

グループ間の発現量に対するANOVA — 集計が重要

「細胞種／グループ間の発現量を比較するF統計量」を求められた場合：

• 各遺伝子／miRNAは1つの観測値です。グループ数がKで、遺伝子数がNの場合、各グループにはN個の値があります（そのグループのサンプルにおける各遺伝子の平均または中央値の発現量）
• 各グループに遺伝子レベルの値がN個ずつある状態で、f_oneway(group1_values, group2_values, ...)を実行します
• サンプルごとにすべての遺伝子を合計しないでください。これは総RNA量になり、まったく異なる値です
• 多くの遺伝子でグループ間に差がない場合、F統計量は通常低く（0.5～2.0）、高く（50以上）はなりません
• F統計量が10を超えている一方で、生物学的な状況が「有意差なし」を示唆する場合は、遺伝子レベルではなくサンプルレベルで集計している可能性があります

グループ間のログ2フォールドチェンジ — 遺伝子ごとに計算してから要約

「グループAとグループB間のログ2フォールドチェンジの中央値」を求められた場合：

• 遺伝子ごとにログ2FCを計算します。各遺伝子について、log2(mean_expr_A / mean_expr_B)を求めます
• 次に、すべての遺伝子の中央値を求めます
• 先にサンプルごとにすべての遺伝子を合計しないでください。これは総発現量の比であり、遺伝子ごとのフォールドチェンジの中央値ではありません
• ログ2FCの中央値が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で3つすべてを当てはめて順位付けする比較ラッパーを使用してください：

   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