「入門はじめての多変量解析」をPythonで写経 Vol.4 ～ ２章「はじめての重回帰分析」③重回帰分析の統計的検定
２章「はじめての重回帰分析」
書籍の著者 石村貞夫 先生、石村光資郎 先生
書籍「入門はじめての多変量解析」２章「はじめての重回帰分析」の Python写経活動記録 です。
多変量解析の入門を Python と一緒に学ぶ写経シリーズです。
この記事は、重回帰分析の統計的検定 を取り扱います。
具体的には 重回帰の分散分析表（回帰の有意性の検定）、偏回帰係数の有意性の検定 に取り組みます。
ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀
はじめに
このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。
書籍の紹介と引用表記はリンク先の記事に掲載しています。
２章 はじめての重回帰分析
この記事は２章の以下のSectionを取り扱います。
2.6 重回帰分析の検定
記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。
この記事で用いるライブラリをインポートします。
### インポート
# 数値計算
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 matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib
分析の準備
テキスト p.35 表 2.1.1 「温度・圧力・配向度のデータ」を引用いたします。
■ データの読み込み
pandas データフレームにデータを設定します。
### 温度・圧力・配向度のデータ p.35 表2.1.1
# データの登録
data1 = pd.DataFrame(
{'配向度': [45, 38, 41, 34, 59, 47, 35, 43, 54, 52],
'温度': [17.5, 17.0, 18.5, 16.0, 19.0, 19.5, 16.0, 18.0, 19.0, 19.5],
'圧力': [30, 25, 20, 30, 45, 35, 25, 35, 35, 40],
}, index=range(1, 11))
data1.index.name = 'サンプルNo.'
# 結果の表示
data1【実行結果】
データの個数（標本サイズ）は 10 です。
目的変数は「配向度」、説明変数の候補は「温度」「圧力」です。
■ 共通設定
説明変数名と目的変数名を設定します。
データフレームの列を特定する際に活躍します。
## 共通設定
# 説明変数
VARS = ['温度', '圧力']
# 目的変数
TARGET = '配向度'【実行結果】なし
■ Python ライブラリで重回帰分析結果を取得
ライブラリ statsmodels で、説明変数に「温度」「圧力」を用いて、重回帰分析を実施しておきます。
🖲️statsmodels
引数 formula で目的変数と説明変数の関係を設定します。
基本文法は「目的変数 ~ 説明変数１ ＋ 説明変数２ ＋ ・・・」です。
## 重回帰分析の実行 statsmodels利用
result1_sm = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1).fit()
result1_sm.summary()【実行結果】
重回帰分析の結果を格納した「result1_sm」は今回記事で頻出します！
重回帰分析のサマリーは以下のようになりました。
統計的検定 p.64
テキストは重回帰分析の統計的検定のうち、次の２つを紹介しています。
・回帰の有意性の検定（重回帰の分散分析表）
・偏回帰係数の有意性の検定（説明変量の検定）
統計的検定は「標本データを用いて母集団に関する仮説が統計的に成り立つかどうかを判断すること」です。
一般に次の手順で実行します（テキストの手順を一部改変）。
母集団の特性値に関する帰無仮説 $${H_0}$$、対立仮説 $${H_1}$$ を設定する
標本データより検定統計量を計算する
検定統計量が棄却域に含まれるとき、帰無仮説 $${H_0}$$ を棄却する
回帰の有意性の検定（重回帰の分散分析表） p.65～
「重回帰モデル全体」が目的変数に対して効果があるかどうか、に関する統計的検定です。
① 仮説
■ テキスト版：
・帰無仮説 $${H_0}$$：重回帰式は予測に役立たない
■ よく見かける版（$${p}$$ 個の母偏回帰係数について）
・帰無仮説 $${H_0}$$：$${\beta_1=\beta_2 = \cdots = \beta_p=0}$$
（すべての母偏回帰係数は０である）
・対立仮説 $${H_1}$$：$${\beta_1 \neq 0\ \text{or}\ \beta_2 \neq 0\ \text{or} \cdots \text{or}\ \beta_p \neq 0}$$
（少なくとも１つの母偏回帰係数は０でない）
② 検定統計量
回帰の有意性検定で用いる検定統計量は $${F}$$ 値です。
テキストは統計解析ソフト SPSS による「分散分析表」で $${F}$$ 値を計算しています。
計算手順付きの分散分析表で $${F}$$ 値の計算方法をチェックしましょう。
$$
\begin{array}{|c|c:c:c|c|}
\hline
\\
変動要因 & 平方和 & 自由度 & 平均平方（分散） & F値 \\
\\
\hline
\\
回帰 & S_R & p & V_R = S_R / p & F = V_R / V_e \\
\\
\hdashline
\\
残差 & S_E & N-p-1 & V_e = S_E / (N-p-1) & - \\
\\
\hline
\\
合計 & S_T & N-1 & - & - \\
\\
\hline
\end{array}
$$
回帰の平方和 $${S_R}$$ と残差平方和 $${S_E}$$、データの個数（標本サイズ） $${N}$$、説明変数の数 $${p}$$ で計算します。
平方和は前々回記事で紹介しました。
そして$${F}$$ 値は「自由度 $${(p,\ N-p-1)}$$ の $${F}$$ 分布」に従います。
🔢 分散分析表の作成・計算
statsmodels の回帰分析結果を用いて、分散分析表を算出してみましょう。
anova_lm 関数に回帰分析結果 result1_sm を渡します。
### 重回帰の分散分析表 F検定 p.65 表2.6.2 ※statsmodels利用
sm.stats.anova_lm(result1_sm, typ=1)【実行結果】
回帰の行が説明変数別に分かれていて、回帰全体の $${F}$$ 値を確認できません…
ということで、statsmodels の回帰分析結果を利用して分散分析表を作成する関数を作りましょう！（DIYです）
### 重回帰の分散分析表作成関数 for statsmodels p.65
def anova_table(result):
## 設定と準備
# 標本サイズn、説明変数の数p
n, p = result.nobs, result.df_model
# 欠損値
NaN = np.nan
## データフレームの作成
# statsmodelsのOLSResultから平均和、自由度、平均平方、F値、P値を設定
df = pd.DataFrame({
'平方和': [result.ess, result.ssr, result.ess + result.ssr],
'自由度': [int(p), int(n - p - 1), int(n - 1)],
'平均平方': [result.mse_model, result.mse_resid, NaN],
'F値': [result.fvalue, NaN, NaN],
'p値': [result.f_pvalue, NaN, NaN]},
index=['回帰', '残差', '全体'])
## 戻り値: 分散分析表データフレーム
return dfstatsmodels の回帰分析結果 result1_sm を引数に与えて関数を実行します。
# 分散分析表の作成
anova_table(result1_sm).round(3)【実行結果】
テキスト p.65 表 2.6.2「SPSS による分散分析表」と同じ顔つきになりました。
$${F}$$ 値は $${21.176}$$ です。
$${p}$$ 値は $${0.001 \leq 0.05}$$ です。
実は分散分析表を作成しなくても、statsmodels の回帰分析結果の「サマリー」を表示することで、 $${F}$$ 値とその $${p}$$ 値を確認できます。
# statsmodels のF値
result1_sm.summary().tables[0]【実行結果】
$${F}$$ 値は F-statistic = $${21.18}$$、$${p}$$ 値は Prob (F-statistic) = $${0.00107}$$ です。
ちなみに $${F}$$ 値だけを取り出すには：
# statsmodels のF値
result1_sm.fvalue【実行結果】
では仮説の判断に進みましょう。
③ 帰無仮説が棄却できるか判断
結論を先に述べると…
$${p}$$ 値 $${\leq}$$ 有意水準 $${5\%}$$ です。
有意水準 $${5\%}$$ で帰無仮説は棄却され、重回帰モデルは有意であると言えます。
「温度」「圧力」を説明変数にした重回帰モデルが、配向度の予測に役立つのです。
重回帰モデルを作った甲斐がありました！！！
テキスト p.65 図 2.6.1「有意確率と有意水準」に相当する可視化を行いましょう。
描画関数を作ります。
### F分布の可視化 for statsmodels
def plot_f_dist(result, alpha=0.05):
## 設定と準備
# 基本色の設定
color = 'tab:blue'
# 標本サイズn、説明変数の数p
n, p = result.nobs, result.df_model
# 重回帰モデルのF値とP値の取得
f_val, p_val = result.fvalue, result.f_pvalue
# F分布の自由度の算出
dfn, dfd = int(p), int(n-p-1)
# F分布の設定
f_dist = stats.f(dfn=dfn, dfd=dfd)
# 棄却限界値の算出
c_val = f_dist.isf(q=alpha)
# 描画のためのx軸の値
x_val1 = np.linspace(0, f_val * 1.1)
x_val2 = np.linspace(c_val, f_val * 1.1)
## 描画
# F分布の確率密度関数の描画
plt.plot(x_val1, f_dist.pdf(x_val1), color=color, label='$F$分布')
# 有意水準alphaの棄却域の塗りつぶし描画
plt.fill_between(x_val2, 0, f_dist.pdf(x_val2), color=color, alpha=0.2,
label='棄却域')
# 棄却限界値の垂直点線（黒）の描画
plt.axvline(c_val, color='black', ls='--', label='棄却限界値')
# F値の垂直点線の描画
plt.axvline(f_val, color='tab:red', ls='--', label='F値')
# 修飾
plt.title(f'自由度{dfn, dfd}の$F$分布\n'
f'棄却限界値:{c_val:.4f}, F値:{f_val:.4f}'
f'\n有意水準{alpha:.1%}, $p$値={p_val:.4f}')
plt.xticks([c_val, f_val])
plt.legend(loc='upper right')
plt.show()では描画しましょう。
### 有意水準の可視化 p.65 図2.6.1
plot_f_dist(result1_sm)【実行結果】
$${F}$$ 値が棄却限界値以上の場合、帰無仮説を棄却します。
このチャートで $${F}$$ 値 $${21.1756}$$が棄却限界値 $${4.7374}$$ を上回っているので、帰無仮説を棄却できることがわかります。
🛸ちょっと寄り道：ChatGPTおすすめチャート 🛸
もっと重回帰モデルの可視化例を知りたいと思い、ChatGPTに訊いてみました。
【プロンプト】
回帰の有意性の検定をグラフで「直感的に」理解したいと思います。
お勧めのグラフを教えてください
回答は・・・🥁
帰無仮説下での「説明優位度（R²）」のシミュレーション分布💹
決定係数 $${R^2}$$ を使って重回帰モデルの有効性をつかめるチャートとのこと。
ひとまず描画しましょう！
### 決定係数 R² の帰無分布をシミュレーション
# R² の帰無分布は、帰無仮説「説明変数と目的変数に何の関係もない」のもとで、
# 偶然でどのくらいの R² が出るかを集めた分布です。
# 観測データで得られた実際の R² が、帰無分布でほとんど現れないほど大きければ、
# 「偶然ではこんなに高い説明力はほぼ起こらない」
# ⇒「モデルは意味がある（有意）」と判断します。
## 設定と準備
# シミュレーション回数
n_sim = 1000
# 重回帰分析の決定係数
observed_r2 = result1_sm.rsquared
# 乱数生成器
rng = np.random.default_rng(seed=123)
## 帰無分析のシミュレーション
# 帰無分布の決定係数を格納するリストの初期化
r2_null = []
# シミュレーションの実行
for _ in range(n_sim):
# 目的変数をランダムに並び替える
y_perm = rng.permutation(data1[TARGET].values)
# 並び替えた目的変数と並び替えていない説明変数で重回帰分析を実行
perm_model = sm.OLS(y_perm, sm.add_constant(data1[VARS])).fit()
# 決定係数を格納
r2_null.append(perm_model.rsquared)
## 描画
# 描画領域の設定
plt.figure(figsize=(6, 4))
# 帰無分布の決定係数のヒストグラムの描画
plt.hist(r2_null, bins=30, density=True, edgecolor='white', alpha=0.7)
plt.axvline(observed_r2, ls='--', lw=2, color='tab:red',
label=f'観測 $R^2$ = {observed_r2:.2f}')
plt.xlabel('帰無分布のもとの $R^2$：帰無仮説「温度・圧力と配向度は無関係」',
fontsize=12)
plt.ylabel('密度', fontsize=12)
plt.title('$R^2$ の帰無分布シミュレーションと観測値')
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()
【実行結果】
このヒストグラムは「帰無仮説：説明変数と目的変数に関係がない」場合の分布を示しています。
具体的には、「ランダムにシャッフルした目的変数」と「元のままのシャッフルしない説明変数」で重回帰分析を行って決定係数を算出する試行を 1000 回繰り返しています。
シャッフルした目的変数と元のままの説明変数は関係がなさそうですよね！
ヒストグラムを読むと、ほとんどのデータの決定係数は $${0.8}$$ 未満であることが分かります。
そして実際のデータから得られた決定係数（観測 $${R^2}$$） は、赤い垂直点線の $${0.86}$$ です。
説明変数と目的変数に関係がない場合の分布では、決定係数 $${0.86}$$ が滅多に起きないことを示しています。
「滅多に起きない ⇒ 説明変数と目的変数に関係がある」と考えられます。
そしてそしてChatGPT は…
$${F}$$ 値と $${R^2}$$ の関係も教えてくれました。
回帰の平方和 $${S_R\ +}$$ 残差平方和 $${S_E\ =\ }$$ 総平方和 $${S_T}$$ の関係から、
$$
\begin{align*}
R^2 = \cfrac{S_R}{S_T} &\Longrightarrow S_R = R^2 \times S_T \\
R^2 = 1-\cfrac{S_E}{S_T} &\Longrightarrow \ S_E = (1-R^2) \times S_T
\end{align*}
$$
となり、$${S_R, S_T}$$ を $${F}$$ 値の計算式に代入すると、
$$
F=\cfrac{\cfrac{S_R}{p}}{\cfrac{S_E}{N-p-1}} = \cfrac{\cfrac{R^2 S_T}{p}}{\cfrac{(1-R^2)S_T}{N-p-1}} = \cfrac{\cfrac{R^2}{p}}{\cfrac{1-R^2}{N-p-1}}
$$
となります。
決定係数 $${R^2}$$ が大きいと $${F}$$ も大きくなる関係が分かりました。
偏回帰係数の有意性の検定（説明変量の検定） p.66～
「個々の説明変数」が目的変数に影響を与えているかどうか、に関する統計的検定です。
① 仮説
特定の説明変数の母偏回帰係数 $${\beta_i}$$ に関する仮説です。
■ テキスト版：
・帰無仮説 $${H_0}$$：説明変数は目的変数に影響を与えない
■ よく見かける版
・帰無仮説 $${H_0}$$：$${\beta_i=0}$$
（母偏回帰係数 $${\beta_i}$$ は０である）
・対立仮説 $${H_1}$$：$${\beta_i \neq 0}$$
（母偏回帰係数 $${\beta_i}$$ は０でない）
② 検定統計量
偏回帰係数の有意性検定で用いる検定統計量は $${t}$$ 値です。
テキストは統計解析ソフト SPSS で $${t}$$ 値を計算しています。
🔢 検定統計量 $${t}$$ 値の計算
statsmodels の回帰分析結果を用いて、$${t}$$ 値を確認しましょう。
### 検定統計量の算出 p.67 表2.6.3 ※statsmodels利用
result1_sm.summary().tables[1]【実行結果】
t 列に偏回帰係数ごとの $${t}$$ 値が表示されています。
ちなみに $${t}$$ 値は
$$
t\ 値 = \cfrac{偏回帰係数の推定値\ \text{coef}}{標準誤差\ \text{std err}}
$$
で計算できます。
🛸ちょっと寄り道：$${\boldsymbol{t}}$$ 値の公式 🛸
テキストに $${t}$$ 値の公式が未掲載だったので気になり…
ChatGPTに教えてもらいました。
なお変数の記号がテキストと異なっていますのでご注意下さい。
まず標準誤差を計算、次に $${t}$$ 値を計算します。
📊 $${j}$$ 番目の偏回帰係数の推定量 $${\hat{\beta}_j}$$ の標準誤差
$$
\begin{align*}
\text{SE}(\hat{\beta_j}) &= \sqrt{s^2 \left[(X^{\top}X)^{-1} \right]_{jj}}\\
\\
s^2 &= \cfrac{\sum_{i=1}^N (y_i - \hat{y}_i)^2}{N-p-1}\\
\end{align*}
$$
【変数・記号の説明】
$$
\begin{array}{l:l}
変数・記号 & 説明 \\
\hline
\\
X & 説明変数（定数項を含む）\\
\\
[(X^{\top}X)^{-1}]_{jj} & [\ ]の行列の\ j\ 行\ j\ 列成分 \\
\\
s^2 & 誤差分散\ \sigma^2\ の推定量 \\
& テキストでは誤差変動の不偏分散\ V_E \\
\\
y_i & 目的変数 \\
& テキストでは\ Y \\
\\
\hat{y}_i & 目的変数の予測値 \\
\\
\sum_{i=1}^N (y_i - \hat{y}_i)^2 & 残差平方和 \\
\\
N & 標本サイズ \\
\\
p & 説明変数の数（定数項を除く） \\
\end{array}
$$
📊 $${j}$$ 番目の偏回帰係数の推定量 $${\hat{\beta}_j}$$の $${t}$$ 値
$$
t = \cfrac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \sim t(N-p-1)
$$
🔢 公式を用いた検定統計量 $${t}$$ 値の計算
上の公式を配向度データに当てはめて、ステップ・バイ・ステップで $${t}$$ 値を 計算しましょう。
配向度データから説明変数と目的変数を取り出して、説明変数に定数項を追加しています。
## 設定と準備
# 説明変数
X = data1[VARS].values
# 目的変数
y = data1[TARGET].values
# 標本サイズN、説明変数の数p
N, p = X.shape
# 説明変数の最初の列に定数項を追加
X_const = np.column_stack([np.ones(N), X])偏回帰係数を算出します。
（せっかくなので偏回帰係数の計算もします）
## 偏回帰係数の推定
beta_hat = np.linalg.inv(X_const.T @ X_const) @ X_const.T @ y
beta_hat【実行結果】
切片、温度、圧力の偏回帰係数の推定値です。
偏回帰係数の標準誤差を算出します。
公式を使っていきます！
## 偏回帰係数の標準誤差の推定
# 誤差分散の推定値
sigma2_hat = sum(result1_sm.resid**2) / (N - p - 1)
# 説明変数の行列積の逆行列
XX_inv = np.linalg.inv(X_const.T @ X_const)
# 偏回帰係数の標準誤差の推定
se_beta_hat = np.sqrt(np.diag(sigma2_hat * XX_inv))
# 結果の表示
print(se_beta_hat)【実行結果】
切片、温度、圧力の標準誤差です。
statsmodels の回帰分析結果と合っています！
偏回帰係数の $${t}$$ 値と $${p}$$ 値（両側検定）を算出します。
## t値、p値の算出
# t値の算出
t_value = beta_hat / se_beta_hat
print('t値: ', t_value)
# t値のp値の算出
p_value = stats.t.sf(abs(t_value), df=N-p-1)*2 # 両側
print('p値: ', p_value)【実行結果】
切片、温度、圧力の順に表示しています。
statsmodels の回帰分析結果と合っています！
statsmodels のように、公式を用いて計算した各値を表にまとめましょう。
## 重回帰分析のサマリー表の作成
pd.DataFrame(
{'係数': beta_hat, '標準誤差': se_beta_hat, 't値': t_value, 'p値': p_value},
index=['切片'] + VARS
).round(3)【実行結果】
表にまとめると、なんだかいい感じですね！
では仮説の判断に進みましょう。
③ 帰無仮説が棄却できるか判断
結論を先に述べると…
切片以外の偏回帰係数について、$${p}$$ 値 $${\leq}$$ 有意水準 $${5\%}$$ です。
有意水準 $${5\%}$$ で帰無仮説は棄却され、切片以外の偏回帰係数の推定値は有意であると言えます。
「温度も圧力も配向度に影響を与えている」ですね。
よかったです！
テキスト p.65 図 2.6.2「温度の場合」、図 2.6.3「圧力の場合」の有意水準と棄却域の関係の図を描画しましょう。
### 温度と圧力の係数の検定 p.67 図2.6.2, 2.6.3
## 設定と準備
alpha = 0.05 # 有意水準（両側検定）
left, right = -4, 4 # グラフのx軸の両端
t_vals = result1_sm.tvalues.iloc[1:].values # t値の取得
p_vals = result1_sm.pvalues.iloc[1:].values # p値の取得
res = result1_sm # statsmodelsの回帰分析結果
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
# 係数ごとに描画を繰り返し処理
for t_val, p_val, col, ax in zip(t_vals, p_vals, VARS, axes.flat):
## 設定と準備
# 基本の色
color = 'tab:blue'
# 標本サイズn、説明変数の数pの取得
n, p = int(res.nobs), int(res.df_model)
# t分布の設定
t_dist = stats.t(df=n-p-1)
# 有意水準alphaの両側検定の棄却限界値
c_val = t_dist.isf(q=alpha/2)
# t分布の確率密度関数を算出するためのx軸の値
x_val1 = np.linspace(left, right) # 確率密度関数用
x_val2 = np.linspace(left, -c_val) # 下側棄却域の塗りつぶし用
x_val3 = np.linspace(c_val, right) # 上側棄却域の塗りつぶし用
## 描画
# t分布の確率密度関数の描画
ax.plot(x_val1, t_dist.pdf(x_val1), color=color, label='$t$分布')
# 下側棄却域の塗りつぶし描画
ax.fill_between(x_val2, t_dist.pdf(x_val2), color=color, alpha=0.2,
label='棄却域')
# 上側棄却域の塗りつぶし描画
ax.fill_between(x_val3, t_dist.pdf(x_val3), color=color, alpha=0.2)
# t値の垂直点線（赤）の描画
ax.axvline(t_val, color='tab:red', ls='--', label='$t$値')
# 棄却限界値の垂直点線（黒）の描画
ax.axvline(-c_val, color='grey', ls='--', label='棄却限界値')
ax.axvline(c_val, color='grey', ls='--')
# 修飾
ax.set_title(f'【{col}】の両側検定の可視化\n'
f'自由度{n-p-1}の$t$分布'''
f'\n棄却限界値:{c_val:.4f}, $t$ 値:{t_val:.4f} \n'
f'有意水準{alpha:.1%}, $p$ 値={p_val:.4f}')
ax.set_xticks([-c_val, 0, c_val])
ax.legend(loc='upper left')
plt.show()【実行結果】
青い曲線は自由度 $${7}$$ の $${t}$$ 分布の確率密度関数です。
曲線の両端の青い領域が有意水準 $${5\%}$$ の棄却域です。
青い領域の確率（面積）は $${0.05}$$ です。
赤い垂直点線が $${t}$$ 値です。
両方の説明変数の $${t}$$ 値は 青い領域に位置しており、帰無仮説は棄却されます。
🛸ちょっと寄り道：ChatGPTおすすめのチャート 🛸
もっともっと重回帰分析の可視化例を知りたくなり、ChatGPTに訊いてみました。
【プロンプト】
回帰係数の有意性の検定の可視化例としてオススメのグラフを教えてください
でてきた答えは…
1. フォレストプロット（Coefficient Plot）
2. ブートストラップ分布のヒストグラム
3. 部分回帰プロット（Added-Variable Plot）
4. ボルケーノプロット（Volcano Plot）
４番に取り組みましょう！
１番・２番は次回記事で紹介します。
３番は前回記事で紹介しています。
📈 ボルケーノプロット（Volcano Plot）
ChatGPTによる解説で、ざっくりボルケーノプロットの概要・読み取り方を見ておきましょう！
🌋 Volcano Plot の原理
X軸：効果量（偏回帰係数や標準化偏回帰係数）
変数がどれだけ目的変数に影響を与えているかを「数値の大きさ」で示します。
Y軸：$${−\log10}$$($${p}$$ 値)
$${p}$$ 値の小ささ（有意性の高さ）を、10を底にした対数で表します。
例えば $${p\ 値=0.01 \rightarrow −\log10(0.01)=2}$$、$${p\ 値=0.001 \rightarrow 3}$$ となり、直感的に「この点はとても有意」という高さがわかります。
🔍️ 読み取り方
右側に大きく離れていてかつ上のほうにある点
→ 正の効果が大きく、かつ統計的にも有意な変数左側に大きく離れていてかつ上のほうにある点
→ 負の効果が大きく、かつ有意な変数左右どちらにも離れていない or 下のほうに位置する点
→ 効果の大きさが小さいか、有意水準を満たしていない（$${p}$$ 値が大きい）可能性が高い変数
💡 ポイント
左右の広がり が「効果量の大きさ」を、上下の高さ が「$${p}$$ 値の小ささ（有意性）」を同時に伝えるので、変数ごとに「影響力」と「信頼度」の両面を一目で比較できます。
火山（Volcano）に見立てると、⛰️ の両脇に飛び出す変数が「強力な有意効果」を持つピークのように見える、というのが名前の由来です😊
この図を使えば、どの説明変数が「大きく」「確かな」影響を与えているかを直感的につかむことができますよ！
横軸（効果）を標準化偏回帰係数、縦軸（有意性）を $${−\log10}$$($${p}$$ 値) にしてボルケーノプロットを描画します。
### Volcano プロット
## 設定と準備
# 偏回帰係数
betas = result1_sm.params[1:]
# 説明変数の標準偏差、p値を取得
std_X = data1[VARS].std(ddof=1)
# 目的変数の標準偏差
std_y = data1[TARGET].std(ddof=1)
# 偏回帰係数のp値
pvals = result1_sm.pvalues[1:]
## 計算
# 標準化偏回帰係数を計算
beta_star = betas * std_X / std_y
# -log10(p値)の計算
neglogp = -np.log10(pvals)
## 描画
# 描画領域の設定
plt.figure(figsize=(6, 4))
# 標準化偏回帰係数と-log10(p値)の散布図
plt.scatter(beta_star, neglogp, s=70, color='tab:blue')
# 各点に変数名を表示
for var, x_val, y_val in zip(beta_star.index, beta_star.values, neglogp.values):
plt.text(x_val, y_val, var, fontsize=12, ha='left', va='bottom')
# 有意水準0.05の水平線の描画
plt.axhline(-np.log10(0.05), color='tab:red', ls='--', label='有意水準 5%')
# 修飾
plt.xlim(0.4, 0.6)
plt.ylim(1, 2)
plt.xlabel('標準化偏回帰係数 ($\\beta^*$)', fontsize=12)
plt.ylabel('-log10($p$値)', fontsize=12)
plt.title('Volcanoプロット (標準化偏回帰係数 vs -log10($p$値))')
plt.grid(True)
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()【実行結果】
説明変数「温度」「圧力」の両方が正の効果が大きく、かつ統計的にも有意な変数と言えそうです。
２変数を比べると、温度のほうが
【影響度】標準化偏回帰係数が大きい
⇒ 目的変数に対する影響が大きい【有意性】$${−\log10}$$($${p}$$ 値) が大きい
⇒ $${p}$$ 値が小さい
と言えるようです。
記事の最後はChatGPTが締めくくります。
今回は日々の成長にたとえて。
📘 ChatGPTのひとこと：
庭に種をまき、芽が出る土の健康をそっと確かめるように、今回はモデル全体の根付き具合（回帰の有意性）と、各成分が本当に効いているか（偏回帰係数の有意性）を統計的検定で優しくチェックしましたね🌱
次回は、その芽がどれくらい伸びるかを予測する「予測値の区間推定」と、成分ごとの働きぶりに確かな目印をつける「偏回帰係数の区間推定」という水やりの目安を学びます。
日々のデータという小さな庭で、芽の動きを感じながら、また一緒に静かな学びの時間を紡いでいきましょう😊✨
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！