Main Content

線形混合効果モデルでのパラメーター推定

線形混合効果モデルは次のような形式になります。

y=Xβfixed+Zbrandom+εerror,

ここで

  • y は n 行 1 列の応答ベクトル、n は観測数です。

  • X は n 行 p 列の固定効果計画行列です。

  • β は p 行 1 列の固定効果ベクトルです。

  • Z は n 行 q 列の変量効果計画行列です。

  • b は q 行 1 列の変量効果ベクトルです。

  • ε は n 行 1 列の観測誤差ベクトルです。

変量効果ベクトル b と誤差ベクトル ε は、次のように独立した事前分布になっていると仮定しています。

b~N(0,σ2D(θ)),ε~N(0,σI2),

ここで、D は分散成分ベクトル θ によってパラメーター表現されている対称な半正定値行列、I は n 行 n 列の単位行列、σ2 は誤差の分散です。

このモデルで推定するパラメーターは、固定効果係数 β および分散成分の θ と σ2 です。線形混合効果モデルにおけるパラメーター推定で最も一般的に使用される 2 つの方法は、最尤法と制限付き最尤法です。

最尤法 (ML)

最尤推定法には、回帰係数と分散成分の両方が含まれており、これらは尤度関数の固定効果と変量効果の項です。

前に定義した線形混合効果モデルでは、与えられた β、b、θ および σ2 に対して、応答変数 y の条件応答は次のようになります。

y|b,β,θ,σ2~N(Xβ+Zb,σ2In).

与えられた β、θ および σ2 に対して y の尤度は次のようになります。

P(y|β,θ,σ2)=P(y|b,β,θ,σ2)P(b|θ,σ2)db,

ここで

P(b|θ,σ2)=1(2πσ2)q21|D(θ)|12exp{12σ2bTD1b}andP(y|b,β,θ,σ2)=1(2πσ2)n2exp{12σ2(yXβZb)T(yXβZb)}.

Λ(θ) を D(θ) の下三角コレスキー因子、Δ(θ) を Λ(θ) の逆数と仮定します。このとき次のようになります。

D(θ)1=Δ(θ)TΔ(θ).

次のように定義します。

r2(β,b,θ)=bTΔ(θ)TΔ(θ)b+(yXβZb)T(yXβZb),

そして、b* が次の式を満たす b の値であると仮定します。

r2(β,b,θ)b|b*=0

(特定の β と θ に対して)。この場合、尤度関数は次のようになります。

P(y|β,θ,σ2)=(2πσ2)n2|D(θ)|12exp{12σ2r2(β,b*(β),θ)}1|ΔTΔ+ZTZ|12.

P(y|β,θ,σ2) は特定の θ に対して、β と σ2 によって最初に最大化されます。これにより最適解 β^(θ) および σ^2(θ) が θ の関数として得られます。これらの解を尤度関数に代入すると、P(y|β^(θ),θ,σ^2(θ)) になります。この式はプロファイル尤度と呼ばれ、β と σ2 がプロファイルされています。P(y|β^(θ),θ,σ^2(θ)) は θ の関数であり、このアルゴリズムでは次にこの式を θ に関して最適化します。θ の最適推定値が見つかると、β と σ2 の推定値は β^(θ) および σ^2(θ) によって与えられます。

ML 法では、分散成分の推定で β を固定の未知の量として扱いますが、固定効果の推定によって失われる自由度は考慮されません。これは ML の推定が小さい分散によってバイアスされる原因になります。ML が REML よりも優れている点の 1 つは、2 つのモデルを固定効果および変量効果の項について比較できることです。一方、REML を使用してパラメーターを推定する場合、同じ固定効果計画をもつ変量効果の項において 2 つのモデルが入れ子になっている場合のみ、モデルを比較できます。

制限付き最尤法 (REML)

制限付き最尤推定法には分散成分のみが含まれます。分散成分とは線形混合効果モデル内の変量効果の項をパラメーター表現するパラメーターです。β は 2 番目のステップで推定されます。β が一様な変則事前分布であると仮定して尤度 P(y|β,θ,σ2) を β について積分したものが制限付き尤度 P(y|θ,σ2) になります。つまり次のようになります。

P(y|θ,σ2)=P(y|β,θ,σ2)P(β)dβ=P(y|β,θ,σ2)dβ.

このアルゴリズムでは、はじめに σ^R2 をプロファイルしてから残りの目的関数を θ に関して最大化することにより θ^R を求めます。その後、σ2 に関して制限付き尤度を最大化して σ^R2 を求めます。そして、次の事後分布に関する期待値を求めることにより β を推定します。

P(β|y,θ^R,σ^R2).

REML は固定効果の推定によって失われる自由度を考慮するため、変量効果の分散を推定するバイアスがやや小さくなります。θ および s2 の推定値は β の値に対して不変であり、ML 推定に比べるとデータの外れ値にあまり影響されません。ただし、REML を使用してパラメーターを推定する場合、同一の固定効果計画行列をもつ変量効果の項において 2 つのモデルが入れ子になっている場合のみ、モデルを比較できます。

参照

[1] Pinherio, J. C., and D. M. Bates. Mixed-Effects Models in S and S-PLUS. Statistics and Computing Series, Springer, 2004.

[2] Hariharan, S. and J. H. Rogers. “Estimation Procedures for Hierarchical Linear Models.” Multilevel Modeling of Educational Data (A. A. Connell and D. B. McCoach, eds.). Charlotte, NC: Information Age Publishing, Inc., 2008.

[3] Raudenbush, S. W. and A. S. Bryk. Hierarchical Linear Models: Applications and Data Analysis Methods, 2nd ed. Thousand Oaks, CA: Sage Publications, 2002.

[4] Hox, J. Multilevel Analysis, Techniques and Applications. Lawrence Erlbaum Associates, Inc, 2002.

[5] Snidjers, T. and R. Bosker. Multilevel Analysis. Thousand Oaks, CA: Sage Publications, 1999.

[6] McCulloch, C.E., R. S. Shayle, and J. M. Neuhaus. Generalized, Linear, and Mixed Models. Wiley, 2008.

参考

| |

関連するトピック