この記事は（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メソッドを実行することで検定結果を取得できます。