「入門はじめての多変量解析」をPythonで写経 Vol.3 ～ ２章「はじめての重回帰分析」②当てはまりの良さ、標準化偏回帰係数

２章「はじめての重回帰分析」

書籍の著者　石村貞夫 先生、石村光資郎 先生

書籍「入門はじめての多変量解析」２章「はじめての重回帰分析」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

この記事は、モデルの評価 と 偏回帰係数の解釈 を取り扱います。
具体的には 決定係数、重相関係数、AIC、標準化偏回帰係数 などをテーマにします。

ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

２章 はじめての重回帰分析

この記事は２章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計
import pingouin as pg
import statsmodels.formula.api as smf
import statsmodels.api as sm

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
from graphviz import Digraph
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],
    '時間': [20, 20, 20, 20, 15, 20, 20, 20, 20, 15],
    }, index=range(1, 11))
data1.index.name = 'サンプルNo.'

# 結果の表示
data1

【実行結果】
データの個数（標本サイズ）は 10 です。

目的変数は「配向度」、説明変数の候補は「温度」「圧力」「時間」です。

■ Python ライブラリで重回帰分析結果を取得
重回帰分析を行えるライブラリ statsmodels、pingouin で、説明変数に「温度」「圧力」を用いて、重回帰分析を実施しておきます。
２つのライブラリは前回記事で扱いました。

🖲️statsmodels

# statsmodels利用
result1_sm = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1).fit()
result1_sm.summary()

【実行結果】
下段の「Intercept」（切片）・「温度」・「圧力」の coef が偏回帰係数です。

🖲️pingouin

# pingouin利用
result1_pg = pg.linear_regression(X=data1[['温度', '圧力']], y=data1['配向度'])
result1_pg

【実行結果】
coef 列が偏回帰係数です。

重回帰式の当てはまりの良さ

データから重回帰式を得たものの、重回帰式は目的変数の予測値をちゃんと推定できるのでしょうか・・・？

「重回帰式の当てはまりの良さ」「モデルの当てはまりの良さ」と呼ばれるものは、ざっくり、分析したデータに関連して「説明変数が目的変数を説明している度合い・指標」で「良さ」を確認します。
テキストは、決定係数とAICを中心に「当てはまりの良さを示す指標」を解説しています。

決定係数 p.50～

決定係数 $${R^2}$$ は、ざっくり「重回帰式が目的変数のばらつきを説明する度合い」的な指標です。
０～１の値をとり、説明度合いが高い＝当てはまりが良いほど、１に近づきます。

テキスト p.50、51 の決定係数の公式をお借りします。

📊 日本語バージョン

１行目の式は分母が「実測値」に関連する量、分子が「予測値」＝「重回帰式」に関連する量になっています。
実測値の何らかの全体に対して、重回帰式の何かで説明できることがイメージできるかと思います。

📊 数式バージョン
「平方和」の正体が明らかになります！

数式は目的変数で構成されています。

• $${y_i}$$：$${i}$$ 番目のデータの目的変数（実測値）

• $${\bar{y}}$$：目的変数の平均

• $${Y_i}$$：$${i}$$ 番目のデータの目的変数の予測値（重回帰式で予測）

各平方和の意味合いはざっくり、こんな感じです。

• 実測値の平方和（総平方和）が示すばらつきは「実測値と平均値の差」の二乗の合計であり、「目的変数のばらつき」（平均からの乖離度合い）です。

• 予測値の平方和（回帰の平方和）は、「重回帰式による予測値と平均値の差」の二乗の合計であり、「説明変数が説明する目的変数のばらつき」です。

• 残差平方和は「実測値と予測値の差＝残差」の二乗の合計であり、「説明変数が説明できていない目的変数のばらつき」です。

🔢 ３つの平方和の計算
テキストの表 2.4.1「３つの平方和」の計算をします。
最初にデータ点ごとの目的変数の「実測値」「予測値」「残差」を計算します。

### ３つの平方和 p.50 表2.4.1

# 実測値・予測値・残差の表の作成
data1_pred = pd.concat([data1['配向度'].rename('実測値'),
                        result1_sm.fittedvalues.rename('予測値'),
                        result1_sm.resid.rename('残差')], axis=1)
data1_pred

【実行結果】

実測値と予測値の乖離具合いを可視化して、直感的にイメージしましょう。

# 実測値と予測値のプロット

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 5))
# 実測値と予測値の散布図の描画
sns.scatterplot(data=data1_pred, x='実測値', y='予測値', s=70, ax=ax)
# 実測値＝予測値となる45度直線の描画
ax.plot([32, 60], [32, 60], color='tab:red', ls='--')
# 修飾
ax.set(xlim=(32, 60), ylim=(32, 60));

【実行結果】
４５度線（赤い点線）は実測値＝予測値を示しています。
この線にデータ点（青い点）が近いほど、予測値の精度が高いと言えます。
この図では、データ点が概ね４５度線付近に点在していますので、そこそこの予測精度のように思えます。

表 2.4.1 の各列の「平均」と３つの平方和を計算します。

# 平均値の表示
data1_pred_mean = data1_pred.mean(axis=0).rename('平均値').to_frame()
data1_pred_mean.round(10).T

【実行結果】
実測値・予測値の平均は「実測値の平均」です。
残差の平均は０です。

この平均値を用いて、３つの平方和を計算します。

# 平方和の算出 

# 上記の平均値を計算に使いやすいように加工
means = data1_pred_mean.values.flatten()

# 平方和の算出
sum_square1 = ((data1_pred.apply(lambda x: x - means, axis=1)**2)
                 .sum(axis=0)
                 .rename('平方和').to_frame().T)
sum_square1.round(2)

【実行結果】
テキスト掲載のとおり、
「実測値の平方和＝予測値の平方和＋残差平方和」です。

🔢 決定係数の計算
実測値の平方和と予測値の平方和を使って決定係数を計算します。

# 決定係数の算出 p.51
print((sum_square1['予測値'] / sum_square1['実測値']).values[0])

【実行結果】
決定係数は $${0.858}$$ です。
この重回帰式の当てはまりは良さそうです。

statsmodels と pingouin の重回帰分析の結果から決定係数を取り出してみましょう。

🖲️statsmodels

# statsmodels利用
result1_sm.rsquared

【実行結果】

🖲️pingouin

# pingouin利用
result1_pg.r2[0]

【実行結果】

自由度調整済み決定係数 p.52

決定係数は説明変数を増やすと値が高くなる性質があります。
説明変数の数の影響を排除できるのが「自由度調整済み決定係数 」$${\widehat{R}^2}$$ です。

テキストの自由度調整済み決定係数 $${\widehat{R}^2}$$ の公式をお借りします。

$${N}$$ はデータの個数（標本サイズ）、$${p}$$ は説明変数の数です。

🔢 自由度調整済み決定係数の計算
さきほど計算した 実測値の平方和（総平方和）$${S_T}$$ と残差平方和 $${S_E}$$ を用いて、自由度調整済み決定係数を計算しましょう。

# 自由度調整済み決定係数の算出 p.53

# 設定と準備
p = 2           # 説明変数の個数
N = len(data1)  # 標本サイズ

# SSR, SSE, SSTの算出
SSR = sum_square1['予測値']
SSE = sum_square1['残差']
SST = SSR + SSE

# 自由度調整済み決定係数の算出
R2_adj = 1 - (SSE/(N - p - 1)) / (SST/(N - 1))
R2_adj.values[0]

【実行結果】
決定係数 $${0.858}$$ より小さい値になりました。

statsmodels と pingouin の２変数モデルの自由度調整済み決定係数を取り出してみましょう。

🖲️statsmodels

# statsmodels利用
result1_sm.rsquared_adj

【実行結果】

🖲️pingouin

# pingouin利用
result1_pg.adj_r2[0]

【実行結果】

🔢 説明変数の数が異なる重回帰分析を比較
テキスト p.53 表 2.4.2「SPSSによる結果」のように、説明変数を
・「温度」「圧力」の２変数
・「温度」「圧力」「時間」の３変数
の２つのモデルで決定係数を比較してみましょう。

# 独立変数に時間を加えない場合と加える場合の決定係数 p.53 表2.4.2

# 独立変数に加える場合の回帰分析 ※statsmodels利用
result1_3_sm = smf.ols(formula='配向度 ~ 温度 + 圧力 + 時間', data=data1).fit()

# 表の作成
data1_r2_df = pd.DataFrame(
    {'温度 圧力': [result1_sm.rsquared, result1_sm.rsquared_adj],
     '温度 圧力 時間': [result1_3_sm.rsquared, result1_3_sm.rsquared_adj]},
     index=['決定係数', '自由度調整済み決定係数'])
data1_r2_df.loc['差'] = data1_r2_df.iloc[0, :] - data1_r2_df.iloc[1, :]
data1_r2_df.columns.name = '説明変数'

# 結果の表示
data1_r2_df.round(3)

【実行結果】
決定係数は、説明変数の数が多い「温度・圧力・時間」のほうが大きな値です。
一方で、自由度調整済み決定係数は、説明変数の数が少ない「温度・圧力」のほうが大きな値です。

重相関係数 p.54

重相関係数 $${R}$$ は目的変数の実測値と予測値の相関係数です。
重相関係数の二乗は決定係数 $${R^2}$$ です。

Excel の分析ツール「回帰分析」の実行結果に含まれる「重相関R」が重相関係数です。

🔢 重相関係数の計算
２変数モデルの重相関係数を計算しましょう。
pandas の corr() を用いて、実測値と予測値の相関係数を算出します。

### 重相関係数 p.54
multi_corr1 = data1_pred[['実測値', '予測値']].corr().iloc[0, 1]
multi_corr1

【実行結果】
重相関係数は 0.926 です。

重相関係数を二乗すると決定係数に一致することを確認しましょう。

### 重相関係数の二乗が決定係数
multi_corr1**2

【実行結果】
先に計算した決定係数 0.858 と一致しました。

AIC（赤池情報量規準） p.55

AICは 値が小さいほど 予測の良いモデルであることを表します。
なおテキストは、AICは当てはまりの悪さを示す統計量、と紹介しています。

テキストの重回帰モデルの場合のAICの公式をお借りします。

$${N}$$ はデータの個数（標本サイズ）、$${p}$$ は説明変数の数です。

🔢 AICの計算
２変数モデルのAICを計算しましょう。

### AICの算出 p.55

# 設定と準備
N = len(data1)  # 標本サイズ
p = 2           # 説明変数の個数

# 残差平方和の算出
SSE = sum((data1_pred['実測値'] - data1_pred['予測値'])**2)

# AICの算出
AIC1 = N * (np.log(2 * np.pi * (SSE/N)) + 1) + 2 * (p + 2)
AIC1

【実行結果】

statsmodels の２変数モデル・３変数モデルのAICを取り出してみましょう。

🖲️statsmodels
こちらは２変数モデルです。

# statsmodels利用 説明変数が２つ
result1_sm.aic

【実行結果】

テキストの公式による値と異なるのは、テキストとは異なる計算式を採用しているからです。
statsmodels 公式によると切片を含むモデルの場合、次式で計算します。

statsmodelsの方が一般的な気がします（個人の見解です）。

３変数モデルのAICです。

# statsmodels利用 説明変数が３つ
result1_3_sm.aic

【実行結果】

AICの値が小さい２変数モデルのほうが当てはまりが良いと評価できます。

（おまけ）分散分析表
テキスト p.55 の末尾に佇む「分散分析表」。
statsmodels の２変数モデルの回帰分析結果を用いて、分散分析表を作ってみましょう。

# 重回帰の分散分析表の例 statsmodels利用
sm.stats.anova_lm(result1_sm)

【実行結果】
回帰による変動は説明変数別に分かれて表示されます。
残差による変動は Residual 行に示されています。

詳しくは次回記事で！

偏回帰係数の意味すること p.56

重回帰式で各説明変数の係数「偏回帰係数」はどんな意味を持つのでしょう・・・？

この意味合いを実際の計算で確かめましょう！

■ 単回帰分析の回帰係数と重回帰分析の偏回帰係数は異なる
テキスト p.56~ の「パス図」を描いて、単回帰分析の回帰係数と重回帰分析の偏回帰係数が異なることを確認していきましょう。
準備として、statsmodels で単回帰分析・重回帰分析を実行しておきます。

### 重回帰分析のパス図 p.56 図2.5.1～図2.5.3

## 係数の推定 by statsmodels

# 説明変数: 温度
result1_sm_temp = smf.ols(formula='配向度 ~ 温度', data=data1).fit()
b1_temp = result1_sm_temp.params.iloc[1]

# 説明変数: 圧力
result1_sm_press = smf.ols(formula='配向度 ~ 圧力', data=data1).fit()
b1_press = result1_sm_press.params.iloc[1]

# 説明変数: 温度、圧力
result1_sm_temp_press = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1).fit()
b1, b2 = result1_sm_temp_press.params.iloc[1:]

【実行結果】なし

目的変数と２つの説明変数の重回帰分析のパス図を描きます。
graphviz ライブラリの 有向グラフ Digraph を利用します。

## 重回帰分析のパス図 図2.5.1

## 設定
# 有向グラフオブジェクトの生成、neatoでnodeの位置調整を実施
g = Digraph(engine='neato')
# nodeの基本属性の設定
g.attr('node', shape='box', fontname='Meiryo UI')

## node：頂点の作成、posで位置固定
g.node('温度x1', pos='0, 1!')
g.node('圧力x2', pos='0, 0!')
g.node('配向度y', pos='2, 0.5!')

## edge：辺の作成
g.edge('温度x1', '配向度y', label=f'b1={b1:.3f}')
g.edge('圧力x2', '配向度y', label=f'b2={b2:.3f}')

## グラフの表示
g

【実行結果】
矢印付近の等式は偏回帰係数です。

続いて温度のみを説明変数にする単回帰分析のパス図です。

## 単回帰分析のパス図 図2.5.2

## 設定
# 有向グラフオブジェクトの生成、neatoでnodeの位置調整を実施
g = Digraph(engine='neato')
# nodeの基本属性の設定
g.attr('node', shape='box', fontname='Meiryo UI')

## node：頂点の作成、posで位置固定
g.node('温度x1', pos='0, 0!')
g.node('配向度y', pos='2, 0!')

## edge：辺の作成
g.edge('温度x1', '配向度y', label=f'b={b1_temp:.3f}')

## グラフの表示
g

【実行結果】
重回帰分析の偏回帰係数 $${b_1}$$ と異なる値になっています。

圧力のみを説明変数にする単回帰分析のパス図です。

## 単回帰分析のパス図 図2.5.3

## 設定
# 有向グラフオブジェクトの生成、neatoでnodeの位置調整を実施
g = Digraph(engine='neato')
# nodeの基本属性の設定
g.attr('node', shape='box', fontname='Meiryo UI')

## node：頂点の作成、posで位置固定
g.node('圧力x2', pos='0, 0!')
g.node('配向度y', pos='2, 0!')

## edge：辺の作成
g.edge('圧力x2', '配向度y', label=f'b={b1_press:.3f}')

## グラフの表示
g

【実行結果】
こちらも重回帰分析の偏回帰係数 $${b_2}$$ と異なる値になっています。

単回帰分析と重回帰分析の係数の違いについて、テキストは「説明変数どうしが何らかの影響を及ぼし合っているから」としています。
温度と圧力の相関係数を確認します。

### 説明変数間の相関係数 p.57 ※pandas利用
data1.corr().loc['温度', '圧力']

【実行結果】
正の相関がありました！

■ 単回帰分析の回帰係数と重回帰分析の偏回帰係数をつなぐ
テキストに沿って、温度の回帰係数・偏回帰係数を通じて、単回帰分析の回帰係数と重回帰分析の偏回帰係数のつながりを確認します。

圧力から温度への影響を単回帰分析で確認します。

### 温度と圧力の単回帰式 p.57 ※statsmodels利用

result1_sm_press2temp = smf.ols(formula='温度 ~ 圧力', data=data1).fit()
b0, b1 = result1_sm_press2temp.params
print(f'温度 = {b1:.4f} x 圧力 + {b0:.4f}')

【実行結果】
回帰係数は $${0.1029}$$ です。

温度から圧力の影響を取り除くことは、温度・圧力の単回帰分析の結果の残差（実測値と予測値の差）を求めることです。
やってみましょう。

### 温度と圧力の単回帰分析における残差V p.58 図2.5.1

# データフレームの作成
data1_V = data1[['温度', '圧力']].copy()
data1_V['予測値'] = result1_sm_press2temp.fittedvalues
data1_V['残差V'] = data1_V['温度'] - data1_V['予測値']

# 結果の表示
data1_V

【実行結果】
残差 $${V}$$ は温度の実測値から予測値を差し引いた値です。
この残差が「圧力の影響を取り除いた温度」になります。

同じように目的変数の配向度から圧力の影響を取り除きましょう。

### 配向度と圧力の単回帰分析における残差W p.58 図2.5.1

# データフレームの作成
data1_W = data1[['配向度', '圧力']].copy()
data1_W['予測値'] = result1_sm_press.fittedvalues
data1_W['残差W'] = data1_W['配向度'] - data1_W['予測値']

# 結果の表示
data1_W

【実行結果】
残差 $${W}$$ は配向度の実測値から予測値を差し引いた値です。
この残差が「圧力の影響を取り除いた配向度」になります。

最後に２つの残差の単回帰分析＝圧力の影響を取り除いた配向度と温度の単回帰分析を行って、係数のつながりを確認します。

### 残差Wと残差Vの単回帰式の回帰係数の算出 p.59 ※statsmodels利用

# 残差Wと残差Vのデータフレームの作成: statsmodels用
result1_WV_df = pd.concat([data1_W['残差W'], data1_V['残差V']], axis=1)
# 単回帰分析の実行
result1_sm_V2W = smf.ols(formula='残差W ~ 残差V', data=result1_WV_df).fit()
# 係数の推定値を取得して表示
b0, b1 = result1_sm_V2W.params
print(f'残差W = {b1:.3f} x 残差V + {b0:.3f}')

【実行結果】
残差 $${V}$$ ＝温度（圧力の影響控除後）の回帰係数 $${3.470}$$ は、２変数の重回帰分析に関する温度の偏回帰係数 $${3.470}$$ と一致しました！

テキストのとおり、「圧力の影響を取り除いたあとの温度から配向度への影響の程度を示す偏回帰係数 $${3.470}$$」は「温度から配向度への単回帰分析の回帰係数 $${3.470}$$」に等しい、ことが分かりました。

🛸ちょっと寄り道：部分回帰プロット 🛸

ここまで圧力の影響を取り除いた「温度と配向度」の関係を計算にて確認しました。
もっとカジュアルにチャートで可視化したいですよね？
「部分回帰プロット」で可視化できます！
１つの説明変数と目的変数の関係を「残りの説明変数の影響を取り除いて」見ることができます。
実はこのプロットの存在を ChatGPT に教えてもらいました。

statsmodels の plot_partregress を簡単に描画できます。

### 部分回帰プロット

# 追加インポート
from statsmodels.graphics.regressionplots import plot_partregress

# 変数の設定
target = '配向度'        # 目的変数名
vars = ['温度', '圧力']  # 説明変数名

# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 説明変数ごとに部分回帰プロット描画を繰り返し処理
for var, ax in zip(vars, axes.flat):
    # 部分回帰プロットの描画
    plot_partregress(
        endog=target,   # 目的変数
        exog_i=var,     # 説明変数
        exog_others=[exc_v for exc_v in vars if exc_v != var], # 除外する説明変数
        data=data1,     # データフレーム
        obs_labels=False,  # 散布図にインデックスを付記するかどうか 
        ax=ax,
    )
    # 修飾
    ax.set_title(f'{var}の部分回帰プロット')
    ax.set_xlabel(f'残差化した{var}', fontsize=12)
    ax.set_ylabel(f'残差化した{target}', fontsize=12)
plt.show()

【実行結果】
左の直線の傾きは温度の偏回帰係数です。
右の直線の傾きは圧力の偏回帰係数です。

statsmodels の ols による回帰分析結果を用いて「もっとシンプルなコード」で部分回帰プロットを描けるコードはこちらです。

### 部分回帰プロット グリッド版

# 追加インポート
from statsmodels.graphics.regressionplots import plot_partregress_grid

## 描画
# 描画領域の設定
fig = plt.figure(figsize=(8, 4))
# 部分回帰プロットの描画
plot_partregress_grid(result1_sm, exog_idx=vars, fig=fig)
# 修飾
plt.tight_layout();

【実行結果】

部分回帰プロットの読み取り方をChatGPTに訊きました！

１．部分回帰プロットの概要

• 「この変数だけがどれだけ効いているか」を見るためのグラフです。

• 手順はざっくり３ステップ：

  1. 他の変数の影響を取り除く

     • 目的変数 y から「残りの説明変数」で説明できる部分を引いて残差を出す

     • 対象の説明変数 X₁ からも同じく「残りの説明変数」で説明できる部分を引いて残差を出す

  2. 残った“残差同士”をプロット

     • 横軸に「X₁ の残差」、縦軸に「y の残差」を並べる

  3. 直線をあてはめ

     • その直線の傾きが、他の変数を固定したときの X₁ の回帰係数にほぼ一致します。

２．読み取り方

• 散らばり具合が小さく直線に沿っていれば
  → 他の変数の影響を除いても、X₁ は y をしっかり説明している

• 点がバラバラだったり外れ値が多いと
  → X₁ 単独では説明力が弱かったり、外れ値の影響が大きい可能性あり

• 直線の傾きの向きと大きさ
  → 「＋」なら X₁ が増えるほど y も増える、「－」なら逆向き
  → 傾きの大きさが、１単位あたりの変化量を表します

標準化偏回帰係数 p.60

配向度に対して、温度の影響と圧力の影響のどちらの影響が大きいのでしょう・・・？
重回帰分析の偏回帰係数は説明変数どうしの単位が異なるため、単純に比較できません。

そこで登場するのが「標準化偏回帰係数」です。
データの標準化済みの目的変数・説明変数を用いた重回帰式の偏回帰係数です。

なおテキストは「標準偏回帰係数」の用語を用いていますが、この記事は個人的に使い慣れている「標準化偏回帰係数」を用いています。

データの標準化は、データからデータの平均を差し引いて、データの標準偏差で割る操作のことです。
標準化後のデータは平均０、分散１になります。

🔢 標準化偏回帰係数の計算
配向度・温度・圧力に対してデータの標準化を行い、標準化偏回帰係数を計算、温度と圧力の影響の大きさを比べてみましょう。

まずはデータの標準化です。
テキスト p.61 表 2.5.7「データの標準化」に相当します。

### データの標準化 p.62 表2.5.7
data1_std = (data1[['配向度', '温度', '圧力']]
             .apply(lambda x: (x - x.mean()) / x.std(ddof=1), axis=0))
data1_std

【実行結果】
各変数は「無単位」になりました。

標準化データを対象にして重回帰分析を実行します。
テキスト p.62 表 2.5.8 に相当します。
statsmodels を利用します。

### 重回帰分析の実行 標準化偏回帰係数の算出 p.62 表2.5.8 ※statsmodels利用
result1_std_sm = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1_std).fit()
result1_std_sm.summary2().tables[1].round(4)

【実行結果】
Coef. が標準化偏回帰係数です。
標準化偏回帰係数を比べると、温度は $${0.5576}$$、圧力は $${0.4836}$$。
温度のほうが配向度に与える影響が大きいです！

🔢 相関係数を用いた標準化偏回帰係数の計算
テキストの数式をお借りします。

📊 分散共分散行列を用いた偏回帰係数の公式

📊 相関行列を用いた標準化偏回帰係数の公式
$${b_1^*, b_2^*}$$ が標準化偏回帰係数です。

相関行列を用いた標準化偏回帰係数の公式を使って計算します。

### 相関行列を用いた標準化偏回帰係数の算出 p.63

# x1とy, x2とyの相関係数の算出
corr_vec = data1.corr().loc['配向度', ['温度', '圧力']].values
print('x1とy, x2とyの相関係数:')
print(corr_vec, '\n')

# x1とx2の相関行列の算出
corr_mtx = data1[['温度', '圧力']].corr().values
print('x1とx2の相関行列:')
print(corr_mtx, '\n')

# 標準化偏回帰係数の算出
b1, b2 = np.linalg.inv(corr_mtx) @ corr_vec
print('標準化偏回帰係数:')
print(f'b1: {b1:.4f},  b2: {b2:.4f}')

【実行結果】
標準化偏回帰係数を算出できました！

🔢 偏回帰係数を使った標準化偏回帰係数の計算
データの標準化をしなくてもＯＫ🙆‍♂️な方法です。

statsmodels の偏回帰係数（標準化していないデータ）を標準化偏回帰係数に変換しましょう。

### 標準化しないデータの偏回帰係数を標準化偏回帰係数に変換

# 説明変数Xと目的変数yの標準偏差を算出
std_X = data1[['温度', '圧力']].std(ddof=1)
std_y = data1['配向度'].std(ddof=1)

# 標準化偏回帰係数 ＝ 通常の偏回帰係数 × Xの標準偏差 ÷ yの標準偏差
(result1_sm.params[1:] * std_X / std_y).rename('標準化偏回帰係数').to_frame()

【実行結果】
標準化偏回帰係数に変身しました！
これでデータの標準化をしなくても標準化偏回帰係数をＧＥＴできますね！

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

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