前回は構造がシンプルな単回帰分析を例に、この分析方法の考え方や残差の取り方、および、なぜそのように取る必要があるか等を説明した。今回は、より実用場面で登場する重回帰分析まで拡張し、 分析の進め方や利用時の注意点を含めて紹介する。 いろいろな場面での応用が可能な手法なので、大いに利用していただきたい。
単回帰分析の時も説明したが、後々の分析(残差プロット)時に、データに欠損値が含まれていると計算できないので、分析対象の変量だけを事前に切り出して、欠損値を除外しておく。
なお、今回もより多くのデータを収録した「StudAll20c.xlsx」 をCSV形式に変換したデータ(StudAll20c.csv)を用いて説明を行うが、 「前のデータのままで良い」という場合は、 以下の「StudAll20c」を「StudAll20a」に読み替えて 実行すれば良い。
# データ切り出し
dim(Student20c)
## [1] 494 8
StuHWC<-na.omit(Student20c[,2:4]) # 欠損値を含むデータを削除
dim(StuHWC)
## [1] 157 3
colnames(StuHWC)
## [1] "Height" "Weight" "Chest"
StuHWC[1:10,]
## Height Weight Chest
## 2 146.7 41.0 85
## 4 148.0 43.0 80
## 7 150.0 46.0 86
## 11 151.7 41.5 80
## 12 152.0 35.0 77
## 17 153.0 46.5 87
## 19 153.0 55.0 78
## 24 154.4 44.0 75
## 26 155.0 48.0 83
## 31 156.0 42.0 85
元々の身長と体重のデータは494サンプルだったが、 欠損値を含むデータが337サンプルあり、残ったのは157サンプルである。 変量が多くなると欠損値を含むサンプルが多くなる傾向にある。 以後は、Student20cの代わりにStudHWCを使って分析を行う。
2変量のときには有用性が理解してもらえないと思ったので紹介しなかったが、 3変量以上の場合に変量間の関係を概観するのに便利な「散布図行列」を紹介しておく。 2変量ごとを組み合わせて散布図を行列形式で配置したものである。 変量感でどの程度関係性がありそうかを視覚的に把握することができて便利である。
# 散布図行列
plot(StuHWC)
式の展開、解法。
直線の方程式:
測定値と予測値のズレ:
ズレの2乗の和を最小に:
# 重回帰分析
RresultHWC<-lm(Weight~Height+Chest, data=StuHWC) # 右辺に説明する変量をプラスマークで追加していく
# 分析結果要約
summary(RresultHWC)
##
## Call:
## lm(formula = Weight ~ Height + Chest, data = StuHWC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.894 -3.868 -1.509 2.849 33.974
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -98.84919 12.07580 -8.186 9.56e-14 ***
## Height 0.77130 0.07213 10.694 < 2e-16 ***
## Chest 0.33565 0.06579 5.101 9.79e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.223 on 154 degrees of freedom
## Multiple R-squared: 0.5384, Adjusted R-squared: 0.5324
## F-statistic: 89.82 on 2 and 154 DF, p-value: < 2.2e-16
# 残差プロット
plot(StuHWC$Height, RresultHWC$residuals) # 対身長
abline(v=seq(140,190,5), lty=3) # 点線を追記
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # ゼロを実線で追記
plot(StuHWC$Chest, RresultHWC$residuals) # 対胸囲
abline(v=seq(30,120,5), lty=3) # 点線を追記
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # ゼロを実線で追記
plot(StuHWC$Weight, RresultHWC$residuals) # 対体重
abline(v=seq(30,100,5), lty=3) # 点線を追記
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # ゼロを実線で追記
残差の正規性を確かめてみよう。
# 残差のヒストグラム
hist(RresultHWC$residuals, right=F)
abline(h=seq(0,150,25), lty=3) # 点線を追記
# 残差の箱ひげ図
boxplot(RresultHWC$residuals, horizontal=T)
abline(v=seq(-10,30,10), lty=3) # 点線を追記
# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHWC$residuals)
qqline(RresultHWC$residuals, lty=2)
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # 実線を追記
[検討事項] 残差が20を超える辺りの4サンプルほどは(少なくとも)分布を乱しているように見える。===> 外れ値か? 吟味が必要。
残差分析の結果、残差が20を超えている4サンプルについて吟味したところ、 一般的な学生の体型と異なっていると判断し、 これらを除外して分析をしてみよう。
元のExcelファイルに戻って除外する手もあるが、Rの中で実現するには、 不等号や否定(!, Not)等をうまく組み合わせて、 残差が20を超えるサンプルを除外する。 それぞれの式の意味している事項が理解できるであろうか? もし理解できないようなら、式の中を分割して出力してみれば理解の促進になるであろう。
# 残差の大きいサンプルを抽出
StuResidOver20<-RresultHWC$residuals>20
RresultHWC$residuals[StuResidOver20] # 20超えの残差一覧
## 231 340 341 428
## 25.21708 33.97434 25.20585 23.96511
StuHWC[StuResidOver20,] # 20超えのサンプル一覧
## Height Weight Chest
## 231 169.3 88.5 94
## 340 173.0 84.0 46
## 341 173.0 90.0 90
## 428 178.0 100.0 112
20以下のサンプルの抽出。157サンプルから4サンプルを除外したので153サンプルになっている。以下では変数名をStuHWC2とした。
StuHWC2<-StuHWC[!StuResidOver20,]
StuHWC2[1:10,]
## Height Weight Chest
## 2 146.7 41.0 85
## 4 148.0 43.0 80
## 7 150.0 46.0 86
## 11 151.7 41.5 80
## 12 152.0 35.0 77
## 17 153.0 46.5 87
## 19 153.0 55.0 78
## 24 154.4 44.0 75
## 26 155.0 48.0 83
## 31 156.0 42.0 85
dim(StuHWC2)
## [1] 153 3
では、抽出したデータ(153サンプル)を用いて改めて重回帰分析を行ってみよう。
# 重回帰分析
RresultHWC2<-lm(Weight~Height+Chest, data=StuHWC2)
# 分析結果要約
summary(RresultHWC2)
##
## Call:
## lm(formula = Weight ~ Height + Chest, data = StuHWC2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.262 -3.502 -1.004 3.010 22.114
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -94.22920 9.66924 -9.745 < 2e-16 ***
## Height 0.71096 0.05774 12.314 < 2e-16 ***
## Chest 0.39223 0.05795 6.769 2.74e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.719 on 150 degrees of freedom
## Multiple R-squared: 0.6365, Adjusted R-squared: 0.6317
## F-statistic: 131.3 on 2 and 150 DF, p-value: < 2.2e-16
# 残差プロット
plot(StuHWC2$Height, RresultHWC2$residuals) # 対身長
abline(v=seq(140,190,5), lty=3) # 点線を追記
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # ゼロを実線で追記
plot(StuHWC2$Chest, RresultHWC2$residuals) # 対胸囲
abline(v=seq(30,120,5), lty=3) # 点線を追記
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # ゼロを実線で追記
plot(StuHWC2$Weight, RresultHWC2$residuals) # 対体重
abline(v=seq(30,100,5), lty=3) # 点線を追記
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # ゼロを実線で追記
残差の正規性を確かめてみよう。
# 残差のヒストグラム
hist(RresultHWC2$residuals, right=F)
abline(h=seq(0,150,25), lty=3) # 点線を追記
# 残差の箱ひげ図
boxplot(RresultHWC2$residuals, horizontal=T)
abline(v=seq(-10,30,10), lty=3) # 点線を追記
# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHWC2$residuals)
qqline(RresultHWC2$residuals, lty=2)
abline(h=seq(-10,30,10), lty=3) # 点線を追記
abline(h=0, lty=1) # 実線を追記
# 分析結果要約 = 再掲
summary(RresultHWC2)
##
## Call:
## lm(formula = Weight ~ Height + Chest, data = StuHWC2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.262 -3.502 -1.004 3.010 22.114
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -94.22920 9.66924 -9.745 < 2e-16 ***
## Height 0.71096 0.05774 12.314 < 2e-16 ***
## Chest 0.39223 0.05795 6.769 2.74e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.719 on 150 degrees of freedom
## Multiple R-squared: 0.6365, Adjusted R-squared: 0.6317
## F-statistic: 131.3 on 2 and 150 DF, p-value: < 2.2e-16
[参考: 出力の読み方の詳細] 統計検定のための、Rの出力結果からわかること(回帰分析編)
このページで取り扱ったプログラムだけを抜き出して以下に列挙しておく。
# ディレクトリの移動。必須ではない。個々人の設定に応じて。 setwd("D:/home_sub3/R_Dir") # ホームディレクトリに移動(Set Working Directory) ## setwd("C:/home/R_Dir") # ホームディレクトリに移動(Set Working Directory) getwd() # 現在のディレクトリ位置を表示 list.files() # ファイル・ディレクトリ一覧を表示 setwd("KougiDS") # ディレクトリを移動 list.files() # ファイル・ディレクトリ一覧を表示 # データの読み込み Ver.c=20年の2つのクラスのデータ Student20c<-read.csv("StudAll20c.csv", skip=6, header=TRUE, na.strings="NA") dim(Student20c) colnames(Student20c) Student20c[1:5,] # データ切り出し dim(Student20c) StuHWC<-na.omit(Student20c[,2:4]) # 欠損値を含むデータを削除 dim(StuHWC) colnames(StuHWC) StuHWC[1:10,] # 散布図行列 plot(StuHWC) # 重回帰分析 RresultHWC<-lm(Weight~Height+Chest, data=StuHWC) # 右辺に説明する変量をプラスマークで追加していく # 分析結果要約 summary(RresultHWC) # 残差プロット plot(StuHWC$Height, RresultHWC$residuals) # 対身長 abline(v=seq(140,190,5), lty=3) # 点線を追記 abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # ゼロを実線で追記 plot(StuHWC$Chest, RresultHWC$residuals) # 対胸囲 abline(v=seq(30,120,5), lty=3) # 点線を追記 abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # ゼロを実線で追記 plot(StuHWC$Weight, RresultHWC$residuals) # 対体重 abline(v=seq(30,100,5), lty=3) # 点線を追記 abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # ゼロを実線で追記 # 残差のヒストグラム hist(RresultHWC$residuals, right=F) abline(h=seq(0,150,25), lty=3) # 点線を追記 # 残差の箱ひげ図 boxplot(RresultHWC$residuals, horizontal=T) abline(v=seq(-10,30,10), lty=3) # 点線を追記 # 残差のQQプロット: 正規確率プロット qqnorm(RresultHWC$residuals) qqline(RresultHWC$residuals, lty=2) abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # 実線を追記 # 残差の大きいサンプルを抽出 StuResidOver20<-RresultHWC$residuals>20 RresultHWC$residuals[StuResidOver20] # 20超えの残差一覧 StuHWC[StuResidOver20,] # 20超えのサンプル一覧 StuHWC2<-StuHWC[!StuResidOver20,] StuHWC2[1:10,] dim(StuHWC2) # 重回帰分析 RresultHWC2<-lm(Weight~Height+Chest, data=StuHWC2) # 分析結果要約 summary(RresultHWC2) # 残差プロット plot(StuHWC2$Height, RresultHWC2$residuals) # 対身長 abline(v=seq(140,190,5), lty=3) # 点線を追記 abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # ゼロを実線で追記 plot(StuHWC2$Chest, RresultHWC2$residuals) # 対胸囲 abline(v=seq(30,120,5), lty=3) # 点線を追記 abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # ゼロを実線で追記 plot(StuHWC2$Weight, RresultHWC2$residuals) # 対体重 abline(v=seq(30,100,5), lty=3) # 点線を追記 abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # ゼロを実線で追記 # 残差のヒストグラム hist(RresultHWC2$residuals, right=F) abline(h=seq(0,150,25), lty=3) # 点線を追記 # 残差の箱ひげ図 boxplot(RresultHWC2$residuals, horizontal=T) abline(v=seq(-10,30,10), lty=3) # 点線を追記 # 残差のQQプロット: 正規確率プロット qqnorm(RresultHWC2$residuals) qqline(RresultHWC2$residuals, lty=2) abline(h=seq(-10,30,10), lty=3) # 点線を追記 abline(h=0, lty=1) # 実線を追記 # 分析結果要約 = 再掲 summary(RresultHWC2)