一般化線形混合効果モデル

一般化線形混合効果モデルとは

GLME (一般化線形混合効果) モデルは、正規分布以外の応答変数分布をもつデータの場合に、1 つ以上のグループ化変数によって変化する可能性のある係数を使用し、応答変数と独立変数との関係を記述します。GLME モデルは、グループに収集および要約されたデータに関する GLM (一般化線形モデル) の拡張と考えられます。あるいは、GLME モデルは、応答変数が正常に分布していないデータに関する LME (線形混合効果モデル) の汎化と考えられます。

混合効果モデルは固定効果および変量効果の項を含みます。固定効果の項は通常、モデルの従来の線形回帰部分です。変量効果の項は、母集団から無作為に抽出された個々の実験単位に関連付けられており、応答に影響を与える可能性のあるグループ間の変動を考慮しています。変量効果には事前分布がありますが、固定効果にはありません。

GLME モデルの式

一般化線形混合効果モデルの標準形式は次のとおりです。

$\begin{matrix} y_{i} | b \end{matrix} \sim D i s t r (μ_{i}, \frac{σ^{2}}{w_{i}})$

$g (μ) = X β + Z b + δ,$

ここで

y は n 行 1 列の応答ベクトルであり、y_i はその i 番目の要素です。
b は変量効果ベクトルです。
Distr は、与えられる b に対する y の指定された条件付き分布です。
μ は与えられた b に対する y の条件付き平均であり、μ_i はその i 番目の要素です。
σ² は分散パラメーターです。
w は有効な観測値の重みベクトルであり、w_i は観測値 i の重みです。
- 二項分布では、有効な観測値の重みは、fitglme において 'Weights' によって名前と値のペアの引数が指定された事前の重みに等しく、'BinomialSize' 名前と値のペアの引数で指定された二項分布のサイズで乗算されています。
- 他のすべての分布の場合、有効な観測値の重みは fitglme における 'Weights' の名前と値のペアの引数によって指定された事前の重みと等しくなります。
g(μ) は、平均応答 μ と予測子の線形結合との関係を定義するリンク関数です。
X は n 行 p 列の固定効果計画行列です。
β は p 行 1 列の固定効果ベクトルです。
Z は n 行 q 列の変量効果計画行列です。
b は q 行 1 列の変量効果ベクトルです。
δ はモデルオフセットベクトルです。

平均応答 μ のモデルは以下です。

$μ = g^{- 1} (η),$

ここで、g^-1 はリンク関数 g(μ) の逆関数であり、 ${\hat{η}}_{M E}$ は次の一般化線形混合効果モデルの固定効果および変量効果の線形予測子です。

$η = X β + Z b + δ .$

GLME モデルは β、θ および σ² でパラメーター表現されています。

一般化線形混合効果モデルの仮定は次のとおりです。

変量効果ベクトル b は以下の事前分布になります。
$b | σ^{2}, θ \sim N (0, σ^{2} D (θ)),$
ここで、σ² は分散パラメーターであり、D は、制約のないパラメーターベクトル θ によりパラメーター表現された、対称な半正定値行列です。
観測値 y_i は与えられた b に対して条件付きで独立しています。

モデル近似用データの準備

GLME モデルをデータに当てはめるには、fitglme を使用します。table データ型を使用して入力データをフォーマットします。表の各行は 1 つの観測値、各列は 1 つの予測子変数を表します。table の作成および使用に関する詳細は、table の作成とその table へのデータの代入を参照してください。

入力データは連続したグループ化変数を含むことができます。fitglme は、以下のデータ型を使用する予測子がカテゴリカルであると仮定します。

logical
categorical
文字ベクトルまたは文字配列
string 配列
文字ベクトルの cell 配列

入力データの表が何らかの NaN 値を含む場合、fitglme はデータの行全体を近似から除外します。さらにデータ行を除外する場合、モデルを当てはめる際に fitglme で 'Exclude' の名前と値のペアの引数を使用できます。

モデルの分布タイプを選択

GLME モデルは、応答データが正規分布に従わないときに使用します。したがって、fitglme を使用してモデルを当てはめる場合、'Distribution' 名前と値のペアの引数を使用して応答分布タイプを指定しなければなりません。多くの場合、応答データの型はモデルの適切な分布タイプを示します。

応答データの型	示される応答分布タイプ
任意の実数	`'Normal'`
任意の正の数	`'Gamma'` または `'InverseGaussian'`
任意の非負の整数	`'Poisson'`
ゼロから n までの整数 (n は正の固定値)	`'Binomial'`

モデルのリンク関数の選択

GLME モデルはリンク関数 g を使用し、平均応答と予測子の線形結合との間の関係をマップします。既定では fitglme は、以下の表に示すように、応答データの指定された分布に基づく、事前定義済みの一般に受け入れられたリンク関数を使用します。ただし、fitglme の 'Link' 名前と値のペアの引数を使用して、事前定義済みの関数のリストから異なるリンク関数を指定することもでき、また独自の関数を定義することもできます。

値	説明
`'comploglog'`	`g(mu) = log(-log(1-mu))`
`'identity'`	`g(mu) = mu` 正規分布の正準リンク。
`'log'`	`g(mu) = log(mu)` ポアソン分布の正準リンク。
`'logit'`	`g(mu) = log(mu/(1-mu))` 二項分布の正準リンク。
`'loglog'`	`g(mu) = log(-log(mu))`
`'probit'`	`g(mu) = norminv(mu)`
`'reciprocal'`	`g(mu) = mu.^(-1)`
スカラー値 `P`	`g(mu) = mu.^P`
構造体 `S`	値が関数ハンドルである 4 つのフィールドを含む構造体は、以下を扱います。 `S.Link` — リンク関数 `S.Derivative` — 導関数 `S.SecondDerivative` — 2 次導関数 `S.Inverse` — リンクの逆関数 `'FitMethod'` が `'MPL'` または `'REMPL'` であるか、`S` の指定された分布に対する正準リンクを表す場合、`S.SecondDerivative` の指定を省略できます。

モデルをデータに当てはめる場合、fitglme は既定で正準リンク関数を使用します。

分布	既定のリンク関数
`'Normal'`	`'identity'`
`'Binomial'`	`'logit'`
`'Poisson'`	`'log'`
`'Gamma'`	`-1`
`'InverseGaussian'`	`-2`

リンク関数 'comploglog'、'loglog' および 'probit' は主に二項モデルで便利です。

モデル式の指定

fitglme のモデル仕様では、'y ~ terms' という形式の文字ベクトルまたは string スカラーであるウィルキンソンの表記法を使用します。ここで y は応答変数名であり、terms は以下の表記法で記述されます。

ウィルキンソンの表記法	標準表記の因子
`1`	定数 (切片) 項
`X^k`、`k` は正の整数	`X`, `X²`, ..., `X^k`
`X1 + X2`	`X1`, `X2`
`X1*X2`	`X1`, `X2`, `X1.*X2 (element-wise multiplication of X1 and X2)`
`X1:X2`	`X1.*X2` のみ
`- X2`	`X2` は含めない
`X1*X2 + X3`	`X1`, `X2`, `X3`, `X1*X2`
`X1 + X2 + X3 + X1:X2`	`X1`, `X2`, `X3`, `X1*X2`
`X1X2X3 - X1:X2:X3`	`X1`, `X2`, `X3`, `X1X2`, `X1X3`, `X2*X3`
`X1*(X2 + X3)`	`X1`, `X2`, `X3`, `X1X2`, `X1X3`

式には既定で定数 (切片) 項が含まれます。モデルから定数項を除外するには、式に –1 を含めます。

一般化線形混合効果モデルについては、式の仕様は 'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)' 形式であり、ここで fixed および random は固定効果項と変量効果項をそれぞれ含みます。

入力データの表は以下を含むと想定します。

応答変数 y
予測子変数 X1, X2, ..., XJ。ここで J は予測子変数の (連続したグループ化変数を含む) 合計数。
グループ化変数 g1, g2, ..., gR。ここで R はグループ化変数の数。

XJ および gR のグループ化変数は、categorical 配列、logical 配列、文字配列、string 配列、または文字ベクトルの cell 配列が可能です。

この場合、'y ~ fixed + (random1|g1) + ... + (randomR|gR)' の形式の式において、項 fixed は固定効果の計画行列 X の仕様に対応し、random1 はグループ化変数 g1 に対応する変量効果の計画行列 Z1 の仕様であり、同様に randomR はグループ化変数 gR に対応する変量効果の計画行列 ZR の仕様です。fixed 項および random 項は、以下のようにウィルキンソンの表記法で表現できます。

式	説明
`'y ~ X1 + X2'`	切片の固定効果: `X1` および `X2`。この式は `'y ~ 1 + X1 + X2'` と等価です。
`'y ~ -1 + X1 + X2'`	切片のない `X1` と `X2` の固定効果。`-1` を含めることによって暗黙的な切片の項は抑制されます。
`'y ~ 1 + (1 \| g1)'`	グループ化変数 `g1` の水準ごとの切片の固定効果と切片の変量効果の和。
`'y ~ X1 + (1 \| g1)'`	固定傾きのランダム切片モデル。
`'y ~ X1 + (X1 \| g1)'`	相関があり得るランダムな切片と傾き。この式は `'y ~ 1 + X1 + (1 + X1\|g1)'` と等価です。
`'y ~ X1 + (1 \| g1) + (-1 + X1 \| g1)'`	切片および傾きに関する独立した変量効果の項。
`'y ~ 1 + (1 \| g1) + (1 \| g2) + (1 \| g1:g2)'`	`g1` と `g2` に対する独立したメイン効果のあるランダムな切片モデル + 独立した交互作用効果。

たとえば、標本データ mfr は世界で 50 の工場を操業する製造企業からのシミュレートされたデータを含みます。各工場は、完成品の生産のためにバッチプロセスを実行します。同社は各バッチの欠陥数を減少させるために新たな製造プロセスを開発しました。新しいプロセスの効果をテストするため、同社は実験に参加させる 20 工場を無作為に選びました。10 工場では新プロセスを実施しますが、残りの 10 工場では旧プロセスの実行を続けます。各 20 工場で、同社は 5 つのバッチ (合計 100 バッチ) を実行し、バッチごとに処理時間 (time_dev)、温度 (temp_dev)、欠陥数 (defects) および原材料の供給業者を示すカテゴリカル変数 (supplier) のデータを記録しました。

新しいプロセス (予測子変数 newprocess で表現されます) が欠陥数を有意に減少させたかどうかを判断するために、newprocess、time_dev、temp_dev および supplier を固定効果予測子として使用し、GLME モデルを当てはめます。factory 別にグループ化された変量効果の切片を含めます。これは工場特有の変動に起因する品質に差がある可能性を考慮するためです。応答変数 defects はポアソン分布になります。

欠陥数はポアソン分布を使用してモデル化できます

$d e f e c t s_{i j} ~ P o i s s o n (μ_{i j})$

これは一般化線形混合効果モデルに対応します

$\begin{array}{l} \log (μ_{i j}) = β_{0} + β_{1} n e w p r o c e s s_{i j} + β_{2} t i m e_d e v_{i j} \\ + β_{3} t e m p_d e v_{i j} + β_{4} s u p p l i e r_C_{i j} + β_{5} s u p p l i e r_B_{i j} + b_{i}, \end{array}$

ここで

defects_ij は、バッチ j (ここで j = 1, 2, ..., 5) 処理中に工場 i (ここで i = 1, 2, ..., 20) で実行されたバッチにおいて観察された欠陥数です。
μ_ij はバッチ j 処理中の工場 i に対応する欠陥の平均数です。
supplier_C_ij および supplier_B_ij はそれぞれ、会社 C または B がバッチ j 処理中に工場 i で実行されたバッチ用の加工化学薬品を供給したかどうかを意味するダミー変数です。
b_i ~ N(0,σ_b²) は、工場特有の品質変動に相当する、各工場 i の変量効果の切片です。

ウィルキンソンの表記法を使用し、このモデルを以下のように指定します。

'defects ~ 1 + newprocess + time_dev + temp_dev + supplier + (1|factory)'

fitglme を使用してモデルをデータに当てはめる際、応答変数のポアソン分布を考慮するには、'Distribution' の名前と値のペアの引数を 'Poisson' として指定します。既定では、fitglme はポアソン分布の応答変数に対数リンク関数を使用します。

モデルの表示

近似関数 fitglme の出力は一般化線形混合効果モデルの情報を提供します。

製造実験データ mfr を使用し、newprocess、time_dev、temp_dev および supplier を固定効果予測子とするモデルを当てはめます。応答分布にポアソンを、リンク関数に対数を、近似メソッドにラプラスを指定します。

load mfr

glme = fitglme(mfr,...
		'defects ~ 1 + newprocess + time_dev + temp_dev + supplier + (1|factory)',...
		'Distribution','Poisson','Link','log','FitMethod','Laplace',...
		'DummyVarCoding','effects')

glme = 


Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations             100
    Fixed effects coefficients           6
    Random effects coefficients         20
    Covariance parameters                1
    Distribution                    Poisson
    Link                            Log   
    FitMethod                       Laplace

Formula:
    defects ~ 1 + newprocess + time_dev + temp_dev + supplier + (1 | factory)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    416.35    434.58    -201.17          402.35  

Fixed effects coefficients (95% CIs):
    Name                 Estimate     SE          tStat       DF    pValue    
    '(Intercept)'           1.4689     0.15988      9.1875    94    9.8194e-15
    'newprocess'          -0.36766     0.17755     -2.0708    94      0.041122
    'time_dev'           -0.094521     0.82849    -0.11409    94       0.90941
    'temp_dev'            -0.28317      0.9617    -0.29444    94       0.76907
    'supplier_C'         -0.071868    0.078024     -0.9211    94       0.35936
    'supplier_B'          0.071072     0.07739     0.91836    94       0.36078


    Lower        Upper    
       1.1515       1.7864
     -0.72019    -0.015134
      -1.7395       1.5505
      -2.1926       1.6263
     -0.22679     0.083051
    -0.082588      0.22473

Random effects covariance parameters:
Group: factory (20 Levels)
    Name1                Name2                Type         Estimate
    '(Intercept)'        '(Intercept)'        'std'        0.31381 

Group: Error
    Name                      Estimate
    'sqrt(Dispersion)'        1

Model information 表は標本データの観測値の合計 (100)、固定効果および変量効果係数の数 (それぞれ 6 および 20)、共分散パラメーターの数 (1) を表示しています。また、応答変数は Poisson 分布であり、リンク関数は Log であり、近似メソッドが Laplace であることもわかります。

Formula はウィルキンソンの表記法によるモデル仕様を示します。

Model fit statistics 表はモデルの適合度の評価に使用された統計を表します。これには赤池情報量基準 (AIC)、ベイズ情報量基準 (BIC) 値、対数尤度 (LogLikelihood) および逸脱度 (Deviance) の値が含まれます。

Fixed effects coefficients 表は、fitglme が 95% の信頼区間を返したことを示します。これには固定効果予測子ごとに 1 行が含まれ、各列にはその予測子に対応する統計が含まれます。列 1 (Name) には各固定効果係数の名前が含まれ、列 2 (Estimate) にはその推定値が含まれ、列 3 (SE) には係数の標準誤差が含まれます。列 4 (tStat) には係数が 0 に等しいという仮説検定のための t 統計量が含まれます。列 5 (DF) と列 6 (pValue) には t 統計量に対応する自由度と p 値がそれぞれ含まれます。最後の 2 列 (Lower および Upper) には、各固定効果係数の 95% 信頼区間の下限と上限がそれぞれ表示されます。

Random effects covariance parameters は各グループ化変数 (ここでは factory のみ) の表を表示します。これには水準の総数 (20)、共分散パラメーターのタイプおよび推定値が含まれます。ここでの std は、工場の予測子に関連付けられている変量効果の標準偏差が fitglme から返されることを示します。この推定値は 0.31381 です。また、誤差パラメーターのタイプ (ここでは分散パラメーターの平方根) およびその推定値 1 を含む表も表示します。

fitglme により生成される標準表示は変量効果パラメーターの信頼区間を指定しません。covarianceParameters を使用して、これらの値を計算し表示します。

モデルの操作

fitglme を用いて GLME モデルを作成した後で、モデルを操作するための追加の関数を使用できます。

係数および信頼区間の検査および検定

固定効果および変量効果係数の推定値、共分散パラメーター、計画行列および関連する統計を抽出するには、次の方法を用います。

fixedEffects は、推定固定効果係数および関連する統計を近似されたモデルから抽出します。関連する統計には、標準誤差、t 統計量、自由度、各パラメーターが 0 に等しいかどうかを仮説検定するための p 値、および信頼区間が含まれます。
randomEffects は推定変量効果係数および関連する統計を近似された GLME モデルから抽出します。関連する統計には、それぞれの変量効果の推定 EBP (経験的ベイズ予測子)、共分散パラメーターおよび応答が与えられた CMSEP (予測の条件付き平均二乗誤差) の平方根、t 統計量、自由度、および各パラメーターが 0 に等しいかどうかを仮説検定するための p 値、ならびに信頼区間が含まれます。
covarianceParameters は推定共分散パラメーターおよび関連する統計を近似された GLME モデルから抽出します。関連する統計には、共分散パラメーターの推定値および信頼区間が含まれます。
designMatrix は固定効果および変量効果の計画行列またはその指定されたサブセットを、近似された GLME モデルから抽出します。

固定効果および変量効果係数の有意性に関して、カスタマイズされた仮説検定を実行し、カスタムの信頼区間を計算するには、以下を行います。

anova は固定効果項に対する周辺 F 検定 (仮説検定) を実施し、固定効果項を表すすべての係数が 0 に等しいかどうかを判定します。anova を使用すると、カテゴリカル予測子の係数の有意性を結合して検定できます。
coefCI は近似 GLME モデルから固定効果および変量効果パラメーターの信頼区間を計算します。既定では、fitglme は信頼区間の 95% を計算します。coefCI を使用し、別の信頼水準での境界を計算します。
coefTest は近似された一般化線形混合効果モデルの固定効果または変量効果ベクトルへのカスタムな仮説検定を実施します。たとえば、対比行列を指定できます。

新たな応答値を生成しモデルを再度当てはめる

近似された GLME モデルに基づき、近似された応答、予測された応答および無作為の応答を含む新たな応答値を生成するには、以下を行います。

fitted は元の予測子の値を使用して近似された応答値を計算し、近似されたモデルから推定係数およびパラメーター値を計算します。
predict は元の予測子の値、または新しい予測子の値のいずれかを使用し、応答の予測された条件付き平均または周辺平均を計算し、近似されたモデルから推定係数およびパラメーター値を計算します。
random は、近似されたモデルから無作為の応答を生成します。
refit は、元のモデルおよび新たな応答ベクトルに基づき、新たな近似された GLME モデルを作成します。

残差の検査および可視化

近似された GLME モデルから残差を抽出および可視化するには、以下を行います。

residuals は、近似されたモデルから生の残差もしくはピアソン残差を抽出します。さらに、条件付き残差または限界残差を計算するかどうかも指定できます。
plotResiduals は、近似されたモデルから生の残差またはピアソン残差により、次のプロットを作成します。
- 残差のヒストグラム
- 残差と近似値の散布図
- 残差とラグ付き残差の散布図

参考

fitglme | GeneralizedLinearMixedModel

トピック

一般化線形混合効果モデルの当てはめ