第8回 回帰分析
はじめに
本シリーズ「Pythonで学ぶデータ分析実践」の第8回では、回帰分析について解説します。
日常生活で遭遇する多くの現象は、複数の要因によって影響されています。例えば「気温が上がるとビールの売上が増える」「広告費を増やすと売上が伸びる」といった関係です。回帰分析は、こうした変数間の関係をモデル化し、未知のデータに対する予測を行うための手法です。
本記事では、単回帰分析から重回帰分析、正規線形モデル、一般化線形モデルまでを解説し、R²やRMSEによるモデル評価の方法を実践的に紹介します。
シリーズ全体の構成（全12回）
| 回 | テーマ |
|---|---|
| 第1回 | データ分析とは |
| 第2回 | データの前処理 |
| 第3回 | データの把握 |
| 第4回 | 相関分析 |
| 第5回 | 統計的推定 |
| 第6回 | 統計的検定 |
| 第7回 | 分散分析 |
| 第8回 | 回帰分析（本記事） |
| 第9回 | 時系列データ分析 |
| 第10回 | クラス分類 |
| 第11回 | クラスタリング |
| 第12回 | 次元削減 |
8-1. 回帰分析とは
回帰分析とは、ある変数（目的変数）が他の変数（説明変数）によってどのように説明されるかをモデル化し、予測や関係性の解明を行う統計手法です。データの関係を直線や曲線で表現し、未知のデータに対する予測を可能にします。
基本用語
| 用語 | 別名 | 説明 | 例 |
|---|---|---|---|
| 目的変数 | 従属変数、被説明変数、y | 予測したい結果 | ビール売上、住宅価格 |
| 説明変数 | 独立変数、特徴量、x | 結果を説明する要因 | 気温、広告費、面積 |
回帰分析の種類
| 手法 | 説明変数の数 | 目的変数の分布 |
|---|---|---|
| 単回帰分析 | 1つ | 正規分布（連続値） |
| 重回帰分析 | 2つ以上 | 正規分布（連続値） |
| ロジスティック回帰 | 1つ以上 | 二項分布（0/1） |
| ポアソン回帰 | 1つ以上 | ポアソン分布（カウント） |
パラメータの推定方法
回帰分析では、モデルの「傾き」や「切片」を決定する作業をパラメータの推定と呼びます。
| 推定方法 | 説明 | 特徴 |
|---|---|---|
| 最小二乗法（OLS） | データと直線のズレ（残差）の2乗和を最小にする | 計算がシンプル、最も広く使われる |
| 最尤推定法（MLE） | データが観測される確率が最大になるようにパラメータを決める | 確率的な考え方、一般化線形モデルで使用 |
線形モデルの前提条件
線形モデルを正しく適用するには、以下の条件を確認する必要があります。
| 条件 | 説明 | 違反時の対処 |
|---|---|---|
| 線形性 | 説明変数と目的変数が線形関係にある | 対数変換、多項式回帰 |
| 残差の正規性 | 残差が正規分布に従う | 一般化線形モデルの使用 |
| 等分散性 | 残差の分散が一定 | 加重最小二乗法、変数変換 |
| 独立性 | 残差が互いに独立 | 時系列分析の手法を使用 |
| 多重共線性がない | 説明変数間に強い相関がない | 変数の除外、正則化 |
8-2. 単回帰分析
単回帰分析は、1つの説明変数と1つの目的変数の関係を分析する最もシンプルな回帰分析です。データ点に最もフィットする直線（回帰直線）を求め、予測を行います。
単回帰式
$$y = ax + b$$
- y: 目的変数（予測値）
- x: 説明変数
- a: 傾き（回帰係数）— xが1単位増えたときのyの変化量
- b: 切片 — x=0のときのyの値
最小二乗法による係数の算出
実際の値と回帰直線での予測値の差を残差と呼び、残差の2乗和を最小にするのが最小二乗法です。
傾き a の計算：
$$a = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n}(x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}}$$
切片 b の計算：
$$b = \bar{y} - a\bar{x}$$
Pythonでの実装
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import japanize_matplotlib
# サンプルデータ（気温とビール売上）
np.random.seed(42)
temperature = np.array([18, 20, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35])
beer_sales = 50 + 8 * temperature + np.random.normal(0, 15, len(temperature))
# 手動で最小二乗法を計算
x_mean = temperature.mean()
y_mean = beer_sales.mean()
# 傾き a
numerator = np.sum((temperature - x_mean) * (beer_sales - y_mean))
denominator = np.sum((temperature - x_mean) ** 2)
a = numerator / denominator
# 切片 b
b = y_mean - a * x_mean
print(f"=== 最小二乗法による単回帰分析 ===")
print(f"傾き a: {a:.4f}")
print(f"切片 b: {b:.4f}")
print(f"回帰式: y = {a:.2f}x + {b:.2f}")
print(f"\n解釈: 気温が1℃上がると、ビール売上は約{a:.1f}本増加する")
statsmodelsによる実装
import statsmodels.formula.api as smf
# DataFrameの作成
df = pd.DataFrame({"temperature": temperature, "beer_sales": beer_sales})
# OLS（最小二乗法）による線形回帰
model = smf.ols("beer_sales ~ temperature", data=df).fit()
# 結果の表示
print(model.summary())
結果の読み方（summary()の主要な項目）
| 項目 | 意味 |
|---|---|
| R-squared | 決定係数（モデルの説明力） |
| Adj. R-squared | 自由度調整済み決定係数 |
| coef（Intercept） | 切片 b |
| coef（temperature） | 傾き a |
| std err | 係数の標準誤差 |
| t | t値（係数が0かどうかの検定） |
| P>|t| | p値（係数が有意かどうか） |
| [0.025, 0.975] | 係数の95%信頼区間 |
回帰直線の可視化
# 予測値の計算
x_pred = np.linspace(temperature.min()-2, temperature.max()+2, 100)
y_pred = model.predict(pd.DataFrame({"temperature": x_pred}))
# 可視化
plt.figure(figsize=(10, 6))
plt.scatter(temperature, beer_sales, s=80, color='steelblue',
edgecolors='black', zorder=5, label='実測値')
plt.plot(x_pred, y_pred, color='red', linewidth=2,
label=f'回帰直線: y = {a:.2f}x + {b:.2f}')
# 残差を可視化
y_hat = model.predict(df)
for i in range(len(temperature)):
plt.plot([temperature[i], temperature[i]], [beer_sales[i], y_hat.iloc[i]],
color='gray', linestyle='--', alpha=0.5)
plt.xlabel("気温（℃）", fontsize=12)
plt.ylabel("ビール売上（本）", fontsize=12)
plt.title(f"単回帰分析: 気温とビール売上\nR² = {model.rsquared:.4f}", fontsize=14)
plt.legend(fontsize=11)
plt.grid(alpha=0.3)
plt.show()
scikit-learnによる実装
scikit-learnのLinearRegressionクラスで線形回帰モデルを構築します。
LinearRegressionの主要パラメータ：
| パラメータ | デフォルト値 | 意味 |
|---|---|---|
| fit_intercept | True | 切片を計算するかどうか |
| copy_X | True | データのコピーを作成するか |
LinearRegressionの属性：
| 属性 | 意味 |
|---|---|
| coef_ | 偏回帰係数 |
| intercept_ | 切片 |
LinearRegressionのメソッド：
| メソッド | 意味 |
|---|---|
| fit(X, y) | モデルの学習を実施 |
| predict(X) | 学習済みモデルで予測を実施 |
| score(X, y) | 決定係数R²を返す |
| get_params() | モデルのパラメータを取得 |
from sklearn.linear_model import LinearRegression
# モデルの構築
X = temperature.reshape(-1, 1) # 2次元配列に変換
y = beer_sales
lr = LinearRegression()
lr.fit(X, y)
print(f"=== scikit-learn による単回帰分析 ===")
print(f"傾き: {lr.coef_[0]:.4f}")
print(f"切片: {lr.intercept_:.4f}")
print(f"R²スコア: {lr.score(X, y):.4f}")
# 予測
new_temp = np.array([[36]]) # 気温36℃のときの売上を予測
predicted = lr.predict(new_temp)
print(f"\n気温36℃のときの予測売上: {predicted[0]:.1f}本")
8-3. 重回帰分析
重回帰分析は、2つ以上の説明変数を用いて目的変数を予測する手法です。現実の多くの現象は複数の要因に影響されるため、実務では重回帰分析の方がよく使われます。
重回帰式
$$y = b_0 + b_1 x_1 + b_2 x_2 + \cdots + b_p x_p$$
- $b_0$: 切片
- $b_1, b_2, \ldots, b_p$: 偏回帰係数（他の変数を固定したときの影響度）
- $x_1, x_2, \ldots, x_p$: 説明変数
例題：住宅価格の予測
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import japanize_matplotlib
# サンプルデータ（住宅データ）
np.random.seed(42)
n = 100
area = np.random.uniform(40, 120, n) # 面積（㎡）
distance = np.random.uniform(1, 30, n) # 駅距離（分）
age = np.random.uniform(0, 40, n) # 築年数
rooms = np.random.randint(1, 5, n) # 部屋数
# 住宅価格（万円）= 面積の影響 + 駅距離の影響 + 築年数の影響 + ノイズ
price = (30 * area - 50 * distance - 20 * age + 200 * rooms
+ 1000 + np.random.normal(0, 300, n))
df = pd.DataFrame({
"price": price, "area": area, "distance": distance,
"age": age, "rooms": rooms
})
print("=== データの概要 ===")
print(df.describe().round(2))
statsmodelsによる重回帰分析
# 重回帰モデルの構築
model = smf.ols("price ~ area + distance + age + rooms", data=df).fit()
print("\n=== 重回帰分析の結果 ===")
print(model.summary())
結果の解釈
print("\n=== 偏回帰係数の解釈 ===")
for var in ["area", "distance", "age", "rooms"]:
coef = model.params[var]
p_val = model.pvalues[var]
sig = "***" if p_val < 0.001 else "**" if p_val < 0.01 else "*" if p_val < 0.05 else ""
print(f" {var:>10}: {coef:>8.2f} (p={p_val:.4f}) {sig}")
print(f"\n R²: {model.rsquared:.4f}")
print(f" 調整済みR²: {model.rsquared_adj:.4f}")
print(f" AIC: {model.aic:.2f}")
print(f"\n=== 解釈の例 ===")
print(f" 面積が1㎡増えると、住宅価格は約{model.params['area']:.0f}万円上昇")
print(f" 駅距離が1分遠くなると、住宅価格は約{abs(model.params['distance']):.0f}万円下落")
多重共線性の確認（VIF）
**多重共線性（マルチコ）**とは、説明変数の中に相関係数が高い組み合わせがあることです。多重共線性が存在すると、高い決定係数が得られても予測がうまくいかないことがあります。
多重共線性を回避するには、説明変数間の相関係数を相関行列で確認し、高い相関がある変数のいずれかを除去します。
相関行列とヒートマップ
import seaborn as sns
# 相関行列の計算
corr_matrix = df[["area", "distance", "age", "rooms"]].corr()
print("=== 相関行列 ===")
print(corr_matrix.round(3))
# ヒートマップで可視化
plt.figure(figsize=(7, 5))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', fmt='.3f',
vmin=-1, vmax=1, square=True)
plt.title("説明変数間の相関行列（ヒートマップ）")
plt.tight_layout()
plt.show()
VIF（分散膨張因子）
VIF（Variance Inflation Factor）は多重共線性を定量的に評価する指標です。
| VIFの値 | 判断 |
|---|---|
| 1 | 多重共線性なし |
| 1〜5 | 軽度（許容範囲） |
| 5〜10 | 要注意 |
| 10以上 | 多重共線性の疑いあり（変数の除去を検討） |
from statsmodels.stats.outliers_influence import variance_inflation_factor
# VIFの計算
X = df[["area", "distance", "age", "rooms"]]
vif_data = pd.DataFrame()
vif_data["変数"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print("\n=== VIF（分散膨張因子） ===")
print(vif_data)
print("\n※ VIF > 10 の場合、多重共線性の疑いあり")
変数選択
すべての変数をモデルに投入するのではなく、重要な変数を選択することが重要です。
# AICによるステップワイズ的な変数選択
from itertools import combinations
variables = ["area", "distance", "age", "rooms"]
best_aic = np.inf
best_model_vars = None
print("\n=== AICによるモデル選択 ===")
for k in range(1, len(variables) + 1):
for combo in combinations(variables, k):
formula = f"price ~ {' + '.join(combo)}"
m = smf.ols(formula, data=df).fit()
if m.aic < best_aic:
best_aic = m.aic
best_model_vars = combo
print(f"最良モデルの変数: {best_model_vars}")
print(f"最良モデルのAIC: {best_aic:.2f}")
8-4. 正規線形モデル
**正規線形モデル（Normal Linear Model）**は、目的変数が正規分布に従うことを前提とした線形モデルです。単回帰・重回帰分析は正規線形モデルの一種です。
正規線形モデルの構造
$$y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + \varepsilon_i$$
- $\varepsilon_i \sim N(0, \sigma^2)$（残差が正規分布に従う）
正規線形モデルの拡張
正規線形モデルは以下のように拡張できます。
| モデル | 説明変数 | 例 |
|---|---|---|
| 単回帰 | 連続変数1つ | y ~ x |
| 重回帰 | 連続変数複数 | y ~ x1 + x2 + x3 |
| 分散分析 | カテゴリ変数 | y ~ C(group) |
| 共分散分析 | 連続 + カテゴリ混合 | y ~ x + C(group) |
| 多項式回帰 | 高次の項 | y ~ x + x² |
| 交互作用モデル | 変数の交互作用 | y ~ x1 * x2 |
多項式回帰
データが直線的でない場合、多項式回帰を使って非線形な関係をモデル化できます。
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import japanize_matplotlib
# 非線形データの生成
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 5 + 3*x - 0.5*x**2 + np.random.normal(0, 3, 50)
df_poly = pd.DataFrame({"x": x, "y": y})
# 線形モデル（1次）
model_1 = smf.ols("y ~ x", data=df_poly).fit()
# 2次多項式モデル
model_2 = smf.ols("y ~ x + I(x**2)", data=df_poly).fit()
# 3次多項式モデル
model_3 = smf.ols("y ~ x + I(x**2) + I(x**3)", data=df_poly).fit()
# AICの比較
print("=== 多項式回帰のモデル比較 ===")
print(f"1次（線形）: R²={model_1.rsquared:.4f}, AIC={model_1.aic:.2f}")
print(f"2次（二次式）: R²={model_2.rsquared:.4f}, AIC={model_2.aic:.2f}")
print(f"3次（三次式）: R²={model_3.rsquared:.4f}, AIC={model_3.aic:.2f}")
# 可視化
plt.figure(figsize=(10, 6))
plt.scatter(x, y, alpha=0.7, color='steelblue', edgecolors='black', label='実測値')
x_plot = np.linspace(0, 10, 100)
df_plot = pd.DataFrame({"x": x_plot})
plt.plot(x_plot, model_1.predict(df_plot), 'g--', linewidth=2, label=f'1次 (R²={model_1.rsquared:.3f})')
plt.plot(x_plot, model_2.predict(df_plot), 'r-', linewidth=2, label=f'2次 (R²={model_2.rsquared:.3f})')
plt.plot(x_plot, model_3.predict(df_plot), 'purple', linewidth=2, linestyle=':', label=f'3次 (R²={model_3.rsquared:.3f})')
plt.xlabel("x", fontsize=12)
plt.ylabel("y", fontsize=12)
plt.title("多項式回帰の比較", fontsize=14)
plt.legend(fontsize=11)
plt.grid(alpha=0.3)
plt.show()
残差診断
モデルの前提条件が満たされているか、残差プロットで確認します。
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# model_2（2次多項式）の残差診断
residuals = model_2.resid
fitted = model_2.fittedvalues
# 1. 残差 vs 予測値（等分散性の確認）
axes[0,0].scatter(fitted, residuals, alpha=0.7, edgecolors='black')
axes[0,0].axhline(0, color='red', linestyle='--')
axes[0,0].set_xlabel("予測値")
axes[0,0].set_ylabel("残差")
axes[0,0].set_title("残差 vs 予測値（等分散性の確認）")
# 2. 残差のヒストグラム（正規性の確認）
axes[0,1].hist(residuals, bins=15, edgecolor='black', alpha=0.7, density=True)
x_norm = np.linspace(residuals.min(), residuals.max(), 100)
axes[0,1].plot(x_norm, stats.norm.pdf(x_norm, residuals.mean(), residuals.std()),
'r-', linewidth=2)
axes[0,1].set_xlabel("残差")
axes[0,1].set_ylabel("密度")
axes[0,1].set_title("残差のヒストグラム（正規性の確認）")
# 3. Q-Qプロット（正規性の確認）
stats.probplot(residuals, plot=axes[1,0])
axes[1,0].set_title("Q-Qプロット（正規性の確認）")
# 4. 残差の自己相関（独立性の確認）
axes[1,1].plot(residuals, 'o-', alpha=0.7, markersize=4)
axes[1,1].axhline(0, color='red', linestyle='--')
axes[1,1].set_xlabel("観測番号")
axes[1,1].set_ylabel("残差")
axes[1,1].set_title("残差の時系列プロット（独立性の確認）")
plt.tight_layout()
plt.show()
8-5. 一般化線形モデル
**一般化線形モデル（GLM: Generalized Linear Model）**は、正規線形モデルを拡張し、目的変数が正規分布以外の分布（二項分布、ポアソン分布など）に従う場合にも対応できるモデルです。
正規線形モデルとの違い
| 項目 | 正規線形モデル | 一般化線形モデル |
|---|---|---|
| 目的変数の分布 | 正規分布のみ | 正規、二項、ポアソン、ガンマなど |
| リンク関数 | 恒等関数（y = Xβ） | 恒等、ロジット、対数など |
| 推定方法 | 最小二乗法（OLS） | 最尤推定法（MLE） |
| 適用例 | 連続値の予測 | 確率予測、カウント予測 |
GLMの3つの構成要素
| 構成要素 | 説明 |
|---|---|
| 確率分布 | 目的変数が従う分布（正規、二項、ポアソンなど） |
| 線形予測子 | $\eta = \beta_0 + \beta_1 x_1 + \cdots$ |
| リンク関数 | 確率分布の期待値と線形予測子を結びつける関数 |
代表的なGLM
| モデル | 確率分布 | リンク関数 | 目的変数の例 |
|---|---|---|---|
| 線形回帰 | 正規分布 | 恒等（identity） | 売上、温度 |
| ロジスティック回帰 | 二項分布 | ロジット（logit） | 購入/非購入、合格/不合格 |
| ポアソン回帰 | ポアソン分布 | 対数（log） | 来店数、事故件数 |
ロジスティック回帰の実装
目的変数が0/1（二値）の場合に使用します。
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report
import matplotlib.pyplot as plt
import japanize_matplotlib
# サンプルデータ（商品購入予測）
np.random.seed(42)
n = 200
age = np.random.uniform(20, 60, n)
income = np.random.uniform(200, 1000, n)
# 購入確率はageとincomeに依存
logit = -5 + 0.05 * age + 0.005 * income
prob = 1 / (1 + np.exp(-logit))
purchase = np.random.binomial(1, prob, n)
df_log = pd.DataFrame({"age": age, "income": income, "purchase": purchase})
# ロジスティック回帰（statsmodels）
model_logit = smf.logit("purchase ~ age + income", data=df_log).fit()
print("=== ロジスティック回帰の結果 ===")
print(model_logit.summary())
# オッズ比の表示
print(f"\n=== オッズ比 ===")
odds_ratios = np.exp(model_logit.params)
for var in ["age", "income"]:
print(f" {var}: {odds_ratios[var]:.4f}")
print(f" → {var}が1単位増えると、購入のオッズが{(odds_ratios[var]-1)*100:.2f}%変化")
ポアソン回帰の実装
目的変数がカウントデータ（0以上の整数）の場合に使用します。
import statsmodels.api as sm
# サンプルデータ（1日の来店客数）
np.random.seed(42)
n = 100
ad_spend = np.random.uniform(10, 100, n) # 広告費（万円）
is_holiday = np.random.binomial(1, 0.3, n) # 休日フラグ
# 来店客数（ポアソン分布に従う）
log_lambda = 3.5 + 0.01 * ad_spend + 0.3 * is_holiday
visitors = np.random.poisson(np.exp(log_lambda))
df_poisson = pd.DataFrame({
"visitors": visitors, "ad_spend": ad_spend, "is_holiday": is_holiday
})
# ポアソン回帰
model_poisson = smf.glm("visitors ~ ad_spend + is_holiday", data=df_poisson,
family=sm.families.Poisson()).fit()
print("\n=== ポアソン回帰の結果 ===")
print(model_poisson.summary())
# 予測例
print(f"\n=== 予測例 ===")
new_data = pd.DataFrame({"ad_spend": [50, 80], "is_holiday": [0, 1]})
predicted = model_poisson.predict(new_data)
for i, row in new_data.iterrows():
print(f" 広告費{row['ad_spend']:.0f}万円, 休日={'Yes' if row['is_holiday'] else 'No'}"
f" → 予測来店数: {predicted[i]:.0f}人")
8-6. モデル評価（R²・RMSE）
構築したモデルが良いモデルかどうかを評価する指標を紹介します。
主要な評価指標
| 指標 | 正式名称 | 範囲 | 良いモデル |
|---|---|---|---|
| R² | 決定係数 | 0〜1 | 1に近いほど良い |
| 調整済みR² | 自由度調整済み決定係数 | 0〜1 | 変数増加による過大評価を補正 |
| RMSE | 平均二乗誤差の平方根 | 0〜∞ | 0に近いほど良い |
| MAE | 平均絶対誤差 | 0〜∞ | 0に近いほど良い |
| AIC | 赤池情報量規準 | -∞〜∞ | 小さいほど良い |
決定係数 R²
モデルがデータの変動をどれだけ説明できるかを表す指標です。
$$R^2 = 1 - \frac{SS_{res}}{SS_{tot}} = 1 - \frac{\sum(y_i - \hat{y}_i)^2}{\sum(y_i - \bar{y})^2}$$
- $SS_{res}$: 残差平方和（予測と実際のズレ）
- $SS_{tot}$: 全変動（平均からの総ズレ）
- R²=1: 完璧な予測
- R²=0: 平均値で予測するのと同程度
調整済み R²
説明変数の数が増えるとR²は必ず上がるため、変数の数を考慮して補正した指標です。
$$R^2_{adj} = 1 - \frac{(1-R^2)(n-1)}{n-p-1}$$
- n: データ数、p: 説明変数の数
RMSE（Root Mean Squared Error）
予測誤差の大きさを目的変数と同じ単位で表す指標です。
$$RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}$$
AIC（赤池情報量規準）
モデルの複雑さ（変数の数）と適合度のバランスを評価する指標です。
$$AIC = -2\ln L + 2k$$
- L: 最大尤度
- k: パラメータ数（自由度）
AICが低いほど「良い」モデルです。変数を増やすと適合度は上がりますが、同時にペナルティ（2k）が加わるため、過学習を防ぐ効果があります。
Pythonでの評価指標の計算
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
import matplotlib.pyplot as plt
import japanize_matplotlib
# 住宅データの再利用
np.random.seed(42)
n = 200
area = np.random.uniform(40, 120, n)
distance = np.random.uniform(1, 30, n)
age = np.random.uniform(0, 40, n)
rooms = np.random.randint(1, 5, n)
price = 30*area - 50*distance - 20*age + 200*rooms + 1000 + np.random.normal(0, 300, n)
X = np.column_stack([area, distance, age, rooms])
y = price
# 学習データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# モデルの構築
model = LinearRegression()
model.fit(X_train, y_train)
# 予測
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)
# 評価指標の計算
print("=== モデル評価指標 ===")
print(f"\n{'指標':<15} {'学習データ':>12} {'テストデータ':>12}")
print("-" * 45)
print(f"{'R²':<15} {r2_score(y_train, y_pred_train):>12.4f} {r2_score(y_test, y_pred_test):>12.4f}")
print(f"{'RMSE':<15} {np.sqrt(mean_squared_error(y_train, y_pred_train)):>12.2f} "
f"{np.sqrt(mean_squared_error(y_test, y_pred_test)):>12.2f}")
print(f"{'MAE':<15} {mean_absolute_error(y_train, y_pred_train):>12.2f} "
f"{mean_absolute_error(y_test, y_pred_test):>12.2f}")
# 過学習の判断
r2_train = r2_score(y_train, y_pred_train)
r2_test = r2_score(y_test, y_pred_test)
if r2_train - r2_test > 0.1:
print("\n⚠️ 学習データとテストデータのR²に大きな差があります（過学習の可能性）")
else:
print("\n✓ 過学習の兆候は見られません")
予測値 vs 実測値プロット
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 学習データ
axes[0].scatter(y_train, y_pred_train, alpha=0.6, edgecolors='black', s=40)
axes[0].plot([y_train.min(), y_train.max()], [y_train.min(), y_train.max()],
'r--', linewidth=2, label='完璧な予測')
axes[0].set_xlabel("実測値（万円）", fontsize=12)
axes[0].set_ylabel("予測値（万円）", fontsize=12)
axes[0].set_title(f"学習データ（R²={r2_train:.4f}）", fontsize=12)
axes[0].legend()
axes[0].grid(alpha=0.3)
# テストデータ
axes[1].scatter(y_test, y_pred_test, alpha=0.6, edgecolors='black', s=40, color='coral')
axes[1].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()],
'r--', linewidth=2, label='完璧な予測')
axes[1].set_xlabel("実測値（万円）", fontsize=12)
axes[1].set_ylabel("予測値（万円）", fontsize=12)
axes[1].set_title(f"テストデータ（R²={r2_test:.4f}）", fontsize=12)
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
交差検証によるモデル評価
ホールドアウト法はシンプルですが、データが少ない場合は評価が不安定になります。**k交差検証法（k-Fold Cross Validation）**を使うことで、データを最大限に活用した信頼性の高い評価が可能です。
k交差検証法の仕組み
データをk個のグループに分割し、そのうち1つを評価用、残りを学習用として使用します。これをk回繰り返し、各グループが1回ずつ評価用として使われるようにします。
KFoldクラスのパラメータ
| パラメータ | デフォルト値 | 意味 |
|---|---|---|
| n_splits | 5 | データセットの分割数 |
| shuffle | False | 分割前にデータをシャッフルするか |
| random_state | None | 乱数シード |
cross_val_score関数のパラメータ
| パラメータ | 意味 |
|---|---|
| estimator | 使用するモデルオブジェクト |
| X | 入力データ（説明変数） |
| y | 正解データ（目的変数） |
| cv | 分割数（整数）またはKFoldオブジェクト |
| scoring | 評価指標（'r2', 'neg_mean_squared_error'など） |
from sklearn.model_selection import cross_val_score, KFold
# KFoldで分割方法を定義
kf = KFold(n_splits=5, shuffle=True, random_state=42)
# 5分割交差検証
cv_scores_r2 = cross_val_score(LinearRegression(), X, y, cv=kf, scoring='r2')
cv_scores_rmse = -cross_val_score(LinearRegression(), X, y, cv=kf,
scoring='neg_root_mean_squared_error')
print("\n=== 5分割交差検証 ===")
print(f"R²: 平均={cv_scores_r2.mean():.4f} ± {cv_scores_r2.std():.4f}")
print(f"RMSE: 平均={cv_scores_rmse.mean():.2f} ± {cv_scores_rmse.std():.2f}")
sklearn.metricsの誤差関数
| 関数 | 指標 | 意味 |
|---|---|---|
mean_absolute_error(y_test, y_pred) |
MAE | 平均絶対誤差 |
mean_squared_error(y_test, y_pred) |
MSE | 平均二乗誤差 |
np.sqrt(mean_squared_error(...)) |
RMSE | 二乗平均平方根誤差 |
r2_score(y_test, y_pred) |
R² | 決定係数 |
補足：回帰分析の実践的なポイント
正則化（過学習対策）
変数が多い場合やデータが少ない場合は、正則化を用いて過学習を防ぎます。
| 手法 | 特徴 | sklearn |
|---|---|---|
| Ridge回帰（L2） | 係数を小さくする（0にはならない） | Ridge() |
| Lasso回帰（L1） | 不要な変数の係数を0にする（変数選択） | Lasso() |
| ElasticNet | RidgeとLassoの混合 | ElasticNet() |
from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
# スケーリング（正則化の前に必須）
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# Ridge回帰
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)
print(f"\n=== 正則化モデルの比較 ===")
print(f"Ridge R²（テスト）: {ridge.score(X_test_scaled, y_test):.4f}")
# Lasso回帰
lasso = Lasso(alpha=1.0)
lasso.fit(X_train_scaled, y_train)
print(f"Lasso R²（テスト）: {lasso.score(X_test_scaled, y_test):.4f}")
print(f"Lasso 係数: {lasso.coef_}")
print(f" → 0になった係数（不要と判断された変数）の数: {(lasso.coef_ == 0).sum()}")
回帰分析の選択ガイド
目的変数のタイプは？
├── 連続値（売上、価格、温度）
│ ├── 説明変数1つ → 単回帰分析
│ ├── 説明変数複数 → 重回帰分析
│ └── 非線形な関係 → 多項式回帰、スプライン回帰
│
├── 二値（購入/非購入、合格/不合格）
│ └── ロジスティック回帰
│
├── カウントデータ（来店数、事故件数）
│ └── ポアソン回帰
│
└── 生存時間（故障までの時間）
└── Cox回帰
まとめ
本記事では、回帰分析として以下の内容を解説しました。
| 項目 | ポイント |
|---|---|
| 回帰分析とは | 変数間の関係をモデル化し予測を行う手法 |
| 単回帰分析 | 1つの説明変数で直線関係をモデル化（y = ax + b） |
| 重回帰分析 | 複数の説明変数で予測。VIFで多重共線性を確認 |
| 正規線形モデル | 目的変数が正規分布に従う前提の線形モデル |
| 一般化線形モデル | 正規分布以外にも対応（ロジスティック、ポアソン） |
| モデル評価 | R²、RMSE、AICで適合度・予測力・モデルの複雑さを評価 |
回帰分析のチェックリスト
- 目的変数と説明変数を正しく設定したか
- データを学習用とテスト用に分割したか
- 前提条件（線形性・正規性・等分散性・独立性）を確認したか
- 多重共線性（VIF）をチェックしたか
- 適切なモデル（線形/GLM）を選択したか
- R²だけでなくRMSEやAICも確認したか
- 残差診断を行ったか
- 過学習していないか（学習/テストの乖離）を確認したか
- 必要に応じて正則化を検討したか
シリーズ記事一覧
- 【第1回】Pythonで学ぶデータ分析実践 〜データ分析とは〜
- 【第2回】データの前処理（欠損値・外れ値・ダミー変数・データスケーリングなど）
- 【第3回】データの把握（基本統計量・ヒストグラム・箱ひげ図・散布図など）
- 【第4回】相関分析（共分散・相関係数・相関比・連関係数など）
- 【第5回】統計的推定（母集団と標本・点推定・区間推定・信頼区間）
- 【第6回】統計的検定（仮説検定・p値・t検定・カイ二乗検定など）
- 【第7回】分散分析（一元分散分析・二元分散分析・多重比較）
- 【第8回】回帰分析（単回帰・重回帰・一般化線形モデル・モデル評価）
- 【第9回】時系列データ分析（トレンド・季節性・ARIMA・SARIMAなど）
- 【第10回】クラス分類（決定木・ランダムフォレスト・ロジスティック回帰など）
- 【第11回】クラスタリング（k-means・階層クラスタリング・エルボー法など）
- 【第12回】次元削減（主成分分析・寄与率・可視化）
- 【発展編】分析精度の向上（特徴量エンジニアリング・ハイパーパラメータ調整・アンサンブル学習など）
- 【発展編】ニューラルネットワーク入門（パーセプトロン・誤差逆伝播法・TensorFlow・Kerasなど）
次回の第9回では、時間とともに変化するデータを分析する時系列データ分析について解説します。トレンド、季節性の分解、自己回帰（AR）、移動平均（MA）、ARIMA、SARIMAモデルによる予測手法を学びます。