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

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

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

とりあえず、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」近い。
ひょっとしたら、モデルに加えない方が良いかもしれない。
あるいは、後円高や前方幅もいらないかもしれない。

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

では、どうするべきか？

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