今度は午後の部，階層線形モデルの話。
ネストされたデータは全部HLMの土俵。反復測定も、個人と集団も。後者の方が得意らしいけど。 ネストされたデータで，HLMを使わないと、サンプルの独立性の検定に違反する。また、平均値が集団の性質を反映していない。後者はカップルデータのようなときに顕著。グループ内の類似性を評価し、それに合わせたモデリングをする。それがHLM，というお話。
ところで，階層回帰分析と階層線形モデリングは，同じ階層という言葉を使っているけど，意味が全然違う。
前者は手続きが順番に行われる，という意味で階層的であり，まず要因A，要因Bを入れ，次に交互作用ABを入れる，という順番でステップを踏む。あくまでも順番であって，あんまり層を積み重ねるという意味じゃないと思う。実際英語としては，Multiple Linear Regression，つまりただ「重回帰分析」と表されるね。重回帰分析は二つ以上の要因がある，という広い意味で使われる。交互作用項を特に【調整変数】と呼んで，センタリング等の適切な処置を経て投入する，というところが特徴。
個人的には「順番にやる回帰分析」「一歩ずつ重回帰」，といった名称にすればよかったんじゃないかと思う。
（実は内心，両者の区別が分からないやつが間違って使い始めたのが定着したのではないかと疑っている）
さて，HLMは本気で階層的。データがレベルを持っている。例えば個人のデータは集団の性質、個人の性質、誤差からなる。場合によっては集団と個人の交互作用も考えることがある。グループ内の類似性が高いことは、集団レベルの情報を多くもっているということ
。
ここでMLRと混乱させられるもう一つの秘密が。それは，HLMでもステップ1，ステップ2という用語があり得るのです。ほんとはレベル1，レベル2というのが正しい。レベル1は個人のデータ。個人レベルの情報。Within(群内)とも呼ぶ。レベル2は集団のデータ。集団レベルの情報で，Betweenとも。変数がどっちレベルで集まっているか，がしっかり把握できていないと混乱するぜ。HLMはレベル2の変数がレベル1の係数を予測する（回帰のパスが下位レベルの係数にささる）こともでき，それを「交互作用」と呼んじゃうからさらにMLRとの混同がおきやすいよね。要注意。
さて，実際にデータを触りながらやった方が分かりやすいかと思うので，Rソースを示し，一歩ずつ解説をしていきます。
サンプルデータは上の記事と同じで，使うデータは清水先生のサイトにあるHLM用サンプルデータをご利用ください。分析資料も同じサイトにあるよ。
RでHLMをするには，lme4というパッケージが必要です。
ソースはこんな感じ。
# #CSVファイルを読み込む sample_data < -readcsv(file.choose(),head=T,na.strings=".") # 読み込んだデータの確認 summary(sample_data) # いちいち書くのがめんどうなので，以下データセット名は自動補完するように仕掛け。 attach(sample_data) # ライブラリの読み込み require(lme4) # センタリング # グループセンタリング sample_data$talk_c <- - talk - ave(talk,Group,na.rm=T) # グランドセンタリング sample_data$per_g <- - per -mean(per,na.rm=T) # 集団平均値の挿入 sample_data$talk_g_m <- - ave(talk,Group,FUN= function(x)mean(x,na.rm=T))-mean(talk,na.rm=T) # Model.1 <- lmer(idt~talk_c+per_g+( 1|Group),sample_data) Model.2 < - lmer(idt~talk_c+talk_g_m+per_g+(1|Group),sample_data) Model.3 < - lmer(idt~talk_c+talk_g_m+per_g+talk_c*per_g+(talk_c|Group),sample_data) # 適合度指標 AIC(Model.1,Model.2,Model.3) BIC(Model.1,Model.2,Model.3)
さて解説。最初の数行は割愛させてください。センタリングの話から。
HLMでもセンタリングをやります。個人レベルの変数に行うセンタリングは，グループごとの平均値で中心化するセンタリングで「グループセンタリング」といいます。集団レベルの変数に行うセンタリングは，全体の平均値を使うやつで「グランドセンタリング」といいます。
それがここ。
# センタリング # グループセンタリング sample_data$talk_c <- talk - ave( talk, Group, na.rm=T) #グランドセンタリング sample_data$per_g <- per -mean( per, na.rm=T)
mean関数は平均値を出す関数。na.rmは欠損値を外せというオプション。もう一つ，グループセンタリングはave関数というのがあって，変数Group毎に平均値を出すという便利な関数。
次の行は，集団平均値をデータとして使う場合があるので，そのための細工。
# 集団平均値の挿 sample_data$talk_g_m<- ave(talk,Group,FUN=function(x)mean(x,na.rm=T))-mean(talk,na.rm=T)
グループごとに平均を出して，さらにそれをセンタリングしたものを入れています。ave関数がちょっとおかしな形になっているけど，これはデフォルトだとna.rmが入らないから。
ave関数は本来，ave(x,id,FUN)と書いて，データに，群ごとに，ある処理FUNをする，というもの。FUNはデフォルトでmeanなんだけど，このmeanのデフォルトがna.rm=Fなので，それを教えてあげないといけない。つまり，FUNはfunction(x)を使うよ，そしてそれはmean(x,na.rm=T)だよ，という二度手間構造(もっといいやり方があったら教えてエロい人)。
ともかく，これで準備オーケー。
まずは群ごとに，idtを従属変数，talk,perを独立変数とした回帰分析をしてみる（清水先生の分析モデル1）。
Model.1<- lmer(idt~talk_c+per_g+(1|Group),sample_data)
関数は，パッケージmlmRevに入っているlmerで，書き方はlmと同じ，チルダの左に従属変数，右に独立変数。独立変数はセンタリングしたtalk_cとper_gを入れる。階層性を表す変数Groupは係数1で関わってきますよ，というのが(1|Group)の意味。
次に，集団平均値を入れたモデル（清水先生資料の分析モデル2）
これはさっきの変数を使うだけだから簡単。
Model.2<- lmer(idt~talk_c+talk_g_m+per_g+(1|Group),sample_data)
Model.1に比べて，talk_g_mという独立変数が増えただけです。
最後に，ランダム係数モデル（清水先生の分析モデル3）。ここでのランダムは乱数という意味ではなく，無作為でもなく，確率変数(random variables)という意味でのランダムね。
Model.3<- lmer(idt~talk_c+talk_g_m+per_g+talk_c*per_g+(talk_c|Group),sample_data)
ここでは，Groupごとにtalk_cの係数が変わってくるよ，ということと，交互作用項talk_c*per_gが増えています。
いずれもHADというエクセルマクロで下準備し，SPSSで処理をするという（資料上の）一連の流れがこれだけで再現できるので，どこにどのような数値が表れているかと確認しながらやってみてください。
級内相関の出し方なんかについては，パッケージを使えば出来るんだけど，それまたちょっと調べてから後日ブログにアップする予定です。
ひとまず。
追記（2012/02/14)
清水先生のスライドだと，まず級内相関をだして階層性を入れる意義を検証しようね，という話がありました。
この級内相関，パッケージで簡単にでるかなとおもったけど，色々調べても手計算している例が多い。
清水氏のHADみたいに，まず各変数において級内相関を（アルファ係数までも！）出すようなかんすうがあればいいんですけどね。
ひとまず，手計算のやり方を書いておきます。
Model.0<- lmer(idt~(1|Group),sample_data) Model.0
とまぁ，このように，まずグループレベルの変数しか入ってない回帰モデルを作る
Linear mixed model fit by REML Formula: idt ~ (1 | Group) Data: sample_data AIC BIC logLik deviance REMLdev 813.5 824.5 -403.7 804.1 807.5 Random effects: Groups Name Variance Std.Dev. Group (Intercept) 0.35004 0.59164 Residual 0.63566 0.79729 Number of obs: 297, groups: Group, 100 Fixed effects: Estimate Std. Error t value (Intercept) 3.43293 0.07514 45.69
と，その結果として集団レベル変数の分散と残差の分散がでるので，集団レベルの分散/(集団レベルの分散＋残差の分散)
を手計算してやる
0.35004/(0.35004+0.63566) [1] 0.3551182
いいってことになります。
以上。