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

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

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



一般的に、擬逆行列の解法は変数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検定という方法を用いることでモデルの適合性を評価するができる。
これに関しては、次回に整理することにする。

文系のための「重回帰分析のモデル」

重回帰分析の基本的な仕組みを理解したところで、
今度は、実際のデータを用いて重回帰モデルの検討方法を整理する。

今回のデータは、某国の某時代の某地域の「古墳」のデータ。
実際のデータを少し加工し、詳細が特定できないようにマスキング。
考古学的な解釈を避けるための処置。興味がある人は自分のデータで。

とりあえず、Google spreadsheet にデータを置いているので、
インターネット経由で直接データを読み込む。
この作業を行うためには、RCurlというライブラリを読み込む必要がある。

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=","))

このデータは、「294」の古墳が「対象」となっている。
属性は、墳長（length）、後円径（bk_radius）、後円高（bk_height）、前方幅（fr_width）、
前方長（fr_length）、前方高（fr_height）、くびれ（constricted）である。

なお、このデータの単位は、全て「メートル」であり、計測値は、10センチ単位まで。
扱うデータの単位を理解しておくことは基本中の基本なのであるが、
単位を無視して分析する人は少なくは無い。

ふむ。どこの国かは、特定できそうな変数であるが...まぁ、それは置いておく。
細かい詮索は無用。とにかく、これは架空のデータとして理解。
ちなみに、かつて、卒業論文研究のために集めたデータを改変したもの。

なお、後の処理のために、今回のデータはマトリクス型では無く、
データフレーム型というデータ型を指定しておく。

まずは、このデータの相関係数行列と偏相関係数行列を確認しておく。
多次元データを扱う際には、最初にデータの全容を掴んでおくことが重要。
やり方は、理解できているハズなので、出力結果のみを表示。

> round(cor(X),3)
            length bk_radius bk_height fr_width fr_length fr_height constricted
length       1.000     0.985     0.933    0.949     0.976     0.873       0.955
bk_radius    0.985     1.000     0.944    0.930     0.939     0.865       0.953
bk_height    0.933     0.944     1.000    0.882     0.887     0.896       0.900
fr_width     0.949     0.930     0.882    1.000     0.956     0.895       0.972
fr_length    0.976     0.939     0.887    0.956     1.000     0.866       0.950
fr_height    0.873     0.865     0.896    0.895     0.866     1.000       0.872
constricted  0.955     0.953     0.900    0.972     0.950     0.872       1.000

小数点以下の桁数が多すぎるので、小数点第三位までで丸め込み。
対象数が294の場合の相関の基準も考えておく。
たしか、相関係数のt検定というのがあった。

# 対象数を変数「n」に格納する。
n <- nrow(X)

# 有意水準1%、5%、10%の場合について検討する
abs(qt((0.01)/2, n-2)/sqrt(n-2+qt((0.01)/2, n-2)^2))
abs(qt((0.05)/2, n-2)/sqrt(n-2+qt((0.05)/2, n-2)^2))
abs(qt((0.10)/2, n-2)/sqrt(n-2+qt((0.10)/2, n-2)^2))

以下が実行結果。

> # 有意水準1%、5%、10%の場合について検討する
> abs(qt((0.01)/2, n-2)/sqrt(n-2+qt((0.01)/2, n-2)^2))
[1] 0.1500135
> abs(qt((0.05)/2, n-2)/sqrt(n-2+qt((0.05)/2, n-2)^2))
[1] 0.1144192
> abs(qt((0.10)/2, n-2)/sqrt(n-2+qt((0.10)/2, n-2)^2))
[1] 0.09611705

以上の結果から、有意水準1%で見ても、相関係数が「0.16」以上であれば、
相関があると行っても構わない。相関係数行列を見てみると、
どの値もこの水準よりも遥かに高い。したがって、全ての組合せには相関がある。

この状態では見難いかもしれない。ついでに、散布図行列も出しておく。

# 「凝った」散布図行列を作るためにライブラリ「psych」を使う
library(psych)
pairs.panels(X, smooth=FALSE, density=FALSE, ellipses=FALSE, scale=TRUE)

なるほど。散布図行列を見てみても、なんらかの関係があるように見える。
ひょっとしたら、見せかけの相関が潜んでいるかもしれない。
ということで、次に、偏相関係数行列も計算しておく。

> library(corpcor)

> round(cor2pcor(cor(X)),3)

       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]
[1,]  1.000  0.835  0.113  0.131  0.816 -0.078 -0.298
[2,]  0.835  1.000  0.249 -0.124 -0.586 -0.045  0.444
[3,]  0.113  0.249  1.000 -0.200 -0.103  0.540  0.073
[4,]  0.131 -0.124 -0.200  1.000  0.129  0.415  0.654
[5,]  0.816 -0.586 -0.103  0.129  1.000  0.068  0.270
[6,] -0.078 -0.045  0.540  0.415  0.068  1.000 -0.134
[7,] -0.298  0.444  0.073  0.654  0.270 -0.134  1.000

Rの偏相関係数行列は、ライブラリ「corpcor」が必要であり、
cor2pcor()関数は、相関係数行列の結果を使って、
偏相関係数を計算するのであった。ラベルは、相関係数行列と同じ。

偏相関を見てみると、墳長と後円径、墳長と前方長の関係が際立って高い。
さて、この偏相関係数の値は、後に重要となるのであるが、
今は、深く踏み込まず、実際に重回帰分析の方法を整理する。

さて、このデータを用いて重回帰モデルを立てるためには、
最初に目的変数を決めなければならないのであるが、
墳長（length） を目的変数とし、その他の変数との関係を調べてみる。

最初にすべきことは、何であったか？そうそう、重回帰係数を求めること。
一応、前回の復習を兼ねて、重回帰モデルの基本形を見てみる。



元のデータのうち、一つの変数が目的変数となり、
切片を計算するために、先頭に「1」が並んだ列が入っている。
このような行列を連立方程式に見たてて解いてやると重回帰係数が出てくるはず。

連立方程式を解くためには、逆行列が必要となるが、
残念ながら逆行列は正方行列にしか使えない。どうれば良かったか？
たしか、擬逆行列を用いるのであった。

二つの方法を紹介したが、ここでは特異値分解を使う。
既に述べたようにn✕p行列とp✕n行列では解き方が異なる。
そのような場合でも、特異値分解を用いれば擬逆行列が計算できる。

重回帰分析を実行できるソフトによっては、変数の数が対象数よりも多い場合に、
デタラメな結果を表示することもあるらしい。注意が必要である。
擬逆行列と重回帰分析の仕組みを知っていれば解ることであるが...。

いずれにせよ、対象数よりも変数の方が多いようなデータに対して、
いわゆる「多変量解析」を行うのは危険極まりないのである。
上記の意味が解らない人は、もう一度、擬逆行列の話を参照すること。

さてさて、特異値分解を用いた方法はどのような方法であったか？
これも、もう一度確認しておくことにする。



Rには、特異値分解を行うためのsvd()関数がある。
この関数を覚えておくと主成分分析など他の分析でも役立つ。

少々、話が逸れてしまったが、上記のことを確認した上で、
実際に、以下の手順で重回帰係数を求めてみる。

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

# 目的変数を元の行列から分離する
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.svd <- svd(X)
X.svd.v <- X.svd$v
X.svd.d <- diag(X.svd$d)
X.svd.u <- X.svd$u

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

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

# 最後に結果を表示する
X.lm.coefficient

実行結果は以下の通り。

> # 最後に結果を表示する
> X.lm.coefficient
            coefficients
Intercept     1.31774995
bk_radius     1.05513862
bk_height     0.41494162
fr_width      0.09120632
fr_length     0.90234517
fr_height    -0.23795594
constricted  -0.33198797

これでめでたく重回帰係数が求まったはずだが...本当に計算できているのか？

Rにおける重回帰分析を行うための関数は、単回帰と同じ。

どのような関数であったか？頑張って思い出してみる。

# Rで重回帰分析を行なってみる。
lm(formula=length~bk_radius+bk_height+fr_width+fr_length+fr_height+constricted, data=O)

実は、単回帰の話で用いたlm()関数で良い。

すでに、Xのデータは、元のデータと変わってしまっているので、

ここでは、元のデータをコピーである行列「O」を用いている。

lm()関数では、少々、見慣れない形でデータが指定されている。

まず、パラメータ「formula」をモデル式と呼び、ここでモデルを定義し、

実際のデータをパラメータ「data」で指定する。

この場合は、「〜」を挟んで右側が目的変数であり、左側が説明変数となる。

今回は、全ての変数を使っているが、使う変数をここで決めることができる。

なお、全ての変数を使う場合には、次のようにも指定できる。

lm(formula=length~., data=O)

もちろん、上記二つの結果は同じになる。
とにかく、実行結果は以下の通り。

> lm(formula=length~., data=O)

Call:
lm(formula = length ~ ., data = O)

Coefficients:
(Intercept)    bk_radius    bk_height     fr_width    fr_length    fr_height  
    1.31775      1.05514      0.41494      0.09121      0.90235     -0.23796  
constricted  
   -0.33199 

重要な部分は、「Coefficients: 」から後の部分。この部分に重回帰係数が並んでいる。
特異値分解を使って解いた場合と結果が同じことを確認できたはず。
重回帰分析の解法を理解できていれば、lm()関数を用いても問題ないだろう。

さて、ここで、重回帰係数の検討に入りたいところであるが、
その前に、もう少し、重回帰モデルそのものについて考えてみたい。

何度も述べているが「回帰」というのは変数間の関係を観察するための方法である。
では、変数間に「全く」関係が無い場合、一体どのようなことになるのか？
ここで、簡単な実験を行なってみる。

# 乱数使って変数間に関係の無いデータを作る
test <- matrix(c(runif(1000),runif(1000)), ncol=2)

# 回帰係数を求める
b0 <- lm(test[,2]~test[,1])$coefficient[1]
b1 <- lm(test[,2]~test[,1])$coefficient[2]

# データをプロットする。
rng <- c(0, 1)
plot(test[,1], test[,2], ylim=rng, xlab="x", ylab="y", pch=16, col="lightblue")
par(new=TRUE)
plot(test[,1], test[,1]*b1+b0, type="l", ylim=rng, xlab="", ylab="", xaxt="n", yaxt="n", col="blue")

# yの平均を表す水平線を引く
abline(h=mean(test[,2]), col="red", lty=2, ltw=2)

さて、この図を見て解るように、関係が無いということは、
傾きは限りなく「0」に近づき、切片のみが反映されていることが解る。

たしか、xの平均値のときにyの平均を通るように調整した際のyの値が
切片であったので、結局は、平均によって説明することと同じなのである。
この状況は、lm()関数のモデル式でも表すことができる。

ここで、もう一度、古墳のデータに戻って確認してみる。

# まずは、行列「O」の「length」の平均値を計算する。
mean(O$length)

# 次に、lm()関数で変数を指定しない場合の係数を確認する。
lm(formula=length~1, data=O)$coefficient

以下がこの実行結果。

> # まずは、行列「O」の「length」の平均値を計算する。
> mean(O$length)
[1] 83.86105
>
> # 次に、lm()関数で変数を指定しない場合の係数を確認する。
> lm(formula=length~1, data=O)$coefficient
(Intercept)
   83.86105

つまり、他の変数との関係を除いてしまえば、平均によって説明することと同じ。
また、作成した回帰モデルが本当に意味があるものであれば、
そのモデルは平均よりも妥当性が高い（説明力が高い）とも言える。

さて、ここからが、今回の話の「本題」である。
いやいや、ここまでの話は、少々長すぎた...かもしれない。

すなわち、重回帰モデルの検討の仕方。ここまでのことが理解できていれば簡単。
とりあえず、重回帰係数の絶対値の大きさと、その符号を観察する。



このモデルの中で、負の符号を取っているのは、前方高とくびれ幅。
したがって、前方高の値が「1メートル」増えると、墳長は「0.238メートル」減り、
同様に、くびれ幅の値が「1メートル」増えると、墳長は「0.332メートル」減る。

逆に、全長が増える場合を考えると、正の係数が最も大きいものが後円径であり、
この値が「1メートル」増えると、墳長は「1.055メートル」増えることになる。
次に値が大きい前方長が「1メートル」増えると、「0.090メートル」墳長は増える。

さて、このモデルでは、小数点第三位まで表しているのであるが、適切であろうか？
元のデータの単位がメートル単位であり、値は10センチ刻みであった。
したがって、ミリ単位というのは細かすぎる。本当は、小数点第二位までで十分。

たまに、コンピュータの出力をそのまま書き出している人がいるが、
計算結果の桁の下の方は、ほとんど計算誤差。
意味のある桁数を理解しておく必要がある。

話が逸れたが、このようにして、ある変数とその他の変数の間の関係を、
それぞれの値の増減によって変化するかを観察するのが重回帰分析である。

ところで、この結果を見て、一種の「違和感」を感じた人もいるかもしれない。
規模の話なので、前方高やくびれ幅が増えると、墳長が減って良いのか？
ここで、もう一度、偏相関係数の結果を見てみる。

> round(cor2pcor(cor(X)),3)
       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]
[1,]  1.000  0.835  0.113  0.131  0.816 -0.078 -0.298
[2,]  0.835  1.000  0.249 -0.124 -0.586 -0.045  0.444
[3,]  0.113  0.249  1.000 -0.200 -0.103  0.540  0.073
[4,]  0.131 -0.124 -0.200  1.000  0.129  0.415  0.654
[5,]  0.816 -0.586 -0.103  0.129  1.000  0.068  0.270
[6,] -0.078 -0.045  0.540  0.415  0.068  1.000 -0.134
[7,] -0.298  0.444  0.073  0.654  0.270 -0.134  1.000

実は、変数間の関係性で言うと、特に、前方径の値は限りなく「0」近い。
ひょっとしたら、モデルに加えない方が良いかもしれない。
あるいは、後円高や前方幅もいらないかもしれない。

また、くびれ部に関しては、全体の大きさに対して、
変化の度合いが鈍感で、そのような状況が影響しているかもしれない。
とにかく、何とも言えない状況がある。

では、どうするべきか？

この問題を考えるためには、重回帰係数の妥当性の検討や、
どの変数を選択すべきであるか、という変数選択が必要になる。
ということで、次回は、重回帰係数の妥当性から整理してみることにする。

文系のための「二変数の関係」（３）

これまでの話から、二つの変数間の関係を見るための方法には、
相関係数以外にも単回帰という方法があることが理解できたハズ。

単回帰は、相関係数に良く似ているが、
一方の変数によってもう一方を説明するという考え方に基づいている。
相関係数は、二つの変数を平等に見ていたので、この点が大きく異なる。

また、重要なこととして、相関係数は二つの変数が平等なので軸は関係無いが、
単回帰の場合は、どちらを説明変数に置くかによって値が異なる。
まぁ、この辺りは、感覚的に理解できていることだろうと信じたい。

さて、今回は、もう少し単回帰について踏み込んでみる。
実は、前回の話では、二つの変数間の関係の強さについては、
十分に検討ができていなかったのである。

相関係数と同様に、二つの変数間の関係の強さを観察するには、
推定された値と実際の値との誤差の部分の解釈が必要となる。
本日は、このことについて整理する。

ということで、前回のデータの読み込みから。

# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",", header=TRUE)

# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",", header=TRUE)

以下が、そのデータ。データの説明に関しては、平均の話を参照。

Oroshi_Kg, Taka, Naka, Yasu
Tai, 1336, 3150, 1217, 53
Chinu, 226, 630, 513, 105
Sawara, 212, 1575, 1259, 840
Yazu, 4834, 840, 315, 179
Suzuki, 618, 1575, 1146, 525
Akou, 21, 4725, 3129, 1575
Kochi, 28, 2100, 874, 210
Okoze, 28, 5250, 2635, 210
Ainame, 29, 3150, 621, 315
Hage, 1205, 1890, 383, 210
Konoshiro, 21, 2100, 1100, 105
Sayori, 12, 5250, 3916, 1260
Anago, 1637, 2625, 1721, 630
Mebaru, 330, 4200, 1680, 315
Tachiuo, 1211, 2310, 930, 105
Hamo, 1622, 3938, 592, 53
Karei, 41, 6300, 3178, 525
Hirame, 540, 2100, 1496, 210
Aji, 5149, 3150, 789, 210
Saba, 5413, 3675, 625, 53
Ika, 2414, 2888, 1083, 105
Tako, 1844, 1890, 1180, 525
Ebi, 560, 7350, 1420, 525
KurumaEbi, 222, 8925, 6508, 2625
Koiwashi, 1316, 1155, 617, 328
ObaIwashi, 2716, 499, 267, 175
Mentai, 678, 1155, 859, 394
Hamachi, 4772, 998, 632, 105
Isaki, 268, 2625, 1465, 840

ふむ。前回と同様に、中値と安値の関係を使うことにするか。
単回帰分析も行うことにしよう。lm()関数があるが今は使わない。

x1 <- X$Naka
x2 <- X$Yasu

# まず、推定値を求める
B <- cov(x1,x2)/var(x1)

# まず、推定値を求める
B0 <- mean(x2) - B*mean(x1)

# x2の推定値を求める
x2.hat <- B0 + B*x1

数式あった方が解り易い。x1によってx2の値を説明するので、
たしか、以下のような式になるハズであった。



ここで、二つ推定値のいうのがそれぞれ、

,  

であった。詳しい説明は既にしてあるから、ここまでは大丈夫であろう。
自信が無い人は、もう一度復習をすること。

さて、現在、この回帰式によって説明されるのは、
あくまで、理屈として推定される値であって、実際の値とは異なる。
理屈上の話だから、直線で推定の直線が描けたのであった。

では、実際の値というのは、推定された値を使ってどのように表現できるか？



この式は、推定値では無くて、実際の値を表しているから、
左辺の「^（ハット）」が取れていることにも注目。

さて、この式の中で重要なことは最後に「ℇ」の値が付いていること。
これを「残差（residuals）」と呼び、この記号は一般的に残差を表す。

つまり、推定された値に誤差を加えたものが実際の値、と考えるのである。
この辺りの考え方は、数学嫌いの人には、少々、理解し難いかもしれない。

では、残差は、どうすれば求めることができるか？
これも中学校までの数学の知識があれば簡単なことであろう。



つまり、実際の値から推定値を引くだけ。混乱することは無い。
では、早速、この誤差の部分を計算してみる。

residuals <- x2 - x2.hat

ふむ。これで誤差が計算できた。ちょっと解りにくいので可視化する。

# 上手く重なるようにx軸とy軸の描画範囲を固定する
x.lim <- c(0,max(x1))
y.lim <- c(0,max(x2))

# 散布図を作成する
plot(x1, x2, xlim=x.lim, ylim=y.lim, pch=16)

# 図を重ねるためのオマジナイ
par(new=TRUE)

# xの値の範囲を最低値と最大値で決める。
x1.range <- seq(min(x1),max(x1))

# x値の範囲に対応する推定値を計算する。
x1.estim <- B0 + B*x1.range

# 回帰直線を引く
plot(x1.range, x1.estim, xlim=x.lim, ylim=y.lim, xaxt="n", yaxt="n", xlab="", ylab="", type="l", col="red")

# 残差を表す線を描画する
for(i in seq(1,length(x1))){
  x <- rep(x1[i],2)                 # x1の値を2つ作る
  y <- c(x2[i], x2[i]-residuals[i]) # x2の値とx2の値から残差を引いた値を作る

  # 残差を表す垂線を引く
  par(new=TRUE)
  plot(x, y, xlim=x.lim, ylim=y.lim, xaxt="n", yaxt="n", xlab="", ylab="", type="l", col="blue")
}

# 最後にグリッドを描いて全体を整える
grid()

この図において、実際の値から回帰直線に引かれた垂線が残差「residuals」である。
この残差というのが、二つの変数間の関係を観察する上で重要な値で、
残差を合わせた値が小さいほど当てはまりが良いことを表す。



二乗しているのは、値の正負が打ち消し合わないようにするため。
これを残差平方和（Sum Squared Error）と呼ぶ。とりあえず、Se と置いておく。
試しに、今回のデータを使って、この値を確認してみる。

(se <- sum(residuals^2))

これを実行すると、以下のようになる。

> (se <- sum(residuals^2))
[1] 2076206

さて、この値が小さいということが重要なのだけれど、これは大きいのか？
この値の大小の程度を見分けるためには、やはり、基準化しないといけない。
原点に帰って、もう一度、この残差というものが何なのかを考えてみる。

データを観察する上で、最も重要なのは、くどいけれど、バラツキ。
そして、個々の値のバラツキを偏差と言い、
偏差を二乗して足し合わせたものが偏差平方和。

とりあえず、目的変数の総平方和（Sum Squared Total Deviation）と呼び、
この値をSt とでも置いておくことにする。式は以下の通り。



これを対象数から1を引いた「n-1」で割ると？そうそう、分散になる。
要するに、分散の式の分子の部分に当たるのであった。

さて、このStとSeの関係に注目してやると、
Stは観測値のバラツキの平方和で、Seは残差のバラツキの平方和だから、
予測値のバラツキとSe を足したものが、St になるはず。

では、予測値のバラツキはどうなるか？というと、以下のようになる。



つまり、St = Se + Sr という関係が成立するはず。
これは、直感的に理解できる。ここからが数学マジック！
両辺をStで割るとどういうことになるのか？



となる。ここで、残差平方和の影響の大きさを見たいので、
この式から、予測値のバラツキの部分を抜いてしまう。
すると、以下のようになる。左辺に何もないと困るので「」を置いておく。



つまり、残差平方和の大きさというのは、このように表すことができる。
この値が「1」に近づくということは、残差の影響が小さいことを意味する。
なかなか、エレガントな数式。私は美しいと思うのだが、どうだろう？

さて、この「」という指標は頻繁に登場する。だからちゃんと名前がある。
一般的に、この指標のことを「決定係数」と呼び、慣例的に「」で表す。

では、この決定係数を実際に計算してみる。

# 一応、三つの値を出しておく。
st <- sum((x2-rep(mean(x2), length(x2)))^2)
sr <- sum((x2.hat-rep(mean(x2), length(x2)))^2)
se <- sum((x2-x2.hat)^2)

# 決定係数の計算
(R2 <- 1-(se/st))

では、確認。以下が実行結果。

> # 決定係数の計算
> (R2 <- 1-(se/st))
[1] 0.7576962

さて、この結果は果たして高いと言えるのかどうか？実は、判断は難しいのであるが、
一方の変数がもう一方を説明する上で、75.8％の説明力を持っていることを示していて、
これは、決して低い値ではない。大体、60〜70％以上であれば、悪くはない。

決定係数に関しては、注意するべきことがある。
相関係数と同様に、この値を盲目的に信用してはならないし、
慎重にデータの解釈をしなければならない。決定係数は、基準の一つでしかない。

例えば、決定係数は、対象数が多くなると必然的に値が高くなるので、
自由度調整済み決定係数と呼ばれる指標を用いたり、
確立統計の手法を用いる場合もある。

さてさて、最後に、Rのlm()関数を使った場合についてもやっておく。
とりあえず、X.reg という変数に、回帰分析の結果を入れて、その要約を確認する。

X.reg <- lm(x2 ~ x1)
summary(X.reg)

以下がこの処理の実行結果。

> X.reg <- lm(x2 ~ x1)
> summary(X.reg)

Call:
lm(formula = x2 ~ x1)

Residuals:
    Min      1Q  Median      3Q     Max 
-678.72 -163.58   -7.29  158.81  506.60 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -69.68447   77.21191  -0.903    0.375    
x1            0.36372    0.03958   9.189 8.45e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 277.3 on 27 degrees of freedom
Multiple R-squared: 0.7577, Adjusted R-squared: 0.7487 
F-statistic: 84.43 on 1 and 27 DF,  p-value: 8.447e-10 

実は、色々と見ないといけないのであるが、
そもそも、今回の話の目的は、二つの変数間の関係を見ることであって、
精密な推定を行うことではない。ということで、決定係数以外は無視。

さて、上記の出力において、Coefficients の部分は係数の部分。
決定係数の部分は、下から二行目の部分。Multiple R-squared の右側。
ちゃんと、手計算で算出した決定係数と一致していることが解る。

回帰分析の解説の多くは、「目的変数を説明変数によって推定する手法」、
と紹介しているのではないかと思うが、変数間の関係を見ることが重要なのである。
これは、相関係数の考え方に非常に良く似ているのである。たしかに、式も似ていた。

相関係数との違いは、一方の変数によってもう一方の変数を説明するという考え方。
したがって、その説明力を「残差」に注目して検討するのである。
残差が小さいほど、二つの変数の間の関係は強いということになる。

他の教科書に書かれているように、推定に注目するのであれば、
より複雑なことが必要になるし、3つ以上の変数間の関係になると結構大変。
残差そのものに関する検討や、複数の変数の当てはまりの高さの検討も必要。

この章で書いてしまいたい話ではあるが、
この辺りの話は、中級編の重回帰分析の話で整理することにする。