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 を使ったほうが簡単です。
では、ボストンのデータセットを１つのデータフレームの形に整えます。
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」としますので、これを「ｙ」とします。
説明変数は、変数「INDUS」「RM」「PTRATIO」「LSTAT」の４つとします。この説明変数の選択には特に大きな意味はありません。多すぎると煩雑になるため、例として選択したと思って下さい。
一番下のコードは、すべての値が、1.0 である列を作成しています。これは、切片を作るために必要な処理です。重回帰分析では切片が必要です。相当特殊な場合は切片を設けないことがあるかもしれませんが、わたしは知りません。
では、各変数を確認してみましょう。
y
出力：
ｙには「PRICE」が入っています。回帰分析の目的変数は連続変数でしたね。
X
出力：
Ｘには説明変数の４つと、「const」という変数があります。「const」は切片を設定するためのもので、データはすべて 1.0 になっています。
機械学習では予測を重視しますので、学習データとテストデータに分割して、学習データのみを回帰分析します。
でも、統計学ではデータの説明を重視するため、すべてのデータを使って回帰分析します。ここは、用途に合わせて分割するしないを判断して下さい。
重回帰分析の実装
model = sm.OLS(y, X_con) # モデルを作成 result = model.fit() # モデルを適応
以上で、重回帰分析ができます。非常に簡単ですね。
１行目で目的変数と説明変数を指定したモデルを作成します。
２行目でそのモデルを適応させて、「result」としています。
結果の確認
result.summary() # 結果を確認
出力：
いっぱい出力結果がありますが、基本的で最低限のものだけ解説します。
Ⓐは、自由度調整済決定係数です。「R²」ですね。これは、１に近いほど良いです。解釈としては、今回の重回帰分析のモデルが実際のデータをどの程度説明できるかを表しています。
Ⓒは、編回帰係数です。Ⓑはそれに対する説明変数です。「INDUS」の編回帰係数は、0.0076 ですね。編回帰係数は、「他の説明変数が一定であった場合に、その変数を１変化させると、目的変数がどの程度変化するか」を表しています。
PRICE = (0.0076×INDUS) + (4.5152×RM) + (-0.9351×PTRATIO) + (-0.5757×LSTAT) +18.6150（切片）
このような計算式で表すことができます。
ここの注意が必要です。先ほどの編回帰係数の説明で、「他の説明変数が一定であった場合」といいましたが、これが曲者です。他の説明変数が一定であり、１つの説明変数のみが変化することがありえるでしょうか。これは、各変数の相関が０ということです。それはあり得ないですね。
この説明変数間の関係性の強さが多重共線性です。多重共線性はあるないではなく、強いか弱いです。たとえば、相関係数が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（回帰係数）の数値だけ、変化します。
標準化は、平均値が０、標準偏差が１にすることなので、説明変数間で比較することができます。
ここでは、「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でできる統計手法｜ロジスティック回帰分析
コメント