回帰モデルの診断・評価(1)：正規性・等分散性・外れ値・独立性

Rothman先生のModern Epidemiology（4th edition）の拾い読みシリーズとして"Chapter 21：Regression Analysis Part II: Model Fitting and Assessment" をまとめようと思ったんですが、もう少し統計学寄りの内容を勉強しようと思い、Applied linear statistical models（5th edition）を拾い読みしました。

この記事では回帰モデルの診断・評価の話だけを扱っています。回帰モデルの選定については下のリンクを参照のこと。

necostat.hatenablog.jp

線形回帰モデルの前提条件
残差とは
- ggfortifyパッケージで残差プロットを描いてみる
正規性
等分散性
- 対処法
外れ値
- 対処法
独立性
おわりに
参考資料

線形回帰モデルの前提条件

「アウトカムである連続変数Yを説明変数Xで回帰する」という比較的シンプルな線形回帰モデルを数式で考えてみます（ちょっとクドイ書き方かも）。下の式で $i$ は $i$ 番目の対象者であることを示す添字です。

$\begin{aligned} E [Y_i | X=X_i ] &= \beta_0 + \beta_1 X_i \\ Y_i &= E [Y_i | X=X_i ]+ \varepsilon_i \\ \varepsilon_i &\sim Normal(0, \sigma^2) \\ \varepsilon_i &\perp \varepsilon_j \end{aligned}$

1番目の式は、アウトカムYの平均値と説明変数Xの関係性を表しています。変数Xが与えられたときに期待されるYの平均値（条件付き期待値）はXの1次式、つまり直線関係になっています。線形回帰モデルというくらいだから直線でないとダメかというとそんなことはなく、別に2次式だろうが、対数変換されていようが構いません。

線形性とは

2つの数値の和を入力したときの出力は、それぞれの値を入力したときの出力の和になる（ $f(x_1 + x_2) = f(x_1) + f(x_2)$ ）
入力をa倍にすると、出力もa倍になる（ $f(ax) = af(x)$ ）

という性質のことで、これを満たしていると「入力に対してどんな出力が得られるか」が理解しやすくなります。

「係数βについて線形独立である」*1（ベクトルとか行列を習ったときに出てきた用語で、「軸」の方向が平行になっていないこと）という言い方もできると思います。モデル式の右辺が変数Xの多項式になることもありうるので、変数Xではなく係数βに注目した説明の方が昔の私には分かりやすかったです。

2番目の式にあるように、実際に観察されるアウトカムYは説明変数Xと係数βで完全に決まるわけではありません。偶然誤差（random error） $\varepsilon_i$ がノイズとして加わって測定されます。

3番目の式は、この偶然誤差が中心が0の正規分布に従っていると仮定していることを示しています。誤差 $\varepsilon_i$ の平均は0なので、アウトカムの期待値（平均値） $E[Y_i$ ] は $\beta_0 + \beta_1 X_i$ となるように想定したモデルであると言えます。偶然誤差 $\varepsilon_i$ のバラつき具合は分散 $\sigma ^2$ で与えられ、 $\sigma ^2$ が大きいほどアウトカムの期待値 $E[Y_i$ ] からのズレが大きくなりえることになります。そして、各対象者における誤差の分散は等しいことが仮定されます。

4番目の式にある $" \perp "$ は、垂直に交わる記号として見覚えがあるかもしれません。ここでは「無関係である・独立である」ことを示しています。つまり、このモデルでは各対象者の誤差同士には関連がないことが仮定されています。

前提条件をまとめてみます。

線形性（linearity）：係数βが互いに線形独立である
正規性（normality）：誤差が正規分布に従う
等分散性（homoscedasticity）：誤差の分散が一定である
独立性（independence）：誤差が互いに独立である

次からは、誤差についての条件（2-4番目）について評価する方法を見ていきましょう。

残差とは

変数Yの真の値が分からないと誤差（error）がいくらなのか分かりません。そこで、観測値と予測値の差、残差（residual）に注目します（残差 = 観測値 - 予測値）。モデルが前提条件を満たしていれば、残差は誤差の満たすべき条件を満たしているはずだからです。

ggfortifyパッケージで残差プロットを描いてみる

Rでお馴染みのirisデータを使って残差分析の方法をいくつか見てみましょう。次の回帰モデルについて評価してみます。

fit <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
          data = iris)

当てはめたモデルを、ggfortifyパッケージのautoplot( )に渡します（実はplot( )に渡すだけでもいいんですが、一度に見栄え良く作図してくれるのでこうしてます）。

library(ggfortify)
autoplot(fit, which=1:6, label.size = 3)

which：表示したいプロットの番号。残差に関するプロットは1〜6まで。
label.size：グラフ内のラベルの大きさ

ちなみに上記に加え、指定したグループ毎に色を塗り分けることもできます。

autoplot(fit, which=1:6, label.size = 3,
         data = iris,
         colour = "Species")

正規性

誤差が正規分布に従うと仮定しているので、残差も大まかに言って正規分布に従っていなければなりません。

正規分布かどうか大まかに確認する方法の1つはヒストグラム（あるいは密度プロット）です。単峰・左右対称な分布ならOKです。

視覚的に確認するもう1つの方法は、Q-Qプロット（quantile-quantile plot）と呼ばれるものです（グラフ(2)）。

Q-Qプロットは2つの分布の中で同じ分位数（パーセンタイル値など）をプロットしたもので、両者の分布が類似していればプロットは直線上に並びます。例えば下の図はそれぞれ、x軸方向・y軸方向にある2つの分布について5つの分位点（1, 25, 50, 75, 99パーセンタイル）の交点をプロットしたものです。

左側は2つとも正規分布の場合です。5つの分位点の相対的な間隔は両分布でだいたい同じようになっているので、それらの交点はだいたい直線上に並んでいます。これに対して、右側は正規分布からは程遠い分布と正規分布を比較しています。こうなると分位点の相対的な間隔が全然違うので、プロットは直線上に並びません。

等分散性

アウトカムYが大きくなるにつれて誤差も大きくなってしまう状況はしばしば遭遇します。そのため、予測値が変化しても誤差が一定かどうかを確認します。

グラフ(1)は、横軸にYの予測値を、縦軸に残差を取ったグラフです。予測値が変化しても残差の中央（青線）が0付近にあって、上下の広がり具合が一定なら等分散性を満たしていると解釈します。

グラフ(3)は、縦軸を残差の代わりにバラツキ具合で補正した残差*2の絶対値の平方根をとったものです。青線は平均的な残差の大きさを示していて、これが水平なら残差の大きさは一定と言えます。

残差をこのように変換すると残差の中央の傾向はわからなくなりますが、残差そのものの上下の広がりを見るよりもバラつきの変化が分かりやすくなります（多分）。

対処法

等分散の仮定を満たしていなくても係数の点推定値に影響はないが、その推定誤差が変わってきてしまいます。対処法としては、

反復再重み付け最小二乗法（Iteratively Reweighted Least Squares, IRLS）：通常どおり推定した残差をもとに、各対象者への重み付けを計算し、それを使って再度残差を計算し、またそれをもとに重み付けを更新し...という方法
一般化推定方程式（Generalized Estimating Equations, GEE）：IRLSを発展させ、対象者間で誤差の相関関係を考慮した方法

などがあります。