「数値シミュレーションで読み解く統計のしくみ」をPythonで写経～第５章5.1、5.2「統計的検定の論理、対応のないｔ検定」

2024年8月31日 11:38

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

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

この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第５章「統計的検定の論理とエラー確率のコントロール」の5.1節「統計的検定の論理」と5.2節「Rによる統計的検定の実際」の Python写経活動 を取り扱います。

今回は対応のないｔ検定の実践を通じて統計的検定の論理に近づきます。

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

はじめに

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

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

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

引用表記

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

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

5.1 統計的検定の論理

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

5.1節では特段、プログラムコードを用いた解説はありません。
そこで、テキストの２つの図の描画に取り組みます。

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

### インポート

# 数値計算
import numpy as np

# 確率・統計
import scipy.stats as stats
import pingouin as pg                   # 統計便利ツール
from statsmodels.stats.weightstats import ttest_ind

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

5.1.2 帰無仮説のデータによる棄却

１つ目の図はテキストの図5.1 ｔ検定の帰無仮説の分布（ｔ分布）の描画です。
ｔ分布の各種数値の算出には scipy.stats の t を用います。
t.pdf で確率密度の取得、t.ppf でパーセンタイル点の取得が可能です。

### 189ページ 図5.1 t検定の帰無分布（t分布）

## 設定
n = 20         # 標本サイズ
nu = n - 2     # 自由度
alpha = 0.05   # 有意水準

## 計算値
lower = stats.t.ppf(q=alpha/2, df=nu)         # 下側の棄却限界値
upper = stats.t.ppf(q=1-alpha/2, df=nu)       # 上側の棄却限界値
line_x = np.linspace(-4, 4, 1001)             # 確率密度関数算出に用いる確率変数X
line_x_95ci = np.linspace(lower, upper, 1001) # 受容域の確率変数X

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 3))
# t分布の確率密度関数の描画
plt.plot(line_x, stats.t.pdf(line_x, df=nu))
# 受容域の塗りつぶしの描画
plt.fill_between(line_x_95ci, 0, stats.t.pdf(line_x_95ci, df=nu), ec='black',
                 alpha=0.2)
# タイトル
plt.title('検定統計量 $t$ が帰無仮説が正しいときに従う分布');

【実行結果】
曲線がｔ分布の確率密度関数、青塗りの領域は帰無仮説が正しいときに検定統計量$${T}$$が大抵は得られると想定される範囲（95%の範囲）です。

２つ目の図はテキストの図5.2 ｔ検定の帰無分布と非心分布です。
非心分布は帰無仮説が正しくないときに検定統計量$${T}$$が従う分布です。
非心ｔ分布の各種数値の算出には scipy.stats の nct を用います。
nct.pdf で確率密度の取得が可能です。

### 189ページ 図5.2 t検定の帰無分布と非心分布

## 設定
n = 20         # 標本サイズ
nu = n - 2     # 自由度
alpha = 0.05   # 有意水準
lam = 3        # 非心度

## 計算値
lower = stats.t.ppf(q=alpha/2, df=nu)         # 下側の棄却限界値
upper = stats.t.ppf(q=1-alpha/2, df=nu)       # 上側の棄却限界値
line_x = np.linspace(-3, 8, 1001)             # 確率密度関数算出に用いる確率変数X
line_x_95ci = np.linspace(lower, upper, 1001) # 受容域の確率変数X

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 3))
# 帰無分布t分布の確率密度関数の描画
plt.plot(line_x, stats.t.pdf(line_x, df=nu))
# 非心分布の確率密度関数の描画
plt.plot(line_x, stats.nct.pdf(line_x, df=nu, nc=lam), color='black', lw=0.7,
         ls='--')
# 受容域の塗りつぶしの描画
plt.fill_between(line_x_95ci, 0, stats.t.pdf(line_x_95ci, df=nu), ec='black',
                 alpha=0.2)
# タイトル
plt.title('検定統計量 $t$ が実際に従う分布（右）');

【実行結果】
点線が非心ｔ分布です。

5.2 Rによる統計的検定の実際

テキストのタイトルは「Rによる」ですが、この記事は「Pythonによる」になります。
5.2節では対応のない$${t}$$検定に取り組みます。

5.2.1 対応のないｔ分布

２群の平均値の差の検定です。

サンプルデータの作成、$${t}$$値の計算と有意性の判断、$${p}$$値の計算と有意性の判断、の流れで進みます。

■ 対応のない２群データの作成
テキストに従い、「設定と準備」の条件で２つの乱数を生成します。
テキストの数式をお借りします。
２データはそれぞれ独立に正規分布に従い、かつ、等分散です。

$$
Y_{1i} \sim \text{Normal}\ (\mu_1,\ \sigma^2) \\
Y_{2i} \sim \text{Normal}\ (\mu_2,\ \sigma^2) \\
$$

テキストより引用

### 192ページ ２つの標本の作成（正規分布乱数） 図5.3

## 設定と準備
n = [20, 25]   # 2つのサンプルサイズ
mu = [6, 5]    # 2つの母平均
sigma = 2      # 母標準偏差（ただしt検定のときは両母集団で等しい）
alpha = 0.05   # 有意水準

## シミュレーション
rng = np.random.default_rng(seed=7) # 1
Y1 = rng.normal(size=n[0], loc=mu[0], scale=sigma)  # 群1のデータ生成
Y2 = rng.normal(size=n[1], loc=mu[1], scale=sigma)  # 群2のデータ生成

## 結果 Y1のヒストグラムの描画
plt.figure(figsize=(6, 3))
plt.hist(Y1, bins=8, ec='white')
plt.xlabel('Y1')
plt.ylabel('Frequency')
plt.title('Histogram of Y1');

【実行結果】
Y1データのヒストグラムです。

■ $${t}$$値の計算
帰無仮説は$${\mu_1 = \mu_2}$$、あるいは、$${\mu_{\delta}=0}$$です。

### 193ページ 検定統計量tの算出

## 設定と準備
# 自由度の計算
df = n[0] + n[1] - 2

## t値の計算
# 群間の標本平均の差
Ybar_d = Y1.mean() - Y2.mean()
# プールされた不偏分散
u2_p = (n[0]*Y1.var(ddof=0) + n[1]*Y2.var(ddof=0)) / (n[0] + n[1] - 2)
# プールされた標本サイズ
n_p = (n[0] * n[1]) / (n[0] + n[1])
# t値
t = (Ybar_d - 0) / np.sqrt(u2_p / n_p)

## 出力
t

【実行結果】

■ 臨界値の計算
臨界値は棄却域の境界です。
ｔ分布の確率密度関数にて棄却域と臨界値を可視化します。

### 193ページ 図5.4 t分布の棄却域

## 計算値
line_x = np.linspace(-4, 4, 1001)            # 確率変数X
lower = stats.t.ppf(q=alpha/2, df=df)        # 下側の棄却限界値
upper = stats.t.ppf(q=1-alpha/2, df=df)      # 上側の棄却限界値
line_x_lower = np.linspace(-4, lower, 1001)  # 下側の棄却域の確率変数X
line_x_upper = np.linspace(upper, 4, 1001)   # 上側の棄却域の確率変数X

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 3))
# t分布の確率密度関数の描画
plt.plot(line_x, stats.t.pdf(line_x, df=df))
# 下側棄却域の塗りつぶしの描画
plt.fill_between(line_x_lower, 0, stats.t.pdf(line_x_lower, df=df), 
                 color='tab:blue', ec='black', alpha=0.2)
# 上側棄却域の塗りつぶしの描画
plt.fill_between(line_x_upper, 0, stats.t.pdf(line_x_upper, df=df),
                 color='tab:blue', ec='black', alpha=0.2)
# タイトル
plt.title('$t$ 分布の棄却域');

【実行結果】
青塗りの領域が棄却域、棄却域の端が臨界値です。

臨界値を計算します。

### 193ページ 臨界値の計算
c = stats.t.ppf(q=1 - alpha/2, df=df)
c

【実行結果】
臨界値は$${t}$$値よりも（絶対値が）小さい結果になりました。

■ 有意性の判断
絶対値に変換の上、$${t}$$値が臨界値を上回る場合、帰無仮説を棄却します。

### 194ページ 有意性の判断
'帰無仮説の棄却' if abs(t) > abs(c) else '帰無仮説の保留'

【実行結果】
絶対値に変換の上、$${t}$$値$${2.16}$$が臨界値$${2.02}$$を上回るので、帰無仮説を棄却します。

■ $${p}$$値の計算
テキストによると「$${p}$$値とは、データから計算された検定統計量と等しいあるいはそれより極端な（$${t}$$分布の場合は絶対値が大きい）値が得られる確率」です。
テキストの図5.5に相当する$${p}$$値の可視化に取り組みます。

### 194ページ 図5.5 t分布の棄却域とp値が示す確率の範囲

## 計算値
line_x = np.linspace(-4, 4, 1001)                   # 確率変数X
line_x_low_p = np.linspace(-4, -abs(t), 1001)       # 下側のp値の確率変数X
line_x_low_c = np.linspace(-abs(t), -abs(c), 1001)  # 下側の棄却域の確率変数X
line_x_up_p = np.linspace(abs(t), 4, 1001)          # 上側のp値の確率変数X
line_x_up_c = np.linspace(abs(c), abs(t), 1001)     # 上側の棄却域の確率変数X

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 3))
# t分布の確率密度関数の描画
plt.plot(line_x, stats.t.pdf(line_x, df=nu))
# 下側p値の塗りつぶしの描画
plt.fill_between(line_x_low_p, 0, stats.t.pdf(line_x_low_p, df=nu), 
                 color='tab:red', ec='black', alpha=0.2, label='$p$ 値')
# 下側棄却域の塗りつぶしの描画
plt.fill_between(line_x_low_c, 0, stats.t.pdf(line_x_low_c, df=nu), 
                 color='tab:blue', ec='black', alpha=0.2)
# 上側p値の塗りつぶしの描画
plt.fill_between(line_x_up_p, 0, stats.t.pdf(line_x_up_p, df=nu), 
                 color='tab:red', ec='black', alpha=0.2)
# 上側棄却域の塗りつぶしの描画
plt.fill_between(line_x_up_c, 0, stats.t.pdf(line_x_up_c, df=nu), 
                 color='tab:blue', ec='black', alpha=0.2)

# タイトル
plt.title('$t$ 分布の棄却域と $p$ 値が示す確率の範囲')
plt.legend();

【実行結果】
赤い領域が$${p}$$値に相当します。

$${p}$$値を用いて統計的検定を行います。
まずは$${p}$$値の算出です。

### 195ページ 累積分布関数を用いたp値の計算
# p値の計算
pvalue = (1 - stats.t.cdf(x=abs(t), df=df)) * 2  # 2倍は両側検定
pvalue

【実行結果】
$${p}$$値はおよそ$${0.036}$$です。

■ 有意性の判断
有意水準を$${5\%}$$とする場合の有意性の判断を実行します。

### 195ページ p値を使った有意性の判断
'帰無仮説の棄却' if pvalue < alpha else '帰無仮説の保留'

【実行結果】
$${p}$$値$${0.036}$$は有意水準$${5\%}$$($${0.05}$$)を下回るので、帰無仮説は棄却されます。

■ Python のライブラリを用いた$${t}$$検定
３つのライブラリで対応のない２標本の母平均の差のｔ検定を実践します。

① scipy.stats の statsのttest_ind を利用

### 195ページ t検定関数で実行
# scipy.statsのttest_ind: 独立した２標本の母平均の差のt検定を使用
stats.ttest_ind(Y1, Y2, equal_var=True)

【実行結果】
statisticは$${t}$$値、pvalueは$${p}$$値、dfは自由度です。

② statsmodels の ttest_ind を利用

### 別解１ statsmodelsの使用 T検定統計量, p値, 自由度
ttest_ind(Y1, Y2, usevar='pooled')

【実行結果】
左から$${t}$$値、$${p}$$値、自由度です。

③ pingouin の ttest を利用

### 別解２ pingouinの使用
pg.ttest(Y1, Y2, correction=False)

【実行結果】
T は$${t}$$値、dof は自由度、alternative の two-sided は両側検定、p-val は$${p}$$値です。

■ Welch の検定
等分散を仮定しない Welch の検定を実践します。
等分散を仮定するｔ検定と同様に３つのライブラリを利用します。

① scipy.stats の statsのttest_ind を利用
引数 equal_var=False とすることで Welch の検定を実行できます。

### 196ページ 等分散を仮定しない二群の平均値の差の検定（welch test）
stats.ttest_ind(Y1, Y2, equal_var=False)

【実行結果】
statisticは$${t}$$値、pvalueは$${p}$$値、dfは自由度です。

② statsmodels の ttest_ind を利用
引数 usevar='pooled' とすることで Welch の検定を実行できます。

### 別解１ statsmodelsの使用 T検定統計量, p値, 自由度
ttest_ind(Y1, Y2, usevar='pooled')

【実行結果】
左から$${t}$$値、$${p}$$値、自由度です。

③ pingouin の ttest を利用
引数 correction=True とすることで Welch の検定を実行できます。

### 別解２ pingouinの使用
pg.ttest(Y1, Y2, correction=True)

【実行結果】
T は$${t}$$値、dof は自由度、alternative の two-sided は両側検定、p-val は$${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の教科書です。
よかったらぜひ、お試しくださいませ。

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