Идея, вопрос

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

Описание, редактирование данных

Начнем с загрузки данных, очистки и описательных характеристик переменных. Для изучения были использованы данные 29-ой волны (репрезентативная выборка, 2020г.) скачанные с сайта rlms.

Загрузим и рассмотрим для поиска регрессии следующие переменные:

  • Вспомните, сколько примерно денег все члены Вашей семьи израсходовали на питание дома и вне дома в течение последних 30 дней? (ye4) Зависимая переменная

  • Каким был денежный доход всей Вашей семьи в течение последних 30 дней? Включите сюда все денежные поступления: заработную плату, пенсии, стипендии, любые другие денежные поступления, в т.ч. и в валюте. (income)

  • Какова общая полезная площадь жилья у Вашей семьи, то есть сумма площадей жилых комнат, кухни, ванной, туалета, прихожей, кладовых и тому подобного в квартире (доме)? (sq_meters)

  • Количество членов семьи (n_members)

  • Ваша семья получала в течение последних 30 дней какие-либо пенсии, включая доплаты от Пенсионного фонда РФ, от государства? (pens_q)

  • Сколько в рублях пенсии, включая доплаты от Пенсионного фонда РФ, от государства получила Ваша семья в течение последних 30 дней? (pens)

  • Тип населенного пункта (place)

  • Если вспомнить последние 3 месяца, сколько рублей в среднем в месяц должна была платить Ваша семья за квартиру, включая аренду, и коммунальные услуги, за вычетом субсидий и льгот, если они у Вас есть? (living_cost)

Объясняемую переменную я подсчитал простым соотношением: расходы на еду к общим доходам. Посмотрим сначала на все наши количественные переменные вместе.

Почти все переменные, кроме количества получаемой пенсии имеют уровень ответа выше 90%. Уже из этой таблицы можно увидеть что доля расходов на еду имеет пропущенные значения - см. медиану и максимальное. На самом деле там значения равные бесконечности. Рассмотрим этот феномен немного ниже.

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

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

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

Видно что есть несколько точек где отношение трат на еду выше единицы (для удобства я провел горизонтальную линию на значении 1), т.е. домохозяйства тратят больше чем получают дохода. Я бы посмотрел на эти точки подробней.

К тому же вызывает вопрос точка с очень большим доходом, справа на графике. Просто из любопытства посмотрим на семью миллионеров. Семья как семья - 4 человека, с доходом 1,943,500 рублей :) в 37.7 кв.м..

Теперь посмотрим на отношения трат на еду больше единицы в виде таблицы.

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

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

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

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

Надо сказать, что прямые линии в этом графике меня довольно сильно смущают. Во-первых видно некоторую “дискретность” данных, т.е. видимо люди округляют числа при опросе, а во вторых видна четкая линейная зависимость (после логарифмического преобразования).

Еще одни данные, рассматриваемые в данной модели - плата за квартиру. Ниже представлен график показывающий разброс значений для этого параметра для разных типов населенных пунктов. Например видно, что в селе медианная плата за квартиру ниже (1800) чем в городе (4000) и областном центре (4530).

Можно провести t-test для проверки значимости разницы среднего значения между выборками (предварительно придется “уравнять” выборки по количеству). Такой тест показывает, что разница статистически значима для разницы между селом и городом.

## 
##  Welch Two Sample t-test
## 
## data:  city and village
## t = 18.286, df = 1981.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1695.501 2102.871
## sample estimates:
## mean of x mean of y 
##  4382.823  2483.637


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

Если подсчитать t-тест между средним доходом в Селе и Областном центре, то это разница статистически значима, p-value сильно меньше 0.05 (при выборе \(\alpha=5\%\))

## 
##  Welch Two Sample t-test
## 
## data:  city and village
## t = 7.388, df = 1999.3, p-value = 2.178e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  10284.5 17717.7
## sample estimates:
## mean of x mean of y 
##   61210.6   47209.5

В то время как разница между средними значениями дохода в Городе и Селе статистически не значима, нулевая гипотеза (что средний доход равен в городе и селе) не отвергается, при уровне значимости 5%.

## 
##  Welch Two Sample t-test
## 
## data:  city and village
## t = 1.2582, df = 1986.6, p-value = 0.2085
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1168.243  5350.670
## sample estimates:
## mean of x mean of y 
##  49300.71  47209.50

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

Начнем с построения простой модели отношение еды/доход к доходу и сравним её с моделью где доход и предсказываемая переменная преобразованы логарифмически. Слева обычная модель, справа модель с лог преобразованием.

  food_ratio log(food_ratio)
Predictors Estimates CI p Estimates CI p
(Intercept) 0.36700518 0.36109876 – 0.37291161 <0.001 2.47987740 2.27570051 – 2.68405428 <0.001
income -0.00000101 -0.00000109 – -0.00000093 <0.001
income [log] -0.35251704 -0.37161043 – -0.33342366 <0.001
Observations 4209 4209
R2 / R2 adjusted 0.134 / 0.134 0.237 / 0.237


В обоих случаях коэффициент статистически значим, хотя сравнивать модели по значению \(R^2\) мы не можем, т.к. предсказываем разные переменные.

Поскольку график выше тоже подсказывал что логарифмическое преобразование дает лучшую зависимость, возьмем такую модель как базовую. Таким образом, модель с логарифмическим преобразованием показывает что при увеличении дохода на один процент отношение трат на еду уменьшается на 0.35 процента.

Построение более сложной модели

Построим модель со всеми рассматриваемыми переменными. Как видно, коэффициент при вопросе о получении пенсии имеет p-value 0.3, что сильно выше чем \(\alpha=5\%\), таким образом я решил исключить эту переменную из регрессии.Переменная отвечающая за жилую площадь имеет значение p-value 0.014, что будет меньше выбранного 5% барьера для ошибки первого рода.

  log(food_ratio)
Predictors Estimates CI p
(Intercept) 3.9872 3.7395 – 4.2349 <0.001
income [log] -0.5150 -0.5401 – -0.4899 <0.001
n_members 0.0873 0.0763 – 0.0984 <0.001
sq_meters -0.0010 -0.0017 – -0.0002 0.014
pens_q [нет] -0.0145 -0.0424 – 0.0133 0.306
place [Обл. центр] 0.1038 0.0725 – 0.1352 <0.001
place [ПГТ] -0.0326 -0.0880 – 0.0228 0.249
place [Село] -0.1734 -0.2112 – -0.1356 <0.001
internet [1] 0.0695 0.0401 – 0.0988 <0.001
Observations 4209
R2 / R2 adjusted 0.319 / 0.318

Ниже показаны параметры окончательной регрессии.
Окончательная модель
  log(food_ratio)
Predictors Estimates std. Error CI p
(Intercept) 3.622787 0.133922 3.360229 – 3.885345 <0.001
income [log] -0.541168 0.013153 -0.566955 – -0.515382 <0.001
n_members 0.086228 0.005585 0.075278 – 0.097179 <0.001
place [Обл. центр] 0.094571 0.015906 0.063386 – 0.125755 <0.001
place [ПГТ] -0.020009 0.028112 -0.075123 – 0.035105 0.477
place [Село] -0.114384 0.020612 -0.154794 – -0.073974 <0.001
sq_meters -0.001141 0.000389 -0.001903 – -0.000379 0.003
internet [1] 0.057043 0.014716 0.028191 – 0.085895 <0.001
living_cost [log] 0.079480 0.010271 0.059343 – 0.099618 <0.001
Observations 4209
R2 / R2 adjusted 0.328 / 0.327


Эта модель имеет все коэффициенты значимыми кроме коэффициента при типе населенного пункта Поселок Городского Типа. Дело может быть в недостаточности данных по этому параметру, он был представлен самыми маленьким набором наблюдений либо действительно в ПГТ зависимость не отличается от города (т.к. город базовая факторная переменная), что мы и увидим ниже.

С помощью F-статистики мы можем сравнить две модели. Проведенный тест показывает, что более сложная модель действительно точнее предсказывает искомую переменную (долю расходов на еду).Н0 в том что более простая модель верна отвергается, т.к. p-value < 0.05.

Res.DfRSSDfSum of SqFPr(>F)
4.21e+03812            
4.2e+03 715796.981.33.78e-111

Тест Рамсея показывает, что в модели видимо есть пропущенные переменные (\(p-value < \alpha=5\%\)), что наверное, не удивительно, учитывая низкий коэффициент корреляции.

## 
##  RESET test
## 
## data:  fit.log.m2
## RESET = 65.047, df1 = 2, df2 = 4198, p-value < 2.2e-16

Графическое представление модели

Ниже графически представлены коэффициенты с их доверительными интервалами.Поскольку мы использовали лог преобразование, то интерпретация коэффициентов может быть следующая - при увеличении дохода на 1% на 0.54 процента уменьшается отношение расходов на еду к общему доходу, для количества людей в семье - при увеличении в семье на одного человека, доля затрат на еду увеличивается на 9 процентов, противоположным образом действует жилая площадь - увеличение на один квадратный метр жилой площади уменьшает долю расходов на еду на 0.1 процент, наличие высокоскоростного интернета увеличивает долю расходов на 6 процентов.


Как же интерпретировать коэффициенты при категорийных регрессорах? На графике ниже видно, что при типе населенного пункта в областном центре кривая проходит выше, значит отношение расходов на еду там в целом выше чем в городе (базовое значение). Это же мы видим и в коэффициентах при этих переменных: положительный при Областном центре и отрицательный при перемной Село. При этом также стоит отметить, что в городе зависимость визуально очень слабо отличается от зависимости в ПГТ, что мы и видели в значениях значимости коэффициента при ПГТ.

Построение прогноза с предиктивным интервалом

На графике снизу показана окончательная модель с предиктивным интервалом. Я построил две кривых из модели: одну для данных в населенном пункте областной центр (красная), другую для населенного пункта село (фиолетовая), т.к. они имеют наибольшее отличие друг от друга. Предиктивный интервал показан только для областного пункта. При построении кривых изменялся только доход, остальные параметры были зафиксированы (кол-во человек в семье - 2; высокоскоростной интернет - да; кол-во жилой площади - 33; квартплата - 4000)

Предсказание для выдуманного наблюдения

Рассмотрим “медианную” семью, проживающую в городе из трех человек. У такой семьи медианный доход будет 57,000, квартплата 4,400 и медианная жилая площадь 34 кв.м..

Если использовать эти данные как вводные, то получим.

Таким образом наша модель предсказывает, что среднее значение доли доходов, которое данная семья будет тратить на еду будет 26%, предиктивный интервал будет равен [11%; 58%].

Скрипт

## ----setup, include=FALSE-------------------------------------------------------------------
knitr::opts_chunk$set(echo = FALSE,message=FALSE,warning = FALSE)
library(tidyverse) # манипуляции с данными
library(skimr) # описательные статистики
library(tidymodels) # удобства для оценивания моделей
library(sjPlot) # визуализация коэффициентов
library(car) # проверка сложных гипотез
library(lmtest) # тесты для моделей
library(texreg) # таблички для сравнения моделей
library(rio) # import / export всех форматов
library(labelled)
library(sjlabelled)
library(cowplot)
library(reshape2)
library(plotly)
library(knitr)
library(kableExtra)
library(DT)
library(rstatix)
library(ggpubr)
library(huxtable)


## ----input, cache=TRUE----------------------------------------------------------------------
df <- import("Data/r29h_os41.sav")
df.na.flt <- df %>% 
  mutate_if(is.numeric,~ifelse(.>99999991, NA_real_,.)) # убираем флаги отвечающие за "нет ответа"


## ----comments filter data-------------------------------------------------------------------
# yf14 Каким был денежный доход всей Вашей семьи в течение последних 30 дней? Включите сюда все денежные поступления: заработную плату, пенсии, стипендии, любые другие денежные поступления, в т.ч. и в валюте

# ye4 Вспомните, сколько примерно денег все члены Вашей семьи израсходовали на питание дома и вне дома в течение последних 30 дней?

# yf10 В течение последних 30 дней Ваша семья получала зарплату или какие-либо другие выплаты по основному или дополнительному месту работы в виде денег, товаров или услуг?

# yc6 Какова общая полезная площадь жилья у Вашей семьи, то есть сумма площадей жилых комнат, кухни, ванной, туалета, прихожей, кладовых и тому подобного в квартире (доме)?

# y_nfm Количество членов семьи

# yf12.1a Ваша семья получала в течение последних 30 дней какие-либо пенсии, включая доплаты от Пенсионного фонда РФ, от государства?

# yf12.1b Сколько в рублях пенсии, включая доплаты от Пенсионного фонда РФ, от государства  получила Ваша семья в течение последних 30 дней?

# status Тип населенного пункта


df.flt <- df.na.flt %>%   
  mutate(food_ratio=ye4/yf14) %>% #cчитаем отношение денег потраченных на еду к доходу
  select(yf14,ye4,y_nfm,status,yf10,food_ratio, yc5, yf12.1a,yf12.1b,yc9.624a,ye12.2) %>%#отбираем нужные переменные 
  rename(income=yf14,food_exp=ye4,n_members=y_nfm,sq_meters=yc5,pension=yf12.1b,internet=yc9.624a,living_cost=ye12.2) %>% #переименовываем переменные
  mutate(add_income=as.factor(if_else(yf10==1,"нет","да"))) %>% #превращаем категорийноую переменную в красивый ответ, для графиков
  mutate(place=as.factor(case_when(
    status == 1 ~ "Обл. центр",
    status == 2 ~ "Город",
    status == 3 ~ "ПГТ",
    status == 4 ~ "Село"  ))) %>%   
  mutate(internet=as.factor(if_else(internet==1,1,0))) %>% 
  mutate(pens_q=as.factor(if_else(yf12.1a==1,"да","нет"))) %>% 
  select(-c(yf10,status,yf12.1a))


## ----summary table--------------------------------------------------------------------------
df.flt %>% 
skim_without_charts() %>% 
  # names()
  dplyr::select(-c(numeric.sd,numeric.p25,numeric.p75)) %>%
  yank("numeric") %>% 
  mutate_if(is.numeric,round,2) %>% 
  datatable(rownames = FALSE,options=list(pageLength=7), autoHideNavigation = TRUE)
  


## ----summary quantity table-----------------------------------------------------------------
df.flt %>% 
skim() %>% 
  # names()
  select(-factor.ordered) %>% 
  yank("factor") %>% 
  mutate_if(is.numeric,round,3) %>%
  datatable(rownames = FALSE,options=list(pageLength=5), autoHideNavigation = TRUE)


## ----places percentage----------------------------------------------------------------------
df.place_ratio <- df.flt %>% 
  group_by(place) %>% 
  summarise(count=n()) %>% 
  mutate(ratio=round(count/sum(count),2)*100)  
  
fig <- plot_ly(df.place_ratio,y=~place,x=~ratio,type="bar",alpha=0.7,color="#2c7bb6",
               hoverinfo="text",
               text=paste(df.place_ratio$place,"<br>", df.place_ratio$ratio,"%"))
fig <- fig %>% layout(xaxis = list(title="Процент, %"),yaxis = list(title="Тип населенного пункта"))
fig


## ----initial scatter income-food_ratio------------------------------------------------------
gg1 <- ggplot(df.flt,aes(y=food_ratio,x=income,
                         color=cut_interval(sq_meters,4),
                         text=paste0("Доход: ", income,
                                     "<br> Еда/доход: ", round(food_ratio,2),
                                     "<br> Кол-во человек: ", n_members,
                                     "<br> Кол-во кв. метров: ", sq_meters)))+
  geom_point(alpha=0.8)+
  scale_color_brewer(palette = "RdYlBu")+
  geom_hline(yintercept = 1, linetype =2,size=0.4,color='red' )+
  theme_cowplot()+
  background_grid()+
  labs(y="Отношение затрат на еду к доходу", x = "Доход", color="Кол-во кв. м")

ggplotly(gg1,tooltip="text")



## ----food_ratio table-----------------------------------------------------------------------
df.flt %>% 
  select(food_exp,income,food_ratio) %>% 
  arrange(-food_ratio) %>% 
  head(30) %>% 
  datatable(head(iris),rownames = FALSE, editable = list(
  target = 'row', disable = list(columns = c(1, 3, 4))))


## ----food_ratio histogram-------------------------------------------------------------------
df.flt.out <- df.flt %>% 
  filter(!is.infinite(food_ratio), food_ratio<1) 
  
fig <- plot_ly(df.flt.out,x=~food_ratio, type="histogram",
               color=~place,
               histnorm="percent",
               colors = "RdYlBu",
               alpha=0.6)
fig <- fig %>% layout(xaxis = list(range = c(0,1),title="Доля расходов на еду"),barmode="relative",bargroupgap=0.1,
                      yaxis = list(title="Проценты"))
fig


## ----food_ratio table per place-------------------------------------------------------------
df.flt.out %>% 
  select(food_ratio,place) %>% 
  group_by(place) %>% 
  skim_without_charts() %>% 
  select(-c(complete_rate,n_missing)) %>% 
  yank("numeric") %>% 
  mutate_if(is.numeric,round,3) %>%
  datatable(rownames = FALSE,options=list(pageLength=4), autoHideNavigation = TRUE)


## ----scatter log-log------------------------------------------------------------------------

gg1 <- ggplot(df.flt.out,aes(y=log(food_ratio),x=log(income),color=place,
                             size=living_cost,
                             text=paste0("Доход: ", income,
                                     "<br> Еда/доход: ", round(food_ratio,2),
                                     "<br> Квартплата: ", round(living_cost,2),
                                     "<br> Кол-во человек: ", n_members,
                                     "<br> Кол-во кв. метров: ", sq_meters,
                                     "<br> Проживание: ", place)))+
  geom_point(alpha=0.4)+
  theme_cowplot()+
  background_grid()+
  scale_color_brewer(palette = "RdYlBu")+
  # coord_cartesian(ylim=c(0,1))+
  # geom_smooth(method="lm")+
  labs(color=paste("Населенный","<br>","пункт"),size="",x="Доход", y="Отношение трат на еду/питание")+
  theme(legend.text=element_text(size=10),legend.title = element_text(size=11),
        plot.margin=unit(c(1,2,1,1),"cm"))
  

ggplotly(gg1,tooltip = "text",autosize = FALSE, hovermode = "x unified")


## ----plot living cost-----------------------------------------------------------------------
fig <- plot_ly(df.flt.out,y=~living_cost, boxpoints="all",
               jitter = 0.4,
               alpha=0.5,
               type="box",color=~place,
               pointpos = -2,
               colors = "RdYlBu")
fig <- fig %>% layout(yaxis = list(title="Плата за квартиру, ком. услуги",range = c(0,20000)),
                      xaxis = list(title="Тип населенного пункта"))
fig


## ----t-test town village living cost--------------------------------------------------------
city <- filter(df.flt.out,place == "Город") %>% select(living_cost) %>% sample_n(1015)
# # dim(city)
village <-  filter(df.flt.out,place == "Село")  %>% select(living_cost)
# # dim(village)
t.test(x=city,y=village,alternative = "two.sided",var.equal = FALSE)





## ----plot income place----------------------------------------------------------------------
fig <- plot_ly(df.flt.out,y=~income, boxpoints="all",
               jitter = 0.4,
               alpha=0.5,
               type="box",color=~place,
               pointpos = -2,
               colors = "RdYlBu")
fig <- fig %>% layout(yaxis = list(title="Доход",range = c(0,3E+5)),
                      xaxis = list(title="Тип населенного пункта"))
fig

## -------------------------------------------------------------------------------------------
df.flt.out %>% 
  select(place,income) %>% 
  group_by(place) %>% 
  skim_without_charts() %>% 
  select(-c(complete_rate,n_missing)) %>% 
  yank("numeric") %>% 
  mutate_if(is.numeric,round,2) %>%
  datatable(rownames = FALSE,options=list(pageLength=4), autoHideNavigation = TRUE)
  # summarize(среднее =  mean(income), медиана =median(income), min=min(income),max=max(income) )


## ----t-test city center village-------------------------------------------------------------
city <- filter(df.flt.out,place == "Обл. центр") %>% select(income) %>% sample_n(1015)
# # dim(city)
village <-  filter(df.flt.out,place == "Село")  %>% select(income)
 # dim(village)
t.test(x=city,y=village,alternative = "two.sided",var.equal = FALSE)



## ----t-test town village income-------------------------------------------------------------
city <- filter(df.flt.out,place == "Город") %>% select(income) %>% sample_n(1015)
# # dim(city)
village <-  filter(df.flt.out,place == "Село")  %>% select(income)
 # dim(village)
t.test(x=city,y=village,alternative = "two.sided",var.equal = FALSE)



## ----simple regression----------------------------------------------------------------------
df.flt.regr <- df.flt.out %>% 
  select(-c(pension,add_income)) %>% 
  na.omit()

fit <- lm(food_ratio ~ income,df.flt.regr)
fit.log <- lm(log(food_ratio) ~ log(income),df.flt.regr)
tab_model(fit,fit.log,digits = 8)


## ----first pass complex regression-----------------------------------------------------------
fit.log.m <- lm(log(food_ratio) ~ log(income)+n_members+sq_meters+pens_q+place+internet,df.flt.regr)
tab_model(fit.log.m,digits=4)


## ----final regression-----------------------------------------------------------------------
fit.log.m2 <- lm(log(food_ratio) ~ log(income)+n_members+place+sq_meters+internet+log(living_cost),df.flt.regr)
tab_model(fit.log.m2, show.se = TRUE,digits=6, show.fstat = TRUE,   title = "Окончательная модель")


## ----anova----------------------------------------------------------------------------------
 anova(fit.log,fit.log.m2)


## ----resettest------------------------------------------------------------------------------
resettest(fit.log.m2)


## ----coeff plot-----------------------------------------------------------------------------
 plot_model(fit.log.m2,show.values = TRUE,vline.color = "darkred")+
    scale_y_continuous()+
  theme_sjplot2()+
  labs(title="Коэффициенты для построенной модели",  y="Значения коэффициентов с доверительным интервалом (95%)")


## ----curves comparison, fig.width=9---------------------------------------------------------
 plot_model(fit.log.m2,type="eff",terms=c("income","place"))+
    scale_y_continuous(limits = c(0,0.5))+
  scale_x_continuous(limits=c(0,250000))+
  theme_sjplot2()+
  labs(title="Сравнение моделей по типу населенного пункта",  x="Доход", y="Отношение расходов на еду к общему доходу", color="Тип населенного пункта")


## ----simulated data-------------------------------------------------------------------------
newx <- tibble(
  income=seq(1500,2e+6,length.out=500),
  n_members=rep(2,500),
  place = as.factor(rep("Обл. центр",500)),
  internet= as.factor(rep(1,500)),
  sq_meters = rep(33,500),
  living_cost = rep(4000,500)
)
# skim(df.flt.out)
# skim(newx)
newy <- predict(fit.log.m2, newdata = newx,interval="prediction") %>%  exp() %>% data.frame()
newxy <- bind_cols(newx,newy)

newyconf <- predict(fit.log.m2, newdata = newx,interval="confidence") %>%  exp() %>% data.frame()
colnames(newyconf) <- c("fit_conf","lwr_conf","upr_conf")
newxyboth <- bind_cols(newxy,newyconf)
# View(newxyboth)


newv <- tibble(
  income=seq(1500,2e+6,length.out=500),
  n_members=rep(2,500),
  place = as.factor(rep("Село",500)),
  internet= as.factor(rep(1,500)),
  sq_meters = rep(33,500),
  living_cost = rep(4000,500)
)
newyv <- predict(fit.log.m2, newdata = newv,interval="prediction") %>%  exp() %>% data.frame()
newxyv <- bind_cols(newv,newyv)



## ----final plot with simulated data, fig.width=9--------------------------------------------
fig <- plot_ly(data=df.flt.out,x=~income,y=~food_ratio,
               type="scatter",
               mode="markers",
               alpha=0.1,
               size=~sq_meters,
               sizes = c(5,100),
               hoverinfo = "text",
               name="Исходные данные",
               color = "#fdae61",
               opacity=0.7,
               text=paste("Доход: ", round(df.flt.out$income,2),
                          "<br>Доля: ", round(df.flt.out$food_ratio,2),
                          "<br>Квартплата: ", round(df.flt.out$living_cost,2),
                          "<br>Место: ", df.flt.out$place)) %>% 

    add_trace(data=newxy,x=~income,y=~fit,type="scatter",mode="line",
              hoverinfo="text",
              alpha=1,
              name = "(Обл. центр) среднее, предсказание",
              text=paste("(Обл. центр) среднее: ", round(newxy$fit,2)),
              line = list(color='#d7191c',width=2, opacity=1)) %>% 
  add_trace(data=newxy,x=~income,y=~lwr,type="scatter",mode="line",
              hoverinfo="text",
              name = paste("(Обл. центр) Нижний предиктивный", "<br>интервал (95%)"),
              text=paste("(Обл. центр) Нижний пред.инт. : ", round(newxy$lwr,2)),
              line = list(color='#2c7bb6',width=2, dash="dot")) %>% 
  add_trace(data=newxy,x=~income,y=~upr,type="scatter",mode="line",
              hoverinfo="text",
            name = paste("(Обл. центр) верхний предиктивный", "<br>интервал (95%)"),
              text=paste("(Обл. центр) Верхний пред.инт. : ", round(newxy$upr,2)),
              line = list(color='#2c7bb6',width=2, dash="dot")) %>% 
  add_trace(data=newxyboth,x=~income,y=~upr_conf,type="scatter",mode="line",
              hoverinfo="text",
            name = paste("(Обл. центр) верхний доверительный", "<br>интервал (95%)"),
              text=paste("(Обл. центр) Верхний дов.инт. : ", round(newxyboth$upr_conf,2)),
              line = list(color='#4f4f4f',width=2, dash="dash")) %>% 
  add_trace(data=newxyboth,x=~income,y=~lwr_conf,type="scatter",mode="line",
              hoverinfo="text",
            name = paste("(Обл. центр) нижний доверительный", "<br>интервал (95%)"),
              text=paste("(Обл. центр) Нижний дов.инт. : ", round(newxyboth$lwr_conf,2)),
              line = list(color='#4f4f4f',width=2, dash="dash")) %>% 
   add_trace(data=newxyv,x=~income,y=~fit,type="scatter",mode="line",
              hoverinfo="text",
              alpha=1,
              name = "(село) среднее, предсказание",
              text=paste("(село) среднее: ", round(newxyv$fit,2)),
              line = list(color='#3300cc',width=2, opacity=1))
    

fig <- fig %>% layout(xaxis = list(title = "Доход",range = c(0,250000)),
                      yaxis = list (title = "Доля расходов на еду",range = c(0,1)),
                      hovermode = "x unified")
fig


## ----median overview------------------------------------------------------------------------
df.flt.regr %>% 
  filter(place == "Город") %>% 
  select(c(income,living_cost,n_members, sq_meters)) %>% 
  group_by(n_members) %>% 
  rename(`Кол-во человек в семье`=n_members) %>% 
  summarize(`Доход (медиана)`=median(income),
            `Квартплата (медиана)`=median(living_cost),
            `Жилая площадь (медиана)`=median(sq_meters)) %>% 
  mutate_if(is.numeric,round,0) %>% 
  datatable(rownames = FALSE,options=list(pageLength=11),autoHideNavigation = TRUE)


## ----one sample prediction------------------------------------------------------------------
family <- tibble(
  income = 57000,
  living_cost = 4400,
  n_members = 3,
  sq_meters=34,
  internet = factor(1),
  place = "Город"
)

family_pred <- predict(fit.log.m2,newdata = family,interval="predict") %>% exp()*100 

family_pred %>%
  data.frame() %>% 
mutate_if(is.numeric,round,0) %>% 
  rename(среднее=fit,
                `нижний предиктивный интервал`=lwr,
                `верхний предиктивный интервал`=upr) %>%
  datatable(rownames = FALSE,options=list(pageLength=1), autoHideNavigation = TRUE)