Created
August 12, 2018 15:38
-
-
Save Gedevan-Aleksizde/3f3331bb0206e67e82823c4c51027358 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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