都内の不動産のデータをhereRを使って可視化・分析する
この記事はHERE Advent Calendar 2023の17日目の記事です。
Rでこんなパッケージを見つけた。試してみたいので記事を書くことにした。
HERE REST APIsをラップしているみたいで、ジオコーディングやら到達圏解析やらできるらしい。 作成者はスイス連邦鉄道のデータサイエンティストの方のようだ。
今回は不動産のデータをWebスクレイプしてhereRパッケージを使いつつ可視化・分析する。
データの取得
以下の記事を参考にし、suumoに掲載されているデータを取得させていただいた。
大量にデータを取得するとサーバに負荷がかかるので、取得するデータは都内のワンルームマンションだけにする。
なお、当該ウェブサイトの利用規約上に著作権に係る規定や禁止行為が明記されているので、当記事を参考にされる方がいれば事前に参照されたい。
以下のステップでWebスクレイプしている。
RでWebスクレイピングする際はrvestを使うと便利だ。 以下のコードでは、rvestに加えてこの後使うパッケージもロードし、スクレイプした各ページのデータをまとめた後に列名を適切なものに変更している。
# import packages to use require(rvest) require(tidyverse) require(ggthemes) require(DataExplorer) require(patchwork) require(hereR) require(tmap) require(tmaptools) require(stringi) require(sf) require(geojsonsf) # Scraping data from SUUMO base_url <- ("https://suumo.jp/jj/bukken/ichiran/JJ010FJ001/?ar=030&bs=011&ta=13&jspIdFlg=patternShikugun&kb=1&kt=9999999&mb=0&mt=9999999&md=0&ekTjCd=&ekTjNm=&tj=0&cnb=0&cn=9999999&srch_navi=1") html <- read_html(base_url) max_page_num <- read_html(base_url) %>% html_elements(".pagination-parts") %>% html_elements("a") %>% html_text2() %>% as.integer() %>% max(na.rm = TRUE) %>% as.character() header <- c("name", "price", "address", "line_station", "area", "room_type", "balcony", "year_built") page_urls <- str_glue(str_c(base_url, "&page={1:", max_page_num, "}")) scraper <- function(url) { read_html(url) %>% html_elements(".dottable--cassette") %>% html_elements("dd") %>% html_text2() %>% matrix(ncol = 8, byrow = TRUE) %>% as_tibble() } dat <- map(page_urls, scraper) %>% list_rbind() %>% rename_with(~ header)
データの整形
取得したデータでは価格に漢数字が使われていたり、最寄り駅の路線と駅名が同じカラムに入っていたりするので整形する。 数値は文字列になっているのでデータ型も変更し、築年数も計算している。
# Data cleansing dat <- dat %>% separate(col = price, sep = "億", into = c("hundread_M", "ten_K"), fill = "left") %>% mutate(hundread_M = as.integer(hundread_M) * 100000000) %>% mutate(ten_K = str_extract(ten_K, "\\d+")) %>% mutate(ten_K = as.integer(ten_K) * 10000) %>% replace_na(replace = list(hundread_M = 0, ten_K = 0)) %>% mutate(price = hundread_M + ten_K) %>% separate(line_station, sep = "「", into = c("line", "station")) %>% separate(station, sep = "」", into = c("station", "walk_from_station")) %>% mutate(walk_from_station = str_extract(walk_from_station, "\\d+")) %>% mutate(walk_from_station = as.integer(walk_from_station)) %>% mutate(area = str_extract(area, "\\d+(?:\\.\\d+)?")) %>% mutate(area = as.numeric(area)) %>% mutate(year_built = str_replace(year_built, pattern = "年", replacement = "-")) %>% mutate(year_built = str_remove(year_built, pattern = "月")) %>% mutate(year_built = as_date(year_built, format = ("%Y-%m"))) %>% mutate(age_of_building = as.integer((Sys.Date() - year_built)/365)) %>% select(-c(hundread_M, ten_K, balcony, room_type))
整形後のしたデータは以下のようなものになる。
データの可視化（チャート）
データは用意できたのでとりあえず可視化していく。 チャートを作成する際はggplot2に加えてDataExplorerを使うと便利だ。
# Data exploration plot_histogram(dat) plot_bar(dat) plot_correlation(dat) g1 <- ggplot(dat, aes(x = area, y = price)) + scale_x_continuous(breaks = seq(0, 160, 10)) + geom_point() + theme_economist() g2 <- ggplot(dat, aes(x = age_of_building, y = price)) + geom_point() + theme_economist() g3 <- ggplot(dat, aes(x = walk_from_station, y = price)) + scale_x_continuous(breaks = seq(0, 20, 5)) + geom_point() + theme_economist() g1 + g2 + g3 + plot_layout(ncol = 1)
ヒストグラムをみると築40年以上の物件がかなり多い。また、ワンルームなのにかなり広い物件が少数あるようだ。おそらく同じ物件だろうが価格も億を超える物件がある。
各変数事の相関係数を見ると価格と専有面積の相関係数が0.73で、築年数と最寄駅からの徒歩時間も若干負の相関がある。
築年数の影響をもう少し確認したいので、グループ化して中央値をとる。
group_by_ages <- dat %>% mutate(age_of_building_class = case_when(age_of_building <= 5 ~ "<=5", age_of_building > 5 & age_of_building <= 10 ~ "<=10", age_of_building > 10 & age_of_building <= 15 ~ "<=15", age_of_building > 15 & age_of_building <= 20 ~ "<=20", age_of_building > 20 & age_of_building <= 25 ~ "<=25", age_of_building > 25 & age_of_building <= 30 ~ "<=30", age_of_building > 30 & age_of_building <= 35 ~ "<=35", age_of_building > 35 & age_of_building <= 40 ~ "<=40", age_of_building > 40 ~ ">40")) %>% group_by(age_of_building_class) %>% summarise(median_price = median(price)) plot_bar(group_by_ages, "median_price")
都内の物件であればある程度価格を維持していると思うがさすがに、25年を超えると下がってくる傾向がありそうだ。また、40年を超えた物件はかなり中央値が下がっている。
データの可視化（マップ）
マップで可視化するので、hereRを使ってジオコーディングする。 マップの表示にはtmapを使う。
ベースマップもカスタムできるので、HEREのラスタータイルを使用しよう。 以下のコードでは試しに「東京都」でジオコーディングして返却されたポイントを可視化している。データはsf形式で返却されるので、そのままtmapの関数にデータを渡すことができる。
# Map visualization key <- "YOUR_API_KEY" set_key(key) basemap_endpoint <- paste0("https://maps.hereapi.com/v3/base/mc/{z}/{x}/{y}/png?style=lite.day&apiKey=", key) tokyo <- geocode(address = "東京都") tmap_mode("view") tm_shape(tokyo) + tm_dots(col = "#CC6666", size = 0.15) + tm_basemap(server = basemap_endpoint)
東京都でジオコーディングすると千代田区役所が代表点として返却されるようだ。
以下のコードでWebスクレイプで取得したデータの住所を元にジオコーディングし、返却された結果を元のデータに結合している。 返却されるデータには結合用にidが付与されるため、元のデータにも事前にidを付与している。 地図での表示には価格を元に比例シンボルで表示している。
# Geocode scraped data dat <- rowid_to_column(dat) %>% rename(id = rowid) geocoded <- geocode(dat$address) %>% right_join(dat, by = "id") %>% rename(geocoded_address = address.x, source_address = address.y) %>% select(id, rank, geocoded_address, type, name, source_address, line, station, walk_from_station, area, year_built, price, age_of_building, geometry) tm_shape(geocoded) + tm_bubbles(size = "price", col = "#CC6666", alpha = 0.7, popup.vars = TRUE) + tm_basemap(server = basemap_endpoint)
当然だが都心部の方が価格が高い物件が多い。
次に、各市区ごとの価格の中央値を地図上で可視化してみる。スタンフォード大学が日本の市区町村の行政区域のポリゴンを公開してくれているので、今回はそれを使った。
Second-level Administrative Divisions, Japan, 2015 | Stanford Digital Repository
GeoJSON形式でDLした行政区域のデータを読み込み、東京だけフィルターして使う。 sfパッケージのst_join関数で空間結合し、各市区町村の不動産価格の中央値を集計し表示している。
# Visualize median price for each city in Tokyo admin_poly <- geojson_sf("share_w_host/here_advent_calendar_2023/stanford-ff095jm2871-geojson.json") %>% filter(name_1 == "Tokyo") tm_shape(admin_poly) + tm_polygons(col = "#CC6666", alpha = 0.5) + tm_basemap(server = basemap_endpoint) admin_medion_price <- st_join(admin_poly, y = geocoded) %>% group_by(name_2) %>% summarise(median_price = median(price)) tm_shape(admin_medion_price) + tm_polygons(col = "median_price", alpha = 0.5, ) + tm_shape(geocoded) + tm_bubbles(size = "price", col = "green", alpha = 0.2, popup.vars = TRUE) + tm_basemap(server = basemap_endpoint)
やはり都心三区を中心に価格が高い傾向があるが、練馬区や足立区、府中なども周辺に比べると高い中央値になっている。
一般化線形モデルによる統計解析
最後に一般化線形モデルで統計解析する。価格を目的変数にして、専有面積、駅からの徒歩時間、築年数、所在する市区町村を都心三区、都心六区、その他に分けたグループを説明変数にする。最寄駅を説明変数にしてもよいが、カテゴリーデータが多くなりすぎるとオーバーフィットしそうなので、上記三つにグループ分けして説明変数にすることにした。 モデルの確率密度関数は非負の連続データなので、とりあえずガンマ関数を使うことにする。
以下のコードで市区町村のグループ分け、一般化線形モデルでの統計解析、サマリの表示を行っている。
# Statistical analysis dat_for_model <- st_join(x = geocoded, y = admin_poly) %>% mutate(region = case_when(nl_name_2 == "千代田区" | nl_name_2 == "中央区" | nl_name_2 == "港区" ~ "都心3区", nl_name_2 == "渋谷区" | nl_name_2 == "新宿区" | nl_name_2 == "文京区" ~ "都心6区", .default = "その他")) model <- glm(price ~ area + walk_from_station + age_of_building + region, data = dat_for_model, family = Gamma(link = "log")) summary(model)
結果は以下の通りだ。
以下のコードでモデルを元に予測した価格を表示し、残差を確認している。
# Prediction predictions <- exp(predict(model, newdata = dat_for_model, type = "link")) # Combine predictions with original data pred_data <- cbind(dat_for_model, pred_price = predictions) # Regression Plot ggplot(pred_data, aes(x = area, y = pred_price)) + geom_point(aes(y = price), color = "#CC6666", alpha = 0.5) + geom_point(aes(y = pred_price), color = "lightgreen", alpha = 0.5) + geom_smooth(method = "glm", formula = y ~ x, se = FALSE, color = "grey", method.args = list(family = "Gamma")) + labs(x = "Area", y = "Predicted Price") + theme_economist_white() # Residual Plot residuals <- residuals(model, type = "deviance") ggplot(pred_data, aes(x = pred_price, y = residuals)) + geom_point(color = "orange", alpha = 0.5) + geom_hline(yintercept = 0, linetype = "dashed", color = "blue") + labs(x = "Predicted Price", y = "Residuals") + theme_economist_white()
予測した価格はピンク色、実際のデータは緑色のドットで表示しているが、専有面積が大きくなると価格が一気に上昇するモデルになっている。おそらくは都心にある外れ値のデータが影響しているのだろう。
残差プロットを見ても外れ値が存在することがわかる。ガンマ関数は外れ値に対してロバストではないので、他の確率密度関数を用いるか外れ値の処理をする方がよさそうだ。
とりあえず今回はここまで。
