Architectural Ergonomics - Statistical Analysis

Statistical Analysis

Step00: Google Colaboratory

Google Colaboratoryにアクセスし、Pythonでコーディングする
ChatGPTやCopilotに要望を伝えれば、とりあえず動くコードを作ってくれる
何をしているか解らないコードも、「＊＊ってなに？」と聞けば教えてくれる

Step00: Visual Studio Code (VSCode)

Visual Studio Code でPython入門【Windows編】

Step01：データの傾向を確認

取得したデータをPythonで取り込み、seabornライブラリでいろんなグラフを作ってみよう
データの可視化については、「指標・特徴量の設計から始めるデータ可視化学入門」を参照すること

GitHubにサンプルコードも公開されていて勉強になる＞https://github.com/tkEzaki/data_visualization

反復測定のダミーデータの生成

スクリプトの練習用データフレームを作成
実験データはすべてロング形式にまとめる

import numpyimport pandas
# サンプルデータの生成（重複測定データ）
participant = numpy.repeat(range(10), 16) # 被験者10名に全組合せcondition_A = numpy.tile(numpy.repeat(['F', 'M'], 80), 1) # 被験者間因子（例えば性別）condition_B = numpy.tile(numpy.repeat(['J', 'K'], 1), 80) # 条件1は交互condition_C = numpy.tile(numpy.repeat(['O', 'P', 'Q', 'R'], 2), 20) # 条件2は2回繰り返しcondition_D = numpy.tile(numpy.repeat(['X', 'Y'], 8), 10) # 条件3は3回繰り返し
numpy.random.seed(0)measurement = numpy.random.normal(900, 150, 160) # 基本的な測定値
measurement[condition_C == 'Q'] += 100 # 条件 'Q' で有意差？measurement[condition_C == 'R'] -= 50 # 条件 'R' で有意差？measurement = numpy.round(measurement, 2) # 小数点以下2桁まで
# データフレームにまとめるdata = pandas.DataFrame({ 'Participant': participant, 'Condition_A': condition_A, 'Condition_B': condition_B, 'Condition_C': condition_C, 'Condition_D': condition_D, 'Measurement': measurement})
# CSV形式のデータとして出力data.to_csv('DATA.csv',index=False)print(data)

CSV形式のデータの取り込み

計測データがある場合はここからはじめる

# CSV形式のデータの取り込みdata = pandas.read_csv('DATA.csv')print(data)

ストリッププロットで可視化

外れ値を含め全てのデータをプロットする
データ数が多いときは色分けやJitterで調整

import matplotlib.pyplot as pltimport pandas # あえて省略（as pd）していませんimport seaborn # あえて省略（as sns）していません
# CSVファイルからデータを読み込む# 上のコードを実行していれば不要data = pandas.read_csv('DATA.csv')
# ストリッププロットの作成plt.figure(figsize=(90/25.4, 90/25.4))plt.rcParams.update({'font.size': 7})seaborn.stripplot(x='Condition_C', y='Measurement', hue='Condition_B', data=data, jitter=True)
plt.title('Measurement by Condition C with Condition B')plt.xlabel('Condition C')plt.ylabel('Measurement')plt.legend(title='Condition B')
plt.savefig("plot.svg", bbox_inches="tight")plt.show()

バイオリンプロットで可視化

データのピークや分布を視覚的に把握する
中央値や四分位範囲など統計量も示せる

import matplotlib.pyplot as pltimport pandas # あえて省略（as pd）していませんimport seaborn # あえて省略（as sns）していません
# CSVファイルからデータを読み込む# 上のコードを実行していれば不要data = pandas.read_csv('DATA.csv')
# バイオリンプロットの作成、線を細く設定plt.figure(figsize=(90/25.4, 90/25.4))plt.rcParams.update({'font.size': 7})seaborn.violinplot(x='Condition_C', y='Measurement', hue='Condition_B', data=data, split=True, linewidth=0.5)
plt.title('Measurement by Condition C with Condition B')plt.xlabel('Condition C')plt.ylabel('Measurement')plt.legend(title='Condition B')
plt.savefig("violin_plot_thin_lines.svg", bbox_inches="tight")plt.show()

スロープグラフで可視化（被験者ごと）

個々のデータの変化のパターンを把握する

import matplotlib.pyplot as pltimport pandas # あえて省略（as pd）していませんimport seaborn # あえて省略（as sns）していません
# CSVファイルからデータを読み込む# 上のコードを実行していれば不要data = pandas.read_csv('DATA.csv')
# Seabornカラーパレットを使用して色を設定palette = seaborn.color_palette("husl", 4)colors = {'O': palette[0], 'P': palette[1], 'Q': palette[2], 'R': palette[3]}
# プロットfig, ax = plt.subplots(figsize=(90/25.4, 90/25.4))plt.rcParams.update({'font.size': 7})
# 各被験者についてループfor participant in data['Participant'].unique(): participant_data = data[data['Participant'] == participant] # Condition Bの水準ごとにプロット for condition in participant_data['Condition_C'].unique(): # Condition Cの水準に応じて色を設定 condition_data = participant_data[participant_data['Condition_C'] == condition] ax.plot(condition_data['Condition_B'], condition_data['Measurement'], '-o', color=colors[condition], label=f'Condition C: {condition}')
# 重複するレジェンドを削除handles, labels = plt.gca().get_legend_handles_labels()by_label = dict(zip(labels, handles))plt.legend(by_label.values(), by_label.keys(), title="Condition C")
plt.title ('Slope Graph for Each Participant across Condition B')plt.xlabel('Condition B')plt.ylabel('Measurement')
plt.savefig("plot.svg", bbox_inches="tight")plt.show()

スロープグラフで可視化（平均の推移）

群としてのデータの挙動を把握する

import matplotlib.pyplot as pltimport pandas # あえて省略（as pd）していませんimport seaborn # あえて省略（as sns）していません
# CSVファイルからデータを読み込む# 上のコードを実行していれば不要# data = pandas.read_csv('DATA.csv')
# Seabornカラーパレットを使用して色を設定palette = seaborn.color_palette("husl", 4)colors = {'O': palette[0], 'P': palette[1], 'Q': palette[2], 'R': palette[3]}
# スロープグラフの作成（横軸はCondition B、Condition Cをレジェンドとして色分け）
# Condition Bごとにデータを分割data_j = data[data['Condition_B'] == 'J']data_k = data[data['Condition_B'] == 'K']
# 条件Cごとの平均値を計算means_j = data_j.groupby('Condition_C')['Measurement'].mean()means_k = data_k.groupby('Condition_C')['Measurement'].mean()
# プロットfig, ax = plt.subplots(figsize=(90/25.4, 90/25.4))plt.rcParams.update({'font.size': 7})
# Condition Cの水準ごとに色分けしてプロットfor condition in ['O', 'P', 'Q', 'R']: plt.plot(['J', 'K'], [means_j[condition], means_k[condition]], marker='o', color=colors[condition], label=condition)
plt.title('Slope Graph for Average Measurement across Condition B')plt.xlabel('Condition B')plt.ylabel('Average Measurement')plt.legend(title='Condition C')
plt.savefig("plot.svg", bbox_inches="tight")plt.show()

スロープグラフで可視化（重ね合わせ）

色分けをしっかりすれば見やすくなる？

import matplotlib.pyplot as pltimport pandas # あえて省略（as pd）していませんimport seaborn # あえて省略（as sns）していません
# CSVファイルからデータを読み込む# 上のコードを実行していれば不要data = pandas.read_csv('DATA.csv')
# 平均の推移は少し濃いめの色設定darker_colors = {cond: color for cond, color in zip(['O', 'P', 'Q', 'R'], seaborn.color_palette("husl", n_colors=4, desat=1.0))}
# 個々の被験者の推移は淡めの色設定light_colors = {cond: color for cond, color in zip(['O', 'P', 'Q', 'R'], seaborn.color_palette("husl", n_colors=4, desat=0.2))}
fig, ax = plt.subplots(figsize=(90/25.4, 90/25.4))plt.rcParams.update({'font.size': 7})
# 各被験者ごとのデータを淡い色でプロットfor participant in data['Participant'].unique(): participant_data = data[data['Participant'] == participant] for condition in ['O', 'P', 'Q', 'R']: condition_data = participant_data[participant_data['Condition_C'] == condition] ax.plot(condition_data['Condition_B'], condition_data['Measurement'], marker='o', linestyle='-', color=light_colors[condition], alpha=0.5)
# 被験者全体の平均値の推移を計算して重ねるfor condition in ['O', 'P', 'Q', 'R']: mean_transition = data[data['Condition_C'] == condition].groupby('Condition_B')['Measurement'].mean() ax.plot(mean_transition.index, mean_transition.values, '-o', label=condition, color=darker_colors[condition], linewidth=2)
plt.title('Individual and Mean Transition across Condition B')plt.xlabel('Condition B')plt.ylabel('Measurement')plt.legend(title='Condition C')
plt.savefig("plot.svg", bbox_inches="tight")plt.show()

Step02：分散分析（Python > R > anovakun）

GoogleColabのPython環境でRを動かします ↓↓↓

PythonでコーディングしながらシームレスにRの関数（anovakunなど）を使いたい
Pythonの分散分析は、高次の分析に対応しておらず、ヘンな結果が出ることもある
分散分析の結果を確認するためだけに、そのつどRやEZRを立ち上げるのはのは面倒
GoogleColabでもR言語を扱うことができるけど、別ファイルを使い分けるのも面倒

一元配置分散分析（対応あり）

rpy2をバージョン指定でインストール

# rpy2をバージョンを指定してインストールし展開!pip install rpy2==3.5.1%load_ext rpy2.ipython

Pythonでワイド型データを生成

# Pythonで「対応あり一元配置」のデータを生成import pandas as pdimport numpy as np
# 乱数のシードを設定np.random.seed(123)
# データフレームを作成participants = np.repeat(np.arange(1, 11), 3) # 各参加者が3条件すべてに参加conditions = np.tile(['条件A', '条件B', '条件C'], 10) # 条件scores = np.concatenate([ np.random.normal(100, 10, 10), # 条件Aのスコア np.random.normal(105, 10, 10), # 条件Bのスコア np.random.normal(110, 10, 10) # 条件Cのスコア])
df = pd.DataFrame({ 'participant': participants, 'condition': conditions, 'score': scores})
# 条件を列として横に展開して参加者ごとにまとめるdata = df.pivot(index='participant', columns='condition', values='score').reset_index(drop=True)
# 最左列の被験者番号を削除data.columns.name = None
data.to_csv('DATA.csv',index=False)print(data)

生成したデータをR領域に取り込む

# Pythonで生成したデータをR領域に取り込む%%R -i datadata

あるいはCSV形式のデータの取り込み

%%R# CSVファイルを読み込むdata <- read.csv("DATA.csv", header = TRUE, sep = ",")
# データの概要を表示head(data)

anovakunで分散分析（1要因参加者内）

%%R # R領域展開
# 下記URLより最新バージョンのanovakunをダウンロード# http://riseki.php.xdomain.jp/index.php?ANOVA君# colabのディレクトリに追加（ドラッグ&ドロップでOK）
# anovakunファイルの読み込みsource("anovakun-489.txt")#　ファイル名はダウンロードした最新のファイル名に合わせる
# anovakunによる分散分析の実行anovakun(wide_df, "sA", 3, auto=T, eta=T)
# 引数は、(データ名、"要因計画"、各要因の水準数、効果量の指定)# 2要因参加者内分散分析：anovakun(dat, "sAB", 2, 3, auto=T, eta=T)# 2要因参加者間分散分析：anovakun(dat, "ABs", 2, 3, eta=T)# 2要因混合分析：anovakun(dat, "AsB", 3, 2, auto=T, eta=T)# https://psycho.hes.kyushu-u.ac.jp/link/howtoanovakun/

一元配置分散分析（対応なし）

rpy2をバージョン指定でインストール

# rpy2をバージョンを指定してインストールし展開!pip install rpy2==3.5.1%load_ext rpy2.ipython

Pythonでロング型データを生成

import numpy as npimport pandas as pd
# 結果の再現性を保証np.random.seed(123)
# 参加者ID、条件、テストスコアを生成participants = np.repeat(np.arange(1, 11), 3) # 10人の参加者が3条件すべてに参加conditions = np.tile(['条件A', '条件B', '条件C'], 10) # 条件scores = np.concatenate([ np.random.normal(100, 10, 10), # 条件Aのスコア np.random.normal(105, 10, 10), # 条件Bのスコア np.random.normal(110, 10, 10) # 条件Cのスコア])
# データフレームを作成data = pd.DataFrame({ 'condition': conditions, 'score': scores})
# データの概要を表示data.to_csv('DATA.csv',index=False)print(data.head())

生成したデータをR領域に取り込む

# Pythonで生成したデータをR領域に取り込む%%R -i datahead(data)

あるいはCSV形式のデータの取り込み

%%R# CSVファイルを読み込むdata <- read.csv("DATA.csv", header = TRUE, sep = ",")
# データの概要を表示head(data)

anovakunで分散分析（1要因参加者間）

%%R # R領域展開
# 下記URLより最新バージョンのanovakunをダウンロード
# Chromeでは接続できないので、EdgeかFireFoxで接続# http://riseki.php.xdomain.jp/index.php?ANOVA君# ドラッグ&ドロップでcolabのディレクトリに追加
# anovakunファイルの読み込みsource("anovakun-489.txt")#　ファイル名はダウンロードした最新のファイル名に合わせる
# anovakunによる分散分析の実行anovakun(data, "As", 3, eta=T)# anovakunの引数は、(データ名、"要因計画"、各要因の水準数、効果量の指定)# 2要因参加者内分散分析：anovakun(dat, "sAB", 2, 3, auto=T, eta=T)# 2要因参加者間分散分析：anovakun(dat, "ABs", 2, 3, eta=T)# https://psycho.hes.kyushu-u.ac.jp/link/howtoanovakun/

多元配置分散分析（対応あり）

anovakunで分散分析（たとえば1要因被験者間・3要因参加者内）
mixed-design ANOVA with one between-subjects factor and three within-subjects factors

%%R # R領域展開
# 下記URLより最新バージョンのanovakunをダウンロード
# Chromeでは接続できないので、EdgeかFireFoxで接続# http://riseki.php.xdomain.jp/index.php?ANOVA君# ドラッグ&ドロップでcolabのディレクトリに追加
# anovakunファイルの読み込みsource("anovakun-489.txt")#　ファイル名はダウンロードした最新のファイル名に合わせる
# anovakunによる分散分析の実行anovakun(data, "AsBCD", 3, eta=T)# anovakunの引数は、(データ名、"要因計画"、各要因の水準数、効果量の指定)# 2要因参加者内分散分析：anovakun(dat, "sAB", 2, 3, auto=T, eta=T)# 2要因参加者間分散分析：anovakun(dat, "ABs", 2, 3, eta=T)# https://psycho.hes.kyushu-u.ac.jp/link/howtoanovakun/

ANOVA君/より高度な入力方式 - 井関龍太のページ（Edgeで接続）4.4.0以降のANOVA君ではロング形式のデータが扱えるようになった。
anovakun(data, "AsB", long = TRUE)
とlong=TRUEと宣言する。ロング形式の場合は水準数の指定は省略できる。データ形式は、最左端の列を参加者ラベル、最右端の列を従属変数とする。中間列は各要因とし、被験者間要因を左に，被験者内要因を右にまとめる。
たとえば、被験者間因子が（A）で、被験者内因子が（B、C、D）ならば
anovakun(data, "AsBCD", long = TRUE)
とし、AとBCDをs（subject）で区切れば、多元配置分散分析が実行される。その他の引数はこれまでどおり使用できるので、下記のような形になる。
anovakun(data, "AsBCD", long = TRUE, mau = T, auto = T )
データに分析に使わない不要な列が含まれているとエラーがでます。

Step03: 発表用グラフの作成

ｍatplot.pyplotによるグラフの生成

サンプルデータの生成（反復測定）

import numpyimport pandas
# サンプルデータの生成（重複測定データ）
participant = numpy.repeat(range(10), 16) # 被験者10名に全組合せcondition_A = numpy.tile(numpy.repeat(['F', 'M'], 80), 1) # 被験者間因子（例えば性別）condition_B = numpy.tile(numpy.repeat(['J', 'K'], 1), 80) # 条件1は交互condition_C = numpy.tile(numpy.repeat(['O', 'P', 'Q', 'R'], 2), 20) # 条件2は2回繰り返しcondition_D = numpy.tile(numpy.repeat(['X', 'Y'], 8), 10) # 条件3は3回繰り返し
numpy.random.seed(0)measurement = numpy.random.normal(900, 150, 160) # 基本的な測定値
measurement[condition_C == 'Q'] += 100 # 条件 'Q' で有意差？measurement[condition_C == 'R'] -= 50 # 条件 'R' で有意差？measurement = numpy.round(measurement, 2) # 小数点以下2桁まで
# データフレームにまとめるdata = pandas.DataFrame({ 'Participant': participant, 'Condition_A': condition_A, 'Condition_B': condition_B, 'Condition_C': condition_C, 'Condition_D': condition_D, 'Measurement': measurement})
# CSV形式のデータとして出力data.to_csv('DATA.csv',index=False)print(data)

CSV形式のデータの取り込み

# CSV形式のデータの取り込みdata = pandas.read_csv('DATA.csv')print(data)

エラーバー付き折れ線グラフの出力

（統計的注釈なし）

import seaborn as snsimport matplotlib.pyplot as plt
# グラフエリアを作成fig = plt.figure(figsize=(90/25.4, 90/25.4))ax = fig.add_subplot(111)ax.spines['top'].set_linewidth(0.5)ax.spines['right'].set_linewidth(0.5)ax.spines['bottom'].set_linewidth(0.5)ax.spines['left'].set_linewidth(0.5)ax.grid(color=(0.8, 0.8, 0.8), linewidth=0.2)
# プロット作成sns.pointplot( data=data, x=data.columns[3], y=data.columns[5], hue=data.columns[4], dodge=0.1, palette="dark:gray", markers="o", markersize=4, linestyle="-", linewidth=0.3, markeredgecolor='black', markeredgewidth=0.2, errorbar="se", capsize=0.05, ax=ax )
# 凡例の設定legend = ax.legend(loc='lower right', bbox_to_anchor=(1, 0), fontsize=6, ncol=2)legend.get_frame().set_linewidth(0)ax.tick_params(axis='both', which='major', labelsize=7)plt.xlabel(data.columns[4], fontsize=7)plt.ylabel(data.columns[5], fontsize=7)
#　科学研究のグラフとして必要な注釈の追加plt.text(1, 1.01, "Error bars denote Mean ± SE.", ha='right', fontsize=6, transform=plt.gca().transAxes)
plt.savefig("plot.svg", bbox_inches="tight")plt.show()

多重比較検定の結果に応じて注釈

エラーバー付き折れ線グラフの出力

（主効果あり・交互作用なし）

# グラフの作成 (主効果あり交互作用なし)import seaborn as snsimport matplotlib.pyplot as plt
# グラフエリアを作成fig = plt.figure(figsize=(90/25.4, 90/25.4))ax = fig.add_subplot(111)ax.spines['top'].set_linewidth(0.5)ax.spines['right'].set_linewidth(0.5)ax.spines['bottom'].set_linewidth(0.5)ax.spines['left'].set_linewidth(0.5)ax.grid(color=(0.8, 0.8, 0.8), linewidth=0.2)
# プロット作成sns.pointplot(    data=data, x=data.columns[3], y=data.columns[5], hue=data.columns[4],    dodge=0.1, palette="dark:gray",    markers="o", markersize=4, linestyle="-", linewidth=0.3,    markeredgecolor='black', markeredgewidth=0.2,    errorbar="se", capsize=0.05,    ax=ax    )
# 凡例の設定legend = ax.legend(loc='lower right', bbox_to_anchor=(1, 0), fontsize=6, ncol=2)legend.get_frame().set_linewidth(0)ax.tick_params(axis='both', which='major', labelsize=7)plt.xlabel(data.columns[4], fontsize=7)plt.ylabel(data.columns[5], fontsize=7)
#　科学研究のグラフとして必要な注釈の追加plt.text(1, 1.01, "Error bars denote Mean ± SE. * indicates p < 0.05 significance.",         ha='right', fontsize=6, transform=plt.gca().transAxes)
# 統計的注釈の追加# 位置調整用の数値の設定dodge = 0.07y_max = ax.get_ylim()[1]y_len = (data[data.columns[5]].max() - data[data.columns[5]].min()) / 120
# plt.ylim(800, 1150)
# 統計的注釈の形状を定義def add_stat_annotation(x1, x2, y, dg):    plt.plot([x1+dg, x1+dg, x2+dg, x2+dg], [y, y + y_len/2, y + y_len/2, y], lw=0.2, c='k')    plt.text((x1 + x2)/2+dg, y, "*", ha='center', va='bottom', fontsize=6)
# 水準間の有意差(交互作用なし)def group_stat_annotation(x1, x2, y, dg):    plt.plot([x1+dg, x1+dg, x2+dg, x2+dg], [y, y + y_len/2, y + y_len/2, y], lw=0.2, c='k')    plt.plot([(x1+x2)/2, (x1+x2)/2], [y+ y_len/2, y + y_len ], lw=0.2, c='k')
levels = data['Condition_C'].unique()level_pairs_g = [(0, 0), (1, 1), (2, 2), (3, 3)]for x1, x2 in level_pairs_g:    group_stat_annotation(x1-0.1, x2+0.1, y_max, 0)
y_max += y_len
significant_pairs =[        (0, 1), (0, 2), (0, 3),        (1, 2), (1, 3),        (2, 3)         ] #有意差が出た水準間のみ残す
for x1, x2 in significant_pairs:    add_stat_annotation(x1, x2, y_max, 0)    y_max += y_len # 次の注釈のためにy_maxを更新

エラーバー付き折れ線グラフの出力

（交互作用あり・単純主効果あり）

# グラフの作成import seaborn as snsimport matplotlib.pyplot as plt
# グラフエリアを作成fig = plt.figure(figsize=(90/25.4, 90/25.4))ax = fig.add_subplot(111)ax.spines['top'].set_linewidth(0.5)ax.spines['right'].set_linewidth(0.5)ax.spines['bottom'].set_linewidth(0.5)ax.spines['left'].set_linewidth(0.5)ax.grid(color=(0.8, 0.8, 0.8), linewidth=0.2)
# プロット作成sns.pointplot(    data=data, x=data.columns[3], y=data.columns[5], hue=data.columns[4],    dodge=0.1, palette="dark:gray",    markers="o", markersize=4, linestyle="-", linewidth=0.3,    markeredgecolor='black', markeredgewidth=0.2,    errorbar="se", capsize=0.05,    ax=ax    )
# 凡例の設定legend = ax.legend(loc='lower right', bbox_to_anchor=(1, 0), fontsize=6, ncol=2)legend.get_frame().set_linewidth(0)ax.tick_params(axis='both', which='major', labelsize=7)plt.xlabel(data.columns[4], fontsize=7)plt.ylabel(data.columns[5], fontsize=7)
#　科学研究のグラフとして必要な注釈の追加plt.text(1, 1.01, "Error bars denote Mean ± SE. * indicates p < 0.05 significance.",         ha='right', fontsize=6, transform=plt.gca().transAxes)
# 統計的注釈の追加# 位置調整用の数値の設定dodge = 0.07y_max = ax.get_ylim()[1]y_len = (data[data.columns[5]].max() - data[data.columns[5]].min()) / 120
# plt.ylim(800, 1150)
# 統計的注釈の形状を定義def add_stat_annotation(x1, x2, y, dg):    plt.plot([x1+dg, x1+dg, x2+dg, x2+dg], [y, y + y_len/2, y + y_len/2, y], lw=0.2, c='k')    plt.text((x1 + x2)/2+dg, y, "*", ha='center', va='bottom', fontsize=6)
# 統計的注釈を総当りで仮置き# 水準内の有意差level_pairs_0 = [                (0, 0), (1, 1), (2, 2), (3, 3)                ] # 有意差が出た水準間のみ残す
for x1, x2 in level_pairs_0:    add_stat_annotation(x1-dodge, x2+dodge, y_max, 0)y_max += y_len # 次の注釈のためにy_maxを更新
# 水準間の有意差(交互作用あり)level_pairs_1 =[        (0, 1), (0, 2), (0, 3),        (1, 2), (1, 3),        (2, 3)         ] #有意差が出た水準間のみ残す
for x1, x2 in level_pairs_1:    add_stat_annotation(x1-dodge, x2-dodge, y_max, 0)    y_max += y_len # 次の注釈のためにy_maxを更新
# 水準間の有意差(交互作用あり)level_pairs_2 =[        (0, 1), (0, 2), (0, 3),        (1, 2), (1, 3),        (2, 3)         ] #有意差が出た水準間のみ残す
for x1, x2 in level_pairs_2:    add_stat_annotation(x1+dodge, x2+dodge, y_max, 0)    y_max += y_len # 次の注釈のためにy_maxを更新
plt.savefig("plot.svg", bbox_inches="tight")

Step04: 論文用テーブルの作成

統計的注釈の自動化（未完成）

主効果あり・交互作用なしの事後検定（significant_pairsの出力）

# 主効果あり交互作用なしの場合の事後検定from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey_result = pairwise_tukeyhsd(data['Measurement'], data['Condition_C'], alpha=0.05)print (tukey_result)
# グラフに統計的注釈を記述するための組み合わせの書き出しsave_data = data['Condition_C'] # 一旦別名で保存data['Condition_C'] = data['Condition_C'].astype('category').cat.codes # 水準名を数値に
tukey_result2 = pairwise_tukeyhsd(data['Measurement'], data['Condition_C'], alpha=0.05)significant_pairs = [tuple(tukey_result2._results_table.data[i+1][:2]) for i, _ in enumerate(tukey_result.reject) if tukey_result2.reject[i]]print ("有意差のあった群間：", significant_pairs)print ("\n")
data['Condition_C']= save_data # 数値を水準名に戻す

グラフの生成（事後検定のsignificant_pairsを拾わせる）

# グラフの作成 (主効果あり交互作用なし)import seaborn as snsimport matplotlib.pyplot as plt
# グラフエリアを作成fig = plt.figure(figsize=(90/25.4, 90/25.4))ax = fig.add_subplot(111)ax.spines['top'].set_linewidth(0.5)ax.spines['right'].set_linewidth(0.5)ax.spines['bottom'].set_linewidth(0.5)ax.spines['left'].set_linewidth(0.5)ax.grid(color=(0.8, 0.8, 0.8), linewidth=0.2)
# プロット作成sns.pointplot(    data=data, x=data.columns[3], y=data.columns[5], hue=data.columns[4],    dodge=0.1, palette="dark:gray",    markers="o", markersize=4, linestyle="-", linewidth=0.3,    markeredgecolor='black', markeredgewidth=0.2,    errorbar="se", capsize=0.05,    ax=ax    )
# 凡例の設定legend = ax.legend(loc='lower right', bbox_to_anchor=(1, 0), fontsize=6, ncol=2)legend.get_frame().set_linewidth(0)ax.tick_params(axis='both', which='major', labelsize=7)plt.xlabel(data.columns[4], fontsize=7)plt.ylabel(data.columns[5], fontsize=7)
#　科学研究のグラフとして必要な注釈の追加plt.text(1, 1.01, "Error bars denote Mean ± SE. * indicates p < 0.05 significance.",         ha='right', fontsize=6, transform=plt.gca().transAxes)
# 統計的注釈の追加# 位置調整用の数値の設定dodge = 0.07y_max = ax.get_ylim()[1]y_len = (data[data.columns[5]].max() - data[data.columns[5]].min()) / 120
# plt.ylim(800, 1150)
# 統計的注釈の形状を定義def add_stat_annotation(x1, x2, y, dg):    plt.plot([x1+dg, x1+dg, x2+dg, x2+dg], [y, y + y_len/2, y + y_len/2, y], lw=0.2, c='k')    plt.text((x1 + x2)/2+dg, y, "*", ha='center', va='bottom', fontsize=6)
# 水準間の有意差(交互作用なし)def group_stat_annotation(x1, x2, y, dg):    plt.plot([x1+dg, x1+dg, x2+dg, x2+dg], [y, y + y_len/2, y + y_len/2, y], lw=0.2, c='k')    plt.plot([(x1+x2)/2, (x1+x2)/2], [y+ y_len/2, y + y_len ], lw=0.2, c='k')
levels = data['Condition_C'].unique()level_pairs_g = [(0, 0), (1, 1), (2, 2), (3, 3)]for x1, x2 in level_pairs_g:    group_stat_annotation(x1-0.1, x2+0.1, y_max, 0)
y_max += y_len
for x1, x2 in significant_pairs:    add_stat_annotation(x1, x2, y_max, 0)    y_max += y_len # 次の注釈のためにy_maxを更新

交互作用あり・単純主効果ありの事後検定（Level_pairs_1,2,3の出力）

# 交互作用ありの場合の事後検定（単純主効果あり）from statsmodels.stats.multicomp import pairwise_tukeyhsd
#################################### Condition_Dの水準に基づいてデータを分割groups = data['Condition_C'].unique()
# 各水準内の群間に対して多重比較検定for group in groups:    subset = data[data['Condition_C'] == group]    # Tukeyの多重比較検定を実施し出力    tukey_result = pairwise_tukeyhsd(subset['Measurement'], subset['Condition_D'], alpha=0.05)    print(f"{group}の水準内での多重比較検定:")    print(tukey_result)
# 統計的注釈を記述するための組み合わせの書き出しlevel_pairs_0 = []save_data = data['Condition_C'] # 一旦別名で保存data['Condition_C'] = data['Condition_C'].astype('category').cat.codes # 水準名を数値に
for group in groups:    subset = data[data['Condition_C'] == group]    # Tukeyの多重比較検定を実施し出力    tukey_result = pairwise_tukeyhsd(subset['Measurement'], subset['Condition_D'], alpha=0.05)    significant_pairs = [tuple(tukey_result._results_table.data[i+1][:2]) for i, _ in                          enumerate(tukey_result.reject) if tukey_result.reject[i]]    level_pairs_0.extend(significant_pairs)print(level_pairs_0)print("\n")
data['Condition_C']= save_data # 数値を水準名に戻す
#################################### Condition_Dの水準に基づいてデータを分割groups = data['Condition_D'].unique()
# 各水準内の群間に対して多重比較検定for group in groups:    subset = data[data['Condition_D'] == group]
    # Tukeyの多重比較検定を実施し出力    tukey_result = pairwise_tukeyhsd(subset['Measurement'], subset['Condition_C'], alpha=0.05)    print(f"{group}の水準内での多重比較検定:")    print(tukey_result)    print("\n")
# 統計的注釈を記述するための組み合わせの書き出しsave_data = data['Condition_C'] # 一旦別名で保存data['Condition_C'] = data['Condition_C'].astype('category').cat.codes # 水準名を数値に
# 有意なペアを保存するための辞書level_pairs_dict = {}for idx, group in enumerate(groups, start=1):    subset = data[data['Condition_D'] == group]
    # 有意な差を持つペアのみを抽出し出力    tukey_result = pairwise_tukeyhsd(subset['Measurement'], subset['Condition_C'], alpha=0.05)    level_pairs = [tuple(tukey_result._results_table.data[i+1][:2]) for i, _ in                         enumerate(tukey_result.reject) if tukey_result.reject[i]]        # globalsを使って動的に変数を作成    globals()[f'level_pairs_{idx}'] = level_pairs        # 辞書に結果を保存    level_pairs_dict[f"level_pairs_{idx}"] = level_pairs
# 例: 第1のグループの有意なペアを参照するにはprint (level_pairs_1)print (level_pairs_2)

グラフの生成（事後検定のLevel_pairs_1と2を拾わせる）

# グラフの作成import seaborn as snsimport matplotlib.pyplot as plt
# グラフエリアを作成fig = plt.figure(figsize=(90/25.4, 90/25.4))ax = fig.add_subplot(111)ax.spines['top'].set_linewidth(0.5)ax.spines['right'].set_linewidth(0.5)ax.spines['bottom'].set_linewidth(0.5)ax.spines['left'].set_linewidth(0.5)ax.grid(color=(0.8, 0.8, 0.8), linewidth=0.2)
# プロット作成sns.pointplot(    data=data, x=data.columns[3], y=data.columns[5], hue=data.columns[4],    dodge=0.1, palette="dark:gray",    markers="o", markersize=4, linestyle="-", linewidth=0.3,    markeredgecolor='black', markeredgewidth=0.2,    errorbar="se", capsize=0.05,    ax=ax    )
# 凡例の設定legend = ax.legend(loc='lower right', bbox_to_anchor=(1, 0), fontsize=6, ncol=2)legend.get_frame().set_linewidth(0)ax.tick_params(axis='both', which='major', labelsize=7)plt.xlabel(data.columns[4], fontsize=7)plt.ylabel(data.columns[5], fontsize=7)
#　科学研究のグラフとして必要な注釈の追加plt.text(1, 1.01, "Error bars denote Mean ± SE. * indicates p < 0.05 significance.",         ha='right', fontsize=6, transform=plt.gca().transAxes)
# 統計的注釈の追加# 位置調整用の数値の設定dodge = 0.07y_max = ax.get_ylim()[1]y_len = (data[data.columns[5]].max() - data[data.columns[5]].min()) / 120
# plt.ylim(800, 1150)
# 統計的注釈の形状を定義def add_stat_annotation(x1, x2, y, dg):    plt.plot([x1+dg, x1+dg, x2+dg, x2+dg], [y, y + y_len/2, y + y_len/2, y], lw=0.2, c='k')    plt.text((x1 + x2)/2+dg, y, "*", ha='center', va='bottom', fontsize=6)
# 水準内の有意差に統計的注釈for x1, x2 in level_pairs_0:    add_stat_annotation(x1-dodge, x2+dodge, y_max, 0)y_max += y_len # 次の注釈のためにy_maxを更新
# 水準間の有意差に統計的注釈for x1, x2 in level_pairs_1:    add_stat_annotation(x1-dodge, x2-dodge, y_max, 0)    y_max += y_len # 次の注釈のためにy_maxを更新
for x1, x2 in level_pairs_2:    add_stat_annotation(x1+dodge, x2+dodge, y_max, 0)    y_max += y_len # 次の注釈のためにy_maxを更新
plt.savefig("plot.svg", bbox_inches="tight")

参考サイト

ANOVA関連

R関連

そのほか

基本解説＞ポストホックテストとしての多重比較検定
SPSS ＞多元配置分散分析（対応あり）
Excelの分析ツールの有効化＞分散分析：エクセル標準とRとPythonの比較
- 実験データ解析再入門：論文を「フェイクニュース」にしないために

多重比較の前に分散分析は必要か？

理論的には必要ありません。むしろ、Bonferroni法やTukeyの多重比較検定は分散分析とは異なる仕組みで解析するため、分散分析を経ずに実施するべきです。
しかし、分散分析を行わずに多重比較検定を適用した結果を認めない立場の研究者もとても多く存在します。特に、論文の査読過程でそうした研究者にあたって修正の指摘を受けると、結果の解釈や結論の大幅な変更をしなくてはならなくなります。そのため、実質的には分散分析を行っておく方が、安全かつ円滑に研究を進めやすいといえます。
それでも分散分析を通さずに多重比較検定を実施する場合は下記のような説明を明示してください。

本研究では、特定の条件間の比較に焦点を当てており、交互作用項にも興味はない。したがって分散分析による全体的な有意性検定は必要ではないと判断した。異なる統計検定量を扱う検定を重ねることによる多重性の問題も考慮した判断である」
＞永田靖(1998)「多重比較法の実際」『応用統計学』27 (2): 93-108 https://doi.org/10.5023/jappstat.27.93.＞永田靖・𠮷田道弘 (1997) 『統計的多重比較法の基礎』サイエンティスト社