Statistical Analysis
Statistical Analysis
Step00: Google Colaboratory
Step00: Google Colaboratory
Google Colaboratoryにアクセスし、Pythonでコーディングする
ChatGPTやCopilotに要望を伝えれば、とりあえず動くコードを作ってくれる
何をしているか解らないコードも、「＊＊ってなに？」と聞けば教えてくれる
Step00: Visual Studio Code (VSCode)
Step00: Visual Studio Code (VSCode)
Step01：データの傾向を確認
Step01：データの傾向を確認
- 取得したデータをPythonで取り込み、seabornライブラリでいろんなグラフを作ってみよう
- データの可視化については、「指標・特徴量の設計から始めるデータ可視化学入門」を参照すること
反復測定のダミーデータの生成
- スクリプトの練習用データフレームを作成
- 実験データはすべてロング形式にまとめる
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)
# サンプルデータの生成（重複測定データ）
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()
# 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()
# 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()
# 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()
# 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()
# 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）
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)
# 乱数のシードを設定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)
# データの概要を表示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/
# 下記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())
# 結果の再現性を保証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)
# データの概要を表示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/
# 下記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/
# 下記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 )
データに分析に使わない不要な列が含まれているとエラーがでます。
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: 発表用グラフの作成
Step03: 発表用グラフの作成
ｍatplot.pyplotによるグラフの生成
ｍ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)
# サンプルデータの生成（重複測定データ）
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()
# グラフエリアを作成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を更新
# グラフエリアを作成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")
# グラフエリアを作成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: 論文用テーブルの作成
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 # 数値を水準名に戻す
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を更新
# グラフエリアを作成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)
#################################### 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")
# グラフエリアを作成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")