1からみなおす線形モデル (1) - 一般線形モデル

はじめに

階層ベイズモデルとMCMC法の理解をしようと思って，Tokyo.Rまわりの資料をいくつかあさってみていたんですが（Rで階層ベイズモデルとかマルコフ連鎖モンテカルロ法入門-１とかマルコフ連鎖モンテカルロ法入門-2とか），@tsuna_1217に下の本を教えてもらったので，この本でも読みながらしばらく復習がてら線形モデルについてまとめていくことにしました．

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (18件) を見る

ということで，まだAmazon注文中なわけですが，最初の線形モデルについては本なしでもさくっとまとめておくことにします．

線形モデル

ここでいう（狭義の）線形モデルってのは，いわゆる回帰分析モデルのことですね．複数のxからyを線形に予測するってやつです．たとえば，身長と足のサイズから体重を予測する，みたいなものですね．

$y=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}$

ほんとうは，上記のモデルが成り立つためには各パラメタが正規分布に従うとか，各サンプルの値が独立であるとか，いくつかの仮定があるわけですが，そのあたりはここで深くは踏み込みません．これをRでやるのは，別に今更感もありますが，lm()を使ってあげればよいわけです．サンプルとしてはじめから入っているwomenデータを読み込んであげて，実際に身長から体重を予測する計算をしてみましょう．

> data(women)
> women.lm <- lm(weight~height, data=women)
> summary(women.lm)

Call:
lm(formula = weight ~ height, data = women)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared: 0.991,	Adjusted R-squared: 0.9903 
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

みてのとおり，調整済み決定係数が0.99で，heightの有意確率も1.091e-14 なので，十分に有意という結果です．基本的に身長と体重は比例するので，当たり前の結果ですね．

一般線形モデル

さて，上でみた数式ですが，これを行列を使ってもう少し抽象的な形に直してみると，以下のようにかけます．これが一般線形モデルです．

$Y=BX+E$

実は，分散分析や共分散分析などは，上記の回帰分析のモデルと同じ行列モデルに落とし込むことができるわけです．というのは，2群の分散分析というのは，先ほどのxに，1か0かの2値変数を当てはめたときになるわけです．さらに3群の分散分析だと，群1か否かの2値変数，群2か否かの2値変数の，2つの独立変数を当てはめた形になります．つまり，上記の行列式の形にすることで，回帰分析と分散分析を同じモデルで扱うことが可能になるわけです．

実際に，先ほどのデータを3群にわけてみて，3群変数を作って分析を行ってみましょう．身長が高い5人をグループ3，それ以外をグループ1と2に分けてみました．このgroupはカテゴリカル変数なので，factor()でカテゴリカル扱いにして独立変数に投入します．Rのlm()では，カテゴリカル変数の場合最初の群を基準として，グループ2か否か，グループ3か否かという独立変数を自動で作成して分析を行います．

結果をみればわかるように，グループ3でのみ係数が有意になっています．つまり，グループ3である場合は，それ以外に比べてweightが25.4高くなるわけです．これで分散分析モデルが回帰分析と同じ形で扱えることがわかるかと思います．

ちなみに，なんでRで回帰分析モデルを行った場合に，上記のような変数の自動変換が行われるかについては，こちらの14-17ページあたりに詳しい解説がありますので，ご一読ください．簡単にいうと，そうしないとパラメタを推定できないからですが...　あとここでは分散分析についてのみ書きましたが，もちろん共分散分析も同じようにlm()で実施することができます．

> women
   height weight
1      58    115
2      59    117
3      60    120
4      61    123
5      62    126
6      63    129
7      64    132
8      65    135
9      66    139
10     67    142
11     68    146
12     69    150
13     70    154
14     71    159
15     72    164
> women$group = ifelse(women$height%%2==1, 1, 2)
> women$group = ifelse(women$height>67, 3, women$group)
> women
   height weight group
1      58    115     2
2      59    117     1
3      60    120     2
4      61    123     1
5      62    126     2
6      63    129     1
7      64    132     2
8      65    135     1
9      66    139     2
10     67    142     1
11     68    146     3
12     69    150     3
13     70    154     3
14     71    159     3
15     72    164     3
> women.lm <- lm(weight~factor(group), data=women)
> summary(women.lm)

Call:
lm(formula = weight ~ factor(group), data = women)

Residuals:
   Min     1Q Median     3Q    Max 
 -12.2   -6.3   -0.4    5.7   12.8 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     129.200      3.977  32.483 4.58e-13 ***
factor(group)2   -2.800      5.625  -0.498 0.627634    
factor(group)3   25.400      5.625   4.516 0.000707 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 8.894 on 12 degrees of freedom
Multiple R-squared: 0.7177,	Adjusted R-squared: 0.6707 
F-statistic: 15.26 on 2 and 12 DF,  p-value: 0.0005056

ということで，今回は単純な回帰分析（=線形モデル）と，その拡張である一般線形モデルについてみてきました．最初の方で軽く触れましたが，一般線形モデルは基本的にすべてのパラメタが正規分布に従うことを仮定しています．そのため，従属変数が正規分布に従わない場合（たとえば従属変数が生きている/死んでいるの2値変数である場合など）には，当然使うべきではないわけです（むりやり当てはめることはできますが，当然意味のある結果とはいえないわけです）．

そういったデータを扱うためには，正規分布以外の分布を仮定したモデルを考えて行く必要があります．それが一般化線形モデルと呼ばれているものです（似た名前でややこしいですが，英語でいうと一般線形モデルがgeneral linear modelで，一般化線形モデルがgeneralized linear modelですね）．ということで，次回は一般化線形モデルについてみていきます．