Матрунич Консалтинг

Анализ данных, визуализация, маркетинговые исследования, язык R

Tag: machine learning

Визуальная диагностика простой линейной регрессии с ggvis

Продолжаю тему с однофакторым регрессионным анализом в R. Теперь мне интересно применить инструменты визуализации, доступные в R-библиотеке ggvis, для диагностики регрессионной модели. Здесь же мы применяем перенаправление с помощью функции %>%.

Напомню, что я использую данные, сопровожающие книгу An Introduction to Statistical Learning with Applications in R. Подготовка данных и построение модели было произведено в предыдущей публикации. Так что будем считать, что мы продолжаем творить в пространстве той публикации.

Нелинейность взаимосвязи

data.frame(predict = predict(autolm),
           resids  = residuals(autolm)) %>%
  ggvis(~predict, ~resids) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения") %>%
  add_axis("y", title = "Остатки") %>%
  layer_smooths(stroke := 'red')
Точечный график со сглаживанием в ggvis

Точечный график со сглаживанием в ggvis

Наблюдается некоторая нелинейная зависимость между остатками и предсказанными значениями. Имеет смысл сделать нелинейное преобразование элементов модели (логарифм, квадратный корень, квадрат).

Корреляция остатков

data.frame(id = seq_len(dim(Auto)[1]),
           resids  = residuals(autolm)) %>%
  ggvis(~id, ~resids) %>% 
  layer_paths() %>%
  add_axis("x", title = "Порядковый номер наблюдения") %>%
  add_axis("y", title = "Остатки")
Визуальная оценка взаимосвязи остатков

Визуальная оценка взаимосвязи остатков

График не выглядит таким, что по нему на основе предыдущего наблюдения можно предсказать следующее. Делаем вывод, что корреляции остатков нет.

Непостоянная дисперсия остатков

data.frame(predict = predict(autolm),
           rstuds  = rstudent(autolm)) %>%
  ggvis(~predict, ~rstuds) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения") %>%
  add_axis("y", title = "Стьюдентизированные остатки")
Стьюдентизированные остатки и предсказанные значения: точечный график в ggvis

Стьюдентизированные остатки и предсказанные значения: точечный график в ggvis

Видно, что с увеличением величины предсказанных значений происходит увеличение разброса остатков: мы наблюдаем изогнутую воронкообразную форму. Это говорит о наличии гетероскедастичности. Отсутствие гетероскедастичности — условие для корректного расчёта стандартной ошибки, доверительных интервалов и тестирования гипотез. В нашем случае можно преобразовать зависимую переменную с помощью какой-нибудь “выгнутой” функции — например, логарифм или корень квадратный.

autolm1 <- update(autolm, log(.) ~ .)
data.frame(predict = predict(autolm1),
           rstuds  = rstudent(autolm1)) %>%
  ggvis(~predict, ~rstuds) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения (log)") %>%
  add_axis("y", title = "Стьюдентизированные остатки")
Точечный график стьюдентизированных остатков и натурального логарифма предсказанных значений

Точечный график стьюдентизированных остатков и натурального логарифма предсказанных значений

Выбросы

Auto %>%
  mutate(rstud = ifelse(abs(rstudent(autolm)) > 3, 'Больше 3', 'Не больше 3'),
         label = row.names(Auto),
         label = ifelse(rstud == 'Больше 3', label, '')) %>% 
  ggvis(~kpl, ~ horsepower, fill = ~rstud) %>%
  layer_points() %>%
  add_axis("x", title = "Километров на литр") %>%
  add_axis("y", title = "Лошадиных сил") %>%
  add_legend("fill", title = "Стьюд.остатки") %>%
  layer_model_predictions(stroke := "red", model = 'lm') %>%
  layer_text(text := ~label, dy := -5, dx := 5)
Визуальное определение предикатов-выбросов

Визуальное определение предикатов-выбросов

Только у двух наблюдений абсолютная величина стьюдентизированных остатков превышает 3.

Точки с избыточным влиянием на модель (High Leverage Points)

В простой (однофакторной) регрессии мы просто можем посмотреть, какие значения независимой переменной слишком выпадают. Но построим график с оценкой уровня влияния точек.

data.frame(leverage = hatvalues(autolm),
           rstuds  = rstudent(autolm),
           label = row.names(Auto)) %>%
  mutate(label = ifelse(leverage > .025, label, '')) %>% 
  ggvis(~leverage, ~rstuds) %>% 
  layer_points() %>%
  add_axis("x", title = "Leverage") %>%
  add_axis("y", title = "Стьюдентизированные остатки") %>%
  layer_text(text := ~label, dy := -5, dx := 5)
Точечный график с оценкой влияния наблюдений на модель

Точечный график с оценкой влияния наблюдений на модель

Обогащённая графика с диагностикой модели

data.frame(predict = predict(autolm),
           resids  = residuals(autolm),
           rstuds  = rstudent(autolm),
           leverage = hatvalues(autolm)) %>%
  ggvis(~predict, ~resids, fill = ~abs(rstuds), size = ~leverage) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения") %>%
  add_axis("y", title = "Остатки")  %>%
  add_legend("size", title = "Влияние на модель") %>% 
  add_legend("fill", title = "Модуль стьюд.ост.", 
             properties = legend_props(legend = list(y = 100)))
График, обобщающий оценку модели: выбросы, влияние отдельных значений, нелинейность взаимосвязи

График, обобщающий оценку модели: выбросы, влияние отдельных значений, нелинейность взаимосвязи

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

Отчётливо видно, что наблюдения с наибольшим влиянием встречаются в том секторе, где наблюдения встречаются реже.

Однофакторный линейный регрессионный анализ в R

На основе книги An Introduction to Statistical Learning with Applications in R авторов Gareth James, Daniela Witten, Trevor Hastie и Robert Tibshirani. Бесплатная копия книги в PDF-формате, или заказ бумажной копии книги на Amazon.

library(ISLR)
library(ggvis)
library(dplyr)

Подготовка данных

Работаем с массивом Auto из пакета ISLR.

# Загружаем массив Auto в пространство пользователя
data(Auto)

# Структура массива
str(Auto)
## 'data.frame':    392 obs. of  9 variables:
##  $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
##  $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
##  $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
##  $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
##  $ weight      : num  3504 3693 3436 3433 3449 ...
##  $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
##  $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
##  $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
# Первые строки массива
head(Auto)
##   mpg cylinders displacement horsepower weight acceleration year origin
## 1  18         8          307        130   3504         12.0   70      1
## 2  15         8          350        165   3693         11.5   70      1
## 3  18         8          318        150   3436         11.0   70      1
## 4  16         8          304        150   3433         12.0   70      1
## 5  17         8          302        140   3449         10.5   70      1
## 6  15         8          429        198   4341         10.0   70      1
##                        name
## 1 chevrolet chevelle malibu
## 2         buick skylark 320
## 3        plymouth satellite
## 4             amc rebel sst
## 5               ford torino
## 6          ford galaxie 500
# Конвертируем миль (1.609 км) на галлон (3.79 л) в км на литр
Auto$kpl <- Auto$mpg * 1.609 / 3.79

Построение модели

Построим однофакторную регрессионную модель, чтобы выявлить, как мощность двигателя влияет на пробег на километр топлива, и ответим на следующие вопросы:

  1. Существует ли взаимосвязь между мощностью и пробегом.
  2. Насколько сильна взаимосвязь между мощностью и пробегом.
  3. Взаимосвязь между мощностью и пробегом прямая или обратная.
  4. Каково прогнозное значение пробега на километр при мощности 98 лошадиных сил. Каковы 95% доверительный интервал и интервал предсказания.
# Задаем модель
autolm <- lm(kpl ~ horsepower, data = Auto)

# Результаты анализа
summary(autolm)
## 
## Call:
## lm(formula = kpl ~ horsepower, data = Auto)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.761 -1.384 -0.146  1.173  7.185 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.95430    0.30461    55.7   <2e-16 ***
## horsepower  -0.06701    0.00274   -24.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.08 on 390 degrees of freedom
## Multiple R-squared:  0.606,  Adjusted R-squared:  0.605 
## F-statistic:  600 on 1 and 390 DF,  p-value: <2e-16

Коэффициент horsepower составляет -0.067, т.е. при увеличении мощности автомобиля на одну лошадиную силу пробег снижается на 0.07 км на литр топлива. Наличие этой взаимосвязи подтверждается p-значением t-статистики: t равен -24.5, а вероятность получить такое значение t при независимых переменных (т.е. p-значение) стремится к нулю (<2e-16).

Прогнозирование

# Прогнозируем пробег при мощности 98 л.с.
# Доверительный интервал
confint <- predict(autolm,
                   newdata = data.frame(horsepower = 98),
                   interval = 'confidence')
confint
##     fit   lwr  upr
## 1 10.39 10.18 10.6

Если мы возьмём несколько разных автомобилей с мощностью двигателя 98 л.с., то средний пробег этих машин на литр топлива с вероятностью 95% попадёт в интервал от 10.2 до 10.6 км.

# Интервал предсказания
predint <- predict(autolm,
                   newdata = data.frame(horsepower = 98),
                   interval = 'prediction')
predint
##     fit   lwr   upr
## 1 10.39 6.287 14.49

Если мы возьмём какой-то конкретный автомобиль с мощностью двигателя 98 л.с., то его пробег на литр топлива с вероятностью 95% попадёт в интервал от 6.3 до 14.5 км.

Auto %>% 
  ggvis(~horsepower, ~kpl) %>%
  layer_points() %>%
  layer_model_predictions(model = 'lm', se = T)

Точечный график и линия тренда

Диагностика

После построения регрессионной модели имеет смысл проверить наличие следующих проблем:

  1. Нелинейная взаимосвязь между независимыми и зависимой переменными.
  2. Корреляция остатков между собой.
  3. Непостоянная дисперсия остатков (гетероскедастичность).
  4. Наличие выбросов.
  5. Избыточное влияние отдельных наблюдений.
  6. Коллинеарность предикатов.
par(mfrow = c(2, 3))
plot(autolm, which = 1:6, caption = 
       list("Прогноз и остатки", "Квантили", "Разброс-Положение",
            "Расстояние Кука", "Влияние и остатки", 
            expression("Влияние и Кук " * h[ii] / (1 - h[ii]))))
Визуальная диагностика регрессионной модели

Визуальная диагностика регрессионной модели

Точечный график предсказанных значений и остатков позволяет обнаружить закономерности распределения последних. График квантилей помогает определить нормальность остатков. График “Разброс-положение” даёт возможность оценить монотонность распределения остатков (гетероскедастичность). С помощью графика “Остатки-влияние” можно определить точки с большим влиянием на модель.

Powered by WordPress & Theme by Anders Norén