В этом примере показано, как подбирать обобщенную линейную модель и анализировать результаты. Типичный рабочий процесс включает эти шаги: импортируйте данные, подбирайте обобщенную линейную модель, протестируйте ее качество, измените модель, чтобы улучшить ее качество и сделать предсказания на основе модели. В этом примере вы используете ирисовые данные Фишера, чтобы вычислить вероятность, что цветок находится в одном из двух классов.
Загрузите ирисовые данные Фишера.
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. Этот результат подразумевает, что ответы в точках высоких рычагов не сопоставимы с ожидаемыми значениями от упрощенной модели.
Используйте mdl3
предсказать вероятность, что цветок со средними измерениями является versicolor. Сгенерируйте доверительные интервалы для предсказания.
[newf,newc] = predict(mdl3,mean(X))
newf = 0.4558
newc = 1×2
0.1234 0.8329
Модель дает почти 46%-ю вероятность, что средний цветок является versicolor с широким доверительным интервалом.
coefCI
| fitglm
| GeneralizedLinearModel
| plotDiagnostics
| predict
| removeTerms
| stepwiseglm