predict

クラス: GeneralizedLinearMixedModel

一般化線形混合効果モデルの応答の予測

構文

ypred = predict(glme)

ypred = predict(glme,tblnew)

ypred = predict(___,Name,Value)

[ypred,ypredCI] = predict(___)

[ypred,ypredCI,DF] = predict(___)

説明

ypred = predict(glme) は、一般化線形混交効果モデル glme の当てはめに使われる元の予測子を使用して、応答の予測された条件付き平均 ypred を返します。

ypred = predict(glme,tblnew) は、tblnew で指定された新たな予測子の値を使用して、予測した条件付き平均を返します。

tblnew のグループ化変数が、元のデータにない水準を保有する場合、そのグループ化変数の変量効果は、グループ化変数に新しい水準が含まれる観測値での 'Conditional' 予測に役立ちません。

ypred = predict(___,Name,Value) は、1 つ以上の Name,Value のペアの引数に指定された追加オプションを使用して、応答の予測した条件付き平均を返します。たとえば、信頼水準、同時信頼限界または固定効果のみからの寄与を指定できます。前の構文の入力引数のいずれかを使用できます。

[ypred,ypredCI] = predict(___) はまた、各予測値の 95% の点別信頼区間 ypredCI を返します。

[ypred,ypredCI,DF] = predict(___) はまた、信頼区間の計算に使用される自由度 DF を返します。

入力引数

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

GeneralizedLinearMixedModel オブジェクト

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

tblnew — 新しい入力データ

テーブル | dataset 配列

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

名前と値の引数

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

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

Alpha — 有意水準

0.05 (既定値) | [0,1] の範囲のスカラー値

有意水準。'Alpha' と [0,1] の範囲にあるスカラー値から構成されるコンマ区切りのペアとして指定します。値が α の場合、信頼水準は 100 × (1 – α)% です。

たとえば、99% の信頼区間の場合は、次のように信頼水準を指定できます。

例: 'Alpha',0.01

データ型: single | double

Conditional — 条件付き予測のインジケーター

true (既定値) | false

条件付き予測のインジケーター。'Conditional' と、以下のいずれかで構成されるコンマ区切りペアとして指定されます。

値	説明
true	固定効果と変量効果の両方からの寄与 (条件付き)
false	固定効果のみからの寄与 (限界)

例: 'Conditional',false

DFMethod — 自由度の近似の計算方法

'residual' (既定値) | 'none'

自由度の近似の計算方法。'DFMethod' と次のいずれかの値で構成されるコンマ区切りのペアとして指定します。

値	説明
'residual'	自由度の値は n - p に等しい定数であると仮定されます。n は観測値の個数、p は固定効果の個数です。
'none'	自由度は無限大に設定されます。

例: 'DFMethod','none'

Offset — モデル オフセット

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

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

Simultaneous — 信頼限界のタイプ

false (既定値) | true

信頼限界のタイプ。'Simultaneous' と false または true のいずれかで構成されるコンマ区切りペアとして指定します。

• 'Simultaneous' が false の場合、predict は非同時信頼限界を計算します。

• 'Simultaneous' が true の場合、predict は同時信頼限界を返します。

例: 'Simultaneous',true

出力引数

ypred — 予測応答

ベクトル

予測応答。ベクトルとして返します。'Conditional' 名前と値のペアの引数が true として指定される場合、ypred には、変量効果が与えられた応答の条件付き平均の予測が含まれます。条件付き予測には、固定効果と変量効果の両方からの寄与が含まれます。限界予測には、固定効果の寄与のみが含まれます。

限界予測を計算するには、predict は条件付き予測を計算しますが、変量効果の EBP (経験的ベイズ予測) の代わりに、ゼロのベクトルを代入します。

ypredCI — 点別信頼区間

2 列の行列

2 列の行列として返される予測値の点別信頼区間。ypredCI の 1 列目には信頼区間の下限が含まれ、2 列目には上限が含まれます。既定の設定では、ypredCI には予測の 95% の非同時信頼区間が含まれます。Alpha 名前と値のペアの引数を使用して信頼水準を変更し、Simultaneous 名前と値のペアの引数を使用して、同時信頼水準にすることができます。

fitglme と最大尤度の近似メソッド ('Laplace' または 'ApproximateLaplace') のいずれかを使用して GLME モデルを当てはめる場合、predict は、推定の共分散パラメーターおよび観測された応答を条件とする CMSEP (予測の条件付き平均二乗誤差) メソッドを使用して、信頼区間を計算します。あるいは、信頼区間とは、推定された共分散パラメーターと観測された応答を条件とする近似のベイズの信頼できる区間として解釈することもできます。

fitglme と疑似尤度の近似メソッド ('MPL' または 'REMPL') の1 つを使用して GLME モデルを当てはめる場合、predict は疑似尤度の最後の反復からの固定線形混合効果モデル計算に基づきます。

DF — 自由度

ベクトル | スカラー値

信頼区間の計算に使用される自由度。ベクトルまたはスカラー値として返されます。

• 'Simultaneous' が false の場合、DF はベクトルです。

• 'Simultaneous' が true の場合、DF はスカラー値です。

例

元の計画値での応答の予測

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

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

元の計画値で応答値を予測します。観測された応答値と共に最初の 10 件の予測を表示します。

ypred = predict(glme);
[ypred(1:10),mfr.defects(1:10)]

ans = 10×2

    4.9883    6.0000
    5.9423    7.0000
    5.1318    6.0000
    5.6295    5.0000
    5.3499    6.0000
    5.2134    5.0000
    4.6430    4.0000
    4.5342    4.0000
    5.3903    9.0000
    4.6529    4.0000

列 1 には、元の計画値における予測応答値が含まれます。列 2 には、観測された応答値が含まれます。

新規テーブルの値での応答の予測

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

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

元の計画値で応答値を予測します。

ypred = predict(glme);

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 を使用して、予測応答値と 99% の非同時信頼区間を計算します。tblnew に基づいた予測値、mfr に基づいた予測値および観測された応答値の最初の 10 行を表示します。

[ypred_new,ypredCI] = predict(glme,tblnew,'Alpha',0.01);
[ypred_new,ypred(1:10),mfr.defects(1:10)]

ans = 10×3

    3.4536    4.9883    6.0000
    4.1142    5.9423    7.0000
    3.5530    5.1318    6.0000
    3.8976    5.6295    5.0000
    3.7040    5.3499    6.0000
    3.6095    5.2134    5.0000
    3.2146    4.6430    4.0000
    3.1393    4.5342    4.0000
    3.7320    5.3903    9.0000
    3.2214    4.6529    4.0000

列 1 には、tblnew のデータ (ここでは newprocess = 1) に基づいた予測応答値が含まれます。列 2 には、mfr の元データに基づいた予測応答値が含まれています。ここで、newprocess = 0 です。列 3 には、mfr の観測された応答値が含まれます。この結果に基づいて、他のすべての予測子が元の値を保持する場合、予測された欠陥の数は新しいプロセスの使用時により少なく見えます。

新しい予測応答値に対応する行 1 から行 10 の 99% の信頼区間を表示します。

ypredCI(1:10,1:2)

ans = 10×2

    1.6983    7.0235
    1.9191    8.8201
    1.8735    6.7380
    2.0149    7.5395
    1.9034    7.2079
    1.8918    6.8871
    1.6776    6.1597
    1.5404    6.3976
    1.9574    7.1154
    1.6892    6.1436

参照

[1] Booth, J.G., and J.P. Hobert. “Standard Errors of Prediction in Generalized Linear Mixed Models.” Journal of the American Statistical Association, Vol. 93, 1998, pp. 262–272.

参考

GeneralizedLinearMixedModel | fitglme | fitted | random