「入門はじめての多変量解析」をPythonで写経 Vol.2 ～ ２章「はじめての重回帰分析」①重回帰式、重回帰モデル、偏回帰係数
２章「はじめての重回帰分析」
書籍の著者 石村貞夫 先生、石村光資郎 先生
書籍「入門はじめての多変量解析」２章「はじめての重回帰分析」の Python写経活動記録 です。
多変量解析の入門を Python と一緒に学ぶ写経シリーズです。
２章は多変量解析の王様（弊社調べ）の「重回帰分析」を堪能する章です。
この記事は、重回帰モデル・重回帰式 と 偏回帰係数の算出 を取り扱います。
ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀
はじめに
このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。
書籍の紹介と引用表記はリンク先の記事に掲載しています。
２章 はじめての重回帰分析
この記事は２章の以下のSectionを取り扱います。
2.1 重回帰分析とは因果なもの？！
2.2 重回帰式と重回帰モデル
2.3 重回帰式を求めよう
記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。
この記事で用いるライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
import pandas as pd
import sympy
# 統計
import pingouin as pg
import statsmodels.formula.api as smf
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
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 です。
このデータは「セラミックス」の製造に関するもののようです。
セラミックスの知識（ドメイン知識）が無いので、ChatGPTに「配向度・温度・圧力の関係」を訊いてみました。
配向度の概要
配向度は、小さな粒や板がどれだけ同じ方向を向いてそろっているかを、０（バラバラ）～１（ピタリとそろう）の数値で示したものです。
配向度でわかること
値が低いほど、あらゆる方向で同じ性質を持つ「均一な」材料
値が高いほど、一方向にだけ強い・熱や電気が流れやすい「向き依存の」材料
配向度・温度・圧力の関係
温度が高いと 粒がやわらかく動きやすくなり、向きを合わせやすくなる
圧力をかけると 粒が押されて同じ向きに並びやすくなる
高温＋高圧 の組み合わせで、よりそろった状態になりやすい（ただし温度が高すぎると粒が大きくなりすぎることも）
テキストの配向度の単位は ％ かもですね！
配向度・温度・圧力の関係を見ると次のような関係がありそうです。
・温度が高い ⇒ 配向度が高い
・圧力が高い ⇒ 配向度が高い
散布図を描いて、配向度と温度・圧力の関係を直感的に確認しましょう。
ヒストグラムを付記できる seaborn の pairplot() を利用します。
# ヒストグラムと散布図の描画
sns.pairplot(
data=data1,
height=2, # グラフ１つの高さ
diag_kws={'ec': 'white', 'alpha': 0.7}, # ヒストグラムの引数
plot_kws={'s': 50} # 散布図の引数
);【実行結果】
配向度の行を見ましょう。
温度と配向度、圧力と配向度は右上がりの傾向が見られます。
正の相関がありそうです。
相関係数を確認しましょう。
pandas の corr() を利用します。
# 相関係数の算出
data1.corr().round(3)【実行結果】
配向度の行を見ましょう。
温度と配向度、圧力と配向度の相関係数は $${0.8}$$ 超です。
強めの正の相関がありました。
重回帰式と重回帰モデル
ChatGPTに訊いた配向度・温度・圧力の関係や、データの散布図から読み取った相関関係から、「温度や圧力を使って配向度を見積もれそう」な気がします！
例えばこんな感じの１次式の関係です。
$$
配向度 = ▢ \times 温度 + △ \times 圧力 + 残りの部分
$$
等式の左辺を「目的変数」と呼びましょう。配向度です。
重回帰分析の目的変数は「量的変数」（数値データ）です。
等式の右辺の温度・圧力は「説明変数」です。
説明変数の係数 ▢、△ を解明したいですよね！
これが重回帰分析の動機になります。
（注）
テキストでは、目的変数を「従属変数」、説明変数を「独立変数」と呼んでいます。
ここでテキスト p.36 図 2.2.1 のグラフを見ておきます。
seaborn の regplot で単回帰直線付きの散布図を描画します。
### 散布図の描画 p.36 図2.2.1
# 目的変数の変数名の設定
target = '配向度'
# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
# 説明変数ごとに単回帰直線付き散布図の描画を繰り返し処理
for col, ax in zip(data1.columns[1:3], ax):
# 説明変数の１つと目的変数の散布図（単回帰直線付き）を描画
sns.regplot(data=data1, x=col, y=target, ci=None, ax=ax,
line_kws={'color': 'tab:red'}, scatter_kws={'s': 70})【実行結果】
温度と配向度の散布図、圧力と配向度の散布図に、単回帰直線（赤い線）を重ねています。
データ点（青い点）は直線上に綺麗に乗っているわけではなく、直線の周辺にバラバラと存在しています。
テキストはこのバラバラと存在する様子を「データは１次式の直線からずれているので、１次式の等号がそのままでは成り立たない」と説明しています。
この流れを受けつつ、重回帰モデル と 重回帰式 を区別して、重回帰分析に進みます。
■ 重回帰モデル p.37
テキストは、直線とのずれを考慮した次の式を「重回帰モデル」と呼んでいます。
$$
y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots \beta_p x_{ip} + \varepsilon_i
$$
$${N}$$ 個のデータ $${i}$$ について、$${y_i}$$ は目的変数、$${x_{i1}}$$ 等は $${p}$$ 個の説明変数、$${\beta_1}$$ 等は $${p+1}$$ 個の母偏回帰係数、$${\varepsilon_i}$$ はデータ $${i}$$ ごとの誤差（ずれ）です。
配向度のデータを重回帰モデルで示してみます。
$$
配向度_i = \beta_0 + \beta_1 \times 温度_i + \beta_2 \times 圧力_i + \varepsilon_i
$$
テキストは、重回帰モデルの誤差 $${\varepsilon_i}$$ に対して、次を仮定するとしています。
・$${\varepsilon_i}$$ は 平均 $${0}$$、分散 $${\sigma^2}$$ の正規分布に従う
・$${\varepsilon_i, \varepsilon_j}$$ は互いに独立である（$${\text{Cov}(\varepsilon_i, \varepsilon_j)=0}$$）
■ 重回帰式 p.37
目的変数の予測値を算出する式です。
誤差（ずれ）を示す項がありません。
$$
Y = b_0 + b_1 x_1 + b_2 x_2 + \cdots + b_p x_p
$$
$${Y}$$ は目的変数の予測値、$${x_{1}}$$ 等は $${p}$$ 個の説明変数、$${b_1}$$ 等は $${p+1}$$ 個の偏回帰係数（母偏回帰係数の推定量）です。
配向度のデータを重回帰式で示してみます。
$$
配向度 = b_0 + b_1 \times 温度 + b_2 \times 圧力
$$
重回帰式を可視化しましょう。
テキスト p.37 図 2.2.2「重回帰式のグラフ」を描きます。
説明変数が２つなので、重回帰式のチャートは「平面」になります。
※説明変数が１つの単回帰式は直線になっていましたね！
目的変数を加えた３つの変数を扱うので３次元グラフになります。
### 重回帰式のグラフ p.37 図2.2.2
## 設定
N = 100 # 2次元平面の縦横の格子の数
b0, b1, b2 = 10, -0.7, -0.1 # 重回帰式の係数 y = b1*x1 + b2*x2 + b0
## 予測値Yの平面の設定
# x1,x2の値の設定
x = np.linspace(-5, 5, N) # x軸・y軸の値の設定
x1, x2 = np.meshgrid(x, x) # 2次元平面の格子データの作成
# Yの垂直方向の値の算出
Y = b1*x1 + b2*x2 + b0 # 重回帰式でYの値を算出
## データ点の設定
# 実測値の設定
x1_obs, x2_obs, y_obs = 2, 0, 17
# 予測値の算出
y_pred = b1*x1_obs + b2*x2_obs + b0
## 描画
# 描画領域の設定 ※3次元指定
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
# 3次元上のYの平面を描画
ax.plot_surface(x1, x2, Y, color='tab:green', alpha=0.3, lw=0, label='重回帰式')
# ３つのデータ点の描画
# 実測値のデータ点（赤）
ax.scatter(x1_obs, x2_obs, y_obs, color='tab:red', s=70, label='実測値')
# 予測値のデータ点（緑）
ax.scatter(x1_obs, x2_obs, y_pred, color='tab:green', s=70, label='予測値')
# 実測値のxy平面上のデータ点（グレイ）
ax.scatter(x1_obs, x2_obs, 0, color='gray', s=70)
# 実測値からxy平面の垂直点線の描画
ax.plot([x1_obs, x1_obs], [x2_obs, x2_obs], [y_obs, 0], color='gray', ls='--')
# 修飾
ax.set_xlabel('$x_1$', labelpad=1, fontsize=12)
ax.set_ylabel('$x_2$', labelpad=1, fontsize=12)
ax.set_zlabel('$Y$', labelpad=-1, fontsize=12)
ax.set_zlim(0, 18)
ax.set_zticks([0, 5, 10, 15])
ax.view_init(elev=20, azim=-70) # 視点の角度(視点の高さ、平面の回転)
ax.legend(bbox_to_anchor=(1.2, 0.8))
plt.tight_layout()
plt.show()【実行結果】
薄緑色の四角形が重回帰式の目的変数 $${Y}$$の平面です。
目的変数の実測値と予測値の例も描いています。
目的変数の予測値の点（緑）は回帰平面上に存在し、実測値の点（赤）は回帰平面から外れて存在します。
予測値と実測値の差を「残差」と呼びます。
重回帰式を求める１「最小二乗法計算」編 p.38～43
目的変数の実測値と予測値の差「残差」が最小になるように偏回帰係数を推定する「最小二乗法」を使います。
計算手順は以下のようになります。
データの残差の二乗和（残差平方和）を算出する
残差平方和を偏回帰係数で偏微分する
偏微分の結果を０と置いて、偏回帰係数について連立方程式を解く
紙とペンで解く方法はテキストをご覧ください。
この記事では sympy ライブラリを用いて、偏微分と連立方程式に取り組みます。
① データの残差の二乗和（残差平方和）を算出する
残差平方和は次の式で求められます。
$$
Q = \sum_{i=1}^N \big( y_i - (b_0 + b_1 x_{i1} + b_2 x_{i2} + \cdots + b_p x_{ip}) \big)^2
$$
説明変数 $${x_1, x_2}$$ と 偏回帰係数 $${b_0, b_1, b_2}$$ を用いて残差平方和の式を作ります。
Q = … の部分で、Python のリスト内包表記を用いて、残差の二乗値をリストに展開して、sum で残差二乗値を足し上げています。
### 重回帰式の求め方（その１） p.38～ sympy利用
# 変数の定義
x1, x2, b0, b1, b2 = sp.symbols('x1 x2 b0 b1 b2')
# 残差平方和の算出
Q = sum([(y - (b1*x1 + b2*x2 + b0))**2 for y, x1, x2 in data1.values])
print('【残差平方和】')
display(Math(f'Q = {sp.latex(Q.expand())}'))【実行結果】
テキスト p.40 の残差平方和 $${Q(b_1, b_2, b_0)}$$ と一致しています。
② 残差平方和を偏回帰係数で偏微分する
diff で 残差平方和 Q を 偏回帰係数 b0、b1、b2 で偏微分します。
# 残差平方和をb0, b1, b2で偏微分
expr_b0 = sp.diff(Q, b0)
expr_b1 = sp.diff(Q, b1)
expr_b2 = sp.diff(Q, b2)
# 偏微分結果を表示
print('【残差平方和を偏微分】')
display(Math('\cfrac{\partial Q}{\partial b_0}= ' + f'{sp.latex(expr_b0)}'))
display(Math('\cfrac{\partial Q}{\partial b_1}= ' + f'{sp.latex(expr_b1)}'))
display(Math('\cfrac{\partial Q}{\partial b_2}= ' + f'{sp.latex(expr_b2)}'))【実行結果】
テキスト p.42 の 式 (1) ～ (3) と一致しています。
③ 偏微分の結果を０と置いて、偏回帰係数について連立方程式を解く
solve() で連立方程式を求解します。
引数に偏微分で得た３つの方程式 expr_b0、expr_b1、expr_b2 をリストにして与えます。
# 偏微分を０とおいて、連立方程式を解く
solved = sp.solve([expr_b0, expr_b1, expr_b2])
# 解を表示
print('【偏回帰係数の推定値：連立方程式の解】')
display(Math(f'{sp.latex(solved)}'))【実行結果】
テキスト p.42 の $${b_0, b_1, b_2}$$ の値と一致しています。
偏回帰係数を用いて、テキスト p.42 にように重回帰式で表現しましょう。
# 重回帰式
display(Math(f'Y = {sp.latex(solved[b1] * x1 + solved[b2] * x2 + solved[b0])}'))【実行結果】
sympy であっさり解けました。
display() の式の表示形式を整える部分は少々複雑ですが、残差平方和算出、偏微分計算、連立方程式求解は、１式１行でコードを書けています。
シンプルでパワフルな sympy ライブラリをぜひご活用ください！
せっかく偏回帰係数を推定して重回帰式を完成できたので・・・
重回帰式の読み取りを行いましょう。
温度 $${x_1}$$ が１℃大きくなると、配向度 $${Y}$$ が $${3.470}$$ 大きくなる関係
圧力 $${x_2}$$ が１単位大きくなると、配向度 $${Y}$$ が $${0.533}$$ 大きくなる関係
と読み取りできそうですね！
３次元のチャートで重回帰式の「平面」を見てみましょう。
### 重回帰式のグラフ
## 設定
N = 100 # 2次元平面の縦横の格子の数
## 偏回帰係数の設定
b0_plot = float(solved[b0])
b1_plot = float(solved[b1])
b2_plot = float(solved[b2])
## 予測値Yの平面の設定
# x1,x2の値の設定
x1 = np.linspace(15, 20, N) # x軸・y軸の値の設定
x2 = np.linspace(10, 50, N) # x軸・y軸の値の設定
x1, x2 = np.meshgrid(x1, x2) # 2次元平面の格子データの作成
# Yの垂直方向の値の算出
Y = b1_plot*x1 + b2_plot*x2 + b0_plot # 重回帰式でYの値を算出
## 描画
# 描画領域の設定 ※3次元指定
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
# 3次元上のYの平面を描画
ax.plot_surface(x1, x2, Y, color='tab:green', alpha=0.3, lw=0, label='重回帰式')
# 修飾
ax.set_xlabel('温度', labelpad=1, rotation=-10, fontsize=12)
ax.set_ylabel('圧力', labelpad=1, rotation=70, fontsize=12)
ax.set_zlabel('配向度', labelpad=3, rotation=90, fontsize=12)
ax.view_init(elev=30, azim=-70) # 視点の角度(視点の高さ、平面の回転)
plt.tight_layout()
plt.show()【実行結果】
薄緑色の平面が重回帰式です。
温度が大きくなるごとに、圧力が大きくなるごとに、配向度の値も大きくなっています。
重回帰式（偏回帰係数）の計算方法にはさまざまあります。
以降で計算方法を確認していきましょう。
重回帰式を求める２「平方和積和行列」編 p.44～45
テキストは別法「平方和積和行列」を利用する方法を掲載しています。
計算手順は以下のようになります。
平方和積和行列を算出する
平方和積和行列から連立方程式を構成する
逆行列を用いて連立方程式を解く
Python で解きましょう。
① 平方和積和行列を算出する
データの偏差行列を作成し、偏差行列の積をとると平方和積和行列になります。
平方和積和行列の「平方和」は行列の対角成分であり、同一の変数の偏差平方和です。
「積和」は対角成分以外の成分であり、２変数の偏差の積の和です。
最初に偏差行列を計算します。
データの中心化とも呼ばれ、データからそのデータの平均を差し引いた値です。
### 平方和積和行列を用いた偏回帰係数の算出 p.44~45
# 偏差の算出
deviation = data1.apply(lambda x: x - x.mean(), axis=0)
deviation【実行結果】
続いて平方和積和行列を作成します。
「偏差行列$${^{\top}}$$ @ 偏差行列」で計算できます（＠は行列積）。
# 平方和積和行列の算出
sum_prod_dev_mtx = deviation.T @ deviation
sum_prod_dev_mtx【実行結果】
② 平方和積和行列から連立方程式を構成する p.45
平方和積和行列の２行目以降で「連立方程式」を構成します。
$$
\begin{cases}
83.5 = 16.0 b_1 + 52.5 b_2 \\
454.0 = 52.5 b_1 + 510.0 b_2 \\
\end{cases} \\
\\
\begin{align*}
\Longleftrightarrow
\begin{bmatrix}83.5 \\ 454.0\end{bmatrix}
&= \begin{bmatrix}16.0 & 52.5 \\ 52.5 & 510.0\end{bmatrix}
\begin{bmatrix}b_1 \\ b_2\end{bmatrix} \\
\\
\Longleftrightarrow
\begin{bmatrix}b_1 \\ b_2\end{bmatrix}
&=\begin{bmatrix}16.0 & 52.5 \\ 52.5 & 510.0\end{bmatrix}^{-1}
\begin{bmatrix}83.5 \\ 454.0\end{bmatrix} \\
\end{align*}
$$
③ 逆行列を用いて連立方程式を解く
３行目の行列の等式を解きます。
# 平方和積和行列を用いた偏回帰係数の算出
np.linalg.inv(sum_prod_dev_mtx.iloc[1:, 1:]) @ sum_prod_dev_mtx.iloc[1:, 0]【実行結果】
偏回帰係数 $${b_1, b_2}$$ が求まりました。
せっかくですので、平方和積和行列を活用する偏回帰係数の算出を関数化しましょう。
引数は、df=データフレーム、cols=[目的変数, 説明変数1, 説明変数2, …] です。
# 関数化 ※pandas, numpy利用
def estimate_coefs_by_sum_prod_dev_mtx(df, cols):
## 偏回帰係数の算出（定数項を除く）
# 偏差の算出
deviation = df[cols].apply(lambda x: x - x.mean(), axis=0).values
# 偏差の平方和積和行列の算出
sum_prod_dev_mtx = deviation.T @ deviation
# 偏回帰係数の算出（定数項を除く）
coefs = np.linalg.inv(sum_prod_dev_mtx[1:, 1:]) @ sum_prod_dev_mtx[1:, 0]
## 定数項の係数の算出
means = df[cols].mean().values
coef0 = means[0] - sum(means[1:] * coefs)
## 戻り値: 係数 [beta_0, beta_1, ..., beta_p]
return [coef0, *coefs]この関数を使って偏回帰係数を算出しましょう。
# 関数の実行
cols = ['配向度', '温度', '圧力'] # 対象列 [y, x1, x2]
coefs = estimate_coefs_by_sum_prod_dev_mtx(data1, cols)
display(Math(f'{coefs}'))【実行結果】
偏回帰係数 $${b_0, b_1, b_2}$$ を推定できました！
重回帰式を求める３「分散共分散行列」編 p.46～48
データの「分散共分散行列」を利用する方法です。
データの最初の列に目的変数を配置することが前提条件です。
計算手順は以下のようになります。
分散共分散行列を算出する
分散共分散行列から連立方程式を構成する
逆行列を用いて連立方程式を解く
２～３番目の手順は「平方和積和行列」編と似ています。
Python で解きましょう。
① 分散共分散行列を算出する
pandas の cov() で分散共分散行列を作成します。
# 分散共分散行列の算出
cov_mtx = data1.cov(ddof=1)
cov_mtx【実行結果】
② 分散共分散行列から連立方程式を構成する p.48
分散共分散行列の２行目以降で「連立方程式」を構成します。
$$
\begin{cases}
9.278 = 1.778 b_1 + 5.833b_2 \\
50.444 = 5.833 b_1 + 56.667 b_2 \\
\end{cases} \\
\\
\begin{align*}
\Longleftrightarrow
\begin{bmatrix}9.278 \\ 50.444\end{bmatrix}
&= \begin{bmatrix}1.788 & 5.833 \\ 5.833 & 56.667\end{bmatrix}
\begin{bmatrix}b_1 \\ b_2\end{bmatrix} \\
\\
\Longleftrightarrow\begin{bmatrix}b_1 \\ b_2\end{bmatrix}
&=\begin{bmatrix}1.788 & 5.833 \\ 5.833 & 56.667\end{bmatrix}^{-1}
\begin{bmatrix}9.278 \\ 50.444\end{bmatrix} \\
\end{align*}
$$
③ 逆行列を用いて連立方程式を解く
３行目の行列の等式を解きます。
# 分散共分散行列を用いた偏回帰係数の算出
np.linalg.inv(cov_mtx.iloc[1:, 1:]) @ cov_mtx.iloc[1:, 0]【実行結果】
偏回帰係数 $${b_1, b_2}$$ が求まりました。
せっかくですので、分散共分散行列を活用する偏回帰係数の算出を関数化しましょう。
引数は、df=データフレーム、cols=[目的変数, 説明変数1, 説明変数2, …] です。
### 分散共分散行列を用いた偏回帰係数の算出 p.46~48
# 関数化
def estimate_coefs_by_cov_mtx(df, cols):
## 偏回帰係数の算出（定数項を除く）
# 分散共分散行列の算出
cov_mtx = df[cols].cov(ddof=1).values
# 偏回帰係数の算出（定数項を除く）
coefs = np.linalg.inv(cov_mtx[1:, 1:]) @ cov_mtx[1:, 0]
## 定数項の係数の算出
means = df[cols].mean().values
coef0 = means[0] - sum(means[1:] * coefs)
## 戻り値: 係数 [beta_0, beta_1, ..., beta_p]
return [coef0, *coefs]この関数を使って偏回帰係数を算出しましょう。
# 関数の実行
cols = ['配向度', '温度', '圧力'] # 対象列 [y, x1, x2]
coefs = estimate_coefs_by_cov_mtx(data1, cols)
display(Math(f'{coefs}'))【実行結果】
偏回帰係数 $${b_0, b_1, b_2}$$ を推定できました！
重回帰式を求める４「正規方程式」編
統計書籍でよく見かける正規方程式から導かれた偏回帰係数の公式を用います。
目的変数ベクトル $${\bm y}$$、説明変数行列 $${\bm X}$$、偏回帰係数ベクトル $${\bm b}$$ として、正規方程式は次のように示されます。
$$
\bm X^{\top} \bm X \bm b = \bm X^{\top} \bm y
$$
行列の$${^{-1}}$$ は逆行列、逆行列の$${^{\top}}$$ は行列の転置です。
$${\bm b}$$ について解くと、次のようになります。
$$
\bm b = (\bm X^{\top} \bm X)^{-1} \bm X^{\top} \bm y
$$
この式を利用して偏回帰係数を求めます。
偏回帰係数 $${b_0}$$ に対応する定数項を説明変数 $${\bm X}$$ に追加しておきます。
### 正規方程式（行列）利用
# 説明変数の設定: [定数項, 変数1, ...]
X = np.hstack([np.ones((len(data1), 1)), data1[['温度', '圧力']].values])
# 目的変数の設定: 縦ベクトル（列がある）
y = data1[['配向度']]
# 偏回帰係数の算出
result1_n = np.linalg.inv(X.T @ X) @ X.T @ y
result1_n.columns = ['偏回帰係数']
# 結果の表示
result1_n【実行結果】
重回帰式を求める５「Pythonライブラリ」編 p.49
テキストは p.49 末尾で統計解析用ソフト SPSS による重回帰式の結果を掲載しています。
Pythonにも重回帰分析ができるライブラリがあります。
SPSS の偏回帰係数、ｔ値、有意確率に相当する統計量を Pythonライブラリでサクッと求めてみましょう。
🖲️statsmodels の ols
ols は最小二乗法 Ordinary Least Squares の頭文字です。
statsmodels には２種類の「表記方法」（API）があります。
この記事では R の formula 記法を使える formula.api （smf）を利用します。
formula 引数に '目的変数 ~ 説明変数1 + 説明変数2 + …' と書きます。
# statsmodels利用
result1_sm = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1).fit()
result1_sm.summary().tables[1]【実行結果】
重回帰分析の結果から偏回帰係数の部分を表示しました。
Intercept（$${b_0}$$、切片）と２つの説明変数が行になっています。
列は、偏回帰係数: coef、標準誤差: std err、検定統計量: t、検定統計量のp値: P>|t|、95%信頼区間: [0.025, 0.975] です。
🖲️pingouin の linear_regression
linear_regression は線形回帰 Linear Regression のことです。
データが pandas データフレームになっていることが前提です。
引数 X に説明変数、Y に目的変数を与えます。
# pingouin利用
result1_pg = pg.linear_regression(X=data1[['温度', '圧力']], y=data1['配向度'])
result1_pg【実行結果】
列は、偏回帰係数の名称: names、偏回帰係数: coef、標準誤差: se、検定統計量: T、検定統計量のp値: pval、自由度調整済み決定係数: adj_r2、95%信頼区間: CI[0.025], CI[0.975] です。
専用ライブラリを使うことで、とてもシンプルに重回帰分析を実践できます！
記事の最後はChatGPTが締めくくります。
今回は料理にたとえて。
📘 ChatGPTのひとこと：
レシピづくり、お疲れさまでした🍳
sympy で最小二乗法の味付けを試し、statsmodels や pingouin で偏回帰係数をしっかり調理して、自分だけのレシピ～重回帰式を完成させましたね。
次回は、出来上がった一皿を味見する時間です。
決定係数、重相関係数、AIC といった「味見スプーン」を使って、モデルの当てはまり具合をチェック。
そして偏回帰係数と標準偏回帰係数を使って、「どの素材がいちばん風味を決めているのか」をじっくり味わいながら解釈していきます😊✨
日常のちょっとした瞬間にもデータの香りを感じながら、次回も一緒に統計のキッチンで学びを深めましょう🔍🎉
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！