一般化線形モデル
一般化線形モデルとは
線形回帰モデルは、応答と 1 つ以上の予測項間の線形関係を記述します。しかし、非線形関係が存在します。非線形回帰では、一般の非線形モデルを説明します。"一般化線形モデル" と呼ばれる非線形モデルの特別なクラスは、線形手法を使用します。
線形モデルが次の特性をもつことを思い出してください。
予測子の値のセットごとに、応答は平均が μ の正規分布になります。
係数ベクトル b は、予測子 X の線形結合 Xb を定義します。
モデルは μ = Xb です。
一般化線形モデルでは、これらの特性は次のように一般化されます。
データの準備
回帰の当てはめを開始するには、データを近似関数に望ましい形式にします。すべての回帰手法は、配列 X
の入力データと独立したベクトル y
の応答データか、テーブルまたはデータセット配列 tbl
内の入力データと tbl
の列としての応答データで始まります。入力データの各行が、1 つの観測値を表します。各列が 1 つの予測子 (変数) を表します。
テーブルまたはデータセット配列 tbl
では、次のように 'ResponseVar'
の名前と値のペアで応答変数を示します。
mdl = fitglm(tbl,'ResponseVar','BloodPressure');
応答変数は既定で最後の列です。
数値の "カテゴリカル" 予測子を使用できます。カテゴリカル予測子は可能性のある固定セットからの値をとります。
数値配列
X
では、'Categorical'
の名前と値のペアでカテゴリカル予測子を示します。たとえば、6 つの予測子から2
と3
が 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
の列を指定します。
入力および応答データのテーブル
ワークスペース変数からテーブルを作成するには、次のようにします。
load carsmall
tbl = table(MPG,Weight);
tbl.Year = ordinal(Model_Year);
入力データの数値配列、応答の数値ベクトル
たとえば、ワークスペース変数から数値配列を作成するには、次のようにします。
load carsmall
X = [Weight Horsepower Cylinders Model_Year];
y = MPG;
Excel® スプレッドシートから数値配列を作成するには、次のようにします。
[X, Xnames] = readmatrix('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 |
| µ = Xb |
| log(µ) = Xb |
| log(µ/(1 – µ)) = Xb |
| log(–log(µ)) = Xb |
'probit' | Φ–1(µ) = Xb、Φ は正規 (ガウス) 累積分布関数 |
'reciprocal' 、分布の既定の設定 'gamma' | µ–1 = Xb |
| µp = Xb |
| ユーザー指定のリンク関数 (カスタム リンク関数を参照) |
既定ではないリンク関数は、主に二項モデルで使用すると便利です。これらの既定ではないリンク関数は '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' | モデルは、切片、線形項、交互作用、二乗項を含みます。 |
'poly | モデルは多項式であり、最初の予測子は次数 i まで、2 番目の予測子は次数 j まで、3 番目以降も同様です。0 から 9 までの数値を使用します。たとえば、'poly2111' には、1 つの定数のほかにすべての線形項と積項があり、また、予測子 1 の二乗の項を含んでいます。 |
項の行列
項行列 T
は、モデル内の項を指定する t 行 (p + 1) 列の行列です。t は項の数、p は予測子変数の数であり、+1 は応答変数に相当します。T(i,j)
の値は、項 i
の変数 j
の指数です。
たとえば、3つの予測子変数 x1
、x2
、x3
と応答変数 y
が x1
、x2
、x3
、y
という順序で入力に含まれていると仮定します。T
の各行は 1 つの項を表します。
[0 0 0 0]
— 定数項 (切片)[0 1 0 0]
—x2
(x1^0 * x2^1 * x3^0
と等価)[1 0 1 0]
—x1*x3
[2 0 0 0]
—x1^2
[0 1 2 0]
—x2*(x3^2)
各項の最後の 0
は、応答変数を表します。一般に、項行列内のゼロの列ベクトルは、応答変数の位置を表します。行列と列ベクトルに予測子と応答変数がある場合、各行の最後の列に応答変数を示す 0
を含めなければなりません。
式
モデル仕様の式は、次のような形式の文字ベクトルまたは string スカラーです。
'
,y
~ terms
'
y
は応答名です。terms
次が含まれます変数名
+
次の変数を含みます-
次の変数を除外します:
項の積である交互作用を定義します*
交互作用とすべての次数の低い項を定義します*
で繰り返されるとおり、^
は予測子をべき乗にし、^
は低い次数の項も含みます()
項をグループ化します
ヒント
式には既定で定数 (切片) 項が含まれます。モデルから定数項を除外するには、式に -1
を含めます。
次に例を示します。
'y ~ x1 + x2 + x3'
は、切片がある 3 変数線形モデルです。
'y ~ x1 + x2 + x3 - 1'
は、切片がない 3 変数線形モデルです。
'y ~ x1 + x2 + x3 + x2^2'
は、切片と x2^2
項がある 3 変数モデルです。
'y ~ x1 + x2^2 + x3'
は、x2^2
に x2
項が含まれるので、前の例と同じです。
'y ~ x1 + x2 + x3 + x1:x2'
には x1*x2
項が含まれています。
'y ~ x1*x2 + x3'
は、x1*x2 = x1 + x2 + x1:x2
なので、前の例と同じです。
'y ~ x1*x2*x3 - x1:x2:x3'
には、3 次交互作用を除く x1
、x2
、x3
間の交互作用がすべてあります。
'y ~ x1*(x2 + x3 + x4)'
には、すべての線形項、および x1
と他の各変数の積があります。
モデルのデータへの当てはめ
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
の関係を以下のように定義します。
次の例では、二項モデルをデータに当てはめます。このデータは、重量が異なる大型車の測定が含まれる carbig.mat
から導出されています。w
の各重量には、total
に対応した自動車数と poor
に対応した燃費の悪い自動車数が含まれています。
total
および w
に依存する成功のパーセンテージでトライアル数を指定すると、poor
の値が二項分布分布に続くと仮定できます。この分布は、リンク関数 で一般化線形モデルを使用することにより、ロジスティック モデルのコンテキスト内で説明されます。このリンク関数は "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')
以下のことが明らかになりました。残差は正規分布していません。代わりに、潜在的なポアソン分布のように裾がより厚くなっています。
予測子の効果とモデルの変更方法を理解するためのプロット
次の例は、各予測子が回帰モデルに与える効果を理解する方法と、モデルを変更して不要な項を削除する方法を示します。
人為的なデータでいくつかの予測子からモデルを作成します。データは 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");
応答のスライス プロットを調査します。これは、各予測子の効果を個別に表示します。
plotSlice(mdl)
最初の予測子のスケールは、プロットを圧倒しています。最初の予測子を無効にするには、[予測子] リストで 2 番目から 5 番目までの予測子のみを選択し、[適用] をクリックします。
2 番目と 3 番目の予測子の効果がほとんどないことは明らかです。青い縦の破線によって表される個々の予測子値をドラッグできます。また、青い影付きの領域で表されている同時信頼限界と非同時信頼限界の間を選択することもできます。予測子のラインをドラッグして、2 番目と 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
メソッドを使用して予測で信頼区間を予測して取得する方法を示します。
人為的なデータでいくつかの予測子からモデルを作成します。データは
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
新しいデータをいくつか生成して、データから予測を評価します。
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
メソッドを使用して平均応答を予測する方法を示しています。
人為的なデータでいくつかの予測子からモデルを作成します。データは
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
新しいデータをいくつか生成して、データから予測を評価します。
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
メソッドを使用して応答をシミュレートする方法を示しています。
人為的なデータでいくつかの予測子からモデルを作成します。データは
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
新しいデータをいくつか生成して、データから予測を評価します。
Xnew = randn(3,5) + repmat([1 2 3 4 5],[3,1]); % new data ysim = random(mdl,Xnew)
ysim = 1111 17121 37457
random
による予測はポアソンの標本であり、整数です。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.