【重回帰モデル】回帰の診断（1）――異常な観測値はないか？
※重回帰モデルについての、小暮厚之先生の講義ノート（2012年）を基にした覚書です。
回帰の診断：回帰モデルを推定した後に・・・
使用したデータ：2004年調査のSCF(Survey of Consumer Finances)から無作為に抽出した有収入の500世帯の中で、生命保険に加入している275世帯のデータ（TermLife.csv）。
分析の目的：生命保険に対する各家計の需要を説明したい。
各家計が加入する生命保険の保険金額（FACE）を説明する。説明変数として、年収（INCOME）、世帯主の教育程度（EDUCATION）、世帯構成人数（NUMHH）を考える。
> # ファイルの読み込み > Term = read.csv("./R/TermLife.csv", header=T) > Term2 = subset(Term, subset=FACE > 0) # Termから保険に加入した世帯（FACE>0）だけを取り出す。 > Term2$LNFACE = with(Term2, log(FACE)) # FACEの対数値 > Term2$LNINCOME = with(Term2, log(INCOME)) # INCOMEの対数値
ここで対数をとっているのは、重回帰モデルが正規分布を前提としたモデルだから。
対数変換することで、だいたいの変数の分布は正規分布に近づく。
> Term1 = subset(Term2, select=c(EDUCATION, LNFACE, LNINCOME, NUMHH)) # EDUCATION, LNFACE, LNINCOME, NUMHHのみからなるデータセットを作成する。 > summary(Term1) EDUCATION LNFACE LNINCOME NUMHH Min. : 2.00 Min. : 6.685 Min. : 5.561 Min. :1.00 1st Qu.:13.00 1st Qu.:10.820 1st Qu.:10.545 1st Qu.:2.00 Median :16.00 Median :11.918 Median :11.082 Median :3.00 Mean :14.52 Mean :11.990 Mean :11.149 Mean :2.96 3rd Qu.:16.50 3rd Qu.:13.288 3rd Qu.:11.695 3rd Qu.:4.00 Max. :17.00 Max. :16.455 Max. :16.118 Max. :9.00 > cor(Term1) # 相関行列 EDUCATION LNFACE LNINCOME NUMHH EDUCATION 1.0000000 0.3828489 0.3427036 -0.0635292 LNFACE 0.3828489 1.0000000 0.4818427 0.2876115 LNINCOME 0.3427036 0.4818427 1.0000000 0.1793354 NUMHH -0.0635292 0.2876115 0.1793354 1.0000000 > Term1.kaiki = lm(LNFACE~LNINCOME+EDUCATION+NUMHH, data=Term1) > summary(Term1.kaiki) Call: lm(formula = LNFACE ~ LNINCOME + EDUCATION + NUMHH, data = Term1) Residuals（残差）: Min 1Q Median 3Q Max -5.7420 -0.8681 0.0549 0.9093 4.7187 Coefficients: Estimate Std. Error t value Pr(>|t|) （回帰係数） （標準誤差） （t値） （p値） (Intercept) 2.58408 0.84643 3.053 0.00249 ** LNINCOME 0.49353 0.07754 6.365 8.32e-10 *** EDUCATION 0.20641 0.03883 5.316 2.22e-07 *** NUMHH 0.30605 0.06333 4.833 2.26e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.525 on 271 degrees of freedom Multiple R-squared（R2乗）: 0.3425, Adjusted R-squared（自由度調整済みR2乗）: 0.3353 F-statistic（F統計量）: 47.07 on 3 and 271 DF, p-value（F統計量に対するp値）: < 2.2e-16
自由度調整済みR2乗値が0.335というのは、小さいのかどうか（回帰モデルが妥当かどうか）は、よくわからない経験則ではなく、F統計量から判断する。
判断方法は、「R二乗値なんて信仰にすぎない！」について（F統計量のお話）に記述してあります。
重回帰分析の仮定
- 誤差項は平均がゼロ、分散が一定値
の正規分布に互いに独立に従う
- 誤差項と各説明変数は無相関
これらは、以下の5つの仮定に分けられる
これらの仮定は、残差を見ることによって事後的に正否を判断できる
> # 残差の検討：外れ値 > e = resid(Term1.kaiki) # 残差 > se = sum(e^2)/271 # 残差標準誤差 > e.std = e/se # 標準化残差 > yhat = predict(Term1.kaiki) # 回帰予測値 > > # 回帰予測値への標準化残差のプロット > plot(yhat, e.std) > abline(h=c(-1.96, 1.96), col="red")
誤差分散が不均一であると、誤差の標準誤差が誤って評価される。その結果、信頼区間や検定が信頼できないものとなる。
赤線の外側にあるものが外れ値。
> # 外れ値の検出 > names(e) # eの各値の名前はFACE<=0のものが抜けた飛び飛びの番号になっている > names(e) = 1:length(e) # eの各値の名前を"1から271"に変更 > e[order(abs(e))] 241 229 172 4.704770582 4.718672358 -5.741985758 > # abs(e)はeの絶対値 > # order(abs(e))は絶対値でのeの順番 > # e[order(abs(e))]は絶対値の順番でeを並び替え
異常値の検出
「今あるデータから主張した」というのではなく、違う場面でも当てはまるものであるべき。ロバストネス大事。異常値によって引っ張られた結果、回帰係数が偶然有意だった、とかいうのはその場限りの頑健性のかけらもない主張。
2つの外れ方
- 外れ値（utlier）：残差が大きい観測値→被説明変数の乖離を表す
- 影響値（influential observation）：レバレジ（影響力）の大きい観測値→説明変数の乖離を表す
レバレジ（影響力）
> library(car) 要求されたパッケージ MASS をロード中です 要求されたパッケージ nnet をロード中です > influencePlot(Term1.kaiki, id.method="identify")
- ハット値は、説明変数の外れ具合
- ステューデント化残差は、被説明変数の外れ具合
- 円の面積はクックの距離であり、2つの外れ具合両方を表すもの
※ここからは、疲れたので端折ってます。後日追記します（多分…）。
クックの距離
クックの距離は、被説明変数と説明変数の両方の外れ具合の指標。
> plot(Term1.kaiki, 4)
> cooks.distance(Term1.kaiki) #クックの距離の数値を取り出す 352 355 356 357 360 361 1.264793e-05 3.727564e-04 1.659352e-03 3.144669e-03 3.005157e-01 9.728116e-06
360番目に対応する観測値は取り除いたほうがよさげなことが分かる。
外れ値の検出
> x = Term1 > rownames(x) = 1:nrow(x) # 行名を"1から275"に変更 > x[rownames(Term1)=="360",] EDUCATION LNFACE LNINCOME NUMHH 197 17 14.57163 5.560682 4 > # 197番目の観測値が外れ値の候補となった！
ダミー変数による外れ値の処理
> x = rep(0, 275) > x[197] = 1 > Term1$d197 = with(Term1, x) # 197番目以外は0のベクトル > Term1.kaiki2 = lm(LNFACE~EDUCATION+LNINCOME+NUMHH+d197, data=Term1) > summary(Term1.kaiki2) Call: lm(formula = LNFACE ~ EDUCATION + LNINCOME + NUMHH + d197, data = Term1) Residuals: Min 1Q Median 3Q Max -5.8399 -0.8279 0.0302 0.8748 4.9129 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.04945 0.84925 2.413 0.01648 * EDUCATION 0.18433 0.03882 4.749 3.33e-06 *** LNINCOME 0.57489 0.08043 7.148 8.25e-12 *** NUMHH 0.28237 0.06272 4.502 1.00e-05 *** d197 5.06235 1.58934 3.185 0.00162 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.5 on 270 degrees of freedom Multiple R-squared: 0.3664, Adjusted R-squared: 0.357 F-statistic: 39.03 on 4 and 270 DF, p-value: < 2.2e-16
d197のt値が有意ということは、197番目の観測値がほかの観測値とは異なる（=外れ値）であるということを意味する！
・・・というわけで、回帰診断その1はここまでです。
