【統計検定準1級】回帰診断

はじめに

この記事では、回帰分析を行うとき、誤差項の仮定が成立しているかどうかを評価する、回帰診断（regression diagnosis）について扱います。統計検定準1級レベルの内容です。

はじめに
誤差項の仮定
回帰診断
Rを使って回帰診断図を描く
回帰診断図の使い道とは？
参考

誤差項の仮定

以下の回帰モデルを仮定します。

$\begin{eqnarray} Y_i &=& α + β X_{i} + u_i \tag{1} \end{eqnarray}$

誤差項の仮定は、

仮定1
$E[u_i]=0$ 　 $, \forall i$
期待値がゼロ

仮定2
$E[u_i|X_i]=0$ 　 $, \forall i$
説明変数 $X_i$ の条件付き期待値がゼロ

仮定3
$V[u_i]=σ^2$ $,\forall i$
分散が等しい（等分散性）

仮定4
$Cov(u_i,u_j)=0$ $,\forall i, i \neq j$ 　
共分散がゼロ（分散共分散行列の対角成分以外はゼロ）

仮定5
$u_i \overset{iid}{\sim} N(0,σ^2)$ $,\forall i$
正規分布に従う（正規性）

これらの仮定を合わせると、誤差項は互いに独立に同一の正規分布に従う（independently and identically distributed）ことが分かります。

回帰診断

①予測値に対する残差のプロット

仮定1 期待値ゼロ、仮定3 等分散性

回帰モデルの線形性が正しい場合には、回帰モデルが観察データに完全に当てはまる（fit）するならば、残差（Residual）は0を中心に一様に分布する。残差が一様ではなく、曲線的な傾向になっている場合は、モデルの線形性の仮定が疑われる。また、残差の0の水平線からのズレが全体的な傾向か、少数の外れ値が原因かを判断することができる。

②残差の正規Q-Qプロット

仮定5 誤差項の正規性

標準化された残差を大きさの順に並べたものの分位点と、標準正規分布の累積分布関数の分位点をプロットしたもの。誤差項が正規分布しているかどうかを判断する。正規性の仮定が正しいとき、残差は傾き1の直線上に並ぶ。直線から外れた値を外れ値として取り除けば、モデルの正規性はおおよそ妥当だと判断できる。

③予測値に対する標準化した残差の絶対値の平方根のプロット

仮定3 等分散性

残差の変動状況を確認するもの。モデルが正しければ、プロットは一様となる。何らかの傾向があればモデルの妥当性や外れ値を疑う。プロットが予測値に対して増加または減少する傾向がある場合は等分散性が成立しないと判断される。

等分散性が成立しない場合に、従属変数や説明変数を変数変換（平方根、自然対数など）することで等分散性に近づく場合がある。

④梃子（てこ）値（leverage）に対する標準化した残差のプロット

外れ値

各観測データが回帰係数の推定値（ $β$ ）へ与える影響度を判断するもの。点線はクック距離（Cook's distance）と呼び、その観測データを含む場合と含まない場合の予測値の変化量を示す。クック距離が大きいと外れ値の可能性があり、0.5を超えると影響力があり、1を超えると異常に影響力があると判断される。

誤差項同士の相関

誤差項同士に相関がある場合、回帰分析によって最良線形推定量を得ることができない。特に時系列データを扱う場合の誤差項の相関は系列相関（serial correlation）と呼ぶ。系列相関を検定するための統計量としてダービン・ワトソン比、系列相関に対処する方法としてコクラン・オーカット法等が知られている。

www.goodnalife.com

Rを使って回帰診断図を描く

Rで回帰診断図（Regression Diagnosis Plots）を描いてみます。Rのデフォルトデータセットである「airquality」を使い、重回帰分析を行います。オゾンの量を、日射量、風力、気温に回帰させます。

オゾンの量（単位：ppb, $Y_i$ ）
日射量（単位：lang, $X_{1i}$ ）
風力（単位：mph, $X_{2i}$ ）
気温（単位：F°, $X_{3i}$ ）

plot関数を使えば回帰診断図を作成できます。

res <- lm(Ozone ~ Solar.R + Wind + Temp, data=airquality)
plot(res)

それでは順番に診断図を見ていきます。

①予測値に対する残差のプロット

曲線的な傾向が見られます。特に観察データ30,62,117が外れ値である可能性が高いことが示唆されます。

②残差の正規Q-Qプロット

観察データ30,62,117などを除けば、概ね回帰モデル（誤差項）の正規性は成立することが分かります。

③予測値に対する標準化した残差の絶対値の平方根のプロット

予測値（fitted value）に対して減少する傾向が見られます。(1)の回帰モデルの説明変数ベクトルを $log(X_i)$ などと変数変換し、等分散に近づけるなどの操作を行うことが考えられます。

④梃子（てこ）値（leverage）に対する標準化した残差のプロット

観測データ9,48のleverageが大きく、回帰係数に影響を与える可能性と、外れ値である可能性があります。しかしいずれもクック距離（Cook's distance）は0.5を超えていません。

【補足】ggplotを使った回帰診断図

plot関数を使えば一瞬で作図はできるものの、自分の勉強のため上述の4枚をggplotを使って作図しました。

4枚目の梃子値と残差のグラフ（Residuals vs Leverage）はうまくクック距離を重ねて表示することができなかったため、残差とクック距離のグラフを描きました。

回帰診断図の使い道とは？

回帰モデルの線形性、誤差項の仮定が妥当かどうかを「診断」しました。それでは診断結果を踏まえてどのような対処ができるのでしょうか。1つはモデルへの当てはまりが悪い外れ値を除外することです。もう1つは、回帰モデルの関数形を見直すことです。両者ともある決まった作法がある訳ではありませんし、特に後者の場合は先行研究における観測データの扱いを踏まえて、妥当なモデルを検討し続けるしかないのだと思います。

回帰診断図より、外れ値である可能性があると考えられる観察データ30,62,117を取り除き、同じ回帰モデルを推定します。