重回帰分析
©2023 Nobuyuki Otake
最終作成日:2024-02-04
1 重回帰分析
重回帰分析は逆行列\((X^tX)^{-1}\)が存在すれば，正規方程式の解\(\hat{β}\)，Yの予測値，残差を求めることができる。
\(\hat{β}=(X^tX)^{-1}X^tY\) # 正規方程式の解
\(\hat{Y}=X\hat{β}=X(X^tX)^{-1}X^tY\) # Yの予測値
\(e=Y_i-\hat{Y}_i\) # 残差
1.1 実際に計算してみる
library(MASS)
data(Cars93)
c <- rep(1,93) # 定数:1のベクトル
x1 <- Cars93$MPG.city # 説明変数 MPG.city
x2 <- Cars93$EngineSize # 説明変数 EngineSize
x3 <- Cars93$Horsepower # 説明変数 Horsepower
y <- Cars93$Price # 目的変数 Price
x <- cbind(c, x1,x2,x3) # 定数，説明変数をまとめておく
# 正規方程式の解を求める:solve()逆行列を求める，上付きtは転置行列 %*%は行列の積
solve(t(x) %*% x) %*% t(x) %*% y ## [,1]
## c 5.6031419
## x1 -0.2115335
## x2 -0.1288022
## x3 0.1319717head(x %*% solve(t(x) %*% x) %*% t(x) %*% y) # 予測値:head()で最初の6サンプル分## [,1]
## [1,] 18.55900
## [2,] 27.77771
## [3,] 23.71096
## [4,] 23.92249
## [5,] 27.94871
## [6,] 15.18293head(y - x %*% solve(t(x) %*% x) %*% t(x) %*% y) # 残差:head()で最初の6サンプル分## [,1]
## [1,] -2.658997
## [2,] 6.122291
## [3,] 5.389044
## [4,] 13.777511
## [5,] 2.051292
## [6,] 0.5170741.2 決定係数を計算
上記重回帰式の精度を判定する時いくつかの指標があるが，ここでは決定係数により判断してみる. Yの全変動（SST），回帰変動（SSR），残差変動（SSE）を以下のように定義する．
\(SST=Σ(Y_i-\bar{Y})^2\)
\(SSR=∑(\hat{Y}_i-\bar{Y})^2\)
\(SSE=∑e_i^2\)
上記は次の関係を満たす．\(SST=SSR+SSE\)
決定係数は次のように定義される．\(R^2=1-\frac{SSR}{SSE}\)
自由度調整済み決定係数は次のように定義される．\(R^{*2}=1-\frac{n-1}{n-p-1}(1-R^2)\)
# 決定係数
sst <- sum((y-mean(y))^2)
ssr <- sum(((x %*% solve(t(x) %*% x) %*% t(x) %*% y)-mean(y))^2)
sse <- sum((y - x %*% solve(t(x) %*% x) %*% t(x) %*% y)^2)
R2 <- 1 - (sse / sst)
R2## [1] 0.6289284# 自由度調整済み決定係数
aR2 <- 1-(((93-1)/(93-3-1))*(1-R2)) # Cars93はn=93
aR2## [1] 0.6164204# 重相関係数；重相関係数は決定係数の平方根
sqrt(R2)## [1] 0.79305012 重回帰分析（Rで解析）
2.1 data読込
実例で使用するデータ：MASSライブラリーのCars93
str()でサンプルサイズ，アイテム数，各アイテムのデータの型，Factor(因子型)の場合カテゴリー数を確認する.
colnames()で列名と列番号を確認する．分析実行時に列番号が分かると便利
library(MASS)
data(Cars93)
str(Cars93)## 'data.frame': 93 obs. of 27 variables:
## $ Manufacturer : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
## $ Model : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
## $ Type : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
## $ Min.Price : num 12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
## $ Price : num 15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
## $ Max.Price : num 18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
## $ MPG.city : int 25 18 20 19 22 22 19 16 19 16 ...
## $ MPG.highway : int 31 25 26 26 30 31 28 25 27 25 ...
## $ AirBags : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
## $ DriveTrain : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
## $ Cylinders : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
## $ EngineSize : num 1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
## $ Horsepower : int 140 200 172 172 208 110 170 180 170 200 ...
## $ RPM : int 6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
## $ Rev.per.mile : int 2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
## $ Man.trans.avail : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
## $ Fuel.tank.capacity: num 13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
## $ Passengers : int 5 5 5 6 4 6 6 6 5 6 ...
## $ Length : int 177 195 180 193 186 189 200 216 198 206 ...
## $ Wheelbase : int 102 115 102 106 109 105 111 116 108 114 ...
## $ Width : int 68 71 67 70 69 69 74 78 73 73 ...
## $ Turn.circle : int 37 38 37 37 39 41 42 45 41 43 ...
## $ Rear.seat.room : num 26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
## $ Luggage.room : int 11 15 14 17 13 16 17 21 14 18 ...
## $ Weight : int 2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
## $ Origin : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
## $ Make : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...colnames(Cars93)## [1] "Manufacturer" "Model" "Type"
## [4] "Min.Price" "Price" "Max.Price"
## [7] "MPG.city" "MPG.highway" "AirBags"
## [10] "DriveTrain" "Cylinders" "EngineSize"
## [13] "Horsepower" "RPM" "Rev.per.mile"
## [16] "Man.trans.avail" "Fuel.tank.capacity" "Passengers"
## [19] "Length" "Wheelbase" "Width"
## [22] "Turn.circle" "Rear.seat.room" "Luggage.room"
## [25] "Weight" "Origin" "Make"3 基礎分析
3.1 基本統計量
summary()で5数要約,平均値を確認
describe()で基本統計量を確認:Rは不偏分散を返す
library(psych)
summary(Cars93[,c(5 ,7 ,12 ,19)]) #アイテムを5,7,12,19列に絞る## Price MPG.city EngineSize Length
## Min. : 7.40 Min. :15.00 Min. :1.000 Min. :141.0
## 1st Qu.:12.20 1st Qu.:18.00 1st Qu.:1.800 1st Qu.:174.0
## Median :17.70 Median :21.00 Median :2.400 Median :183.0
## Mean :19.51 Mean :22.37 Mean :2.668 Mean :183.2
## 3rd Qu.:23.30 3rd Qu.:25.00 3rd Qu.:3.300 3rd Qu.:192.0
## Max. :61.90 Max. :46.00 Max. :5.700 Max. :219.0describe(Cars93[,c(5 ,7 ,12 ,19)])## vars n mean sd median trimmed mad min max range skew
## Price 1 93 19.51 9.66 17.7 18.29 8.30 7.4 61.9 54.5 1.48
## MPG.city 2 93 22.37 5.62 21.0 21.61 4.45 15.0 46.0 31.0 1.65
## EngineSize 3 93 2.67 1.04 2.4 2.56 0.89 1.0 5.7 4.7 0.83
## Length 4 93 183.20 14.60 183.0 183.27 13.34 141.0 219.0 78.0 -0.09
## kurtosis se
## Price 3.05 1.00
## MPG.city 3.58 0.58
## EngineSize 0.23 0.11
## Length 0.29 1.51trimmed;トリム平均(Rでは10％ずつ外している)
mean(dat.x$test.a , trim = 0.1) ; mean関数ではtrim=で片側10%ずつ外すことができる．
MAD(平均絶対偏差;mean absolute deviation)：Rでは，MAD関数は1.4826倍している．正規分布の場合に標準偏差と比較しやすくするために，MADの定義を1.4826倍しておくことがある．Rの mad(x) で求められる値は1.4826倍されている(補正しないMADは0.8)
mad(dat.x$test．a)
1.4826*median(abs(dat.x\(test.a-median(dat.x\)test.A)))
標準偏差のほうが数学的処理が便利であるのに加えて、平均値が「偏差の絶対値の総和」ではなく「偏差の２乗和」を最小にする定数 （偏差の絶対値の総和を最小にする定数は中央値）
3.2 データの視覚化:散布図行列
pairs.panels(Cars93[,c(5 ,7 ,12 ,19)],
cor = T ,stars = T , #上側三角に相関係数,有意水準(星)
hist.col = "red", #ヒストグラムを赤
rug = F, #ヒストグラムの床なし
ellipses = F, #相関楕円なし
lm = T, #回帰直線を引く
density = F, #確率密度曲線なし
cex.main = 1.0 , main = "散布図行列") 
Price,MPG.city,EngineSizeの3アイテムは若干右に歪んだ分布であり，外れ値とみられる値がある．
4 重回帰分析:lm()関数
引数(y~x,data名)
s.out <- lm(Price ~ MPG.city + EngineSize + Horsepower , data = Cars93)
summary(s.out)##
## Call:
## lm(formula = Price ~ MPG.city + EngineSize + Horsepower, data = Cars93)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.201 -2.765 -1.117 1.753 32.090
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.60314 5.99706 0.934 0.353
## MPG.city -0.21153 0.16627 -1.272 0.207
## EngineSize -0.12880 0.97850 -0.132 0.896
## Horsepower 0.13197 0.01844 7.155 2.26e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.982 on 89 degrees of freedom
## Multiple R-squared: 0.6289, Adjusted R-squared: 0.6164
## F-statistic: 50.28 on 3 and 89 DF, p-value: < 2.2e-16call:は分析したモデル
Residuals:は残差の5数要約
Coefficients:は係数
Estimateは偏回帰係数，Std.Errorは標準誤差，t valueはt値(帰無仮説:母偏回帰係数の値は0である)，Prは有意水準(Horsepowerだけが有意である)
Multiple R-squaredは決定係数，Adjusted R-squaredは自由度調整済み決定係数であり，0.62と目的変数（従属変数）の分散の62％が説明されている．
最終行は分散分析の結果である(帰無仮説:すべての母偏回帰係数の値は0である．対立仮説：いずれかの母偏回帰係数の値は0ではない)．p-valueは0であり帰無仮説は棄却される．
library(pander)
pander(summary(s.out)) # pander()関数で出力を変える| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 5.603 | 5.997 | 0.9343 | 0.3527 |
| MPG.city | -0.2115 | 0.1663 | -1.272 | 0.2066 |
| EngineSize | -0.1288 | 0.9785 | -0.1316 | 0.8956 |
| Horsepower | 0.132 | 0.01844 | 7.155 | 2.262e-10 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 93 | 5.982 | 0.6289 | 0.6164 |
library(equatiomatic) ### 以下はRのコードではない ###
extract_eq(s.out) #理論モデル\[ \operatorname{Price} = \alpha + \beta_{1}(\operatorname{MPG.city}) + \beta_{2}(\operatorname{EngineSize}) + \beta_{3}(\operatorname{Horsepower}) + \epsilon \]
extract_eq(s.out, use_coefs = TRUE) #理論モデルに係数を入れる\[ \operatorname{\widehat{Price}} = 5.6 - 0.21(\operatorname{MPG.city}) - 0.13(\operatorname{EngineSize}) + 0.13(\operatorname{Horsepower}) \]
names(s.out) #分析結果で何が出力されているか## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"head(s.out$fitted.values) # fitted.valuesは予測値:head()で最初の6サンプル分## 1 2 3 4 5 6
## 18.55900 27.77771 23.71096 23.92249 27.94871 15.18293head(s.out$residuals) # residualsは残差:head()で最初の6サンプル分## 1 2 3 4 5 6
## -2.658997 6.122291 5.389044 13.777511 2.051292 0.5170744.1 多重共線性診断
VIFは,１つの説明変数\(x_i\)を目的変数とし他の説明変数を説明変数とした重回帰分析の重相関係数\(R_i^2\)を用いて以下のように定義される.
\(VIF_i=\frac{1}{1-R_i^2}\)
VIF(variance inflation factor)は多重共線性が生じているか否かを判断する指標である．VIFが10以上であれば多重共線性が認められる(VIFが5以上で多重共線性が疑われる).
# MPG.city EngineSize Horsepowerの相関係数行列を算出し逆行列を求める（対角部分がVIF）
car <- Cars93[c(7, 12, 13)] # MPG.city EngineSize Horsepowerだけのdata行列
solve(cor(car)) # 逆行列を求める## MPG.city EngineSize Horsepower
## MPG.city 2.2443768 1.052306 0.7392348
## EngineSize 1.0523062 2.648557 -1.2312417
## Horsepower 0.7392348 -1.231242 2.3986524いずれの変数も5以下であり，多重共線性は疑われない．
4.1.1 vif()関数
VIFはVIF()関数を使うことで求めることできる.
library(car)
vif(s.out)## MPG.city EngineSize Horsepower
## 2.244377 2.648557 2.3986524.1.2 多重共線性診断:カテゴリカル変数がある場合
AirBogsを説明変数に追加：AirBogsは3カテゴリー
vif(lm(Price ~ MPG.city + EngineSize + Horsepower + AirBags , data = Cars93)) ## GVIF Df GVIF^(1/(2*Df))
## MPG.city 2.245722 1 1.498573
## EngineSize 2.674076 1 1.635260
## Horsepower 2.615469 1 1.617241
## AirBags 1.307617 2 1.069351GVIF：一般化されたVIF
GVIF^(1/(2*Df))：自由度で調整されたVIF
ただし，自由度が2を超えるカテゴリカル変数の場合，自由度で調整されたVIFを使用する．判断は平方根を取る必要があり，\(\sqrt{5}\)（2.2）は多重共線性が疑われ，\(\sqrt{10}\)（3.2）は多重共線性が認められる．自由度で調整されたVIFを2乗して元の経験的判断(VIF5,10以上)と比較することもでもよい．
vif(lm(Price ~ MPG.city + EngineSize + Horsepower + Origin , data = Cars93))## MPG.city EngineSize Horsepower Origin
## 2.260691 3.306490 2.751688 1.366619Originは2カテゴリ：カテゴリ数が2の場合自由度は1，連続量は自由度1である．
4.2 回帰診断
par(mfrow=c(2,2)) #一画面2行2列の作図
plot(s.out)
par(mfrow=c(1,1)) #一画面1行1列に戻す回帰診断のグラフは，
残差プロット（左上）は，誤差が独立で同一の分布に従うという仮定を確認するために，予測値を横軸，残差を縦軸にして描画。誤差が互いに独立で正規分布に従うならば，予測値と残差は独立になる。予測値の大小にかかわらず残差が０の水平軸の周りにランダムにプロットされてれば仮定を満たしていると判断される。 残差プロットみると，予測値が小さくなると残差のばらつきは小さくなる傾向を持っている。58，28は残差が大きいことがわかる。それを除くと水平軸の周りに均一にプロットされている。
規準化残差（Scale LocationPlot）（左下）は，予測値を横軸，規準化残差の絶対値を縦軸にして描画。残差プロットと概ね同じ目的の散布図，誤差分散の均一性はこの図のほうが見やすい。 規準化残差のプロットみると予測値が小さいと規準化残差の絶対値の平方根のちらばりが小さい傾向がある。58，28は規準化残差の絶対値が大きい。（規準化残差の絶対値が2の平方根（1.41）を超えている）
規準化残差の正規QQプロット（右上）は，誤差が正規分布に従うという仮定を確認するためのグラフ。仮定が満たされている場合概ね直線上に並ぶ。58，28は直線上から外れている。
影響プロット（右下）は，横軸にてこ比，縦軸に規準化残差をとったグラフ。さらにクックの距離が0.5と1のところに破線が引かれる。てこ比は各観測値が予測値にどの程度影響しているかを表している。説明変数の数（ｐ）とサンプルサイズ（ｎ）から2(p+1)/n，3(p+1)/nより大ならばその観測値は予測値に影響を与えているということになる。p=3，n=93であるから2(3+1)/93=0.086，3(3+1)/93=0.129となる。この値を超える観測値が6つあることがわかる。クックの距離は，てこ比とスチューデント化残差から各観測値の影響度を測っている。クックの距離が0.5以上ならば回帰係数の推定に影響が大きいと判定できる。28はクックの距離0.5上にある。
回帰診断をみてみると，残差の絶対値ではsample59,28が大きいことが分かる．Normal Q-Q(残差の正規性)をみると，やはり59,28が残差が大きくかつ正規性から外れている．Residuals vs Leverageをみると，赤い点線はクックの距離であり，28は若干大きいことが分かる．（クックの距離:そのsampleが入るか入らないかで回帰係数の値が大きく変化する）
4.3 influencePlot関数による回帰診断
library(car)
influencePlot(s.out, id=TRUE)
## StudRes Hat CookD
## 28 -2.9456438 0.19176624 0.47380898
## 39 1.0418977 0.24935393 0.09006456
## 57 -0.5839244 0.27164645 0.03202897
## 59 6.6756859 0.03830895 0.29796023影響プロットは，横軸にてこ比，縦軸にスチューデント化残差，バブルのブルーの濃さがクックの距離を表すグラフである。てこ比は，説明変数の数（ｐ）とサンプルサイズ（ｎ）から2(p+1)/n，3(p+1)/nより大ならばその観測値は予測値に影響を与えているということになる。p=3，n=93であるから2(3+1)/93=0.086，3(3+1)/93=0.129となり，その位置に破線が引かれている。，縦軸にスチューデント化残差，2と-2の位置に破線が引かれている。ブルーの濃さでクックの距離を表している。 2(3+1)/93=0.086，3(3+1)/93=0.129を超える値を取る観測値が6つある。スチューデント化残差が2，-2を超える観測値が5つある。
4.4 autoplot関数による回帰診断
library("ggfortify")
autoplot(lm(Price ~ Horsepower + Passengers + Wheelbase , data = Cars93),
which = 1:6, ncol = 2,label.size = 3 )
4.4.1 層別回帰診断
autoplot(lm(Price ~ Horsepower + Passengers + Wheelbase , data = Cars93),
which = 1:6, ncol = 2,label.size = 3 ,
data = Cars93, colour = "Origin") # colour = 層別変数
4.5 予測値:predict()関数
Pred.m <- predict(s.out) #推定された回帰モデルでyの値を推定
head(Pred.m)## 1 2 3 4 5 6
## 18.55900 27.77771 23.71096 23.92249 27.94871 15.182934.6 データを標準化して回帰分析
scale()関数を使ってデータを標準化し標準化偏回帰係数を求める．
z.out <- lm(scale(Price) ~ scale(MPG.city) +
scale(EngineSize) + scale(Horsepower) , data = Cars93)
summary(z.out)##
## Call:
## lm(formula = scale(Price) ~ scale(MPG.city) + scale(EngineSize) +
## scale(Horsepower), data = Cars93)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5737 -0.2862 -0.1156 0.1814 3.3222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.861e-17 6.422e-02 0.000 1.000
## scale(MPG.city) -1.231e-01 9.673e-02 -1.272 0.207
## scale(EngineSize) -1.383e-02 1.051e-01 -0.132 0.896
## scale(Horsepower) 7.156e-01 1.000e-01 7.155 2.26e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6193 on 89 degrees of freedom
## Multiple R-squared: 0.6289, Adjusted R-squared: 0.6164
## F-statistic: 50.28 on 3 and 89 DF, p-value: < 2.2e-164.7 標準化偏回帰係数算出:lm.beta()関数
library(QuantPsyc)
lm.beta(s.out)## MPG.city EngineSize Horsepower
## -0.12306919 -0.01383256 0.71556388Rは偏回帰係数を返してくる．標準化偏回帰係数を計算したい場合は,
lm.beta()関数で計算する．ただしdummy変数を含む場合errorになる．
計算結果は，データを標準化して分析した係数と一致している．
次の式で計算をする．\(標準化偏回帰係数=偏回帰係数\sqrt{\frac{独立変数の分散}{従属変数の分散}}\)
5 参考
5.1 Modelの比較
重回帰分析は独立変数(説明変数)を選択する必要がある．選択する場合，統計的に選択する(有意な変数だけを残す)，分析目的を考えて変数を選択する．といった方法が考えられる．
library(memisc)out1 <- lm(Price ~ MPG.city + EngineSize + Horsepower, data = Cars93)
out2 <- update(out1 , ~. + Type)
out3 <- update(out1 , ~. + Origin)
out4 <- update(out1 , ~. + Type + Origin)out <- mtable(out1 , out2 , out3 , out4 ,#モデル名を指定
summary.stats=c("R-squared" ,"F" , "p" , "N",
"AIC" , "BIC" , "Deviance"))
print(out)##
## Calls:
## out1: lm(formula = Price ~ MPG.city + EngineSize + Horsepower, data = Cars93)
## out2: lm(formula = Price ~ MPG.city + EngineSize + Horsepower + Type,
## data = Cars93)
## out3: lm(formula = Price ~ MPG.city + EngineSize + Horsepower + Origin,
## data = Cars93)
## out4: lm(formula = Price ~ MPG.city + EngineSize + Horsepower + Type +
## Origin, data = Cars93)
##
## =============================================================================
## out1 out2 out3 out4
## -----------------------------------------------------------------------------
## (Intercept) 5.603 4.904 3.430 3.878
## (5.997) (6.368) (5.697) (6.041)
## MPG.city -0.212 -0.095 -0.258 -0.173
## (0.166) (0.199) (0.158) (0.190)
## EngineSize -0.129 -0.408 1.457 1.083
## (0.978) (1.178) (1.032) (1.206)
## Horsepower 0.132*** 0.125*** 0.109*** 0.101***
## (0.018) (0.019) (0.019) (0.020)
## Type: Large/Compact 0.376 0.886
## (2.841) (2.696)
## Type: Midsize/Compact 3.734 3.489
## (2.019) (1.914)
## Type: Small/Compact -2.658 -2.419
## (2.267) (2.149)
## Type: Sporty/Compact -2.489 -1.509
## (2.141) (2.051)
## Type: Van/Compact -1.606 -2.384
## (2.649) (2.521)
## Origin: non-USA 4.718*** 4.418**
## (1.370) (1.358)
## -----------------------------------------------------------------------------
## R-squared 0.629 0.680 0.673 0.716
## F 50.282 22.291 45.277 23.253
## p 0.000 0.000 0.000 0.000
## N 93 93 93 93
## AIC 602.556 598.846 592.799 589.679
## BIC 615.219 624.172 607.995 617.538
## Deviance 3185.287 2748.660 2807.011 2437.671
## =============================================================================
## Significance: *** = p < 0.001; ** = p < 0.01; * = p < 0.05モデルの選択;AIC,BIC,Deviance(値が小さい方が良い)
モデルをいくつか考え，それらを分析することでモデルを選択する．
ここでは，Price ~ MPG.city + EngineSize + Horsepowerというベースモデルを考え，そのモデルにType，Originを加えるという4モデルを作成し，解析結果からどのモデルを採用するかを考える．
R-squared(決定係数)，AIC,BIC,Devianceをみるとout4というモデルのfitが最も良い．ここでも，統計指標だけから判断するのではなく，分析目的，解釈などを考慮して選択する必要がある．
5.2 変数選択
5.2.1 step()関数
方法：“both”（ステップワイズ法:変数減増法）, “backward”（変数減少法,） “forward”（変数増加法）
情報規準：“AIC”, “AICc”, “AICb1”, “AICb2”, “BIC”, “KIC”, “KICc”, “KICb1”, “KICb2”
デフォルト：backward，AIC
Cars93から5,6,8,12,13,14,19,21列を抽出しstep関数により分析してみる．
dat <- Cars93[c(5,6,8,12,13,14,19,21)]
step(lm(Price ~. , data = dat), direction = "both")## Start: AIC=89.55
## Price ~ Max.Price + MPG.highway + EngineSize + Horsepower + RPM +
## Length + Width
##
## Df Sum of Sq RSS AIC
## - RPM 1 0.20 205.29 87.64
## - MPG.highway 1 2.73 207.82 88.78
## <none> 205.09 89.55
## - EngineSize 1 5.14 210.23 89.85
## - Horsepower 1 12.50 217.59 93.05
## - Length 1 19.14 224.23 95.85
## - Width 1 21.48 226.57 96.81
## - Max.Price 1 2463.54 2668.63 326.18
##
## Step: AIC=87.64
## Price ~ Max.Price + MPG.highway + EngineSize + Horsepower + Length +
## Width
##
## Df Sum of Sq RSS AIC
## - MPG.highway 1 3.28 208.57 87.11
## <none> 205.29 87.64
## + RPM 1 0.20 205.09 89.55
## - EngineSize 1 11.87 217.16 90.87
## - Length 1 19.13 224.42 93.92
## - Width 1 22.16 227.45 95.17
## - Horsepower 1 23.67 228.96 95.79
## - Max.Price 1 2491.32 2696.61 325.14
##
## Step: AIC=87.11
## Price ~ Max.Price + EngineSize + Horsepower + Length + Width
##
## Df Sum of Sq RSS AIC
## <none> 208.57 87.11
## + MPG.highway 1 3.28 205.29 87.64
## + RPM 1 0.75 207.82 88.78
## - EngineSize 1 12.04 220.61 90.33
## - Length 1 18.71 227.28 93.10
## - Width 1 19.24 227.81 93.32
## - Horsepower 1 26.37 234.94 96.19
## - Max.Price 1 2621.11 2829.68 327.62##
## Call:
## lm(formula = Price ~ Max.Price + EngineSize + Horsepower + Length +
## Width, data = dat)
##
## Coefficients:
## (Intercept) Max.Price EngineSize Horsepower Length Width
## 7.28795 0.75598 0.81268 0.01953 0.05733 -0.28561step(lm(Price ~. , data = dat), direction = "forward")## Start: AIC=89.55
## Price ~ Max.Price + MPG.highway + EngineSize + Horsepower + RPM +
## Length + Width##
## Call:
## lm(formula = Price ~ Max.Price + MPG.highway + EngineSize + Horsepower +
## RPM + Length + Width, data = dat)
##
## Coefficients:
## (Intercept) Max.Price MPG.highway EngineSize Horsepower RPM
## 12.2540190 0.7501384 -0.0471711 0.7134974 0.0206103 -0.0001595
## Length Width
## 0.0591826 -0.3267359変数選択法では，一番最後に出力されているモデルが採用すべきモデルである．
変数増減法と変数増加法では結果が異なっている。変数増加法は有意な変数をひとつずつ追加していくことから，モデルで考慮されない変数がある，一方変数減少法は最初にすべての変数をモデルに取り込み，有意ではない変数をひとつずつ減らしていく。Rのステップワイズ法は変数減増法である．
5.2.2 変数減増法と変数増減法の比較
model.1 <- lm(Price ~ 1 , data = dat) # 定数項だけのモデル
model.all <- lm(Price ~ . , data = dat) # 最大モデル
step(model.1,
scope = list(upper = model.all, lower = model.1 )) # 変数増減法 ## Start: AIC=422.83
## Price ~ 1
##
## Df Sum of Sq RSS AIC
## + Max.Price 1 8270.7 313.3 116.96
## + Horsepower 1 5333.1 3250.9 334.53
## + EngineSize 1 3063.8 5520.2 383.77
## + MPG.highway 1 2698.5 5885.5 389.73
## + Length 1 2177.3 6406.8 397.62
## + Width 1 1785.1 6798.9 403.15
## <none> 8584.0 422.83
## + RPM 1 0.2 8583.8 424.83
##
## Step: AIC=116.96
## Price ~ Max.Price
##
## Df Sum of Sq RSS AIC
## + Horsepower 1 63.6 249.7 97.85
## + EngineSize 1 62.8 250.5 98.15
## + Length 1 50.6 262.7 102.57
## + Width 1 31.3 282.0 109.17
## + MPG.highway 1 26.9 286.4 110.61
## + RPM 1 7.5 305.8 116.71
## <none> 313.3 116.96
## - Max.Price 1 8270.7 8584.0 422.83
##
## Step: AIC=97.85
## Price ~ Max.Price + Horsepower
##
## Df Sum of Sq RSS AIC
## + Length 1 20.12 229.6 92.03
## + EngineSize 1 15.80 233.9 93.77
## + RPM 1 8.71 241.0 96.54
## <none> 249.7 97.85
## + MPG.highway 1 4.62 245.1 98.11
## + Width 1 1.89 247.8 99.14
## - Horsepower 1 63.63 313.3 116.96
## - Max.Price 1 3001.19 3250.9 334.53
##
## Step: AIC=92.03
## Price ~ Max.Price + Horsepower + Length
##
## Df Sum of Sq RSS AIC
## + Width 1 8.96 220.6 90.33
## <none> 229.6 92.03
## + EngineSize 1 1.76 227.8 93.32
## + MPG.highway 1 0.68 228.9 93.76
## + RPM 1 0.31 229.2 93.91
## - Length 1 20.12 249.7 97.85
## - Horsepower 1 33.13 262.7 102.57
## - Max.Price 1 2961.92 3191.5 334.81
##
## Step: AIC=90.33
## Price ~ Max.Price + Horsepower + Length + Width
##
## Df Sum of Sq RSS AIC
## + EngineSize 1 12.04 208.57 87.11
## + RPM 1 8.61 212.00 88.63
## <none> 220.61 90.33
## + MPG.highway 1 3.45 217.16 90.87
## - Width 1 8.96 229.56 92.03
## - Length 1 27.19 247.79 99.14
## - Horsepower 1 42.08 262.69 104.57
## - Max.Price 1 2664.18 2884.79 327.42
##
## Step: AIC=87.11
## Price ~ Max.Price + Horsepower + Length + Width + EngineSize
##
## Df Sum of Sq RSS AIC
## <none> 208.57 87.11
## + MPG.highway 1 3.28 205.29 87.64
## + RPM 1 0.75 207.82 88.78
## - EngineSize 1 12.04 220.61 90.33
## - Length 1 18.71 227.28 93.10
## - Width 1 19.24 227.81 93.32
## - Horsepower 1 26.37 234.94 96.19
## - Max.Price 1 2621.11 2829.68 327.62##
## Call:
## lm(formula = Price ~ Max.Price + Horsepower + Length + Width +
## EngineSize, data = dat)
##
## Coefficients:
## (Intercept) Max.Price Horsepower Length Width EngineSize
## 7.28795 0.75598 0.01953 0.05733 -0.28561 0.81268step(model.all,
scope = list(upper = model.all, lower = model.1 )) # 変数減増法 ## Start: AIC=89.55
## Price ~ Max.Price + MPG.highway + EngineSize + Horsepower + RPM +
## Length + Width
##
## Df Sum of Sq RSS AIC
## - RPM 1 0.20 205.29 87.64
## - MPG.highway 1 2.73 207.82 88.78
## <none> 205.09 89.55
## - EngineSize 1 5.14 210.23 89.85
## - Horsepower 1 12.50 217.59 93.05
## - Length 1 19.14 224.23 95.85
## - Width 1 21.48 226.57 96.81
## - Max.Price 1 2463.54 2668.63 326.18
##
## Step: AIC=87.64
## Price ~ Max.Price + MPG.highway + EngineSize + Horsepower + Length +
## Width
##
## Df Sum of Sq RSS AIC
## - MPG.highway 1 3.28 208.57 87.11
## <none> 205.29 87.64
## + RPM 1 0.20 205.09 89.55
## - EngineSize 1 11.87 217.16 90.87
## - Length 1 19.13 224.42 93.92
## - Width 1 22.16 227.45 95.17
## - Horsepower 1 23.67 228.96 95.79
## - Max.Price 1 2491.32 2696.61 325.14
##
## Step: AIC=87.11
## Price ~ Max.Price + EngineSize + Horsepower + Length + Width
##
## Df Sum of Sq RSS AIC
## <none> 208.57 87.11
## + MPG.highway 1 3.28 205.29 87.64
## + RPM 1 0.75 207.82 88.78
## - EngineSize 1 12.04 220.61 90.33
## - Length 1 18.71 227.28 93.10
## - Width 1 19.24 227.81 93.32
## - Horsepower 1 26.37 234.94 96.19
## - Max.Price 1 2621.11 2829.68 327.62##
## Call:
## lm(formula = Price ~ Max.Price + EngineSize + Horsepower + Length +
## Width, data = dat)
##
## Coefficients:
## (Intercept) Max.Price EngineSize Horsepower Length Width
## 7.28795 0.75598 0.81268 0.01953 0.05733 -0.28561変数増減法はStartのAIC=422.83 変数減増法はStartのAIC=89.55
5.2.3 総当り法
bestglm関数はy(目的変数)を最終列にする必要がある．select関数で並び替え．
情報規準：“AIC”, “BIC”, “BICg”, “BICq”, “LOOCV”, “CV”
library(dplyr)
dat2 <- select(dat, Max.Price,MPG.highway,
EngineSize,Horsepower,RPM,Length,Width,Price,everything())
names(dat2)## [1] "Max.Price" "MPG.highway" "EngineSize" "Horsepower" "RPM"
## [6] "Length" "Width" "Price"library(bestglm)
out.aic <- bestglm(dat2, IC="AIC")
out.aic## AIC
## BICq equivalent for q in (0.509145549018386, 0.821889676136107)
## Best Model:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.28794573 5.340327985 1.364700 1.758667e-01
## Max.Price 0.75598352 0.022863182 33.065543 4.831775e-51
## EngineSize 0.81268176 0.362696498 2.240666 2.759451e-02
## Horsepower 0.01953361 0.005889715 3.316563 1.331000e-03
## Length 0.05733156 0.020522207 2.793635 6.409383e-03
## Width -0.28560725 0.100826905 -2.832649 5.735725e-03out.bic <- bestglm(dat2, IC="BIC")
out.bic## BIC
## BICq equivalent for q in (0.162301641925943, 0.509145549018386)
## Best Model:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.08998008 2.271639455 -3.121085 2.429816e-03
## Max.Price 0.77176240 0.022774764 33.886735 1.186899e-52
## Horsepower 0.01846543 0.005152635 3.583688 5.523013e-04
## Length 0.03844367 0.013763191 2.793223 6.388673e-03