random

クラス: GeneralizedLinearMixedModel

当てはめた一般化線形混合効果モデルからのランダム応答の生成

構文

ysim = random(glme)

ysim = random(glme,tblnew)

ysim = random(___,Name,Value)

説明

ysim = random(glme) は、元の計画点で当てはめた一般化線形混合効果モデル glme から、シミュレートされた応答 ysim を返します。

ysim = random(glme,tblnew) は、テーブルまたはデータセット配列 tblnew で指定された新規の入力値を使用して、シミュレートされた応答を返します。

ysim = random(___,Name,Value) は、前述の構文のいずれかを使用する 1 つ以上の Name,Value のペアの引数により指定された追加オプションによって、シミュレートされた応答を返します。たとえば、観測値の重み、二項分布のサイズまたはモデルのオフセットを指定できます。

入力引数

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

GeneralizedLinearMixedModel オブジェクト

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

tblnew — 新しい入力データ

table | データセット配列

応答変数、予測変数およびグループ化変数が含まれる新規入力データ。table またはデータセット配列として指定します。予測変数は連続変数またはグループ化変数にすることができます。tblnew は、一般化線形混合効果モデル glme の当てはめに使用される元のテーブルまたはデータセット配列 tbl と同じ変数をもっていなければなりません。

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

BinomialSize — 二項分布の試行回数

ones(m,1) (既定値) | 正の整数値の m 行 1 列のベクトル

二項分布の試行回数。'BinomialSize' と正の整数値の m 行 1 列のベクトルから構成されるコンマ区切りペアとして指定します。ここで、m は tblnew の行の数になります。'BinomialSize' 名前と値のペアは、二項分布にのみ適用します。値は、ランダム応答値を生成する場合に二項の試行数を指定します。

データ型: single | double

Offset — モデル オフセット

zeros(m,1) (既定値) | スカラー値のベクトル

モデル オフセット。長さ m のスカラー値のベクトルとして指定されます。ここで、m は tblnew の行の数です。オフセットは追加の予測子として使用され、1 で固定した係数値をもちます。

Weights — 観測値の重み

非負のスカラー値の m 行 1 列のベクトル

観測値の重み。'Weights' と、非負のスカラー値の m 行 1 列のベクトルで構成されるコンマ区切りのペアとして指定します。ここで、m は、tblnew の行数です。応答の分布が二項またはポアソンである場合、'Weights' は正の整数のベクトルでなければなりません。

データ型: single | double

出力引数

ysim — シミュレートされた応答値

m 行 1 列のベクトル

シミュレートされた応答値。m 行 1 列のベクトルとして返され、m は tblnew の行数です。random は、当てはめた事前分布に基づいて変量効果のベクトルを最初に生成することにより ysim を作成します。次に、random は変量効果が与えられた近似条件付き分布から ysim を生成します。random は、fitglme を使用してモデルを当てはめるとき、指定された観測値の重みの効果があれば、これを考慮します。

例

GLME モデルからのランダム応答のシミュレート

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

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');

random を使用して、当てはめたモデルから新しい応答ベクトルをシミュレートします。

rng(0,'twister');  % For reproducibility
ynew = random(glme);

シミュレートされた応答ベクトルの最初の 10 行を表示します。

ynew(1:10)

ans = 10×1

     3
     3
     1
     7
     5
     8
     7
     9
     5
     9

新しい入力値を使用して、新しい応答ベクトルをシミュレートします。mfr の最初の 10 件を tblnew にコピーすることにより新規テーブルを作成します。

tblnew = mfr(1:10,:);

mfr の先頭 10 行には、工場 1 および 2 の検査 1 ～ 5 から収集したデータが含まれています。どちらの工場でも実験時のすべての検査に古いプロセスを使用したので、10 個の観測値はすべて newprocess = 0 です。

newprocess の値を tblnew の観測値の 1 に変更します。

tblnew.newprocess = ones(height(tblnew),1);

tblnew の新しい入力値を使用して、新しい応答をシミュレートします。

ynew2 = random(glme,tblnew)

ynew2 = 10×1

     2
     3
     5
     4
     2
     2
     2
     1
     2
     0

詳細

条件付き分布メソッド

random は、以下のように当てはめた一般化線形混合効果モデルから乱数データを生成します。

• $$b_{sim}\sim P\left(b|\widehat{\theta },{\widehat{\sigma }}^2\right)$$ をサンプリングします。ここで、$$P\left(b|\widehat{\theta },{\widehat{\sigma }}^2\right)$$ は変量効果について推定した事前分布、$$\widehat{\theta }$$ は推定した共分散パラメーターのベクトル、$${\widehat{\sigma }}^2$$ は推定した分散パラメーターです。

• i= 1 ～ m について与えられた b_sim に対して $$y_{sim\_i}\sim P\left(y_{new\_i}|b_{sim},\widehat{\beta },\widehat{\theta },{\widehat{\sigma }}^2\right)$$ をサンプリングします。ここで $$P\left(y_{new\_i}|b_{sim},\widehat{\beta },\widehat{\theta },{\widehat{\sigma }}^2\right)$$ は、与えられた b_sim およびモデル パラメーターに対する i 番目の新しい応答 y_{new_i} の条件付き分布です。

参考

GeneralizedLinearMixedModel | fitglme | predict | fitted