Рабочий процесс линейной регрессии

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

Шаг 1. Импортируйте данные в таблицу.

hospital.xls - электронная таблица Excel ®, содержащая имена пациентов, пол, возраст, вес, артериальное давление и даты лечения в экспериментальном протоколе. Сначала считайте данные в таблицу .

patients = readtable('hospital.xls','ReadRowNames',true);

Исследуйте пять строк данных.

patients(1:5,:)
ans=5×11 table
                   name         sex     age    wgt    smoke    sys    dia    trial1    trial2    trial3    trial4
               ____________    _____    ___    ___    _____    ___    ___    ______    ______    ______    ______

    YPL-320    {'SMITH'   }    {'m'}    38     176      1      124    93       18       -99       -99       -99  
    GLI-532    {'JOHNSON' }    {'m'}    43     163      0      109    77       11        13        22       -99  
    PNI-258    {'WILLIAMS'}    {'f'}    38     131      0      125    83      -99       -99       -99       -99  
    MIJ-579    {'JONES'   }    {'f'}    40     133      0      117    75        6        12       -99       -99  
    XLK-030    {'BROWN'   }    {'f'}    49     119      0      122    80       14        23       -99       -99  

The sex и smoke поля, по-видимому, имеют два выбора каждый. Поэтому смените эти поля на категориальные.

patients.smoke = categorical(patients.smoke,0:1,{'No','Yes'});
patients.sex = categorical(patients.sex);

Шаг 2. Создайте подобранную модель.

Ваша цель - смоделировать систематическое давление как функцию возраста, веса, пола и статуса курения пациента. Создайте линейную формулу для 'sys' как функцию 'age', 'wgt', 'sex', и 'smoke' .

modelspec = 'sys ~ age + wgt + sex + smoke';
mdl = fitlm(patients,modelspec)
mdl = 
Linear regression model:
    sys ~ 1 + sex + age + wgt + smoke

Estimated Coefficients:
                   Estimate        SE        tStat        pValue  
                   _________    ________    ________    __________

    (Intercept)       118.28      7.6291      15.504    9.1557e-28
    sex_m            0.88162      2.9473     0.29913       0.76549
    age              0.08602     0.06731       1.278       0.20438
    wgt            -0.016685    0.055714    -0.29947       0.76524
    smoke_Yes          9.884      1.0406       9.498    1.9546e-15


Number of observations: 100, Error degrees of freedom: 95
Root Mean Squared Error: 4.81
R-squared: 0.508,  Adjusted R-Squared: 0.487
F-statistic vs. constant model: 24.5, p-value = 5.99e-14

Предикторы пола, возраста и веса имеют довольно высокие p- значения, указывающие, что некоторые из этих предикторов могут быть ненужными.

Шаг 3. Найдите и удалите выбросы.

Проверьте, есть ли выбросы в данных, которые должны быть исключены из подгонки. Постройте график невязок.

plotResiduals(mdl)

Figure contains an axes. The axes with title Histogram of residuals contains an object of type patch.

Существует один возможный выброс со значением более 12. Вероятно, это не действительно выбросы. Для демонстрации вот как ее найти и удалить.

Найдите выбросы.

outlier = mdl.Residuals.Raw > 12;
find(outlier)
ans = 84

Удалите выбросы.

mdl = fitlm(patients,modelspec,...
    'Exclude',84);

mdl.ObservationInfo(84,:)
ans=1×4 table
               Weights    Excluded    Missing    Subset
               _______    ________    _______    ______

    WXM-486       1        true        false     false 

Наблюдение 84 больше не находится в модели.

Шаг 4. Упростите модель.

Попытайтесь получить более простую модель, с меньшим количеством предикторов, но с той же прогнозирующей точностью. step ищет лучшую модель путем добавления или удаления одного термина за раз. Разрешить step сделать до 10 шагов.

mdl1 = step(mdl,'NSteps',10)
1. Removing wgt, FStat = 4.6001e-05, pValue = 0.9946
2. Removing sex, FStat = 0.063241, pValue = 0.80199
mdl1 = 
Linear regression model:
    sys ~ 1 + age + smoke

Estimated Coefficients:
                   Estimate       SE       tStat       pValue  
                   ________    ________    ______    __________

    (Intercept)     115.11       2.5364    45.383    1.1407e-66
    age            0.10782     0.064844    1.6628       0.09962
    smoke_Yes       10.054      0.97696    10.291    3.5276e-17


Number of observations: 99, Error degrees of freedom: 96
Root Mean Squared Error: 4.61
R-squared: 0.536,  Adjusted R-Squared: 0.526
F-statistic vs. constant model: 55.4, p-value = 1.02e-16

step сделал два шага. Это означает, что он не мог улучшить модель дальше путем добавления или вычитания одного термина.

Постройте график эффективности более простой модели на обучающих данных.

plotResiduals(mdl1)

Figure contains an axes. The axes with title Histogram of residuals contains an object of type patch.

Невязки выглядят примерно такими же маленькими, как и у исходной модели.

Шаг 5. Прогнозируйте ответы на новые данные.

Предположим, у вас есть четыре новых человека в возрасте 25, 30, 40 и 65, и первый и третий дым. Предсказать их сестрическое давление можно используя mdl1.

ages = [25;30;40;65];
smoker = {'Yes';'No';'Yes';'No'};
systolicnew = feval(mdl1,ages,smoker)
systolicnew = 4×1

  127.8561
  118.3412
  129.4734
  122.1149

Чтобы делать предсказания, вам нужны только переменные, которые mdl1 использует.

Шаг 6. Поделитесь моделью.

Вы можете захотеть, чтобы другие могли использовать вашу модель для предсказания. Доступ к терминам в линейной модели.

coefnames = mdl1.CoefficientNames
coefnames = 1x3 cell
    {'(Intercept)'}    {'age'}    {'smoke_Yes'}

Просмотрите формулу модели.

mdl1.Formula
ans = 
sys ~ 1 + age + smoke

Доступ к коэффициентам членов.

coefvals = mdl1.Coefficients(:,1).Estimate
coefvals = 3×1

  115.1066
    0.1078
   10.0540

Модель sys = 115.1066 + 0.1078*age + 10.0540*smoke, где smoke является 1 для курильщика и 0 в противном случае.

См. также

| | | |

Похожие темы