CSVだったりデータベースだったりエクセルだったりデータセットは色々ありますが、何かを目的に分析をする上で変数の種類は非常に重要になります。
具体的に変数の種類は以下のようになります
| 説明変数 | 目的変数 | 分析方法 |
|---|---|---|
| 質的変数 | 質的変数 | クロス集計表 |
| 質的変数 | 量的変数 | 平均値の差や要約統計量 |
| 量的変数 | 質的変数 | 分類 |
| 量的変数 | 量的変数 | 回帰 |
| 量的変数と時間 | 質的変数 | ハザード分析 |
では実際に変数別の分析の一例を考えてみましょう
説明変数が質的変数で目的変数が質的変数
ここでは喫煙の有無と肺癌の有無についてです。データセットは『基礎から学ぶ 楽しい疫学 第4版』にある表から作成しました。
ライブラリのインポート
今回はクロス集計表を作るだけでなく、有意かどうかも見るためscipyも使います。
import pandas as pd
import scipy.stats as stats
データの読み込み
df = pd.read_csv("LungCanser.csv")
df.head()
ここで1は有で0は無です。
クロス集計表
crs = pd.crosstab(df["Smoke"], df["LungCanser"])
crs
リスク比の計算
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
この結果から喫煙者は非喫煙者と比べて2.3倍肺癌になりやすいと解釈できます。
有意かどうか
p値を実際に算出して帰無仮説である理論値と実測値が同じになる確率を計算して有意な結果なのかを考えます。
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に限りなく近いため有意であることが分かります。
理論値
この集計表において理論値はどうなっていたかを確認します。
df_exp = pd.DataFrame(expected)
df_exp.columns = crs.columns
df_exp.index = crs.index
df_exp
当たり前ですが全て50になります。
実測値と理論値の差
この差の大きさや正負によって実際理論値とどのように違うか考察します。
df_diff = pd.DataFrame(crs-expected)
df_diff
この結果から喫煙しない人は本来肺癌になりにくく、喫煙者は肺がんになりやすい事が分かります。
説明変数が質的変数で目的変数が量的変数
説明変数が質的変数という事は場合分けして目的変数の値を考察する事が出来ます。そこで、効果検証を使ってみようと思います。
また、しっかりと解析するにはしっかりと先述したp値も使います。
ライブラリのインポート
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import scipy.stats as stats
因果推論の材料を作る関数
過去に因果推論を効率的にやるプログラムを作りましたが、複数回やるなら関数化した方が効率的なので関数にします。
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で読み込みます。
df = pd.read_csv("rossi2.csv")
df.head()
関数の利用と図示
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()
この結果を見るとwexp(労働経験の有無)がとarrest(逮捕)の有無が大きく影響している事が分かります。
要約統計量
箱ひげ図は要約統計量を図示しているだけなので実際の数値を見てみます。
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
仮説検証
ここでは厳密にやりたいのでF検定を行い分散を等しいとするかどうかも考慮してt検定を行います。
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
結果から逮捕と労働が大きく関係していることが分かります。
(補足)条件付きでやる場合
先ほど因果推論の材料を作りましたので、これを利用して条件を付けて同じことをやってみようと思います。
条件としては逮捕の有無で分布がどう変わるかを見てみます。
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()
逮捕されている場合だとmar(結婚の有無)に少し差ができて労働の有無がより差が出る事が分かります。
説明変数が量的変数で目的変数が質的変数
この場合は分類問題が適用できます。
ライブラリのインポート
from sklearn.metrics import classification_report
import statsmodels.api as sm
import pandas as pd
import numpy as np
データの読み込み
ここでは乳がんのデータセットを読み込みます。
df = pd.read_csv("breast_cancer.csv")
df.head()
ブラウザの都合上見切れていますがこのようなデータです。
モデルの適合
ロジスティック回帰の適合させます。
y = df["y"]
x = sm.add_constant(df.drop("y", axis=1))
logit = sm.Logit(y, x).fit_regularized()
logit.summary()
適合した結果はこうなりました。
精度の確認
ここではしきい値を0.5にします。
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
説明変数が量的変数で目的変数が量的変数
どちらとも量的変数である場合は回帰を使うことができます。
ライブラリのインポート
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
データの読み込み
ここでは糖尿病のデータセットを読み込みます。
df = pd.read_csv("diabetes.csv")
df.head()
モデルの適合
モデルには最小二乗法をここでは用います。
y = df["y"]
x = sm.add_constant(df.drop("y", axis=1))
ols = sm.OLS(y, x).fit()
ols.summary()
誤差を確認
決定係数は0.518であることが分かったが実際の誤差がどの程度かも検証する。
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
間違え方の確認
ここまで精度が出ましたが、ではどういう風に間違えたか、異常値や傾向が無いか確認します。
plt.scatter(y_pred, y)
plt.ylabel("real")
plt.xlabel("pred")
plt.title("corr = %.2f"%(np.corrcoef(y, y_pred)[0][1]))
plt.show()
特に傾向も無ければ異常な値も無いと思います。
説明変数が量的で時間もあり目的変数が質的変数の場合
この場合イベント発生におけるハザード比を算出します。
ライブラリのインポート
from lifelines import CoxPHFitter
import matplotlib.pyplot as plt
import pandas as pd
データの読み込み
ここでは白血病のデータセットを使います。
df = pd.read_csv("leukemia.csv")
df.head()
要約統計量
後にハザード比の影響を確認する上で使用する値を確認するために出力します。
df.describe()
モデルの適合
ハザード比の算出にCox比例ハザードモデルを使用します。
cox = CoxPHFitter()
cox.fit(df, duration_col="t",event_col="status")
cox.print_summary()
この結果からlogWBCとRxが値が大きくなるほどハザード比に影響しやすいことが分かります。
変数別ハザード比のプロット
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()
この結果からRx(投与の有無)とlogWBC(白血球数のlog)が影響力が強いことが分かります。
まとめ
データの種類によって分析方法を考えましょう。