Загрузка данных

Загрузите ирисовые данные Фишера.

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%-ю вероятность, что средний цветок является versicolor с широким доверительным интервалом.