変数種別データ分析の方針

CSVだったりデータベースだったりエクセルだったりデータセットは色々ありますが、何かを目的に分析をする上で変数の種類は非常に重要になります。
具体的に変数の種類は以下のようになります
| 説明変数 | 目的変数 | 分析方法 |
|:-:|:-:|:-:|
| 質的変数  | 質的変数  | クロス集計表  |
|  質的変数 |  量的変数 | 平均値の差や要約統計量  |
|  量的変数 | 質的変数  | 分類  |
| 量的変数  | 量的変数  | 回帰  |
|  量的変数と時間 | 質的変数  | ハザード分析  |

では実際に変数別の分析の一例を考えてみましょう
# 説明変数が質的変数で目的変数が質的変数
ここでは喫煙の有無と肺癌の有無についてです。データセットは『基礎から学ぶ 楽しい疫学 第4版』にある表から作成しました。
## ライブラリのインポート
今回はクロス集計表を作るだけでなく、有意かどうかも見るためscipyも使います。
```Python3
import pandas as pd
import scipy.stats as stats
```
## データの読み込み
```Python3
df = pd.read_csv("LungCanser.csv")
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/dc56fc8d-7960-b93b-b060-ebe08b2dee61.png)
ここで1は有で0は無です。
## クロス集計表
```Python3
crs = pd.crosstab(df["Smoke"], df["LungCanser"])
crs
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/4a1bd53d-b715-244b-25e7-19b018ea278c.png)

## リスク比の計算
```Python3
val = crs.values
risk = []
risk.append([val[0][1]/sum(val[0]), val[1][1]/sum(val[1])])
df_risk = pd.DataFrame(risk).T
rr = []
rr.append([risk[0][1]/risk[0][0]])
df_rr = pd.DataFrame(rr)
df_risk = pd.concat([df_risk, df_rr], axis=1)
df_risk.columns = ["risk", "risk_ratio"]
crs = pd.concat([crs, df_risk], axis=1)
crs
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/a18a5b15-c368-c3b3-18f9-0310b3f2d559.png)
この結果から喫煙者は非喫煙者と比べて2.3倍肺癌になりやすいと解釈できます。
## 有意かどうか
p値を実際に算出して帰無仮説である理論値と実測値が同じになる確率を計算して有意な結果なのかを考えます。
```Python3
x2, p, df, expected = stats.chi2_contingency(crs)
print("X2\t=%f\np_val\t=%f\ndf\t=%d"%(x2, p, df))
```
```
X2      =30.420000
p_val	=0.000000
df	    =1
```
この結果からカイ二乗値は30.42で自由度が1、そしてp値は0に限りなく近いため有意であることが分かります。
## 理論値
この集計表において理論値はどうなっていたかを確認します。
```Python3
df_exp = pd.DataFrame(expected)
df_exp.columns = crs.columns
df_exp.index = crs.index
df_exp
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/c81bf0dc-01f9-f49c-67bb-99b9f632ca23.png)

当たり前ですが全て50になります。
## 実測値と理論値の差
この差の大きさや正負によって実際理論値とどのように違うか考察します。
```Python3
df_diff = pd.DataFrame(crs-expected)
df_diff
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/fffb2ee3-c248-84ea-b7a7-b50e37b07515.png)
この結果から喫煙しない人は本来肺癌になりにくく、喫煙者は肺がんになりやすい事が分かります。

# 説明変数が質的変数で目的変数が量的変数
説明変数が質的変数という事は場合分けして目的変数の値を考察する事が出来ます。そこで、効果検証を使ってみようと思います。
また、しっかりと解析するにはしっかりと先述したp値も使います。
## ライブラリのインポート
```Python3
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import scipy.stats as stats
```
## 因果推論の材料を作る関数
過去に因果推論を効率的にやるプログラムを作りましたが、複数回やるなら関数化した方が効率的なので関数にします。

https://qiita.com/murasaki1994/items/1b3bc37ee01203363edc

```Python3
def effecttest(df, columns, y_name):
    ave = []
    dtr = []
    pos = []
    lab = []
    x = 0
    for col in columns:
        values = list(set(df[col].values))
        tmp_ave = []
        tmp_dtr = []
        tmp_pos = []
        tmp_lab = []
        for val in values:
            df_tmp = df[df[col]==val]
            tmp_ave.append(df_tmp[y_name].mean())
            tmp_dtr.append(df_tmp[y_name])
            tmp_pos.append(x)
            tmp_lab.append(col+"_"+str(val))
            x = x + 1
        ave.append(tmp_ave)
        dtr.append(tmp_dtr)
        pos.append(tmp_pos)
        lab.append(tmp_lab)
    return ave, dtr, pos , lab
```
ではここから効果検証を様々とやってみましょう。
## データの読み込み
lifelineというライブラリにあるrossiというデータからprio(逮捕回数)以外を除いたデータセットをCSVで読み込みます。
```Python3
df = pd.read_csv("rossi2.csv")
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/54ed2be5-fd29-b7d6-fb8d-01ebcb66206b.png)

## 関数の利用と図示
```Python3
columns = []
for col in df.columns:
    if col != "prio":
        columns.append(col)
ave, dtr, pos, lab = effecttest(df, columns, "prio")
for i in range(len(dtr)):
    plt.boxplot(dtr[i], positions=pos[i], labels=lab[i])
    plt.plot(pos[i], ave[i], marker="x")
plt.xticks(rotation=90)
plt.show()
```
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/712ee438-39af-19f9-bbc5-e22b7f6a9606.png)
この結果を見るとwexp(労働経験の有無)がとarrest(逮捕)の有無が大きく影響している事が分かります。
## 要約統計量
箱ひげ図は要約統計量を図示しているだけなので実際の数値を見てみます。
```Python3
df_dcb = pd.DataFrame(dtr[0][0].describe())
for i in range(len(dtr)):
    for j in range(len(dtr[i])):
        if i != 0 or j!=0:
            df_dcb = pd.concat([df_dcb, pd.DataFrame(dtr[i][j].describe())], axis=1)
df_dcb.columns = np.array(lab).flatten()
df_dcb
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/3f1612b4-7bb6-489f-ae48-3d5400f3f0fd.png)

## 仮説検証
ここでは厳密にやりたいのでF検定を行い分散を等しいとするかどうかも考慮してt検定を行います。
```Python3
for i in range(len(dtr)):
    print(columns[i])
    f, p = stats.bartlett(dtr[i][0], dtr[i][1])
    if (2 * p) <= 0.05:
        t, p = stats.ttest_ind(dtr[i][0], dtr[i][1], equal_var=False)
    else:
        t, p = stats.ttest_ind(dtr[i][0], dtr[i][1], equal_var=True)
    print("t   = %f"%(t))
    print("p   = %f"%(p))
    print("val = %f"%(ave[i][0]-ave[i][1]))
    print()
```
```
arrest
t   = -2.931893
p   = 0.003877
val = -1.070672

fin
t   = 0.016594
p   = 0.986768
val = 0.004630

race
t   = 1.871865
p   = 0.061904
val = 0.792702

wexp
t   = 5.252460
p   = 0.000000
val = 1.522070

mar
t   = 0.614361
p   = 0.539302
val = 0.261114

paro
t   = 2.479255
p   = 0.013753
val = 0.751856
```
結果から逮捕と労働が大きく関係していることが分かります。
## (補足)条件付きでやる場合
先ほど因果推論の材料を作りましたので、これを利用して条件を付けて同じことをやってみようと思います。
条件としては逮捕の有無で分布がどう変わるかを見てみます。
```Python3
df_a0 = df[df["arrest"]==0]
df_a1 = df[df["arrest"]==1]
ave0, dtr0, pos0, lab0 = effecttest(df_a0, columns, "prio")
ave1, dtr1, pos1, lab1 = effecttest(df_a1, columns, "prio")
for i in range(len(dtr)):
    plt.boxplot(dtr0[i], positions=pos0[i], labels=lab0[i])
    plt.plot(pos0[i], ave0[i], marker="x")
plt.xticks(rotation=90)
plt.show()
for i in range(len(dtr)):
    plt.boxplot(dtr1[i], positions=pos1[i], labels=lab1[i])
    plt.plot(pos1[i], ave1[i], marker="x")
plt.xticks(rotation=90)
plt.show()
```
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/436083ee-1a98-09f5-c605-3a766173f639.png)
![Untitled-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/7164c80b-9aa0-9524-2916-48ab9b9f61c1.png)

逮捕されている場合だとmar(結婚の有無)に少し差ができて労働の有無がより差が出る事が分かります。

# 説明変数が量的変数で目的変数が質的変数
この場合は分類問題が適用できます。
## ライブラリのインポート
```Python3
from sklearn.metrics import classification_report
import statsmodels.api as sm
import pandas as pd
import numpy as np
```
## データの読み込み
ここでは乳がんのデータセットを読み込みます。
```Python3
df = pd.read_csv("breast_cancer.csv")
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/2d2a79c4-1763-5888-c2dc-7ef5e64a4436.png)

ブラウザの都合上見切れていますがこのようなデータです。
## モデルの適合
ロジスティック回帰の適合させます。
```Python3
y = df["y"]
x = sm.add_constant(df.drop("y", axis=1))
logit = sm.Logit(y, x).fit_regularized()
logit.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/3fe617dd-74b0-73f9-7020-d0bb6db5a7ed.png)
適合した結果はこうなりました。

## 精度の確認
ここではしきい値を0.5にします。
```Python3
y_pred_proba = logit.predict(x)
y_pred = np.where(y_pred_proba >= 0.5, 1, 0)
print(classification_report(y, y_pred))
```
```
              precision    recall  f1-score   support

         0.0       0.99      0.98      0.98       212
         1.0       0.99      0.99      0.99       357

    accuracy                           0.99       569
   macro avg       0.99      0.99      0.99       569
weighted avg       0.99      0.99      0.99       569

```

# 説明変数が量的変数で目的変数が量的変数
どちらとも量的変数である場合は回帰を使うことができます。
## ライブラリのインポート
```Python3
from sklearn.metrics import mean_squared_error as mse
from sklearn.metrics import mean_absolute_error as mae
import statsmodels.api as sm
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
```
## データの読み込み
ここでは糖尿病のデータセットを読み込みます。
```Python3
df = pd.read_csv("diabetes.csv")
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/6561b0cd-e19e-4846-09f4-980af4c8efcb.png)


## モデルの適合
モデルには最小二乗法をここでは用います。
```Python3
y = df["y"]
x = sm.add_constant(df.drop("y", axis=1))
ols = sm.OLS(y, x).fit()
ols.summary()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/9b9b5d09-c001-47e1-bd76-de8626529e04.png)

## 誤差を確認
決定係数は0.518であることが分かったが実際の誤差がどの程度かも検証する。
```Python3
y_pred = ols.predict(x)
print("MAE  = %f"%mae(y, y_pred))
print("RMSE = %f"%np.sqrt(mse(y, y_pred)))
```
```
MAE  = 43.277452
RMSE = 53.476129
```
## 間違え方の確認
ここまで精度が出ましたが、ではどういう風に間違えたか、異常値や傾向が無いか確認します。
```Python3
plt.scatter(y_pred, y)
plt.ylabel("real")
plt.xlabel("pred")
plt.title("corr = %.2f"%(np.corrcoef(y, y_pred)[0][1]))
plt.show()
```
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/182e51b6-0ea2-0e2b-6ed2-d6a14db1c5ae.png)

特に傾向も無ければ異常な値も無いと思います。

# 説明変数が量的で時間もあり目的変数が質的変数の場合
この場合イベント発生におけるハザード比を算出します。
## ライブラリのインポート
```Python3
from lifelines import CoxPHFitter
import matplotlib.pyplot as plt
import pandas as pd
```
## データの読み込み
ここでは白血病のデータセットを使います。
```Python3
df = pd.read_csv("leukemia.csv")
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/e97aee3b-cbaa-0910-59ac-e16c19b1d953.png)

## 要約統計量
後にハザード比の影響を確認する上で使用する値を確認するために出力します。
```Python3
df.describe()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/cb718973-def1-3c39-b338-8a0625cb5a81.png)

## モデルの適合
ハザード比の算出にCox比例ハザードモデルを使用します。
```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/31267225-f361-c1c3-e856-73394a03bd1f.png)
この結果からlogWBCとRxが値が大きくなるほどハザード比に影響しやすいことが分かります。
## 変数別ハザード比のプロット
```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/3b789d71-dfef-5720-c305-d336fd23b6fa.png)
![Untitled-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/b8385081-86fd-77df-9530-be1501813ba9.png)
![Untitled.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/a6ab700b-d40b-2d9c-63c6-a6ec677ee97f.png)
この結果からRx(投与の有無)とlogWBC(白血球数のlog)が影響力が強いことが分かります。

# まとめ
データの種類によって分析方法を考えましょう。

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