ABテストの一連の流れをやってみる

Pythonのライブラリの中でCox比例ハザードモデルを扱う「lifelne」というライブラリには疫学データが多くあります。
で、この疫学データというのが非常に大事で、疫学はデータの集め方や原因と結果を繋げる学問なので一般的なトイデータのような「症例対象研究(簡単に言えば結果論)」だけでなく原因から結果を見る「介入研究」や「コホート研究(観察)」のデータがあります。
ここでABテストは介入データでランダム化されている必要があります。

そこで今回はその中で介入研究である白血病のデータセット「leukemia」を使ってABテストをやってみます。

ライブラリのインポート

from scipy import stats
import pandas as pd

データセットの読み込み

df = pd.read_csv("leukemia.csv")
df.head()

ここで変数は以下のようになっています(朧気なので間違っているかもしれません)。

変数名	説明
t	生存時間・打ち切りインジゲータ
status	効果(治ったか)
sex	性別
logWBC	白血球数の対数
Rx	投与

投与と効果のABテスト

ここでは両者質的変数であるためクロス集計表を使用します。

df_crs = pd.crosstab([df["sex"], df["Rx"]], df["status"])
df_crs

このクロス集計表から性別ごとのリスク比とカイ二乗検定を行います。

リスク比

crs = df_crs.values
risk = []
risk_ratio = []
for i in range(len(crs)):
    risk.append(crs[i][1]/sum(crs[i]))
    if i % 2 == 1:
        risk_ratio.append(risk[i]/risk[i-1])
    else:
        risk_ratio.append("")
rr = ["", "", "", risk_ratio[3]/risk_ratio[1]]
df_risk = pd.DataFrame(risk)
df_rr = pd.DataFrame(risk_ratio)
df_rr = pd.concat([df_risk, df_rr], axis=1)
df_rr_s = pd.DataFrame(rr)
df_rr = pd.concat([df_rr, df_rr_s], axis=1)
df_rr.columns = ["risk", "risk_ratio", "sex"]
df_rr.index = df_crs.index
df_rr

この結果から投与していない群とした群で2倍違うことから効果がある事が分かります。
また、性別は1に近いためほぼ関係ないことが分かります。

カイ二乗検定

ここまでリスク比を出しましたが、この結果が有意かどうかを算出します。

x2, p, freedeg, expected = stats.chi2_contingency(df_crs)
print("x2 = %f"%(x2))
print("p  = %f"%(p))
df_exp = pd.DataFrame(expected)
df_exp.index = df_crs.index
df_exp.columns = df_crs.columns
df_exp

x2 = 16.876364
p  = 0.000749

サンプル数が少ない中で有意確率が低く非常に有意な結果が得られました。

この結果から投与と効果に因果関係がある事が推測できます。

投与と白血球数のABテスト

次に投与群とプラセボ群で白血球数に違いがあるかを見ていきます。

Rx_1 = df[df["Rx"]==1]["logWBC"]
Rx_0 = df[df["Rx"]==0]["logWBC"]
t, p = stats.ttest_ind(Rx_1, Rx_0)
print("t-val = %f"%(t))
print("p-val = %f"%(p))
print("mean  = %f"%(Rx_1.mean()-Rx_0.mean()))

t-val = 2.168723
p-val = 0.036107
mean  = 0.588095

この結果から有意水準を5%とした場合、「投与群とプラセボ群で白血球数が同じ」という帰無仮説は棄却できることが分かります。

性別を加味する

次に性別が交絡になっていないかを見てみます。

R1s1 = df[(df["Rx"]==1) & (df["sex"]==1)]["logWBC"]
R0s0 = df[(df["Rx"]==0) & (df["sex"]==0)]["logWBC"]
R1s0 = df[(df["Rx"]==1) & (df["sex"]==0)]["logWBC"]
R0s1 = df[(df["Rx"]==0) & (df["sex"]==1)]["logWBC"]
t1, p1 = stats.ttest_ind(R1s1, R0s1)
t2, p2 = stats.ttest_ind(R1s0, R0s0)
t3, p3 = stats.ttest_ind(R1s1, R1s0)
t4, p4 = stats.ttest_ind(R0s0, R0s1)
print("Rx sex = 1")
print("t-val = %f"%(t1))
print("p-val = %f"%(p1))
print("mean  = %f"%(R1s1.mean()-R0s1.mean()))
print()
print("Rx sex = 0")
print("t-val = %f"%(t2))
print("p-val = %f"%(p2))
print("mean  = %f"%(R1s0.mean()-R0s0.mean()))
print()
print("sex Rx = 1")
print("t-val = %f"%(t1))
print("p-val = %f"%(p1))
print("mean  = %f"%(R1s1.mean()-R1s0.mean()))
print()
print("sex Rx = 0")
print("t-val = %f"%(t2))
print("p-val = %f"%(p2))
print("mean  = %f"%(R0s1.mean()-R0s0.mean()))

Rx sex = 1
t-val = 2.963402
p-val = 0.008321
mean  = 1.271000

Rx sex = 0
t-val = -0.111589
p-val = 0.912261
mean  = -0.032727

sex Rx = 1
t-val = 2.183527
p-val = 0.041743
mean  = 0.850909

sex Rx = 0
t-val = 1.368048
p-val = 0.187260
mean  = 0.452818

基本的にはp値は低いですが2つ目の性別を0にして投与したかどうかだと帰無仮説が棄却できないp値になっています。性別が0の性質を調べる必要があります。

データの分布

これら性別を加味した時にデータの分布がどのようになっているか確認します。

import matplotlib.pyplot as plt
plt.boxplot([R1s0, R0s0, R1s1, R0s1], labels=["R1s0", "R0s0", "R1s1", "R0s1"], positions=[1, 2, 3, 4])
plt.scatter([1, 2, 3, 4], [R1s0.mean(), R0s0.mean(), R1s1.mean(), R0s1.mean()], color="#000000", marker="x")
plt.show()
dcb = pd.concat([pd.DataFrame(R1s0.describe()), pd.DataFrame(R0s0.describe()), pd.DataFrame(R1s1.describe()), pd.DataFrame(R0s1.describe())], axis=1)
dcb.columns = ["R1s0", "R0s0", "R1s1", "R0s1"]
dcb

このp値が1に近いデータは四分位と平均値が近いことが分かります。また、外れ値の影響もあってか標準偏差も近くなった結果性別が0の場合は投与に白血球数の違いが見られなかったと考えられます。

実際にt検定でも関連する信頼区間を見てみましょう。

i = 0
for data in [R1s0, R0s0, R1s1, R0s1]:
    trange = stats.t(loc=data.mean(), scale=np.sqrt(data.var()/len(data.values)), df=len(data.values)-1)
    low, high = trange.interval(0.95)
    plt.plot([i, i], [low, high], marker="_", color="#000000")
    plt.scatter([i], [data.mean()], marker="x", color="#000000")
    i = i + 1
plt.xticks([0, 1, 2, 3], ["R1s0", "R0s0", "R1s1", "R0s1"])
plt.show()

こうやってみ見てみるとs0(左2つ)の時においては投与の有無に限らず分布と平均値が類似(近似)していることが分かります。
なのでsexが0ではlogWBCに影響が出ず、sexが1の時にlogWBCに影響が出ていることが分かります。
そのためこの薬品はsexが1の人に効果がある事が考えられます。

まとめ

ABテストは単純に比や差を求めるだけでなくp値を求める事でより説得力の増す結果になります。
また、交絡も加味するとより詳細な結果になり、それが結果として有意ではないことも有り得るということもありますし、どこかだけ傾向が違う(有意、有意でない)となったらそれを分析する事で因果関係がよりはっきりします。