Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
You can’t perform that action at this time.