階層重回帰モデルの実例

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

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

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

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

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

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

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

平均値を出すのは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として入れなさい(< -)という意味です。

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

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

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

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

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

ひとまず。

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

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

コメントは受け付けていません。