library(MASS) #一般化線型混合モデルのための関数(glmmPQL)はここ attach(bacteria)# データの詳細は ?bacteria でどうぞ #治療と2週間以上かどうかを考慮した条件付確立分布が二項分布を仮定したモデル fit<-glmmPQL(y ~ trt + I(week > 2), random = ~1 | ID, family = binomial) #random の部分でどの変数にレベル2でのランダムネスを許可するか?を指定 #|以下がレベル2のユニット #I()は条件式で帰ってくるファクターをそのままつかえということ library(boot)#inv.logit というlogitの逆関数がbootライブラリにあります。 p<-inv.logit(predict(fit)) plot(trt,p,m="インフルエンザ菌の存在確率",,xla="治療群",yla="確率")
以下はモデルから予測された確率をプロットしたもの。
それぞれの生徒が卒業した小学校、いま在学している中学校のデータがあるときに、両者がネスト構造になっていない、クロス分類になっている。
XC<-read.table("http://www.okada.jp.org/RWiki/index.php?plugin=attach&openfile=pupcross.dat&refer=%5B%5B%A5%DE%A5%EB%A5%C1%A5%EC%A5%D9%A5%EB%CA%AC%C0%CF%5D%5D",h=T) cons<-rep(1,length(XC)) XC<-cbind(XC,cons) XCdata <-groupedData(ACHIEV~PUPSEX | cons,data=XC) fit<-lme(ACHIEV~PUPSEX,random=pdBlocked(list(pdIdent(~PSCHOOL-1),pdIdent(~SSCHOOL-1))),data=XCdata) summary(fit)
結果
Random effects: Composite Structure: Blocked Block 1: PSCHOOL Formula: ~PSCHOOL - 1 | cons PSCHOOL StdDev: 0.0003159525 Block 2: SSCHOOL Formula: ~SSCHOOL - 1 | cons SSCHOOL Residual StdDev: 0.003773505 0.8593165 Fixed effects: ACHIEV ~ PUPSEX Value Std.Error DF t-value p-value (Intercept) 6.186725 0.05365852 998 115.29810 <.0001 PUPSEX 0.238468 0.05443146 998 4.38107 <.0001 Correlation: (Intr) PUPSEX -0.487
小学校間、中学校間ともに小さい標準偏差しかなく、通常の線型モデルで十分?
Paterson, L. (1991). Socio economic status and educational attainment: a multidimensional and multilevel study. Evaluation and Research in Education 5: 97-121.
の方が意味があるかも?
lme(score~-1+wtn+cwk+I(wtn*gender)+I(cwk*gender),random=~-1+wtn+cwk|school,weights=varIdent(form=~1| wtn),corr=corCompSymm(form=~1 |school/student),data=gcseex)
様々な(マルチレベルモデル用の)データセットもダウンロードできる。