Загрузить данные

Загрузите данные радужки Фишера.

load fisheriris

Извлеките строки 51-150, которые имеют классификацию versicolor или virginica.

X = meas(51:end,:);

Создание логических переменных ответа true для versicolor и false для virginica.

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

Подгонка обобщенной линейной модели

Подгонка биномиальной обобщенной линейной модели к данным.

mdl = fitglm(X,y,'linear','Distribution','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

В соответствии с отображением модели, некоторые значения 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

Только четвертый предиктор x3 имеет коэффициент, доверительный интервал которого не включает 0.

Коэффициенты x1 и x2 имеют большие значения p, и их 95% доверительные интервалы включают 0. Проверьте, могут ли оба коэффициента быть нулевыми. Укажите матрицу гипотез для выбора коэффициентов x1 и x2.

M = [0 1 0 0 0     
     0 0 1 0 0];   
p = coefTest(mdl,M)

p = 0.1442

Значение p приблизительно равно 0,14, что не мало. Удалить x1 и x2 из модели.

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

Значение p (pValue) для x2 в таблице коэффициентов больше 0,05, но stepwiseglm включает x2 в модели, поскольку значение p (PValue) для добавления x2 меньше 0,05. stepwiseglm вычисляет функцию PValue использование посадок с и без x2, тогда как функция вычисляет pValue на основе приблизительной стандартной ошибки, вычисленной только из окончательной модели. Поэтому PValue является более надежным, чем pValue.

Определение отклонений

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

plotDiagnostics(mdl2,'leverage')

Наблюдение можно считать отклонением, если его рычаги влияния существенно превышают p/n, где p - количество коэффициентов и n - количество наблюдений. Пунктирная опорная линия является рекомендуемым порогом, вычисляемым по 2*p/n, что соответствует 0,08 на этом сюжете. Некоторые наблюдения имеют значения, превышающие 10*p/n (то есть 0,40). Определите эти точки наблюдения.

idxOutliers = find(mdl2.Diagnostics.Leverage > 10*mdl2.NumCoefficients/mdl2.NumObservations)

idxOutliers = 4×1

    19
    21
    57
    85

Проверьте, изменяются ли коэффициенты модели при подгонке модели, исключая эти точки.

oldCoeffs = mdl2.Coefficients.Estimate;
mdl3 = fitglm(X,y,'linear','Distribution','binomial', ...
    'PredictorVars',2:4,'Exclude',idxOutliers);
newCoeffs = mdl3.Coefficients.Estimate;
disp([oldCoeffs newCoeffs])

   50.5268   44.0085
    8.3761    5.6361
   -7.8745   -6.1145
  -21.4296  -18.1236

Коэффициенты модели в mdl3 отличаются от таковых в mdl2. Этот результат подразумевает, что реакции в точках с высоким рычагом не согласуются с прогнозируемыми значениями из уменьшенной модели.

Прогнозирование вероятности существования Versicolor

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

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

newf = 0.4558

newc = 1×2

    0.1234    0.8329

Модель даёт почти 46% вероятность того, что средний цветок версиколор, с широким доверительным интервалом.