階層重回帰モデルの実例

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

階層重回帰モデルについて,東大の深谷先生にご講演いただく。話を聞きながら,横で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.csvcsvファイルを読み込みます。環境によって,どこにサンプルファイルをおいているか分からないので,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()とかの関数がそれかな,と思ったけどエラーが出てうまくうごかなーい。

日記
父の日

昨日は父の日でしたね。 これまでやったこともやられたこともなかったので、意識したこともなかったが…… …

日記
蕁麻疹て結局なんも分からんのやないか

皮膚科に行きました。 朝起きたら、例によって下半身があちこち噛まれたような跡になってて痒い。ただ、以 …

日記
布団の外科手術

夜中に足が痒くて目が覚めるんですよ。 見たら、ものすごく「何か」に咬まれてる。蚊にかまれたみたいな跡 …