とーけい部

現在、このブログは工事中です。

統計学 統計検定 数学

線形モデルと回帰診断図

統計学

修正：2020/12/01　標準化残差の式が間違っていた(いやまあそういう定義もあるんだけど、R内での定義とは異なっていたことに気づかんかった)ので修正しました。申し訳ないです。

外れ値を含むデータ

以下の仮想データを見てみましょう

set.seed(1)
X <- c( rnorm(n = 10 ,mean = 0 , sd = 1) , 5 , 0.5)
Y <- 3*X + rnorm(n = 12 , mean = 0 , sd = 2)
Y[12] <- Y[12] + 20
plot(X , Y , xlim = c(-2 , 7) , ylim = c(-4 , 26) , col = ifelse(Y<10 , "black" , "red") , pch = 16 , cex =1, main = "仮想データ1")

f:id:hasegawa-yuta289:20201115024731p:plain

これは、以下のデータ

DATA <- data.frame(X,Y)
DATA
##             X          Y
## 1  -0.6264538  1.1442009
## 2   0.1836433  1.3306164
## 3  -0.8356286 -3.7493670
## 4   1.5952808  0.3564426
## 5   0.3295078  3.2383852
## 6  -0.8204684 -2.5512724
## 7   0.4874291  1.4299066
## 8   0.7383247  4.1026465
## 9   0.5757814  3.3697864
## 10 -0.3053884  0.2716375
## 11  5.0000000 16.8379547
## 12  0.5000000 23.0642726

の散布図になっています。
まあパッと見てわかるのは、このデータは線形に分布しているみたいだな～ってことですね。
しかし、明らかにデータの11番目の観察と12番目の観察は、データ全体からみると外れた値、外れ値のように見えます。

この二つのデータ点（赤点で示しました）は、果たして適切なんだろうか、
もしかしたら、データの打ち間違えか何かなんじゃないのか……というのが、今回のテーマになります。

とりあえず線形モデルを当てはめてみよう

上記の仮想データについてとりあえず線形モデルで単回帰してみましょう。
つまり、以下

\displaystyle{ \large \begin{eqnarray} y_i &=& \beta_0 + \beta_1 x_i + \epsilon_i \\ \epsilon_i &\overset{iid}\sim& N(0 , \sigma^2)\\ (i &=& 1,2,\cdots , 12) \end{eqnarray} }

というモデルを考えて、このパラメータ\displaystyle{ \beta _ 0 , \beta _ 1 }を最小二乗法で推定してみます。

このモデルでは、以下の仮定を置いていることに注意してください。

1.各観察iは、互いに独立である(言い換えれば、各誤差\displaystyle{ \epsilon _ i }は独立である)
2.各誤差iは、期待値が0(\displaystyle{ E(\epsilon _ i) = 0 })であり、かつ分散が等しい(\displaystyle{ V(\epsilon _ i) = \sigma^ 2 })
3.各観察iは、正規分布している(言い換えれば、各誤差\displaystyle{ \epsilon _ i }は正規分布に従う)

実際の推定は、Rでlm()関数を使うだけで簡単に終わります。

LM <- lm(data = DATA , Y ~ X)
LM
## 
## Call:
## lm(formula = Y ~ X, data = DATA)
## 
## Coefficients:
## (Intercept)            X  
##       2.389        2.957

この線形モデルの切片は2.389、傾きは2.957と推定されました。

実際のデータに推定された線形モデルを図示してみましょう。

plot(DATA , xlim = c(-2 , 7) , ylim = c(-4 , 26) , col = ifelse(Y<10 , "black" , "red") , pch = 16 , cex =1 , main = "仮想データ1")
abline(LM)

f:id:hasegawa-yuta289:20201115024755p:plain

う～ん、なんかあんまり上手くフィットしてないように見えますね。
観察1~10(黒点)よりも上の方を回帰直線が通過してますね。
あれ？ でも、11番目の観察(右上の赤点)にはよく当てはまっているように見えます。
回帰直線よりもだいぶ上側に位置している12番目の観察にはもう全然フィットしてませんね。

まあこの例の場合はさすがにあからさま過ぎるんで言っちゃいますが、これはデータの12番目の観察に引っ張られることで推定された回帰式が上に引っ張られてしまい、うまくデータにフィットしなくなってしまったわけです。

では、この12番目のデータ点は、果たして異常値として排除してしまって良いのでしょうか？
この判断を、散布図をパッと見て「よし、これがおかしいから外そう！」でやってるんじゃ困りますよね。
人によっては「いやいや、この点ってそんなにおかしいかな？」とか言い出しそうです。
そのような相手に対して、「よし、どう見てもおかしいから排除！」というやり方では、この12番目の観察が異常値として排除するに足るだけの十分な理由を説明できません。

そこで登場するのが回帰診断です。
回帰診断とはデータに対してフィットさせた回帰モデルが、実際どれくらい手元のデータによく当てはまっているのか、どこかおかしい点はないのかなどを調べる方法です。
(どちらかというと、データそのものが、フィットさせた回帰モデルにどれくらい当てはまっているのか、前提をきちんと満たしていそうかどうか、のチェックする方法という呼び方の方が正しいかも)

LMに対する回帰診断

まあ一口に回帰診断とか言っても色々あるんで詳しく知りたい人は別途調べてほしいんですが、今回はRにplot(LM)と打ち込むことで(LMはlm()関数の出力)表示される4枚の回帰診断図について説明していきます。

plot(LM)

f:id:hasegawa-yuta289:20201115024823p:plain

f:id:hasegawa-yuta289:20201115024834p:plain

f:id:hasegawa-yuta289:20201115024847p:plain

f:id:hasegawa-yuta289:20201115024900p:plain

なんか色々出てきました。
一つずつ見ていきましょう。

残差vs当てはめ値

f:id:hasegawa-yuta289:20201115024823p:plain

Residualsってのは残差ですね。
fittedってのは当てはめ値のことです。ちなみに当てはめ値\displaystyle{ \hat{y} }のことを推定値と呼ぶこともありますが、yという値そのものは単に観測された値でしかない(真の値を知ることができないパラメータなどとは違って、真の値はすでにそこに"ある"ものである)ので、まあ個人的には推定値という呼び方は良くないのかなとか思います(まあyhatのことを期待値μの推定としてとらえるなら推定値って呼び方もできるんで、まあどっちでもいいが)

\displaystyle{ 残差_i = 実測値_i - 当てはめ値_i\\ つまり、y_i - \hat{y}_iのこと }
oldpar <- par(no.readonly = TRUE)
par(mfrow=c(2,1));par(oma=c(0,0,0,0));par(mar=c(4,4,1.5,1))

plot(DATA , xlim = c(-2 , 7) , ylim = c(-4 , 26) , col = ifelse(Y<10 , "black" , "red") , pch = 16 , cex =1 , main = "仮想データ1")
text(c("11" , "12") , x = c(5 , 0.5) , y = c(19 , 25))
abline(LM)
plot(LM,1)
par(oldpar)

f:id:hasegawa-yuta289:20201115025020p:plain

この図はまあ、fittedに対して残差がどのくらいか見てるわけですね。
LMにおける仮定では全ての誤差は独立同一にとある正規分布に従うとしていましたので、
この図で見て当てはめ値に対して残差が系統的に変化するようでしたら、そもそもこの1次式による線形モデルを利用すること自体をちょっとやめた方がいいんじゃないの？ということになります。

たとえば以下のような仮想データに対する例を見てみれば明らかでしょう。

set.seed(1)
X2 <- rnorm(n = 10 , mean = 0 , sd = 3)
Y2 <- X2^2 + rnorm(n=10 , mean = 0 , sd = 1)
DATA2 <- data.frame(X2 , Y2)
oldpar <- par(no.readonly = TRUE)
par(mfrow=c(2,1));par(oma=c(0,0,0,0));par(mar=c(4,4,1.5,1))
plot(DATA2 , main = "仮想データ2")
curve(add=T , x^2 , xlim = c(-3 , 6) , col = "blue")
LM2 <- lm(data = DATA2 , Y2 ~ X2)
abline(LM2)
plot(LM2 , 1)
par(oldpar)

f:id:hasegawa-yuta289:20201115025047p:plain

実際のデータ点が当てはめ値（線形式）よりも上に位置していれば、残差は正。
実際のデータ点が当てはめ値（線形式）よりも下に位置していれば、残差は負。

この仮想データでは、データが中心から離れれば離れるほど残差が正に大きくなっていく傾向が見て取れ、逆にデータが中心に近づけば近づくほど残差が負に大きくなっていく傾向がありますね。

このように、残差に傾向が存在するかしないか？を残差vs当てはめ値プロットでは確認できます。
この図で残差に何かしらの傾向が存在するならば、線形モデルの仮定である「誤差分布の独立同一性」が非常に怪しいですよね。
まあ、このような使い方をします。

正規QQプロット

QQプロットって知ってますかね？知ってるならこの項目は飛ばして良いです。

QQプロットとは、ざっくり言うと「確率変数Aの分布」と、「確率変数Bの分布」がどれくらい違うのかをチェックする方法です。

具体的には、確率変数Aの標本から求めた分位点と、Bの標本から求めた分位点とを、各々対応している分位どうしでプロットしていきます。（まあ要するに分位点の散布図がQQプロットになります）
二つの分布が等しいなら、各分位点は、(当たり前ですが)同じ値になりますよね？
なので、二つの分布が等しいならQQプロットは理論上は完全な45度線(y=x上の点)をとります。実際のデータにはばらつきがあるため、二つの分布が等しいならおおむね直線を描きます。

set.seed(1)
A <- rnorm(n=10000)
B <- rnorm(n=10000)
qqplot(A,B , cex = 0.5 , main = "QQプロット")
abline(0,1 , col="red" , lwd = 2)

f:id:hasegawa-yuta289:20201115025115p:plain

正規QQプロットとは、このQQプロットを利用して、「片方は正規分布」と仮定し、もう片方は実際のデータを利用することで、理論上の正規分布の分位点の値からどれだけずれているのかをチェックします。

ある確率変数Zがあったとします。

set.seed(1)
Z <- rt(n = 60 , df = 3 )
hist(Z , main = "Z" , xlab = "Z")

f:id:hasegawa-yuta289:20201115025139p:plain

う～ん、このZは正規分布してるんですかね？
正規QQプロットでチェックしてみます。

qqnorm(Z)
abline(0,1)

f:id:hasegawa-yuta289:20201115025205p:plain

あらあら、中心部はともかく、端っこの方は随分と直線関係から外れてしまいました。
これはZが実際には正規分布ではなくt分布（自由度3）に従っていることによります。

正規分布N(20,100)に従う確率変数Vを考えてみましょう。

set.seed(1)
V <- rnorm(n=1000 , mean = 20 , sd = 10)
qqnorm(V)

f:id:hasegawa-yuta289:20201115025222p:plain

正規QQプロットでチェックするのは標準正規分布との当てはまりなので、今回これは45度線にはなっていませんが(縦軸と横軸の縮尺に注意！)、正規分布は位置パラメータ（平均）と尺度パラメータ(分散)しかもたないので、分布の形自体は変わらないために直線形状にはなってくれます。正規分布に対してQQプロットが上手くいく理由がわかりますね。
ちなみに縦軸の単位が大きいのは、確率変数Vの分散が標準正規分布(横軸)と比べて大きいからです。

最後に、ガンマ分布Ga(2,5)に従う確率変数Dと、ガンマ分布Ga(6,8)に従う確率変数Eを考えてみましょう。
Ga分布の最初のパラメータは形状パラメータとよばれ、この値が変わると分布の形状が大きく変化してしまいます。

curve(dgamma(x , shape = 2 , scale = 5) , xlim = c(0 , 100) , col = "red" , main = "ガンマ分布")
curve(dgamma(x , shape = 6 , scale = 8) , xlim = c(0 , 100) , col = "blue" , add = T)
legend("topright" , legend = c("形状パラメータ=2","形状パラメータ=6") , col = c("red" , "blue") , lty = 1 , cex = 1.5 )

f:id:hasegawa-yuta289:20201115025242p:plain

DとEを1000個発生させて、QQプロットを見ましょう。

set.seed(1)
D <- rgamma(n = 1000 , shape = 2 , scale = 5)
E <- rgamma(n = 1000 , shape = 6 , scale = 8)
qqplot(D,E,main="QQプロット")

f:id:hasegawa-yuta289:20201115025255p:plain

これは直線関係にはなってないですよね。
QQプロットは一応どのような分布にも適用できますが、ガンマ分布のようにパラメータによってそもそも形状が変わってしまうような分布について適用しても、(形状パラメータが一致していない限り)あまり意味がなさそうですね。

ちなみに、形状パラメータを合わせてみると以下のようになります(Ga(2,4)とGa(2,8)の比較)。

set.seed(1)
D <- rgamma(n = 1000 , shape = 2 , scale = 4)
E <- rgamma(n = 1000 , shape = 2 , scale = 8)
qqplot(D,E,main="QQプロット")

f:id:hasegawa-yuta289:20201115025312p:plain

まあ～おおむね直線関係になりました。

先ほどのLMの回帰診断図のQQプロットに戻ってみましょう。

plot(LM,2)

f:id:hasegawa-yuta289:20201115025331p:plain

縦軸を見ればわかりますが、これは横軸の理論上の正規分布に対して、縦軸を残差としています。
正規分布に対する独立同一性の仮定が正しければ、残差は(おおむね)正規分布するのでこれは直線を描くはずですね。実際ほぼ直線を描いていますが、やはり観察12だけはおかしな値になっています。こいつだけ正規分布に従っていなさそう、と考えられます。

標準化残差のプロット

standardized residualっていうのは文字通り標準化残差です。

標準化っていうのは標準正規分布を考えてもらえばわかると思います。N(5,100)に従うような確率変数Xに対して、これを標準化したものは(X-5)/√100 ですが、これは標準正規分布に従いますよね。

標準化残差というのは、残差を標準化したものになります。具体的には以下の計算式によって求めた残差の変換値のことを言います。

\displaystyle{ \Large{ 線形モデル\ y = X\beta + \epsilon\ に対して、\\ 標準化残差_i = \frac{残差_i}{\sqrt{\hat{\sigma}^2 (1-h_{ii})}} \\ ただし、\hat{\sigma}^2 = \frac{\sum_i( y_i - \hat{y_i})^2}{n-rankX}であり、\\ h_{ii}はi番目の観察のてこ比である } }

この式の\displaystyle{ \hat{\sigma}^ 2 }は、誤差\displaystyle{ \epsilon }の分散の不偏推定量になってます。
なんで？て思った人は、まあ例えば『自然科学の統計学』の2章とか読めば載ってるんで読んでみてください。
てこ比については次の節で紹介しているので、ここではまあ雰囲気で流してください。

なお、この値のことを標準化残差ではなくスチューデント化残差と呼ぶこともあります。その場合、このスチューデント化残差はインターナルなスチューデント化残差(internally studentized residual)とよばれます。Rで標準化残差を算出するための関数rstandard(LM)は上の標準化残差、もしくはインターナルなスチューデント化残差を計算しています。ちなみに、エクスターナルなスチューデント化残差(あるいは、単にスチューデント化残差)はRの関数rstudent(LM)によって計算できます。参考:Studentized residual

正規分布というものは大体の値が3σ範囲内に入るので、この標準化残差が3を超えているような値は何かおかしいんじゃねえの～？ということが言えるわけです。
(正規分布の仮定が正しければ標準化残差は漸近的に標準正規分布に従うので)

もう一度LMの回帰診断図(標準化残差プロット)を見てみましょう。

plot(LM,3)

f:id:hasegawa-yuta289:20201115025429p:plain

理由はよくわかりませんが、標準化残差の絶対値のルートを縦軸に取ってますね。横軸は当てはめ値です。
√3=1.73205ぐらいらしいので、まあ1.5を超えているような場合には要注意でしょう。これを見るとやはり観察12はおかしな値をとっているように見えます。

てこ比プロットとクックの距離

最後の図はてこ比プロットです。もう一度見てみましょう。

plot(LM,5)

f:id:hasegawa-yuta289:20201115025453p:plain

まあ縦軸は良いですよね。これは標準化残差です。
横軸のLeverage（レバレッジ：てこ比）とは何でしょうか？

レバレッジについて説明する前に、まず思い出してほしいのですが、最小二乗推定量というのはyの線形推定量、つまりyの各々の要素を線形結合したような推定量になってました。(知らなかった人は豆知識として覚えておくとお得です)

めんどくさいので行列表現してしまいますが、以下のようになっています。

\displaystyle{ 線形モデル\ y = X\beta + \epsilon に対して、\\ 最小二乗推定量\hat{\beta}は、以下の正規方程式\\ X'X\beta = X'y\\ の解であるから、Xが正則とすれば、\\ \hat{\beta} = (X'X)^{-1}X'y\ より、\\ 当てはめ値\ \hat{y} = X\hat{\beta}=X(X'X)^{-1}X'y\\ したがって、\hat{y}の各要素は、\\yの各々の要素の線形結合になっている。 }

ということは、yのk番目の観察の当てはめ値は、以下のように表現できますよね。

\displaystyle{ \Large \hat{y_k} = h_{k1}y_1+h_{k2}y_2+\cdots+h_{kk}y_k+\cdots+h_{kn}y_n }

これについて、てこ比とは以下のように定義します。

\displaystyle{ \Large k番目の観察のてこ比とは、h_{kk}のことである }

このてこ比(レバレッジ)が横軸の値になっているわけですね。観察11のレバレッジが0.8程度と、他と比べて格段に高いのがてこ比プロットから見て取れます。

てこ比の定義を見ればわかりますが、
てこ比が高いということは、その値がその値自身を説明する能力が高いということになります。
したがって、てこ比が高ければ高いほど、その観察の当てはめ値はデータ全体について、その観察以外の値に対してあまり影響を受けないということになります。
散布図を思い出してみましょう。

plot(DATA , xlim = c(-2 , 7) , ylim = c(-4 , 26) , col = ifelse(Y<10 , "black" , "red") , pch = 16 , cex =1 , main = "仮想データ1")
abline(LM)

f:id:hasegawa-yuta289:20201115025522p:plain

黒点(レバレッジは0~0.2程度)については回帰直線の当てはまりが良くありませんが、右上の赤点(観察11)にはよく当てはまっているように見えます。これは観察11のレバレッジが0.8程度と非常に高い値を取っているからです。

先ほどの説明の通り、レバレッジが高い観察に対しては推定されたモデルの当てはまりが良いので、残差が(標準化残差が)低くなる傾向にあります。
このことを逆手に取ると、同じ標準化残差をもっているとしても、レバレッジの値が高ければその残差が持つ意味は大きいということがわかります。てこ比が高いにも関わらずその観察の標準化残差が高ければ、その値は非常に「おかしい」値ということになります。

ところで、てこ比プロットにて赤い点線はCook's distance(クックの距離)と書いてありますが、これはなんでしょうか？

クックの距離について説明していきます。
簡単に言えば、「観察kのクックの距離」とは、「観察kを引っこ抜いたデータから推定されたあてはめ値たちと、実際の当てはめ値たちとの距離」になります。
クックの距離が高いほど、その観察がモデル推定に対して与えている影響が絶大である、ということになりますね。
まあ詳しくはね、Wikipediaとかあるんでね、読むとね、良い可能性が無きにしもあらずCook's distance

てこ比プロットを見れば、レバレッジが高いほど、小さな標準化残差に対してクックの距離が大きくなっていくことがわかります。まあ、そりゃそうだろうねという感じがします。

実際に計算してみましょう。観察12のクックの距離は0.5より低いくらいとのことなので、観察12を引っこ抜いて線形モデルを推定し、本当にそうなのかチェックしてみます。

DATA3 <- DATA[-12,]           #観察12を抹消する
LM3 <- lm(data = DATA3 , Y~X)
X12.pred.LM <- fitted(LM)  #観察12を含むデータからのあてはめ値
X12.pred.LM3 <- predict(LM3 , newdata = data.frame(X = DATA$X)) #観察12を除くデータからのあてはめ値
sigma2hat <- sum(resid(LM)^2)/(12-2)
sum((X12.pred.LM - X12.pred.LM3)^2) / (2 * sigma2hat)
## [1] 0.4117405

クックの距離は0.41ぐらいとわかりました。

ちなみに、plot(LM)では出てこない図ですが、plot(LM , 4)と打ち込めばもっとわかりやすくクックの距離を出してくれます。

plot(LM , 4)

f:id:hasegawa-yuta289:20201115025546p:plain

まあ露骨に観察12のクックの距離が大きいですね。
この値があまりに大きいと、その値がモデル推定に対して与えている影響が大き過ぎるということなので、その値が本当に分析に含めるべきまっとうなデータなのか？があやしくなります。

なお、Rでクックの距離を計算したければcooks.distance(LM)という関数があるのでそれを使えば一発です。

観察12はおかしいよね

というわけで、今までの回帰診断図を見る限り、まあ観察12はおかしい値に思えますね。逆に、外れ値っぽい観察11は、それほどおかしな値ではありませんでした。

観察12を異常値として排除し、回帰直線を引いてみましょう。

plot(DATA , xlim = c(-2 , 7) , ylim = c(-4 , 26) , col = ifelse(Y<10 , "black" , "red") , pch = 16 , cex =1 , main = "仮想データ1")
abline(LM3)

f:id:hasegawa-yuta289:20201115025558p:plain

当てはまりが直感的にも良くなりました。やったぜ。(おわり)

参考文献

買いましょう

これだけ読めばOK

Generalized Linear Models With Examples in R (Springer Texts in Statistics)

Generalized Linear Models With Examples in R (Springer Texts in Statistics)

Amazon

多変量解析の本だけどてこ比や標準化残差について載ってる

多変量解析法入門 (ライブラリ新数学大系)

多変量解析法入門 (ライブラリ新数学大系)

Amazon

分散の不偏推定量～の下りがわからなかったり、そもそも行列を使って表現した線形モデルについて詳しくない人はこれ読むと良いかも

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

Amazon

そもそも行列？線形代数？なにそれわかんね～～～～って人は(高校数学わかるなら)これとかおすすめかも、ジョルダン標準形とかあんま使わないから最後まで読まんで良いかもしれん

1冊でマスター 大学の線形代数

1冊でマスター 大学の線形代数

Amazon

Rを使った社会科学系の実践的なデータ分析がやってみたい・実データを通してRのコードの書き方を学びたい人はこれ読むと良いかも。日本語訳がかなり微妙なので英語得意なら原著の方が良いかも

社会科学のためのデータ分析入門(上)

社会科学のためのデータ分析入門(上)

Amazon

あと僕はまだ読んだことないんだけど、回帰診断に特化した書籍だとこんなのがあったり

回帰診断 (統計ライブラリー)

回帰診断 (統計ライブラリー)

Amazon