Metadata-Version: 2.4
Name: affinityfit
Version: 0.1.0
Summary: A library for saturation curve fitting and estimation reliability diagnostics
Keywords: curve-fitting,binding-affinity,dose-response,dissociation-constant,langmuir-isotherm,michaelis-menten,hill-equation,confidence-interval,nonlinear-regression
Author: matsudan
Author-email: matsudan <daaamatsun@gmail.com>
License-Expression: MIT
License-File: LICENSE
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Typing :: Typed
Requires-Dist: numpy>=2.5.1
Requires-Dist: scipy>=1.18.0
Requires-Python: >=3.12
Project-URL: Homepage, https://github.com/matsudan/affinityfit
Description-Content-Type: text/markdown

# affinityfit

滴定データ（濃度とシグナルの 2 列）から **Kd** などを求め、その推定値を信用してよいかを
診断する Python ライブラリ。

観測量が結合率と線形の関係にあれば、測定手段を問わず同じ数式（飽和曲線）が使えます。
NMR のピーク強度、SPR の定常状態応答、酵素反応の初速度がその例です。

```
signal = baseline + Bmax * [L] / (Kd + [L])
```

複数のデータセットを同時にフィッティングし、パラメータをデータセット間で共有したり
定数に固定したりできます。測定範囲が Kd を挟んでいないために単独では決定不能な曲線
でも、共有によって推定が成立します。

## Usage

Python >= 3.12

```bash
uv add affinityfit
```

```python
from affinityfit import fit, load_csv

conc, signal = load_csv("titration.csv")
res = fit(conc, signal, unit="nM")

print(res.report())
print(res.params["kd"], res.intervals["kd"].format("nM"))
for warning in res.warnings:
    print(warning)
```

入力 CSV は 1 列目が濃度、2 列目がシグナル。ヘッダ行・コメント行・空行は読み飛ばします。
同じ濃度の行が複数あればレプリケートとしてそのまま使われます。

欠測（空欄）、`nan`、`inf` は行番号つきのエラーになります。空欄の行を黙って捨てると
点数が変わり、自由度と診断が静かに変わってしまうためです。

```csv
concentration_nM,signal
0,0.0249
1,0.0963
3,0.2493
10,0.5090
30,0.7657
100,0.9249
300,0.9758
1000,1.0035
```

`report()` の出力:

```
model    : langmuir  (1:1 binding: signal = baseline + Bmax * L / (Kd + L))
interval : profile
Kd       = 10.1 +/- 0.8 nM
Bmax     = 0.995 +/- 0.017
baseline = 0.018 +/- 0.014
R^2      = 0.9998   (n = 8)
AICc     = -71.12   (AIC = -77.12)

診断チェック: 問題は検出されませんでした。
```

固定側（受容体・レクチンなど）の濃度がわかっている場合は渡してください。
リガンド枯渇のチェックが有効になります。

```python
res = fit(conc, signal, receptor_conc=1.0, unit="nM")
```

モデルを変えるには `model=` を渡します。

```python
from affinityfit import hill

res = fit(conc, signal, model=hill)  # 協同性を検定する
print(res.intervals["n"].contains(1.0))
```

複数データセットを同時にフィッティングし、パラメータを共有・固定する場合:

```python
from affinityfit import Dataset, fit_global

res = fit_global(
    [Dataset("oxidized", conc, sig_ox), Dataset("reduced", conc, sig_red)],
    shared=["bmax"],  # 全データセットで 1 つの値を推定
    fixed={"baseline": 0.0},  # 定数に固定して推定しない
    unit="mM",
)
print(res.report())

# 1 データセット分を取り出す。そのデータセットの警告・注記も一緒に付いてくる。
sub = res.result_for("oxidized")
print(sub.report())
x, y = sub.curve()
```

`warnings` は問題、`notes` は「共有したおかげで推定できている」のような補足です。
`GlobalFitResult.warnings` / `.notes` は `[データセット名]` の接頭辞つきで全件を返し、
`result_for()` はそのデータセット分（およびフィッティング全体に関わるもの）を接頭辞
なしで返します。

## Weighting

既定では全点を等価値に扱います。これは「誤差の大きさがどこでも同じ」という仮定で、
SPR の RU のような系では妥当ですが、蛍光・発光・吸光では誤差が信号に比例するため
成り立ちません。飽和曲線は信号が baseline から baseline+Bmax まで動くので、飽和側の
点の絶対誤差が低濃度側の数十倍になります。

```python
res = fit(conc, signal, sigma=0.01 + 0.10 * signal)  # 点ごとの標準偏差
```

`sigma` は相対値として扱われ、全体のスケールは残差から推定されます。定数倍しても
結果は変わりません。グローバルフィットでは `Dataset(..., sigma=...)` でデータセット
ごとに渡せ、データセット間の相対的な重み（片方が 10 倍ノイジー等）も反映されます。

sigma は**全データセットに与えるか、どこにも与えないか**のどちらかです。一部にだけ
与えると、sigma のないデータセットの重みが暗黙に 1 になり、両者の重み比が sigma の
絶対スケール（分数か % か ppm か）で決まってしまうためです。混在はエラーになります。

重みを渡さずに不均一分散のデータをフィッティングすると、精度が落ちるだけでなく
**信頼区間が狭く出ます**。Kd=10 に対し 30% の比例誤差、Kd が測定範囲上端にある設計
での実測です。

| | Kd の四分位範囲 | 95%CI の被覆率 |
|---|---|---|
| 重みなし | 19.7 | 79% |
| sigma 指定 | 11.3 | 94% |

残差の大きさがフィッティング値とともに増えている場合は診断が警告します（誤報率 2%、
比例誤差の検出率 56〜71%）。

## Confidence Intervals

`ci=` で 3 つの方法を選べます。

| 方法 | 内容 |
|---|---|
| `profile`（既定） | Kd を固定して他を再フィットし、残差平方和が有意に増える点を境界とする（F 検定）。区間は非対称になりうる。片側が決まらない場合は「Kd > 1.6」という形で返す |
| `asymptotic` | 共分散行列の曲率から区間を出す。最も高速。Kd のような濃度定数は対数スケールで区間を作るので下限が負になることはないが、曲面が閉じていなくても有限の両側区間を主張してしまう（識別不能なデータで 2 桁以上に広がった区間が出る） |
| `bootstrap` | データを再抽出して推定値の分布を直接得る。反復測定があれば反復を、なければ残差を再抽出する |

反復実験から誤差を出す慣行はこの分野で広く使われています。反復測定を渡せばその
流儀に揃えられます。

```python
res = fit(conc, signal, ci="bootstrap", replicates=reps, n_boot=2000)
```

`Dataset` では反復測定だけを渡せます（`signal` は平均が使われます）。

```python
Dataset("oxidized", conc, replicates=reps)
```

`n_boot` は百分位区間に必要な下限（100）以上でなければエラーになります。下回った値を
渡すと区間が作れないためです。

報告される桁数は不確かさに合わせて丸められます。相対誤差 17% の測定から
`Kd = 4.70e-8` と 3 桁で書くのは過剰なので、`(4.7 +/- 0.8)e-08` になります。

## Model

| モデル | パラメータ | 用途 |
|---|---|---|
| `langmuir`（既定） | kd, bmax, baseline | 1:1 結合。`signal = baseline + Bmax·L/(Kd+L)` |
| `hill` | kd, bmax, baseline, n | 協同性。n の信頼区間が 1 を含むかで判定 |
| `michaelis` | km, vmax, baseline | 酵素反応速度。式は langmuir と同一だが Km は親和性ではない |

`langmuir` と `hill` の `bmax` は符号を制限していません。結合によって**減少する**観測量
（蛍光クエンチ、強度の減少）でも Kd の意味は変わらず、負の応答係数として同じ式で表せます。
`michaelis` の `vmax` は反応速度なので非負に留めているため、減少するデータには
`langmuir` か `hill` を使ってください。噛み合わないモデルを指定した場合は、それらしい
数値を返さずに警告します。

`Km = (k_off + k_cat)/k_on` であり、`k_cat` が無視できるときだけ Kd に一致します。

Kd や Km は対数スケールで最適化されます。線形に扱うとフィッティング結果が濃度の単位に
依存し、M 単位のまま pM 親和性（Kd = 1e-12 M）を扱うと相対誤差が 6800% に達します。
現在は fM から 100 M まで 18 桁にわたって相対誤差 1e-6 以内で回収し、同じ実験を
M / nM / pM のどの単位で書いても同じ答えになります。
モデルの追加は `models.py` に `Model` を 1 つ定義するだけで、フィッティングと診断の
コードには手を入れません。

## Model Selection

モデルやパラメータ共有の比較には `aicc`（小標本補正つき赤池情報量規準）を使ってください。
`aic` も参照用に持っていますが、補正なしの AIC は漸近的にしか妥当でなく、滴定実験の
規模（点数 6〜15、パラメータ 3〜6 で n/k が 2〜5）ではパラメータの多い側に偏ります。

実測: 振幅が真に共通なデータ（12 点、共有で 3 パラメータ / 非共有で 4）40 通りで、
補正なし AIC が共有を選んだのは 32 回、AICc は 38 回でした。10 点・5 対 6 パラメータの
構成では 12 通り中 7 通りで結論が食い違い、いずれも AIC が過剰パラメータ側を選びました。

```python
shared.aicc < free.aicc
fit_global(ds, model=hill).aicc < fit_global(ds, model=langmuir).aicc
```

標本が小さすぎて補正が定義できない場合（n − k − 1 ≤ 0）、`aicc` は無限大になります。
`report()` は AICc を先に、参照用の AIC を括弧内に表示します。

## Diagnostics

R² が 0.99 を超えていても Kd が決まっていないことは珍しくありません。以下を自動で警告します。
Kd / Km のようにモデルごとに呼び名が変わる量は、そのモデルの表記で表示されます。

診断は2種類あります。**フィッティング自体が健全か**と、**測定点の配置が適切か**です。

| 条件 | 意味 |
|---|---|
| 振幅がほぼ 0 に潰れた | フィッティングは実質的に水平線。モデルがデータの形を表現できていない |
| R² < 0.5 | モデルがデータの傾向を捉えていない。Kd の値に意味はない |
| R² < 0.9 | 飽和曲線としては低め。ノイズかモデル選択の問題 |
| 残差の符号が系統的 | 決定係数が高くてもモデルの形が機構に合っていない |
| パラメータが境界に張り付いた | その値は推定結果ではなく制約の産物。報告できない |
| 最高濃度 < 3 × Kd | 飽和未達。Kd と Bmax が数学的に分離できず、「Kd > 最高濃度」としか言えない |
| 最高濃度 < 10 × Kd | Bmax の推定が不安定で、Kd の信頼区間も広がる |
| データ点 < 2 × 推定パラメータ数 | 情報が不足。信頼区間は参考値（固定したパラメータは数えない） |
| Kd 近傍（Kd/3〜3Kd）に点が 0〜1 個 | 曲線の変曲点が決まらない。ここに点を足すのが最も効く |
| 最低濃度 > Kd | 全点が飽和側。Kd は外挿で決まっており、有効数字を増やして報告できない |
| Kd/10 以下の点がない | baseline が曲線と一緒に推定され、Bmax がずれうる |
| 受容体濃度 > Kd/10 | 結合でリガンドが減るため Kd を過大評価。tight-binding（二次式）モデルが必要 |
| Hill 係数 n の CI が 1 を含む | 協同性を主張できない |
| 残差の大きさがフィッティング値に比例 | 誤差が不均一。`sigma` を渡さないと区間が狭く出る |

フィッティング自体が壊れている場合（R² < 0.5、振幅が潰れた）、「共有したおかげで
推定できている」という NOTE は出しません。矛盾した助言になるためです。
| 自由度が 0（点数 ≤ パラメータ数） | 曲線は定義上すべての点を通る。信頼区間は算出できない |
| ヤコビ行列がランク落ち | パラメータが区別できず値が一意に決まらない。濃度の種類を増やす |
| 信頼区間の片側が決定できない | 点推定値を報告してはいけない。共有か測定範囲の拡張が必要 |

グローバルフィットでは、共有・固定によって前提が消えた指摘は抑制されます。たとえば
`bmax` を共有していれば「Kd と Bmax を分離できない」は当てはまらないので、代わりに
「共有しているため推定できている」という NOTE が出ます。

`examples/titration_unsaturated.csv` は R² = 0.9936 と当てはまりが良く見えても、最高濃度が
Kd の 1/3 しかないため上限が決まらず `Kd > 76 nM` としか言えない例です。

信頼区間が算出できない場合は数値を作りません。自由度が 0 のとき、ヤコビ行列が
ランク落ちしているとき、プロファイル探索の内側の再フィットが収束しなかったとき、
ブートストラップの成功回数が足りないときは、区間を「決定不能」として返し、理由を
警告で説明します。点推定値そのものは返します。

ブートストラップは収束しなかった再抽出の回数を数えて報告します。収束する再抽出は
フィッティングしやすいものに偏るため、黙って捨てると区間が系統的に狭くなります。

残差が厳密に 0 のフィッティング（合成データなど）では区間幅も 0 になりますが、裸の
数値ではなく `10.000 +/- 0 (no residual scatter)` と理由つきで表示します。

残差の検定にはラン検定（Wald–Wolfowitz）と隣接残差の自己相関の両方を使い、どちらかが
有意なら警告します。15 点の滴定を 300 通りのノイズで試した実測で、誤ったモデル
（協同性データを 1:1 でフィッティングする、R² は 0.94〜0.98）の検出率は 100%、
正しいモデルでの誤検出率は 0.3〜1.7% でした。

### 統計量

警告の判定は固定の有意水準（例: 残差検定は片側 2.5%、不均一分散は片側 1%）で行われて
います。これは 1 件のフィッティングには適切ですが、複数のデータセットや複数の論文図に
わたる検定を積み重ねる場合には、その有意水準は検定の数だけ増えた偽陽性率を勝手に
受け入れていることになります。何本の検定になるか、それをどう補正すべきかはこちら側では
決められないので、判定の元になった統計量と p 値をそのまま返します。

```python
res = fit(conc, signal)
for s in res.statistics:
    print(s.name, s.statistic, s.p_value)
```

複数データセットの場合は `GlobalFitResult.statistics_per["データセット名"]` で
データセットごとに取得できます。集めた p 値に対して、`scipy.stats.false_discovery_control`
や `statsmodels.stats.multitest.multipletests` で Bonferroni や Benjamini–Hochberg
補正をかけられます。

```python
from scipy.stats import false_discovery_control

p_values = [
    s.p_value
    for name, stats in res.statistics_per.items()
    for s in stats
    if s.name == "heteroscedasticity" and s.p_value is not None
]
q_values = false_discovery_control(p_values, method="bh")
```

`residual_autocorrelation` は固定閾値（0.3）による判定で、p 値を持ちません
（`p_value` は `None`）。

## Plotting

図を返す API はありません。論文用の図はフォント・色・パネル構成を自分で制御する必要が
あるので、`affinityfit` は描画に必要な数値だけを返します。

| メソッド | 返すもの |
|---|---|
| `res.curve(conc_min, conc_max, n)` | 対数間隔のフィッティング曲線 `(x, y)` |
| `res.predict(conc)` | 任意の濃度におけるフィッティング値 |
| `res.residuals(conc, signal)` | 残差（観測値 − フィッティング値） |

```python
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(*res.curve(), "-")
ax.plot(conc, signal, "o")
ax.set_xscale("log")
```

濃度軸は対数にしてください。線形軸だと低濃度側が潰れ、Kd 近傍の曲率と「測定点が
飽和領域に偏っている」といった設計上の問題が見えなくなります。

`examples/plot_fit.py` に、測定点・フィッティング曲線・残差パネルを 1 枚に描くサンプル
があります。配布物には含まれません。

```bash
uv run python examples/plot_fit.py examples/titration_good.csv --unit nM
```

`examples/titration_good.csv` にこのサンプルを実行した結果:

![langmuir fit of examples/titration_good.csv](examples/fit_good.png)

## Development

```bash
uv run pytest -q
uv run ruff check .
uv run ruff format --check .
uv run ty check
```

テストは既知のパラメータ（Kd は 1e-16 〜 1e2 の 18 桁、Hill 係数は 0.4 〜 3）を回収
できるか、各警告が意図した条件でのみ発火するか、共有パラメータが識別不能なデータセット
を救うかを検証しています。信頼区間については、良条件のデータで 3 手法が一致すること、
プロファイル区間の端が F 検定の閾値に厳密に一致すること、識別不能なデータで片側限界に
なることを確認しています。`examples/` のサンプルもテストから実際に呼んでいるので、
API を変えれば落ちます。

## Roadmap

未実装のもの。上の各節に書かれている機能はすべて動きます。

**モデル**

- 時間経過モデル `A(1 - exp(-k·t))`。プラトーに達していないと振幅と速度定数が
  分離できないため、その判定つきで
- 競合阻害モデル（阻害剤存在下の見かけの Kd から真の Ki を出す）
- tight-binding（二次式）モデル。受容体濃度が Kd と同程度の場合に必要
- 山型のプロファイル（金属濃度やイオン強度の最適値、pH 最適）。単調でないため
  飽和曲線とは別の族になる

**誤差モデル**

- `sigma` を自前で用意できない場合の分散モデル（`proportional` / `poisson`）。
  当てはめ値から重みを作るため反復再重み付けが必要で、当てはめ値が 0 の点で
  重みが発散する。σ の下限を含めた定式化が前提
- 反復測定から `sigma` を推定する。反復数が少ないと標準偏差自体がばらつくので、
  そのばらつきを織り込む必要がある

**入出力**

- 複数データセットの結果表を CSV / JSON で出力するヘルパ
- 装置固有のエクスポート形式（SPR のセンサーグラムなど）のインポータ
- LOD / LOQ の算出（IUPAC / CLSI 準拠）
