混合効果モデルと固定効果モデルの比較

前回の続きです。

変量効果モデルで推定した場合と、固定効果モデルで交互作用として推定した場合の違いを絵でみるとこのようになります。

成長曲線のデータで、8～14才を追跡し身長が伸びた分を記録しているそうです。

以下のグラフは9人分の図。

• 変量効果

• 固定効果

まぁ、あんまり変わらないね。。。

変量効果の方が欠測値を扱いやすいとかっていう利点があった（固定効果でも推定はできるけど）。

あまりに欠測が多かったりカテゴリが多いとMCMCで推定する必要が出てきます。

結果はあまり変わらないけど、変量効果モデルを使うと「お、統計を知ってるな」って思われるかも。

プログラムはこちら↓

#------変量効果と固定効果の比較
#---変量効果
library(nlme)
formula(Orthodont)

#交互作用モデルではFactorの値が使われるので文字にしておく
Orthodont$Subject <- as.character(Orthodont$Subject)

fm1      <- lme(distance ~ age, data = Orthodont)	#subject毎にageと切片の推定値を算出
IntMix   <- fm1$coefficients$fixed[1] + fm1$coefficients$random$Subject[, 1]		#混合効果モデルでの切片
SlopeMix <- fm1$coefficients$fixed[2] + fm1$coefficients$random$Subject[, 2]		#混合効果モデルでの傾き

MixedPlot <- function(Group){
	plot(Orthodont[Orthodont[ ,3] == Group, c(2, 1)], main=Group)
	abline(a=IntMix[names(IntMix) == Group], b=SlopeMix[names(SlopeMix) == Group])
}

par(mfrow=c(3, 3), pty="m")
MixedPlot("M01")
MixedPlot("M02")
MixedPlot("M03")
MixedPlot("M04")
MixedPlot("M05")
MixedPlot("M06")
MixedPlot("M07")
MixedPlot("M08")
MixedPlot("M09")


#---固定効果
fm2      <- lm(distance ~ age + age * Subject, data = Orthodont)	#固定効果で交互作用を指定したモデル
IntFix   <- fm2$coefficients[1] + fm2$coefficients[3:28]
SlopeFix <- fm2$coefficients[2] + fm2$coefficients[29:54]

FixedPlot <- function(Group){
	plot(Orthodont[Orthodont[ ,3] == Group, c(2, 1)], main=Group)
	#---交互作用モデルでは推定パラメータ名に変数名Subjectが付いているので名前をあわせる
	SGroup <- paste("Subject", Group, sep="")
	SAGroup <- paste("age:Subject", Group, sep="")
	abline(a=IntFix[names(IntFix) == SGroup], b=SlopeFix[names(SlopeFix) == SAGroup])
}

par(mfrow=c(3, 3), pty="m")

#M01を基底にして交互作用を推定しているので、M01だけはそのまま描く
plot(Orthodont[Orthodont[ ,3] == "M01", c(2, 1)], main="M01")
SGroup <- paste("Subject", "M01", sep="")
SAGroup <- paste("age:Subject", "M01", sep="")
abline(a=IntFix[names(IntFix) == SGroup], b=SlopeFix[names(SlopeFix) == SAGroup])

FixedPlot("M02")
FixedPlot("M03")
FixedPlot("M04")
FixedPlot("M05")
FixedPlot("M06")
FixedPlot("M07")
FixedPlot("M08")
FixedPlot("M09")