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

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

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

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

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

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

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

# 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つ以上の変数間の関係になると結構大変。
残差そのものに関する検討や、複数の変数の当てはまりの高さの検討も必要。

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

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

データを観察する上で最も重要なことは、何度でも言うけれども、
「データの中心」と「データのバラツキ」を把握すること。
あらゆる統計的な手法はこの二つの点が全ての原点になっている。

そう言えば、数量化理論を創り上げた林知己夫氏が、
以下のように語っていることを思い出した。

「科学は所詮平均値の議論である。・・・（中略）・・・
平均や構造をみることで大きな知見を得ることができ、
さらに個々を介して新たに発展の道をたどることになる。」

この言は、日本行動計量学会35年記念誌『行動計量学』2008年9月号の、
「論文　ーその世俗的な話しー」と題したエッセイに書かれてある。
思わず、苦笑いしたくなる名言が詰まっていて、是非、研究者には読んでもらいたい。

初っ端から話が逸れてしまったが、重要なことは、あくまでデータの中心であって、
その中心からのバラツキを測った指標として、分散や標準偏差があった。

さらに、このバラツキを二つの変数間で観察するための指標として、
共分散というものがあり、それを基準化したものとして相関係数があった。
ふむふむ。では、相関係数とはどのようなものだったか？

確か、二つの変数の関係の強さを-1〜1の範囲で表し、

一方が増加するともう一方が減少する場合を負の相関と呼び、
一方が増加するともう一方も増加する場合を正の相関と呼んだ。

そうそう。そういった関連性の強さを表すものであった。

実を言うと、このように二つの変数の関連の強さを観察する方法というのは、
相関係数だけではない。別の考え方に基づく方法がある。
今回の話は、その方法についての話。

とりあえず、相関の話で使ったデータの読み込みから。

# 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

最初にするべきことは、やはり相関係数を出すこと。
話を解りやすく進めるため、今回は「中値」と「安値」で見ていくことにする。
この二つの変数をそれぞれ、x1とx2という変数に入れ直して、相関係数を見る。

x1 <- X$Naka
x2 <- X$Yasu
cor(x1, x2)

実行結果は以下のようになる。

> x1 <- X$Naka
> x2 <- X$Yasu
> cor(x1, x2)
[1] 0.8704575

相関係数の符号はプラスで、値は、0.8と1に近い値になっている。
ところで、相関係数の算出式はどようなものであったか？
今一度、思い出してみる。たしか、以下のような式だった。



要するに、分母が二つの変数の標準偏差に相当し、分子が共分散となっている。
二つの変数の値が全く同じ場合は、共分散の部分が分散となり、
結果として、相関係数が「1」になるのであった。忘れた人は、もう一度復習を。

さて、この式は、二つ変数の実際のバラツキの連動性を、
理屈上のバラツキの連動性で基準化しているということになる。
見方を変えると、この式は二つの変数を対等に見ているとも言える。

はて？「対等」ということはどのようなことか？
また、逆に「対等でない」ということはどのようなことか？

要するに、一方の変数のバラツキから、
もう一方の変数のバラツキを「説明する」という考え方もできる。



β の上に「＾（ハット）」が付いているが、これは推定値という意味。
数学の世界では慣例的に、ある数式によって導かれた値をこのように表す。

なお、相関係数の式では、二つの標準偏差の積を使うことで、
分散の値に揃えようとしていたのだけれど、ここでは一方のみで良いので、
分母の部分は分散に相当するようになっている。分子はそのまま。

さて、この式は、一体何を表しているか？

要するに、二つの共分散を一方の変数のみで基準化すると、
もう一方の「理屈上」の分散を推定できる、と考えているのである。

さらに突き詰めて考えてみると、一方の値が分かれば、
理屈上のもう一方の値の「増減率」を推定することができる。

さて、言葉で考えても解りにくいので、
とりあえず、散布図を描きながら考えてみることにする。

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

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

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

# まずは、普通に散布図を描く
plot(x1, x2, ylim=c(0,max(x2)), xaxt="n", yaxt="n", xlab="", ylab="")

# 図を重ね合わすためのオマジナイ。
par(new=TRUE)

# 推定され得る値を直線で表す
plot(x1.range, x1.estim, col="red", type="l", ylim=c(0,max(x2)), xlab="", ylab="")

# x軸とy軸の軸ラベルを付ける
title(xlab="x1", ylab="x2")

# 変数x2の平均（横線）と変数x1の平均（縦線）を描く
abline(h=mean(x2), v=mean(x1), lty=2, col="blue")

このように、相関係数とは異なった方法で二つの変数の関係を表すことができる。
今回の例の場合は、「それとなく」線と散布図がフィットしているように見えるが、
この赤い線によって推定されているのは、二つの変数の増減率の関係を表す「傾き」だけ。

値そのものを推定したいのであれば、二つの値の中心に原点を合わす必要がある。
要するに、上の図で言うと、二つの青い破線の交点に赤い線を通すようにすれば良い。
どのようにすれば良いかは、頭で考えても解るハズ。中学校までの知識で十分。

うん？という人のために補足説明。

推定の赤い線が、x1の平均値の時に、x2の平均値を取るようすれば良いので、
以下のような補正項を考えれば良いことになる。



x1の推定された値が、実際の平均よりもy軸方向にズレているのが問題なので、
そのズレの部分を計算しているだけ。大して難しい話では無い。
さて、この部分を先ほどの増減率の式と一緒にしてみる。



つまり、最終的にこのようになる。推定されるべきx2の値は、
補正項の部分とx1の推定された増減率の値を足したものからなる。

これを「単回帰」と呼ぶ。

実は、このブログでは、相関係数から単回帰についての説明を導いたが、
歴史的には、相関係数よりも単回帰の方が先に発見された。
ちなみに、発見した人物はナイチンゲールの親戚でもあるフランシス・ゴールトン。

少し話が逸れたが、この回帰直線と相関係数を同時に表示すると、
さらに、二つの変数間の関係がよく解るのである。
そのため、散布図を描いた際に、回帰直線を一緒に描くことがよくある。

そう言えば、先程の図は、まだ補正項の部分が修正されていなかった。
折角なので、補正した図を作成してみる。

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

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

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

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

# まずは、普通に散布図を描く
plot(x1, x2, ylim=c(0,max(x2)), xaxt="n", yaxt="n", xlab="", ylab="")

# 図を重ね合わすためのオマジナイ。
par(new=TRUE)

# 推定され得る値を直線で表す
plot(x1.range, x1.estim, col="red", type="l", ylim=c(0,max(x2)), xlab="", ylab="")

# x軸とy軸の軸ラベルを付ける
title(xlab="x1", ylab="x2")

# 変数x2の平均（横線）と変数x1の平均（縦線）を描く
abline(h=mean(x2), v=mean(x1), lty=2, col="blue")

今回は、仕組みを理解するために手計算で回帰直線を引いたが、
もちろん、Rには簡単に回帰分析を行うための関数がある。以下は、実行結果。

> lm(x2 ~ x1)
Call:
lm(formula = x2 ~ x1)
Coefficients:
(Intercept)           x1
   -69.6845       0.3637

Rでは、lm()関数というものを用いる。括弧の中身が解りにくいかもしれない。
これは、モデル式と呼ばれるもので、「〜」の左側に求めたい変数を置き、
右側に求めたい値を「説明するための変数」を置く。

統計学では、説明するための変数のことを「説明変数」と呼び、
その説明変数によって説明される変数のことを「目的変数」と呼ぶ。

ちなみに、

数学では、説明するための変数のことを「独立変数」と呼び、
その説明変数によって説明される変数のことを「従属変数」と呼ぶ。

文系出身者の多くには、二つの分野が同じように見えるかもしれないが、
前者は工学系の分野で、後者は理学系の分野。
論理の立て方が異なるので、与えられる用語は（記号も）微妙に異なる。

実際には同じなのだけれど、用語の使い方に惑わされないように。

さてさて。出力結果で見ないといけないのは、Coefficientsと書かれている部分。
(Intercept) というのは、日本語では「切片」と訳される。要するに、 の部分。
そして、x1と書かれている部分が「傾き」の係数部分になる。

lm()関数は、単回帰だけでなく、重回帰分析でも使用する関数なので、
今の間に覚えておいた方が良いかもしれない。


ふむ。確かに、推定するための線があった方が散布図は見易そうである。
ところで、この推定された値は、実際にはどの程度上手く当てはまっているのだろうか？
これを検討するには、推定された値と、実際の値とのズレを分析する必要があるのだが…。

今回も、話が少し長くなったので、この「ズレ」の検討については次回にしよう。