require(tidyverse) # ver. 1.2.1 | |
require(CARBayes) # ver. 5.0 | |
require(sf) # ver. 0.6-3 | |
require(ggplot2) # ver. 3.0.0 | |
require(ggthemes) # ver. 4.0.0 | |
require(ggmcmc) # ver. 1.1 | |
rmse <- function(y, pred){ | |
return(sqrt(mean((y - pred)^2))) | |
} | |
# voronoi <- readShapePoly("out.shp") | |
voronoi <- st_read("out.shp") %>% mutate(obs_id=1:NROW(voronoi)) | |
W <- st_intersects(voronoi, voronoi, sparse=F) * 1 | |
ggplot(voronoi) + geom_sf(aes(fill=log(price)), size=.1) + theme_tufte() | |
summary(voronoi) | |
ggplot(voronoi, aes(x=price)) + geom_histogram() + scale_x_log10() + theme_tufte() + labs(title="地価分布(対数)") | |
# intrinsic CAR | |
fit <- S.CARleroux( | |
log(price) ~ is_hosp + is_shop + is_plant + region_2 + size, data=voronoi, | |
family = "gaussian", W=W, fix.rho=T, rho=1, | |
burnin = 2000, n.sample = 32000, thin=10) | |
params <- c("beta") | |
for(name in params){ | |
print(name) | |
ggmcmc(ggs(fit$samples[[name]]), | |
file=NULL, plot=c("traceplot", "histogram", "autocorrelation")) | |
} | |
coef(fit) | |
voronoi <- voronoi %>% mutate(predict = exp(fit$fitted.values)) %>% | |
mutate(residual=price - predict) | |
print(rmse(voronoi$price, voronoi$predict)) | |
summary(voronoi) | |
ggplot(voronoi, aes(x=residual)) + geom_histogram() + theme_tufte() + labs(title="CAR モデル予測残差") | |
ggplot(voronoi) + geom_sf(aes(fill=residual), size=.1) + theme_tufte() + labs(title="CAR モデル予測残差") | |
ggplot(voronoi %>% filter(grepl("区$", region_2)), aes(x=residual)) + geom_histogram() + | |
theme_tufte() + labs(title="CAR モデル予測残差 (23区)") | |
ggplot(voronoi %>% filter(grepl("区$", region_2))) + geom_sf(aes(fill=residual), size=.1) + | |
theme_tufte() + labs(title="CAR モデル予測残差 (23区)") | |
rail <- read_sf("N02-17_GML/N02-17_Station.shp", | |
options = "ENCODING=CP932") %>% | |
filter(N02_004 %in% c("東日本旅客鉄道")) %>% | |
st_intersection(y=voronoi %>% filter(grepl("区$", region_2))) | |
# how to display labels | |
# https://notchained.hatenablog.com/entry/2018/05/28/003910 | |
# rail_coords <- st_transform(rail %>% dplyr::select(N02_005), st_crs(rail)) | |
# rail_coords <- cbind(rail_coords, st_point_on_surface(rail_coords$geometry) %>% st_coordinates) | |
ggplot(voronoi %>% filter(grepl("区$", region_2))) + | |
geom_sf(aes(fill=residual), size=.1) + | |
geom_sf(data=rail, color="white") + | |
# geom_label(aes(X, Y, label=N02_005), data=rail_coords) + | |
theme_tufte() | |
# add station information | |
rail <- read_sf("N02-17_GML/N02-17_Station.shp", | |
options = "ENCODING=CP932") %>% | |
filter(N02_004 %in% c("東日本旅客鉄道")) %>% | |
st_intersection(y=voronoi) | |
voronoi <- st_intersection( | |
voronoi, rail_coords) %>% mutate(has_station=T) %>% | |
as.data.frame %>% dplyr::select(obs_id, has_station) %>% unique %>% left_join(x=voronoi, y=., on=id) %>% | |
mutate(has_station=if_else(is.na(has_station), F, T)) | |
voronoi <- st_intersection( | |
voronoi, rail_coords %>% filter(N02_005 %in% c("東京", "池袋", "渋谷", "新宿"))) %>% mutate(has_hub_station=T) %>% | |
as.data.frame %>% dplyr::select(obs_id, has_hub_station) %>% unique %>% left_join(x=voronoi, y=., on=id) %>% | |
mutate(has_hub_station=if_else(is.na(has_hub_station), F, T)) | |
# reestimate by CAR model | |
fit2 <- S.CARleroux( | |
log(price) ~ has_station + has_hub_station+ is_hosp + is_shop + is_plant + region + size, data=voronoi, | |
family = "gaussian", W=W, fix.rho=T, rho=1, | |
burnin = 2000, n.sample = 32000, thin=10) | |
params <- c("beta") | |
for(name in params){ | |
print(name) | |
ggmcmc(ggs(fit2$samples[[name]]), | |
file=NULL, plot=c("traceplot", "histogram", "autocorrelation")) | |
} | |
coef(fit2) | |
voronoi <- voronoi %>% mutate(predict = exp(fit2$fitted.values)) %>% | |
mutate(residual=price - predict) | |
print(rmse(voronoi$price, voronoi$predict)) | |
ggplot(voronoi, aes(x=residual)) + geom_histogram() + theme_tufte() + labs(title="CAR モデル予測残差 (with 駅情報)") | |
ggplot(voronoi) + geom_sf(aes(fill=residual), size=.1) + theme_tufte() + labs(title="CAR モデル予測残差 (with 駅情報)") | |
ggplot(voronoi %>% filter(grepl("区$", region_2)), aes(x=residual)) + geom_histogram() + | |
theme_tufte() + labs(title="CAR モデル予測残差 (23区)") | |
ggplot(voronoi %>% filter(grepl("区$", region_2))) + geom_sf(aes(fill=residual), size=.1) + | |
geom_sf(data=rail, color="white") + | |
theme_tufte() + labs(title="CAR モデル予測残差 (23区)") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment