「数値シミュレーションで読み解く統計のしくみ」をPythonで写経～第６章6.5「非心分布を使わないサンプルサイズ設計のシミュレーション」

2024年10月19日 11:07

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

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

この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第６章「適切な検定のためのサンプルサイズ設計」の6.5節「非心分布を使わないサンプルサイズ設計のシミュレーション」の Python写経活動 を取り扱います。

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

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

はじめに

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

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

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

引用表記

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

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

6.5 非心分布を使わないサンプルサイズ設計のシミュレーション

この節では、複雑なモデルの場合に非心分布でタイプⅡエラー確率の計算ができないことを想定して、タイプⅡエラー確率自体をシミュレーションして必要サンプルサイズを計算する方法を学びます。

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

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

### インポート

# 数値計算
import numpy as np

# 確率・統計
import scipy.stats as stats
from scipy.special import expit         # ロジスティック関数（シグモイド関数）
import statsmodels.formula.api as smf   # 回帰分析等の統計処理

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

テキストはシミュレーションの方法を次のように説明しています。

データ生成に必要なパラメータを決めて、そのパラメータに基づいてデータを何度も発生させ、検定を繰り返し実施します。
ここで、母数（パラメータ）が０ではないモデルからデータを発生させているのに有意にならなければタイプⅡエラーに該当するので、有意にならない回数をシミュレーション回数で割ったものがタイプⅡエラー確率の推定値になります。

テキストの文章を一部改変して引用

■ 対応のないｔ検定のシミュレーション
差$${\delta}$$、各群の標準偏差$${\sigma}$$を設定して、iter_t2e回シミュレーションを行って、タイプⅡエラー確率を推定する関数を定義します。
関数中の単回帰分析には scipy.stats の linregress を使用しています。

### 251ページ 対応のない2標本の差のt検定のタイプⅡエラー確率シミュレーション関数
# 単回帰分析で対応のないt検定のタイプⅡエラー確率が算出されるのが不思議･･･

def t2e_ttest(alpha, delta, sigma, n, iter_t2e, seed=1234):

    ## 設定と準備
    # Xの作成: 0をn個、1をn個
    X = np.repeat(a=[0, 1], repeats=n)
    # p値を格納する変数の初期化
    pvalue = np.zeros(iter_t2e)
    # 乱数生成器・乱数シードの設定
    rng = np.random.default_rng(seed=seed)

    ## シミュレーション
    for i in range(iter_t2e):
        # 変数Yの作成: 正規分布乱数（x=0のとき平均0, x=1のとき平均delta）
        Y = np.hstack([rng.normal(size=n, loc=0, scale=sigma),
                       rng.normal(size=n, loc=delta, scale=sigma)])
        # 'Y ~ 1 + X'の回帰分析を実行 scipy.statsのlinregress()を使用
        result = stats.linregress(x=X, y=Y)
        # Xの回帰係数のp値を取得
        pvalue[i] = result.pvalue

    ## 有意でない（p値が有意水準α以上の）割合を算出
    t2e = (alpha <= pvalue).mean()

    return t2e

20000回のシミュレーションでタイプⅡエラー確率を推定します。
各標本のサンプルサイズは$${50}$$です。

### 251ページ 対応のない2標本の差のt検定のタイプⅡエラー確率シミュレーション実行

## 設定と準備
alpha = 0.05   # 有意水準
delta = 1      # 母平均の差
sigma = 2      # 各標本の母標準偏差
n = 50         # 各標本のサンプルサイズ
iter = 20000   # シミュレーション回数

## シミュレーション
t2e_ttest(alpha, delta, sigma, n, iter_t2e=iter, seed=123)

【実行結果】
各標本のサンプルサイズが$${50}$$のとき、タイプⅡエラー確率は$${29.8\%}$$と推定されました。

■ サンプルサイズを変化させて$${\beta=20\%}$$に近づける
各標本のサンプルサイズを$${100}$$と$${64}$$にしてシミュレーションします。
まずはサンプルサイズ$${100}$$のケースです。

### 252ページ 続き サンプルサイズ100でシミュレーション
n = 100
t2e_ttest(alpha, delta, sigma, n, iter_t2e=iter, seed=123)

【実行結果】

続いてサンプルサイズ$${64}$$のケースです。

### 252ページ 続き サンプルサイズ64でシミュレーション
n = 64
t2e_ttest(alpha, delta, sigma, n, iter_t2e=iter, seed=123)

【実行結果】
$${\beta=20}$$にかなり近づきました。

前々回の記事にて、本シミュレーションと同じ条件で、非心ｔ分布を用いて、対応のないｔ分布の必要サンプルサイズ（理論値）を算出したところ、各群$${64}$$人、タイプⅡエラー確率$${19.9\%}$$でした。
本シミュレーションのタイプⅡエラー確率は、理論値とかなり近いことが分かりました。

非心分布が分からなくてもデータ生成ができれば、シミュレーションの方法を用いることでサンプルサイズ設計ができる、ということです！

■ ロジスティック回帰のサンプルサイズ設計
テキストに合わせて、ロジスティック回帰の回帰係数に関するタイプⅡエラー確率計算を関数化します。

### 252ページ ロジスティック回帰分析のタイプⅡエラー確率シミュレーション関数

def t2e_logistic(alpha, b0, b1, n, iter_t2e, seed=123):
    
    ## 設定と準備
    # 説明変数の生成。ここでは2群のダミー変数（非確率変数）を想定
    X = np.repeat(a=[0, 1], repeats=n)
    # p値を格納する変数の初期化
    pvalue = np.zeros(iter_t2e)
    # 乱数生成器・乱数シードの設定
    rng = np.random.default_rng(seed=seed)
    
    ## シミュレーション
    for i in range(iter_t2e):
        # ベルヌーイ分布乱数を２種類生成 expit(x)はロジスティック関数
        Y = np.hstack([rng.binomial(size=n, n=1, p=expit(b0)),
                       rng.binomial(size=n, n=1, p=expit(b0 + b1))])
        # ロジスティック回帰の実行
        result = smf.logit('Y ~ X', data=dict(X=X, Y=Y)).fit(disp=False)
        # Xの係数のp値を取得
        pvalue[i] = result.pvalues.X

    ## エラー確率の計算 ：p値が有意水準α以上の割合
    t2e = (alpha <= pvalue).mean()

    return t2e

20000回のシミュレーションを実行します。

### 253ページ 続き ロジスティック回帰分析のタイプⅡエラー確率シミュレーションの実行
# 課題：statsmodelsのログを抑制する方法 → fit(disp=False)で解決 40秒

## 設定と準備
b0 = -1.5     # 線形予測子の切片
b1 = 0.8      # 線形予測子の傾き
n = 98        # 各群のサンプルサイズ
iter = 20000  # シミュレーション回数

## シミュレーション
seed = 0
t2e_logistic(alpha, b0, b1, n=n, iter_t2e=iter, seed=seed)

【実行結果】
タイプⅡエラー確率は約$${33\%}$$です。

6.6 演習問題

6.6.1 演習問題１
非心分布を使わないで、相関係数のサンプルサイズ設計をシミュレーションするコードを書きます。
まずはシミュレーション関数の定義です。

### 253ページ 6.6.1 演習問題１
# 相関係数のサンプルサイズ設計を非心分布を使わずに計算する

## 分散共分散行列を作る関数
def Vcov(sigma, rho):
    size = len(sigma)
    mat = np.diag(sigma**2).astype(float)
    for i in range(size - 1):
        for j in range(i+1, size):
            mat[i, j] = mat[j, i] = sigma[i] * sigma[j] * rho
    return mat

## 相関係数のタイプⅡエラー確率シミュレーション関数
def t2error_corr(n, mu, sigma, rho, alpha, iter=10000, seed=123):

    ## 設定と準備
    # 有意でない場合１を立てる変数の初期化
    reject_result = np.zeros(iter)
    # 乱数生成器・乱数シードの設定
    rng = np.random.default_rng(seed=seed)

    ## シミュレーション
    for i in range(iter):
        # ２変量正規分布乱数の生成
        sim = rng.multivariate_normal(size=n, mean=mu, cov=Vcov(sigma, rho))
        # 相関係数の検定の実行 scipy.statsのpearsonr()を使用
        result = stats.pearsonr(sim[:, 0], sim[:, 1])
        # 相関係数が有意でない（p値が有意水準α以上の）場合、1を設定
        reject_result[i] = np.where(alpha <= result.pvalue, 1, 0)

    ## 有意でない（p値が有意水準α以上の）割合を算出
    type2_error = reject_result.mean()

    return type2_error

必要サンプルサイズのシミュレーションを実行します。
シミュレーション回数は 10000 回です。

### 253ページ 続き シミュレーションの実行 45秒

## 設定と準備
alpha = 0.05               # 有意水準
beta = 0.20                # タイプⅡエラー確率の閾値（使っていない変数）
mu = [5, 5]                # 母平均
sigma = np.array([2, 3])   # 母標準偏差
rho = 0.3                  # 母相関係数
iter = 10000               # シミュレーション回数
n_ptn = range(80, 86)      # チェックするサンプルサイズのパターン

## シミュレーション
seed = 1
for n in n_ptn:
    type2_error = t2error_corr(n, mu, sigma, rho, alpha, iter, seed)
    print(f'Sample size {n:3d},  Error ratio {type2_error:.5f}')

【実行結果】
$${\beta=0.20}$$未満となる必要サンプルサイズは$${84}$$です。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

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

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

この記事が参加している募集

#私のこだわり旅

with ことりっぷ

6,073件

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