mvregress

多変量線形回帰

構文

beta = mvregress(X,Y)

beta = mvregress(X,Y,Name,Value)

[beta,Sigma] = mvregress(___)

[beta,Sigma,E,CovB,logL] = mvregress(___)

説明

beta = mvregress(X,Y) は、X の計画行列に対応する Y について、d 次元応答の多変量正規回帰の推定係数を返します。

beta = mvregress(X,Y,Name,Value) では、1 つ以上の名前と値のペア引数で指定された追加オプションを使用して推定係数を返します。たとえば、推定アルゴリズム、初期推定値または回帰の最大反復回数を指定できます。

[beta,Sigma] = mvregress(___) は、前の構文の入力引数のいずれかを使用して、d 行 d 列の Y の推定分散共分散行列も返します。

[beta,Sigma,E,CovB,logL] = mvregress(___) は、残差の行列 E、回帰係数の推定分散共分散行列 CovB、最後の反復後の対数尤度目的関数の値 logL も返します。

例

異なる切片をもつパネル データの多変量回帰モデル

多変量回帰モデルをパネル データに近似する場合、切片が異なっていても勾配は共通であると仮定します。

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

load('flu')

データセット配列 flu には、米国疾病対策予防センター (CDC) による米国全土のインフルエンザ発生件数の推定値と、Google® 検索データに基づく 9 つの地域別の推定値が含まれています。

応答データと予測子データを抽出します。

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

Y の応答は 9 つの地域別インフルエンザ発生件数の推定値です。1 年のすべての週に観測値が存在するので、$$n$$ = 52 です。応答の次元は地域に対応するので、$$d$$ = 9 です。x 内の予測子は、1 週間における全国のインフルエンザ発生件数の推定値です。

地域別にインフルエンザ データをプロットします。

figure;
regions = flu.Properties.VarNames(2:end-1);
plot(x,Y,'x')
legend(regions,'Location','NorthWest')

多変量回帰モデル $$y_{ij}={\alpha }_j+\beta x_{ij}+{ϵ}_{ij}$$ を当てはめます。ここで、$$i=1,\dots ,n$$ および $$j=1,\dots ,d$$ であり、地域間の同時相関は $$COV({ϵ}_{ij},{ϵ}_{ij})={\sigma }_{jj}$$ です。

$$K$$ = 10 個の回帰係数を推定します。9 つが切片の項、1 つが共通の勾配です。入力引数 X は、$$d$$ 行 $$K$$ 列の計画行列が含まれている $$n$$ 要素の cell 配列でなければなりません。

X = cell(n,1);
for i = 1:n
	X{i} = [eye(d) repmat(x(i),d,1)];
end
[beta,Sigma] = mvregress(X,Y);

beta には、$$K$$ 次元係数ベクトル $$({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_9,\beta {)}^{\prime }$$ の推定値が格納されます。

Sigma には、地域間の同時相関に関する $$d$$ 行 $$d$$ 列の分散共分散行列 $${({\sigma }_{ij})}_{d\times d}$$ ($$i,j=1,\dots ,d$$) の推定値が格納されます。

近似した回帰モデルをプロットします。

B = [beta(1:d)';repmat(beta(end),1,d)];
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;

figure;
h = plot(x,Y,'x',xx,fits,'-');
for i = 1:d
	set(h(d+i),'color',get(h(i),'color'));
end
legend(regions,'Location','NorthWest');

このプロットは、各回帰線の切片は異なっていますが、勾配が同じであることを示しています。視覚的には、いくつかの回帰線は他の回帰線よりもデータを適切に近似しているように見えます。

異なる勾配をもつパネル データの多変量回帰

切片と勾配が異なると仮定して、最小二乗法を使用して多変量回帰モデルをパネル データに近似します。

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

load('flu');

データセット配列 flu には、米国疾病対策予防センター (CDC) による米国全土のインフルエンザ発生件数の推定値と、Google® 検索に基づく 9 つの地域別の推定値が含まれています。

応答データと予測子データを抽出します。

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

Y の応答は 9 つの地域別インフルエンザ発生件数の推定値です。1 年のすべての週に観測値が存在するので、$$n$$ = 52 です。応答の次元は地域に対応するので、$$d$$ = 9 です。x 内の予測子は、1 週間における全国のインフルエンザ発生件数の推定値です。

多変量回帰モデル $$y_{ij}={\alpha }_j+{\beta }_j{x}_{ij}+{ϵ}_{ij}$$ を当てはめます。ここで、$$i=1,\dots ,n$$ および $$j=1,\dots ,d$$ であり、地域間の同時相関は $$COV({ϵ}_{ij},{ϵ}_{ij})={\sigma }_{jj}$$ です。

$$K$$ = 18 個の回帰係数を推定します。9 つが切片の項、9 つが勾配の項です。X は、$$d$$ 行 $$K$$ 列の計画行列が含まれている $$n$$ 要素の cell 配列です。

X = cell(n,1);
for i = 1:n
    X{i} = [eye(d) x(i)*eye(d)];
end
[beta,Sigma] = mvregress(X,Y,'algorithm','cwls');

beta には、$$K$$ 次元係数ベクトル $$({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_9,{\beta }_1,{\beta }_2,\dots ,{\beta }_9{)}^{\prime }$$ の推定値が格納されます。

近似した回帰モデルをプロットします。

B = [beta(1:d)';beta(d+1:end)'];
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;

figure;
h = plot(x,Y,'x',xx,fits,'-');
for i = 1:d
	set(h(d+i),'color',get(h(i),'color'));
end

regions = flu.Properties.VarNames(2:end-1);
legend(regions,'Location','NorthWest');

このプロットは、各回帰線は切片と勾配が異なっていることを示しています。

単一の計画行列を使用した多変量回帰

すべての応答次元について、単一の $$n$$ 行 $$P$$ 列の計画行列を使用して多変量回帰モデルを当てはめます。

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

load('flu')

データセット配列 flu には、米国疾病対策予防センター (CDC) による米国全土のインフルエンザ発生件数の推定値と、Google® 検索に基づく 9 つの地域別の推定値が含まれています。

応答データと予測子データを抽出します。

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

Y の応答は 9 つの地域別インフルエンザ発生件数の推定値です。1 年のすべての週に観測値が存在するので、$$n$$ = 52 です。応答の次元は地域に対応するので、$$d$$ = 9 です。x 内の予測子は、1 週間における全国のインフルエンザ発生件数の推定値です。

$$n$$ 行 $$P$$ 列の計画行列 X を作成します。定数項を回帰に含めるため、1 の列を追加します。

X = [ones(size(x)),x];

多変量回帰モデルを当てはめます。

$$y_{ij}={\alpha }_j+{\beta }_j{x}_{ij}+{ϵ}_{ij},$$

ここで、$$i=1,\dots ,n$$ および $$j=1,\dots ,d$$ であり、地域間の同時相関は次のようになります。

$$COV({ϵ}_{ij},{ϵ}_{ij})={\sigma }_{jj}.$$

推定する回帰係数は 18 個で、9 個は切片の項、9 個は勾配の項です。

[beta,Sigma,E,CovB,logL] = mvregress(X,Y);

beta には $$P$$ 行 $$d$$ 列の係数行列の推定値が格納されます。Sigma には地域間の同時相関に関する $$d$$ 行 $$d$$ 列の分散共分散行列の推定値が格納されます。E は残差の行列です。CovB は回帰係数の推定分散共分散行列です。logL は最後の反復後の対数尤度目的関数の値です。

近似した回帰モデルをプロットします。

B = beta;
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;

figure
h = plot(x,Y,'x', xx,fits,'-');
for i = 1:d
    set(h(d+i),'color',get(h(i),'color'))
end

regions = flu.Properties.VarNames(2:end-1);
legend(regions,'Location','NorthWest')

このプロットは、各回帰線は切片と勾配が異なっていることを示しています。

入力引数

X — 計画行列

行列 | 行列の cell 配列

多変量回帰の計画行列。行列または行列の cell 配列として指定します。ここで、n はデータ内の観測の数、K は推定する回帰係数の数、p は予測子変数の数、d は応答変数行列 Y 内の次元の数です。

• d = 1 である場合、X を 1 つの n 行 K 列の計画行列として指定します。

• d > 1 でありすべての d 次元が同じ計画行列をもつ場合は、X を、1 つの n 行 p 列の計画行列 (cell 配列内に含まれていない) として指定できます。

• d > 1 でありすべての n 観測が同じ計画行列をもつ場合は、X を、1 つの d 行 K 列の計画行列を含む cell 配列として指定できます。

• d > 1 でありすべての n 観測が同じ計画行列をもつ場合は、X を、d 行 K 列の計画行列を含む長さ n の cell 配列として指定します。

回帰モデルに定数項を含めるには、各計画行列に 1 の列を 1 つ含めなければなりません。

mvregress は、X の NaN 値を欠損値として扱い、欠損値を含む X の行を無視します。

データ型: single | double | cell

Y — 応答変数

行列

n 行 d 列の行列として指定される応答変数。n はデータ内の観測値の数、d は応答の次元の数です。d = 1 である場合、mvregress は Y の値を n 個の独立応答値のように扱います。

mvregress は、Y の NaN 値を欠損値として扱い、名前と値のペア引数 algorithm で指定された推定アルゴリズムに従って処理します。

データ型: single | double

名前と値の引数

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

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

例: 'algorithm','cwls','covar0',C は、共分散行列 C を使用して、共分散の重み付き最小二乗推定法を使用するよう指定します。

algorithm — 推定アルゴリズム

'mvn' | 'ecm' | 'cwls'

'algorithm' および以下のいずれかで構成されるコンマ区切りペアとして指定される推定アルゴリズム。

'mvn'	標準の多変量正規最尤推定法。
'ecm'	ECM アルゴリズムを用いた最尤推定法。
'cwls'	共分散の重み付き最小二乗推定法。

既定のアルゴリズムは、欠損値の有無によって異なります。

• 完全なデータの場合、既定の設定は 'mvn' です。

• 欠損値 (つまり NaN) がある場合、標本のサイズが全パラメーターの推定に十分であれば、既定の設定は 'ecm' になります。それ以外の場合、既定のアルゴリズムは 'cwls' です。

メモ

algorithm の値が 'mvn' の場合、mvregress は推定を行う前に、欠落した応答値をもつ観測値を削除します。

例: 'algorithm','ecm'

beta0 — 回帰係数の初期推定

ベクトル

回帰係数の初期推定。'beta0' および K 個の要素のベクトルで構成されるコンマ区切りのペアとして指定します。既定値は 0 のベクトルです。

引数 beta0 は、推定アルゴリズム algorithm が 'mvn' の場合は使用されません。

covar0 — 分散共分散行列の初期推定

行列

分散共分散行列 Sigma の初期推定。コンマ区切りのペアとして指定し、'covar0' および対称の d 行 d 列の正定値行列で構成されます。既定の値は単位行列です。

推定アルゴリズム algorithm が 'cwls' の場合、mvregress は各反復において covar0 を重み行列としてその値を変えずに使用します。

covtype — 分散共分散行列のタイプ

'full' (既定値) | 'diagonal'

Y について推定する分散共分散行列のタイプ。'covtype' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

'full'	すべての d(d + 1)/2 分散共分散要素を推定します。
'diagonal'	分散共分散行列の d 対角要素のみを推定します。

例: 'covtype','diagonal'

maxiter — 最大反復回数

100 (既定値) | 正の整数

推定アルゴリズムの最大反復回数。'maxiter' と正の整数で構成されるコンマ区切りのペアとして指定します。

推定が tolbeta および tolobj の収束許容誤差内になるか、maxiter で指定された最大反復回数に達するまで、反復が繰り返されます。tolbeta と tolobj がともに 0 の場合、mvregress は収束テストなしで maxiter 回の反復を実行します。

例: 'maxiter',50

outputfcn — 各反復で評価を行う関数

関数ハンドル

各反復において評価を行う関数。'outputfcn' と関数ハンドルで構成されるコンマ区切りのペアとして指定します。関数は true または false の論理値を返さなければなりません。mvregress は各反復においてこの関数の評価を行います。結果が true の場合、反復は停止します。それ以外の場合、反復は継続されます。たとえば、現在の反復結果をプロットまたは表示して、その図が閉じられると true を返すような関数を指定できます。

この関数は、以下の順序で 3 つの入力引数を受け入れなければなりません。

• 現在の係数推定値のベクトル

• 以下の 3 つのフィールドをもつ構造体:

  Covar	分散共分散行列の現在の値
  iteration	現在の反復番号
  fval	対数尤度目的関数の現在の値

• 次の 3 つの値をとるテキスト。

  'init'	初期化中に関数が呼び出された場合
  'iter'	反復の後で関数が呼び出された場合
  'done'	完了後に関数が呼び出された場合

tolbeta — 回帰係数の収束許容誤差

sqrt(eps) (既定値) | 正のスカラー値

回帰係数の収束許容誤差。'tolbeta' と正のスカラー値で構成されるコンマ区切りのペアとして指定します。

$$b^t$$ が反復 t における係数ベクトルの推定値を表し、$${\tau }_{\beta }$$ が tolbeta で指定された許容誤差であるとします。回帰係数推定の収束基準は次のようになります。

$$\Vert b^t-b^{t-1}\Vert <{\tau }_{\beta }\sqrt{K}\left(1+\Vert b^t\Vert \right),$$

ここで、K は $$b^t$$ の長さ、$$\Vert v\Vert$$ はベクトル $$v.$$ のノルムです。

推定が tolbeta および tolobj の収束許容誤差内になるか、maxiter で指定された最大反復回数に達するまで、反復が繰り返されます。tolbeta と tolobj がともに 0 の場合、mvregress は収束テストなしで maxiter 回の反復を実行します。

例: 'tolbeta',1e-5

tolobj — 対数尤度目的関数の収束許容誤差

eps^(3/4) (既定値) | 正のスカラー値

対数尤度目的関数の収束許容誤差。'tolobj' と正のスカラー値で構成されるコンマ区切りのペアとして指定します。

$$L^t$$ が反復 t における対数尤度目的関数の値を表し、$${\tau }_{\ell }$$ が tolobj で指定された許容誤差であるとします。目的関数の収束基準は次のようになります。

$$\left|L^t-L^{t-1}\right|<{\tau }_{\ell }\left(1+\left|L^t\right|\right).$$

推定が tolbeta および tolobj の収束許容誤差内になるか、maxiter で指定された最大反復回数に達するまで、反復が繰り返されます。tolbeta と tolobj がともに 0 の場合、mvregress は収束テストなしで maxiter 回の反復を実行します。

例: 'tolobj',1e-5

varformat — 分散共分散行列のパラメーター推定の形式

'beta' (既定値) | 'full'

パラメーター推定の分散共分散行列 CovB の形式。コンマ区切りのペアとして指定し、'varformat' および以下のいずれか 1 つで構成されます。

'beta'	回帰係数の推定 beta のみの分散共分散行列を返します。
'full'	回帰係数の推定 beta および分散共分散行列の推定 Sigma の両方の分散共分散行列を返します。

例: 'varformat','full'

vartype — パラメーター推定の分散共分散行列のタイプ

'hessian' (既定値) | 'fisher'

パラメーター推定のための分散共分散行列のタイプ。'vartype' と、'hessian' または 'fisher' のいずれかで構成される、コンマ区切りのペアとして指定します。

• 値が 'hessian' である場合、mvregress はヘッシアン、つまり観測された情報を使用して CovB を計算します。

• 値が 'fisher' である場合、mvregress は完全なデータのフィッシャー行列、つまり予想される情報を使用して CovB を計算します。

'hessian' メソッドでは、データの欠損に伴う不確実性の増大を考慮しますが、'fisher' メソッドではこれが考慮されません。

例: 'vartype','fisher'

出力引数

beta — 推定された回帰係数

列ベクトル | 行列

列ベクトルまたは行列として返される、推定回帰係数。

• X を 1 つの n 行 K 列の計画行列として指定した場合、mvregress は beta を長さ K の列ベクトルとして返します。たとえば、X が 20 行 5 列の計画行列の場合、beta は 5 行 1 列の列ベクトルです。

• X を 1 つ以上の d 行 K 列の計画行列を含む cell 配列として指定した場合、mvregress は beta を長さ K の列ベクトルとして返します。たとえば、X が 2 行 10 列の計画行列を含む cell 配列の場合、beta は 10 行 1 列の列ベクトルです。

• X を 1 つの n 行 p 列の計画行列 (cell 配列内に含まれていない) として指定し、Y が次元 d > 1 をもつ場合、mvregress は beta を p 行 d 列の行列として返します。たとえば、X が 20 行 5 列の計画行列で、Y が d = 2 となるような 2 つの次元をもつ場合、beta は 5 行 2 列の行列で、近似された Y 値は X × beta です。

Sigma — 推定分散共分散行列

正方行列

Y の応答の推定された分散共分散行列。d 行 d 列の正方行列として返されます。

メモ

推定分散共分散行列 Sigma は、残差行列 E の標本共分散行列ではありません。

E — 残差

行列

近似された回帰モデルの残差。n 行 d 列の行列として返されます。

algorithm の値が 'ecm' または 'cwls' の場合、mvregress は Y の欠損値に対応する残差を条件付き補定値と当てはめた値との差として計算します。

メモ

algorithm の値が 'mvn' の場合、mvregress は推定を行う前に、欠落した応答値をもつ観測値を削除します。

CovB — パラメーター推定分散共分散行列

正方行列

パラメーター推定の分散共分散行列。正方行列として返します。

• varformat の値が 'beta' (既定の設定) である場合、CovB は beta にある係数推定の推定された分散共分散行列になります。

• varformat の値が 'full' である場合、CovB は beta と Sigma にある推定の組み合わせの推定分散共分散行列になります。

logL — 対数尤度目的関数の値

スカラー値

スカラー値として返される、最後の反復後の対数尤度目的関数の値。

詳細

多変量正規回帰

多変量正規回帰は、正規分布の誤差がある、予測子変数の計画行列に対する d 次元応答の回帰です。誤差は不等分散である可能性および相関をもつ可能性があります。

モデルは次のようになります。

$$y_i=X_i\beta +e_i,i=1,\dots ,n,$$

ここで

• $$y_i$$ は、応答が含まれている d 次元ベクトルです。

• $$X_i$$ は、予測子変数の計画行列です。

• $$\beta$$ は、回帰係数が含まれているベクトルまたは行列です。

• $$e_i$$ は、誤差項が含まれている d 次元ベクトルです。次のような多変量正規分布になります。

  $$e_i~MV{N}_d(0,\Sigma ).$$

条件付き補定値

期待値/条件付き最大化 ('ecm') および共分散の重み付き最小二乗 ('cwls') の推定アルゴリズムには、欠けている応答値の補定が含まれています。

$$\tilde{y}$$ が欠損観測値を表すとします。条件付き補定値は、観測データが与えられた場合の欠損観測値の期待値 $${\rm E}\left(\tilde{y}|y\right)$$ です。

欠損している応答と観測された応答の同時分布は、次のように多変量正規分布になります。

$$\left(\begin{array}{l}\tilde{y}\\ y\end{array}\right)~MVN\left\{\left(\begin{array}{l}\tilde{X}\beta \\ X\beta \end{array}\right),\left(\begin{array}{cc}{\Sigma }_{\tilde{y}}& {\Sigma }_{\tilde{y}y}\\ {\Sigma }_{y\tilde{y}}& {\Sigma }_y\end{array}\right)\right\}.$$

多変量正規分布の性質を使用すると、補定された条件付き期待値は次で与えられます。

$${\rm E}\left(\tilde{y}|y\right)=\tilde{X}\beta +{\Sigma }_{\tilde{y}y}{\Sigma }_y^{-1}(y-X\beta ).$$

メモ

mvregress は、欠けている応答値のみを補定します。計画行列では、欠損値をもつ観測が削除されます。

参照

[1] Little, Roderick J. A., and Donald B. Rubin. Statistical Analysis with Missing Data. 2nd ed., Hoboken, NJ: John Wiley & Sons, Inc., 2002.

[2] Meng, Xiao-Li, and Donald B. Rubin. “Maximum Likelihood Estimation via the ECM Algorithm.” Biometrika. Vol. 80, No. 2, 1993, pp. 267–278.

[3] Sexton, Joe, and A. R. Swensen. “ECM Algorithms that Converge at the Rate of EM.” Biometrika. Vol. 87, No. 3, 2000, pp. 651–662.

[4] Dempster, A. P., N. M. Laird, and D. B. Rubin. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society. Series B, Vol. 39, No. 1, 1977, pp. 1–37.

バージョン履歴

R2006b で導入

参考

manova1 | mvregresslike

トピック

• 多変量回帰の問題の設定
• 多変量一般線形モデル
• 同時相関を使用する固定効果のパネル モデル
• 縦方向解析
• 多変量線形回帰
• 多変量回帰モデルの推定