Павел Сивохин, “Введение в пространственное моделирование. Тест Морана”
Подготовка
Загружаем пакеты.
#-------------------------Package Installer--------------------------
# load packages and install if missing
# thanks to Richard Schwinn for the code, http://stackoverflow.com/a/33876492
# list the packages you need
p <- c('data.table', 'dplyr', 'sf', 'tmap', 'mapview', 'spdep', 'osmdata')
# this is a package loading function
loadpacks <- function(package.list = p) {
new.packages <- setdiff(package.list, installed.packages()[,'Package'])
if(length(new.packages) != 0) {
install.packages(new.packages, repos = "https://cran.rstudio.com")
}
lapply(package.list, require, character.only = TRUE)
}
loadpacks(p) # calling function to load and/or install packages
rm(loadpacks, p) # cleanup namespace
#----------------------End of Package Installer----------------------Загружаем файлы
df <- fread("EKB_Cian_avito_flats_rent_2021_01_01_2022_12_22.csv", encoding = 'UTF-8')Посмотрим их структуру
df[1, ]## V1 param_12723 avitoid city person_type
## 1: 1 90 2295567152 Свердловская область, Екатеринбург Агентство
## source metro param_1943
## 1: avito.ru Геологическая Сдам
## url
## 1: https://www.avito.ru/ekaterinburg/kvartiry/3-k._kvartira_110m_1216et._2295567152
## cat1_id param_2315
## 1: 1 12
## description
## 1: Квартира после косметического ремонта, в зале поставят диван, остальные подробности по телефону!\n\nВашему вниманию предлагается трёхкомнатная квартира в центре Екатеринбурга. Очень светлая, просторная, качественный ремонт, мебель, посуда, бытовая техника (стиральная машина, пылесос, утюг с парогенератором, холодильник, посудомоечная машина, варочная панель, вытяжка). Квартира расположена на пересечении улиц Радищева и Сакко и Ванцетти. Рядом крупные торговые центры `Гринвич`, `Алатырь`, парк `Зелёная роща` и дендропарк, цирк. Гимназия N5, лицей N173, школа N10, детские сады, поликлиника `УГМК-Здоровье` в шаговой доступности. Закрытый двор, детская площадка. Окна на две стороны - во двор и на пешеходную зону, благодаря чему не слышен шум центральных улиц, В квартире 3 санузла, на фото нет санузла с душевой кабиной, могу отправить отдельно. Все остальные подробности по телефону!\n\nУсловия: 55 000 + ку ( лето 6 000 , зима 10 000 ) , депозит 55 000 , по факту подписания договора комиссия 27 500, минимальный срок аренды от 6 месяцев. Квартира освободится в течение декабря, просмотры возможны уже сейчас.\n\nЛучший специалист аренды в Екатеринбурге, работаю как с выездом, так и дистанционно во время пандемии. Комиссия оплачивается только после заключения договора с собственником квартиры (возможен безналичный расчет по картам Visа/МаstеrСаrd/Маеstrо любого банка), звоните в любое удобное для вас время!
## nedvigimost_type price cat2 param_6862 contactname phone_protected
## 1: Сдам 55000 Квартиры Есть NA 0
## cat1 id km_do_metro param_6866 param_2019
## 1: Недвижимость 565176488 0.9 1 3
## person address param_2415 cat2_id
## 1: Аренда элитных квартир ул. Радищева, 33 16 2
## time title param_2016
## 1: 2022-04-22 23:32:22 3-к. квартира, 110 м², 12/16 эт. На длительный срок
## param_6865 param_2515 phone person_type_id nedvigimost_type_id
## 1: 50 110 89584900928 2 2
## source_id region city1 phone_operator
## 1: 1 Свердловская область Екатеринбург МТТ
## phone_region count_ads_same_phone param_1945 param_12724 lat
## 1: Свердловская область 6 <NA> NA 56.82849
## lng param_2078
## 1: 60.58855 <NA>Фильтруем по цене и месяцам. Выбираем только летние (туристические) месяцы.
df1 <- df %>% select(avitoid, price, lat, lng, param_2515, time)
lower_bound <- quantile(df1$price, 0.025)
upper_bound <- quantile(df1$price, 0.975)
lower_bound <- quantile(df1$param_2515, 0.025)
upper_bound <- quantile(df1$param_2515, 0.975)
df2 <- df1 %>%
filter(price >= 12000) %>%
filter(price <= 58000) %>%
filter(param_2515 >= 22) %>%
filter(param_2515 <= 100)
df2$time <- as.Date(df2$time)
df3 <- df2 %>%
filter(as.numeric(strftime(df2$time, '%m')) %in% 6:8)Создаем sf объект.
Пока в классической проекции WGS84. Функция st_difference убирает наблюдения, которые имеют одинаковые координаты. В некоторых случаях это можно использовать, но лучше попытаться как-то агрегировать. Например выгрузить полигоны домов и агрегировать несколько квартир на уровне дома.
df_sf <- st_as_sf(df3, coords = c("lng", "lat"), crs = 4326)
df_sf <- st_difference(df_sf)## although coordinates are longitude/latitude, st_difference assumes that they
## are planarmapview(df_sf)
Moran test with points
Тест Морана - это проверка на пространственную автокорреляцию в данных. О том, что это такое и почему эта тема должна быть среди первых при моделировании на пространственных данных можно почитать здесь.
Matrix W - Distance.
Перед проверкой данных на пространственную автокорреляцию необходимо определить, какие точки являются соседями (neighbours). Это можно сделать разными способами. В данном случае мы используем расстояние между двумя точками. Сначала используем проекцию Екатеринбурга. Можете попробовать посчитать в WGS84 и посмотреть что выйдет. Более того, в проекции геометрия представлена не в lat и long, а в метрах. Для построения матрицы мы берем координаты точек из sf объекта и с помощью функции dnearneigh заайем параметры соседей (точки находящиеся на расстоянии < 500 метров). Так как наши данные перепроецированы, мы указываем, что геометрия задана не в формате long и lat. С помощью функции nb2listw создаем ту самую матрицу в формате листа, так как тест Морана принимает на вход именно такой формат. Style задает параметр стандартизации, в данном случае W — row strandartization, zero.policy — параметр, который задает, добавляем ли мы точки у который нет соседей. В данном случае да.
df_sf1 <- df_sf %>% st_transform(32641)
df_coords <- st_coordinates(df_sf1)
df_nb <- dnearneigh(df_coords, 0, 500, longlat = F)
df_W <- nb2listw(df_nb, style = 'W', zero.policy = T)
df_W$neighbours## Neighbour list object:
## Number of regions: 7230
## Number of nonzero links: 440806
## Percentage nonzero weights: 0.8432783
## Average number of links: 60.96902
## 18 regions with no links:
## 18 1803 2073 3198 3222 3564 3596 5155 5545 5648 5817 5987 6290 6321
## 6336 6620 6635 6999Matrix W - K nearest.
Делаем похожую на предыдущий чанк операцию, но в этот раз строим матрицу не на основе дистанции, а на основе определенного количества ближайших точек. В данном случае это K = 10. То есть соседями являются 10 ближайших точек от данной.
list_k <- knearneigh(df_coords, k = 10)
neighbour_knear <- knn2nb(knn = list_k)
knear_W <- nb2mat(neighbours = neighbour_knear, style = "W")
plot(knn2nb(list_k), df_coords)
Matrix W - triangulatio.
Вариант для точечных слоев. С помощью триангуляции выделяем точки, полигоны которых (по сути Voronoi polygons) соприкасаются.
tri_neighbours <- tri2nb(st_geometry(df_sf1))Global Moran test.
Фукнция построения глобального индекса Морана. На вход берет sf объект, лист матрицы - df_W и тут можно задать параметр о включении точек без связей. Permutation — проверка на то, насколько полученные результаты не артефакт. Мы 999 раз (nsim) перемешиваем данные: меняем координаты точек, при этом сохраняя их параметр (цена) и смотрим, насколько изменяются параметры тестов (если сохраняются — это знак того, что найденный эффект не является следствием пространственного расположения объектов).
global_moran <- moran.test(df_sf1$price, df_W, zero.policy = T)
global_moran##
## Moran I test under randomisation
##
## data: df_sf1$price
## weights: df_W n reduced by no-neighbour observations
##
##
## Moran I statistic standard deviate = 89.167, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 2.341261e-01 -1.386770e-04 6.902485e-06permutation <- moran.mc(df_sf1$price, df_W, zero.policy = T, nsim = 999)
permutation##
## Monte-Carlo simulation of Moran I
##
## data: df_sf1$price
## weights: df_W
## number of simulations + 1: 1000
##
## statistic = 0.23413, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greaterplot(permutation)
На последнем графике смотрим на то, пересекает ли вертикальная линия (значение Moran теста на наших данных) линию распределения слева (значения теста, полученные в результате случайных перестановок). Если не пересекает (как в нашем случае), то можно сделать вывод, что найденный пространственный эффект не мог получиться случайно при прочих равных.
Residuals Moran test.
Алтернативой классическому Морану является проверка распределения остатков. Мы строим простеньку модель и смотрим есть ли пространственная корреляция в остатках. Если да, то значит, что наши данные пространственно скоррелированы, но не все параметры добавлены в модель.
price.ols <- lm(price ~ param_2515, data = df_sf1)
lm_moran <- lm.morantest(model = price.ols, listw = df_W, zero.policy = TRUE)
lm_moran##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = price ~ param_2515, data = df_sf1)
## weights: df_W
##
## Moran I statistic standard deviate = 95.83, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 2.516449e-01 -1.494024e-04 6.903758e-06Polygons.
Функция для построения полигонов Воронова. Один из вариантов работы с точечными слоями, при котором вы превращаете их в полигоны. Чтобы понять в чём суть метода, предлагаю обратиться к замечательному посту на Хабре и к двум видео-иллюстрациям ( 1 и 2)
df_sf2 <- df_sf1 %>%
st_set_geometry(value = st_voronoi(do.call(c, st_geometry(df_sf1))) %>%
st_collection_extract()) %>%
st_set_crs(32641)mapview(df_sf2)Обрежем получившуюся данные по границам Екатеринбурга
admin_osm <- opq(bbox = 'Ekaterinburg') %>%
add_osm_feature(key = "admin_level", value = '6') %>%
osmdata_sf()admin_osm <- admin_osm[["osm_multipolygons"]] %>%
dplyr::filter(name.en == 'Yekaterinburg Municipality') %>%
st_union() %>%
st_transform(32641)## Error in `stopifnot()`:
## ℹ In argument: `name.en == "Yekaterinburg Municipality"`.
## Caused by error:
## ! object 'name.en' not founddf_sf2 <- st_intersection(df_sf2, admin_osm)## Error in geos_op2_geom("intersection", x, y, ...): st_crs(x) == st_crs(y) is not TRUEmapview(df_sf2)