Многофакторная линейная регрессионная модель зависимости продаж мороженного от погоды
Содержание
- Аннотация
- Подготовка данных
- Разведочный анализ
- Построение модели 4.1. ТОП лето 4.2 ТОП весна
- Выводы
1. Аннотация
Напомню вопросы, которые поставлены перед исследованием:
* Какова степень влияния погодных условий (температура, облачность, влажность, давление, осадки и т.п.)?
* Как повышение средней температуры месяца на 1 градус влияет на изменение продаж?
* Что важнее: температура воздуха или ясная погода?
* Как влияют на продажи другие метеорологические показатели - влажность воздуха, атмосферное давление, скорость и направление ветра?
В отчете исследована зависимость продаж мороженного (т.н. “импульсного” в упаковках до 200гр.) от метеорологических факторов - температуры, влажности, давления, скорости ветра. В целом можно сказать, что погода значительно влияет на продажи, при этом летом это влияние не так заметно, гораздо заметнее связь в межсезонье, зимой связь совсем слабая.(см. таблицу ниже) Также хочется отметить, что основным фактором с подавляющим перевесом является температура воздуха, остальные факторы влияют гораздо менее значимо, например летом дополнительное влияние оказывают влажность и давление, а весной или осенью - скорость ветра и облачность. Однако, усложнение модели с введением этих факторов не оправдано по причине их низкой значимости.
Таблица с параметрами моделей по сезонам.
| Бренд | Сезон | Adj R-squared | Основной предиктор | Дополнительные предикторы |
|---|---|---|---|---|
| ТОП | Лето | 0.43 | tem | humifity, ppp |
| ТОП | Весна | 0.66 | tem | winds, clowdly |
| ТОП | Осень | 0.64 | tem | winds, clowdly |
| ТОП | Зима | 0.13 | tem | - |
2. Подготовка данных
Подключаем нужные библиотеки
library("dplyr")
library("ggplot2")
library("lattice")
library("xtable")
library("car")Загрузка обработанных данных
retail <- dget(file = "./retailtidy")
weather <- dget(file = "./weathertidy")Объединение данных о продажах и погоде
retweat <- merge(x=retail, y=weather, by = "date", all.x = T, all.y = F)Выборка данных для бренда ТОП
retweattop <- retweat[retweat$BRAND=="ТОП",]в дальнейшем работаем с датасетом retweattop
3. Разведочный анализ
предварительный анализ показал, что характер зависимости продаж от погоды значительно различается для различных времен года. смотрите для примера характер зависимости продаж от температуры воздуха (обратите внимание на различие шкал). Поэтому было принято решение строить отдельные модели для различных времен года для каждого из брендов (Топ, Юкки и Soletto)
xyplot(saleskg~tem|season, data = retweattop, layout = c(1, 4), scales
= "free", breaks = 50)
Далее, посмотрим, как взаимосвязаны различные параметры в датасете.
3.1. Диаграмма рассеивания Влажность ~ Температура, цветом обозначены уровни облачности
ggplot(data = retweat, aes(x = tem, y=humidity, col=clowdlyf)) + geom_point()+
xlab("Температура воздуха (гр.Цельсия)")+ylab("Влажность (%)")
3.2. Боксплот Температура ~ Направление ветра
ggplot(data = retweat, aes(x = windd, y=tem)) + geom_boxplot()+
xlab("Направление ветра")+ylab("Температура (гр.Цельсия)")
3.3. Диаграмма рассеивания Давление ~ Температура, цветом обозначены направления ветра
ggplot(data = retweat, aes(x = tem, y=ppp, col=windd)) + geom_point()+
xlab("Температура воздуха (гр.Цельсия)")+ylab("Давление (гПа)")
3.4. Диаграмма рассеивания Давление ~ Влажность, цветом обозначены уровни температуры
ggplot(data = retweat, aes(x = humidity, y=ppp, col=tem)) + geom_point()+
xlab("Влажность воздуха (%)")+ylab("Давление (гПа)")
как показывает разведочный анализ, нет сильно-выраженных зависимостей между метеорологическими показателями, поэтому все их будем включать в регрессионный анализ.
Вывод: на основании предварительного анализа принято решение строить модели отдельно для каждого сезона используя все возможные предикторы.
4.1 Построение модели - бренд: ТОП, сезон: ЛЕТО
tmpseasonstr <- "summer"
tmpseason <- retweattop[retweattop$season==tmpseasonstr,]Подбор модели
Используем функцию step для автоподбора (берем только числовые переменные)
model_full <- lm(saleskg ~ I(tem^2)+tem*clowdly*winds*humidity*ppp,
data = tmpseason[,c("saleskg", "tem", "clowdly", "winds", "humidity", "ppp")])
model_null <- lm(saleskg ~ 1, data = tmpseason[,c("saleskg", "tem", "clowdly", "winds", "humidity", "ppp")])
ideal_model <- step(object = model_full, scope = list(lower = model_null, upper = model_full), direction = 'backward')Результат для наилучшей модели ideal-model
## Residual standard error: 223.7 on 261 degrees of freedom
## Multiple R-squared: 0.4083 Adjusted R-squared: 0.3766
## F-statistic: 12.87 on 14 and 261 DF, p-value: < 2.22e-16| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -281948.4356 | 71120.7431 | -3.96 | 0.0001 |
| I(tem^2) | 3.2647 | 1.1122 | 2.94 | 0.0036 |
| tem | 15010.7442 | 3343.0715 | 4.49 | 0.0000 |
| clowdly | -73.3839 | 33.9305 | -2.16 | 0.0315 |
| winds | 64291.2917 | 25154.7789 | 2.56 | 0.0112 |
| humidity | 3337.7939 | 1030.4499 | 3.24 | 0.0014 |
| ppp | 280.9210 | 70.3252 | 3.99 | 0.0001 |
| tem:winds | -3274.7395 | 1186.1410 | -2.76 | 0.0062 |
| tem:humidity | -188.6255 | 52.5078 | -3.59 | 0.0004 |
| clowdly:humidity | 1.3585 | 0.6952 | 1.95 | 0.0517 |
| tem:ppp | -14.9668 | 3.3073 | -4.53 | 0.0000 |
| winds:ppp | -63.5383 | 24.8338 | -2.56 | 0.0111 |
| humidity:ppp | -3.3237 | 1.0190 | -3.26 | 0.0013 |
| tem:winds:ppp | 3.2341 | 1.1703 | 2.76 | 0.0061 |
| tem:humidity:ppp | 0.1873 | 0.0519 | 3.61 | 0.0004 |
Мы видим значимые различия найденной модели ideal_model по сравнению с моделью с максимальным количеством факторов model_full
xt <- anova(model_full, ideal_model)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) | |
|---|---|---|---|---|---|---|
| 1 | 243 | 12606436.21 | ||||
| 2 | 261 | 13061934.90 | -18 | -455498.70 | 0.49 | 0.9619 |
Но ручной подбор позволяет упростить модель, сохранив при этом значение коэффициента детерминации Adjusted R-squared и значимость коэффициентов
fit <- lm(saleskg^(1/3) ~ I(tem^2)+tem*humidity*ppp, data = tmpseason) ## Residual standard error: 1.068 on 267 degrees of freedom
## Multiple R-squared: 0.3855 Adjusted R-squared: 0.3671
## F-statistic: 20.94 on 8 and 267 DF, p-value: < 2.22e-16| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -705.0332 | 281.9774 | -2.50 | 0.0130 |
| I(tem^2) | 0.0102 | 0.0052 | 1.96 | 0.0506 |
| tem | 36.5901 | 12.9073 | 2.83 | 0.0049 |
| humidity | 12.3805 | 4.7400 | 2.61 | 0.0095 |
| ppp | 0.7064 | 0.2786 | 2.54 | 0.0118 |
| tem:humidity | -0.6822 | 0.2339 | -2.92 | 0.0038 |
| tem:ppp | -0.0365 | 0.0128 | -2.86 | 0.0046 |
| humidity:ppp | -0.0123 | 0.0047 | -2.62 | 0.0092 |
| tem:humidity:ppp | 0.0007 | 0.0002 | 2.93 | 0.0037 |
Проверка модели (подобранной вручную)
1.1 нормальность остатков (тест Шапиро-Вилка)
test1 <- shapiro.test(fit$residuals)cat(ifelse(test1$p.value < 0.05, "Остатки не распределены нормально", "Остатки распределены нормально"), ", p-значение: ", test1$p.value)## Остатки не распределены нормально , p-значение: 0.047194211.2 визуальная оценка нормальности остатков (QQ-plot)
qqPlot(fit, labels=row.names(tmpseason), id.method="identify",
simulate=TRUE, main="Q-Q Plot")
- независимость остатков (тест Дарбина-Уотсона)
test2 <- durbinWatsonTest(fit)cat(ifelse(test2$p < 0.05, "Есть зависимость", "Нет зависимости"), ", p-значение: ", test2$p)## Есть зависимость , p-значение: 0.0323.1 гомоскедастичность остатков
test3 <- ncvTest(fit) cat(ifelse(test3$p < 0.05, "Гомоскедастичность НЕ выполняется", "Гомоскедастичность выполняется"), ", p-значение: ", test3$p)## Гомоскедастичность выполняется , p-значение: 0.17892093.2 Визуальная оценка гомоскедастичности
spreadLevelPlot(fit)
Вывод: можно сказать что модель соответствует требованиям, те небольшие отклонения от нормальности можно нивелировать при подгонке модели.
Подгонка модели
Схема этапа подгонки - выбросы находим с помощью теста Бонферони, и расчета метрики Кука для поиска влиятельных значений, реже визуально на графиках, (как правило удалял немного - около 4-7 наблюдений.) Также на этапе подгонки, в случае невыполнения требования гомоскедастичности или нормальности, рассчитываются показатели степени для преобразования зависимой переменной (используются функции spreadLevelPlot и powerTransform соответственно)
Итак, ищем выбросы с помощью теста Бонферони
tt <- outlierTest(fit)
outn <- which(row.names(tmpseason)==names(tt$p), arr.ind = T)
tmpseason <- tmpseason[-c(outn),]Пересчитаем параметры модели после удаления выбросов и увидим, что все требования выполняются
fit <- lm(saleskg^(1/3) ~ I(tem^2)+tem*humidity*ppp, data = tmpseason) | Тест | Результат исходный | Результат после удаления выбросов |
|---|---|---|
| Нормальность остатков | 0.047 (NO) |
0.294 (YES) |
| Независимость остатков | 0.032 (NO) |
0.07 (YES) |
| Гомоскедастичность остатков | 0.179 (YES) |
0.142 (YES) |
Визуализация результатов
Кривая зависимости продаж от температуры при различных комбинациях давления и влажности
Вывод: На данном графике хочется отметить (см. верхний), что при повышенном давлении продажи выше, если влажность выше. В целом, попытка посмотреть на найденную зависимость при различных срезах предикторов наводит на мысль, что переменные Влажность и Давление не так явно влияют на продажи, но при этом значительно усложняют модель и ее интерпретацию.
Далее, рассмотрим еще модель для бренда ТОП для весны, будем смотреть зависимость только от температуры.
4.2 Построение модели - бренд: ТОП, сезон: ВЕСНА
tmpseasonstr <- "spring"
tmpseason <- retweattop[retweattop$season==tmpseasonstr,]Подбор модели
Здесь подбирал коэффициенты вручную
fit <- lm(saleskg^(1/3) ~ tem+I(tem^2), data = tmpseason) ## Residual standard error: 0.9642 on 273 degrees of freedom
## Multiple R-squared: 0.6638 Adjusted R-squared: 0.6613
## F-statistic: 269.52 on 2 and 273 DF, p-value: < 2.22e-16| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 4.3871 | 0.0971 | 45.19 | 0.0000 |
| tem | 0.0971 | 0.0168 | 5.77 | 0.0000 |
| I(tem^2) | 0.0035 | 0.0008 | 4.49 | 0.0000 |
как видим модель гораздо более полно описывает продажи, чем предыдущая, доля объясненной дисперсии в весенней модели гораздо выше - 0.66 против 0.37 летом. Но это особенность сезона - по всем брендам летом влияние погоды (а точнее я смотрел температуру) не так значительно как весной или осенью. Зимой метео условия вообще мало влияют на продажи.
Проверка модели
1.1 нормальность остатков (тест Шапиро-Вилка)
test1 <- shapiro.test(fit$residuals)## Остатки не распределены нормально , p-значение: 6.171487e-051.2 визуальная оценка нормальности остатков (QQ-plot) 
- независимость остатков (тест Дарбина-Уотсона)
test2 <- durbinWatsonTest(fit)## Есть зависимость , p-значение: 0.0023.1 гомоскедастичность остатков
test3 <- ncvTest(fit) ## Гомоскедастичность НЕ выполняется , p-значение: 4.80107e-083.2 Визуальная оценка гомоскедастичности
spreadLevelPlot(fit)
Вывод: есть значительные отклонения от требований.
Подгонка модели
Ищем выбросы с помощью теста Бонферони и удаляем их
for (i in 1:6) {
fit <- lm(saleskg^(1/3) ~ tem+I(tem^2), data = tmpseason)
tt <- outlierTest(fit)
outn <- which(row.names(tmpseason)==names(tt$p), arr.ind = T)
tmpseason <- tmpseason[-c(outn),]
}Пересчитаем параметры модели после удаления выбросов
fit <- lm(saleskg^(1/3) ~ tem+I(tem^2), data = tmpseason) | Тест | Результат исходный | Результат после удаления выбросов |
|---|---|---|
| Нормальность остатков | 10^{-4} (NO) |
0.0048 (NO) |
| Независимость остатков | 0.002 (NO) |
0.006 (NO) |
| Гомоскедастичность остатков | 0 (NO) |
0 (NO) |
подгонка модели не помогла полностью привести модель в соответствие с требованиями.
Попробуем трансформировать зависимую переменную
tt <- summary(powerTransform(tmpseason$saleskg))
round(tt$result[1],2)## [1] -0.05Как видим функция предлагает нам достаточно радикальное преобразование зависимой переменной z=y^(-0.05). На мой взгляд такое преобразование затруднит интерпретацию модели. Вопрос о выполнении требований к линейной модели пока остается открытым и требует дополнительного исследования.
Визуализация результатов
Кривая зависимости продаж от температуры
ggplot(data = tmpseason, aes(x=tem, y=saleskg)) + geom_point(col="grey", size=4, alpha = 1/2) +
geom_line(aes(x = tem, y = (fit$fitted.values)^3), col = 'blue', lwd=1)+
geom_line(aes(x = tem, y = (fit$fitted.values+1.645*sd(fit$residuals))^3), col = 'red', linetype = 3)+
geom_line(aes(x = tem, y = (fit$fitted.values-1.645*sd(fit$residuals))^3), col = 'red', linetype = 3)+xlab("Температура воздуха (градусов Цельсия)")+ylab("Продажи суточные (кг)")
Выводы
- Какова степень влияния погодных условий (температура, облачность, влажность, давление, осадки и т.п.)?
В целом можно сказать, что погода значительно влияет на продажи, Также хочется отметить, что основным фактором с подавляющим перевесом является температура воздуха, остальные факторы влияют гораздо менее значимо, но усложняют модель. При этом летом это влияние не так заметно, гораздо заметнее связь в межсезонье, зимой связь совсем слабая. >* Как повышение средней температуры месяца на 1 градус влияет на изменение продаж?
Летом прирост температуры на 1 градус влечет рост продаж на 15-18% (расчет выполнен по летней модели при фиксированных Влажности = 70%, давлении = 1000 гПа) Весной рост температуры на 1 градус влечет рост продаж на 7-10%. При этом рост температуры с 0 до 1 градуса дает прирост продаж около 6-7%, а рост температуры с 15 до 16 градусов дает прирост продаж около 9-10%.
- Что важнее: температура воздуха или ясная погода?
Гораздо важнее температура воздуха, если говорить о лете и межсезонье, зимой же ни температура ни ясность не влияют значительно на продажи.
- Как влияют на продажи другие метеорологические показатели - влажность воздуха, атмосферное давление, скорость и направление ветра?
основным фактором с подавляющим перевесом является температура воздуха, остальные факторы влияют гораздо менее значимо, например летом дополнительное влияние оказывают влажность и давление, а весной или осенью - скорость ветра и облачность.