ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

一般化線形モデル

一般化線形モデルとは

線形回帰モデルは、応答と 1 つまたは複数の予測項間の線形関係を記述します。しかし、非線形関係が存在します。非線形回帰では、一般の非線形モデルを説明します。"一般化線形モデル" と呼ばれる非線形モデルの特別なクラスは、線形手法を使用します。

線形モデルが次の特性をもつことを思い出してください。

  • 予測子の値のセットごとに、応答は平均が μ の正規分布になります。

  • 係数ベクトル b は、予測子 X の線形結合 Xb を定義します。

  • モデルは μ = Xb です。

一般化線形モデルでは、これらの特性は次のように一般化されます。

  • 予測子の値のセットごとに、応答には分布が含まれます。これは、正規二項ポアソンガンマ逆ガウス分布 のいずれかで、平均 μ を含むパラメーターをもちます。

  • 係数ベクトル b は、予測子 X の線形結合 Xb を定義します。

  • "リンク関数" f は、モデルを f(μ) = Xb として定義します。

データの準備

回帰の近似を開始するには、データを近似関数に望ましい形式にします。すべての回帰手法は、配列 X の入力データと独立したベクトル y の応答データか、テーブルまたはデータセット配列 tbl 内の入力データと tbl の列としての応答データで始まります。入力データの各行が、1 つの観測値を表します。各列が 1 つの予測子 (変数) を表します。

テーブルまたはデータセット配列 tbl では、次のように 'ResponseVar' の名前と値のペアで応答変数を示します。

mdl = fitglm(tbl,'ResponseVar','BloodPressure');

応答変数は既定で最後の列です。

数値の "カテゴリカル" 予測子を使用できます。カテゴリカル予測子は可能性のある固定セットからの値をとります。

  • 数値配列 X では、'Categorical' の名前と値のペアでカテゴリカル予測子を示します。たとえば、6 つの予測子から 23 が categorical であることを示すには、次のようにします。

    mdl = fitglm(X,y,'Categorical',[2,3]);
    % or equivalently
    mdl = fitglm(X,y,'Categorical',logical([0 1 1 0 0 0]));
  • テーブルまたはデータセット配列 tbl では、近似関数はこれらのデータ型が categorical であることを想定しています。

    • logical ベクトル

    • categorical ベクトル

    • 文字配列

    • string 配列

    数値予測値が categorical であることを示すには、'Categorical' 名前と値のペアを使用します。

欠損数値データは NaN で表されています。他のデータ型用の欠損データを表すには、グループ化変数の欠損値を参照してください。

  • データ行列 X をもつ 'binomial' モデルの応答 y は、次のいずれかです。

    • バイナリ列ベクトル — 各エントリは成功 (1) または失敗 (0) を示します。

    • 整数の 2 列行列 — 1 列目は各観測での成功回数、2 列目はその観測での試行回数を示します。

  • tbl テーブルまたはデータセットをもつ 'binomial' モデルの場合、次の手順を実行します。

    • ResponseVar 名前と値のペアを使用して、各観測での成功回数を示す tbl の列を指定します。

    • BinomialSize 名前と値のペアを使用して、各観測での試行回数を示す tbl の列を指定します。

入力と応答データのデータセット配列

たとえば、Excel® スプレッドシートからデータセット配列を作成するには、次のようにします。

ds = dataset('XLSFile','hospital.xls',...
    'ReadObsNames',true);

ワークスペース変数からデータセット配列を作成するには、次のようにします。

load carsmall
ds = dataset(MPG,Weight);
ds.Year = ordinal(Model_Year);

入力および応答データのテーブル

ワークスペース変数からテーブルを作成するには、次のようにします。

load carsmall
tbl = table(MPG,Weight);
tbl.Year = ordinal(Model_Year);

入力データの数値配列、応答の数値ベクトル

たとえば、ワークスペース変数から数値配列を作成するには、次のようにします。

load carsmall
X = [Weight Horsepower Cylinders Model_Year];
y = MPG;

Excel スプレッドシートから数値配列を作成するには、次のようにします。

[X, Xnames] = xlsread('hospital.xls');
y = X(:,4); % response y is systolic pressure
X(:,4) = []; % remove y from the X matrix

sex などの数値以外のエントリは X には表示されません。

一般化線形モデルとリンク関数の選択

データから一般化線形モデルの分布タイプがわかる場合がよくあります。

応答データのタイプ推奨されるモデルの分布タイプ
任意の実数'normal'
任意の正の数'gamma' または 'inverse gaussian'
任意の非負の整数'poisson'
ゼロから n までの整数 (n は正の固定値)'binomial'

Distribution 名前と値のペアを使用してモデル分布タイプを設定します。モデル タイプを選択した後、平均 µ と線形予測子 Xb 間でマップするためのリンク関数を選択します。

説明
'comploglog'

log(–log((1 – µ))) = Xb

'identity'、分布の既定の設定 'normal'

µ = Xb

'log'、分布の既定の設定 'poisson'

log(µ) = Xb

'logit'、分布の既定の設定 'binomial'

log(µ/(1 – µ)) = Xb

'loglog'

log(–log(µ)) = Xb

'probit'

Φ–1(µ) = Xb、Φ は正規 (ガウス) 累積分布関数

'reciprocal'、分布の既定の設定 'gamma'

µ–1 = Xb

p (数値)、分布 'inverse gaussian' の既定の設定 (p = -2 の場合)

µp = Xb

リンク (FL)、リンクの導関数 (FD)、逆リンク (FI) を定義する 3 つの関数ハンドル (@ を使用して作成) が含まれている、{FL FD FI} という形式の cell 配列。または、FL を含むフィールド LinkFD を含むフィールド Derivative、および FI を含むフィールド Inverse をもつ、関数ハンドルの構造体。

ユーザー指定のリンク関数 (カスタム リンク関数を参照)

既定ではないリンク関数は、主に二項モデルで使用すると便利です。これらの既定ではないリンク関数は 'comploglog''loglog' および 'probit' です。

カスタム リンク関数

リンク関数は、平均応答 µ と予測子の線形結合である Xb = X*b の関係 f( µ) = Xb を定義します。組み込みリンク関数の 1 つを選択するか、リンク関数 FL、リンク関数の導関数 FD およびリンク関数の逆関数 FI を指定してユーザーが独自に定義することができます。

  • リンク関数 FL は f(µ) を計算します。

  • リンク関数 FD の導関数は df(µ)/dµ を計算します。

  • リンク関数 FI の逆関数は g(Xb) = µ を計算します。

カスタム リンク関数は次の 2 つの方法のいずれかで指定できます。どちらの方法にも、µ または Xb を表す値の単一の配列を受け入れて同じサイズの配列を返す、関数ハンドルが含まれています。この関数ハンドルは、以下の cell 配列または構造体のいずれかです。

  • リンク (FL)、リンクの微分 (FD)、逆リンク (FI) を定義する 3 つの関数ハンドル (@ を使用して作成) が含まれている、{FL FD FI} という形式の cell 配列。

  • 3 つのフィールドをもつ構造体 s で、各フィールドには @ を使用して作成された関数ハンドルが含まれます。

    • s.Link — リンク関数

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

    • s.Inverse — リンク関数の逆関数

たとえば、'probit' リンク関数を使用してモデルを近似するには、以下の手順に従います。

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]';
g = fitglm(x,[y n],...
    'linear','distr','binomial','link','probit')
g = 

Generalized Linear regression model:
    probit(y) ~ 1 + x1
    Distribution = Binomial

Estimated Coefficients:
                   Estimate     SE            tStat     pValue    
    (Intercept)      -7.3628       0.66815    -11.02    3.0701e-28
    x1             0.0023039    0.00021352     10.79    3.8274e-27

12 observations, 10 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 241, p-value = 2.25e-54

リンク関数 'probit' と同一の機能のカスタム リンク関数を使用して、同じ近似を実行できます。

s = {@norminv,@(x)1./normpdf(norminv(x)),@normcdf};
g = fitglm(x,[y n],...
    'linear','distr','binomial','link',s)
g = 

Generalized Linear regression model:
    link(y) ~ 1 + x1
    Distribution = Binomial

Estimated Coefficients:
                   Estimate     SE            tStat     pValue    
    (Intercept)      -7.3628       0.66815    -11.02    3.0701e-28
    x1             0.0023039    0.00021352     10.79    3.8274e-27

12 observations, 10 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 241, p-value = 2.25e-54

この 2 つのモデルは同じです。

どちらも同様に、関数ハンドルをもつ cell 配列ではなく、構造体として s を作成できます。

s.Link = @norminv;
s.Derivative = @(x) 1./normpdf(norminv(x));
s.Inverse = @normcdf;
g = fitglm(x,[y n],...
    'linear','distr','binomial','link',s)
g = 

Generalized Linear regression model:
    link(y) ~ 1 + x1
    Distribution = Binomial

Estimated Coefficients:
                   Estimate     SE            tStat     pValue    
    (Intercept)      -7.3628       0.66815    -11.02    3.0701e-28
    x1             0.0023039    0.00021352     10.79    3.8274e-27

12 observations, 10 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 241, p-value = 2.25e-54

近似手法およびモデルの選択

近似モデルを作成するには、次の 2 つの方法があります。

  • 一般化線形モデルの概念をよく理解しているか、あとでモデルを調整して特定の項を追加または除外する場合は、fitglm を使用します。

  • ステップワイズ回帰を使用してモデルを近似する場合は、stepwiseglm を使用します。stepwiseglm は、定数などの 1 つのモデルから開始し、一度に 1 つずつ項を増減する貪欲法で、それ以上改良できなくなるまで、最適な項を毎回選択します。ステップワイズ近似を使用して、有効な項のみを含む適正なモデルを見つけます。

    この結果は、開始モデルに依存します。通常、定数モデルで開始すると、小さいモデルになります。さらに多くの項で開始すると、複雑なモデルになりますが、平均二乗誤差はより小さくなります。

どちらの場合も、近似関数にモデルを指定します (これは stepwiseglm の開始モデルです)。

次のいずれかの方法でモデルを指定します。

短縮されたモデル名

名前モデル タイプ
'constant'モデルは定数 (切片) 項だけを含みます。
'linear'モデルは各予測子に対して切片と線形項を含みます。
'interactions'モデルは、切片、線形項、異なる予測子のペアのすべての積 (二乗項なし) を含みます。
'purequadratic'モデルは、切片、線形項、二乗項を含みます。
'quadratic'モデルは、切片、線形項、交互作用、二乗項を含みます。
'polyijk'モデルは多項式であり、最初の予測子は次数 i まで、2 番目の予測子は次数 j まで、3 番目以降もすべて同様に 0 から 9 までの数値を使用します。たとえば、'poly2111' には、1 つの定数のほかにすべての線形項と積項があり、また、予測子 1 の二乗の項を含んでいます。

項の行列

項行列T は、モデル内の項を指定する t 行 (p + 1) 列の行列です。t は項の個数、p は予測子変数の個数であり、+ 1 は応答変数に相当します。T(i,j) の値は、項 i の変数 j の指数です。

たとえば、3つの予測子変数 ABC と応答変数 YABCY という順序で入力に含まれていると仮定します。T の各行は 1 つの項を表します。

  • [0 0 0 0] — 定数項 (切片)

  • [0 1 0 0]B (A^0 * B^1 * C^0 と等価)

  • [1 0 1 0]A*C

  • [2 0 0 0]A^2

  • [0 1 2 0]B*(C^2)

各項の最後の 0 は、応答変数を表します。一般に、項行列内のゼロの列ベクトルは、応答変数の位置を表します。行列と列ベクトルに予測子と応答変数がある場合、各行の最後の列に応答変数を示す 0 を含めなければなりません。

モデル仕様の式は、次のような形式の文字ベクトルまたは string スカラーです。

'Y ~ terms'

  • Y は応答名です。

  • terms 次が含まれます

    • 変数名

    • + 次の変数を含みます

    • - 次の変数を除外します

    • : 項の積である交互作用を定義します

    • * 交互作用とすべての次数の低い項を定義します

    • * で繰り返されるとおり、^ は予測子をべき乗にし、^ は低い次数の項も含みます

    • () 項をグループ化します

ヒント

式には既定で定数 (切片) 項が含まれます。モデルから定数項を除外するには、式に -1 を含めます。

次に例を示します。

'Y ~ A + B + C' は、切片がある 3 変数線形モデルです。
'Y ~ A + B + C - 1' は、切片がない 3 変数線形モデルです。
'Y ~ A + B + C + B^2' は、切片と B^2 項がある 3 変数モデルです。
'Y ~ A + B^2 + C' は、B^2B 項が含まれるので、前の例と同じです。
'Y ~ A + B + C + A:B' には A*B 項が含まれています。
'Y ~ A*B + C' は、A*B = A + B + A:B なので、前の例と同じです。
'Y ~ A*B*C - A:B:C' には、3 次交互作用を除く ABC 間の交互作用がすべてあります。
'Y ~ A*(B + C + D)' には、すべての線形項、および A と他の各変数の積があります。

モデルのデータへのあてはめ

fitglm または stepwiseglm を使用して近似モデルを作成します。近似手法およびモデルの選択 で説明した方法で、どちらかを選択します。一般化線形モデルでは、正規分布をもつモデルを除き、Distribution 名前と値のペアを 一般化線形モデルとリンク関数の選択 と同じように指定します。たとえば、以下のようにします。

mdl = fitglm(X,y,'linear','Distribution','poisson')
% or
mdl = fitglm(X,y,'quadratic',...
         'Distribution','binomial')

品質の調査と近似モデルの調整

モデルで近似した後に、結果を確認します。

モデル表示

線形回帰モデルの名前または disp(mdl) を入力すると、いくつかの診断情報が表示されます。この表示は、近似モデルがデータを適切にあらわすかどうかを確認するために基本情報のいくつかを提供します。

たとえば、ポアソン モデルを 5 つの予測子のうち 2 つが応答に影響せず、切片の項がないデータにあてはめるには、次のようにします。

rng('default') % for reproducibility
X = randn(100,5);
mu = exp(X(:,[1 4 5])*[.4;.2;.3]);
y = poissrnd(mu);
mdl = fitglm(X,y,...
    'linear','Distribution','poisson')
mdl = 


Generalized Linear regression model:
    log(y) ~ 1 + x1 + x2 + x3 + x4 + x5
    Distribution = Poisson

Estimated Coefficients:
                   Estimate     SE          tStat       pValue    
    (Intercept)     0.039829     0.10793     0.36901       0.71212
    x1               0.38551    0.076116      5.0647    4.0895e-07
    x2             -0.034905    0.086685    -0.40266        0.6872
    x3              -0.17826    0.093552     -1.9054      0.056722
    x4               0.21929     0.09357      2.3436      0.019097
    x5               0.28918      0.1094      2.6432     0.0082126


100 observations, 94 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 44.9, p-value = 1.55e-08

次の点に注意してください。

  • 表示には Estimate 列に各係数の推定値が含まれます。これらの値は真の値 [0;.4;0;0;.2;.3] にかなり近い値ですが、x3 の係数だけは 0 にさほど近くない可能性があります。

  • 係数推定には標準誤差列があります。

  • 予測子 1、4 および 5 に対してレポートされた pValue (標準誤差と仮定して t 統計から導出された) は小さい値です。これらは、応答データ y を作成するために使われた 3 つの予測子です。

  • (Intercept)x2 および x3 に対する pValue は 0.01 より大きい値です。これらの 3 つの予測子は、応答データ y を作成するためには使われませんでした。x3 に対応する pValue.05 をわずかに超える値であるため、有意性が認識される可能性もあります。

  • 表示にはカイ二乗統計量が含まれます。

診断プロット

診断プロットによって外れ値を特定でき、モデルや近似で他の問題を確認できます。これらのプロットについて説明するために、ロジスティック リンク関数をもつ二項回帰について考えます。

"ロジスティック モデル" は比率データに役立ちます。比率 p と重量 w の関係を以下のように定義します。

log[p/(1 – p)] = b1 + b2w

次の例では、二項モデルをデータにあてはめます。このデータは、重量が異なる大型車の測定が含まれる carbig.mat から導出されています。w の各重量には、total に対応した自動車数と poor に対応した燃費の悪い自動車数が含まれています。

total および w に依存する成功のパーセンテージでトライアル数を指定すると、poor の値が二項分布に続くと仮定できます。この分布は、リンク関数 log(µ/(1–µ)) = Xb によって一般化された線形モデルを使用することにより、ロジスティック モデルのコンテキスト内で説明されます。このリンク関数は 'logit' と呼ばれます。

w = [2100 2300 2500 2700 2900 3100 ...
     3300 3500 3700 3900 4100 4300]';
total = [48 42 31 34 31 21 23 23 21 16 17 21]';
poor = [1 2 0 3 8 8 14 17 19 15 17 21]';
mdl = fitglm(w,[poor total],...
    'linear','Distribution','binomial','link','logit')
mdl = 

Generalized Linear regression model:
    logit(y) ~ 1 + x1
    Distribution = Binomial

Estimated Coefficients:
                   Estimate     SE            tStat      pValue    
    (Intercept)       -13.38         1.394    -9.5986    8.1019e-22
    x1             0.0041812    0.00044258     9.4474    3.4739e-21

12 observations, 10 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 242, p-value = 1.3e-54

このモデルがどれだけ適切にデータを近似しているかを確認します。

plotSlice(mdl)

この近似は信頼限界がかなり広く、適正と思われます。

さらに詳細を調べるには、てこ比のプロットを作成します。

plotDiagnostics(mdl)

これは予測子変数によって順序付けられた点をもつ、典型的な回帰です。近似の各点のてこ比は、比較的極端な予測子の値 (どちらの方向でも) をもつ点で高く、平均的な予測子の値をもつ点では低くなります。複数の予測子があり、点が予測子の値で順序付けられていない場合は、てこ比の高い観測値が予測子の値で測定された外れ値となるため、見つけやすくなります。

残差 — 学習データのモデル品質

モデルまたはデータ内の誤差、外れ値または相関を検出できるいくつかの残差プロットが存在します。最もシンプルな残差プロットは既定のヒストグラム プロットであり、これは残差の範囲と頻度を示します。また、確率プロットは残差の分布が正規分布と一致する分散を比較する方法を示します。

次の例は、近似ポアソン モデルの残差プロットを示しています。このデータ構造は、5 つのうち 2 つが応答に影響しない予測子をもち、切片の項をもちません。

rng('default') % for reproducibility
X = randn(100,5);
mu = exp(X(:,[1 4 5])*[2;1;.5]);
y = poissrnd(mu);
mdl = fitglm(X,y,...
    'linear','Distribution','poisson');

残差の検査

plotResiduals(mdl)

ほとんどの残差クラスターは 0 に近くなっていますが、いくつかは ±18 に近い値であるため、別の残差プロットを調べます。

plotResiduals(mdl,'fitted')

大きな残差は近似値のサイズとあまり関連性がないようです。

おそらく、確率プロットがより有益でしょう。

plotResiduals(mdl,'probability')

以下のことが明らかになりました。残差は正規分布していません。代わりに、潜在的なポアソン分布のように裾がより厚くなっています。

予測子の効果とモデルの変更方法を理解するためのプロット

次の例は、各予測子が回帰モデルに与える効果を理解する方法と、モデルを変更して不要な項を削除する方法を示します。

  1. 人為的なデータでいくつかの予測子からモデルを作成します。データは X の 2 番目と 3 番目の列を使用しません。そのため、モデルがこれらの予測子に強く依存していないことを予想します。

    rng('default') % for reproducibility
    X = randn(100,5);
    mu = exp(X(:,[1 4 5])*[2;1;.5]);
    y = poissrnd(mu);
    mdl = fitglm(X,y,...
        'linear','Distribution','poisson');
  2. 応答のスライス プロットを調査します。これは、各予測子の効果を個別に表示します。

    plotSlice(mdl)

    最初の予測子のスケールは、プロットを圧倒しています。[予測子] メニューを使用して、この予測子を無効にします。

    2 番目と 3 番目の予測子の効果がほとんどないことは明らかです。

    青い縦の破線によって表される個々の予測子値をドラッグできます。また、赤い破線の曲線で表されている、同時信頼限界と非同時信頼限界を選択することもできます。予測子のラインをドラッグして、2 番目と 3 番目の予測子に効果がほとんどないことを確認します。

  3. 不要な予測子は、removeTerms または step のいずれかを使用して削除します。1 つの項を削除してから、予期せずに別の項が重要になるケースもあるため、step を使用するほうがより安全です。ただし、step の処理が続行しない場合に、removeTerms が効果的な場合があります。この場合、どちらの方法でも結果は同じです。

    mdl1 = removeTerms(mdl,'x2 + x3')
    mdl1 = 
    
    
    Generalized Linear regression model:
        log(y) ~ 1 + x1 + x4 + x5
        Distribution = Poisson
    
    Estimated Coefficients:
                       Estimate    SE          tStat     pValue     
        (Intercept)    0.17604     0.062215    2.8295       0.004662
        x1              1.9122     0.024638    77.614              0
        x4             0.98521     0.026393    37.328    5.6696e-305
        x5             0.61321     0.038435    15.955     2.6473e-57
    
    
    100 observations, 96 error degrees of freedom
    Dispersion: 1
    Chi^2-statistic vs. constant model: 4.97e+04, p-value = 0
    mdl1 = step(mdl,'NSteps',5,'Upper','linear')
    1. Removing x3, Deviance = 93.856, Chi2Stat = 0.00075551, PValue = 0.97807
    2. Removing x2, Deviance = 96.333, Chi2Stat = 2.4769, PValue = 0.11553
    
    mdl1 = 
    
    
    Generalized Linear regression model:
        log(y) ~ 1 + x1 + x4 + x5
        Distribution = Poisson
    
    Estimated Coefficients:
                       Estimate    SE          tStat     pValue     
        (Intercept)    0.17604     0.062215    2.8295       0.004662
        x1              1.9122     0.024638    77.614              0
        x4             0.98521     0.026393    37.328    5.6696e-305
        x5             0.61321     0.038435    15.955     2.6473e-57
    
    
    100 observations, 96 error degrees of freedom
    Dispersion: 1
    Chi^2-statistic vs. constant model: 4.97e+04, p-value = 0

新しいデータに対する応答を予測またはシミュレート

新しいデータに対する応答を予測するために、3 つの手法を使用できます。

predict

predict メソッドは平均応答の予測と、必要な場合には信頼限界を提供します。

この例では、predict メソッドを使用して予測で信頼区間を予測して取得する方法を示します。

  1. 人為的なデータでいくつかの予測子からモデルを作成します。データは X の 2 番目と 3 番目の列を使用しません。2 番目と 3 番目の予測子の効果がほとんどないことは明らかです。関連する予測子を自動的に含むようにモデルをステップワイズに構築します。

    rng('default') % for reproducibility
    X = randn(100,5);
    mu = exp(X(:,[1 4 5])*[2;1;.5]);
    y = poissrnd(mu);
    mdl = stepwiseglm(X,y,...
        'constant','upper','linear','Distribution','poisson');
    1. Adding x1, Deviance = 2515.02869, Chi2Stat = 47242.9622, PValue = 0
    2. Adding x4, Deviance = 328.39679, Chi2Stat = 2186.6319, PValue = 0
    3. Adding x5, Deviance = 96.3326, Chi2Stat = 232.0642, PValue = 2.114384e-52
  2. 新しいデータをいくつか生成して、データから予測を評価します。

    Xnew = randn(3,5) + repmat([1 2 3 4 5],[3,1]); % new data
    [ynew,ynewci] = predict(mdl,Xnew)
    ynew =
    
       1.0e+04 *
    
        0.1130
        1.7375
        3.7471
    
    
    ynewci =
    
       1.0e+04 *
    
        0.0821    0.1555
        1.2167    2.4811
        2.8419    4.9407

feval

テーブルまたはデータセット配列からモデルを構築するときに、多くの場合、feval のほうが predict よりも平均応答を予測するのに便利です。ただし、feval は信頼限界を提供しません。

この例は feval メソッドを使用して平均応答を予測する方法を示しています。

  1. 人為的なデータでいくつかの予測子からモデルを作成します。データは X の 2 番目と 3 番目の列を使用しません。2 番目と 3 番目の予測子の効果がほとんどないことは明らかです。関連する予測子を自動的に含むようにモデルをステップワイズに構築します。

    rng('default') % for reproducibility
    X = randn(100,5);
    mu = exp(X(:,[1 4 5])*[2;1;.5]);
    y = poissrnd(mu);
    
    X = array2table(X); % create data table
    y = array2table(y);
    tbl = [X y];
    
    mdl = stepwiseglm(tbl,...
        'constant','upper','linear','Distribution','poisson');
    1. Adding x1, Deviance = 2515.02869, Chi2Stat = 47242.9622, PValue = 0
    2. Adding x4, Deviance = 328.39679, Chi2Stat = 2186.6319, PValue = 0
    3. Adding x5, Deviance = 96.3326, Chi2Stat = 232.0642, PValue = 2.114384e-52
  2. 新しいデータをいくつか生成して、データから予測を評価します。

    Xnew = randn(3,5) + repmat([1 2 3 4 5],[3,1]); % new data
    ynew = feval(mdl,Xnew(:,1),Xnew(:,4),Xnew(:,5)) % only need predictors 1,4,5
    ynew =
    
       1.0e+04 *
    
        0.1130
        1.7375
        3.7471

    同様に、

    ynew = feval(mdl,Xnew(:,[1 4 5])) % only need predictors 1,4,5
    ynew =
    
       1.0e+04 *
    
        0.1130
        1.7375
        3.7471

random

random メソッドは、指定された予測子の値に対して新しいランダムな応答値を生成します。応答値の分布は、モデルで使用される分布です。random は予測子、推定係数およびリンク関数から分布の平均を計算します。正規分布などの分布の場合、モデルは応答の分散の推定値も提供します。二項分布およびポアソン分布の場合、応答の分散を平均によって判別します。random は別個の "分布" 推定値を使用しません。

この例は random メソッドを使用して応答をシミュレートする方法を示しています。

  1. 人為的なデータでいくつかの予測子からモデルを作成します。データは X の 2 番目と 3 番目の列を使用しません。2 番目と 3 番目の予測子の効果がほとんどないことは明らかです。関連する予測子を自動的に含むようにモデルをステップワイズに構築します。

    rng('default') % for reproducibility
    X = randn(100,5);
    mu = exp(X(:,[1 4 5])*[2;1;.5]);
    y = poissrnd(mu);
    mdl = stepwiseglm(X,y,...
        'constant','upper','linear','Distribution','poisson');
    1. Adding x1, Deviance = 2515.02869, Chi2Stat = 47242.9622, PValue = 0
    2. Adding x4, Deviance = 328.39679, Chi2Stat = 2186.6319, PValue = 0
    3. Adding x5, Deviance = 96.3326, Chi2Stat = 232.0642, PValue = 2.114384e-52
  2. 新しいデータをいくつか生成して、データから予測を評価します。

    Xnew = randn(3,5) + repmat([1 2 3 4 5],[3,1]); % new data
    ysim = random(mdl,Xnew)
    ysim =
    
            1111
           17121
           37457

    random による予測はポアソンの標本であり、整数です。

  3. random メソッドを再度評価すると、異なる結果が得られます。

    ysim = random(mdl,Xnew)
    ysim =
    
            1175
           17320
           37126

近似モデルの共有

モデルの表示には、他のユーザーが理論的にモデルを再作成するために必要な情報が含まれています。たとえば、以下のようにします。

rng('default') % for reproducibility
X = randn(100,5);
mu = exp(X(:,[1 4 5])*[2;1;.5]);
y = poissrnd(mu);
mdl = stepwiseglm(X,y,...
    'constant','upper','linear','Distribution','poisson')
1. Adding x1, Deviance = 2515.02869, Chi2Stat = 47242.9622, PValue = 0
2. Adding x4, Deviance = 328.39679, Chi2Stat = 2186.6319, PValue = 0
3. Adding x5, Deviance = 96.3326, Chi2Stat = 232.0642, PValue = 2.114384e-52

mdl = 


Generalized Linear regression model:
    log(y) ~ 1 + x1 + x4 + x5
    Distribution = Poisson

Estimated Coefficients:
                   Estimate    SE          tStat     pValue     
    (Intercept)    0.17604     0.062215    2.8295       0.004662
    x1              1.9122     0.024638    77.614              0
    x4             0.98521     0.026393    37.328    5.6696e-305
    x5             0.61321     0.038435    15.955     2.6473e-57


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

モデル記述にはプログラムを使用してアクセスすることもできます。たとえば、以下のようにします。

mdl.Coefficients.Estimate
ans =

    0.1760
    1.9122
    0.9852
    0.6132
mdl.Formula
ans = 

log(y) ~ 1 + x1 + x4 + x5

参照

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

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

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

[4] Neter, J., M. H. Kutner, C. J. Nachtsheim, and W. Wasserman. Applied Linear Statistical Models, Fourth Edition. Irwin, Chicago, 1996.