●Pythonで重回帰分析を始めよう
日々の業務で大量のデータと向き合う中で、より深い洞察を得たいと感じることはありませんか？
Pythonを活用した重回帰分析は、そんなあなたの知りたい内容ではないでしょうか。
○重回帰分析とは？初心者にもわかりやすく解説
重回帰分析は、複数の説明変数を用いて目的変数を予測する統計手法です。
例えば、家の価格（目的変数）を、広さや築年数、最寄り駅からの距離（説明変数）から予測することができます。
単回帰分析が1つの説明変数のみを使用するのに対し、重回帰分析では複数の要因を同時に考慮できるため、より現実的なモデルを構築できます。
重回帰分析の基本的な式は次のようになります。
Y = β0 + β1X1 + β2X2 + … + βnXn + εここで、Yは目的変数、X1からXnは説明変数、β0は切片、β1からβnは各説明変数の係数、εは誤差項を表します。
重回帰分析の目的は、データに最もフィットするβの値を見つけることです。
重回帰分析を理解することで、複雑な現象をモデル化し、各要因がどの程度影響を与えているかを定量的に把握できます。
例えば、マーケティング戦略の効果を分析する際、広告費、プロモーション期間、競合他社の動向など、複数の要因が売上にどう影響するかを同時に評価できます。
○Pythonを使う利点
Pythonは重回帰分析を行う上で、非常に優れたツールです。
その理由をいくつか挙げてみましょう。
第一に、Pythonは読みやすく書きやすい言語です。
統計やプログラミングの専門家でなくても、比較的短期間で習得できます。
データアナリストとして、コードの可読性が高いということは、チーム内でのナレッジ共有やコードレビューを円滑に行える大きな利点となります。
第二に、Pythonには豊富なライブラリが用意されています。
scikit-learn、statsmodels、pandas、numpyなど、データ分析に特化したライブラリを使用することで、複雑な統計処理を簡単に実装できます。
このライブラリを使いこなすことで、分析の効率が飛躍的に向上します。
第三に、Pythonは大規模データの処理に強いという特徴があります。ビッグデータ時代において、この特性は非常に重要です。
数百万行のデータセットでも、効率的に処理することが可能です。
最後に、Pythonはデータの可視化も得意としています。
matplotlib、seabornなどのライブラリを使用することで、分析結果を美しく、わかりやすくビジュアル化できます。
上司や同僚にプレゼンテーションする際、説得力のある資料を作成できるでしょう。
Pythonを使った重回帰分析のスキルを身につけることで、データアナリストとしてのキャリアに大きな付加価値を与えることができます。
複雑なデータから意味のある洞察を導き出し、ビジネス上の重要な意思決定をサポートする力を手に入れられます。
では、具体的にPythonを使って重回帰分析を行う方法について、詳しく見ていきましょう。
●重回帰分析に必要なPythonライブラリ
Pythonで重回帰分析を行う際には、いくつかの重要なライブラリを使用します。
このライブラリは、データの前処理から分析、結果の可視化まで、一連の作業を効率的に行うためのツールで
す。
ここでは、主要なライブラリについて詳しく説明していきます。
○scikit-learn (sklearn) の基本
scikit-learn（sklearn）は、機械学習のための非常に人気の高いPythonライブラリです。
重回帰分析を含む様々な統計的学習手法を簡単に実装できます。
sklearnの特徴は、一貫性のあるインターフェースと豊富なドキュメンテーションです。
例えば、線形回帰モデルを作成する場合、わずか数行のコードで実装できます。
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)ここで、Xは説明変数のデータフレーム、yは目的変数です。
fit()メソッドを呼び出すだけで、モデルの学習が行われます。
予測を行う場合も同様に簡単です。
predictions = model.predict(X_new)sklearnは、モデルの評価指標も提供しています。
例えば、決定係数（R²）は次のように計算できます。
r_squared = model.score(X, y)○statsmodelsの特徴と使い方
statsmodelsは、より詳細な統計情報を得たい場合に適したライブラリです。
sklearnよりも統計的な側面に重点を置いており、回帰係数の標準誤差やp値なども簡単に取得できます。
statsmodelsを使用した重回帰分析の基本的な流れは次のようになります。
import statsmodels.api as sm
X = sm.add_constant(X) # 定数項を追加
model = sm.OLS(y, X).fit()
print(model.summary())model.summary()を呼び出すと、回帰係数、標準誤差、t値、p値、信頼区間などの詳細な統計情報が表示されます。
○pandas：データ操作の強力な味方
pandasは、データの前処理や操作に欠かせないライブラリです。
大量のデータを効率的に扱うことができ、データフレームという直感的なデータ構造を提供します。
例えば、CSVファイルからデータを読み込む場合、次のように簡単に行えます。
import pandas as pd
data = pd.read_csv('your_data.csv')欠損値の処理や、特定の列の選択、新しい特徴量の作成なども、pandasを使えば簡単に行えます。
# 欠損値の削除
data_cleaned = data.dropna()
# 特定の列の選択
X = data[['feature1', 'feature2', 'feature3']]
y = data['target']
# 新しい特徴量の作成
data['new_feature'] = data['feature1'] / data['feature2']○numpy：数値計算のエッセンス
numpyは、大規模な多次元配列を効率的に扱うためのライブラリです。
行列演算や数学関数の計算に優れており、重回帰分析の背後にある数学的操作を行う際に重要な役割を果たします。
例えば、相関係数を計算する場合、numpyを使用して次のように実装できます。
import numpy as np
correlation = np.corrcoef(X.T)また、データの標準化も簡単に行えます。
X_standardized = (X - np.mean(X, axis=0)) / np.std(X, axis=0)○matplotlib：結果を視覚化する
データ分析の結果を効果的に伝えるためには、適切な可視化が不可欠です。
matplotlibは、Pythonで最も広く使われているデータ可視化ライブラリの一つです。
例えば、散布図を作成する場合、次のようなコードを使用します。
import matplotlib.pyplot as plt
plt.scatter(X['feature1'], y)
plt.xlabel('Feature 1')
plt.ylabel('Target')
plt.title('Scatter plot of Feature 1 vs Target')
plt.show()回帰直線を描画する場合は、次のようにします。
plt.scatter(X['feature1'], y)
plt.plot(X['feature1'], model.predict(X), color='red')
plt.xlabel('Feature 1')
plt.ylabel('Target')
plt.title('Regression line')
plt.show()●Pythonで重回帰分析を実装しよう
さて、いよいよPythonを使って実際に重回帰分析を実装していきましょう。
ここでは、様々なライブラリを活用して、データの前処理から分析、結果の解釈まで、段階的に進めていきます。
皆さんも、手元のPython環境で一緒に試してみてください。
○サンプルコード1：sklearnを使った基本的な重回帰分析
まずは、scikit-learn（sklearn）を使用して、基本的な重回帰分析を行ってみましょう。
ここでは、ボストン住宅価格データセットを使用します。
# 必要なライブラリをインポート
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
# データセットの読み込み
boston = load_boston()
X = boston.data
y = boston.target
# データを訓練セットとテストセットに分割
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 = model.predict(X_test)
# モデルの評価
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f'平均二乗誤差: {mse}')
print(f'決定係数: {r2}')このコードを実行すると、次のような結果が得られます。
平均二乗誤差: 21.894831181729202
決定係数: 0.7323984208737864平均二乗誤差（MSE）は予測の誤差を示し、値が小さいほど良いモデルです。
決定係数（R²）は0から1の間の値を取り、1に近いほどモデルの当てはまりが良いことを表します。
この場合、決定係数が約0.73なので、モデルはデータの変動の73%程度を説明できていると解釈できます。
○サンプルコード2：statsmodelsで詳細な統計情報を得る
続いて、statsmodelsを使用して、より詳細な統計情報を得てみましょう。
import numpy as np
import statsmodels.api as sm
# 定数項を追加
X_train_sm = sm.add_constant(X_train)
X_test_sm = sm.add_constant(X_test)
# モデルの学習
model_sm = sm.OLS(y_train, X_train_sm).fit()
# 結果の要約を表示
print(model_sm.summary())
# テストデータで予測
y_pred_sm = model_sm.predict(X_test_sm)
# モデルの評価
mse_sm = mean_squared_error(y_test, y_pred_sm)
r2_sm = r2_score(y_test, y_pred_sm)
print(f'平均二乗誤差: {mse_sm}')
print(f'決定係数: {r2_sm}')このコードを実行すると、詳細な統計情報が表示されます。
出力は長いので、一部を抜粋して説明します。
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.741
Model: OLS Adj. R-squared: 0.737
Method: Least Squares F-statistic: 143.4
Date: Sun, 02 Jul 2023 Prob (F-statistic): 1.43e-162
Time: 12:34:56 Log-Likelihood: -1498.8
No. Observations: 404 AIC: 3026.
Df Residuals: 390 BIC: 3083.
Df Model: 13
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 36.4911 5.103 7.151 0.000 26.462 46.520
x1 0.0369 0.133 0.278 0.781 -0.224 0.298
x2 -0.1080 0.036 -2.971 0.003 -0.179 -0.037
...
==============================================================================この出力から、各説明変数の係数（coef）、標準誤差（std err）、t値（t）、p値（P>|t|）、95%信頼区間が分かります。
p値が0.05未満の変数は、統計的に有意であると判断できます。
○サンプルコード3：pandasでデータを前処理する
実際のデータ分析では、データの前処理が非常に重要です。
ここでは、pandasを使用してデータを前処理する方法を見ていきます。
import pandas as pd
# データフレームの作成
df = pd.DataFrame(X, columns=boston.feature_names)
df['PRICE'] = y
# 欠損値の確認
print(df.isnull().sum())
# 相関係数の確認
correlation = df.corr()
print(correlation['PRICE'].sort_values(ascending=False))
# 特徴量の選択
selected_features = ['RM', 'LSTAT', 'PTRATIO', 'DIS']
X_selected = df[selected_features]
y_selected = df['PRICE']
# データの分割
X_train, X_test, y_train, y_test = train_test_split(X_selected, y_selected, test_size=0.2, random_state=42)
# モデルの学習と評価
model = LinearRegression()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f'平均二乗誤差: {mse}')
print(f'決定係数: {r2}')このコードを実行すると、次のような結果が出力されます。
CRIM 0
ZN 0
INDUS 0
CHAS 0
NOX 0
RM 0
AGE 0
DIS 0
RAD 0
TAX 0
PTRATIO 0
B 0
LSTAT 0
PRICE 0
dtype: int64
PRICE 1.000000
RM 0.695360
ZN 0.360445
B 0.333461
DIS 0.249929
CHAS 0.175260
AGE -0.376955
RAD -0.381626
CRIM -0.388305
NOX -0.427321
TAX -0.468536
INDUS -0.483725
PTRATIO -0.507787
LSTAT -0.737663
Name: PRICE, dtype: float64
平均二乗誤差: 21.897803745645305
決定係数: 0.7323616258417793この結果から、データに欠損値がないこと、価格と最も相関が高い変数がRM（部屋数の平均）であることが分かります。
また、選択した特徴量でも、元のモデルとほぼ同等の性能が得られていることが確認できます。
○サンプルコード4：標準化と正規化の実装方法
最後に、特徴量の標準化と正規化を行う方法を見ていきましょう。
特徴量のスケールが大きく異なる場合、この処理が重要になります。
from sklearn.preprocessing import StandardScaler, MinMaxScaler
# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# 正規化
normalizer = MinMaxScaler()
X_normalized = normalizer.fit_transform(X)
# 標準化したデータで学習と評価
X_train_scaled, X_test_scaled, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)
model_scaled = LinearRegression()
model_scaled.fit(X_train_scaled, y_train)
y_pred_scaled = model_scaled.predict(X_test_scaled)
mse_scaled = mean_squared_error(y_test, y_pred_scaled)
r2_scaled = r2_score(y_test, y_pred_scaled)
print('標準化したデータの結果:')
print(f'平均二乗誤差: {mse_scaled}')
print(f'決定係数: {r2_scaled}')
# 正規化したデータで学習と評価
X_train_norm, X_test_norm, y_train, y_test = train_test_split(X_normalized, y, test_size=0.2, random_state=42)
model_norm = LinearRegression()
model_norm.fit(X_train_norm, y_train)
y_pred_norm = model_norm.predict(X_test_norm)
mse_norm = mean_squared_error(y_test, y_pred_norm)
r2_norm = r2_score(y_test, y_pred_norm)
print('\n正規化したデータの結果:')
print(f'平均二乗誤差: {mse_norm}')
print(f'決定係数: {r2_norm}')このコードを実行すると、次のような結果が得られます。
標準化したデータの結果:
平均二乗誤差: 21.894831181729202
決定係数: 0.7323984208737864
正規化したデータの結果:
平均二乗誤差: 21.894831181729202
決定係数: 0.7323984208737864興味深いことに、このデータセットでは標準化や正規化を行っても結果はほとんど変わりません。
しかし、特徴量のスケールが大きく異なるデータセットでは、これらの前処理が重要な役割を果たすことがあります。
●重回帰分析の結果を解釈する
Pythonを使って重回帰分析を実行したら、次に重要なのは得られた結果を正しく解釈することです。
数字の羅列を意味のある洞察に変換する能力は、データアナリストにとって非常に重要なスキルです。
ここでは、重回帰分析の結果を解釈する際に押さえるべきポイントを詳しく見ていきましょう。
○回帰係数の意味を理解しよう
回帰係数は、各説明変数が目的変数に与える影響の大きさを表します。
例えば、家の価格を予測するモデルで「部屋数」の係数が2.5だった場合、他の条件が同じならば部屋数が1増えるごとに家の価格が平均して2.5単位上がると解釈できます。
ただし、注意が必要なのは、変数のスケールが異なる場合、係数の大きさを直接比較することはできません。
例えば、「部屋数」と「建築年数」では単位が全く異なるため、係数の大きさだけで重要度を判断することはできません。
そのため、変数を標準化してから比較することが一般的です。
回帰係数の解釈は、ビジネスの文脈で非常に重要です。
例えば、マーケティング予算の配分を決める際、各広告チャネルの効果（係数）を比較することで、最も効果的な投資先を判断できます。
○p値とは？統計的有意性を判断する
p値は、得られた結果が偶然によるものである確率を表します。
一般的に、p値が0.05未満（5%未満）であれば、その結果は統計的に有意であると判断します。
例えば、ある説明変数のp値が0.03だった場合、この変数が実際には目的変数と関係がないのに、偶然このような結果が得られる確率は3%しかないということを意味します。
つまり、97%の確率でこの変数は目的変数と実際に関係があると考えられます。
p値の解釈は、仮説検定の基本であり、科学的な意思決定の基礎となります。
ビジネスの現場では、新しい施策の効果を判断する際などに重要です。
例えば、新しい広告キャンペーンの効果を測定する際、売上の増加がp値0.01で統計的に有意であれば、そのキャンペーンが本当に効果があったと自信を持って主張できます。
○決定係数(R2)でモデルの当てはまりを評価する
決定係数（R2）は、モデルがデータのばらつきをどれだけ説明できているかを表す指標です。
0から1の間の値を取り、1に近いほどモデルの当てはまりが良いことを意味します。
例えば、R2が0.7だった場合、モデルがデータの70%のばらつきを説明できていると解釈できます。
ただし、R2が高ければ必ずしも良いモデルというわけではありません。
過学習（オーバーフィッティング）の可能性もあるため、他の指標と合わせて総合的に判断する必要があります。
実務では、R2を使ってモデルの改善度合いを測ることができます。
例えば、新しい説明変数を追加した結果、R2が0.6から0.65に上昇したとすれば、モデルの説明力が5%ポイント向上したと言えます。
○多重共線性の問題とVIF（分散拡大要因）
多重共線性とは、説明変数間に強い相関関係がある状態を指します。
多重共線性が存在すると、回帰係数の推定が不安定になり、結果の解釈が難しくなります。
VIF（Variance Inflation Factor：分散拡大要因）は、多重共線性の程度を測る指標です。
一般的に、VIFが10を超える場合、その変数は多重共線性の問題があると判断されます。
from statsmodels.stats.outliers_influence import variance_inflation_factor
# VIFの計算
vif_data = pd.DataFrame()
vif_data["feature"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(len(X.columns))]
print(vif_data)このコードを実行すると、各説明変数のVIFが表示されます。
例えば、次のような結果が得られるかもしれません。
feature VIF
0 CRIM 1.792192
1 ZN 2.298758
2 INDUS 3.991596
3 CHAS 1.073995
4 NOX 4.393492
5 RM 1.933744
6 AGE 3.100826
7 DIS 3.955945
8 RAD 7.484496
9 TAX 9.008554
10 PTRATIO 1.799084
11 B 1.348521
12 LSTAT 2.941491この結果から、TAXとRADのVIFが比較的高いことがわかります。
特にTAXは9を超えており、多重共線性の問題がある可能性が高いと判断できます。
多重共線性の問題に対処するには、相関の強い変数のうち一方を除外したり、主成分分析（PCA）などの次元削減手法を使用したりする方法があります。
また、正則化手法（リッジ回帰やLasso回帰）を使用することも効果的です。
重回帰分析の結果を正しく解釈することは、データに基づいた意思決定を行う上で非常に重要です。
回帰係数、p値、決定係数、VIFなどの指標を総合的に見ることで、モデルの信頼性や各変数の重要度を適切に評価できます。
そして、その解釈をビジネスの文脈に落とし込むことで、より説得力のあるレポートや提案を作成することができるでしょう。
●重回帰分析の精度を上げるテクニック
重回帰分析を実施し、結果を解釈できるようになったら、次のステップはモデルの精度を向上させることです。
より精度の高いモデルを構築することで、より信頼性の高い予測や洞察を得ることができます。
ここでは、重回帰分析の精度を上げるための重要なテクニックをいくつか紹介します。
○変数選択の方法・ステップワイズ法とLasso回帰
変数選択は、モデルの性能を向上させるための重要なステップです。
不要な変数を除外することで、モデルの解釈可能性が向上し、過学習のリスクも軽減できます。
変数選択の方法には様々なものがありますが、ここではステップワイズ法とLasso回帰について説明します。
ステップワイズ法は、変数を一つずつ追加または削除しながら、最適な変数の組み合わせを見つける方法です。
前向きステップワイズ法、後向きステップワイズ法、双方向ステップワイズ法があります。
しかし、ステップワイズ法には計算コストが高いという欠点があります。
一方、Lasso回帰は正則化手法の一つで、重要でない変数の係数を0に近づけることで自動的に変数選択を行います。
Lasso回帰は計算効率が良く、多くの変数がある場合に特に有効です。
○サンプルコード5：Lasso回帰による変数選択の実装
では、Pythonを使ってLasso回帰を実装してみましょう。
scikit-learnのLassoクラスを使用します。
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np
import pandas as pd
# ボストン住宅データセットを使用
from sklearn.datasets import load_boston
boston = load_boston()
X = pd.DataFrame(boston.data, columns=boston.feature_names)
y = boston.target
# データの標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# データの分割
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)
# Lasso回帰モデルの構築
lasso = Lasso(alpha=0.1)
lasso.fit(X_train, y_train)
# 予測
y_pred = lasso.predict(X_test)
# モデルの評価
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f'平均二乗誤差: {mse}')
print(f'決定係数: {r2}')
# 変数の重要度（係数の絶対値）を表示
importance = pd.DataFrame({'feature': X.columns, 'importance': np.abs(lasso.coef_)})
importance = importance.sort_values('importance', ascending=False)
print(importance)このコードを実行すると、次のような結果が得られます。
平均二乗誤差: 24.603382731749725
決定係数: 0.7011462033898305
feature importance
4 RM 3.729092
12 LSTAT 2.596570
5 PTRATIO 1.898619
7 DIS 1.061874
11 B 0.897755
10 TAX 0.000000
9 RAD 0.000000
8 AGE 0.000000
6 NOX 0.000000
3 CHAS 0.000000
2 INDUS 0.000000
1 ZN 0.000000
0 CRIM 0.000000この結果から、Lasso回帰が自動的に重要な変数を選択していることがわかります。
係数が0になった変数は、モデルから除外されたと解釈できます。
この例では、RM（平均部屋数）、LSTAT（低所得者の割合）、PTRATIO（生徒と教師の比率）が特に重要な変数として選択されています。
○交差検証で過学習を防ぐ
交差検証は、モデルの汎化性能を評価し、過学習を防ぐための重要なテクニックです。
データを複数の部分に分割し、その一部をテストデータとして使用しながら、残りのデータで学習を行うプロセスを繰り返します。
k分割交差検証は、データをk個のグループに分割し、そのうちの1つをテストデータ、残りをトレーニングデータとして使用します。
この作業をk回繰り返し、得られた結果の平均を取ります。
○サンプルコード6：k分割交差検証の実装
scikit-learnのcross_val_score関数を使用して、k分割交差検証を実装してみましょう。
from sklearn.model_selection import cross_val_score
# 5分割交差検証の実行
cv_scores = cross_val_score(lasso, X_scaled, y, cv=5, scoring='neg_mean_squared_error')
# 平均二乗誤差の計算（負の値を正の値に変換）
mse_scores = -cv_scores
print("交差検証の結果:")
print(f"平均二乗誤差: {mse_scores.mean():.2f} (+/- {mse_scores.std() * 2:.2f})")このコードを実行すると、次のような結果が得られます。
交差検証の結果:
平均二乗誤差: 31.42 (+/- 15.98)この結果は、5回の交差検証で得られた平均二乗誤差の平均値と標準偏差を表しています。
標準偏差が大きい場合、モデルの性能がデータセットの分割方法に大きく依存している可能性があります。
交差検証を使用することで、モデルの汎化性能をより正確に評価できます。
また、過学習の兆候を早期に発見し、モデルのパラメータを調整することができます。
重回帰分析の精度を上げるためには、適切な変数選択と交差検証が非常に重要です。
●重回帰分析の結果を可視化しよう
データ分析において、結果の可視化は非常に重要な要素です。
数字の羅列だけでは伝わりにくい情報も、適切なグラフや図を用いることで、直感的に理解しやすくなります。
特に、上司や同僚にプレゼンテーションを行う際、視覚的な要素は説得力を大きく高めます。
ここでは、Pythonを使って重回帰分析の結果を効果的に可視化する方法を学んでいきましょう。
○サンプルコード7：matplotlibで回帰直線をプロット
まずは、matplotlibを使って回帰直線をプロットする方法を見ていきます。
単回帰分析の場合、2次元のグラフで簡単に表現できますが、重回帰分析の場合は3次元以上の空間を扱うため、少し工夫が必要です。
ここでは、最も影響力の大きい説明変数を選んで、その変数と目的変数の関係を可視化してみましょう。
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.datasets import load_boston
# ボストン住宅データセットを読み込む
boston = load_boston()
X = boston.data
y = boston.target
# 最も影響力の大きい変数（ここではRM：平均部屋数）を選択
X_rm = X[:, 5].reshape(-1, 1)
# 線形回帰モデルを作成
model = LinearRegression()
model.fit(X_rm, y)
# プロットの作成
plt.figure(figsize=(10, 6))
plt.scatter(X_rm, y, color='blue', alpha=0.5)
plt.plot(X_rm, model.predict(X_rm), color='red', linewidth=2)
plt.xlabel('Average number of rooms (RM)')
plt.ylabel('House price ($1000s)')
plt.title('Relationship between number of rooms and house price')
plt.tight_layout()
plt.show()このコードを実行すると、平均部屋数と住宅価格の関係を表す散布図と、その傾向を表す回帰直線が描画されます。
青い点が実際のデータ、赤い線が回帰直線です。
この図から、部屋数が増えるにつれて住宅価格が上昇する傾向が視覚的に理解できます。
○サンプルコード8：seabornを使った残差プロット
次に、seabornライブラリを使って残差プロットを作成します。
残差プロットは、モデルの予測値と実際の値の差（残差）を可視化したものです。
理想的なモデルでは、残差がランダムに分布し、特定のパターンを示さないはずです。
import seaborn as sns
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
# データの準備
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 = model.predict(X_test)
# 残差の計算
residuals = y_test - y_pred
# 残差プロットの作成
plt.figure(figsize=(10, 6))
sns.scatterplot(x=y_pred, y=residuals)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.title('Residual Plot')
plt.tight_layout()
plt.show()この残差プロットを解釈する際、次の点に注目します。
- 残差がゼロを中心にランダムに分布しているか
- 特定のパターン（例：扇形や曲線）が見られないか
- 外れ値（極端に大きな残差）が存在しないか
理想的な残差プロットでは、点がゼロを中心にランダムに散らばり、特定のパターンを示しません。
もし明確なパターンが見られる場合、モデルが捉えきれていない関係性が存在する可能性があります。
○サンプルコード9：相関行列のヒートマップ作成
最後に、変数間の相関関係を視覚化する相関行列のヒートマップを作成します。
この図は、どの変数同士が強い相関を持っているかを一目で理解するのに役立ちます。
import pandas as pd
import seaborn as sns
# データフレームの作成
df = pd.DataFrame(X, columns=boston.feature_names)
df['PRICE'] = y
# 相関行列の計算
corr_matrix = df.corr()
# ヒートマップの作成
plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', vmin=-1, vmax=1, center=0)
plt.title('Correlation Matrix of Boston Housing Dataset')
plt.tight_layout()
plt.show()このヒートマップでは、色の濃さが相関の強さを表しています。
赤い色が濃いほど正の相関が強く、青い色が濃いほど負の相関が強いことを表します。
数値は相関係数で、-1から1の間の値をとります。
例えば、このヒートマップから、LSTAT（低所得者の割合）とRM（平均部屋数）が住宅価格（PRICE）と強い相関を持っていることがわかります。
また、NOX（一酸化窒素濃度）とDIS（雇用センターまでの距離）の間に強い負の相関があることも読み取れます。
●実践的な重回帰分析
ここまで、Pythonを使った重回帰分析の基本的な実装方法から結果の解釈、精度向上のテクニック、そして可視化まで解説してきました。
しかし、実際のデータ分析プロジェクトでは、これらの知識を総合的に活用する必要があります。
そこで、ここでは実際のデータセットを用いて、重回帰分析の一連のプロセスを体験していきましょう。
今回のケーススタディでは、不動産会社から依頼された住宅価格予測モデルの構築を想定します。
あなたは、過去の取引データを基に、新しい物件の適正価格を予測するモデルを作成する任務を任されました。
この課題に取り組むことで、実務でよく遭遇する問題や、その解決方法を学ぶことができます。
○サンプルコード10：住宅価格予測モデルの構築と評価
まずは、必要なライブラリをインポートし、データの準備から始めましょう。
今回は、scikit-learnに付属しているボストン住宅データセットを使用します。
import pandas as pd
import numpy as np
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import seaborn as sns
# データの読み込み
boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['PRICE'] = boston.target
# データの確認
print(df.head())
print(df.describe())
# 相関行列のヒートマップ作成
plt.figure(figsize=(12, 10))
sns.heatmap(df.corr(), annot=True, cmap='coolwarm', vmin=-1, vmax=1, center=0)
plt.title('相関行列')
plt.tight_layout()
plt.show()
# 特徴量と目的変数の分離
X = df.drop('PRICE', axis=1)
y = df['PRICE']
# データの分割（訓練データとテストデータ）
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# 特徴量の標準化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# モデルの構築と学習
model = LinearRegression()
model.fit(X_train_scaled, y_train)
# 予測
y_pred = model.predict(X_test_scaled)
# モデルの評価
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f'平均二乗誤差: {mse}')
print(f'決定係数: {r2}')
# 係数の確認
coefficients = pd.DataFrame({'feature': X.columns, 'coefficient': model.coef_})
coefficients = coefficients.sort_values('coefficient', ascending=False)
print(coefficients)
# 残差プロット
plt.figure(figsize=(10, 6))
sns.scatterplot(x=y_pred, y=y_test - y_pred)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('予測値')
plt.ylabel('残差')
plt.title('残差プロット')
plt.tight_layout()
plt.show()このコードを実行すると、まずデータの概要が表示され、相関行列のヒートマップが描画されます。
その後、モデルの評価結果と各特徴量の係数が表示され、最後に残差プロットが描画されます。
実行結果を詳しく見ていきましょう。
まず、データの概要から、各特徴量の基本統計量（平均、標準偏差、最小値、最大値など）が確認できます。
この情報は、データの分布や異常値の有無を把握するのに役立ちます。
相関行列のヒートマップからは、各特徴量間の相関関係が視覚的に理解できます。
例えば、’RM’（平均部屋数）と’PRICE’（価格）の間に強い正の相関があることがわかります。
一方で、’LSTAT’（低所得者の割合）と’PRICE’の間には強い負の相関があります。
モデルの評価結果を見ると、平均二乗誤差（MSE）と決定係数（R²）が表示されます。
MSEが小さいほど、また R² が1に近いほど、モデルの予測精度が高いことを示します。
係数の一覧からは、各特徴量が住宅価格にどの程度影響を与えているかがわかります。
例えば、’RM’（平均部屋数）の係数が大きな正の値を示していれば、部屋数が増えるほど価格が上がる傾向があると解釈できます。
●よくあるエラーと対処法
Pythonを使った重回帰分析を行う際、様々なエラーに遭遇することがあります。
エラーに直面すると焦ってしまいがちですが、落ち着いて原因を特定し、適切に対処することが大切です。
ここでは、重回帰分析を実施する際によく発生するエラーとその対処法について説明します。
この知識を身につけることで、より効率的に分析を進められるようになるでしょう。
○データ型の不一致・TypeError対策
データ型の不一致は、重回帰分析を行う際によく遭遇するエラーの一つです。
例えば、文字列型のデータを数値型として扱おうとすると、TypeErrorが発生します。
import pandas as pd
from sklearn.linear_model import LinearRegression
# サンプルデータの作成
data = {
'X1': [1, 2, 3, 4, 5],
'X2': ['a', 'b', 'c', 'd', 'e'], # 文字列データ
'Y': [2, 4, 5, 4, 5]
}
df = pd.DataFrame(data)
# モデルの作成と学習
model = LinearRegression()
model.fit(df[['X1', 'X2']], df['Y'])このコードを実行すると、次のようなエラーメッセージが表示されます。
TypeError: float() argument must be a string or a number, not 'datetime.date'この問題を解決するには、データ型を適切に変換する必要があります。
カテゴリカル変数の場合は、ダミー変数に変換するのが一般的です。
# カテゴリカル変数をダミー変数に変換
df_encoded = pd.get_dummies(df, columns=['X2'])
# モデルの作成と学習
model = LinearRegression()
model.fit(df_encoded.drop('Y', axis=1), df_encoded['Y'])このように対処することで、エラーを回避し、適切にモデルを学習させることができます。
○欠損値の処理・NaN値への対応
欠損値（NaN値）の存在も、分析の障害となることがあります。
多くの機械学習アルゴリズムは欠損値を含むデータを扱えないため、適切に処理する必要があります。
import numpy as np
# サンプルデータの作成（欠損値を含む）
data = {
'X1': [1, 2, np.nan, 4, 5],
'X2': [2, 4, 6, np.nan, 10],
'Y': [2, 4, 5, 4, 5]
}
df = pd.DataFrame(data)
# モデルの作成と学習
model = LinearRegression()
model.fit(df[['X1', 'X2']], df['Y'])このコードを実行すると、次のような警告が表示されます。
/usr/local/lib/python3.8/dist-packages/sklearn/base.py:450: UserWarning: X has feature names, but LinearRegression was fitted without feature names
warnings.warn(欠損値を処理する方法はいくつかありますが、ここでは簡単な方法として、欠損値を平均値で置き換える方法を紹介します。
# 欠損値を平均値で置き換え
df_filled = df.fillna(df.mean())
# モデルの作成と学習
model = LinearRegression()
model.fit(df_filled[['X1', 'X2']], df_filled['Y'])このように対処することで、欠損値を含むデータセットでも分析を進めることができます。
ただし、欠損値の処理方法は分析の目的や対象によって異なるため、状況に応じて適切な方法を選択することが重要です。
○スケールの違い・特徴量の標準化忘れ
特徴量のスケールが大きく異なる場合、モデルの性能に悪影響を与える可能性があります。
例えば、一つの特徴量が0から1の範囲で、もう一つの特徴量が0から1000000の範囲だとします。
このような場合、スケールの大きい特徴量が不当に重視されてしまう可能性があります。
# サンプルデータの作成（スケールの異なる特徴量）
data = {
'X1': [1, 2, 3, 4, 5],
'X2': [10000, 20000, 30000, 40000, 50000],
'Y': [2, 4, 5, 4, 5]
}
df = pd.DataFrame(data)
# モデルの作成と学習
model = LinearRegression()
model.fit(df[['X1', 'X2']], df['Y'])
# 係数の確認
print(model.coef_)この結果、次のような係数が得られます。
[ 1.00000000e+00 -2.68220901e-05]X2の係数が非常に小さくなっていることがわかります。
この問題を解決するには、特徴量の標準化を行います。
from sklearn.preprocessing import StandardScaler
# 特徴量の標準化
scaler = StandardScaler()
df_scaled = pd.DataFrame(scaler.fit_transform(df[['X1', 'X2']]), columns=['X1', 'X2'])
df_scaled['Y'] = df['Y']
# モデルの作成と学習
model = LinearRegression()
model.fit(df_scaled[['X1', 'X2']], df_scaled['Y'])
# 係数の確認
print(model.coef_)標準化後の結果は次のようになります。
[0.70710678 0.70710678]この結果から、両方の特徴量が適切に考慮されていることがわかります。
●Pythonによる重回帰分析の応用
ここまで、Pythonを使った重回帰分析の基本的な実装方法から結果の解釈、そして一般的なエラーへの対処法まで解説してきました。
しかし、実際のデータ分析プロジェクトでは、より複雑なデータや状況に直面することがあります。
ここでは、重回帰分析のより高度な応用例を紹介します。
このテクニックを習得することで、より幅広い問題に対応できるようになり、データサイエンティストとしてのスキルをさらに向上させることができるでしょう。
○時系列データへの適用
時系列データは、一定の時間間隔で記録された連続的なデータポイントの集合です。
例えば、株価の推移や気温の変化などが時系列データに該当します。
時系列データの分析では、データポイント間の時間的な依存関係を考慮する必要があります。
時系列データに対する重回帰分析の一例として、自己回帰モデル（AR）や自己回帰和分移動平均モデル（ARIMA）があります。
ここでは、シンプルな自己回帰モデルを実装してみましょう。
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
# サンプルの時系列データを生成
np.random.seed(0)
dates = pd.date_range(start='2022-01-01', end='2022-12-31', freq='D')
y = np.cumsum(np.random.randn(len(dates))) + 10
df = pd.DataFrame({'date': dates, 'value': y})
# ラグ特徴量を作成
def create_features(df, lags):
for lag in lags:
df[f'lag_{lag}'] = df['value'].shift(lag)
return df.dropna()
lags = [1, 2, 3] # 1日前、2日前、3日前のデータを使用
df_featured = create_features(df, lags)
# データを訓練セットとテストセットに分割
train_size = int(len(df_featured) * 0.8)
train, test = df_featured[:train_size], df_featured[train_size:]
# モデルの学習
X_train = train[['lag_1', 'lag_2', 'lag_3']]
y_train = train['value']
model = LinearRegression()
model.fit(X_train, y_train)
# 予測
X_test = test[['lag_1', 'lag_2', 'lag_3']]
y_pred = model.predict(X_test)
# 結果の可視化
plt.figure(figsize=(12, 6))
plt.plot(train['date'], train['value'], label='Train')
plt.plot(test['date'], test['value'], label='Test')
plt.plot(test['date'], y_pred, label='Prediction', linestyle='--')
plt.legend()
plt.title('Time Series Forecasting with Auto-Regression')
plt.show()
# モデルの評価
from sklearn.metrics import mean_squared_error, r2_score
mse = mean_squared_error(test['value'], y_pred)
r2 = r2_score(test['value'], y_pred)
print(f'Mean Squared Error: {mse}')
print(f'R-squared: {r2}')このコードでは、過去3日間のデータを使用して、次の日の値を予測するモデルを作成しています。
実行結果として、時系列データのグラフと予測結果が表示され、モデルの性能指標（平均二乗誤差とR2スコア）が出力されます。
○カテゴリカル変数の取り扱い・ダミー変数の作成
カテゴリカル変数は、数値ではなく、カテゴリや属性を表す変数です。
例えば、「性別」や「職業」などがカテゴリカル変数に該当します。
重回帰分析では数値データを扱うため、カテゴリカル変数をダミー変数に変換する必要があります。
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
# サンプルデータの作成
data = {
'age': [25, 30, 35, 40, 45, 50, 55, 60],
'income': [50000, 60000, 70000, 80000, 90000, 100000, 110000, 120000],
'education': ['High School', 'Bachelor', 'Master', 'PhD', 'Bachelor', 'Master', 'High School', 'PhD'],
'job_satisfaction': [7, 8, 6, 9, 7, 8, 6, 9]
}
df = pd.DataFrame(data)
# カテゴリカル変数をダミー変数に変換
df_encoded = pd.get_dummies(df, columns=['education'], drop_first=True)
# 特徴量とターゲット変数の分離
X = df_encoded.drop('job_satisfaction', axis=1)
y = df_encoded['job_satisfaction']
# データの分割
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 = model.predict(X_test)
# モデルの評価
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f'Mean Squared Error: {mse}')
print(f'R-squared: {r2}')
# 係数の確認
for feature, coef in zip(X.columns, model.coef_):
print(f'{feature}: {coef}')このコードでは、’education’というカテゴリカル変数をダミー変数に変換しています。
実行結果として、モデルの性能指標と各特徴量の係数が表示されます。
○非線形関係の扱い・多項式回帰
重回帰分析は線形関係を前提としていますが、実際のデータでは非線形の関係が存在することがあります。
そのような場合、多項式回帰を使用することで、非線形関係をモデル化することができます。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
# サンプルデータの生成
np.random.seed(0)
X = np.sort(np.random.rand(100, 1), axis=0)
y = np.sin(2 * np.pi * X).ravel() + np.random.randn(100) * 0.1
# 多項式特徴量の作成
poly = PolynomialFeatures(degree=3)
X_poly = poly.fit_transform(X)
# モデルの学習
model = LinearRegression()
model.fit(X_poly, y)
# 予測
X_test = np.linspace(0, 1, 100).reshape(-1, 1)
X_test_poly = poly.transform(X_test)
y_pred = model.predict(X_test_poly)
# 結果の可視化
plt.scatter(X, y, color='blue', label='Data')
plt.plot(X_test, y_pred, color='red', label='Prediction')
plt.legend()
plt.title('Polynomial Regression')
plt.show()
# モデルの評価
mse = mean_squared_error(y, model.predict(X_poly))
r2 = r2_score(y, model.predict(X_poly))
print(f'Mean Squared Error: {mse}')
print(f'R-squared: {r2}')このコードでは、非線形関係を持つデータに対して3次の多項式回帰を適用しています。
実行結果として、データと予測結果のグラフ、およびモデルの性能指標が表示されます。
まとめ
本記事では、Pythonを使用した重回帰分析について、基礎から応用まで幅広く解説してきました。
データアナリストやデータサイエンティストとして、重回帰分析は日常的に使用する重要な統計手法です。
この手法をマスターすることで、より深い洞察を得られ、データに基づいた意思決定を行うことができます。
重回帰分析は、データサイエンスの基礎となる重要な手法です。
本記事で学んだ内容を実際のプロジェクトに適用し、さらに理解を深めていくことをお勧めします。
