重回帰分析
被験者諸属性と試験結果(架空データ)を用いる。
ある集団から無作為に選んだ326人に学力テストを行い、その点数を記録し、さらに受験者のさまざまな属性と合わせて結果を分析する。得点と身長・年齢の関係について回帰分析を行う。
説明変数が複数になる。
回帰直線
回帰分析における近似直線を回帰直線と呼ぶ。回帰直線は傾きと切片によりあたえられる。独立変数が複数取れるので直線の式は一般に
y = a + b1x1 + b2x2 + … + bnxn
であたえられる。
data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/exam-result.csv", fileEncoding = "utf-8")
summary(data)
#サンプルサイズ
N <- length(data$no)
x.colnum <- c(3:4) #説明変数は3列目と4列目
y.colnum <- 2 #目的変数は2列目
#説明変数が先、目的変数がその後に来るようにデータを作り替える(Rでの作業の都合上)
data.tmp <- data.frame(c(data[x.colnum],data[y.colnum]))
x.colnum <- c(1:length(x.colnum))
y.colnum <- length(x.colnum) + 1
#説明変数の偏差積和行列
Sxx <- combn(x.colnum,1,function(x2)combn(x.colnum,1,function(x1)sum ((data.tmp[[x1]] - mean(data.tmp[[x1]]))*(data.tmp[[x2]] - mean(data.tmp[[x2]]))) ))
#説明変数と目的変数の偏差積和行列
Sxy <- combn(x.colnum,1,function(x)sum ((data.tmp[[x]] - mean(data.tmp[[x]]))*(data.tmp[[y.colnum]] - mean(data.tmp[[y.colnum]]))) )
#偏回帰係数
b <- Sxy %*% solve(Sxx)
#切片
a <- mean(data.tmp[[y.colnum]]) - b %*% combn(x.colnum,1,function(x) mean(data.tmp[[x]]))
#回帰式
expected.y <- apply(combn(x.colnum,1,function(x)data.tmp[[x]]*b[,x]),1,function(x)sum(x) + a)
- 傾きbi…説明変数の偏差積和逆行列*説明変数と目的変数の偏差積和行列
- 切片a…目的変数の平均値-傾きb1*説明変数1の平均-傾きb2*説明変数2の平均…
※説明変数が複数になると行列計算が必要になる。おまけにRで一般的な形で式を立てようとすると組み合わせ(Combine)が必要になり、式自体はかなり複雑になってしまう。ただ実質的に行っていることは単回帰分析の時とさほど変わらない。
> b
[,1] [,2]
[1,] 0.003497438 6.641647
> a
[,1]
[1,] -5.548134
- 身長の係数が正の値なので身長が高ければ得点も上がる(身長1cmに付き得点は0.0035点上がる)。
- 年齢の係数が正の値なので年齢が高ければ得点も上がる(年齢1歳に付き得点は6.64点上がる)。
※身長1cmと年齢1歳の得点にあたえる影響について単純な大小比較は出来ない。
分散分析
#分散分析表 #変動 r.SS <- sum((expected.y - mean(data.tmp[[y.colnum]]))^2) e.SS <- sum((data.tmp[y.colnum] - expected.y)^2) t.SS <- sum((data.tmp[y.colnum] - mean(data.tmp[[y.colnum]]))^2) #自由度 r.df <- length(x.colnum) e.df <- (length(data.tmp[[y.colnum]]) - 1 - length(x.colnum)) t.df <- (length(data.tmp[[y.colnum]]) - 1) #分散 r.MS <- r.SS / r.df e.MS <- e.SS / e.df t.MS <- t.SS / t.df #F検定 f.value <- r.MS / e.MS p.value <- 1 - pf(f.value,r.df,e.df)
- 回帰自由度は説明変数の数
> p.value [1] 0
- 分散分析より説明変数の変動により目的変数も変動すると言える。つまりこの回帰モデルは有意である。
回帰係数の有意性検定
#回帰係数の標準誤差 b.se <- combn(x.colnum,1,function(x)sqrt(e.MS*abs(solve(Sxx)))[x,x]) #回帰係数の有意性検定 t.value <- b[1,]/b.se b.p.value <- (1 - pt(abs(t.value),e.df))*2 #回帰係数のF値 b.f.value <- t.value^2
- 回帰係数ごとに有意性が出る。
- 回帰係数のF値はt値を二乗したもの。F値が2より大きいものは変数選択の際に残す。
> b.p.value [1] 0.9154989 0.0000000 > b.f.value [1] 0.01127596 952.18491566
- 身長の回帰係数のp値…0.9154989←有意ではない
- 年齢の回帰係数のp値…0.0000000←1%水準で有意である
- 身長の回帰係数のF値…0.01←2より小さいので変数選択の際は外す
- 年齢の回帰係数のF値…952.18←2より大きいので変数選択の際は残す
標準化回帰係数
#標準化回帰係数 x.sd <- combn(x.colnum,1,function(x)sd(data.tmp[[x]])) beta <- b[1,]*x.sd/sd(data.tmp[[y.colnum]])
標準化回帰係数は大小比較が可能である。
> beta [1] 0.003351413 0.973894245
- 身長の標準化回帰係数…0.0034
- 年齢の標準化回帰係数…0.97
得点にあたえる影響は「年齢」 > 「身長」
相関係数
#重相関係数(目的変数と目的変数の期待値との相関係数） r <- cor(data$marks,expected.y) #決定係数 R2 <- 1 - e.SS / t.SS #自由度調整済み決定係数 adjusted.R2 <- 1 - e.MS / t.MS
- 重相関係数…目的変数と目的変数の期待値との相関係数
- 決定係数(寄与率)…回帰モデルの当てはまりの良さを示す(説明変数が多いと当てはまりは良くなる一方でモデルが複雑になる)
- 自由度調整済み決定係数…説明変数の数で調整をして、モデルの複雑さと当てはまりの良さのバランスを取る
> r [1] 0.9770006 > R2 [1] 0.9545302 > adjusted.R2 [1] 0.9542487
多重共線性とVIF
#説明変数ごとの重相関係数(自分自身を目的変数、自分以外の説明変数を説明変数とする重相関係数） r.age <- cor(data$age,data$height) r.height <- cor(data$height,data$age) #VIF(Variance Inflation Factor) vif.age <- 1/(1-r.age^2) vif.height <- 1/(1-r.height^2)
- VIFが10を超えると説明変数同士がぶつかり合う「多重共線性」が発生していることが疑われる。
※VIFの計算は説明変数が2つだと簡単だが、3つ以上になると改めて重回帰分析を繰り返さなければならないので、結構面倒。今回は変数の決め打ちで式を立てている。
> vif.age [1] 7.075908 > vif.height [1] 7.075908
少し値は大きいが、10未満なので問題なしとする。
Rでの重回帰分析
data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/exam-result.csv", fileEncoding = "utf-8")
summary(data)
#得点を目的変数、身長と年齢を説明変数とする重回帰分析
lm.result <- lm(marks~height+age,data)
#結果詳細を見る
summary(lm.result)
Call:
lm(formula = marks ~ height + age, data = data)
Residuals:
Min 1Q Median 3Q Max
-6.0221 -2.6579 0.1528 2.6190 6.9551
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.548134 2.609547 -2.126 0.0343 *
height 0.003497 0.032936 0.106 0.9155
age 6.641647 0.215236 30.857 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.013 on 323 degrees of freedom
Multiple R-squared: 0.9545, Adjusted R-squared: 0.9542
F-statistic: 3390 on 2 and 323 DF, p-value: < 2.2e-16