designMatrix

クラス: GeneralizedLinearMixedModel

固定効果と変量効果の計画行列

構文

D = designMatrix(glme)

D = designMatrix(glme,'Fixed')

D = designMatrix(glme,'Random')

Dsub = designMatrix(glme,'Random',gnumbers)

[Dsub,gnames] = designMatrix(glme,'Random',gnumbers)

説明

D = designMatrix(glme) または D = designMatrix(glme,'Fixed') は、一般化線形混合効果モデル glme の固定効果計画行列を返します。

D = designMatrix(glme,'Random') は、一般化線形混合効果モデル glme の変量効果計画行列を返します。

Dsub = designMatrix(glme,'Random',gnumbers) は、gnumbers で示されているグループ化変数に対応する一般化線形混合効果モデル glme の変量効果計画行列のサブセットを返します。

[Dsub,gnames] = designMatrix(glme,'Random',gnumbers) も gnumbers に対応するグループ化変数名を返します。

入力引数

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

GeneralizedLinearMixedModel オブジェクト

一般化線形混合効果モデル。GeneralizedLinearMixedModel オブジェクトとして指定します。このオブジェクトのプロパティとメソッドについては、GeneralizedLinearMixedModel を参照してください。

gnumbers — グループ化変数の数

整数値の配列

グループ化変数の数。範囲 [1,R] の要素を含む整数値の配列として指定します。ここで、R は、一般化線形混合効果モデル glme のグループ化変数を含む cell 配列の長さです。

たとえば、グループ化変数 g₁、g₃ および g_r は [1,3,r] として指定できます。

データ型: single | double

出力引数

D — 計画行列

行列

一般化線形混合効果モデル glme の計画行列。次のいずれかとして返されます。

• 固定効果の計画行列 — glme の固定効果計画行列で構成される n 行 p 列の行列。ここで、n は観測値の数、p は固定効果の項の数です。D 内の固定効果項の次数は、GeneralizedLinearMixedModel オブジェクト glme の CoefficientNames プロパティにおける項の次数と一致します。

• 変量効果の計画行列 ― glme の変量効果計画行列で構成される n 行 k 列の行列。ここで、k は length(B) に等しく、B は一般化線形混合効果モデル glme の変量効果の係数ベクトルです。変量効果計画行列はスパース行列として返されます。詳細は、スパース行列を参照してください。

  glme には R 個のグループ化変数 g₁、g₂、...、g_R (水準はそれぞれ m₁、m₂、...、m_R) が含まれ、q₁、q₂、...、q_R はそれぞれ g₁、g₂、...、g_R と関連付けられる変量効果ベクトルの長さである場合、B は長さが q₁*m₁ + q₂*m₂ + ... + q_R*m_R の列ベクトルです。

  B は、各グループ化変数の各水準に対応する変量効果ベクトルの経験的ベイズ予測子を、次のように連結することで作成されます。[g1level1; g1level2; ...; g1levelm1; g2level1; g2level2; ...; g2levelm2; ...; gRlevel1; gRlevel2; ...; gRlevelmR]' として連結することによって作成されます。

データ型: single | double

Dsub — 変量効果の計画行列の部分行列

行列

gnumbers で指定されるグループ化変数に対応する変量効果計画行列の部分行列。n 行 k 列の行列として返されます。ここで、k は列ベクトル Bsub の長さです。

Bsub に含まれる連結された変量効果ベクトルの経験的ベイズ予測子は、gnumbers によって指定されたグループ化変数の各水準に対応します。

たとえば、gnumbers が [1,3,r] である場合、これはグループ化変数 g₁、g₃、g_r に対応します。このとき、Bsub には、グループ化変数 g₁、g₃、g_r の各水準に対応する変量効果ベクトルの経験的ベイズ予測子が含まれ、次のようになります。

[g1level1; g1level2; ...; g1levelm1; g3level1; g3level2; ...; g3levelm3; grlevel1; grlevel2; ...; grlevelmr]'.

したがって、Dsub*Bsub は、グループ化変数 g₁、g₃、g_r に対応するすべての変量効果の、glme の応答に対する寄与を表します。

gnumbers が空の場合は、Dsub は変量効果計画行列の全体になります。

データ型: single | double

gnames — グループ化変数の名前

k 行 1 列の cell 配列

設計タイプが 'Random' の場合は、gnumbers の整数に対応するグループ化変数の名前。k 行 1 列の cell 配列として返されます。設計タイプが 'Fixed' の場合は、gnames は空の行列 [] です。

データ型: cell

例

固定効果と変量効果の計画行列の取得

標本データを読み込みます。

load mfr

このシミュレーションされたデータは、世界中で 50 の工場を操業している製造企業から取得しており、各工場が完成品の生産のためにバッチ処理を実行しています。同社は各バッチの欠陥数を減少させるために新たな製造プロセスを開発しました。新しいプロセスの効果をテストするため、同社は実験に参加させる 20 工場を無作為に選びました。10 工場では新プロセスを実施しますが、残りの 10 工場では旧プロセスの実行を続けます。各 20 工場で、同社は 5 つのバッチ (合計 100 バッチ) を実行し以下のデータを記録しました。

• 新しいプロセスがバッチに使用されたかどうかを示すフラグ (newprocess)

• 各バッチの処理時間。時間単位 (time)

• バッチの温度。摂氏 (temp)

• バッチで使用する化学薬品の供給業者 (A、B または C) を示すカテゴリカル変数 (supplier)

• バッチ内の欠陥数 (defects)

またデータに含まれる time_dev と temp_dev は、摂氏 20 度で 3 時間の標準プロセスから得られる時間と温度の絶対偏差をそれぞれ表します。

固定効果予測子として newprocess、time_dev、temp_dev および supplier を使用して一般化線形混合効果モデルを当てはめます。工場特有の変動に起因して品質に差がある可能性を考慮するために、factory 別にグループ化された切片の変量効果項を含めます。応答変数 defects はポアソン分布であり、このモデルの適切なリンク関数は対数です。係数の予測にラプラス近似メソッドを使用します。ダミー変数エンコードを 'effects' として指定すると、ダミー変数の係数の合計が 0 になります。

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

$${\text{defects}}_{ij}\sim \text{Poisson}({\mu }_{ij}).$$

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

$$\log ({\mu }_{ij})={\beta }_0+{\beta }_1{\text{newprocess}}_{ij}+{\beta }_2{\text{time}\text{\_}\text{dev}}_{ij}+{\beta }_3{\text{temp}\text{\_}\text{dev}}_{ij}+{\beta }_4{\text{supplier}\text{\_}\text{C}}_{ij}+{\beta }_5{\text{supplier}\text{\_}\text{B}}_{ij}+b_i,$$

ここで

• $${\text{defects}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ で実行されたバッチで観測された欠陥数です。

• $${\mu }_{ij}$$ は、バッチ $$j$$ ($$j=1,2,...,5$$) 処理中の工場 $$i$$ ($$i=1,2,...,20$$) に対応する欠陥の平均数です。

• $${\text{newprocess}}_{ij}$$、$${\text{time}\text{\_}\text{dev}}_{ij}$$ および $${\text{temp}\text{\_}\text{dev}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ に対応する各変数の測定値です。たとえば $${\text{newprocess}}_{ij}$$ は、工場 $$i$$ で実行されたバッチ $$j$$ 処理中に新プロセスが使用されたかどうかを示します。

• $${\text{supplier}\text{\_}\text{C}}_{ij}$$ および $${\text{supplier}\text{\_}\text{B}}_{ij}$$ はエフェクト (ゼロサム) コーディングを使用するダミー変数であり、バッチ $$j$$ 処理中に工場 $$i$$ で実行されたバッチに対して、それぞれ会社 C または B が加工化学薬品を供給したかどうかを示します。

• $$b_i\sim N(0,{\sigma }_b^2)$$ は、工場特有の品質変動に相当する、各工場 $$i$$ の変量効果の切片です。

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

固定効果の計画行列を抽出し、行 1 から 行 10 を表示します。

Dfe = designMatrix(glme,'Fixed');
disp(Dfe(1:10,:))

    1.0000         0    0.1834    0.2259    1.0000         0
    1.0000         0    0.3035    0.0725         0    1.0000
    1.0000         0    0.0717    0.1630    1.0000         0
    1.0000         0    0.1069    0.0809   -1.0000   -1.0000
    1.0000         0    0.0241    0.0319    1.0000         0
    1.0000         0    0.1214    0.1114         0    1.0000
    1.0000         0    0.0033    0.0553    1.0000         0
    1.0000         0    0.2350    0.0616    1.0000         0
    1.0000         0    0.0488    0.0177         0    1.0000
    1.0000         0    0.1148    0.0105    1.0000         0

固定効果の計画行列 Dfe の列 1 には定数項があります。列 2、3 および 4 にはそれぞれ newprocess、time_dev および temp_dev 項があります。列 5 と 6 にはそれぞれ supplier_C と supplier_B のダミー変数があります。

変量効果の計画行列を抽出し行 1 から 行 10 を表示します。

Dre = designMatrix(glme,'Random');
disp(Dre(1:10,:))

   (1,1)        1
   (2,1)        1
   (3,1)        1
   (4,1)        1
   (5,1)        1
   (6,2)        1
   (7,2)        1
   (8,2)        1
   (9,2)        1
  (10,2)        1

スパース行列 Dre を非スパース行列に変換し、行 1 から 行 10 を表示します。

full(Dre(1:10,:))

ans = 10×20

     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0

各列は、グループ化変数 factory の水準に対応します。

参考

GeneralizedLinearMixedModel | fitglme | fitted | residuals | response