疫学の教科書に載っていた数理モデル

前に因果関係で疫学について書きましたが

https://qiita.com/murasaki1994/items/267a200ad9075fbcc613

疫学は原因と結果を繋げる学問で、そこに多変量解析が載っているので紹介してみたいと思います。
医学書院にある『基礎から学ぶ楽しい疫学第4版』の第10章では多変量解析について取り扱われています。
![108378.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/d59b1001-9b22-29de-ea3f-762961f0f2db.jpeg)

第10章は「疫学に必要な統計」で、その中の多変量解析の中身を今回取り扱います。
使用する数理モデルは以下の4つです。
- 重回帰分析(OLS)
数量データをそのまま従属変数に用いたもの
- ポアソン回帰
層化を用いた多変量解析に利用される
- ロジスティック回帰
2値データ(疾患の発生/非発生, 症例/対照, など)を従属変数に用いるもの
- Cox比例ハザードモデル
2群間の生存曲線を比較する手法


になります。

# コーディング
まずライブラリをインポートします。
※本当なら<a href="https://qiita.com/murasaki1994/items/b7f76390752d4423846c">VIF</a>をしっかり考慮して変数選択してください
※生データにおける係数を今回見るため標準化はしていませんが、重回帰分析とポアソン回帰は元々標準化されているデータです。

https://qiita.com/murasaki1994/items/b7f76390752d4423846c

```Python3
from lifelines import CoxPHFitter
import matplotlib.pyplot as plt
import statsmodels.api as sm
import pandas as pd
```
## 重回帰分析
### データの読み込み
```Python3
df = pd.read_csv("diabetes.csv")
df.head()
```
### 説明変数と目的変数を分ける
```Python3
x = sm.add_constant(df.drop("y", axis=1))
y = df["y"]
```
### モデル作成と諸々結果の表示
```Python3
ols = sm.OLS(y, x).fit()
ols.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/6a349863-5325-e4fd-a985-4fedebae3b26.png)
## ポアソン回帰
```Python3
poi = sm.Poisson(y, x).fit()
poi.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/ad76118a-ae19-f11d-d822-3e1800fc1844.png)

この結果から考えると決定係数では重回帰分析に分がありますが、一方で有意にな係数が多いのはポアソン回帰であることが分かります。

## ロジスティック回帰
### データの読み込み
```Python3
df = pd.read_csv("breast_cancer.csv")
df.head()
```
### 説明変数と目的変数を分ける
```Python3
x = sm.add_constant(df.drop("y", axis=1))
y = df["y"]
```
### モデル作成と諸々結果の表示
```Python3
logit = sm.Logit(y, x).fit_regularized()
logit.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/62207036-2497-9f5b-37ba-c275b2f5bda2.png)
この時z値が多いいほど分類において変数の値が大きいと陽性(1)になりやすくなります。
## Cox比例ハザード
あまり聞きなれないモデルだと思いますが、簡単に言いますと何らかの介入からイベントが発生するまでの時間や打ち切りまでの時間を含めて適合させることで時間と共にイベントがどの程度発生しやすいかを表すロジスティック回帰のような数理モデルです(多分厳密には違う)。
### データの読み込み
```Python3
df = pd.read_csv("leukemia.csv")
df.head()
```

### 適合と結果の表示
```Python3
cox = CoxPHFitter()
cox.fit(df, duration_col="t", event_col="status")
cox.print_summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/c49ba3d6-c8d3-6b74-ada3-4ca423c1326b.png)
で、ここで見るべき項目としてはp値とz値で、性別がp値で有意ではないのに対して白血球数と投与は有意であることが分かります。
これらを基に生存曲線を考察してみましょう。
### 要約統計量
要約統計量に出てきた値と平均値を基に生存曲線を出力します。
```Python3
df.describe()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/4e453267-127b-b494-3a04-e85767de6dd1.png)

### 生存曲線の出力
```Python3
cox.plot_partial_effects_on_outcome(covariates="sex", values=[0, 1], cmap='coolwarm')
cox.plot_partial_effects_on_outcome(covariates="logWBC", values=[1.45, 2.3025, 2.8, 2.930238, 3.49, 5], cmap='coolwarm')
cox.plot_partial_effects_on_outcome(covariates="Rx", values=[0, 1], cmap='coolwarm')
plt.show()
```
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/ed34eed3-e061-4517-c4f3-88b3a7ac7954.png)
![Untitled-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/794c62e5-dffe-f78d-084d-fab361dc23ab.png)
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/d18ba49a-e48f-b23d-5343-adedf20a1f67.png)
性別は1の場合が時間と共にイベント発生しやすいことが分かりますが、それほど差が無い事が分かります。実際p値も0.49で有意とは言えません。
しかしlogWBCは値が高くなるほど顕著に0になるまで早く、投与群も投与した方が早いことが分かります。

# まとめ
回帰分析も最小二乗法やポアソン回帰など色々ありますが用途に合わせて適切に使用しましょう。
2値分類においても時間があるかないかや使い道でロジスティック回帰を使うかCox比例ハザードモデルを使うか分けましょう。

# 参考文献
<a href="https://www.igaku-shoin.co.jp/book/detail/108378">基礎から学ぶ 楽しい疫学　第4版</a>