ドキュメンテーション

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

一般化線形モデルのワークフロー

この例では、一般化線形モデルの近似方法と結果の解析方法について説明します。典型的なワークフローは、データのインポート、一般化線形モデルの近似、品質テスト、品質改善のための変更、そしてモデルに基づく予測という流れです。フィッシャーのアヤメのデータをもとに、花が 2 つのクラスのいずれかである確率を計算します。

手順 1. データを読み込む。

フィッシャーのアヤメのデータを読み込みます。versicolor または virginica の分類をもつ行を抽出します。これに該当するのは 51 行目から 150 行目です。versicolor の花の場合に true である論理応答変数を作成します。

load fisheriris
X = meas(51:end,:); % versicolor and virginica
y = strcmp('versicolor',species(51:end));

手順 2. 一般化線形モデルで近似する。

二項一般化線形モデルをデータにあてはめます。

mdl = fitglm(X,y,'linear',...
    'distr','binomial')
mdl = 
Generalized linear regression model:
    logit(y) ~ 1 + x1 + x2 + x3 + x4
    Distribution = Binomial

Estimated Coefficients:
                   Estimate      SE       tStat      pValue 
                   ________    ______    _______    ________

    (Intercept)     42.638     25.708     1.6586    0.097204
    x1              2.4652     2.3943     1.0296     0.30319
    x2              6.6809     4.4796     1.4914     0.13585
    x3             -9.4294     4.7372    -1.9905    0.046537
    x4             -18.286     9.7426    -1.8769    0.060529


100 observations, 95 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 127, p-value = 1.95e-26

手順 3. 結果を調べ、別のモデルを検討する。

pValue 列には、あまり小さくない p 値が含まれています。おそらく、このモデルは簡素化できます。

係数の中で 95% の信頼区間に 0 が含まれているものがあるかどうかを確認します。0 が含まれている場合、おそらくその項をモデルから削除できます。

confint = coefCI(mdl)
confint = 5×2

   -8.3984   93.6740
   -2.2881    7.2185
   -2.2122   15.5739
  -18.8339   -0.0248
  -37.6277    1.0554

信頼区間に 0 を含まない係数をもつ予測子は 2 つのみです。

'x1''x2' の係数は p 値が最大になっています。両方の係数を 0 にできるかどうかを検定します。

M = [0 1 0 0 0     % picks out coefficient for column 1
     0 0 1 0 0];   % picks out coefficient for column 2
p = coefTest(mdl,M)
p = 0.1442

約 0.14 という p 値はあまり小さい値ではありません。これらの項をモデルから削除してください。

mdl1 = removeTerms(mdl,'x1 + x2')
mdl1 = 
Generalized linear regression model:
    logit(y) ~ 1 + x3 + x4
    Distribution = Binomial

Estimated Coefficients:
                   Estimate      SE       tStat       pValue  
                   ________    ______    _______    __________

    (Intercept)     45.272     13.612      3.326    0.00088103
    x3             -5.7545     2.3059    -2.4956      0.012576
    x4             -10.447     3.7557    -2.7816     0.0054092


100 observations, 97 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 118, p-value = 2.3e-26

おそらく、stepwiseglm を使用して最初にモデルを特定したほうが、より良い結果を得られたでしょう。

mdl2 = stepwiseglm(X,y,...
    'constant','Distribution','binomial','upper','linear')
1. Adding x4, Deviance = 33.4208, Chi2Stat = 105.2086, PValue = 1.099298e-24
2. Adding x3, Deviance = 20.5635, Chi2Stat = 12.8573, PValue = 0.000336166
3. Adding x2, Deviance = 13.2658, Chi2Stat = 7.29767, PValue = 0.00690441
mdl2 = 
Generalized linear regression model:
    logit(y) ~ 1 + x2 + x3 + x4
    Distribution = Binomial

Estimated Coefficients:
                   Estimate      SE       tStat      pValue 
                   ________    ______    _______    ________

    (Intercept)     50.527     23.995     2.1057    0.035227
    x2              8.3761     4.7612     1.7592    0.078536
    x3             -7.8745     3.8407    -2.0503    0.040334
    x4              -21.43     10.707    -2.0014     0.04535


100 observations, 96 error degrees of freedom
Dispersion: 1
Chi^2-statistic vs. constant model: 125, p-value = 5.4e-27

stepwiseglm では p 値が 0.05 と 0.10 の間にある項を追加も削除もしないので、'x2' がモデルに含まれています。

手順 4. 外れ値を探して除外する。

てこ比のプロットを調べて、影響のある外れ値を探します。

plotDiagnostics(mdl2,'leverage')

てこ比が 1 に近い観測が 1 つあります。データ カーソルを使用して点をクリックし、インデックス 69 の観測を見つけます。

この点を削除したモデルを近似したときに、モデル係数が変わるかどうかを確認します。

oldCoeffs = mdl2.Coefficients.Estimate;
mdl3 = fitglm(X,y,'linear',...
    'distr','binomial','pred',2:4,'exclude',69);
newCoeffs = mdl3.Coefficients.Estimate;
disp([oldCoeffs newCoeffs])
   50.5268   50.5268
    8.3761    8.3761
   -7.8745   -7.8745
  -21.4296  -21.4296

モデル係数は変化せず、てこ比の高い点の応答が、次元削減されたモデルでの予測値と整合していることを示します。

手順 5. 新しい花が versicolor である確率を予測する。

平均的な測定値の花が versicolor である確率を予測するには、mdl2 を使用します。この予測に対する信頼区間を生成します。

[newf,newc] = predict(mdl2,mean(X))
newf = 0.5086
newc = 1×2

    0.1863    0.8239

このモデルは、この推定に関する広い信頼区間もち、平均的な花が vesicolor である確率を約 50% としています。