Введение

В файле assumptions.xlsx содержатся данные переменных x, y1, y2.

library(readxl) # считывание данных из файлов Excel
library(tidyverse) # манипулирование данными 
library(ggplot2) # визуализация 
library(ggfortify) # визуализация диагностических графиков
library(modelr) # вспомогательные функции для работы с моделями
library(broom) # преобразование результатов моделирования в табличный вид
library(GGally) # построение матрицы диаграмм рассеяния 
library(car) # функции для степенных преобразований

Задания

Задение 1.

Загрузите данные и постройте визуализации зависимостей y1(x) и y2(x).

assumptions <- read_excel('data/assumptions.xlsx')
ggplot(assumptions, aes(x, y1)) + 
  geom_point() + 
  geom_smooth(method = lm, se = F, col = 'blue') +
  labs(title = 'Визуализация зависимости y1(x)')

ggplot(assumptions, aes(x, y2)) + 
  geom_point() + 
  geom_smooth(method = lm, se = F, col = 'blue') +
  labs(title = 'Визуализация зависимости y2(x)')

Задение 2.

Постройте линейные модели y1(x) и y2(x) и проверьте, выполняются ли допущения линейной модели.

Построение модели y1(x) и проверка на выполнение допущений линейной модели:

m1 <- lm(y1 ~ x, assumptions) 
autoplot(m1)

Первый график показывает зависимость остатков модели от предсказанных значений: Если допущения линейной регрессии верны, то не должно быть зависимости остатков от предскзанных значений. Искривленная форма облака точек на этом графике означает, что зависимость между x и y1нелинейная и необходимо либо включать в модель нелинейные члены, либо использовать линеаризующее преобразование данных. На этом графике прослеживается непостоянство дисперсии остатков. Можно предполодить, что результаты проверки гипотез о значимости модели и коэффициентов будут некорректными.

Второй график позволяет проверить предположение о нормальности остатков: Остатки на графики близки к прямой линии, что свидетельствует о справедливости допущений.

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

Четвертый график помогает выявить выбросы и влиятельные наблюдения.

Модель плохая: гетероскедастичность.

Построение модели y2(x) и проверка на выполнение допущений линейной модели:

m2 <- lm(y2 ~ x, assumptions) 
autoplot(m2)

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

Второй график позволяет проверить предположение о нормальности остатков: отклонение от нормального распределения в области положительных остатков более заметно.

Третий график позволяет обнаружить гетероскедастичность: присутствует изменение дисперсии остатков.

Задение 3.

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

Рассмотрим первую модель:

t1 <- powerTransform(y1 ~ x, assumptions) # Эта функция подбирает параметр λ таким образом, чтобы распределение остатков линейной модели в результате было наиболее близко к нормальному.

m1_2 <- lm(y1^-0.2 ~ x, assumptions)
coef(m1_2)
##  (Intercept)            x 
##  0.611124375 -0.004365662

Уравнение регрессии: \(\frac{1}{\sqrt[5]{y1}} = 0.611−0.0044*x\)

Диагностика остатков:

autoplot(m1_2)

Рассмотрим вторую модель:

t2 <- powerTransform(y2 ~ x, assumptions) # Эта функция подбирает параметр λ таким образом, чтобы распределение остатков линейной модели в результате было наиболее близко к нормальному.

m2_2 <- lm(y2^-0.143 ~ x, assumptions)
coef(m2_2)
##  (Intercept)            x 
##  0.600139449 -0.001116676

Уравнение регрессии: \(\frac{1}{\sqrt[7]{y2}} = 0.6−0.0011*x\)

Диагностика остатков:

autoplot(m2_2)

К сожалению, устранить проблемы в полной мере не удалось.

Анализ влиятельных наблюдений:

autoplot(m1_2, which = c(1, 4:6))

autoplot(m2_2, which = c(1, 4:6))

Для модели y2(x) влиятельные наблюдения являются более значительными и выделяющимися на фоне всего ряда данных. Связано это может быть с тем, что в модели y2(x) проявляется большая корреляция.

Задение 4.

Постройте прогноз и доверительные интервалы прогноза на основе полученных моделей.

# Создаем точки для получения нового прогноза
newdata_ht_1 <- data.frame(x = seq_range(assumptions$x, by = 0.1))

# Интервальный прогноз:
pred_ht_transformed_1 <- predict(m1_2, newdata = newdata_ht_1, 
                   interval = 'prediction',
                   level = .95)

# Обратное преобразование:
pred_ht_1 <- pred_ht_transformed_1^(-5)

# Соединяем все в одну таблицу:
newdata_ht_1 <- bind_cols(newdata_ht_1, as.data.frame(pred_ht_1))

# График
ggplot(data = newdata_ht_1, mapping = aes(x)) +
  geom_ribbon(aes(ymin=lwr, ymax=upr), fill='lightskyblue', alpha=1/2) +
  geom_line(aes(y = fit), colour='blue') +
  geom_point(data = assumptions, aes(x, y1)) +
  labs(title='Интервальный прогноз модели y1(x) (95% интервал)',
       y = 'y1')

# Создаем точки для получения нового прогноза
newdata_ht_2 <- data.frame(x = seq_range(assumptions$x, by = 0.1))

# Интервальный прогноз:
pred_ht_transformed_2 <- predict(m2_2, newdata = newdata_ht_2, 
                   interval = 'prediction',
                   level = .95)

# Обратное преобразование:
pred_ht_2 <- pred_ht_transformed_2^(-7)

# Соединяем все в одну таблицу:
newdata_ht_2 <- bind_cols(newdata_ht_2, as.data.frame(pred_ht_2))

# График
ggplot(data = newdata_ht_2, mapping = aes(x)) +
  geom_ribbon(aes(ymin=lwr, ymax=upr), fill='lightskyblue', alpha=1/2) +
  geom_line(aes(y = fit), colour='blue') +
  geom_point(data = assumptions, aes(x, y2)) +
  labs(title='Интервальный прогноз модели y2(x) (95% интервал)',
       y = 'y2')

Границы доверительного интервала расходятся, т.к. к ним применили нелинейное преобразование. Это означает, что уровень неопределенности возрастает при увеличении y. Этого следовало ожидать, т.к. гетероскедастичность - по определению рост дисперсии.

Визуально наиболее достоверным прогнозом кажется интервальный прогноз модели y2(x) (95% интервал).