【初心者用】Pythonでできる統計手法|重回帰分析

スポンサーリンク
Python
【2021年】ゼロから統計学を独学したい人が読むべき書籍10冊
データサイエンスについて興味があり、統計学を学びたいけど、どんな学習方法がいいか分からない 統計学を独学するために必要な書籍が知りたい 統計...
【初心者】Rと統計学をいっぺんに学ぶ最初の5冊
Rとはオープンソースの統計解析用のプログラミング言語です。 誰でも無料で使えます。統計解析にはSPSSやSASなどが有名ですが、いずれも有料です。大学や研究機関でないと使えないですね。個人でも統計解析できた方が、職場に依存せずにス...

Pythonはデータ分析を行う上で、必要なプログラミング言語ですが、主には機械学習で用いられます。しかし、データ分析には統計学も重要です。特に、重回帰分析は非常に使い勝手が良く、ビジネスでもよく使われています。機械学習でも重回帰分析(線形回帰分析)は使われますが、その用途は予測に使われることが多いです。

重回帰分析|初心者や文系素人でも直感的に理解できる

【初心者】コピペでできる線形回帰

統計学ではデータの解釈や説明に使われます。
今回は、統計学からの重回帰分析をPythonを使って実装して理解を深めていきましょう。

必要なライブラリをインポート

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
from sklearn.datasets import load_boston
import statsmodels.api as sm

# 警告メッセージを表示しない
import warnings
warnings.filterwarnings('ignore')

データは scikit-learn のデータセットであるボストンの住宅価格(boston)を使って説明していきます。

線形回帰は scikit-learn にもありますが、統計に特化するなら statsmodels を使ったほうが簡単です。

では、ボストンのデータセットを1つのデータフレームの形に整えます。

boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['PRICE'] = boston.target
df

出力:

14列、506行のデータフレームになっていることが確認できます。

データの前処理

では、このデータフレームを重回帰分析ができるように前処理していきます。
前処理は適切に行わないと、分析結果が間違うことがあるので、慎重に行わなければなりませんが、今回は、割愛し重回帰分析のみに焦点を当てて説明します。

y = df['PRICE']
X = df[["INDUS", "RM", "PTRATIO", "LSTAT"]]

X_con = sm.add_constant(X)

目的変数は、住宅価格である変数「PRICE」としますので、これを「y」とします。

説明変数は、変数「INDUS」「RM」「PTRATIO」「LSTAT」の4つとします。この説明変数の選択には特に大きな意味はありません。多すぎると煩雑になるため、例として選択したと思って下さい。

一番下のコードは、すべての値が、1.0 である列を作成しています。これは、切片を作るために必要な処理です。重回帰分析では切片が必要です。相当特殊な場合は切片を設けないことがあるかもしれませんが、わたしは知りません。

では、各変数を確認してみましょう。

y

出力:

yには「PRICE」が入っています。回帰分析の目的変数は連続変数でしたね。

X

出力:

Xには説明変数の4つと、「const」という変数があります。「const」は切片を設定するためのもので、データはすべて 1.0 になっています。

機械学習では予測を重視しますので、学習データとテストデータに分割して、学習データのみを回帰分析します。

でも、統計学ではデータの説明を重視するため、すべてのデータを使って回帰分析します。ここは、用途に合わせて分割するしないを判断して下さい。

重回帰分析の実装

model = sm.OLS(y, X_con)  # モデルを作成

result = model.fit()      # モデルを適応

以上で、重回帰分析ができます。非常に簡単ですね。

1行目で目的変数と説明変数を指定したモデルを作成します。

2行目でそのモデルを適応させて、「result」としています。

結果の確認

result.summary()          # 結果を確認

出力:

いっぱい出力結果がありますが、基本的で最低限のものだけ解説します。

Ⓐは、自由度調整済決定係数です。「R²」ですね。これは、1に近いほど良いです。解釈としては、今回の重回帰分析のモデルが実際のデータをどの程度説明できるかを表しています。

Ⓒは、編回帰係数です。Ⓑはそれに対する説明変数です。「INDUS」の編回帰係数は、0.0076 ですね。編回帰係数は、「他の説明変数が一定であった場合に、その変数を1変化させると、目的変数がどの程度変化するか」を表しています。

PRICE = (0.0076×INDUS) + (4.5152×RM) + (-0.9351×PTRATIO) + (-0.5757×LSTAT) +18.6150(切片)

このような計算式で表すことができます。

ここの注意が必要です。先ほどの編回帰係数の説明で、「他の説明変数が一定であった場合」といいましたが、これが曲者です。他の説明変数が一定であり、1つの説明変数のみが変化することがありえるでしょうか。これは、各変数の相関が0ということです。それはあり得ないですね。

この説明変数間の関係性の強さが多重共線性です。多重共線性はあるないではなく、強いか弱いです。たとえば、相関係数が0.3以下なら問題ない、などと判断します。または、VIFが10以下。VIFはPythonでの計算方法がわかりませんでした。すみません。Rでは簡単です。

Ⓓは、P値です。ふつうは、分析前に有意水準を決めておき、その有意水準より小さければ有意差があると判断します。一般的には、0.05(5%)を設定しますが、あまり気にしなくてもいいです。今回の結果のような、「INDUS」が 0.862 であれば、「INDUS」は目的変数に関係ないかもね、と判断してもいいかもしれません。

しかし、P値の解釈は意外と難しいので、データ数が多ければ気にしなくてもよいです。

標準化

どの説明変数の影響が大きいかを知りたいときは、標準化します。

標準化した後は、同じコードで大丈夫です。

from sklearn.preprocessing import StandardScaler
####### 標準化 ########################################################
scaler = StandardScaler() # インスタンス
scaler.fit(np.array(df)) # 標準化を適応
df_std = scaler.transform(np.array(df)) # 標準化する
df_std = pd.DataFrame(df_std, columns=df.columns) # データフレームに変換
######################################################################
X = df_std[["INDUS", "RM", "PTRATIO", "LSTAT"]]
y = df_std['PRICE']

X_con = sm.add_constant(X)

model = sm.OLS(y, X_con) # モデルを作成

result = model.fit() # モデルを適応

result.summary() # 結果を確認

出力:

coef(回帰係数)の数値だけ、変化します。

標準化は、平均値が0、標準偏差が1にすることなので、説明変数間で比較することができます。

ここでは、「LSTAT」の数値が一番大きい(絶対値で)ので、目的変数に大きな影響があると言えます。

標準化した回帰係数を標準化回帰係数といいます。この標準化回帰係数が一番活用できる数値だと思います。

つまり、「標準化回帰係数が大きい説明変数を変化させることで、目的変数を改善できる」と言えます。厳密にはこれは因果関係がはっきりしないと言えませんが。。。

残差を確認

回帰分析は残差が正規分布になっていることが望ましいです。

残差とは、モデルから導かれた目的変数と、実際の目的変数の差です。以下のコードでヒストグラムを確認してみましょう。

sns.distplot(result.resid,bins=30)

出力:

若干、右へ伸びていますね。

ちなみに、下のコードでも同じヒストグラムを描けます。

y_pred = result.predict(X_con)  # 予測値を作成

sns.distplot(y-y_pred, bins=30)

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

モデルの結果と実際の結果の散布図を確認してみましょう。

plt.figure(figsize=(5,5),dpi=150)
plt.scatter(y, y_pred, s=15)
plt.xlabel('True', fontsize=20)
plt.ylabel('Pred', fontsize=20)
plt.plot([0, 50], [0, 50], c="g")
plt.xlim(0,52)
plt.ylim(0,52)
plt.show()

出力:

緑の線が予測値と実測値が一致していることを意味しますが、まあまあばらつきがありますね。

scikit-learnの線形回帰

では、scikit-learnを使った回帰分析を行って、結果を比べてみましょう。

from sklearn.linear_model import LinearRegression

reg_model = LinearRegression(fit_intercept=True)
reg_result = reg_model.fit(X, y)

print(" r2: ", reg_model.score(X, y))
print(X.columns)
print(" Coef: ", reg_model.coef_)
print(" intercept: ", reg_model.intercept_)

出力:

あたりまえですが、ほぼ同じような結果になりましたね。

重回帰分析|初心者や文系素人でも直感的に理解できる

【初心者用】Pythonでできる統計手法|ロジスティック回帰分析

【初心者】コピペでできる線形回帰

【初心者】コピペでできるRidge・Lasso回帰

【初心者】コピペでできる決定木回帰

【初心者】コピペでできるランダムフォレス回帰

【初心者】コピペでできるXGBoost回帰



コメント

タイトルとURLをコピーしました