「数値シミュレーションで読み解く統計のしくみ」をPythonで写経～第５章5.4「一元配置分散分析のシミュレーション」

2024年9月14日 11:06

第５章「統計的検定の論理とエラー確率のコントロール」

書籍の著者　小杉考司先生、紀ノ定保礼先生、清水裕士先生

この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第５章「統計的検定の論理とエラー確率のコントロール」の5.4節「一元配置分散分析のシミュレーション」の Python写経活動 を取り扱います。

今回は一元配置分散分析のシミュレーションです。

R の基本的なプログラム記法はざっくり Python の記法と近いです。
コードの文字の細部をなぞって、R と Python の両方に接近してみましょう！
ではテキストを開いてシミュレーションの旅に出発です🚀

はじめに

テキスト「数値シミュレーションで読み解く統計のしくみ」のご紹介

このシリーズは書籍「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」（技術評論社、「テキスト」と呼びます）の Python写経です。

テキストは、2023年9月に発売され、副題「Rでためしてわかる心理統計」のとおり、統計処理に定評のある R の具体的なコードを用いて、心理統計の理解に役立つ数値シミュレーションを実践する素晴らしい書籍です。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」初版第１刷、著者小杉考司・紀ノ定保礼・清水裕士、技術評論社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

5.4 一元配置分散分析のシミュレーション

テキスト冒頭の分散分析の概説を引用いたします。

・分散分析は３群以上の平均値の差を同時に検定する手法
・データに対応がある場合の分散分析を参加者間計画と呼ぶ
・データに対応がない場合の分散分析を参加者内計画と呼ぶ

テキストより引用（一部文章を改変）

一元配置分散分析は一要因参加者間計画の分散分析です。

この記事は Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。
概ね確率分布の特性値算出には scipy.stats を用い、乱数生成には numpy.random.generator を用いるようにしています。

主に利用するライブラリをインポートします。

### インポート

# 数値計算
import numpy as np

# データフレーム
import pandas as pd

# 確率・統計
import scipy.stats as stats
import statsmodels.api as sm            # 分散分析関連
import statsmodels.formula.api as smf   # 分散分析関連
import pingouin as pg                   # 統計便利ツール

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

5.4.1 検定の多重性

テキストによると一元配置分散分析は複数の群の平均値の差を検定するための方法です。
以前の記事では、２群の平均値の差の検定に「ｔ検定」を用いました。

テキストは、ｔ検定を繰り返して複数の群の平均値の検定をすることに問題があることを、シミュレーションで明らかにします。

■ ｔ検定を繰り返して３群の平均値の差を検定するときのタイプⅠエラー確率をシミュレーション
３群のペア「群１と群２」「群１と群３」「群２と群３」ごとにｔ検定を実行し、少なくとも１つの検定がエラーになる確率を有意水準$${\alpha=0.05}$$のもとで計算します。
10000回の検定をシミュレーションします。
scipy.stats の ttest_ind() でｔ検定を行います。

### 205ページ 検定の多重性、エラー確率のシミュレーション

## 設定と準備
alpha = 0.05       # 有意水準
k = 3              # 群の数
n = [20, 20, 20]   # サンプルサイズ
mu = [5, 5, 5]     # 群ごとの母平均
sigma = [2, 2, 2]  # 群ごとの母標準偏差
iter = 10000       # シミュレーション回数

# 結果を格納するオブジェクト
pvalue12 = np.zeros(iter)
pvalue13 = np.zeros(iter)
pvalue23 = np.zeros(iter)

## シミュレーション
rng = np.random.default_rng(seed=1234)
for i in range(iter):
    Y1 = rng.normal(size=n[0], loc=mu[0], scale=sigma[0])
    Y2 = rng.normal(size=n[1], loc=mu[1], scale=sigma[1])
    Y3 = rng.normal(size=n[2], loc=mu[2], scale=sigma[2])

    pvalue12[i] = stats.ttest_ind(Y1, Y2, equal_var=False).pvalue
    pvalue13[i] = stats.ttest_ind(Y1, Y3, equal_var=False).pvalue
    pvalue23[i] = stats.ttest_ind(Y2, Y3, equal_var=False).pvalue

# 少なくとも1つの検定が間違えている確率を計算
type1_error = np.where(
    (pvalue12 < alpha) | (pvalue13 < alpha) | (pvalue23 < alpha), 1, 0).mean()

# 出力
type1_error

【実行結果】
少なくとも１つの検定がエラーになる確率～タイプⅠエラー確率はおよそ$${12\%}$$になりました。

テキストによると、このシミュレーションで使ったデータの場合、タイプⅠエラー確率は数学的に$${1-(1-\alpha)^3}$$で求められる、とのこと。
では計算してみましょう。

### 206ページ どれか１つ間違えてしまうタイプⅠエラーの計算(タイプⅠエラーの理論値)
round(1 - (1 - alpha)**3, 10)

【実行結果】
タイプⅠエラー確率の理論値は$${14.26\%}$$です。

検定を複数回繰り返して、そのうち有意な結果を報告することによって、タイプⅠエラー確率が有意水準$${\alpha}$$を超えてしまう問題を「検定の多重性の問題」と呼ぶそうです。
テキストは検定の多重性の対策として「平均値の差の検定の文脈では、分散分析と呼ばれる方法によって部分的に解決されています」と書いています。

5.4.2 一元配置分散分析のデータ生成

一元配置分散分析による３群の平均値の差の検定に進みます。
３群の母平均を$${\mu_1, \mu_2, \mu_3}$$とすると、一元配置分散分析の帰無仮説は「$${\mu_1=\mu_2=\mu_3}$$」となり、対立仮説は「３つの母平均の少なくとも１つは他と異なる」となります。

■ データ作成① 群ラベルの作成
１～３の群を識別するラベルを作成します。

### 207ページ 一元配置分散分析用の3群データの生成１ ラベル作成

## 設定と準備
k = 3              # 群の数
n = [20, 20, 20]   # サンプルサイズ
mu = [6, 5, 4]     # 群ごとの母平均
sigma = [2, 2, 2]  # 群ごとの母標準偏差、ただし分散分析では共通

# 各群について、1,2,3という数列を20ずつ作る ※factor型ではない
x = np.hstack([np.repeat(a=1, repeats=n[0]),
               np.repeat(a=2, repeats=n[1]),
               np.repeat(a=3, repeats=n[2])])
# 中身を確認
x

【実行結果】
３つの群をそれぞれ 20個作成しました。

■ データ作成② 正規分布乱数で各群のデータを作成
前のコードで設定した群ごとの母平均と母標準偏差を用いて正規分布乱数を生成します。
テキストと異なる乱数を生成したので、一元配置分散分析の分析結果はテキストと異なることをご了承ください。

### 207 ページ　一元配置分散分析用の3群データの生成２ 3つの群のデータ（乱数）を作成
rng = np.random.default_rng(seed=1)
Y1 = rng.normal(size=n[0], loc=mu[0], scale=sigma[0])
Y2 = rng.normal(size=n[1], loc=mu[1], scale=sigma[1])
Y3 = rng.normal(size=n[2], loc=mu[2], scale=sigma[2])
Y = np.hstack([Y1, Y2, Y3])

# データフレーム化 xをカテゴリ型に変換
df = pd.DataFrame({'x': pd.Categorical(x), 'Y': Y})
print('df.shape:', df.shape)
display(df.head())

【実行結果】
$${x}$$列は群のラベル（どのグループに所属するか）、$${y}$$列は分析対象の変数です。

■ 一元配置分散分析の実行
Pythonには一元配置分散分析に使えるライブラリ等が複数あるので、私チョイスで４つの一元配置分散分析を実践します！

① statsmodels の anova_lm
まず線形回帰 ols を行い、線形回帰の結果を anova_lm に渡します。

### 207ページ 一元配置分散分析の実行 statsmodelsのanova_lmの使用
result_lm = smf.ols(formula='Y ~ x', data=df).fit()
sm.stats.anova_lm(result_lm, typ=1)

【実行結果】
$${F}$$値の$${p}$$値（PR(>F)）に注目します。
有意水準を$${5\%}$$とする場合、$${p}$$値$${0.000371}$$は有意水準を下回っているので帰無仮説は棄却され、対立仮説「３つの母平均の少なくとも１つは他と異なる」が採択されることになります。

ちなみに$${F}$$値は$${F}$$分布に従う変数であり、$${F}$$分布の自由度は$${(2, 57)}$$です。
df列の自由度を使っています。

② statsmodels の anova_oneway
線形回帰を行う必要はありません。

### 別解１statsmodelsのanova_onewayの使用
sm.stats.anova_oneway(data=Y, groups=x, use_var='equal').tuple

【実行結果】
カッコ内の左が$${F}$$値、右が$${p}$$値です。
最初の一元配置分散分析と同じ結果を得られました。

③ pingouin の anova
私の一推しのライブラリ pingouin を使います🐧

### 別解２ pingouinのanovaの使用
pg.anova(dv='Y', between='x', data=df, detailed=True)

【実行結果】
F列が$${F}$$値、p-unc列が$${p}$$値です。
最初の一元配置分散分析と同じ結果を得られました。

④ scipy.stats の f_oneway
群別にデータを分けてから f_oneway を実行します。

### 別解３ scipy.statsのf_onewayの使用

# GroupByオブジェクトの生成
groups = df.groupby(['x'], sort=False, observed=False)

# 水準別データの作成
data_by_category = []
for category in list(groups.groups):
    data_by_category.append(groups.get_group((category,))['Y'])

# 一元配置分散分析の実行
stats.f_oneway(*data_by_category)

【実行結果】
statisticが$${F}$$値、pvalue列が$${p}$$値です。
最初の一元配置分散分析と同じ結果を得られました。

■ 一元配置分散分析の$${F}$$値に迫る！
自由度$${(2, 57)}$$の$${F}$$分布の確率密度関数を描画し、有意水準$${5\%}$$を示す棄却域を描きましょう。

### 208ページ 自由度(2, 57)のF分布の確率密度関数の描画 209ページ図5.11

## 描画領域の設定
plt.figure(figsize=(6, 3))

## F分布の確率密度関数の描画
# x軸の値（Xの描画範囲）
line_x = np.linspace(0, 8)
# F分布の確率密度関数の描画：stats.f.pdf()で確率密度を取得
plt.plot(line_x, stats.f.pdf(line_x, dfn=2, dfd=57))

## 有意水準5%の棄却域の塗りつぶし描画
# 棄却限界値Critical valueの算出 上側5%点=95%点
cv = stats.f.ppf(q=0.95, dfn=2, dfd=57)
# x軸の値（棄却域のXの描画範囲）
line_x_cv = np.linspace(cv, 8)
# 棄却域の塗りつぶし描画
plt.fill_between(line_x_cv, 0, stats.f.pdf(line_x_cv, dfn=2, dfd=57),
                 ec='gray', alpha=0.3, label='棄却域')

## ラベル、タイトル、凡例の設定
plt.xlabel('x')
plt.ylabel('Density')
plt.title('自由度(2, 57)のF分布の確率密度関数')
plt.legend();

【実行結果】

ここからは、最初に作ったデータを使って$${F}$$値と$${p}$$値を算出します！

■ $${F}$$値算出関数の定義
テキストの関数をPython化します。

### 209ページ F値を計算するための関数を定義
# 各群の要素数（テキストのn）を関数内で算出するように改造
# Pythonのリスト内包表記を用いたサンプルコードをコメントで書いてます

def Fvalue_cul(Y, x):
    # 群ラベル(1, 2, 3)と各群の要素数(20, 20, 20)を取得
    groups, num_group = np.unique(x, return_counts=True)
    
    # 群の数
    k = groups.size

    # プールされた不偏分散
    # テキストのコード
    u2_p = 0
    for j in range(k):
        u2_p = u2_p + num_group[j] * Y[x==groups[j]].var(ddof=0)
    u2_p = u2_p / (num_group.sum() - k)

    # # リスト内包表記
    # u2_p = sum(np.array([n * Y[x==g].var(ddof=0)
    #                      for (g, n) in zip(groups, num_group)])) \
    #        / (num_group.sum() - k)

    # F値の計算
    # テキストのコード
    Fvalue = 0
    for j in range(k):
        Fvalue = Fvalue + (Y[x==groups[j]].mean() - Y.mean())**2 \
                 / (u2_p / num_group[j])
    Fvalue = Fvalue / (k - 1)

    # # リスト内包表記
    # Fvalue = sum(np.array([(Y[x==g].mean() - Y.mean())**2 / (u2_p / n)
    #                        for (g, n) in zip(groups, num_group)])) \
    #          / (k - 1)
    
    return Fvalue

■ $${F}$$値・$${p}$$値の算出
この関数にデータを与えて$${F}$$値を計算します。

### 210ページ F値の計算
Fvalue = Fvalue_cul(Y, x)
Fvalue

【実行結果】
一元配置分散分析の$${F}$$値とほぼ同じ結果を得られました。

続いて$${F}$$値に関する$${p}$$値を算出します。
scipy.stats を利用します。

### 210ページ p値の計算 f.sf()で上側確率を算出
stats.f.sf(x=Fvalue, dfn=k - 1, dfd=sum(n) - k)

【実行結果】
一元配置分散分析の$${p}$$値とほぼ同じ結果を得られました。

5.4.3 分散分析のタイプⅠエラー確率のシミュレーション

テキストは「$${F}$$分布を用いた一元配置分散分析の検定が、タイプⅠエラー確率を有意水準$${\alpha}$$未満にコントロールできているか確認するシミュレーション」に突入します！

以前の記事でも使った$${p}$$値のヒストグラムを描画する関数を定義します。

### p値のヒストグラム描画関数の定義

def plot_hist_pvalue(pvalue, bins=20):
    # 描画領域の設定
    plt.figure(figsize=(6, 3))
    # ヒストグラムの描画
    plt.hist(pvalue, bins=bins, ec='white', alpha=0.7)
    # x軸ラベル、y軸ラベル、タイトルの設定
    plt.xlabel('p value')
    plt.ylabel('Frequency')
    plt.title('Histogram of pvalue')
    
    plt.show()

■ 分散分析のタイプⅠエラー確率のシミュレーション
群の数＝３、各群の標本サイズ＝20で正規分布乱数を生成して、有意水準$${\alpha=0.05}$$、シミュレーション回数$${10000}$$で$${F}$$値と$${p}$$値を算出し、有意になった割合をシミュレーションします。

### 210ページ 

## 設定と準備
alpha = 0.05        # 有意水準
k = 3               # 群の数
n = [20, 20, 20]    # 群ごとのサンプルサイズ
mu = [5, 5, 5]      # 群ごとの母平均（群間で差が無い）
sigma = [2, 2, 2]   # 群ごとの母標準偏差
iter = 10000        # シミュレーション回数
x = np.hstack([np.repeat(a=1, repeats=n[0]),  # factor型ではない
               np.repeat(a=2, repeats=n[1]),
               np.repeat(a=3, repeats=n[2])])

# 結果を格納するオブジェクト
pvalue = np.zeros(iter)

## シミュレーション
rng = np.random.default_rng(seed=0)
for i in range(iter):
    Y1 = rng.normal(size=n[0], loc=mu[0], scale=sigma[0])
    Y2 = rng.normal(size=n[1], loc=mu[1], scale=sigma[1])
    Y3 = rng.normal(size=n[2], loc=mu[2], scale=sigma[2])
    Y = np.hstack([Y1, Y2, Y3])
    Fvalue = Fvalue_cul(Y, x)
    pvalue[i] = stats.f.sf(x=Fvalue, dfn=k - 1, dfd=sum(n) - k)

## 結果
# 有意になった割合
(pvalue < alpha).mean()

【実行結果】
タイプⅠエラー確率は有意水準$${\alpha=0.05}$$にとどまっていることを確認できました。

$${p}$$値の分布を確認します。

### 211ページ p値の分布 図5.12
plot_hist_pvalue(pvalue)

【実行結果】
$${p}$$値の分布は一様分布になっていることを確認できました。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルはＲ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

この記事が気に入ったらサポートをしてみませんか？