混合一般化線形モデルの解析テスト¶
前置き¶
同じ形式のデータがいくつかのグループに分かれている場合、全体を一つの回帰モデルで分析するのでは不十分で、グループごとの特性を踏まえるほうが妥当な場合がある。
たとえば、
- プレテストと本テストの関係データで、学校間の差異を調べてみる必要がある場合 (教育学分野)
- 生態学の調査で、複数の調査地点のデータがあり、地点別の特性を踏まえるべき場合 (生態学)
などである。
方法としては同じであるが、教育学分野では「階層線形モデル」(Hierarchical Linear Model: HLM)、生態学分野では「一般化混合モデル」(Generalized Linear Mixed Model: GLMM)と呼ばれる。
データは、「回帰分析入門」（豊田秀樹、東京図書）９章より
この豊田本では、mlmRevというライブラリを用いている。mlmRevは、評価(review)用データの集合体のようで、統計ライブラリとしてはlme4が用いられているので、実践的には、lme4を使えばよいだろう。
https://rpubs.com/hoxo_m/19257
https://cran.r-project.org/web/packages/mlmRev/mlmRev.pdf
豊田本のホームページからプログラムはダウンロードできる（短いので手打ちでもよいが）。下記"SR09.R"は９章で用いるプログラムである。
モジュールlatticeは、プロットに用いる高次グラフ作成ライブラリである。（次のように、それを読み込むと、前提としてlme4とMatrixもよみこまれるようだ。）
生態学分野の人の文献¶
多くの一般化線形モデルの実践的な解説は、生態学者によるものが多いようだ。（フィールドデータの分析に必要なためだろう。）
- 有名なのは、「データ解析のための統計モデリング入門」（久保拓弥、岩波） ... これは通称「緑本」と呼ばれ、人気があるもののようだ。
たくさん出張講義を行っていてそのスライド（PDF)が http://hosho.ees.hokudai.ac.jp/~kubo/ce/InvitedLecture.html に載っている。
2015年の社会心理学会でのセミナーのものがよさそう。（ただ、生態学でのロジスティック関数でフィッティングできる例を用いているのだが、これはちょっと合わない。） - 同じ生態学者によるものだが、北大の隅田さんのものが、シンプルでとっつきやすい。 こちらを勧める http://www.lowtem.hokudai.ac.jp/plantecol/akihiro/sumida-index.html
library("lattice")
source("S09.R") # 豊田本のホームページよりゲットする
mlmRevにあるデータを使った階層線形モデルのテスト¶
豊田本９章にある解析をなぞる。
下では、豊田本のサイトからとったcsvファイルを用いてるが、mlmRevに含まれているExamというデータを使ってみてもよい。
# 使用するライブラリの読み込み
library(mlmRev)
# 表9.1 「9校のデータ」の読み込み
test<-read.csv("dat/nine-school-testdata.csv", header=T)
head(test) # 最初の部分のみ表示 (単に、testとすると全部表示される)
つぎのような２つの列の単回帰を、階層モデルで分析する。
| 列名 | 意味 |
|---|---|
| standLRT | プレテスト |
| normexam | ポストテスト |
ちなみに、次に示すExamデータを使ってみてもよい。（ここでは行わない。）
# mlmRevに含まれている上記と同様のデータでもっと大きなデータ
head(Exam)
# データ全体の単回帰係数
reg.result <- lm(normexam~standLRT,data=test)
reg.result
# 散布図と単回帰直線
coef(reg.result)
plot(test$standLRT,test$normexam)
abline(reg.result)
# 散布図の描画
xyplot(normexam~standLRT|test[,2], type=c("p"),data=test,layout=c(3,3),
xlim=c(-3,3),ylim=c(-3,3),xlab="プリテスト得点",ylab="ポストテスト得点")
# 最小モデル結果
(null<-lmer(normexam~1+(1|school),data=test))
# ランダム切片モデルの結果
(RI<-lmer(normexam ~standLRT+(1|school),data=test))
# ランダム切片・係数モデルの結果（ランダム切片とランダム係数の相関含む）
(RIC<-lmer(normexam ~standLRT+(standLRT+1|school),data=test))
# 学校ごとの切片・傾きと回帰直線（ランダム切片モデル）
coef(RI)
xyplot(fitted(RI)~standLRT,group=school,type=c("l")
,data=test,xlab="プリテスト得点",ylab="ポストテスト得点")
# 学校ごとの切片・傾きと回帰直線（ランダム切片・係数モデル）
coef(RIC) # 切片、傾きの計算結果出力
xyplot(fitted(RIC)~standLRT,group=school,type=c("l"),
data=test,xlab="プリテスト得点",ylab="ポストテスト得点")
lme4のlmerの利用¶
参考： http://www.lowtem.hokudai.ac.jp/plantecol/akihiro/sumida-index.html
上記にならって混合モデル解析lmerを使ってみる。（同じことをやっているのだが。）
モデルは 従属変数 ~ 固定因子 + (ランダム因子) のように書く。
mixedM <- lmer(normexam~standLRT + (standLRT+1|school),data=XYdata)
のように書くと、切片、傾き両方にランダム因子を入れることになる。（が、+1はなくても同じ結果になる。）
# 豊田本「9校のデータ」の読み込み
XYdata<-read.csv("dat/nine-school-testdata.csv", header=T)
library("lme4")
mixedM <- lmer(normexam~standLRT + (standLRT|school),data=XYdata)
summary(mixedM)
RI<-lmer(normexam ~standLRT+(1|school),data=test) # 切片だけをランダム項とした場合
LM<-lm(normexam ~standLRT,data=XYdata) # グループ分けしない単回帰
AIC(LM, RI, mixedM) # AICで妥当性をチェック
BIC(LM, RI, mixedM)
モデルの妥当性のチェックのためには、赤池情報基準(AIC), ベイジアン情報基準(BIC)が用いられる。
それぞれ、小さいほうがモデルとして妥当と言える。（ただし、数値の絶対値に意味はないので、異なるデータ間の比較はしない。）
# データと回帰直線を重ね書きするには、latticeによる描画(xyplotなど)の場合はlatticeExtraに含まれるas.layer関数を利用する
library(latticeExtra)
coef(mixedM)
fit_lines <- xyplot(fitted(mixedM)~standLRT,group=school,type=c("l")
,data=test,xlab="プリテスト得点",ylab="ポストテスト得点")
data_plot <- xyplot(XYdata$normexam~XYdata$standLRT, group=XYdata$school, type=c("p"), add=T)
fit_lines + as.layer(data_plot)
