「数値シミュレーションで読み解く統計のしくみ」をPythonで写経～第６章6.2「タイプⅡエラー確率のコントロールとサンプルサイズ設計」

2024年9月28日 11:16

第６章「適切な検定のためのサンプルサイズ設計」

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

この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第６章「適切な検定のためのサンプルサイズ設計」の6.2節「タイプⅡエラー確率のコントロールとサンプルサイズ設計」の Python写経活動 を取り扱います。

第６章はシミュレーションを通じて、調査・実験の事前段階でサンプルサイズを設計する方法を学びます。
今回はタイプⅡエラー確率の計算を学びます。
非心分布の再降臨です！

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

はじめに

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

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

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

引用表記

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

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

6.2 タイプⅡエラー確率のコントロールとサンプルサイズ設計

タイプⅡエラーは帰無仮説が偽の場合に帰無仮説を保留する誤りのことです。
本当は効果があるのに誤って効果がないといってしまうことです。

テキストは「追試の際、タイプⅡエラー確率がとても高い実験をしてしまうと、本当は効果があるのに追試失敗の烙印を押されること」、「許容できるエラー確率$${\beta}$$を極端に小さくすることは、社会の貴重な資源や負担を過剰に高めてしまう問題がある」と言う点を指摘し、タイプⅡエラー確率の適切なコントロールの必要性を提示しています。

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

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

### インポート

# 数値計算
import numpy as np

# 確率・統計
import scipy.stats as stats

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

6.2.1 非心分布

テキストによると、非心分布は「実際に検定統計量が従う分布」です。
「非心」とは母数（パラメータ）で中心化されていない量が従う分布とのことです。
イメージとしては、母数の真の値を使わないで、帰無仮説が想定する母数の仮説値を用いて検定統計量を計算するので、非心になるようです。

非心ｔ分布の場合、非心度$${\lambda}$$は、母平均$${\mu}$$、帰無仮説の母平均（想定値）$${\mu_0}$$、母標準偏差$${\sigma}$$、サンプルサイズ$${n}$$を用いて、次式のように表せるそうです。

$$
\lambda = \cfrac{\mu - \mu_0}{\sigma / \sqrt{n}}
$$

テキストより引用

Python では例えば、scipy.stats の nct で非心ｔ分布を実装できます。

ｔ分布と非心ｔ分布の確率密度関数を描画します。

### 234ページ t分布と非心t分布 235ページ図6.3

# 描画領域の設定
plt.figure(figsize=(6, 3))
# x軸の値
line_x = np.linspace(-4, 8, 1001)
# t分布の確率密度関数の描画 scipy.statsのtを使用
plt.plot(line_x, stats.t.pdf(line_x, df=25), color='tab:orange', ls='--',
         label='t 分布')
# 非心t分布の確率密度関数の描画 scipy.statsのnctを使用
plt.plot(line_x, stats.nct.pdf(line_x, df=25, nc=3), color='tab:blue',
         label='非心 t 分布')
# x軸ラベル、y軸ラベル、タイトル、凡例の設定
plt.xlabel('x')
plt.ylabel('Density')
plt.title('t 分布と非心 t 分布')
plt.legend();

【実行結果】

6.2.2 タイプⅡエラー確率の計算

非心度$${\lambda=3}$$のタイプⅡエラー確率を可視化して確かめます。
最初に描画関数を定義します。

### タイプⅠエラー確率とタイプⅡエラー確率の描画関数

def plot_statistical_error(lambda_, df, alpha, i, n=None):

    # 描画領域の設定
    plt.figure(figsize=(6, 3))
    # t分布の描画
    x = np.linspace(-3, 7, 1001)
    plt.plot(x, stats.t.pdf(x, df=df), color='tab:orange', linestyle='--')

    # 帰無分布の描画
    # 臨界値から8までの値を200区切りで用意
    xx = np.linspace(stats.t.ppf(q=1 - alpha / 2, df=df), 7, i)
    yy = stats.t.pdf(xx, df=df) # xxに対応したt分布の密度を得る
    # タイプⅠエラー確率を色付け
    plt.fill_between(xx, yy, color='tab:orange', ec='gray', alpha=0.2,
                    label='タイプⅠエラー確率')

    # 非心分布の描画
    plt.plot(x, stats.nct.pdf(x, df=df, nc=lambda_), color='tab:blue')
    # -3から臨界値までの値を200区切りで用意
    xx = np.linspace(-3, stats.t.ppf(q=1 - alpha / 2, df=df), i)
    yy = stats.nct.pdf(xx, df=df, nc=lambda_) # xxに対応した非心t分布の密度を得る
    # タイプⅡエラー確率を色付け
    plt.fill_between(xx, yy, color='tab:blue', ec='gray', alpha=0.2,
                    label='タイプⅡエラー確率')

    # タイトルの設定
    if n is None:
        title = f'非心度={lambda_}のときの\nタイプⅠエラー確率とタイプⅡエラー確率'
    else:
        title = f'n={n}のときのタイプⅡエラー確率'
    plt.title(title)
    # x軸ラベル、y軸ラベル、凡例の設定
    plt.xlabel('x')
    plt.ylabel('Density')
    plt.legend()
    plt.show()

では描画してみます。
有意水準$${\alpha=0.05}$$です。
薄い青色の部分がタイプⅡエラー確率です。

### 235ページ 非心度λ=3の場合のタイプⅡエラー確率 236ページ図6.4

# 設定と準備
lambda_ = 3    # 非心度
df = 25        # 自由度
alpha = 0.05   # 有意水準
i = 200        # 描画の細かさ

# 描画
plot_statistical_error(lambda_, df, alpha, i)

【実行結果】

続いて非心度$${\lambda=5}$$のケースを可視化します。

### 236ページ 図6.5 非心度λ=5の場合のタイプⅡエラー確率

# 設定と準備
lambda_ = 5    # 非心度
df = 25        # 自由度
alpha = 0.05   # 有意水準
i = 200        # 描画の細かさ

# 描画
plot_statistical_error(lambda_, df, alpha, i)

【実行結果】
非心度を大きくすることで、タイプⅡエラー確率はかなり小さくなりました。

テキストにならって、母効果量$${\delta_0}$$を用いて非心度を書き換えます。

$$
\delta_0 = \cfrac{\mu - \mu_0}{\sigma} \\
\\
\lambda = \delta_0 \sqrt{n}
$$

テキストより引用

テキストはタイプⅡエラー確率をコントロールするためには、母効果量とサンプルサイズを適切に設定する必要がある、と解説しています。
また、サンプルサイズを設計する場合には、母効果量を決めておく必要があるとしています。

6.2.3 母効果量の見積もりとサンプルサイズ設計

母効果量を見積もる際、テキストは想定する母効果量の値よりも小さく見積もることで、非心度を小さくしておく例を掲載しています。
その小さな非心度よりも実際の非心度が大きくなった場合であっても、非心度が大きいほどタイプⅡエラー確率は小さくなるので、想定する母効果量で定めたタイプⅡエラー確率$${\beta}$$未満に収まる、という考え方です。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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