文系のための「自由度調整済み決定係数」

ここまでの話では、重回帰分析の仕組みとモデルの読み方を整理してきた。
要するに、重回帰係数について、連立方程式を解き、
重回帰係数の値の正負と値の大小によって変数間の関係を観察する。

そうそう。特異値分解による擬逆行列から連立方程式を解く場合には、
以下のような式によって、重回帰係数を算出することができた。
少々、クドいかもしれないが、とりあえず再掲しておく。



一般的に、擬逆行列の解法は変数pが対象数nよりも小さい場合と（n > p）、
変数pが対象数nよりも大きい場合（n < p）で異なっているが、
特異値分解を用いた場合には近似的に計算することができる。 

とにかく、この式によって導出されたベクトルが重回帰係数であり、
この重回帰係数を用いて、次のような重回帰モデルを導出することができる。



さて、基本的な話は、これで収まっているのであるが、 問題となるのは、
このモデルが本当に正しいのか？という問題である。
まずは、このモデルの適合の度合いというか、説明力を知りたい。

最も単純な方法は、残差から重回帰モデルの適合性を評価する決定係数を用いる方法。
 予測値のバラツキ（回帰の平方和）「Sr」、残差のバラツキ（残差平方和）「Se」、
そして、目的変数の総平方和 「St」から以下の式によって表すことができる。



この式を変形し、以下のようにしたものが決定係数と呼ばれるものであった。



詳細に関しては、残差と決定係数の話ですでに整理済みの話である。
決定係数は、モデルで説明できない残差の大きさを反映した指標である。

ところで、この指標は直感的に理解しやすいのであるが、少々、問題がある。
実は、変数の数が増えれば増えるほど決定係数の値は高くなってしまうのである。
そこで登場するのが、「自由度修正済み決定係数」と呼ばれるものである。



ここでは、目的変数の総平方和「St」に対して対象数「n-1」で基準化し、
Seは、Stと同様に対象数で基準化すると同時に、説明変数の数「p」で基準化をしている。
つまり、全体では一つの対象における一つの変数当たりの残差でモデルを評価する。

このような調整を行うことで、対象数と変数の数の影響が除去される。
とにかく、前回のデータを用いて、自由度調整済み決定係数とやらを計算してみる。
まず、前回と同様に、インターネット経由で直接データを読み込む。
今回の演習で使用するデータも前回と同じ。

library(RCurl) 

ライブラリを読み込めたら、
以下のコマンドを実行する。

# Google Spreadsheet からデータをダウンロードする
data <- getURL("https://docs.google.com/spreadsheet/pub?key=0AtOtIs5BjRVhdDdyTmgtNm9DcEVJb3VLbnh3ZFlNb2c&single=true&gid=0&range=A1%3AH295&output=csv")

# ダウンロードしたデータを読み込む
X <- as.data.frame(read.table(textConnection(data), header=TRUE, row.names=1, sep=","))

次に、重回帰分析を実行する。前回と同様に、「墳長」を目的変数とし、
その他の変数から重回帰モデルを立てる。

# 念の為に最初のデータの状態を保存しておく。
O <- X

# 実際の目的変数も保存しておく。
y <- O[, 1]

 # 目的変数を元の行列から分離する
y <- X[, 1]
X <- as.matrix(X[, -1])

# 切片の計算のために先頭行を「1」の行列を付け加える。
Intercept <- matrix(rep(1, nrow(X)), ncol=1)
X <- cbind(Intercept, X)

# 特異値分解を行う
X.svd <- svd(X)
X.svd.v <- X.svd$v
X.svd.d <- diag(X.svd$d)
X.svd.u <- X.svd$u

 # 重回帰係数を求める
X.lm <- ((X.svd.v %*% solve(X.svd.d)) %*% t(X.svd.u)) %*% y 

# 解りやすいように名前を付ける。
rownames(X.lm) <- c("Intercept",colnames( X[,2:ncol(X)]))
colnames(X.lm) <- "coefficients" 

# 最後に重回帰モデルから予測する。
# なお、以下の計算はベクトルの成分の掛け算。
y_hat <- X.lm[1] + X.lm[2]*X[,2] + X.lm[3]*X[,3] + X.lm[4]*X[,4] + X.lm[5]*X[,5] + X.lm[6]*X[,6] + X.lm[7]*X[,7]

# まずは、StとSeを求めて、決定係数を算出する。
Se <- sum((y-y_hat)^2)
St <- sum((y-mean(y))^2)

# 以下が標準的な決定係数。
1-(Se/St) 

# 次に、自由度調整済み決定係数を求める
p <- ncol(X)-1
n <- length(y)
1-(Se/(n-p-1))/(St/(n-1)) 

実行結果は以下の通り。決定係数と自由度調整済み決定係数はかなり高い。 

> # 以下が標準的な決定係数。
> 1-(Se/St)  
[1] 0.9926216
> 
> # 次に、自由度調整済み決定係数を求める
> p <- ncol(X)-1
> n <- length(y)
> 1-(Se/(n-p-1))/(St/(n-1))  
[1] 0.9924674

決定係数は、モデルの当てはまりの良さを表す指標であるので、
今回の結果の場合は、モデル全体では99%以上の説明力を持っていることが解る。
したがって、今回のモデルは、非常に当てはまりの良いモデルだと言える。

さて、残差の検討もしておく。これも、必ず確認しておく必要がある。

たしか、残差プロットは、まんべんなく散らばっていた方が良かった。

残差に何らかの傾向が現れているということは、重回帰モデルとして問題がある。

# 残差を標準化してバラツキ具合を観察する。
plot(scale(y-y_hat), ylab="Residuals")
abline(h=0)
abline(h=c(-2, 2), col="gray", lty=2)

# 外れ値を検出する
(r.lower <- which(scale(y-y_hat) < -2)) # −2よりも小さいもの
(r.upper <- which(scale(y-y_hat) > 2))  # +2よりも大きいもの
1-((length(r.upper)+length(r.lower))/n)

そして、実行結果は以下の通り。

# 外れ値を検出する
> (r.lower <- which(scale(y-y_hat) < -2)) # −2よりも小さいもの
[1] 91 98 137 243 258
> (r.upper <- which(scale(y-y_hat) > 2)) # +2よりも大きいもの
[1] 53 94 121 127 131 162 164 263 264 267
> 1-((length(r.upper)+length(r.lower))/n)
[1] 0.9489796 

一般的に重回帰分析における残差は正規分布に従うとされている。
そして、標準正規分布に従うのであれば、±2σ に95.4%が入ることになるので、
簡易な方法として、標準化した値の±2の外側を「外れ値」とする。 

今回の場合、全体の約94.9%が、±2σの範囲内に収まっていて、
また、散布図を見る限りは、何の傾向も見えてこない。
したがって、やはり、このモデルは非常に良い状況を表していることになる。 

自由度調整済み決定係数を用いることで、モデルの適合の度合いを知ることができ、
さらに、F検定という方法を用いることでモデルの適合性を評価するができる。
これに関しては、次回に整理することにする。