CAR_20180813.R

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区)")