3-1_地震リスクと地価
Mizuki Kawabata & Keisuke Takano
ステップ0：パッケージの読み込み
library(sf)
library(tidyverse)
library(ggplot2)
library(stargazer)
library(lmtest)ステップ1：国土数値情報の「地価公示」データのダウンロード
省略。
ステップ2：地価（住宅地）データの作成
#地価データ（東京都）
#文字コードをCP932で、shapefile内部の変数をfactor化せずに読み込む
rchika <- sf::st_read(dsn="L01-21_13.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932")
#変数citycodeの作成
rchika <- rchika %>%
#変数citycodeを、変数L01_021を数値型に変換したもので定義
dplyr::mutate(citycode=as.numeric(L01_021))
#島嶼部を除く観測地点に絞る
rchika <- rchika %>%
#変数citycodeが13308以下の観測地点に絞る
dplyr::filter(citycode<=13308)
#平面直角座標系第9系（JGD2011）に投影変換
rchika <- rchika %>%
#平面直角座標系第9系（JGD2011）のEPSGは6677
sf::st_transform(crs=sf::st_crs(6677))
#住宅地の観測地点に絞る
rchika <- rchika %>%
#変数L01_001が000に等しい観測地点に絞る
dplyr::filter(L01_001=="000")
#観測地点固有IDを追加する
rchika <- rchika %>%
#ID変数（id_lp）を追加する
dplyr::mutate(id_lp=dplyr::row_number())#上の処理をパイプで繋ぐことで、一気に実行できる
rchika <- sf::st_read(dsn="L01-21_13.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932") %>%
#変数citycodeを、変数L01_021を数値型に変換したもので定義
dplyr::mutate(citycode=as.numeric(L01_021)) %>%
#変数citycodeが13308以下の観測地点に絞る
dplyr::filter(citycode<=13308) %>%
#平面直角座標系第9系（JGD2011）のEPSGは6677
sf::st_transform(crs=sf::st_crs(6677)) %>%
#変数L01_001が000に等しい観測地点に絞る
dplyr::filter(L01_001=="000") %>%
#ID変数（id_lp）を追加する
dplyr::mutate(id_lp=dplyr::row_number())#作成した地価ポイントをプロット
ggplot2::ggplot() +
ggplot2::geom_sf(data=rchika)
ステップ3：使用データの追加
#市区町村ポリゴン
municipality_tokyo <- sf::st_read(dsn="municipality_tokyo.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=UTF-8")
#鉄道駅ライン
station <- sf::st_read(dsn="station.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=UTF-8")
#地震による地域危険度ポリゴン
kikendo <- sf::st_read(dsn="kikendo.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=UTF-8") %>%
sf::st_transform(crs=sf::st_crs(6677))#読み込んだデータを重ねてプロット
ggplot2::ggplot() +
#市区町村ポリゴン（白で描画）
ggplot2::geom_sf(data=municipality_tokyo,fill="white") +
#駅ライン
ggplot2::geom_sf(data=station) +
#地価ポイント（青で描画）
ggplot2::geom_sf(data=rchika,color="blue")
ステップ4：都心主要駅の作成
#駅ラインに駅IDを変数id_stnとして追加
station <- station %>%
#通し番号を変数として追加
dplyr::mutate(id_stn=dplyr::row_number())
#都心主要駅を駅ラインから抽出
major4stn <- station %>%
dplyr::filter((N02_005=="新宿"&N02_003=="山手線"&id_stn==750)|
(N02_005=="渋谷"&N02_003=="山手線")|
(N02_005=="池袋"&N02_003=="山手線")|
(N02_005=="東京"&N02_003=="総武線"))ggplot2::ggplot() +
ggplot2::geom_sf(data=municipality_tokyo) +
ggplot2::geom_sf(data=station) +
#都心主要駅ラインを赤＆太線（引数sizeで指定）で描画
ggplot2::geom_sf(data=major4stn,color="red",size=3)
ステップ5：地価データに都心主要駅までの距離を追加
#各都心主要駅ラインの重心ポイントを作成
major4stn_cent <- major4stn %>%
sf::st_centroid()
#重心ポイントの座標値を取り出し
major4stn_cent_coord <- sf::st_coordinates(major4stn_cent)
#重心ポイントに座標値を変数として追加
major4stn_cent <- major4stn_cent %>%
dplyr::mutate(x_major4stn=major4stn_cent_coord[,"X"],
y_major4stn=major4stn_cent_coord[,"Y"])#重心ポイント作成〜座標値の追加もパイプを用いて一気に実行できる
major4stn_cent <- major4stn %>%
sf::st_centroid() %>%
dplyr::mutate(x_major4stn=sf::st_coordinates(.)[,"X"],
y_major4stn=sf::st_coordinates(.)[,"Y"])#地価ポイントについても座標値を変数として追加
rchika <- rchika %>%
dplyr::mutate(x_rchika=sf::st_coordinates(.)[,"X"],
y_rchika=sf::st_coordinates(.)[,"Y"])#地価ポイントに都心主要駅重心ポイントを空間結合
rchika2 <- rchika %>%
sf::st_join(y=major4stn_cent,join=st_nearest_feature)#空間結合された最寄主要駅の座標値と、地価ポイントの座標値を用いて、最寄主要駅までの直線距離を計算
rchika2 <- rchika2 %>%
dplyr::mutate(dis2m4stn=sqrt((x_rchika-x_major4stn)^2+(y_rchika-y_major4stn)^2))ステップ6：地価データに最寄り駅までの距離を追加
#最寄り駅までの距離変数を文字列型から数値型に変換し、新たな変数dis2stnとして追加
rchika2 <- rchika2 %>%
dplyr::mutate(dis2stn=as.numeric(L01_046))#変数dis2stnの値が0の場合、40に置換する
rchika2 <- rchika2 %>%
dplyr::mutate(dis2stn=ifelse(dis2stn==0,40,dis2stn))ステップ7：建物倒壊危険度の地図作成
ggplot2::ggplot() +
#市区町村ポリゴン
ggplot2::geom_sf(data=municipality_tokyo) +
#危険度ポリゴンを、変数B_rankに基づき色分けしてプロット
ggplot2::geom_sf(data=kikendo,aes(fill=B_rank),lty=0) +
#変数B_rankの最小値では白、最大値では赤で塗りつぶすようにグラデーションを設定
ggplot2::scale_fill_gradient(low="white",high="red",na.value=NA)
ステップ8：地価データに地域危険度データを追加
#地価ポイントに危険度ポリゴンの情報を空間結合
rchika3 <- rchika2 %>%
sf::st_join(y=kikendo,join=st_intersects)ステップ9：地域危険度のないデータの削除
#変数B_rankが0でない観測地点のみ残す
rchika4 <- rchika3 %>%
dplyr::filter(B_rank!=0)ステップ10：モデル1の変数の準備
rchika4 <- rchika4 %>%
#各種変数を数値型に変換した上で、新たな変数として追加
dplyr::mutate(rchika=as.numeric(L01_006),
chiseki=as.numeric(L01_024),
FAR=as.numeric(L01_053)) %>%
#変数B_rankの値に応じたダミー変数を追加
dplyr::mutate(BriskL=ifelse(B_rank==1|B_rank==2,1,0),
BriskM=ifelse(B_rank==3,1,0),
BriskH=ifelse(B_rank==4|B_rank==5,1,0))ステップ11：モデル1の推定
model1 <- rchika ~ chiseki + FAR + dis2m4stn + dis2stn + BriskM + BriskH
lm_model1 <- lm(model1,data=rchika4)
summary(lm_model1)##
## Call:
## lm(formula = model1, data = rchika4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1235047 -113168 -22490 79064 3634037
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.518e+05 2.345e+04 19.267 < 2e-16 ***
## chiseki 9.338e+01 1.420e+01 6.574 6.5e-11 ***
## FAR 1.772e+03 9.301e+01 19.053 < 2e-16 ***
## dis2m4stn -1.656e+01 7.756e-01 -21.351 < 2e-16 ***
## dis2stn -2.124e+01 8.230e+00 -2.581 0.00994 **
## BriskM -2.730e+05 2.172e+04 -12.566 < 2e-16 ***
## BriskH -3.629e+05 3.300e+04 -10.998 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277900 on 1683 degrees of freedom
## Multiple R-squared: 0.5318, Adjusted R-squared: 0.5302
## F-statistic: 318.6 on 6 and 1683 DF, p-value: < 2.2e-16ステップ12：散布図の作成
ggplot2::ggplot(data=rchika4,aes(x=dis2m4stn,y=rchika)) +
ggplot2::geom_point()
ggplot2::ggplot(data=rchika4,aes(x=log(dis2m4stn),y=log(rchika))) +
ggplot2::geom_point()
ステップ13：自然対数を用いたモデル2の推定
model2 <- log(rchika) ~ log(chiseki) + FAR + log(dis2m4stn) + log(dis2stn) + BriskM + BriskH
lm_model2 <- lm(model2,data=rchika4)
summary(lm_model2)##
## Call:
## lm(formula = model2, data = rchika4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41855 -0.18206 0.01104 0.19790 1.26524
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.3232526 0.1479193 130.634 < 2e-16 ***
## log(chiseki) 0.1259213 0.0139435 9.031 < 2e-16 ***
## FAR 0.0005656 0.0001188 4.761 2.09e-06 ***
## log(dis2m4stn) -0.6127471 0.0119311 -51.357 < 2e-16 ***
## log(dis2stn) -0.2426033 0.0131856 -18.399 < 2e-16 ***
## BriskM -0.1405952 0.0256531 -5.481 4.88e-08 ***
## BriskH -0.1753882 0.0390919 -4.487 7.73e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3239 on 1683 degrees of freedom
## Multiple R-squared: 0.8121, Adjusted R-squared: 0.8114
## F-statistic: 1212 on 6 and 1683 DF, p-value: < 2.2e-16ステップ14：残差の可視化
#model2の残差を変数m2residとしてデータセットに追加
rchika4m2 <- rchika4 %>%
dplyr::mutate(m2resid=lm_model2$residuals)
ggplot2::ggplot() +
ggplot2::geom_sf(data=municipality_tokyo,fill="white") +
#残差が結合された地価ポイントをプロット（残差の色分け階級は、最小値から最大値まで標準偏差で分級）
ggplot2::geom_sf(data=rchika4m2,
aes(color=cut(m2resid,breaks=c(min(m2resid),-2*sd(m2resid),-sd(m2resid),0,
sd(m2resid),2*sd(m2resid),max(m2resid)))),
lty=0) +
#階級の色分けには「Spectral」というパレットを使う
ggplot2::scale_color_brewer(palette="Spectral",direction=-1) +
#図の凡例のタイトルを手動で指定する
ggplot2::guides(color=ggplot2::guide_legend(title="m2resid"))
ステップ15：地域ダミー変数の作成
rchika4m2 <- rchika4m2 %>%
#変数city_nameに基づくダミー変数の作成
dplyr::mutate(chiyoda=ifelse(city_name=="千代田区",1,0),
minato=ifelse(city_name=="港区",1,0),
toshima=ifelse(city_name=="豊島区",1,0),
oume=ifelse(city_name=="青梅市",1,0))吉祥寺駅周辺ダミーの作成。
kichijo <- station %>%
#駅ラインから吉祥寺駅を取り出し
dplyr::filter(N02_005=="吉祥寺") %>%
#重心ポイントを作成
sf::st_centroid()
#重心ポイントから1500mバッファを作成
kichijo1500m <- sf::st_buffer(x=kichijo,dist=1500)
#地価ポイントを1500mバッファでクリップ
rchika4m2_kichijo1500m <- sf::st_intersection(x=rchika4m2,
y=kichijo1500m)
rchika4m2_kichijo1500m <- rchika4m2_kichijo1500m %>%
#クリップされた地価ポイントのgeometry情報を削除
sf::st_set_geometry(NULL) %>%
#観測地点ID（id_lp）のみを変数として残した上で、値を一意にする
dplyr::distinct(id_lp)
rchika4m2_kichijo1500m <- rchika4m2_kichijo1500m %>%
#全てのレコードで1を取る変数kichijoを追加
dplyr::mutate(kichijo=1)上と同じ要領で、国立駅周辺ダミーを作成。パイプを使うことで、シンプルに書くことが可能。
kunitachi <- station %>%
dplyr::filter(N02_005=="国立") %>%
sf::st_centroid()
kunitachi1500m <- sf::st_buffer(x=kunitachi,dist=1500)
rchika4m2_kunitachi1500m <- sf::st_intersection(x=rchika4m2,
y=kunitachi1500m) %>%
sf::st_set_geometry(NULL) %>%
dplyr::distinct(id_lp) %>%
dplyr::mutate(kunitachi=1)rchika4m3 <- rchika4m2 %>%
#吉祥寺駅周辺ダミーkichijoを地価ポイントに左結合
dplyr::left_join(y=rchika4m2_kichijo1500m,by="id_lp") %>%
#国立駅周辺ダミーkunitachiを地価ポイントに左結合
dplyr::left_join(y=rchika4m2_kunitachi1500m,by="id_lp") %>%
#変数kichijo,kunitachiについて、値がNAならば0に置換、NAでなければ元の値を保持
dplyr::mutate(across(c(kichijo,kunitachi),~ifelse(is.na(.),0,.)))ステップ16：地域ダミー変数を追加したモデル3の推定
model3 <- log(rchika) ~ log(chiseki) + FAR + log(dis2m4stn) + log(dis2stn) + BriskM + BriskH + chiyoda + minato + kichijo + kunitachi + toshima + oume
lm_model3 <- lm(model3,data=rchika4m3)
summary(lm_model3)##
## Call:
## lm(formula = model3, data = rchika4m3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10396 -0.17413 0.00896 0.18541 0.81656
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.4421893 0.1326590 146.558 < 2e-16 ***
## log(chiseki) 0.0942566 0.0122521 7.693 2.43e-14 ***
## FAR 0.0005118 0.0001064 4.812 1.63e-06 ***
## log(dis2m4stn) -0.6149235 0.0109605 -56.103 < 2e-16 ***
## log(dis2stn) -0.2336642 0.0114984 -20.321 < 2e-16 ***
## BriskM -0.0981469 0.0225909 -4.345 1.48e-05 ***
## BriskH -0.1640730 0.0343548 -4.776 1.95e-06 ***
## chiyoda 0.6844991 0.1098282 6.232 5.80e-10 ***
## minato 0.6240759 0.0547225 11.404 < 2e-16 ***
## kichijo 0.6753265 0.0711852 9.487 < 2e-16 ***
## kunitachi 0.5709381 0.0820303 6.960 4.86e-12 ***
## toshima -0.7524259 0.0620837 -12.120 < 2e-16 ***
## oume -0.4423210 0.0560350 -7.894 5.26e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.282 on 1677 degrees of freedom
## Multiple R-squared: 0.858, Adjusted R-squared: 0.857
## F-statistic: 844.5 on 12 and 1677 DF, p-value: < 2.2e-16#推定結果を表形式にまとめる
stargazer::stargazer(lm_model1,lm_model2,lm_model3,type="text",aligh=TRUE)##
## ====================================================================================================
## Dependent variable:
## --------------------------------------------------------------------------------
## rchika log(rchika)
## (1) (2) (3)
## ----------------------------------------------------------------------------------------------------
## chiseki 93.378***
## (14.203)
##
## log(chiseki) 0.126*** 0.094***
## (0.014) (0.012)
##
## FAR 1,772.128*** 0.001*** 0.001***
## (93.009) (0.0001) (0.0001)
##
## dis2m4stn -16.560***
## (0.776)
##
## dis2stn -21.240***
## (8.230)
##
## log(dis2m4stn) -0.613*** -0.615***
## (0.012) (0.011)
##
## log(dis2stn) -0.243*** -0.234***
## (0.013) (0.011)
##
## BriskM -272,966.500*** -0.141*** -0.098***
## (21,722.040) (0.026) (0.023)
##
## BriskH -362,938.800*** -0.175*** -0.164***
## (32,999.930) (0.039) (0.034)
##
## chiyoda 0.684***
## (0.110)
##
## minato 0.624***
## (0.055)
##
## kichijo 0.675***
## (0.071)
##
## kunitachi 0.571***
## (0.082)
##
## toshima -0.752***
## (0.062)
##
## oume -0.442***
## (0.056)
##
## Constant 451,832.600*** 19.323*** 19.442***
## (23,451.620) (0.148) (0.133)
##
## ----------------------------------------------------------------------------------------------------
## Observations 1,690 1,690 1,690
## R2 0.532 0.812 0.858
## Adjusted R2 0.530 0.811 0.857
## Residual Std. Error 277,903.700 (df = 1683) 0.324 (df = 1683) 0.282 (df = 1677)
## F Statistic 318.635*** (df = 6; 1683) 1,212.351*** (df = 6; 1683) 844.538*** (df = 12; 1677)
## ====================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
##
## ====
## TRUE
## ----#99%信頼区間
confint(lm_model3,level=0.99)## 0.5 % 99.5 %
## (Intercept) 19.1000930649 19.7842855941
## log(chiseki) 0.0626612095 0.1258519835
## FAR 0.0002375676 0.0007861213
## log(dis2m4stn) -0.6431881116 -0.5866588019
## log(dis2stn) -0.2633160205 -0.2040124592
## BriskM -0.1564035170 -0.0398903579
## BriskH -0.2526659098 -0.0754800258
## chiyoda 0.4012781905 0.9677200578
## minato 0.4829595044 0.7651922014
## kichijo 0.4917566100 0.8588963323
## kunitachi 0.3594013273 0.7824748788
## toshima -0.9125250705 -0.5923266951
## oume -0.5868219092 -0.2978200319#BP検定
bptest(lm_model3)##
## studentized Breusch-Pagan test
##
## data: lm_model3
## BP = 54.96, df = 12, p-value = 1.84e-07