Обобщенный линейный рабочий процесс модели

В этом примере показано, как подбирать обобщенную линейную модель и анализировать результаты. Типичный рабочий процесс включает следующее: импортируйте данные, подбирайте обобщенную линейную модель, протестируйте ее качество, измените их, чтобы улучшить качество и сделать прогнозы на основе модели. Это вычисляет вероятность, что цветок находится в одном из двух классов, на основе ирисовых данных Фишера.

Шаг 1. Загрузите данные.

Загрузите ирисовые данные Фишера. Извлеките строки, которые имеют классификацию versicolor или virginica. Это строки 51 - 150. Создайте логические переменные отклика, которые являются true для versicolor цветов.

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. Исследуйте результат, рассмотрите альтернативные модели.

Некоторые $p$ - значения в pValue столбец очень не мал. Возможно, модель может быть упрощена.

Смотрите, включают ли приблизительно 95% доверительных интервалов для коэффициентов 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.

Коэффициенты 'x1' и 'x2' имейте самое большое $p$ Значения. Протестируйте, могли ли оба коэффициента быть нулем.

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

$p$ - значение приблизительно 0,14 очень не мало. Исключите те условия из модели.

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 включенный 'x2' в модели, потому что это не добавляет и не удаляет условия с $p$ - значения между 0,05 и 0.10.

Шаг 4. Ищите выбросы и исключите их.

Исследуйте график рычагов искать влиятельные выбросы.

plotDiagnostics(mdl2,'leverage')

Существует одно наблюдение с рычагами близко к одному. Используя Data Cursor, кликните по точке и найдите, что это имеет индекс 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.

Используйте mdl2 предсказать вероятность, что цветок со средними измерениями является versicolor. Сгенерируйте доверительные интервалы для своего прогноза.

[newf,newc] = predict(mdl2,mean(X))

newf = 0.5086

newc = 1×2

    0.1863    0.8239

Модель дает почти 50%-ю вероятность, что средний цветок является versicolor с широким доверительным интервалом об этой оценке.

Документация