【統計】一般化線形混合モデル (GLMM)


一般化線形混合モデル (Generalized Linear Mixed Model, GLMM) のメモ。
マルチレベル分析とも言う。

【目次】


・線形回帰モデルでは、誤差が互いに独立に正規分布に従うことを前提としている。
クラスター構造を持つデータでは、データの属する集団内での相関が考えられるため、誤差が互いに独立とは考えにくい。
・GLMMではクラスター毎の相関構造をモデリングできる。
クラスター構造を持つデータは次のようなものが考えられる。A),B)はG行列のランダム効果、C)は反復効果。
 A) 組織をクラスター単位としたデータ(組織毎の性格・性質によるズレ)
 B) 評価者をクラスター単位としたデータ(評価者毎の採点の偏り。甘辛採点)
 C) 個人をクラスター単位とした反復測定データ/パネルデータ


計算式


固定効果モデルと混合効果モデルをそれぞれ記載。
β=固定効果, γ=G行列のランダム効果, G=回帰係数γの分散共分散行列, R=誤差εの分散共分散行列。

固定効果モデルではβとσを推定する。
混合効果モデルではβとσに加え、γ, G, Rを推定する。

【固定効果モデル】
 ■ LM (Linear Model)

 \qquad Y=X' \beta+ \varepsilon   , \quad \varepsilon \sim N(0,\sigma^{2})

 \qquad E(Y) = X' \beta

 \qquad V(Y) = V(\varepsilon)


 ■ GLM (Generalized Linear Model)

 \qquad E(Y) = g^{-1}(X' \beta)

 【統計】一般化線形モデル(GLM) - こちにぃるの日記


【混合効果モデル】
 ■ LMM (Linear Mixed Model)

 \qquad Y_{i} = X_{i}' \beta + Z_{i}' \gamma_{i} + \varepsilon_{i} , \quad \gamma\sim N(0,G),\quad \varepsilon\sim N(0,R_{i})

 \qquad E(Y_{i}) = X_{i}' \beta

 \qquad V(Y_{i}) = V ( Z_{i}' \gamma_{i} ) + V (\varepsilon_{i}) = Z_{i} \ G \ Z_{i}' + R_{i}


  クラスター単位(i番目の水準)の期待値・分散

 \qquad E(Y_{i}|\gamma_{i}) = X_{i}' \beta

 \qquad E(Y_{i}|\gamma_{i}) = V (\varepsilon_{i}) = R_{i}


 ■ GLMM (Generalized Linear Mixed Model)

 \qquad E(Y_{i}|\gamma_{i}) = g^{-1}(X_{i}' \beta + Z_{i}' \gamma_{i})


プログラムコード


*-- ランダム効果を含めたモデル例 ;
proc mixed data=ADS method=reml;
    class Factor Cluster ;
    model Y = Factor / s ddfm=kr ; /* 固定効果 */
    random Intercept Cluster / type=&type. g gcorr ; /* ランダム効果(ランダム切片+ランダム傾き) */
run;
  
*-- 反復効果を含めたモデル例 ;
proc mixed data=ADS method=reml;
    class Group Time SubjectId ;
    model Y = Group Time Group*Time Covariate / s ddfm=kr ; /* 固定効果 */
    repeated Time / subject=SubjectId type=&type. r rcorr ; /* 反復効果 */
    lsmeans Group*Time / pdiff cl ;
run;

 

SASのMIXEDプロシジャで混合効果モデルが作成できる。
・Modelステートメントで「固定効果」を指定する。
・Randomステートメントで「ランダム効果」を指定する。TYPE=でG行列の構造を指定する。オプションにggcorrを入れるとG行列とその相関行列も出力できる。
・Repeatedステートメントで「反復効果」を指定する。TYPE=でR行列の構造を指定する。ARARMAなど前時点の影響度合いを考慮した構造も指定可能。オプションにrrcorrを入れるとG行列とその相関行列も出力できる。
・推定方法のデフォルトは制限付き最尤法 (REstricted Maximum Likelihood, REML) で、最尤法 (Maximum Likelihood, ML) への切り替えもできる。

正規分布以外の分布はGENMODプロシジャやGLIMMIXプロシジャで扱えるようである。


#-- nlme
library(nlme)
fit1 <- nlme::lme(
  data=ads
  ,fixed=Y~Group*Time
  ,random=~1|SubjectId
  ,corr=nlme::corCompSymm(form=~1|SubjectId)
  ,method="REML"
)
summary(fit1)
anova(fit1)  #-- ANOVAテーブル
nlme::VarCorr(fit1)  #-- 共分散パラメータの推定
  
#-- lmer
library(lme4)
library(lmerTest)
fit1 <- lmer(Y~ Factor + (1 | Subject), ADS)) #-- ランダム切片
fit2 <- lmer(Y~ Factor + (1 + Factor | Subject), ADS)) #-- ランダム切片+ランダム傾き
summary(fit1)
anova(fit1, ddf="Kenward-Roger")  #-- ANOVAテーブル
lme4::ranef(fit1)  #-- random effect
lme4::fixef(fit1)  #-- fixed effect

nlmelmerパッケージで混合効果モデルが作成できる。
nlme ではG行列を random、R行列を cor/weights で指定するようだが未勉強。
RPubs - Covariance structures in R


G行列の構造


VCかUNが一般的。SASでのデフォルトはVC(分散成分)。
クラスター3水準(A, B, C)とした場合のG行列:


R行列の構造(反復測定の誤差構造)


誤差構造の一部まとめ。誤差の分散が一定か時点毎に異なるか、誤差の共分散が一定か時点間で異なるかをR行列でデザインする。
R行列はSubject毎(A, B, Cさん3名の例):

 R=\begin{bmatrix} R_{A} & 0 & 0 \\ 0 & R_{B} & 0 \\ 0 & 0 & R_{C} \end{bmatrix}


1被験者5時点とした場合の R_{i}の分散共分散構造。





SAS Help Center
JMP Help
共分散構造


参考


SAS Help Center
http://www.sigmath.es.osaka-u.ac.jp/~kano/research/application/gasshuku02/sas1/MIXED.pdf
https://content.sph.harvard.edu/fitzmaur/ala/lectures.pdf
https://www.mwsug.org/proceedings/2007/stats/MWSUG-2007-S02.pdf
https://agriknowledge.affrc.go.jp/RN/2030812281.pdf
https://www.jstage.jst.go.jp/article/jjb/36/Special_Issue/36_S33/_pdf/-char/ja
https://www.jstage.jst.go.jp/article/jappstat/46/2/46_53/_pdf/-char/ja
一般化線形モデルと一般化線形混合モデル - Qiita
Introduction to SAS® Proc Mixed - ppt download
統計手法アラカルト Mixed Model ~混合モデル~ - ppt download
第4章 MIXED Model 4.1 MIXED Model とは 4.2 反復測定データの分析1 分割法タイプのデータ - ppt download

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。