Анализ выборочных совокупностей

# Настройки R и knitr ----------------------------
# опции R
options(width = 70, scipen = 16, digits = 3) 
# опции чанков по умолчанию
knitr::opts_chunk$set(tidy = FALSE, # сохранять форматирование кода
               fig.width = 4.5, # ширина рисунков по умолчанию
               fig.height = 3)  # высота рисунков по умолчанию

{r include=FALSE}
# Загрузим необходимые пакеты
library(ggplot2)
library(car)
library(lmtest)
library(dplyr)
library(MLmetrics)
library(forecast)
library(cowplot)
library(agricolae)
library(sandwich)
# Установим тему для оформления графиков
theme_set(theme_bw(base_size = 12))

Введение

В 2016 году рынок жилой недвижимости в России был на пике перегрева. Цены на все типы жилья в столице, а также индексы стоимости упали относительно 2015 года, продолжая общий нисходящий тренд. Эксперты сравнивали рынок жилой недвижимости с рынком нефти 2013-2014 годов, когда цены показывали отрицательную динамику в сравнении с предыдущим периодом, при этом производители (в случае рынка жилой недвижимости - застройщики) продолжали наращивание предложения [1]. На конец 2022 году аналитики прогнозируют наступление ситуации на рынке жилой недвижимости, схожей с ситуацией в 2016 году, вследствие окончания действия программы низких ипотечных ставок [2]. Вследствие предположения о том, что рынок недвижимости в среднем стабилен, было решено исследовать регрессоры, влияющие на стоимость квартиры. При этом было решено ориентироваться на крупный город России в силу отсутствия возможных данных для всех показателей во всех городах, к примеру, влияние расстояния квартиры от метро может иметь влияние, однако для городов без метрополитена данное влияние будет неактуально. Следовательно, невозможно в рамках текущей осведомленностью базами данных построить репрезентативную модель по выборке для всей страны. Поэтому было решено построить модели, описывающие влияние на цену квартиры в Москве, и подобрать наиболее оптимальную. Также помимо этого было решено сформулировать дополнительны гипотезы, которые помогут наиболее лучше понять ситуацию на рынке жилой недвижимости в Москве.

Постановка задач

1. Определить число и тип наблюдений;

2. Определить, есть ли отсутствующие значения в ячейках;

3. Проанализировать максимальные, медианные, средние и минимальные значения для количественных переменных;

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

5. Определить, существует ли разница в ценах квартир в зависимости от географии расположения;

6. Построить диаграмму рассеяния и оценить коэффициенты корреляции между переменными

7. Построить первоначальную модель, оценить валидность модели. Найти наиболее оптимальную и интерпретируемую модель

Выбор и описание данных

После формулирования гипотез следует второй наиболее трудозатратный этап - поиск данных. Под наши гипотезы соответственно необходимыми параметрами являются цены на жилье в Москве, а также возможные регрессоры, которые описывали бы переменную-отклик. Для работы с гипотезами были выбраны данные с официального портала открытых данных правительства Москвы [3] за 2016 год. Данный датасет содержит информацию о следующих параметрах:

1. n – номер квартиры по порядку

2. price – цена квартиры в $1000

3. totsp – общая площадь квартиры, кв.м.

4. livesp жилая площадь квартиры, кв.м.

5. kitsp – площадь кухни, кв.м.

6. dist – расстояние от центра в км.

7. metrdist – расстояние до метро в минутах

8. walk – 1 – пешком от метро, 0 – на транспорте

9. brick 1 – кирпичный, монолит ж/б, 0 – другой

10. floor 1 – этаж кроме первого и последнего, 0 – иначе.

11. code – число от 1 до 8, при помощи которого мы группируем наблюдения по подвыборкам:

11.1. Наблюдения сгруппированы на севере, вокруг Калужско-Рижской линии метрополитена

11.2. Север, вокруг Серпуховско-Тимирязевской линии метрополитена

11.3. Северо-запад, вокруг Замоскворецкой линии метрополитена

11.4. Северо-запад, вокруг Таганско-Краснопресненской линии метрополитена

11.5. Юго-восток, вокруг Люблинской линии метрополитена

11.6. Юго-восток, вокруг Таганско-Краснопресненской линии метрополитена

11.7. Восток, вокруг Калиниской линии метрополитена

11.8. Восток, вокруг Арбатско-Покровской линии метрополитена

Этап анализа датасета

# Установим директорию
setwd("C:/Users/User/Desktop")
# Присвоим данным переменную "flats".
flats <- read.table('flats_moscow.txt', sep='\t', header=T)
# Удалим 1 столбец в дата фрейме, поскольку он содержит лишь нумерацию данных по порядку.
flats <- flats[,-1]

Задание 1

nrow(flats)
str(flats)

Выборка имеет 2040 наблюдений и 10 переменных, при этом исходя из описания данных, можно заметить, что 4 переменные: walk, brick, floor, code являются факторными. С помощью функции string было выявлено, что исходный тип ранее упомянутых переменных является целочисленное значение (integer). Для дальнейшего исследования корректно трансформировать данные переменные в факторные.

flats_f <- flats
flats_f$walk <- factor(flats_f$walk)
flats_f$brick <- factor(flats_f$brick)
flats_f$floor <- factor(flats_f$floor)
flats_f$code <- factor(flats_f$code)
str(flats_f)

Задание 2

colSums(is.na(flats))

Пропущенных значений в ячейках не наблюдается.

Задание 3

summary(flats[c(1:6)])

Можно заметить, что практически для всех количественных переменных справедливо можно отметить равенство средних и медианных значений, что наталкивает на мысль о нормальности распределения этих параметров, за исключением переменной price. Для переменной price также можно заметить, что 75% квартир находятся примерно в одном ценовом сегменты, находясь между значениями в 50 тыс. и 142 тыс; при этом 25% квартир имеют цену от 142 тыс. до 730 тыс., что предварительно говорит либо о высокоценовом районе расположения квартир с наибольшей стоимостью, либо о влияния определенных параметров на общую цены квартиры в Москве, либо об этих предпосылках в совокупности. Для остальных переменных сильных различий не замечено.

plot_grid(
  ggplot(data=flats, aes(x=price)) + geom_histogram(fill='powderblue', colour='black'),
  ggplot(data=flats, aes(x=totsp)) + geom_histogram(fill='powderblue', colour='black'),
  ggplot(data=flats, aes(x=livesp)) + geom_histogram(fill='powderblue', colour='black'),
  ggplot(data=flats, aes(x=kitsp)) + geom_histogram(fill='powderblue', colour='black'),
  ggplot(data=flats, aes(x=dist)) + geom_histogram(fill='powderblue', colour='black'),
  ggplot(data=flats, aes(x=metrdist)) + geom_histogram(fill='powderblue', colour='black')
)

Рис. 1. Плотности распределения переменных

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

Задание 4

outlier_tukey <- subset(flats_f$price, flats_f$price < quantile(flats_f$price, 0.25) - 1.5*IQR(flats_f$price) |
                          flats_f$price > quantile(flats_f$price, 0.75) + 1.5*IQR(flats_f$price))
tidy_flats <- flats_f[!flats_f$price %in% outlier_tukey, ]

По методу Тьюки было выявлено 98 значений-выбросов для переменной цены, что неудивительно, поскольку при изначальном анализе гистограммы цены на квартиры в Москве был замечен тяжелый правый хвост, свидетельствующий о необычно дорогих квартирах. Следовательно, в новую переменную tidy_flats были отправлены очищенные от выбросов данные.

Задание 5

Median.test(tidy_flats$price, tidy_flats$code, correct=T)

С помощью функции median.test, которая проверяет гипотезу H0: нет значимой разницы в медианных значениях квартир в разных географических областях, было выяснено, что существует значимая разница хотя бы для одного медианного значения, поскольу p.value оказалось меньше 0.05.

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

ggplot(tidy_flats, aes(x = code, y = price)) +
  geom_boxplot() +
  labs(x = 'Код расположения', y = 'Цена квартиры')

Рис. 2. График "ящик с усами" для цены квартиры в зависимости от места расположения

Исходя из графика также заметно значительное отличие медианных значений цены в Москве по местоположению, в особенности выделяются Северо-запад, вокруг Замоскворецкой линии метрополитена, и Северо-запад, вокруг Таганско-Краснопресненской линии метрополитена

Задание 6

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

plot(tidy_flats[1:6])

Рис. 3. Диаграмма рассеивания

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

Оценим коэффициенты корреляции

cor(flats_f[1:6])

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

С помощью графика "ящик с усами" оценим влияние дамми переменных на цену квартиры.

plot_grid(
  ggplot(tidy_flats, aes(x = walk, y = price)) +
    geom_boxplot() + labs(x = 'Пешая доступность от метро', y = 'Цена квартиры') +
    scale_x_discrete(labels = c('Нет', 'Да')),
  ggplot(tidy_flats, aes(x = brick, y = price)) +
    geom_boxplot() + labs(x = 'Материал стен', y = 'Цена квартиры') +
    scale_x_discrete(labels = c('Другой', 'Кирпичный монолит')),
  ggplot(tidy_flats, aes(x = floor, y = price)) +
    geom_boxplot() + labs(x = 'Этаж', y = 'Цена квартиры') +
    scale_x_discrete(labels = c('Иначе', 'Кроме первого и последнего'))
)

Рис. 4. Графики "ящик с усами" для дамми переменных

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

Задание 7

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

model1 <- lm(data=tidy_flats, price~.-dist - metrdist - walk - code)
summary(model1)
MAPE(model1$fitted.values, tidy_flats$price)

Модель (model1) оказалось значимой, при этом коэфиициенты в модели также оказались значимыми на уровне 5%, за исключением переменной жилой площади. Доля объясненных скорректированных значений оказалась равна 55.6%, а средняя ошибка апроксимации - 12,7%. В принципе предварительная модель оказалась довольно хорошо предсказывающей. Однако, как было отмечено ранее, возможно стоить включить дамми переменную пешой доступности от метро.

model2 <- lm(data=tidy_flats, price~.-dist -metrdist - code)
summary(model2)
MAPE(model2$fitted.values, tidy_flats$price)

Новая модель (model2) также осталась значимой, как и ее коэффиценты на уровне 5%, однако коэффициент перед переменной жилой площади оказался значимым лишь на уровне 1%. Тем не менее, увеличилась доля скорректированных объясненных значений до 57.9%, а средняя ошибка апроксимации упала до уровня 12,4%. В целом можно опираться от второй модели в нахождении наиболее валидной, поскольку она оказалась более лучше предсказывающей значения.

Далее оценим наличие мультиколлинеарности между регрессорами во 2 модели с помощью коэффициента раздутости дисперсии.

vif(model2)

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

Оценим коэффиент раздутости дисперсии третьей модели (model3), удалив из второй модели переменную общей площади квартиры.

model3 <- lm(data=tidy_flats, price~.-dist -metrdist - code - totsp)
vif(model3)

Показатели для текущих переменных оказались меньше 2, что говорит о допостимости новой модели вследствие наличия лишь слабой мультиколлинеарности. Возможно, следовало ввести дополнительный штраф для модели, построив РИДЖ или ЛАССО регрессию, однако для этого необходим детальный анализ коэффициента лямбда на основе кросс-валидации.

Оценим новую третью модель.

vif(model3)
summary(model3)
MAPE(model3$fitted.values, tidy_flats$price)

Модель оказалась все еще значимой, при этом все коэффиценты также значимы на уровне 5%. Доля скорректированных объясненных значений снизиласть незначительно до 53.8%, а средняя ошибка апроксимации осталась на допустимом уровне в 13%.

Теперь проанализируем остатки для оценки применимости модели.

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

model3_diag <- data.frame(fortify(model3), tidy_flats)
D_cut_off <- 4/(nrow(model3_diag)-2)

ggplot(model3_diag, aes(x = 1:nrow(model3_diag), y = .cooksd)) +
  geom_bar(stat = 'identity') + coord_cartesian(ylim = c(0, 0.05)) +
  geom_hline(yintercept = D_cut_off, linetype = 2) +
  labs(x='Номера наблюдений', y='Расстояние Кука')

Рис.5. График расстояния Кука

Можно отметить, что влиятельных наблюдений не выявлено, поскольку не оказалось расстояний Кука, превышаюших стандартный уровень - 1. Однако относительно уровня, учитавающего число наблюдений в модели, который оказался равен 0.002, присутствуют влиятельные наблюдения. Тем не менее, их общее влияние на модель можно считать незначимым.

Далее проверим разброс дисперсии остатков.

ggplot(data = model3_diag, aes(x = .fitted, y = .stdresid)) +
  geom_point() + geom_hline(yintercept = 0) + geom_smooth() + 
  labs(x='Предсказанные значения', y='Стандартизованные остатки')

Рис. 6. График остатков от предсказанных значений

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

Однако также проведем тесты Голдфелда-Квандта и Бройша-Пагана для оценки наличия гетероскедастичности

gqtest(model3, alternative='t')
bptest(model3, studentize=F)

Исходя из результатов, выявлено, что тест Голдфелда-Квандта выдал значение p.value > 0.05, то есть нет оснований отвергать H0 о наличии гомоскедастичности. Однако тест Бройша-Пагана выдал значение p.value < 0.05, вследствие чего H0 о наличии гомоскедастичности должна быть отвергнута. Основываясь на выводах из графика и тестах, можно сделать замечание, что остатки обладают гетероскедастичностью. Вследствие чего следует оценить коэффициенты с помощью робастных ошибок, устойчивых к гетероскедастичности.

coeftest(model3, vcov.=vcovHC(model3))

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

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

coeftest(model3, vcov.=vcovHC(model3))
conftable <- coeftest(model3, vcov.=vcovHC(model3))
ci <- data.frame(estimate=conftable[,1],
                 se_HC=conftable[,2])
ci <- mutate(ci, left_ci=estimate-1.96*se_HC,
             right_ci=estimate+1.96*se_HC)
ci
confint(model3)

ci
confint(model3)

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

Далее проверим остатки на нормальность.

qqPlot(model3, xlab='Квантили', ylab='Остатки модели')

Рис. 7. График Q-Q plot

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

Проведем тест Шапиро-Уилка для проверки гипотезы H0 о нормальном распределение остатков.

shapiro.test(model3$residuals)

По рузультатам теста заметно, что p.value оказалось меньше 5%, что действительно говорит о ненормальном распределении остатков. Нарушение данной предпоссылки довольно существенно и предположительно вызвано либо нелинейностью зависимости переменной отклика от регрессорова, либо наличием гетероскедастичности остатков.

Далее следует проверить предпосылку о равенстве математического ожидания остатков 0, однако в силу ненормальности их распределения это будет не корректно. Тем не менее, с помощью t.test проверим гипотезу H0 о равенстве математического ожидания остатков 0.

t.test(model3$residuals)

Исходя из результатов теста, можно уведить, что p.value оказалось более 5%, вследствие чего нет оснований отвергать гипотезу H0.

Последним шагом проверим предпосылку о наличии автокорреляции между остатками с помощью теста Дарбина-Уотсона.

dwtest(model3)

По результатам теста p.value оказалось более 5%, следовательно, нет оснований отвергнать H0 об отсутствии автокорреляции между остатками.

Вывод

В итоге мы получили линейную модель следующего вида, построенную на основе МНК:

price = -19.92 + 1.83livesp + 4.67kitsp + 9.2walk1 + 13.09brick1 + 6.79*floor1

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

Можно выделить наиболее влиятельные коэффициенты для визуализации упрощенной линейной модели. Для это стандартизуем количественные переменные в модели.

lm(data=tidy_flats, price~scale(livesp) + scale(kitsp) + walk + brick + floor)

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

ggplot(tidy_flats, aes(x = livesp, y = price, colour=brick)) +
  geom_point() + labs(x = 'Жилая площадь', y = 'Цена квартиры') + geom_smooth(method = 'lm')

Рис. 8. Упрощенная линейная модель

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

Ссылки на использованные источники

[1] Индикаторы рынка недвижимости // Обзор рынка недвижимости по итогам 2016 года [Электронный ресурс]. URL: https://www.irn.ru/news/112490.html(дата обращения:19.12.2021), режим доступа - свободный.

[2] Журнал стратегия // Клуб экспертов: прогнозы по рынку недвижимости на 2021-2022 годы [Электронный ресурс]. URL: https://strategyjournal.ru/ekonomika-i-biznes/klub-ekspertov-prognozy-po-rynku-nedvizhimosti-na-2021-2022-gody/(дата обращения:19.12.2021), режим доступа - свободный.

[3] Портал открытых данных правительства Москвы [Электронный ресурс]. URL: https://data.mos.ru/(дата обращения: 10.06.2021), режим доступа - свободный.