次のデータは、「Step1. 基礎編　27-3. 重回帰分析」で用いたデータです。このデータを使って、年日照時間（）を目的変数とし、年平均気温（）、人口密度（）、持ち家比率（）、降水量（）を説明変数とする重回帰式を求めてみましょう。

ファイルのダウンロードはこちらから

• （日本語）linear_regression2.csv
• （英語）linear_regression2_eng.csv

■Excelによる解析例

1. データ > データ分析 > 回帰分析を選択 > [OK]をクリック

※「データ分析」が見当たらない方はこちらを参考に設定を行ってください。



2. 「入力元」と「出力オプション」を図のように設定

3. 回帰分析の結果が出力される

※この結果を見ると、推定された単回帰式はであることが分かります。




4. 多重共線性の確認

重回帰分析を行う際には、多重共線性を確認する必要があります。多重共線性とは、いくつかの説明変数の間で線形関係が認められる状態を言います。このような場合、偏回帰係数が計算できなかったり、計算できても値の信頼性が低くなります。

ここでは、説明変数間の多重共線性を検出するための指標の1つであるVIF（Variance Inflation Factor）を計算してみます。VIFは説明変数間の相関係数行列の逆行列の対角要素であり、値が大きい場合（よく使われる基準は10以上）はその変数を分析から除いた方がよいと考えられます。

①説明変数の相関行列を求める





②相関行列の空欄を埋める



③相関行列の逆行列を求める

※相関行列と同じ大きさの領域を選択（B8からF12まで）→関数の挿入エリアに「=MINVERSE(B2:F6)」と入力→「control」＋「shift」＋「Enter」を押す




算出された逆行列のうち、赤い文字で示した部分がVIFになります。例えば年平均気温のVIFは1.33です。VIFは全て10以下であることから、除外すべき変数は無いと考えられます。

■Rによる解析例

dat <- read.csv("linear_regression2_eng.csv", header=T, row.names=1)  #データの読み込み

#線形回帰
res <- lm(formula = Daylight_h ~ Temperature + Pop_density + Homeownership_rate + Precipitation, data = dat)  #Daylight_hを目的変数、それ以外の変数を説明変数、読み込むデータをdatとして線形回帰を行う
summary(res)  #線形回帰の結果を表示

#VIF
library(car)
vif(res)  #説明変数のVIFを出力

#回帰診断
par(mfrow = c(2, 2))  #2x2の図を表示するように設定
plot(res)  #回帰診断の結果を表示

par(mfrow=c(1,1))  #1x1の図を表示するように設定（初期設定に戻す）

結果の見方

• Residuals：残差の分布
• Coefficients：回帰分析により推定された偏回帰係数、標準誤差、t値、P値
• Residual standard error：残差の標準誤差（ばらつき）と自由度
• Multiple R-squared：決定係数
• Adjusted R-squared：自由度調整済み決定係数
• F-statistic：回帰式の有意性の検定結果（分散分析の結果；F値、自由度、P値）

結果の見方

• Residuals vs Fitted：x軸は回帰直線により予測されたyの値、y軸は予測値と観測値との差（残差）
• Normal Q-Q：残差のQ-Qプロット（残差の正規性を確認するためのプロット）
• Scale-Location：x軸は回帰直線により予測されたyの値、y軸は標準化した残差の絶対値の平方根を取った値
• Residuals vs Leverage：Cookの距離（Cookの距離が0.5を超えると回帰式への影響が大きい＝外れ値の可能性が高い）

■Pythonによる解析例

import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

dat = pd.read_csv("linear_regression2_eng.csv", index_col=0)  #データの読み込み

#線形回帰
x = dat[["Temperature", "Pop_density", "Homeownership_rate", "Precipitation"]]  #4つの変数を説明変数とする
x2 = sm.add_constant(x)  #切片を計算するために説明変数のデータの中に全て1からなる列を追加
y = dat["Daylight_h"]  #Daylight_hを目的変数とする
res = sm.OLS(y, x2).fit()  #y、x2を用いて線形回帰を行う
print(res.summary())  #線形回帰の結果を表示

#VIF
from statsmodels.stats.outliers_influence import variance_inflation_factor

vif = pd.DataFrame()  #結果を入れるためのDataFrameを作成
vif["VIF"] = [variance_inflation_factor(x2.values, i) for i in range(1, x2.shape[1])]  #VIFを算出
vif["Variable"] = x2.columns[1:5]  #各VIFに変数名を付与
print(vif)  #VIFを表示

#回帰診断
plt.plot(res.predict(x2), y - res.predict(x2), "o")  #yの予測値と残差の値をプロット     
plt.hlines(y = 0, xmin = 1500, xmax = 2500, color = 'red')  #y = 0の直線を重ねる
plt.xlabel('Fitted values')  #x軸ラベルを設定                            
plt.ylabel('Residuals')  #y軸ラベルを設定                                   

plt.show() 

結果の見方

• R-squared：決定係数
• Adj. R-squared：自由度調整済み決定係数
• F-statistic、Prob(F-statistic)：回帰式の有意性の検定結果
• Log-Likelihood：対数尤度（対数尤度が大きいほど当てはまりが良いモデル）
• AIC：赤池情報量規準（AICが小さいほど当てはまりが良いモデル）
• BIC：ベイズ情報量基準（対数尤度が小さいほど当てはまりが良いモデル）

• Coef：回帰分析により推定された偏回帰係数
• std err：標準誤差
• t：t値
• P>|t|：P値
• [0.025、0.975]：偏回帰係数の95%信頼区間

• Omnibus、Prob(Omnibus)：残差に対する歪度と尖度によるオムニバス検定（正規性の検定の1つ）とP値
• Skew：残差の歪度
• Kurtosis：残差の尖度
• Durbin-Watson：残差間の自己相関
• Jarque-Bera(JB)、Prob(JB)：残差に対するジャーク=ベラ検定（正規性の検定の1つ）とP値
• Condo. No.：多重共線性の指標