8章_回帰分析
1 ファイルの読み込み
dat <- read.csv("data/data_ch8-1.csv", header = TRUE, row.names = 1)#row.namesは行名の指定。row.names=1だと1列目が列名になる
head(dat)| mid | end | quiz | |
|---|---|---|---|
| S001 | 77 | NA | 36 |
| S002 | 54 | 47 | 28 |
| S003 | 73 | 65 | 31 |
| S004 | 63 | 60 | 26 |
| S005 | 44 | 62 | 24 |
| S006 | 42 | 58 | 28 |
dim(dat)## [1] 100 31.1 欠損値（NA）の有無を確認
complete.cases(dat) ## [1] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [37] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [73] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [97] TRUE TRUE TRUE TRUEsubset(dat, complete.cases(dat) == FALSE)| mid | end | quiz | |
|---|---|---|---|
| S001 | 77 | NA | 36 |
1.2 NAを含む行を削除
dat.2 <- na.omit(dat[, 1 : 2])
dim(dat.2)## [1] 99 2head(dat.2)| mid | end | |
|---|---|---|
| S002 | 54 | 47 |
| S003 | 73 | 65 |
| S004 | 63 | 60 |
| S005 | 44 | 62 |
| S006 | 42 | 58 |
| S007 | 57 | 45 |
2 Mahalanobisの距離
重心からの距離.
マハラノビス距離による外れ値の検索[R]
apply(dat.2, 2, mean)#列平均の計算（２は列の平均を意味する、１は行）## mid end
## 56.70707 54.21212d <- mahalanobis(dat.2, apply(dat.2, 2, mean), cov(dat.2))#平均からの距離の計算
head(d)## S002 S003 S004 S005 S006 S007
## 0.2224415 1.1729144 0.2066332 1.8564335 1.6609444 0.49547452.1 データの行数と列数の計算
n <- nrow(dat.2)
v <- ncol(dat.2)2.2 外れ値の検出
総務省統計局 和田かず美 様 「多変量外れ値の検出」によると、
\(F=\dfrac{(n-p)n}{(n^2-1)p}D^2(x)\sim
F(p,n-p)\).
マハラノビスの距離が大きくなるとF値が大きくなることを用いて外れ値を検出する。
out <- n * (n - v) / ((n ^ 2 - 1) * v) * d > qf(0.9, n, v)#qfは分布関数
head(out)## S002 S003 S004 S005 S006 S007
## FALSE FALSE FALSE FALSE FALSE FALSE2.3 散布図に外れ値を表示
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(dat.2, pch = ifelse(out, 16, 21), xlim = c(0, 100), ylim = c(0, 100), xlab = "中間試験", ylab = "期末試験")
#プロット点の形を指定するパラメータ pch 16黒円 21白丸
#text(dat.2[out, ] - 3, labels = paste(which(out), ":", rownames(dat.2)[out]))
# 外れ値の行番号とラベルを表示
dat.2[out, ] | mid | end | |
|---|---|---|
| S035 | 88 | 5 |
| S069 | 92 | 11 |
2.4 外れ値の除外
dat.3 <- dat.2[-which(out == TRUE), ]3 単回帰モデル
3.1 結果変数 ~ 予測変数
model.1 <- lm(end ~ mid, data = dat.3)
summary(model.1)##
## Call:
## lm(formula = end ~ mid, data = dat.3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.2474 -5.4599 -0.7961 5.0888 26.6908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.5133 3.7623 3.858 0.000208 ***
## mid 0.7257 0.0650 11.164 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.314 on 95 degrees of freedom
## Multiple R-squared: 0.5674, Adjusted R-squared: 0.5629
## F-statistic: 124.6 on 1 and 95 DF, p-value: < 2.2e-16#回帰診断
par(mfrow=c(2,2)) #par(mfrow=c(2,2),oma=c(2,2,2,1),mar=c(5,4,3,2))
plot(model.1) # plot(taikei.lm,cex=1.5,pch=21,bg="blue",col="blue") 
3.2 回帰直線の可視化
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(dat.3, xlim = c(0, 100), ylim = c(0, 100), xlab = "中間試験", ylab = "期末試験")
abline(model.1)
3.3 傾きと切片の95%信頼区間
(confint(model.1, level = 0.95))## 2.5 % 97.5 %
## (Intercept) 7.0442767 21.982321
## mid 0.5966093 0.8547013.4 95%信頼区間の可視化
# 新しいデータ（0点から100点まで1点刻み）を用意
new <- data.frame(mid = seq(0, 100, 1))
confidence <- predict(model.1, newdata = new, interval ='confidence', level = 0.95)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(dat.3, xlim = c(0, 100), ylim = c(0, 100), xlab = "中間試験", ylab = "期末試験")
abline(model.1)
lines(new$mid, confidence[, 2], lty = 3)
lines(new$mid, confidence[, 3], lty = 3)
str(confidence[, 2])## Named num [1:101] 7.04 7.89 8.75 9.6 10.45 ...
## - attr(*, "names")= chr [1:101] "1" "2" "3" "4" ...str(new$mid)## num [1:101] 0 1 2 3 4 5 6 7 8 9 ...3.5 点数の予測
pred <- predict(model.1, newdata = new, interval = 'prediction', level = 0.95)
head(pred)## fit lwr upr
## 1 14.51330 -5.428320 34.45492
## 2 15.23895 -4.656247 35.13415
## 3 15.96461 -3.884904 35.81412
## 4 16.69026 -3.114297 36.49483
## 5 17.41592 -2.344430 37.17627
## 6 18.14157 -1.575309 37.85846head(sort(dat.3$ mid))## [1] 16 20 26 31 31 33pred[6, ]## fit lwr upr
## 18.141574 -1.575309 37.8584574 予測区間の可視化
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(dat.3, xlim = c(0, 100), ylim = c(0, 100), xlab = "中間試験", ylab = "期末試験")
abline(model.1)
lines(new$mid, pred[, 2], lty = 2)
lines(new$mid, pred[, 3], lty = 2)
5 残差の等分散性と正規性のプロット
- 予測値に対する残差のプロット（誤差の独立性、誤差分散の均一性）
- 残差の正規 Q-Q プロット （誤差の正規性）
par(mfcol = c(1, 2))
plot(model.1, which = c(1, 2))
6 重回帰モデル
6.1 ファイル
library("pequod")## 要求されたパッケージ car をロード中です## 要求されたパッケージ carData をロード中です##
## 次のパッケージを付け加えます: 'car'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## recode## 以下のオブジェクトは 'package:purrr' からマスクされています:
##
## somehead(dat)| mid | end | quiz | |
|---|---|---|---|
| S001 | 77 | NA | 36 |
| S002 | 54 | 47 | 28 |
| S003 | 73 | 65 | 31 |
| S004 | 63 | 60 | 26 |
| S005 | 44 | 62 | 24 |
| S006 | 42 | 58 | 28 |
dat.4 <- na.omit(dat)6.2 外れ値の除外
d.2 <- mahalanobis(dat.4, apply(dat.4, 2, mean), cov(dat.4))
n.2 <- nrow(dat.4)
v.2 <- ncol(dat.4)
out.2 <- n.2 * (n.2 - v.2) / ((n.2 ^ 2 - 1) * v.2) * d.2 > qf(0.9, n.2, v.2)
dat.4[out.2, ]| mid | end | quiz | |
|---|---|---|---|
| S035 | 88 | 5 | 45 |
| S069 | 92 | 11 | 35 |
dat.5 <- dat.4[-which(out.2 == TRUE), ]
head(dat.5)| mid | end | quiz | |
|---|---|---|---|
| S002 | 54 | 47 | 28 |
| S003 | 73 | 65 | 31 |
| S004 | 63 | 60 | 26 |
| S005 | 44 | 62 | 24 |
| S006 | 42 | 58 | 28 |
| S007 | 57 | 45 | 22 |
write.csv(dat.5, file = "data/dat.5.csv", row.names = TRUE)6.3 重回帰モデル
Rで重回帰分析で交互作用を検討する方法(pequodパッケージ)
model.2 <- lmres(end ~ mid + quiz, data = dat.5)
summary(model.2)## Formula:
## end ~ mid + quiz
## <environment: 0x7f88a7b1ced8>
##
## Models
## R R^2 Adj. R^2 F df1 df2 p.value
## Model 0.802 0.642 0.635 84.466 2.000 94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -19.703 -4.910 0.029 0.000 4.757 18.907
##
## Coefficients
## Estimate StdErr t.value beta p.value
## (Intercept) 14.11264 3.43967 4.10291 9e-05 ***
## mid 0.47964 0.08122 5.90549 0.4979 <2e-16 ***
## quiz 0.54056 0.12169 4.44216 0.3745 2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Collinearity
## VIF Tolerance
## mid 1.869 0.535
## quiz 1.869 0.535betaは標準化回帰係数で、注目する変数以外の説明変数の値を0に固定した上で，注目する説明変数が1標準偏差変化した時に従属変数が何標準偏差変化するかを表します。 単回帰分析の場合、標準化回帰係数は相関係数と同値である。
Collinearityは多重共線性の診断結果
VIFが10以上またはTolerance 0.1以下の場合、多重共線性がある。
この場合、変数を除外することになる。
交互作用項(x1x2)はx1とx2で作られた変数であるため、x1あ るいはx2 と相関が非常に高くなってしまう(→多重共線性の可能性)
6.4 回帰平面の作成
library(scatterplot3d)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot3d <- scatterplot3d(dat.4$mid, dat.4$quiz, dat.4$end, type = "n",
xlim = c(0, 100), ylim = c(0, 50), zlim = c(0, 100), angle = 25,
xlab = "中間試験", ylab = "期末試験", zlab = "小テスト")
plot3d$plane3d(lm(end ~ mid + quiz, data = dat.4),
lty.box = "solid", col = "grey", draw_polygon = TRUE, draw_lines = TRUE,
polygon_args = list(col = rgb(0.97, 0.97, 0.97, 0.97)))
orig <- plot3d$xyz.convert(dat.4$mid, dat.4$quiz, dat.4$end)
plane <- plot3d$xyz.convert(dat.4$mid, dat.4$quiz, fitted(lm(end ~ mid + quiz, data = dat.4)))
i.negpos <- 1 + (resid(lm(end ~ mid + quiz, data = dat.4)) > 0)
segments(orig$x, orig$y, plane$x, plane$y,
lty = c(2, 1)[i.negpos])
wh1 <- resid(lm(end ~ mid + quiz, data = dat.4)) > 0
wh2 <- resid(lm(end ~ mid + quiz, data = dat.4)) < 0
segments(orig$x[wh1], orig$y[wh1], plane$x[wh1], plane$y[wh1], lty = 1)
plot3d$points3d(dat.4$mid[wh1], dat.4$quiz[wh1], dat.4$end[wh1], pch = 19)
plot3d$points3d(dat.4$mid[wh2], dat.4$quiz[wh2], dat.4$end[wh2], pch = 21)
6.5 交互作用の検討
model.3 <- lmres(end ~ mid + quiz + mid : quiz, centered = c("mid", "quiz"), data = dat.5)
summary(model.3)## Formula:
## end ~ mid + quiz + mid.XX.quiz
## <environment: 0x7f88c16cb568>
##
## Models
## R R^2 Adj. R^2 F df1 df2 p.value
## Model 0.813 0.660 0.649 60.289 3.000 93 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -20.1458 -5.6955 -0.1822 0.0000 5.0320 20.9739
##
## Coefficients
## Estimate StdErr t.value beta p.value
## (Intercept) 56.26278 0.98115 57.34392 < 2e-16 ***
## mid 0.43949 0.08162 5.38476 0.4562 < 2e-16 ***
## quiz 0.57447 0.12021 4.77878 0.3980 1e-05 ***
## mid.XX.quiz -0.01140 0.00514 -2.21574 -0.1374 0.02915 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Collinearity
## VIF Tolerance
## mid 1.9659 0.5087
## quiz 1.8998 0.5264
## mid.XX.quiz 1.0532 0.94956.6 単純傾斜の検定
eff.1 <- simpleSlope(model.3, pred = "mid", mod1 = "quiz")
summary(eff.1)##
## ** Estimated points of end **
##
## Low mid (-1 SD) High mid (+1 SD)
## Low quiz (-1 SD) 42.602 58.709
## High quiz (+1 SD) 57.070 66.670
##
##
##
## ** Simple Slopes analysis ( df= 93 ) **
##
## simple slope standard error t-value p.value
## Low quiz (-1 SD) 0.5507 0.0858 6.42 <2e-16 ***
## High quiz (+1 SD) 0.3283 0.1049 3.13 0.0023 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
##
## ** Bauer & Curran 95% CI **
##
## lower CI upper CI
## quiz 16.311 363.6eff.2 <- simpleSlope(model.3, pred = "quiz", mod1 = "mid")
summary(eff.2)##
## ** Estimated points of end **
##
## Low quiz (-1 SD) High quiz (+1 SD)
## Low mid (-1 SD) 42.602 57.07
## High mid (+1 SD) 58.709 66.67
##
##
##
## ** Simple Slopes analysis ( df= 93 ) **
##
## simple slope standard error t-value p.value
## Low mid (-1 SD) 0.741 0.150 4.95 <2e-16 ***
## High mid (+1 SD) 0.408 0.133 3.06 0.0029 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
##
## ** Bauer & Curran 95% CI **
##
## lower CI upper CI
## mid 23.195 418.866.7 単純傾斜のプロット
PlotSlope(eff.1)
PlotSlope(eff.2)
6.8 検定によるモデル選択
anova(model.1, model.2$Stepfin, model.3$Stepfin)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 95 | 8240.790 | NA | NA | NA | NA |
| 94 | 6811.001 | 1 | 1429.7893 | 20.553504 | 0.0000172 |
| 93 | 6469.476 | 1 | 341.5248 | 4.909487 | 0.0291498 |
6.9 情報量基準によるモデル選択
AIC(model.1, model.2$Stepfin, model.3$Stepfin)| df | AIC | |
|---|---|---|
| model.1 | 3 | 712.1617 |
| model.2\(Stepfin | 4| 695.6777| |model.3\)Stepfin | 5 | 692.6876 |
6.10 3D
#3D表示
options(rgl.printRglwidget = TRUE)
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
df <- data.frame(x=dat.5$mid,y=dat.5$quiz,z=dat.5$end)
head(df)| x | y | z |
|---|---|---|
| 54 | 28 | 47 |
| 73 | 31 | 65 |
| 63 | 26 | 60 |
| 44 | 24 | 62 |
| 42 | 28 | 58 |
| 57 | 22 | 45 |
attach(df)
res <- aov(z ~ x * y)
coefs <- coef(res)
a <- coefs["x"];a## x
## 0.7385204b <- coefs["y"];b## y
## 1.212942c <- coefs["x:y"];c## x:y
## -0.01139717X <- seq(min(df$x),max(df$x),len=100)#定義域方眼
Y <- seq(min(df$y),max(df$y),len=100)#定義域方眼
plot.df <- expand.grid(x=X,y=Y)#定義域方眼
#定義域方眼の座標を(A,B,C)とする
A <- plot.df$x
B <- plot.df$y
C <- a*A+b*B+c*A*B#z座標の列を追加
plot3d(A,B,C,col="lightblue",
alpha = 0.2,xlab="x",ylab="y",zlab="z")
plot3d(x, y, z, type = "s", col = "blue", size = 1,add=T)curve(a*x+b*mean(y)+c*x*mean(y),xlim=c(10,max(df$x)),ylim=c(40,70),col="blue")
curve(a*x+b*(mean(y)+10)+(mean(y)+10)*c*x,add=T,col="pink")
curve(a*x+b*(mean(y)-10)+(mean(y)-10)*c*x,add=T,col="green")
7 マハラノビスの距離とF分布
\(F=\dfrac{(n-p)n}{(n^2-1)p}D^2(x)\sim
F(p,n-p)\).
このF値がF分布に従うことを具体例で確かめておく。
library(MASS)##
## 次のパッケージを付け加えます: 'MASS'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## selectmu <- c(0,0,0)
Sigma <- matrix(c(1, 0,0.6,0,1,0.5,0.6,0.5, 1), 3, 3)
data <- mvrnorm(100, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
attach(d)## 以下のオブジェクトは .GlobalEnv によってマスクされています:
##
## X, Ycov(data)## [,1] [,2] [,3]
## [1,] 1.0003765 0.0223299 0.6714040
## [2,] 0.0223299 1.0282425 0.6023764
## [3,] 0.6714040 0.6023764 1.1328329colMeans(d)## X Y Z
## 0.06321200 0.02006497 -0.04697340ff <- mahalanobis(d, colMeans(d), cov(data))*97*100/(100^2-1)/3
hist(ff,prob=T)
curve(df(x, 3, 97),add=T)