Main Content

predict

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

説明

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

ypred = predict(lme,tblnew) は、新しいテーブルまたはデータセット配列 tblnew の値で近似線形混合効果モデル lme から条件付きの予測された応答 ypred のベクトルを返します。モデル lme の近似にテーブルまたはデータセット配列を使用している場合は predict のテーブルまたはデータセット配列を使用します。

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

ypred = predict(lme,Xnew,Znew) は、新しい固定効果および変量効果の計画行列 XnewZnew それぞれの値での近似線形混合効果モデル lme から条件付き予測応答 ypred のベクトルを返します。Znew は行列の cell 配列にすることもできます。この場合、グループ化変数 Gones(n,1) になります。n は近似で使用される観測値の数です。

モデル lme の近似に計画行列を使用している場合、predict に行列形式を使用します。

ypred = predict(lme,Xnew,Znew,Gnew) は、新しい固定効果および変量効果の計画行列 XnewZnew それぞれの値およびグループ化変数 Gnew の値での近似線形混合効果モデル lme から条件付き予測応答 ypred のベクトルを返します。

Znew および Gnew は、それぞれ行列とグループ化変数の cell 配列にすることもできます。

ypred = predict(___,Name,Value) は、1 つ以上の Name,Value のペアの引数で指定された追加のオプションにより、近似線形混合効果モデル lme から予測応答 ypred のベクトルを返します。

たとえば、信頼水準、同時信頼限界または固定効果のみからの寄与を指定できます。

[ypred,ypredCI] = predict(___) は、前の構文のいずれかの入力引数に対する予測 ypred の信頼区間 ypredCI も返します。

[ypred,ypredCI,DF] = predict(___) は、前の構文の任意の入力引数についての信頼区間の計算に使用される自由度 DFも返します。

すべて折りたたむ

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

load('fertilizer.mat');

このデータセット配列には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。

実用目的でこのデータを ds という名前のデータセット配列に保存し、TomatoSoil および Fertilizer をカテゴリカル変数として定義します。

ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);

線形混合効果モデルを当てはめます。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。

lme = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)');

元の計画値で応答値を予測します。観測された応答値をもつ最初の 5 つの観測値を表示します。

yhat = predict(lme);
[yhat(1:5) ds.Yield(1:5)]
ans = 5×2

  115.4788  104.0000
  135.1455  136.0000
  152.8121  158.0000
  160.4788  174.0000
   58.0839   57.0000

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

load carsmall

Weight に対する固定効果と、Model_Year でグループ化されたランダム切片で、線形混合効果モデルを当てはめます。まず、データをテーブルに保存します。

tbl = table(MPG,Weight,Model_Year);
lme = fitlme(tbl,'MPG ~ Weight + (1|Model_Year)');

データに対する予測された応答を作成します。

yhat = predict(lme,tbl);

元の応答と予測された応答をプロットして、相違点を確認します。それらをモデル年別にグループ化します。

figure()
gscatter(Weight,MPG,Model_Year)
hold on
gscatter(Weight,yhat,Model_Year,[],'o+x')
legend('70-data','76-data','82-data','70-pred','76-pred','82-pred')
hold off

Figure contains an axes object. The axes object with xlabel Weight, ylabel yhat contains 6 objects of type line. One or more of the lines displays its values using only markers These objects represent 70-data, 76-data, 82-data, 70-pred, 76-pred, 82-pred.

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

load('fertilizer.mat');

このデータセット配列には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。

実用目的でこのデータを ds という名前のデータセット配列に保存し、TomatoSoil および Fertilizer をカテゴリカル変数として定義します。

ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);

線形混合効果モデルを当てはめます。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。

lme = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)');

計画値をもつ新しいデータセット配列を作成します。新しいデータセット配列は、モデル lme の近似に使用する元のデータセット配列と同じ変数をもたなければなりません。

dsnew = dataset();
dsnew.Soil = nominal({'Sandy';'Silty'});
dsnew.Tomato = nominal({'Cherry';'Vine'});
dsnew. Fertilizer = nominal([2;2]);

元の計画点の条件付き応答と限界応答を予測します。

yhatC = predict(lme,dsnew);
yhatM = predict(lme,dsnew,'Conditional',false);
[yhatC yhatM]
ans = 2×2

   92.7505  111.6667
   87.5891   82.6667

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

load carbig

ガロンあたりの走行マイル数 (MPG) の線形混合効果モデルを当てはめます。加速度、馬力、気筒数は固定効果で、モデル年によってグループ化される切片と加速度については相関された変量効果の可能性があります。

最初に、線形混合効果モデルを当てはめるための計画行列を準備します。

X = [ones(406,1) Acceleration Horsepower];
Z = [ones(406,1) Acceleration];
Model_Year = nominal(Model_Year);
G = Model_Year;

次に、定義した計画行列とグループ化変数で fitlmematrix を使用してモデルを当てはめます。

lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Horsepower'},'RandomEffectPredictors',...
{{'Intercept','Acceleration'}},'RandomEffectGroups',{'Model_Year'});

応答値を予測する対象データを含む計画行列を作成します。XnewX のように 3 つの列をもたなければなりません。最初の列は 1 の列でなければなりません。また、最後の 2 列の値はそれぞれ AccelerationHorsepower に一致しなければなりません。Znew の 1 列目は 1 の列でなければならず、2 列目には Xnew と同じ Acceleration 値が含まれなければなりません。G の元のグループ化変数はモデル年です。そのため、Gnew にはモデル年の値が含まれなければなりません。Gnew にはノミナル値が含まれていなければなりません。

Xnew = [1,13.5,185; 1,17,205; 1,21.2,193];
Znew = [1,13.5; 1,17; 1,21.2]; % alternatively Znew = Xnew(:,1:2);
Gnew = nominal([73 77 82]);

新しい計画行列のデータに対する応答を予測します。

yhat = predict(lme,Xnew,Znew,Gnew)
yhat = 3×1

    8.7063
    5.4423
   12.5384

次に、切片と加速度について変量効果が無相関な線形混合モデルに、同じ手順を繰り返します。最初に、元の変量効果計画と変量効果グループ化変数を変更します。次に、モデルを再度当てはめます。

Z = {ones(406,1),Acceleration};
G = {Model_Year,Model_Year};

lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Horsepower'},'RandomEffectPredictors',...
{{'Intercept'},{'Acceleration'}},'RandomEffectGroups',{'Model_Year','Model_Year'});

ここで、新しい変量効果計画 Znew およびグループ化変数計画 Gnew を再作成します。これらのどちらかを使用して応答値を予測します。

Znew = {[1;1;1],[13.5;17;21.2]};
MY = nominal([73 77 82]);
Gnew = {MY,MY};

新しい計画行列を使用して応答を予測します。

yhat = predict(lme,Xnew,Znew,Gnew)
yhat = 3×1

    8.6365
    5.9199
   12.1247

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

load carbig

ガロンあたりの走行マイル数 (MPG) の線形混合効果モデルを当てはめます。加速度、馬力、気筒数は固定効果で、モデル年によってグループ化される切片と加速度については相関された変量効果の可能性があります。まず、変数をテーブルに格納します。

tbl = table(MPG,Acceleration,Horsepower,Model_Year);

次に、定義した計画行列とグループ化変数で fitlme を使用してモデルを当てはめます。

lme = fitlme(tbl,'MPG ~ Acceleration + Horsepower + (Acceleration|Model_Year)');

新しいデータを作成し、新しいテーブルに格納します。

tblnew = table();
tblnew.Acceleration = linspace(8,25)';
tblnew.Horsepower = linspace(min(Horsepower),max(Horsepower))';
tblnew.Model_Year = repmat(70,100,1);

linspace は、入力値の下限と上限の間で 100 個の等間隔の値を作成します。Model_Year は 70 で固定されています。これは、任意のモデル年について繰り返すことができます。

予測値および 95% の信頼限界 (非同時) を計算およびプロットします。

[ypred,yCI,DF] = predict(lme,tblnew);
figure(); 
h1 = line(tblnew.Acceleration,ypred);
hold on;
h2 = plot(tblnew.Acceleration,yCI,'g-.');

Figure contains an axes object. The axes object contains 3 objects of type line.

自由度を表示します。

DF(1)
ans = 389

同時信頼限界を計算およびプロットします。

[ypred,yCI,DF] = predict(lme,tblnew,'Simultaneous',true);
h3 = plot(tblnew.Acceleration,yCI,'r--');

Figure contains an axes object. The axes object contains 5 objects of type line.

自由度を表示します。

DF
DF = 389

サタースウェイト法を使用して同時信頼限界を計算し、自由度を計算します。

[ypred,yCI,DF] = predict(lme,tblnew,'Simultaneous',true,'DFMethod','satterthwaite');
h4 = plot(tblnew.Acceleration,yCI,'k:');
hold off
xlabel('Acceleration')
ylabel('Response')
ylim([-50,60])
xlim([8,25])
legend([h1,h2(1),h3(1),h4(1)],'Predicted response','95%','95% Sim',...
'95% Sim-Satt','Location','Best')

Figure contains an axes object. The axes object with xlabel Acceleration, ylabel Response contains 7 objects of type line. These objects represent Predicted response, 95%, 95% Sim, 95% Sim-Satt.

自由度を表示します。

DF
DF = 3.6001

入力引数

すべて折りたたむ

線形混合効果モデル。fitlme または fitlmematrix を使用して構築した LinearMixedModel オブジェクトとして指定します。

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

n 行 p 列の行列で指定される新しい固定効果の計画行列。ここで n は観測値の数、p は固定予測子変数の数です。X の各行は 1 つの観測値に対応し、X の各列は 1 つの変数に対応します。

データ型: single | double

n 行 q 列の行列または R 計画行列 Z{r} の cell 配列として指定される新しい変量効果計画。ここで、r = 1, 2, ..., R となります。Znew が cell 配列の場合、各 Z{r} は n 行 q(r) 列の行列になります。ここで n は観測値の数、q(r) は無作為な予測子変数の数です。

データ型: single | double | cell

新しいグループ化変数または変数。線形混合効果モデル lme の当てはめには、元のグループ化変数と同じ水準またはグループをもつグループ化変数の長さ R のベクトルまたは cell 配列として指定します。

データ型: single | double | categorical | logical | char | string | cell

名前と値の引数

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

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

例: ypred = predict(lme,'Conditional',false);

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

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

例: 'Alpha',0.01

データ型: single | double

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

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

例: 'Conditional,false

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

'residual'既定の設定。自由度は定数で n – p に等しいと仮定されます。ここで n は観測値の数、p は固定効果の数です。
'satterthwaite'サタースウェイトの近似法。
'none'すべての自由度は無限大に設定されます。

たとえば、次のようにサタースウェイトの近似法を指定できます。

例: 'DFMethod','satterthwaite'

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

false既定の設定。非同期区間。
true同期区間。

例: 'Simultaneous',true

予測のタイプ。'Prediction' と以下のいずれかで構成される、コンマ区切りペアとして指定します。

'curve'既定の設定。近似関数に基づく予測の信頼限界です。
'observation'新しい観測値の観測誤差に起因する変動性も、信頼限界の計算に含まれるので、区間が広くなります。

例: 'Prediction','observation'

出力引数

すべて折りたたむ

ベクトルとして返される予測応答。ypred には、'Conditional' の名前と値のペア引数の値の選択に応じて、条件付き応答または限界応答を含むことができます。条件付き予測には、固定効果と変量効果の両方からの寄与が含まれます。

2 列の行列として返される予測値の点別信頼区間。yCI の 1 列目には信頼区間の下限が含まれ、2 列目には上限が含まれます。既定では、yCI には予測の 95% の信頼区間が含まれます。Alpha の名前と値のペア引数を使用して信頼水準を変更したり、Simultaneous の名前と値のペア引数を使用して信頼水準を同時にしたり、Prediction の名前と値のペア引数を使用して、曲線ではなく新しい観測値用にすることもできます。

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

  • 'Simultaneous' の名前と値のペア引数が false である場合、DF はベクトルになります。

  • 'Simultaneous' の名前と値のペア引数が true である場合、DF はスカラー値になります。

詳細

すべて折りたたむ

条件付き予測と限界予測

条件付き予測には固定効果と変量効果の両方からの寄与が含まれますが、限界モデルには固定効果のみからの寄与が含まれます。

線形混合効果モデル lme には、n 行 p 列の固定効果の計画行列 X と、n 行 q 列の変量効果の計画行列 Z があるものとします。また、推定した固定効果が含まれている p 行 1 列のベクトルが β^、変量効果について推定した最良線形不偏予測量 (BLUP) が含まれている q 行 1 列のベクトルが b^ であるとします。予測される条件付き応答は次のようになります。

y^Cond=Xβ^+Zb^,

これは、名前と値のペアの引数 'Conditional','true' に対応します。

予測される限界応答は次のようになります。

y^Mar=Xβ^,

これは、名前と値のペアの引数 'Conditional','false' に対応します。

予測を行う場合、特定のグループ化変数に新しい水準 (元のデータ内に存在していない 1) がある場合、そのグループ化変数に対する変量効果は、グループ化変数に新しい水準が含まれる観測値での 'Conditional' 予測に寄与しません。

バージョン履歴

R2013b で導入