階層重回帰モデルの実例

広島二日目。前日に聞いた「坂の途中で乗り捨てられるスクーター」というサイジョーク（西条ジョーク）の坂はこの辺かな，とかいいながら広島大学に移動。

階層重回帰モデルについて，東大の深谷先生にご講演いただく。話を聞きながら，横でRでどのように実現させるのか，やってみました。以下，参考までにソース。使うデータは清水先生のサイトにあるHLM用サンプルデータをご利用ください。

<span class="synComment"># 　CSVファイルを読み込む</span>
sample_data <span class="synStatement">&lt;-</span> read.<a class="keyword" href="http://d.hatena.ne.jp/keyword/csv">csv</a><span class="synSpecial">(</span>file.choose<span class="synSpecial">(),</span>head=T<span class="synSpecial">,</span>na.strings=<span class="synConstant">&quot;.&quot;</span><span class="synSpecial">)</span>
<span class="synComment">#　読み込んだデータの確認</span>
summary<span class="synSpecial">(</span>sample_data<span class="synSpecial">)</span>
<span class="synComment">#　普通の線形回帰</span>
result.lm <span class="synStatement">&lt;-</span> lm<span class="synSpecial">(</span>idt~talk+per<span class="synSpecial">,</span>data=sample_data<span class="synSpecial">)</span>
summary<span class="synSpecial">(</span>result.lm<span class="synSpecial">)</span>
<span class="synComment">#</span>
<span class="synComment"># 　センタリング</span>
sample_data$talk_c <span class="synStatement">&lt;-</span> sample_data$talk - mean<span class="synSpecial">(</span>sample_data$talk<span class="synSpecial">,</span>na.rm=T<span class="synSpecial">)</span>
sample_data$per_c <span class="synStatement">&lt;-</span> sample_data$per - mean<span class="synSpecial">(</span>sample_data$per<span class="synSpecial">,</span>na.rm=T<span class="synSpecial">)</span>
<span class="synComment"># model A</span>
model.A <span class="synStatement">&lt;-</span> lm<span class="synSpecial">(</span>idt~talk_c+per_c<span class="synSpecial">,</span>data=sample_data<span class="synSpecial">)</span>
summary<span class="synSpecial">(</span>model.A<span class="synSpecial">)</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/AIC">AIC</a><span class="synSpecial">(</span>model.A<span class="synSpecial">)</span>
BIC<span class="synSpecial">(</span>model.A<span class="synSpecial">)</span>
<span class="synComment"># model B</span>
model.B <span class="synStatement">&lt;-</span> lm<span class="synSpecial">(</span>idt~talk_c+per_c+talk_c*per_c<span class="synSpecial">,</span>data=sample_data<span class="synSpecial">)</span>
summary<span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/AIC">AIC</a><span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>
BIC<span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>

以下，順に解説。まずここ。

<span class="synComment"># CSVファイルを読み込む</span>
sample_data <span class="synStatement">&lt;-</span> read.<a class="keyword" href="http://d.hatena.ne.jp/keyword/csv">csv</a><span class="synSpecial">(</span>file.choose<span class="synSpecial">(),</span>head=T<span class="synSpecial">,</span>na.strings=<span class="synConstant">&quot;.&quot;</span><span class="synSpecial">)</span>
<span class="synComment">#読み込んだデータの確認</span>
summary<span class="synSpecial">(</span>sample_data<span class="synSpecial">)</span>

read.csvでcsvファイルを読み込みます。環境によって，どこにサンプルファイルをおいているか分からないので，file.choose()関数にしました。この行を実行してもらうと，ファイル選択ダイアログが開くので，サンプルファイルを指定すると読み込まれます。

モデルは基本的に，idtが従属変数，talk,perが独立変数。Rの書き方に従って，~(チルダ）の左側に従属変数，右側に独立変数を書く。独立変数が複数ある場合は，+でつなげていくという感じ。この変数間関係を回帰分析＝線形モデル＝linear modelにいれると

result.lm <span class="synStatement">&lt;-</span> lm<span class="synSpecial">(</span>idt~talk+per<span class="synSpecial">,</span>data=sample_data<span class="synSpecial">)</span>
summary<span class="synSpecial">(</span>result.lm<span class="synSpecial">)</span>

となります。lm関数の結果を，result.lmという入れ物の中にしまいなさい（左向きの矢印を「小なりハイフン（< -)」で表している），と書いて，その入れ物の中身を要約して示せ(summary)という書き方。summaryなしで直接書き出してもいいのだけど，lmの場合は要約した方が情報が増えてグッド。

さて，センタリングがHRMのミソ。独立変数をセンタリングするために，平均値を引きます。それがここ。

<span class="synComment"># 　センタリング</span>
sample_data$talk_c <span class="synStatement">&lt;-</span> sample_data$talk - mean<span class="synSpecial">(</span>sample_data$talk<span class="synSpecial">,</span>na.rm=T<span class="synSpecial">)</span>
sample_data$per_c <span class="synStatement">&lt;-</span> sample_data$per - mean<span class="synSpecial">(</span>sample_data$per<span class="synSpecial">,</span>na.rm=T<span class="synSpecial">)</span>

平均値を出すのはmean関数。sample_dataの中のtalk変数，をsample$talkで指定しています。いちいちsample_dataって書くのが面倒な人は，attach(sample_data)と書いておけば，以下の実行は全部sample_dataというオブジェクトの中の変数かな，とチェックしてくれます。
mean関数の後ろのオプション(na.rm=T)は欠損値(na)を取り外せ(リムーブ，rm)という意味で，欠損値の取り外し(na.rm)にチェック！（=TRUE)という意味。これを入れていないと，欠損値の含まれる変数の平均値は算出されません。
その計算結果を，sample_dataの中に新しくtalk_cとして入れなさい（< -）という意味です。

で，こうしてセンタリングが終わった変数で，改めて普通の線形回帰をしてみる。

<span class="synComment"># model A</span>
model.A <span class="synStatement">&lt;-</span> lm<span class="synSpecial">(</span>idt~talk_c+per_c<span class="synSpecial">,</span>data=sample_data<span class="synSpecial">)</span>
summary<span class="synSpecial">(</span>model.A<span class="synSpecial">)</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/AIC">AIC</a><span class="synSpecial">(</span>model.A<span class="synSpecial">)</span>
BIC<span class="synSpecial">(</span>model.A<span class="synSpecial">)</span>

切片は当然違うけど，傾きの係数は同じのままね。
ついでに，情報量を出す関数AIC,BICを入れてみました。後ほどの比較のためです。

最後に，センタリングをふまえての交互作用項（調整項）を入れる。交互作用はかけ算（アスタリスク，*)で表される。

model.B <span class="synStatement">&lt;-</span> lm<span class="synSpecial">(</span>idt~talk_c+per_c+talk_c*per_c<span class="synSpecial">,</span>data=sample_data<span class="synSpecial">)</span>
summary<span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>
<a class="keyword" href="http://d.hatena.ne.jp/keyword/AIC">AIC</a><span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>
BIC<span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>

とまぁ，こうすることで交互作用項がはいるという寸法です

下位検定については，Preacherのサイトを使うといいそうです。が，Preacherのサイトを見るとMBESSというパッケージを作っているみたいなので，この辺の使い方が分かったらまた書きます。

ひとまず。

追記）プリちゃんのサイトでひつような，係数の分散共分散行列は，次のようにすれば出せます。

vcov<span class="synSpecial">(</span>model.B<span class="synSpecial">)</span>

信頼帯とかも書けるので，いいサイトだとは思うのですが，Rのパッケージ上でやれるようにしてほしいもんです。
intr.plot(),intr.plot.2d()とかの関数がそれかな，と思ったけどエラーが出てうまくうごかなーい。