単回帰分析

単回帰分析とは

テキストでは  

・データの分布の中心に直線を求める手法

補足  

・直線関係のありそうな二組のデータの関係を数式として明らかにすることで関係性を見出す分析手法
・yi=a+bxi

単回帰分析の手順  

Step1. 散布図を書く

Step2. 相関係数を計算する

Step3. 単回帰直線にあてはめ,決定係数を求める

ではテキストの分析をやっていきましょう

データの場所を指定

setwd(“ファイルの場所”)

ファイルの読み込み

dat<-read.csv("data_ch8-1.csv", header = TRUE, row.names = 1)

欠損値(NA)の有無を確認

subset(dat,complete.cases(dat)==FALSE)
##      mid end quiz
## S001  77  NA   36
dat.2 <- na.omit(dat[, 1 : 2])
#NAを含む行を削除したデータの確認
head(dat.2)
##      mid end
## S002  54  47
## S003  73  65
## S004  63  60
## S005  44  62
## S006  42  58
## S007  57  45

Step1. 散布図を書く

関係を知りたいデータをx軸,y軸に従ってプロットする (family=“font”の部分はwindowsでは書く必要なし)

plot(dat.2, xlim = c(0, 100), ylim = c(0, 100), xlab ="中間試験", ylab = "期末試験",family="HiraKakuProN-W3")

散布図から読み取ること

1.データの分布が直線に近いかどうか

データの分布が直線に近い→相関がある

データの分布が直線に近くない→相関は低い

2.直線の傾きの方向性

右肩上がりなら正の相関   右肩下がりなら負の相関

テキストに戻ります。プロットの右すみにあるデータは?

外れ値かどうかを判断する

mahalanobisの距離

d <- mahalanobis(dat.2, apply(dat.2, 2, mean), cov(dat.2))

データの行数と列数の計算

n<-nrow(dat.2)
v<-ncol(dat.2)

外れ値の検出

out <- n*(n-v) / ((n ^ 2-1)* v)*d > qf(0.9, n,v)

散布図に外れ値を表示

plot(dat.2, pch = ifelse(out, 16, 21), xlim = c(0, 100), ylim = c(0, 100), 
     xlab = "中間試験", ylab = "期末試験", family = "HiraKakuProN-W3")

# 外れ値の行番号とラベルを表示
text(dat.2[out, ]-3,labels=paste(which(out), ":",rownames(dat.2)[out]))

外れ値の除外

dat.3<-dat.2[-which(out==TRUE),]

(テキストにはなかったが相関係数を求める)

(理由)2グループの間に相関があることを保証する(単回帰分析の前提となる)
式はcor(x,y)

mid<-dat.3$mid
end<-dat.3$end
cor(mid,end)
## [1] 0.7532885

相関の考え方

絶対値として0から1の間の数値として1に近いほど相関が高く,低いほど相関が低い (標本サイズが多ければノイズが入るため多少数値が小さくても相関があると考えることもできる)

(参考)相関係数を判断する目安(Rによるデータ分析のレシピより)
-1〜-0.3:強い負の相関
-0.7〜-0.3:負の相関あり
-0.3〜0.3:相関なし
0.3〜0.7:正の相関あり
0.7〜1:強い正の相関

単回帰分析で出るが相関係数から決定係数を求めることができる 決定係数=相関係数の2乗

#rに代入
r<-cor(mid,end)
r^2
## [1] 0.5674435

単回帰モデル(結果変数~予測変数)

中間試験の結果から期末試験の成績を予測する単回帰分析を行う

確認したいこと

結果変数/目的変数
予測変数/説明変数

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

データの見方

14.5133は切片,0.7257は傾き
Multiple R-squaredは調整済みの決定係数R二乗
R二乗が1に近いほど,当てはまりがよい。

最小二乗法・・・

yの残差の二乗の総和が最小となるようにaとbを決定する。

単回帰分析におけるp値

傾きが0である(つまりxがいくらであってもyに影響しない)という帰無仮説に対するもの

(確認したいこと)

偏回帰係数は、重回帰分析の用語ではないか? 重回帰分析は説明変数が複数あるため,傾きが複数になり,それらの傾きのことを偏回帰係数と呼ぶ。

→テキストでは,偏回帰係数とかかれているが,ここは回帰係数が正解

回帰直線の可視化

plot(dat.3, xlim = c(0, 100), ylim = c(0, 100), xlab = "中間試験", ylab = "期末試験",
     family = "HiraKakuProN-W3")
abline(model.1)

#  傾きと切片の95%の信頼区間
confint(model.1, level = 0.95)
##                 2.5 %    97.5 %
## (Intercept) 7.0442767 21.982321
## mid         0.5966093  0.854701
# 新しいデータ(0点から100点まで1点刻み)を用意
new<-data.frame(mid=seq(0,100,1))

# 95%信頼区間の可視化
confidence <- predict(model.1, newdata = new, interval = 'confidence', level = 0.95)
lines(new$mid,confidence[,2], lty=3)
lines(new$mid,confidence[,3], lty=3)

信頼区間とは,同様の標本と回帰分析を100回行ったら,100回のうち95回ほどこの範囲に母集団のモデルが含むということ意味する。

(確認ポイント) new<-data.frame(mid=seq(0,100,1))はどういった意味か

中間試験が60点の場合の期末試験の点数の予測

pred <- predict(model.1, newdata = new, interval = 'prediction', level = 0.95)
pred[60, ]
##      fit      lwr      upr 
## 57.32695 38.73786 75.91604

予測区間の可視化

得られるデータの95%が含まれる範囲を示す。これにより,新しく得られるデータを予測することができる。

信頼区間は母集団についての範囲
予測区間は得られたデータについての範囲,即ちデータ自体の特徴を示す

plot(dat.3, xlim = c(0, 100), ylim = c(0, 100), xlab = "中間試験", ylab = "期末試験",
     family = "HiraKakuProN-W3")
abline(model.1)
lines(new$mid, pred[, 2], lty = 2)
lines(new$mid, pred[, 3], lty = 2)

残差(誤差)の等分散性と正規性のプロット

残差とは回帰分析において各データと予測値との差のこと。つまり直線と実データとの距離

残差の等分散性はy=0を中心にデータが均等に散らばっているかを示す。 残差の正規性は各データがQ-Qプロットから離れすぎていないかを検討する。

par(mfcol = c(1, 2))
plot(model.1, which = c(1, 2))