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

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

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

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

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

では、どうするべきか？

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

文系のための「擬逆行列」（２）

逆行列というのは計算がとにかく厄介で、手計算で計算するのは大変。

それ故に、通常は、コンピュータの力を借りることになる。
しかしながら、逆行列の性質を知らなくても良いということでは無い。

これまで、行列の足し算、引き算、掛け算を行なってきたが、
行列の「割り算」そのものは存在しない。その代わりに、存在するのが逆行列である。
逆行列は、スカラーにおける逆数に相当し、行列を基準化するために用いる。

つまり、逆行列は、厳密には「割り算」では無いけれど、
逆行列を掛ける（あるいは掛けられる）ということは、
スカラーにおいて逆数を掛けるような意味を持つ。

では、実際にどのような時に役立つのか？

実際には、様々なところで逆行列は登場するのであるが、
特に、逆行列を用いることで連立方程式が解ける、
というのは、逆行列の重要な性質の一つである。

私は辛党なのだが、ケーキ屋さんになったつもりになって、
洋菓子を作るための分量を考えてみることにする。
まず、以下のようなレシピとそれに対応する分量が与えられたとする。

	pancake	doughnut	sponge cake
cake flour (gram per piece)	50.0	7.0	75.0
egg (per piece)	0.5	0.2	3.0
milk (millilitre per piece)	50.0	0.0	20.0
sugar (gram per piece)	5.0	6.0	50.0

これは、ホットケーキ、ドーナツ、カステラのレシピであり、
それぞれ、一個辺りの薄力粉の量、卵の量、牛乳の量、砂糖の量が書いてある。
この通りに作れば、秘伝の洋菓子を作れる...かもしれない。

なお、これまでの話では、対象が行で変数が列となる「n×p行列」を扱ってきたが、
今回に関しては、対象が列で変数が行となる「p×n行列」となっている。
これは、後の話を理解しやすくするためであり、例外的なデータ構造となっている。

まずは、上記のレシピの分量を行列としてRに入れてみる。

# 各レシピにおける材料の分量をベクトルとして準備する
pancake <- c(50, 0.5, 50, 5.0)
doughnut <- c(7.0, 0.2, 0.0, 6.0)
spongecake <- c(75, 3.0, 20.0, 50.0)

# 作成したベクトルを結合して行列を作成する。
A <- cbind(pancake, doughnut, spongecake)

# 対称名を付ける
rownames(A) <- c("flour", "egg", "milk", "sugar")

ここで、cbind()関数というのは、ベクトルあるいは行列を列に追加するための関数。
この関数も非常に良く使う。列を追加（結合）するためにはrbind()関数を使う。

上手くできると、以下のような行列が出来るはず。

> A
      pancake doughnut spongecake
flour    50.0      7.0         75
egg       0.5      0.2          3
milk     50.0      0.0         20
sugar     5.0      6.0         50

まずは、復習から。とりあえず、行列の掛け算を用いて、
ホットケーキ2枚、ドーナッツ10個、カステラ1個
を作ったときの総分量はどのようになるのか？これを考えてみる。

とりあえず、小麦粉の総量例を考えてみると、
ホットケーキの小麦粉、ドーナツの小麦粉、カステラの小麦粉のそれぞれの分量と
作りたい個数を掛けて足せば良いのだから、この場合は、

小麦粉の総分量= (50.0 ✕ 2)＋(7.0 ✕ 10) ＋(75 ✕ 1)

となる。これは行列の掛け算で簡単に計算することができる。
行列の掛け算のルールは、「横✕縦」が原則であった。
つまり、以下のような式を考えれば良いことになる。



となるはず。Rで実際に計算して確かめてみる。
たしか、行列の「掛け算」のときには「%*%」を使うのだった。

# 作りたい個数を変数にセットする。
pieces <- c(2, 10, 1)

# 次に、行列の掛け算を実行
A %*% pieces

以下がこの行列の計算結果。

> # 次に、行列の掛け算を実行
> A %*% pieces
      [,1]
flour  245
egg      6
milk   120
sugar  120

このように、総分量を計算するのは至って簡単。
上手くやれば、洋菓子作りの分量計算を自動的にやってくれるソフトも作れる。
まぁ、現実問題として、そういったソフトの需要があるかという問題はあるが。

さて、以上のことが確認できたら、今度は、逆の場合を考えてみる。
すなわち、総分量があったとして、それで作れるお菓子の量を知りたい、
という状況を考えてみる。

薄力粉565g、卵15個、牛乳260ml、砂糖290gを使ったときに、
ホットケーキ、ドーナッツ、カステラはそれぞれいくつ作れるのか？
今度は、これを考えてみる。まずは、式を考えなくてはならない。



知りたい場所が変わったので、今度はこのような式を考えることができる。
おそらく、ここまでは、想像できたはず。知りたい場所が変わっただけ。
計算はしないが、試しに、スカラーの計算に展開してみる。



そうそう。このようなものを連立方程式と呼ぶのであった。
我々が中学校までに習った連立方程式は、実は行列で表すことができる。
計算の順番を確認してみると、ちゃんと、横✕縦の形式になっていて対応がとれるハズ。

では、いよいよ、この連立方程式を解くことを考えてみる。
ここで、先程の行列の式を、よく観察してみる。
すると、あることに気づく。

左辺にxのベクトルを残し、右辺を基準化したら答えが解るのではないか？
ここで、スカラーの世界では、両辺を同じ数で割れば良かったはずだが、
残念ながら、行列の世界には「割り算」は存在しない。そこで、逆行列が登場する。



思わず割り算をしたくなるかもしれない。でも、割り算をしてはいけない。
移動したレシピの行列を良く見てみる。右上に「-1」が付いている。
この右上の記号が逆行列を意味するのであった。

さて、もう一つ思い出してもらいたい。通常の逆行列の条件は何だったか？
既に述べたように、逆行列は正方行列にしか対応しなかった。
つまり、ここでは、擬逆行列を使う必要が出てくる。

今の状況では、少々解りにくい。記号を使って式を般化してみる。



元の行列の転置行列と元の行列を掛けると正方行列になる。
この正方行列の逆行列は、転置行列の分だけ余分に基準化しているので、
さらに、転置行列を掛けて、元の状態に対応する行列に戻す。

また、特異値分解から、次のようにして計算することも可能。



詳しくは、特異値分解の話に書いてあるが、
スカラーにおける因数分解のように、行列においても分解というものが存在する。
要するに、行列の状態の見通しを良くするのであるが、これを使うと擬逆行列も求まる。

とりあえず、この二つの方法で、先程の問題を解決してみる。
まずは、転置行列を使った方法。

# 材料の総量
total <- c(565, 15, 260, 290)

# 連立方程式を解く
solve((t(A)%*%A))%*%t(A)%*%total

ここで、solve()関数というのが逆行列を求める方法であった。
以下がこの計算の結果。

> # 連立方程式を解く
> solve((t(A)%*%A))%*%t(A)%*%total
           [,1]
pancake       4
doughnut     20
spongecake    3

次に、特異値分解を用いた方法。

# 材料の総量
total <- c(565, 15, 260, 290)

# 特異値分解を行う
A.svd <- svd(A)

# 解りやすいように変数を入れなおす。
D <- diag(A.svd$d)
V <- A.svd$v
U <- A.svd$u

# 連立方程式を解く
V %*% solve(D)%*%t(U)%*%total

そして、この結果。

> # 連立方程式を解く
> V %*% solve(D)%*%t(U)%*%total
     [,1]
[1,]    4
[2,]   20
[3,]    3

そうそう。逆行列の関数ginv()関数も存在していた。一応、それも比較。

# ライブラリを読み込む
library(MASS)

# 連立方程式を解く
ginv(A)%*%total

そして、この結果を以下に。

> # 連立方程式を解く
> ginv(A)%*%total
     [,1]
[1,]    4
[2,]   20
[3,]    3

以上の結果から、いずれの方法も同じ結果を導いていることが確認できる。
結論としては、ホットケーキ4枚と、ドーナツ20個、カステラ3個を作ることができる。

ところで、今回の場合は、きちんと計算することができたが、
かならず、キレイな計算結果が出てくるとは限らない。
つまり、そもそも、連立方程式が解けない場合が考えられる。

実は、そのような時でも、擬逆行列を使うと近似的に最適な答えを見つけてくれる。
最適な答えというのは、要するに、誤差の二乗が最も小さくなるような答えである。
この性質から、一連方法は、「最小二乗法」ととも呼ばれる。

このことは、後に中級編で解説する「重回帰分析」とも密接な関わりがある。