この記事は（MATLABやRではなく）Pythonで論文用の図を作成したい人に向けたチュートリアルです。
データに標準誤差の範囲を付けてプロットしたり、統計検定をしたりします。「***」や「n.s.」などの検定有意差の表示については、オリジナルの関数を作成したので、そちらを活用してください。
必要なライブラリの読み込み
本記事では以下のライブラリを使用します。インストールされていないものがある場合は、pip installコマンドを使用してライブラリを使えるようにしてください。
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.anova import AnovaRM
from statsmodels.stats.multicomp import pairwise_tukeyhsd仮想データの生成
10人の被験者に対し、3種類の薬を1日ごとに投与して、scoreを計測することを想定したデータを生成します。薬を投与した結果のscoreは、平均がそれぞれ 0.2, -1.0, 0.5 の正規分布を取るものとします。
drug_name = ["A", "B", "C"]
np.random.seed(1)
data = np.array([np.random.normal(loc=i, scale=1, size=10) for i in [0.2, -1, 0.5]]).T
ns, nd = data.shapeデータのプロット
平均±標準誤差のプロット
scipy.stats.sem関数で標準誤差を計算し、平均に「ヒゲ」を付けてプロットします。そのためには、matplotlib.pyplot.errorbar関数を使用します。
白黒の図の場合（すべての点が同じ色の場合）は3点を同時にプロットすることもできますが、今回はわかりやすさを重視して、for文を使って薬ごとに点の色を変えました。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k")
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(3), drug_name)
plt.ylabel("score")
plt.show()
生データの表示
点を１列に並べる
最近の論文誌では、平均±標準誤差という要約データだけでなく、全データ点を表示することを義務付けている場合も少なくありません。
その時はmatplotlib.pyplot.scatter関数を使用して、全データ点を平均±標準誤差に重ねて表示しましょう。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k")
plt.scatter(np.tile(np.arange(nd), ns), data.reshape(-1), s=20, edgecolors=["k" for _ in range(nd)], c="w")
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(nd), drug_name)
plt.ylabel("score")
plt.show()
（補足）プロットの上下関係を調整する
ちなみに、Matplotlibでは、プロットの実行順ではなくプロットの種類ごとに重なりの上下関係が決まっています。先ほどはerrorbarの後にscatterを実行したにも関わらず、scatterの上にerrorbarが重なっていたのはそのためです。
この順番が気に入らない場合には、zorder（Z軸方向の順番）オプションを利用して重なりの順番を入れ替えます。zorderで指定した値が大きいプロットの方が上に配置されます。
以下はzorderでプロットの順番を入れ替えた（失敗）例です。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k", zorder=0)
plt.scatter(np.tile(np.arange(nd), ns), data.reshape(-1), s=20, edgecolors=["k" for _ in range(nd)], c="w", zorder=1)
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(nd), drug_name)
plt.ylabel("score")
plt.show()
点をバラつかせて重なりを回避する
データ点を白抜きの丸で表示する他にも、各色の点で表現することもできます。この場合は点の重なりを避けるため、numpy.random.normal関数を使用し、X軸方向に微小なバラツキを与えると良いでしょう。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k")
plt.scatter(np.random.normal(loc=0, scale=0.05, size=data.shape[0]) + i, data[:, i], s=3)
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(3), drug_name)
plt.ylabel("score")
plt.show()
対応のある点を線でつなぐ
また、今回は3種類の薬を同じ被験者に与えることを想定しているため、それぞれの点には対応があります。
この関係を表示する際にはmatplotlib.pyplot.plot関数を使用し、対応する点を繋ぎます。線が多いと見にくくなるため、ここではalphaパラメータを0.2に指定して線を透過しています。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k")
plt.plot(data.T, c="k", alpha=0.2)
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(3), drug_name)
plt.ylabel("score")
plt.show()
統計検定の実行
検定用のデータ構造を作る
投与した薬によって、scoreに差が出るかを検定します。
そのための準備として、データをpandas.DataFrameにまとめておきます。ここでIDは被験者に割り振られた番号、drugは投与した薬を意味します。
このデータフレームを作るとき、numpy.repeat関数とnumpy.tile関数を使うと非常に便利です。[0, 1, 2]という配列にそれぞれを適用すると、repeatは[0, ..., 0, 1, ..., 1, 2, ..., 2]、tileは[0, 1, 2, 0, 1, 2, ...]という配列を作ります。
df = pd.DataFrame({"ID": np.repeat(np.arange(ns), nd),
"drug": np.tile(drug_name, ns),
"score": data.reshape(-1)})
df[:6] # 上から６つ目までを表示
対応のある一元配置分散分析
今回は被験者間に対応があるので、対応のある一元配置分散分析を使用して薬ごとの効果の違いを検定します。
statsmodels.stats.anova.AnovaRMクラスにdata（データフレーム）、depvar（従属変数＝計測値）、subject（サンプルの対応関係）、within（影響を比較したい因子群）を与え、fitメソッドを実行することで検定結果を取得できます。
入力
anova = AnovaRM(data=df, depvar="score", subject="ID", within=["drug"]).fit()
print(anova)出力
Anova
==================================
F Value Num DF Den DF Pr > F
----------------------------------
drug 9.0438 2.0000 18.0000 0.0019
==================================Tukey法で多重検定（post-hoc検定）
分散分析の結果に有意差があったので、Tukey法で多重検定を行い、投与した薬のグループ間有意差を検証します（post-hoc test）。
statsmodels.stats.multicomp.pairwise_tukeyhsdにendog（計測値の配列）とgroups（影響を比較したい因子の配列）を与えると、検定結果を取得できます。
alphaは有意差の閾値です。表のp-adjの値がこれを下回ると、rejectがTrueになります。
入力
tukey = pairwise_tukeyhsd(endog=df["score"], groups=df["drug"], alpha=0.05)
print(tukey)出力
Multiple Comparison of Means - Tukey HSD, FWER=0.05
====================================================
group1 group2 meandiff p-adj lower upper reject
----------------------------------------------------
A B -1.2724 0.0314 -2.4451 -0.0998 True
A C 0.4841 0.5688 -0.6886 1.6567 False
B C 1.7565 0.0026 0.5838 2.9292 True
----------------------------------------------------検定結果（有意差）のプロット
「***」や「n.s.」を表示するための関数
統計検定の結果に基づいて、プロットに有意差の情報を足します。そのために、significance関数を作成しました。
また、検定のp値によって表示する記号が決まっている場合は、以下のp_markのような関数を用意しておくと便利です。
def significance(text, y, xmin, xmax, bottom_left=None, bottom_right=None):
plt.hlines(y=y, xmin=xmin, xmax=xmax, color="k", lw=1)
if bottom_left is not None:
plt.vlines(x=xmin, ymin=bottom_left, ymax=y, color="k", lw=1)
if bottom_right is not None:
plt.vlines(x=xmax, ymin=bottom_right, ymax=y, color="k", lw=1)
plt.text((xmax + xmin)/2, y, text, ha="center", va="bottom")
def p_mark(pvalue):
if pvalue < 0.001:
return "***"
elif pvalue < 0.01:
return "**"
elif pvalue < 0.05:
return "*"
else:
return "n.s."| 記号 | p値 |
|---|---|
| *** | p< 0.001 |
| ** | p< 0.01 |
| * | p< 0.05 |
| n.s. | それ以上 |
多重検定の結果を表示（足あり）
significance関数の各引数は上図のように指定します。ここでは、多重検定の結果をプロットに重ねて表示します。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k")
plt.plot(data.T, c="k", alpha=0.2)
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(3), drug_name)
plt.ylabel("score")
significance(p_mark(tukey.pvalues[0]), 1.2, 0.02, 0.98, 0.7, -0.6)
significance(p_mark(tukey.pvalues[1]), 1.5, 1.02, 1.98, -0.6, 1)
significance(p_mark(tukey.pvalues[2]), 2.2, -0.02, 2.02, 0.7, 1)
plt.ylim(-3.5, 3)
plt.show()
分散分析の結果を表示（足なし）
bottom_left, bottom_rightは指定しなくても構いません。その場合は対応する足が描画されません。
ここでは例として、分散分析の結果を表示します。
余談ですが、AnovaRMで計算した結果のp値はAnovaResults.anova_table["Pr > F"].valuesで取得できます。
plt.figure(figsize=(4, 3))
mu = data.mean(axis=0)
sem = stats.sem(data)
for i in range(3):
plt.errorbar(i, mu[i], yerr=sem[i], fmt="o", markersize=10, capsize=5, ecolor="k", markeredgecolor="k")
plt.plot(data.T, c="k", alpha=0.2)
plt.xlim(-0.5, 2.5)
plt.xticks(np.arange(3), drug_name)
plt.ylabel("score")
significance(p_mark(anova.anova_table["Pr > F"].values[0]), 1.5, 0, 2)
plt.ylim(-3.5, 3)
plt.show()
もっと知りたいこと、感想を教えてください！