このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

mnrfit

多項分布ロジスティック回帰

構文

B = mnrfit(X,Y)

B = mnrfit(X,Y,Name,Value)

[B,dev,stats] = mnrfit(___)

説明

B = mnrfit(X,Y) は、X の予測子における Y のノミナル応答の多項分布ロジスティック回帰に対する係数推定値の行列 B を返します。

B = mnrfit(X,Y,Name,Value) は、1 つまたは複数の Name,Value ペア引数で指定された追加オプションを使用して、多項分布モデル近似の係数推定値の行列 B を返します。

たとえば、ノミナル モデル、順序モデル、階層モデルを近似するか、リンク関数を変更できます。

[B,dev,stats] = mnrfit(___) は、あてはめの逸脱度 dev と、前の入力引数のいずれかに対する stats 構造体も返します。stats には、自由度、係数推定値の標準誤差、残差などのモデルの統計量が含まれます。

例

ノミナル応答の多項分布回帰

ノミナル結果の多項分布回帰を近似して、結果を解釈します。

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

load fisheriris

列ベクトル species は、3 種類のアヤメ setosa、versicolor、virginica で構成されています。double 行列 meas は、花に関する 4 種類の測定値、がく片の長さと幅 (cm) と花弁の長さと幅 (cm) で構成されています。

categorical 配列を使用してノミナル応答変数を定義します。

sp = categorical(species);

多項分布回帰モデルを近似し、測定値を使用して種類を予測します。

[B,dev,stats] = mnrfit(meas,sp);
B

B = 5×2
103 ×

    1.9078    0.0426
    0.6371    0.0025
   -0.5375    0.0067
   -0.4879   -0.0094
   -2.6110   -0.0183

これは、4 つのすべての予測子において個別の勾配、つまり meas の各カテゴリをもつ、応答カテゴリ相対リスクのノミナル モデルです。B の最初の行には、最初の 2 つの応答カテゴリ (setosa および versicolor と基準カテゴリ virginica) の相対リスクの切片項が含まれます。最後の 4 つの行には、最初の 2 つのカテゴリのモデルに対する勾配が含まれます。mnrfit は 3 番目のカテゴリを基準カテゴリとして受け入れます。

アヤメの花が種類 2 (versicolor) である場合と種類 3 (virginica) である場合の相対リスクは、2 つの確率 (種類 2 である確率と種類 3 である確率) の比率です。相対リスクのモデルは次のようになります。

$$\ln \left(\frac{{\pi }_{versicolor}}{{\pi }_{virginica}}\right)=42.6+2.5X_1+6.7X_2-9.4X_3-18.3X_4.$$

係数は、相対リスクに対する予測子変数の影響と、あるカテゴリに含まれる確率と基準カテゴリに含まれる確率の対数オッズの両方を表します。たとえば、他の条件がすべて等しい場合、推定された係数 2.5 は、1 番目の測定値 $$X_1$$ が 1 単位増加すると、種類 3 (virginica) になる場合に対する種類 2 (versicolor) になる場合の相対リスクが exp(2.5) 倍になることを示します。他の条件がすべて等しい場合、$$X_1$$ が 1 単位増加すると、virginica ではなく versicolor になる相対的な対数オッズは 2.5 倍になります。

係数が無限大または負の無限大に収束する場合、推定される係数はオペレーティング システムによってわずかに異なる可能性があります。

モデル係数の統計的な有意性を確認します。

stats.p

ans = 5×2

         0    0.0000
         0    0.0281
         0    0.0000
         0    0.0000
         0    0.0000

小さい $$p$$ 値は、virginica ではなく setosa になる (種類 3 と比較した種類 1 の) 相対的なリスクおよび virginica ではなく versicolor になる (種類 3 と比較した種類 2 の) 相対的なリスクに関してすべての尺度が有意であることを示します。

係数推定値の標準誤差を要求します。

stats.se

ans = 5×2

   12.4038    5.2719
    3.5783    1.1228
    3.1760    1.4789
    3.5403    1.2934
    7.1203    2.0967

係数に対する 95% の信頼限界を計算します。

LL = stats.beta - 1.96.*stats.se;
UL = stats.beta + 1.96.*stats.se;

setosa である場合と virginica である場合の相対リスクのモデルの係数に対する信頼区間 (B の係数の最初の列) を表示します。

[LL(:,1) UL(:,1)]

ans = 5×2
103 ×

    1.8835    1.9321
    0.6301    0.6441
   -0.5437   -0.5313
   -0.4949   -0.4810
   -2.6249   -2.5970

versicolor である場合と virginica である場合の相対リスクのモデルの係数に対する信頼区間 (B の係数の 2 番目の列) を見つけます。

[LL(:,2) UL(:,2)]

ans = 5×2

   32.3049   52.9707
    0.2645    4.6660
    3.7823    9.5795
  -11.9644   -6.8944
  -22.3957  -14.1766

順序応答の多項分布回帰

カテゴリ間に自然な順序があるカテゴリカル応答について、多項分布回帰モデルを近似します。

標本データを読み込み、予測子変数を定義します。

load carbig
X = [Acceleration Displacement Horsepower Weight];

予測子変数は、自動車の速度、エンジン排気量、馬力および重量です。応答変数は、ガロンあたりの走行マイル数 (mpg) です。

応答値 9 ～ 19 を 1、20 ～ 29 を 2、30 ～ 39 を 3 および 40 ～ 48 のように 4 の範囲でラベル付けすることで、MPG を 9 から 48 までの 4 つのレベルに分類する順序応答変数を作成します。

miles = ordinal(MPG,{'1','2','3','4'},[],[9,19,29,39,48]);

応答変数 miles の順序応答モデルを近似します。

[B,dev,stats] = mnrfit(X,miles,'model','ordinal');
B

B = 7×1

  -16.6895
  -11.7208
   -8.0606
    0.1048
    0.0103
    0.0645
    0.0017

B の最初の 3 つの要素はモデルの切片項で、B の最後の 4 つの要素は共変量の係数です。これらはすべてのカテゴリに共通であると仮定します。このモデルは、"並列回帰" ("比例オッズ" モデルとも呼ばれます) に対応します。モデルのカテゴリ間の切片は異なりますが、勾配は共通しています。これは、順序モデルの既定値である 'interactions','off' の名前と値のペア引数を使用して指定できます。

[B(1:3)'; repmat(B(4:end),1,3)]

ans = 5×3

  -16.6895  -11.7208   -8.0606
    0.1048    0.1048    0.1048
    0.0103    0.0103    0.0103
    0.0645    0.0645    0.0645
    0.0017    0.0017    0.0017

モデルのリンク関数は、順序モデルの既定値である logit ('link','logit') です。係数は、自動車の mpg が特定の値以下である場合とその値よりも大きい場合の相対リスクまたは対数オッズを表します。

この例の比例オッズ モデルは以下のようになります。

$$\begin{array}{l}\ln \left(\frac{P\left(mpg\le 19\right)}{P\left(mpg>19\right)}\right)=-16.6895+0.1048X_A+0.0103X_D+0.0645X_H+0.0017X_W\\ \\ \ln \left(\frac{P\left(mpg\le 29\right)}{P\left(mpg>29\right)}\right)=-11.7208+0.1048X_A+0.0103X_D+0.0645X_H+0.0017X_W\\ \\ \ln \left(\frac{P\left(mpg\le 39\right)}{P\left(mpg>39\right)}\right)=-8.0606+0.1048X_A+0.0103X_D+0.0645X_H+0.0017X_W\end{array}$$

たとえば、係数推定値 0.1048 は、他のすべてが一定であるとして、速度での 1 単位の変更が、自動車の mpg が 19 以下である場合と 19 を超える場合、29 以下である場合と 29 を超える場合、39 以下である場合と 39 を超える場合のオッズに exp(0.01048) の係数で影響することを示します。

係数の有意性を評価します。

stats.p

ans = 7×1

    0.0000
    0.0000
    0.0000
    0.1899
    0.0350
    0.0000
    0.0118

エンジン排気量、馬力および自動車の重量についての $$p$$ 値 0.035、0.0000 および 0.0118 は、自動車の mpg が特定の値より大きくなるのではなくその値以下になるオッズに関して、これらの因子が有意であることを示しています。

階層型多項分布回帰モデル

階層型多項分布回帰モデルを近似します。

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

load('smoking.mat');

データセット smoking には、次の 5 つの変数があります。性別、年齢、重量、収縮期血圧、拡張期血圧です。性別はバイナリ変数で、1 は女性患者、0 は男性患者を示します。

応答変数を定義します。

Y = categorical(smoking.Smoker);

Smoker のデータには次の 4 つのカテゴリが含まれています。

• 0:非喫煙者、1 日に 0 本

• 1:喫煙者、1 日に 1 ～ 5 本

• 2:喫煙者、1 日に 6 ～ 10 本

• 3:喫煙者、1 日に 11 本以上

予測子変数を定義します。

X = [smoking.Sex smoking.Age smoking.Weight...
    smoking.SystolicBP smoking.DiastolicBP];

階層型多項分布モデルを近似します。

[B,dev,stats] = mnrfit(X,Y,'model','hierarchical');
B

B = 6×3

   43.8148    5.9571   44.0712
    1.8709   -0.0230    0.0662
    0.0188    0.0625    0.1335
    0.0046   -0.0072   -0.0130
   -0.2170    0.0416   -0.0324
   -0.2273   -0.1449   -0.4824

B の最初の列には、非喫煙者と喫煙者である場合を比較した相対リスクのモデルについての切片と係数推定値が含まれます。2 番目の列には、対象者を喫煙者として、1 日 1 ～ 5 本吸う場合と、1 日 6 本以上吸う場合の対数オッズをモデル化するパラメーター推定が含まれています。最後の 3 番目の列には、対象者が 1 日 6 本以上吸う喫煙者とすると、1 日 6 ～ 10 本以上吸う場合と、1 日 11 本以上吸う場合の対数オッズをモデル化するパラメーター推定が含まれています。

係数はカテゴリ間で異なります。これは、階層モデルの既定値である 'interactions','on' の名前と値のペア引数を使用して指定できます。したがって、この例のモデルは以下のようになります。

$$\ln \left(\frac{P\left(y=0\right)}{P\left(y>0\right)}\right)=43.8148+1.8709X_S+0.0188X_A+0.0046X_W-0.2170X_{SBP}-0.2273X_{DBP}$$

$$\ln \left(\frac{P\left(1\le y\le 5\right)}{P\left(y>5\right)}\right)=5.9571-0.0230X_S+0.0625X_A-0.0072X_W+0.0416X_{SBP}-0.1449X_{DBP}$$

$$\ln \left(\frac{P\left(6\le y\le 10\right)}{P\left(y>10\right)}\right)=44.0712+0.0662X_S+0.1335X_A-0.0130X_W-0.0324X_{SBP}-0.4824X_{DBP}$$

たとえば、1.8709 の係数推定値は、他のすべてが一定であることを条件に、喫煙者である場合と非喫煙者である場合の尤度が、性別が女性から男性へ変わると exp(1.8709) = 6.49 倍ずつ増えることを示します。

項の統計的な有意性を評価します。

stats.p

ans = 6×3

    0.0000    0.5363    0.2149
    0.3549    0.9912    0.9835
    0.6850    0.2676    0.2313
    0.9032    0.8523    0.8514
    0.0009    0.5187    0.8165
    0.0004    0.0483    0.0545

性別、年齢、体重は、どのレベルにおいても有意性がありません。0.0009 および 0.0004 という $$p$$ 値は、喫煙者である場合と非喫煙者である場合の相対リスクにおいて両方のタイプの血圧が有意であることを示しています。0.0483 という $$p$$ 値は、1 日に 0 ～ 5 本吸う人と 1 日に 5 本より多く吸う人のオッズに対して拡張期血圧のみが有意であることを示しています。同様に、0.0545 という $$p$$ 値は、1 日に 6 ～ 10 本吸う人と 1 日に 10 本より多く吸う人のオッズに対して拡張期血圧が有意であることを示しています。

有意でない因子が相関するかどうかを確認します。性別でグループ分けされた年齢と体重の散布図を描画します。

figure()
gscatter(smoking.Age,smoking.Weight,smoking.Sex)
legend('Male','Female')
xlabel('Age')
ylabel('Weight')

個人の体重の範囲は、性別によって異なるように見えます。年齢は、性別や体重との明らかな相関はないようです。年齢は有意ではありませんが、体重は性別と相関しているようなので、両方を排除して、モデルを再構築できます。

モデルから年齢と体重を排除し、予測子変数として性別、収縮期血圧、拡張期血圧を使用する階層モデルを近似します。

X = double([smoking.Sex smoking.SystolicBP...
smoking.DiastolicBP]);
[B,dev,stats] = mnrfit(X,Y,'model','hierarchical');
B

B = 4×3

   44.8456    5.3230   25.0248
    1.6045    0.2330    0.4982
   -0.2161    0.0497    0.0179
   -0.2222   -0.1358   -0.3092

ここで、1.6045 の係数推定値は、非喫煙者である場合と喫煙者である場合の尤度が、性別が男性から女性へ変わると exp(1.6045) = 4.97 倍増えることを示します。収縮期血圧における単位の増加は、非喫煙者である場合と喫煙者ある場合の尤度が exp(–.2161) = 0.8056 減少することを示します。同様に、拡張期血圧における単位の増加は、非喫煙者である場合と喫煙者ある場合の相対リスクが exp(–.2222) = 0.8007 減少することを示します。

項の統計的な有意性を評価します。

stats.p

ans = 4×3

    0.0000    0.4715    0.2325
    0.0210    0.7488    0.6362
    0.0010    0.4107    0.8899
    0.0003    0.0483    0.0718

0.0210、0.0010、0.0003 という $$p$$ 値は、モデルの他の項が与えられた場合に、非喫煙者と喫煙者の相対リスクに対して性別および両方のタイプの血圧の項が有意であることを示しています。0.0483 という $$p$$ 値に基づくと、対象者が喫煙者であることを条件にして、1 日に 1 ～ 5 本吸う場合と 1 日に 5 本より多く吸う場合の相対リスクに対して拡張期血圧は有意であるように見えます。3 番目の列の $$p$$ 値はどれも 0.05 未満ではないため、この対象者が 1 日に 5 本より多く吸うことを条件にして、6 ～ 10 本吸う場合と 10 本より多く吸う場合の相対リスクに対して、どの変数も統計的に有意ではないことがわかります。

入力引数

X — 予測子変数に関する観測

n 行 p 列の行列

予測子変数に関する観測値。n 行 p 列の行列として指定します。X には p 個の予測子に対する n 個の観測値が含まれます。

メモ

mnrfit は自動的に定数項 (切片) をすべてのモデルに含めます。X に 1 の列を入力しないでください。

データ型: single | double

Y — 応答値

n 行 k 列の行列 | n 行 1 列の列ベクトル

応答値。列ベクトルまたは行列として指定します。Y は次のいずれかになります。

• n 行 k 列の行列。Y(i,j) は、X(i,:) で与えられる予測子の組み合わせに対する多項カテゴリ j の結果の数です。この場合、各予測子の組み合わせごとに観測値数が作成されます。

• 各観測に対する応答値を示す 1 ～ k の整数スカラー の n 行 1 列の列ベクトル。この例ではすべての標本サイズが 1 です。

• 各観測に対する応答のノミナル値または順序値を表す、n 行 1 列の categorical 配列。この例ではすべての標本サイズが 1 です。

名前と値のペアの引数

オプションの Name,Value 引数のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

例: 'Model','ordinal','Link','probit' は、プロビット リンク関数をもつ順序モデルを指定します。

Model — 近似するモデルのタイプ

'nominal' (既定値) | 'ordinal' | 'hierarchical'

近似するモデルのタイプ。'Model' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

'nominal'	既定の設定。応答カテゴリに順序はありません。
'ordinal'	応答カテゴリは自然な順序になっています。
'hierarchical'	応答カテゴリの選択は、逐次的/入れ子です。

例: 'Model','ordinal'

Interactions — 多項分布カテゴリと係数間の交互作用を表すインジケーター

'on' | 'off'

多項分布カテゴリと係数間の交互作用を表すインジケーター。'Interactions' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

'on'	ノミナル モデルと階層モデルの既定値。カテゴリ間で異なる係数をもつモデルを近似します。
'off'	順序モデルの既定値。すべての多項分布カテゴリ間で予測子変数の係数の共通セットをもつモデルを近似します。多くの場合、これは "並列回帰" または "比例オッズ モデル" と呼ばれます。

いずれの場合も、モデルはカテゴリ間で異なる切片をもちます。'Interactions' の選択により、出力配列 B の次元が決まります。

例: 'Interactions','off'

Link — リンク関数

'logit' (既定値) | 'probit' | 'comploglog' | 'loglog'

順序モデルと階層モデルに使用するリンク関数。'Link' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

'logit'	既定。f(γ) = ln(γ/(1 –γ))
'probit'	f(γ) = Φ-1(γ) — 誤差項は分散 1 で正規分布されます。
'comploglog'	相補 log-log f(γ) = ln(–ln(1 – γ))
'loglog'	f(γ) = ln(–ln(γ))

リンク関数は、応答確率と、予測子の線形結合である Xβ との関係を定義します。リンク関数は、モデルが順序応答または逐次的/入れ子応答のどちらに対応するかに応じて、累積確率または条件付き確率の関数になります。たとえば、順序モデルの場合、γ はカテゴリ 1 ～ j になる累積確率と、ロジット リンク関数をもつモデルを次のように表します。

$$\ln \left(\frac{\gamma }{1-\gamma }\right)=\ln \left(\frac{{\pi }_1+{\pi }_2+\cdots +{\pi }_j}{{\pi }_{j+1}+\cdots +{\pi }_k}\right)={\beta }_{0j}+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p,$$

ここで、k は最後のカテゴリを表します。

ノミナル モデルの場合は 'Link' パラメーターを指定できません。このモデルでは常に多項ロジット リンクを使用します。

$$\ln \left(\frac{{\pi }_j}{{\pi }_r}\right)={\beta }_{j0}+{\beta }_{j1}{X}_{j1}+{\beta }_{j2}{X}_{j2}+\cdots +{\beta }_{jp}{X}_{jp},j=1,\dots ,k-1,$$

ここで、π はカテゴリカル確率を表し、r は基準カテゴリに対応します。mnrfit では、ノミナル モデルの基準カテゴリとして最後のカテゴリを使用します。

例: 'Link','loglog'

EstDisp — 分散パラメーターを推定するインジケーター

'off' (既定値) | 'on'

分散パラメーターを推定するインジケーター。'EstDisp' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

'off'	既定の設定。理論分散値 1 を使用します。
'on'	標準誤差の計算における多項分布の分散パラメーターを推定します。

例: 'EstDisp','on'

出力引数

B — 係数推定値

ベクトル | 行列

Y の応答の多項分布ロジスティック回帰の係数推定値。ベクトルまたは行列として返されます。

• 'Interaction' が 'off' の場合、B は k – 1 + p ベクトルになります。B の最初の k – 1 行は、k – 1 多項分布カテゴリごとに 1 つずつ、切片項に対応し、残りの p 行は最初の k – 1 カテゴリすべてに共通の予測子係数に対応します。

• 'Interaction' が 'on' の場合、B は (p + 1) 行 (k – 1) 列の行列になります。B の各列は、最初の k – 1 多項分布カテゴリに 1 つずつ、推定された切片項と予測子係数に対応します。

mnrfit では最後のカテゴリを基準カテゴリとして取るため、k 番目のカテゴリの推定はゼロとして取られます。

dev — あてはめの逸脱度

スカラー値

あてはめの逸脱度。スカラー値として返されます。これは達成可能な最大対数尤度と近似されたモデルで達成された対数尤度との間の差の 2 倍になります。これは、逸脱度残差の合計に対応します。

$$dev=2*{\displaystyle \sum _i^n{\displaystyle \sum _j^k{y}_{ij}*\log \left(\frac{y_{ij}}{{\widehat{\pi }}_{ij}*m_i}\right)}}={\displaystyle \sum _i^{n}r{d}_i},$$

ここで、rd_i は逸脱度残差です。逸脱度残差の詳細は、「stats」を参照してください。

stats — モデルの統計量

構造体

モデルの統計量。以下のフィールドを含む構造体として返されます。

beta	係数推定値。これらは B と同じです。
dfe	誤差に対する自由度 'Interactions' が 'off' の場合、自由度は n*(k – 1) – (k – 1 + p) です。 'Interactions' が 'on' の場合、自由度は (n – p + 1)*(k – 1) です。
sfit	推定された分散パラメーター。
s	理論的または推定された分散パラメーター。 'Estdisp' が 'off' の場合、s は理論的分散パラメーターを表す 1 です。 'Estdisp' が 'on' の場合、s は推定された分散パラメーター sfit と等しくなります。
estdisp	理論的または推定された分散パラメーターのインジケーター。
se	係数推定値 B の標準誤差。
coeffcorr	B の推定相関行列。
covb	B の推定共分散行列。
t	B の t 統計量。
p	B の p 値。
resid	生の残差。観測値から当てはめ値を減算した値。 $$r_{ij}=y_{ij}-{\widehat{\pi }}_{ij}*m_i,\{\begin{array}{c}i=1,\cdots ,n\\ j=1,\cdots ,k\end{array},$$ ここで、πij はカテゴリカル確率、累積確率または条件付き確率であり、mi は対応する標本のサイズです。
residp	ピアソン残差。これは、推定標準偏差でスケーリングした生の残差です。 $$r{p}_{ij}=\frac{r_{ij}}{{\widehat{\sigma }}_{ij}}=\frac{y_{ij}-{\widehat{\pi }}_{ij}*m_i}{\sqrt{{\widehat{\pi }}_{ij}*\left(1-{\widehat{\pi }}_{ij}\right)*m_i}},\{\begin{array}{c}i=1,\cdots ,n\\ j=1,\cdots ,k\end{array},$$ ここで、πij はカテゴリカル確率、累積確率または条件付き確率であり、mi は対応する標本のサイズです。
residd	逸脱度残差。 $$r{d}_i=2*{\displaystyle {\sum }_j^k{y}_{ij}*\log \left(\frac{y_{ij}}{{\widehat{\pi }}_{ij}*m_i}\right)},i=1,\cdots ,n.$$ ここで、πij はカテゴリカル確率、累積確率または条件付き確率であり、mi は対応する標本のサイズです。

アルゴリズム

mnrfit は、X または Y 内の NaN を欠損値として扱い、無視します。

参照

[1] McCullagh, P., and J. A. Nelder. Generalized Linear Models. New York: Chapman & Hall, 1990.

[2] Long, J. S. Regression Models for Categorical and Limited Dependent Variables. Sage Publications, 1997.

[3] Dobson, A. J., and A. G. Barnett. An Introduction to Generalized Linear Models. Chapman and Hall/CRC. Taylor & Francis Group, 2008.

参考

fitglm | mnrval | glmfit | glmval

トピック

• 多項分布
• ノミナル応答の多項モデル
• 順序応答の多項モデル
• 階層型多項モデル

R2006b で導入

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい