混合効果モデル（変量効果モデル、mixed effect model）について

混合効果モデルについてです。

論文1：Random-Effects Models for Longitudinal Data. Laird NM and Ware JH. 1982. Biometrics. 38. 963–74.

論文2：Approximate Inference in Generalized Linear Mixed Models. Breslow NE and Clayton DG. JASA, 1993. 88;421. 9-25.

http://www.stat.ubc.ca/~ruben/papers/BreslowClaytonPQLpaper.pdf

論文2は論文1の混合効果モデルを一般化した論文で、このBreslow先生とClayton先生はかなり著名な先生です。混合効果モデルは機械学習の教科書には載ってなく、統計の基礎的な教科書にもあまり載ってません。Armitageの教科書（以下）にかろうじて載っていましたが、生物統計以外の分野にはあまり普及してないのかもしれません（追記参照）。

この本は日本語の訳本もあり、医療統計学の最も有名な本の1つです。混合効果モデルの教科書は次のものがあります。

教科書：Linear Mixed Models for Longitudinal Data. Verbeke G and Molenberghs G. Springer. 2009.

東京大学の松山裕先生がこの本を訳されていてとても分かりやすい良本でしたが、残念ながら絶版です（本の名前は「医学統計のための線型混合モデル」）。

「変量効果の推定とBLUP法」という日本語の教科書もありますが、難解です。

以下、説明です。

【混合効果モデル】

• 固定効果モデル
  • 通常の線形モデル
  • y=βx+ε
  • パラメータβは固定値（分布しない）

• 変量効果モデル
  • 混合効果モデル
  • y=βx+uz+ε
  • パラメータuは変量効果（ランダムに分布する）
    • zはクラスターを示す
    • クラスター毎にパラメータuが推定される
    • uが全体として正規分布に従う
  • u～正規分布、y～正規分布という2段階のモデルになっている
    • 階層モデリング（hierarchical modeling、マルチレベルモデル）とも呼ばれる
    • ベイズ的階層モデルとも関連が深い
  • パラメータが分布するという概念は頻度論者としてはとても気持ちが悪い

• クラスターが存在するようなデータで使われる
  • 個人経時データ（個人を変量効果に入れる）
  • 多施設データ（施設を変量効果に入れる）
  • メタアナリシス

【推定方法】

• 固定効果、分散成分
  • 最尤推定法（maximum likelihood、ML）と制限付き最尤推定法（restricted maximum likelihood、REML）を反復して計算
  • 分散成分の初期値を決定→固定効果をMLで計算→分散成分をREMLで再計算→…収束するまで反復

• 変量効果
  • ベイズの定理を用いて事後分布の事後平均（期待値）を計算する
  • 最良線形不偏予測（best linear unbiased prediction、BLUP）になっている
  • 経験ベイズ推定量と解釈できる

【一般化線形混合モデル】

• generalized linear mixed model、GLMM
• 線形混合モデルの結果変数をリンク関数で変換したもの
• 式が複雑になり積分計算できないので次のような方法で推定する
  • EMアルゴリズム
  • ギブスサンプリング
  • PQL(Penalized Quasi-Likelihood) method
    • 階層モデルで変量効果の推定に興味
  • MQL(Marginal Quasi-Likelihood) method
    • 共変量と結果変数間の関係に興味

• Subject Specific（SS）モデル
  • 階層モデル
  • 変量効果の推定に重点を置く

• Population Average（PA）モデル
  • 変量効果に重点を置かず固定効果に大きな興味がある

【Rでの例】

• パッケージ：nlme
  • 著者にR Core Teamが入っているのでこのパッケージで問題ないであろう
• サンプルデータ：Orthodont
  • 27人の子供のOrthodont length（歯の矯正？）を測定
  • 8、10、12、14歳の4回
  • lengthを予測するモデルをいくつか試してみて、固定効果モデルと混合効果モデルの違いを確認する

• groupedData関数で変量効果を指定した後でlme関数で推定する

• lme4
  • 計算速度が速くなった？

次のムービーは以下のモデル1～6を当てはめたときの推定値の変化を表しています。これは分散分析（analysis of variance、ANOVA）の平方和の分解にも相当します。

1. y = mu
   • 平均値のみのモデル
2. y = mu + age
   • モデル1+年齢
3. y = mu + subject
   • モデル1+対象者
4. y = mu + age + subject
   • モデル2+対象者
5. y = mu + age + subject + age*subject
   • モデル4+年齢・対象者交互作用
6. y = mu + age + subject(mixed)
   • モデル2+対象者変量効果

• モデル6が混合効果モデルであり、固定効果モデル5と同等である
• 混合モデルと固定モデルの違いは、推定値が縮小（shrink）していること
  • 混合効果のパラメータは経験ベイズ推定値であるため、このような性質になっている

固定効果の精度が高くなる（有意になりやすくなる）傾向があること。ちなみに、クラスターの数が多いと推定に結構な時間がかかります。

追記：@phosphor_mさんとの議論

• 社会科学の分野では良く使われる（マルチレベルモデル）
• 国際比較をする際に調整目的で変量効果を入れる
• 交互作用モデルと混合効果モデルは目的に応じて使い分けるべき？
  • 交互作用モデル：交互作用の値自体に興味がある
  • 混合効果モデル：確認したいのは主効果のみで、クラスターの効果は調整したい
• 結果が大きく変わることは無いが、ニュアンスが違う

library(nlme)
str(Orthodont)

Orthodont <- groupedData(distance ~ age | Subject, data = Orthodont)


lm1 <- lm(distance  ~ age, data = Orthodont)
lm2 <- lm(distance  ~ age + Subject, data = Orthodont)
lm3 <- lm(distance  ~ age + Subject + age*Subject, data = Orthodont)
lm4 <- lm(distance  ~ Subject, data = Orthodont)
fm1 <- lme(distance ~ age, data = Orthodont)


summary(lm1)
summary(lm2)
summary(lm3)
summary(fm1)


plot(Orthodont$distance, rep(mean(Orthodont$distance), nrow(Orthodont)), 
	lwd = 2, main = "y = mu", xlab = "observed length", ylab = "predicted length", 
	xlim = c(15, 35), ylim = c(15, 35), pch = 19)

#---lm1
plot(Orthodont$distance, predict(lm1, Orthodont), 
	lwd = 2, main = "y = mu + age", xlab = "observed length", ylab = "predicted length", 
	xlim = c(15, 35), ylim = c(15, 35), col = Orthodont$age, pch = 19)
legend("topleft", names(table(Orthodont$age)), col = names(table(Orthodont$age)), pch = 19)


#---lm4
plot(Orthodont$distance, predict(lm4, Orthodont), 
	lwd = 2, main = "y = mu + subject", xlab = "observed length", ylab = "predicted length", 
	xlim = c(15, 35), ylim = c(15, 35), col = Orthodont$age, pch = as.numeric(Orthodont$Subject))
legend("topleft", names(table(Orthodont$age)), col = names(table(Orthodont$age)), pch = 19)
legend("bottomright", names(table(Orthodont$Subject)), pch = 1:length(names(table(Orthodont$Subject))))


#---lm2
plot(Orthodont$distance, predict(lm2, Orthodont), 
	lwd = 2, main = "y = mu + age + subject", xlab = "observed length", ylab = "predicted length", 
	xlim = c(15, 35), ylim = c(15, 35), col = Orthodont$age, pch = as.numeric(Orthodont$Subject))
legend("topleft", names(table(Orthodont$age)), col = names(table(Orthodont$age)), pch = 19)
legend("bottomright", names(table(Orthodont$Subject)), pch = 1:length(names(table(Orthodont$Subject))))


#---lm3
plot(Orthodont$distance, predict(lm3, Orthodont), 
	lwd = 2, main = "y = mu + age + subject + age*subject", xlab = "observed length", ylab = "predicted length", 
	xlim = c(15, 35), ylim = c(15, 35), col = Orthodont$age, pch = as.numeric(Orthodont$Subject))
legend("topleft", names(table(Orthodont$age)), col = names(table(Orthodont$age)), pch = 19)
legend("bottomright", names(table(Orthodont$Subject)), pch = 1:length(names(table(Orthodont$Subject))))


#---fm1
plot(Orthodont$distance, predict(fm1, Orthodont), 
	lwd = 2, main = "y = mu + age + subject(mixed)", xlab = "observed length", ylab = "predicted length", 
	xlim = c(15, 35), ylim = c(15, 35), col = Orthodont$age, pch = as.numeric(Orthodont$Subject))
legend("topleft", names(table(Orthodont$age)), col = names(table(Orthodont$age)), pch = 19)
legend("bottomright", names(table(Orthodont$Subject)), pch = 1:length(names(table(Orthodont$Subject))))