【R】回帰診断図(Regression Diagnosis Plots)をggplot2で表示する

こんにちは，@PKです.

最近，統計の勉強をRでやり直しています．
Rではlm関数を使うと，簡単に線形回帰分析を行うことができます．
lm function | R Documentation
このlm関数で回帰分析を行った際に，回帰モデルに加えて，残差を視覚的に評価するための回帰診断図(Regression Diagnosis Plots)を，plot(lm)で簡単に表示することができます．
この回帰診断図をきれいな図が描けるggplot2で表示させたいと思います．

サンプルデータと線形回帰分析

サンプルデータや解析の流れはこちらのサイトを参考にさせていただいています．
データ解析・マイニングとR言語

まず，身長と体重のデータを作成し，ggplot2で散布図と回帰直線を引きます．

> library(ggplot2)
> taikei = matrix(0,10,2) 
> taikei[,1]<-c(165,170,172,175,170,172,183,187,180,185) 
> taikei[,2]<-c(50,60,65,65,70,75,80,85,90, 95) 
> taikei
      [,1] [,2]
 [1,]  165   50
 [2,]  170   60
 [3,]  172   65
 [4,]  175   65
 [5,]  170   70
 [6,]  172   75
 [7,]  183   80
 [8,]  187   85
 [9,]  180   90
[10,]  185   95
> colnames(taikei)<- c("Height","Weight")
> taikei <- data.frame(taikei)
> plot <- ggplot(taikei, aes(x = Height, y = Weight)) +
+   geom_point() + 
+   geom_abline(slope = taikei.lm$coefficients[2], intercept = taikei.lm$coefficients[1])
> plot

再掲ですが，ggplot2の書き方はこのあたりの記事が本当にわかりやすいです．
ggplot2による可視化入門
ggplot2: きれいなグラフを簡単に合理的に - Heavy Watal

この回帰分析の結果をSummary関数で表示します． 以下のように回帰係数やt値，p値，決定係数を取得することができます．

> taikei.F<- data.frame(taikei) 
> taikei.lm<- lm(Weight~Height,data=taikei.F) 
> summary(taikei.lm)

Call:
lm(formula = Weight ~ Height, data = taikei.F)

Residuals:
   Min     1Q Median     3Q    Max 
-7.158 -5.370 -2.764  6.364  9.608 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -222.1647    56.7579  -3.914 0.004454 ** 
Height         1.6809     0.3224   5.213 0.000809 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.158 on 8 degrees of freedom
Multiple R-squared:  0.7726, Adjusted R-squared:  0.7442 
F-statistic: 27.18 on 1 and 8 DF,  p-value: 0.0008091

回帰診断図(Regression Diagnosis Plots)

回帰診断図はplot(lm)で簡単に作成することができます．
以下のように，「残差と予測値の散布図」，「標準化した残差の絶対値の平方根と予測値の散布図」，「残差のnormal Q-Q プロット」，「Cookの距離」の図が表示されます．

par(mfrow=c(2,2))
plot(taikei.lm)　

plot()だけでこのような残差の可視化をしてくれるので非常に便利なのですが，勉強も兼ねてこの図をggplot2で書いてみます．

Residuals VS Fitted

残差のばらつきを確認します．
以下に残差と予測値の散布図を表示させるコードと結果を示します．

> #Residuals VS Fitted
> library(data.table)
> taikei.lm.F <- data.frame(resid(taikei.lm))
> setDT(taikei.lm.F, keep.rownames = "index")
> class(taikei.lm.F$index) #Classが"character"なので変換する必要がある
[1] "character" 
> taikei.lm.F$index <- as.numeric(unlist(taikei.lm.F$index))
> class(taikei.lm.F$index)
[1] "numeric"
> taikei.res <- ggplot(taikei.lm.F, aes(x = index, y = resid.taikei.lm.)) +
+   geom_point() +
+   geom_hline(yintercept=0, linetype="dashed", color = "red")
> taikei.res

上記コードを書く際に参考にした記事をメモします．

• data frameの行名を列に追加する．
  r - Convert row names into first column - Stack Overflow
• ggplot2でy=0の線を引く．
  ggplot2 add straight lines to a plot : horizontal, vertical and regression lines - Easy Guides - Wiki - STHDA
• リスト内の文字列を数値に変換する
  http://www.programmingr.com/r-error-messages/list-object-cannot-be-coerced-to-type-double/

scale_location

続いて，標準化した残差の絶対値の平方根と予測値の散布図を示します． 残差を標準化しており，外れ値を確認することができます．
上の「Residuals VS Fitted」とコードはほぼ同じですが，残差の絶対値の平方根を計算します．

> #scale_location
> pred <-predict(taikei.lm) 
> resid <-rstandard(taikei.lm) 
> taikei.Pred.F <- data.frame(taikei.F,pred,resid)
> taikei.Pred.F$resid <- sqrt(abs(taikei.Pred.F$resid))
> names(taikei.Pred.F)[ which( names(taikei.Pred.F)=="resid" ) ] <- "sqrt(abs(resid))"
> taikei.scale <- ggplot(taikei.Pred.F, aes(x = pred, y = sqrt(abs(resid)))) +
+   geom_point() +
+   labs(title = "scale-location")+
+   scale_y_continuous(limits = c(0, NA))
> taikei.scale

normal Q-Q plot

Q-Qプロットはデータの正規性の確認で使用されます．

続いてnormal Q-Q plotですが，ggplot2でQQplotとQQlineを表示する記事があったので，こちらをそっくり参考にしました．
r - qqnorm and qqline in ggplot2 - Stack Overflow

QQ plotの基礎的なことはこちらの記事がわかりやすかったです．
Normal Q-Q プロットを理解する｜hanaori｜note

> #QQ plot
> ggQQ <- function(LM) # argument: a linear model
+ {
+   y <- quantile(LM$resid[!is.na(LM$resid)], c(0.25, 0.75))
+   x <- qnorm(c(0.25, 0.75))
+   slope <- diff(y)/diff(x)
+   int <- y[1L] - slope * x[1L]
+   p <- ggplot(LM, aes(sample=.resid)) +
+     stat_qq(alpha = 0.5) +
+     geom_abline(slope = slope, intercept = int, color="blue")
+   
+   return(p)
+ }
> QQ.plot <- ggQQ(taikei.lm)
> QQ.plot 

Cook's distance

Cookの距離とは、全サンプルを用いて推定した場合の予測値 と、1サンプルだけ除いて推定した場合の予測値 との差を表しており，Cookの距離が大きいと，その値が回帰モデルへの影響が大きく，外れ値の可能性があります．

> #Cook's distance
> library(data.table)
> taikei.lm.F <- data.frame(cooks.distance(taikei.lm))
> setDT(taikei.lm.F, keep.rownames = "index")
> class(taikei.lm.F$index)
[1] "character"
> taikei.lm.F$index <- as.numeric(unlist(taikei.lm.F$index))
> class(taikei.lm.F$index)
[1] "numeric"
> taikei.cook <- ggplot(taikei.lm.F, aes(x = index, y = cooks.distance.taikei.lm.)) +
+   geom_point() +
+   geom_hline(yintercept=0.5, linetype="dashed", color = "red")
> taikei.cook

まとめ

最後に，これまで示した４つのggplot2の図を４枚まとめて表示させます．
plot()であれば，par(mfrow=c(2,2))でよいのですが，ggplot2の場合はgridExtraというライブラリを使用します．

ggplot2 ライブラリで作った複数のグラフを1枚にまとめるには、gridExtra ライブラリが使える｜Colorless Green Ideas

> library(gridExtra)
>grid.arrange(taikei.res, QQ.plot, taikei.scale, taikei.cook, ncol = 2, nrow = 2)

lm関数のplot()のように，回帰診断図を表示させることができました．