このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

glmfit

一般化線形回帰モデルの当てはめ

構文

b = glmfit(X,y,distr)

b = glmfit(X,y,distr,Name,Value)

[b,dev] = glmfit(___)

[b,dev,stats] = glmfit(___)

説明

b = glmfit(X,y,distr) は、分布 distr を使用して、X 内の予測子に対する y 内の応答の一般化線形回帰モデルにおける係数推定値が含まれているベクトル b を返します。

b = glmfit(X,y,distr,Name,Value) では、1 つ以上の名前と値の引数を使用して追加オプションを指定します。たとえば、'Constant','off' を指定するとモデルから定数項を省略できます。

[b,dev] = glmfit(___) は、当てはめの逸脱度の値 dev も返します。

[b,dev,stats] = glmfit(___) は、モデルの統計量 stats も返します。

例

プロビット リンクによる一般化線形モデルの当てはめ

一般化線形回帰モデルを当てはめ、当てはめたモデルを使用して予測子データの予測値 (推定値) を計算します。

標本データ セットを作成します。

x = [2100 2300 2500 2700 2900 3100 ...
     3300 3500 3700 3900 4100 4300]';
n = [48 42 31 34 31 21 23 23 21 16 17 21]';
y = [1 2 0 3 8 8 14 17 19 15 17 21]';

x には予測子変数の値が格納されています。y の各値は対応する n の試行回数に対する成功回数です。

x における y についてのプロビット回帰モデルを当てはめます。

b = glmfit(x,[y n],'binomial','Link','probit');

推定成功回数を計算します。

yfit = glmval(b,x,'probit','Size',n);

観測された成功率と推定された成功率を x の値に対してプロットします。

plot(x,y./n,'o',x,yfit./n,'-')

カスタム リンク関数を使用した一般化線形モデルの当てはめ

カスタム リンク関数を定義し、それを使用して一般化線形回帰モデルを当てはめます。

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

load fisheriris

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

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

X = meas(51:end,:);
y = strcmp('versicolor',species(51:end));

ロジット リンク関数のカスタム リンク関数を定義します。リンク関数、リンク関数の導関数、逆リンク関数を定義する 3 つの関数ハンドルを定義します。これらを cell 配列に格納します。

link = @(mu) log(mu./(1-mu));
derlink = @(mu) 1./(mu.*(1-mu));
invlink = @(resp) 1./(1+exp(-resp));
F = {link,derlink,invlink};

glmfit とカスタム リンク関数を使用してロジスティック回帰モデルを当てはめます。

b = glmfit(X,y,'binomial','link',F)

b = 5×1

   42.6378
    2.4652
    6.6809
   -9.4294
  -18.2861

組み込みのリンク関数 logit を使用して一般化線形モデルを当てはめ、結果を比較します。

b = glmfit(X,y,'binomial','link','logit')

b = 5×1

   42.6378
    2.4652
    6.6809
   -9.4294
  -18.2861

逸脱度検定の実行

切片と各予測子の線形項を含む一般化線形回帰モデルを当てはめます。このモデルが定数モデルよりも大幅に良好な近似をしているかどうかを判別する逸脱度検定を実行します。

基となる 2 つの予測子 X(:,1) および X(:,2) のポアソン乱数を使って標本データを生成します。

rng('default') % For reproducibility
rndvars = randn(100,2);
X = [2 + rndvars(:,1),rndvars(:,2)];
mu = exp(1 + X*[1;2]);
y = poissrnd(mu);

切片と各予測子の線形項を含む一般化線形回帰モデルを当てはめます。

[b,dev] = glmfit(X,y,'poisson');

2 番目の出力引数 dev は、当てはめの逸脱度です。

切片のみを含む一般化線形回帰モデルを当てはめます。glmfit でモデルに定数項を含めないように、予測子変数を 1 の列として指定し、'Constant' を 'off' として指定します。

[~,dev_noconstant] = glmfit(ones(100,1),y,'poisson','Constant','off');

dev_constant と dev の差を計算します。

D = dev_noconstant - dev

D = 2.9533e+05

D は自由度が 2 のカイ二乗分布となります。自由度は、dev に対応するモデル内の推定パラメーターの数と定数モデル内の推定パラメーターの数の差に等しくなります。逸脱度検定の p 値を求めます。

p = 1 - chi2cdf(D,2)

p = 0

小さい p 値は、モデルが定数と大きく異なっていることを示します。

代わりに、関数fitglmを使用してポアソン データの一般化線形回帰モデルを作成することもできます。モデルの表示に統計 (Chi^2-statistic vs. constant model) と p 値が含まれます。

mdl = fitglm(X,y,'y ~ x1 + x2','Distribution','poisson')

mdl = 
Generalized linear regression model:
    log(y) ~ 1 + x1 + x2
    Distribution = Poisson

Estimated Coefficients:
                   Estimate       SE        tStat     pValue
                   ________    _________    ______    ______

    (Intercept)     1.0405      0.022122    47.034      0   
    x1              0.9968      0.003362    296.49      0   
    x2               1.987     0.0063433    313.24      0   


100 observations, 97 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 2.95e+05, p-value = 0

当てはめたモデル オブジェクトで関数devianceTestを使用することもできます。

devianceTest(mdl)

ans=2×4 table
                             Deviance     DFE     chi2Stat     pValue
                            __________    ___    __________    ______

    log(y) ~ 1              2.9544e+05    99                         
    log(y) ~ 1 + x1 + x2         107.4    97     2.9533e+05       0  

入力引数

X — 予測子変数

数値行列

予測子変数。n 行 p 列の数値行列として指定します。ここで、n は観測値の数、p は予測子変数の数です。X の各列が 1 つの変数を表し、各行が 1 つの観測値を表します。

既定では、glmfit はモデルに定数項を含めます。X に 1 の列を直接追加しないでください。glmfit の既定の動作を変更するには、名前と値の引数 'Constant' を指定します。

データ型: single | double

y — 応答変数

ベクトル | 行列

応答変数。ベクトルまたは行列として指定します。

• distr が 'binomial' ではない場合、y は、n 行 1 列のベクトルでなければなりません。ここで、n は観測値の数です。y の各エントリは X の対応する行に対する応答です。データ型は single または double でなければなりません。

• distr が 'binomial' である場合、y は、観測ごとの成功や失敗を示す n 行 1 列のベクトルか、1 列目が観測ごとの成功数を示し、2 列目が観測ごとの試行回数を示す n 行 2 列のベクトルになります。

データ型: single | double | logical | categorical

distr — 応答変数の分布

'normal' (既定値) | 'binomial' | 'poisson' | 'gamma' | 'inverse gaussian'

応答変数の分布。次の表のいずれかの値として指定します。

値	説明
'normal'	正規分布 (既定の設定)
'binomial'	二項分布
'poisson'	ポアソン分布
'gamma'	ガンマ分布
'inverse gaussian'	逆ガウス分布

名前と値の引数

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

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

例: b = glmfit(X,y,'normal','link','probit') は応答の分布が正規であることを指定し、プロビット リンク関数を使用するように glmfit に指示します。

B0 — 係数推定値の初期値

数値ベクトル

係数推定値の初期値。数値ベクトルとして指定します。既定値は、入力データから派生する初期近似値です。

データ型: single | double

Constant — 定数項のインジケーター

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

当てはめにおける定数項 (切片) のインジケーター。モデルに定数項を含める場合は 'on'、モデルから定数項を削除する場合は 'off' を指定します。

• 'on' (既定の設定) — glmfit はモデルに定数項を含め、(p + 1) 行 1 列の係数推定値 b のベクトルを返します。ここで、p は X 内の予測子の個数です。定数項の係数は b の最初の要素です。

• 'off' — glmfit は定数項を省略し、p 行 1 列の係数推定値 b のベクトルを返します。

例: 'Constant','off'

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

'off' ('binomial' および 'poisson' 分布の場合) (既定値) | 'on'

'binomial' および 'poisson' 分布の分散パラメーターを計算するインジケーター。'on' または 'off' として指定します。

値	説明
'on'	標準誤差を計算するときに分散パラメーターを推定します。分散パラメーターの推定値は、ピアソン残差の二乗和を誤差の自由度 (DFE) で除算した値です。
'off'	標準誤差を計算するときに理論値 1 を使用します (既定の設定)。

近似関数は常に他の分布の分散を予測します。

例: 'EstDisp','on'

Link — リンク関数

正準リンク関数 (既定値) | スカラー値 | カスタム リンク関数の構造体または cell 配列

正準リンク関数の代わりに使用するリンク関数。次の表のいずれかの組み込みリンク関数、またはカスタム リンク関数として指定します。

リンク関数名	リンク関数	平均 (逆) 関数
'identity' ('normal' 分布の既定の設定)	f(μ) = μ	μ = Xb
'log' ('poisson' 分布の既定の設定)	f(μ) = log(μ)	μ = exp(Xb)
'logit' ('binomial' 分布の既定の設定)	f(μ) = log(μ/(1 – μ))	μ = exp(Xb) / (1 + exp(Xb))
'probit'	f(μ) = Φ–1(μ)。Φ は標準正規分布の累積分布関数です。	μ = Φ(Xb)
'loglog'	f(μ) = log(–log(μ))	μ = exp(–exp(Xb))
'comploglog'	f(μ) = log(–log(1 – μ))	μ = 1 – exp(–exp(Xb))
'reciprocal' ('gamma' 分布の既定の設定)	f(μ) = 1/μ	μ = 1/(Xb)
p (数値、'inverse gaussian' 分布の既定の設定、p = –2)	f(μ) = μp	μ = Xb1/p

'Link' の既定値は、引数 distr で指定された応答変数の分布に基づく正準リンク関数です。

カスタム リンク関数は構造体または cell 配列を使用して指定できます。

• 3 つのフィールドをもつ構造体。構造体 (S など) の各フィールドは、入力のベクトルを受け入れ、同じサイズのベクトルを返す関数ハンドルを保持します。

  • S.Link — リンク関数、f(μ) = S.Link(μ)

  • S.Derivative — リンク関数の導関数

  • S.Inverse — 逆リンク関数、μ = S.Inverse(Xb)

• リンク関数 (FL(mu))、リンク関数の導関数 (FD = dFL(mu)/dmu)、逆リンク関数 (FI = FL^(–1)) を定義する {FL FD FI} の形式の cell 配列。各エントリは、入力のベクトルを受け入れ、同じサイズのベクトルを返す関数ハンドルを保持します。

リンク関数は、平均応答 μ と、予測子の線形結合 X*b の関係である f(μ) = X*b を定義します。

例: 'Link','probit'

データ型: single | double | char | string | struct | cell

Offset — オフセット変数

[ ] (既定値) | 数値ベクトル

近似のオフセット変数。応答 y と同じ長さの数値ベクトルとして指定します。

glmfit は、係数値を 1 で固定した追加の予測子として Offset を使用します。つまり、当てはめの式は次のようになります。

f(μ) = Offset + X*b,

ここで、f はリンク関数、μ は平均応答、X*b は予測子 X の線形結合です。予測子 Offset の係数は 1 です。

たとえば、ポアソン回帰モデルを検討してください。理論上の理由から、カウントの数が予測子 A に比例していると仮定します。log リンク関数を使用し、オフセットに log(A) を指定することにより、この理論上の制約を満たすことをモデルに強制できます。

データ型: single | double

Options — 最適化オプション

statset('glmfit') (既定値) | 構造体

最適化オプション。構造体を指定します。この引数は、glmfit が使用する反復アルゴリズムの制御パラメーターを決定します。

'Options' の値を作成するには、関数 statset を使用するか、次の表に記載されているフィールドと値が含まれている構造体配列を作成します。

フィールド名	値	既定値
Display	アルゴリズムで表示される情報量 'off' — 情報表示なし 'final' — 最終出力を表示	'off'
MaxIter	許容される最大反復回数。正の整数として指定	100
TolX	パラメーターの終了許容誤差。正のスカラーとして指定	1e-6

コマンド ウィンドウで「statset('glmfit')」と入力して、glmfit が名前と値の引数 'Options' で受け入れるフィールドの名前と既定値を表示することもできます。

例: 'Options',statset('Display','final','MaxIter',1000) は、反復アルゴリズムの結果の最終情報を表示し、許容される最大反復回数を 1000 に変更するよう指定します。

データ型: struct

Weights — 観測値の重み

ones(n,1) (既定値) | 非負のスカラー値の n 行 1 列のベクトル

観測値の重み。非負のスカラー値の n 行 1 列のベクトル (n は観測値の数) として指定します。

データ型: single | double

出力引数

b — 係数推定値

数値ベクトル

係数推定値。数値ベクトルとして返されます。

• 'Constant' が 'on' (既定の設定) の場合、glmfit はモデルに定数項を含め、(p + 1) 行 1 列の係数推定値 b のベクトルを返します。ここで、p は X 内の予測子の個数です。定数項の係数は b の最初の要素です。

• 'Constant' が 'off' の場合、glmfit は定数項を省略し、p 行 1 列の係数推定値 b のベクトルを返します。

dev — 当てはめの逸脱度

数値

当てはめの逸脱度。数値として返されます。逸脱度は、一方のモデルが他方のモデルの特別なケースである 2 つのモデルを比較するために役立ちます。2 つのモデルの逸脱度の差異は、2 つのモデル間に推定されるパラメーター数の差異と等しい自由度をもつカイ二乗分布になります。

詳細については、逸脱度を参照してください。

stats — モデルの統計量

構造体

モデルの統計量。次のフィールドをもつ構造体として返されます。

• beta — 係数推定値 b

• dfe — 誤差に対する自由度

• sfit — 推定された分散パラメーター

• s — 理論的または推定された分散パラメーター

• estdisp — 'EstDisp' が 'off' の場合は 0、'EstDisp' が 'on' の場合は 1

• covb — b の推定共分散行列

• se — 係数推定値b の標準誤差のベクトル

• coeffcorr — b の相関行列

• t — b の t 統計量

• p — b の p 値

• resid — 残差のベクトル

• residp — ピアソン残差のベクトル

• residd — 逸脱度残差のベクトル

• resida — Anscombe 残差のベクトル

二項分布またはポアソン分布の分散パラメーターを推定する場合、stats.s は stats.sfit と等しくなります。また、stats.se の要素とそれらの理論値には、要因 stats.s の分だけ差異が生じます。

詳細

逸脱度

逸脱度は、残差二乗和を汎化したものです。飽和モデルと比較した適合度を測定します。

モデル M₁ の逸脱度は、モデル M₁ の対数尤度と飽和モデル M_s の対数尤度の差の 2 倍です。飽和モデルは、推定可能な最大数のパラメーターを含むモデルです。

たとえば、n 個の観測値 (y_i, i = 1, 2, ..., n) があり、X_i^Tβ の値が異なる可能性がある場合、n 個のパラメーターを含む飽和モデルを定義できます。L(b,y) がパラメーター b をもつモデルの尤度関数の最大値を示しているものとします。この場合、モデル M₁ の逸脱度は以下のようになります。

$$-2\left(\log L\left(b_1,y\right)-\log L\left(b_S,y\right)\right),$$

ここで、b₁ および b_s には、それぞれモデル M₁ および飽和モデルの推定パラメーターが含まれます。逸脱度は自由度 n – p のカイ二乗分布になります。ここで、n は飽和モデルのパラメーター数、p はモデル M₁ のパラメーター数です。

2 つの異なる一般化線形回帰モデル M₁ および M₂ があり、M₁ は M₂ 内の項のサブセットをもつと仮定します。モデルの当てはめは、それらの逸脱度 D₁ と D₂ を比較することによって評価できます。逸脱度の差異は以下のとおりです。

$$\begin{array}{l}D=D_2-D_1=-2\left(\log L\left(b_2,y\right)-\log L\left(b_S,y\right)\right)+2\left(\log L\left(b_1,y\right)-\log L\left(b_S,y\right)\right)\\ =-2\left(\log L\left(b_2,y\right)-\log L\left(b_1,y\right)\right).\end{array}$$

差異 D は漸近的にカイ二乗分布となりますが、その自由度 v は、M₁ および M₂ で推定されたパラメーター数の差異に等しくなります。この検定の p 値は 1 — chi2cdf(D,v) を使用して取得できます。

通常は、定数項をもち予測子をもたないモデル M₂ を使用して D を調べます。そのため、D は自由度が p – 1 のカイ二乗分布となります。分散を推定する場合、推定された分散で除算した差異は、分子の自由度が p – 1、分母の自由度が n – p の F 分布となります。

ヒント

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

代替機能

glmfit は、関数の出力引数のみが必要である場合、またはループ内でモデルの当てはめを複数回繰り返す場合に便利です。当てはめたモデルをさらに調べる必要がある場合は、fitglm または stepwiseglm を使用して一般化線形回帰モデル オブジェクト GeneralizedLinearModel を作成します。GeneralizedLinearModel オブジェクトは、glmfit より多くの機能を提供します。

• 当てはめたモデルを調べるには、GeneralizedLinearModel のプロパティを使用します。オブジェクト プロパティには、係数推定値、要約統計量、当てはめ手法および入力データに関する情報が含まれています。

• 応答の予測と、一般化線形回帰モデルの修正、評価および可視化を行うには、GeneralizedLinearModel のオブジェクト関数を使用します。

• glmfit の出力に含まれる情報は、GeneralizedLinearModel のプロパティとオブジェクト関数を使用して取得できます。

  glmfit の出力	GeneralizedLinearModel の同等の値
  b	Coefficients プロパティの Estimate 列を参照してください。
  dev	Deviance プロパティを参照してください。
  stats	コマンド ウィンドウのモデル表示を参照してください。統計量はモデルのプロパティ (CoefficientCovariance、Coefficients、Dispersion、DispersionEstimated、および Residuals) から取得できます。 glmfit の stats.s の分散パラメーターは係数の標準誤差のスケール係数であるのに対し、一般化線形モデルの Dispersion プロパティの分散パラメーターは応答の分散のスケール係数です。したがって、stats.s は Dispersion の値の平方根です。

参照

[1] Dobson, A. J. An Introduction to Generalized Linear Models. New York: Chapman & Hall, 1990.

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

[3] Collett, D. Modeling Binary Data. New York: Chapman & Hall, 2002.

拡張機能

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細は、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

glmval | regress | regstats | GeneralizedLinearModel | fitglm | stepwiseglm

トピック

• 一般化線形モデルによるデータの近似
• 一般化線形モデル