効果測定関数を改良した

以前にいろんな<a href="https://qiita.com/murasaki1994/items/7b116747a022aa9805f0">関数集の記事</a>を書いたんですが、その中で「包括的な効果測定」の関数について改良したのでリスク差の測定をする際に是非使ってください。
改良した内容は平均値とボックスプロットだけでなく信頼区間を計測して図示するようになっています。
関数に入れるデータは

- df・・・データフレーム
- columns・・・質的変数
- y_name・・・結果系変数(量的変数)
- auto・・・分析結果の可視化と数値出力
- per・・・信頼区間

# 関数
```Python3
import scipy.stats as stats
import matplotlib.pyplot as plt
import numpy as np
def effecttest3(df, columns, y_name, auto=False, per=0.95):
    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]
            if auto:
                low1, high1 = stats.t.interval(per, len(df_tmp)-1, 
                                               loc=df_tmp[y_name].mean(), 
                                               scale=np.sqrt(df_tmp[y_name].var()/len(df[y_name].values)))
                plt.plot([x, x], [low1, high1], color="#000000", marker="_")
                print(col+"_"+str(val), ":", low1, "<=", "t", "<=", high1)
            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)
    if auto:
        for i in range(len(dtr)):
            for j in range(len(dtr[i])):
                for k in range(j, len(dtr[i])):
                    if j != k:
                        f, p = stats.bartlett(dtr[i][0], dtr[i][1])
                        if (2 * p) <= 0.05:
                            t, p = stats.ttest_ind(dtr[i][j], dtr[i][k], equal_var=False)
                        else:
                            t, p = stats.ttest_ind(dtr[i][j], dtr[i][k], equal_var=True)
                        print(lab[i][j], lab[i][k])
                        print("t   = %f"%(t))
                        print("p   = %f"%(p))
                        print("val = %f"%(ave[i][j]-ave[i][k]))
                        print()
        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()
    return ave, dtr, pos, lab
```

# 利用例
白血病データセットにおいて白血球数について原因と結果の変数と交絡因子を質的変数として入力して結果を出力します。
```Python3
import pandas as pd
df = pd.read_csv("leukemia.csv")
df.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/67d3dd3b-aa26-7baa-c9c9-3dc3919e44d2.png)

このデータを関数に入れます。
```Python3
ave, dtr, pos ,lab = effecttest3(df, ["sex", "Rx", "status"], "logWBC", auto=True, per=0.95)
```
```
sex_0 : 2.619994321512767 <= t <= 3.050914769396324
sex_1 : 2.6667672879020703 <= t <= 3.4022327120979297
Rx_0 : 2.3871019745403443 <= t <= 2.8852789778406085
Rx_1 : 2.9113372789768297 <= t <= 3.5372341495945987
status_0 : 2.0630662436551264 <= t <= 2.42860042301154
status_1 : 2.9201498667301653 <= t <= 3.487850133269835
sex_0 sex_1
t   = -0.681453
p   = 0.500781
val = -0.199045

Rx_0 Rx_1
t   = -2.168723
p   = 0.036107
val = -0.588095

status_0 status_1
t   = -3.436990
p   = 0.001386
val = -0.958167
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2782889/9e0a304d-a38c-6733-fe66-ede20b44c68e.png)

例えばこの結果だとsex＝1は平均値の外れ値が大きく95%信頼区間の下限が中央値になっていることがわかる。また、性別のp値は0.5と高いが信頼区間がほとんど被っていることが理由として考えられる。そのため性別による差は考えにくい。また、Rx(投与)について着目すると投与群は優位に白血球数が多いことが分かる。また、結果から見たとき、statusは信頼区間が一切被っていなく平均値が同じ値になる確率も0.1%となっている。

といった感じに利用できます。
材質系とかで使えるらしいのでぜひご活用ください。