ラベル 単回帰 の投稿を表示しています。 すべての投稿を表示
ラベル 単回帰 の投稿を表示しています。 すべての投稿を表示

2012/10/04

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

変数の関係を観察するための基本的方法に「相関」と「回帰」の考え方がある。
相関は、変数間の結びつきの「連動性」を表す方法であり、
回帰は、変数間の結びつきの「関係性」を表す方法であった。

ここで「うん...?」となった人は、もう一度、相関の話回帰の話を復習。

相関の考え方の場合では、複数の変数の場合相関係数行列を用いたが、
回帰の考え方の場合には、どのように表すのか?というのが、今回のテーマである。

ところで、一般的な解説では「説明する」という言葉に重きを置いて、
回帰分析を解説することが多いが、回帰分析は変数間の関係を観察する手法である。
「説明変数」と「目的変数」という言葉に惑わされてはいけない。

まずは、2変数の場合における回帰、すなわち、単回帰の話を思い出してみる。
たしか、以下のような式によって、表されたのであった。



この式では、一方の変数のみで、もう一方の変数を変数を表そうとしていることが解る。
一般的な言い方をすると、一方の変数でもう一方の変数の「説明」を試みている。
ここで、説明される側の変数を目的変数説明する側の変数を説明変数と呼んだ。

この式において左辺の目的変数には「^(ハット)」が付いている。
したがって、説明変数によって導かれているのは推定値だと解る。
ここで、実際の値を表現するのであれば、以下のようになる。



左辺の目的変数の「^(ハット)」が取れて、右辺の最後に「」が付くだけ。
最後に付け足した「」は、説明しきれない「残差(誤差)を表すのであった。
説明変数によって、説明しきれなかった部分を最後に足してやれば元の値が解るハズ。

おそらく、ここまでの話は、理解できている...と信じたい。

さて、変数がさらに多くなるとどうなるのか?
p個の変数からなる多次元データがあったとして、
一つの変数を目的変数としたとき、次のように表すことができる。



p-1となっているのは、元の変数の内の1つが目的変数とな っているため。
p個の変数のうち、説明変数を抜いた「p-1」個の変数が説明変数となる。
この式では、目的変数推定値であるが、実際の値は次のように表すことができる。



これが、「重回帰分析」と呼ばれるもののモデルである。
添字が並んでいるので、混乱するかもしれないが、
要するに、変数が増えただけ。モデルとしては、基本的に単回帰と同じ。

これを実際に計算するとなると、少々面倒な手続きが必要となる。

まずは、目的変数を推定する式をと眺めてみる。気づくことはないか?
このままだと、少々、解りにくいかもしれない。
この式はベクトルの式で表されているので行列の式で表してみる。



行列の「掛け算」を理解していれば、難しくは無いハズ。
ここで知りたいのは、「β」のベクトルの部分。回帰係数のベクトルである。
これで、意味が解った人はスゴイ。解らない人のために式を変形してみる

たしか、行列には割り算が存在しない代わりに、
スカラー逆数に相当する、逆行列というのがあった。
まずは、これを使って両辺を基準化してみる。すると以下のようになる。



逆行列の性質から、元の行列に逆行列を掛けると、単位行列が出てくるので、
要するに、上の式は、実際には以下のようになっている。



単位行列に別の行列を掛けると、掛けた行列が出てくる。
つまり、何も起きないのであった。したがって、次のようになる。



最後に、左辺と右辺を入れ替えると次のようになる。



さて、ここまで変形してみて意味が解ればと合格。
要するに、連立方程式の形にもってきただけのこと。
逆行列を使えば連立方程式が解ける

もう少し、一般的な式に書きなおしてみると、以下のようになる。



本当にこの式で大丈夫だろうか?
逆行列は正方行列にしか使えないので、擬逆行列にしないといけない。
つまり、以下のようになる。擬逆行列の話を思い出してみる。



あるいは、特異値分解による逆行列の近似によって、以下のようにも表すことができる。
もちろん、計算結果は上の式と同じ。詳しくは、擬逆行列の話でも述べている。



以上のようにして、複数の変数における回帰係数を求めることができる。
すでに、ここまでの話を、一つずつ積み重ねてきた人は、
Rで重回帰分析における回帰係数を計算することができるハズ。

次の話では、実際のデータを用いてもう少し、
重回帰分析意味モデルの見方について整理したい。

2012/09/26

文系のための「多次元データの要約」(2)

分散共分散行列相関係数行列を用いることで、様々な定量的手法が可能となる。
実際の手法については「中級レベル」で解説するとして、要するに、
これら行列の意味を理解することは極めて重要である。

そもそも、この二つの行列は、どのようなものであったか?

たしか、分散共分散行列は、各変数の分散が対角成分に並び、
それ以外の成分には「共分散が入っているような「対称行列」であった。

一方、相関係数行列は、対角成分に「1」が並び、
それ以外の部分には「相関係数が入っているような行列であった。
相関係数行列は、変数間の関係を一目で確認できるので、色々と重宝するのだった。

ふむふむ。そのような話を確かにした。何か「問題」でもあるのか?

実は、2変数の相関係数の話では気づかなかったのであるが、
多変数の相関係数行列にしてみると、ある種の疑惑が生まれる。

本当に、二つの変数の関係「のみ」を表しているのか?
実は、他の変数の影響を受けているのではないか?
という疑惑である。

例えば、3つの変数、「」、「」、「」の3つの変数が存在した場合、
この3つの変数の相関係数行列は以下のようになるはず。



ふむ。何も不審な点は見当たらないように見えるが、
ひょっとすると、「」と「」の関係の中には、
」の影響が「潜在的」に存在するかもしれない。

2変数の間の相関係数の話では、意識する必要は無いが、
たしかに、多変数となる相関係数行列では、そういった懸念が生じる。
では、そのような別の変数の影響を取り除くにはどうすれば良いか?

この問題が、本日のメインテーマ。

とりあえず、いつものように、いつもの如く、チュートリアルデータを準備する。
今回のデータは、総務省統計局のデータ(e-stat)のデータを用いて、
チュートリアル用の「擬似」データを作成したもの。
今回のデータは、居住世帯の有無一ヶ月辺りの平均家賃の関係を見るためのデータ。
家賃のデータは、19階級に区分し家賃を集計したものであったので、
各階級の平均を一ヶ月の家賃とし、世帯数を掛けたもので都道府県別の平均家賃を算出。

かなり「大雑把な加工」をしているので、正式な分析には用いてはならない
あくまで、手法を学ぶための「擬似データ」であるので、
ここから実際の解釈を導いても、「深い意味」はない。

一応、注意事項を整理したところで、いつもどおりにデータを読み込む。

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

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

今回のチュートリアルのために用意したデータは以下の通り。

OCU,UNO,PRI
Hokkaido,2340300,390100,40687
Aomori-ken,493500,87400,111566
Iwate-ken,470700,78800,17002
Miyagi-ken,869700,144200,16084
Akita-ken,380300,57100,41254
Yamagata-ken,383000,49700,10523
Fukushima-ken,699700,108600,11879
Ibaraki-ken,1036200,187600,25207
Tochigi-ken,708700,131300,36003
Gumma-ken,725300,130400,27494
Saitama-ken,2688000,341100,22045
Chiba-ken,2344500,373100,94198
Tokyo-to,5939900,840500,60863
Kanagawa-ken,3612200,455600,246996
Niigata-ken,810700,119000,113336
Toyama-ken,368800,55500,24592
Ishikawa-ken,421600,76400,10485
Fukui-ken,259700,49000,15203
Yamanashi-ken,314600,83600,6005
Nagano-ken,758300,188000,10698
Gifu-ken,712600,123100,25608
Shizuoka-ken,1359400,238500,22521
Aichi-ken,2764400,368400,54336
Mie-ken,680900,110100,122595
Shiga-ken,491300,76300,17786
Kyoto-fu,1086800,183400,14849
Osaka-fu,3685100,660900,39764
Hyogo-ken,2169400,351300,150179
Nara-ken,502500,90100,68843
Wakayama-ken,382100,85700,12267
Tottori-ken,208600,38600,10868
Shimane-ken,249900,45900,6686
Okayama-ken,734700,131900,7459
Hiroshima-ken,1147600,208700,27600
Yamaguchi-ken,584100,107600,48448
Tokushima-ken,297000,58600,21579
Kagawa-ken,372700,73700,10305
Ehime-ken,574000,107100,12877
Kochi-ken,312800,64900,22953
Fukuoka-ken,2034000,340800,11205
Saga-ken,286100,36800,103313
Nagasaki-ken,539200,91900,9725
Kumamoto-ken,663800,105700,20633
Oita-ken,467200,79300,29546
Miyazaki-ken,443800,65900,22566
Kagoshima-ken,718200,133100,17694
Okinawa-ken,504400,62100,29919

このデータにおいて、
  • OCU居住世帯
  • UNO非居住世帯
  • PRI一ヶ月の平均家賃
まずは、分散共分散行列相関係数行列を出してみる。

cov(X)    # 分散共分散行列
cor(X)    # 相関係数行列

以下が、その結果。

> cov(X)    # 分散共分散行列
             OCU          UNO         PRI
OCU 1.297147e+12 185407240994 21146986138
UNO 1.854072e+11  27363809917  2614407571
PRI 2.114699e+10   2614407571  2151559501

> cor(X)    # 相関係数行列
          OCU       UNO       PRI
OCU 1.0000000 0.9841103 0.4002922
UNO 0.9841103 1.0000000 0.3407284
PRI 0.4002922 0.3407284 1.0000000

さて、相関係数行列を見てみると、
居住世帯数OCU)と非居住世帯数UNO)の相関が非常に高く
平均家賃PRI)と他の変数の関係は強く無さそうに見える。

相関の有無は、ひとまずは置いておいて、
平均家賃に注目してやると、居住世帯数が多いと賃料が高く
同様に、非居住世帯数が多くて「も」平均賃料が高くなる傾向がある。

うん?非居住世帯が多くなるほど、月の平均賃料が高くなるというのは、
直感的に合わないような気がする。
居住・非居住の世帯数の関係が影響しているかも。

問題は、2変数間の関係に、別の変数の影響が存在している可能性があること。
では、この「潜在的」な影響は、どのようにしたら見えてくるのか?

ここで、相関係数とは別の方法で2変数の関係を考えてみる。
別の方法とは、すなわち、単回帰と呼ばれる方法である。
何のことか解らない人は、もう一度、やり直し

たしか、単回帰は、2つの変数の関係について、
一方の変数で、もう一方の変数を「説明」するという考え方だった。
この方法を用いることで、他の変数の潜在的な影響を知ることができるかもしれない。

OCUで「説明できない部分」のPRIの値と、
OCUで「説明できない部分」のUNOの値との、
関係を比較」すれば良いのではないか?とにかくやってみる。

手を動かした方が解りやすいので、まずは、普通に回帰分析を行なってみる。
Rでは、単回帰を行うためのlm()関数があるので、今回は、この関数を使ってみる。

# 平均家賃を固定し、単回帰を行う
X.lm.3_1 <- lm(X[,3]~X[,1])    # PRIをOCUで説明する
X.lm.2_1 <- lm(X[,2]~X[,1])    # UNOをOCUで説明する

本当は、PRIとUNOを比較したいが、OCUの影響が気になる。
したがって、PRIUNOをそれぞれ「y軸」に固定し、OCUからの説明を考える。
# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))

# PRIをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,3], xlab="OCU", ylab="PRI", ylim=c(0, max(X[,3])))

# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.3_1$fitted.value, type="l", ylim=c(0, max(X[,3])), xaxt="n", yaxt="n", xlab="", ylab="")

# UNOをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,2], xlab="OCU", ylab="UNO", ylim=c(0, max(X[,2])))

# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.2_1$fitted.value, type="l", ylim=c(0, max(X[,2])), xaxt="n", yaxt="n", xlab="", ylab="")

次に、各々の「残差の状況」を確認してみる。
# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))

# PRIとOCUの残差の検討(標準化残差)
plot(X[,1], scale(X.lm.3_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="PRI")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")

# UNOとOCUの残差の検討(標準化残差)
plot(X[,1], scale(X.lm.2_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="UNO")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")

この図では、残差の状況を見やすくするために残差を標準化している。
これは「標準化残差」と呼ばれる。「±2」の線は、おおよその「外れ値」の基準線。
詳しい話は、正規分布の話で述べる。

さてさて、たしか、理屈上は、説明できない部分が「残差」と呼ばれるものだった。
つまり、OCUの影響を取り除くには、「残差同士の相関」を見れば良い。

# 残差同士の相関係数
cor(X.lm.3_1$residuals,X.lm.2_1$residuals)

以下は実行結果。

> # 残差同士の相関係数
> cor(X.lm.3_1$residuals,X.lm.2_1$residuals)
[1] -0.3269775

なるほど。OCUの影響を抜いた状況相関係数を見てみると、
月の平均賃料非居住世帯数の関係は、「負の関係」となった。
つまり、非世帯数が増加するほど、平均賃料は低くなる傾向がある。

この結果は、影響を抜かない相関係数よりも直感に即しているように思える。
すなわち、人が住んでいない住宅が多いということは、
需要に対する過剰供給が背後に存在し得るので賃料が低くなる

試しに、残差同士の関係を散布図で可視化すると以下のようになる。
# とりあえず、変数を定義し直す
pri.res <- X.lm.3_1$residuals
uno.res <- X.lm.2_1$residuals

# 残差同士の単回帰(偏回帰)
res.lm <- lm(uno.res ~ pri.res)

# 残差同士の関係を散布図で表す
plot(pri.res, uno.res, ylim=c(min(uno.res),max(uno.res)), xlab="", ylab="")

# 一応、回帰直線も追加
par(new=TRUE)
plot(pri.res, res.lm$fitted.value, type="l", xlab="", ylab="",  xaxt="n", yaxt="n",  ylim=c(min(uno.res),max(uno.res)), col="red")

# 最後に図を整える
pcor <-  round(cor(pri.res, uno.res),3)
plm <- round(res.lm$coefficient[2],3)
leg1 <- paste("Partial Correlation", pcor, sep=": ")
leg2 <- paste("Partial Regression Coefficient", plm, sep=": ")
legend("topright", legend=c(leg1, leg2))

title(xlab="Residuals for building rental price", ylab="Residuals for unoccupied buildings")
abline(h=0, v=0)

さて、このように、ある2つの変数の関係を観察するために、
別の変数の影響を取り除いた相関係数のことを、
偏相関係数(Partial Correlation)」と呼ぶ。

また、ここでは、残差同士の単回帰を行い、その回帰直線も描いているが、
残差同士の単回帰の係数のことを「偏回帰係数(Partial Regression Coefficients)」と呼ぶ。
ちなみに、残差同士の関係なので切片は「0」となる。

さて、今回は、単回帰の残差の関係から偏相関を導いたが、
一般的に、偏相関係数は元の相関係数から、以下のように計算できる。



なお、は、の影響を抜いた と との偏相関を意味し、
 はそれぞれ、元の相関係数を意味している。
」は、ギリシア文字で「ロー」と発音する。「ピー」では無い。

さて、要するに、知りたい2つの変数の相関係数から、
潜在的に影響しているかもしれない別の1つの変数の影響を取り除き
結果の値が、「±1」の間に収まるように分母の部分で基準化している。

今回のように、3変数の場合は、直感的に理解しやすいのであるが、
4変数以上になると計算が非常に厄介になってくる。行列を使えば、
式自体はコンパクトになるが...余裕の無い人は以下の式の話は飛ばしても良い。

まず、元の相関係数行列における、変数「i」と「jとの相関係数を「」 とし、
その相関係数行列の逆行列におけるi」と「j」の関係を 「」する。
ちなみに、逆行列の添字を右下ではなく「右上」に置く。ちょっと、注意。

このとき、4変数以上の偏相関係数行列は、次のように求める。



逆行列の対応する要素を、2つの対角要素の平方根で割って基準化し、
さらに、「符号を反転」させている逆行列を使った数学の魔法!!
この方法で、偏相関係数を計算することができるのである。

まぁ、実際には、手計算で偏相関係数を求める人は少ないであろう。

Rでは、この偏相関係数を簡単に計算できるため、通常はこれを用いるのであるが、
残念ながら、偏相関行列を計算してくれるパッケージは最初からは入っていないので、
「corpcor」パッケージをインストールする。

# corpcor パッケージのインストール
install.packages("corpcor", dep=TRUE)

上手くインストールできたら、実際に計算してみる。

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

# 偏相関係数行列を計算する。
X.pcor <- cor2pcor(cor(X))

# 変数名が付かないので付け直す。
rownames(X.pcor) <- rownames(cor(X))
colnames(X.pcor) <- colnames(cor(X))

以下が、この実行結果。

> X.pcor
          OCU        UNO        PRI
OCU 1.0000000  0.9839439  0.3892439
UNO 0.9839439  1.0000000 -0.3269775
PRI 0.3892439 -0.3269775  1.0000000

PRIUNOの関係を見てみると、確かに、残差同士の相関係数と一致している。

相関係数行列は、変数間の関係を観察する上で非常に重要であるが、
さらに、細かい状況を把握するためには、
偏相関係数行列で見た方が、状況が良く解る場合もある。

2012/09/19

文系のための「二変数の関係」(3)

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

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

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

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

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

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

# 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」で割ると?そうそう、分散になる。
要するに、分散の式の分子の部分に当たるのであった。

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

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