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

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

load fisheriris

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

X = meas(51:end,:);

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

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. The 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. Этот результат подразумевает, что отклики в точках с высоким плечом не соответствуют предсказанным значениям из уменьшенной модели.

Предсказание вероятности того, что это версиколор

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

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

newf = 0.4558

newc = 1×2

    0.1234    0.8329

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