3-2_水害リスクと地価
Shohei Yasuda & Keisuke Takano
ステップ0：パッケージの読み込み
library(sf)
library(tidyverse)
library(fixest)
library(ggplot2)
library(magrittr)
library(dummy)ステップ1：国土数値情報のデータのダウンロード
省略。
ステップ2：座標系の設定
作業ディレクトリに設置した地価データと浸水想定区域データを読み込む。その際、JGD2000からJGD2011への投影変換と、JGD2011から平面直角座標第9系への投影変換も同時に行う。地価データにはID変数FIDを付与する。
#地価データ（茨城県）
#文字コードをCP932で、shapefile内部の変数をfactor化せずに読み込む
chika_ibaraki <- sf::st_read(dsn="L01-21_08.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932") %>%
#JGD2011への投影変換
sf::st_transform(crs=sf::st_crs(6668)) %>%
#平面直角座標第9系への投影変換
sf::st_transform(crs=sf::st_crs(6677)) %>%
#各地価観測地点にIDを付与
dplyr::mutate(FID=dplyr::row_number())
#浸水想定区域データ（茨城県）
#文字コードをCP932で、shapefile内部の変数をfactor化せずに読み込む
shinsui <- sf::st_read(dsn="A31-12_08.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932") %>%
#JGD2011への投影変換
sf::st_transform(crs=sf::st_crs(6668)) %>%
#平面直角座標第9系への投影変換
sf::st_transform(crs=sf::st_crs(6677)) %>%
#不整地物の整形
sf::st_make_valid(shinsui) %>%
#【重要】1つの地点が複数の浸水想定区域に指定されている場合があるので、そのような重複を削除する
#ポリゴンの融合を行うことによって、重なったポリゴンの重複を削除できる
dplyr::mutate(d=1) %>%
dplyr::group_by(d) %>%
dplyr::summarise()読み込んだデータの可視化。地価ポイントが最上位に来るように描画。
ggplot2::ggplot() +
ggplot2::geom_sf(data=shinsui) +
ggplot2::geom_sf(data=chika_ibaraki)
ステップ3：浸水想定区域のダミー変数の作成
以下の手順で浸水想定区域ダミーを作成する。
- 浸水想定区域ポリゴンと地価ポイントの共通部分を取る。これにより、浸水想定区域に含まれる地価観測地点が特定できる。
- 共通部分として抽出されたポイントからgeometry属性を削除し、浸水想定区域ダミーd_shinsuiを作成した上で、不要な変数を削除する。
- FIDをキー変数として、d_shinsuiを地価ポイントに左結合する。
#浸水想定区域に含まれる地価ポイントの抽出
chika_ibaraki_shinsui <- sf::st_intersection(chika_ibaraki,shinsui) %>%
#geometry属性を削除
sf::st_set_geometry(value=NULL) %>%
#浸水想定区域ダミーd_shinsuiを作成
dplyr::mutate(d_shinsui=1) %>%
#不要な変数を削除
dplyr::select(FID,d_shinsui)
#d_shinsuiを地価ポイントに左結合する
chika_ibaraki <- chika_ibaraki %>%
dplyr::left_join(chika_ibaraki_shinsui,by="FID") %>%
#d_shinsuiが結合されなかった（＝区域外）地価ポイントについて、d_shinsuiの値を0に置換
dplyr::mutate(d_shinsui=ifelse(!is.na(d_shinsui),d_shinsui,0))ステップ4：浸水被害データの加工
「常総地区の推定浸水範囲_201509101800.kml」は中身がRで読み込めない状態（部分的に破損）していた為、それ以外の3つのkmlファイルのみを読み込み、加工を行う。なお、同じ処理を何度も繰り返しコーディングするのを避ける為、読み込み・投影変換・ポリゴン化の一連の処理を関数化した上で処理を実行する。
#浸水被害データの読み込み・投影変換・ポリゴン化を一括して行う為の関数を定義
read_polygonize <- function(x){
#文字コードをCP932で、kml内部の変数をfactor化せずに読み込む
dat <- sf::st_read(dsn=x,
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932") %>%
#JGD2011への投影変換
sf::st_transform(crs=sf::st_crs(6668)) %>%
#平面直角座標第9系への投影変換
sf::st_transform(crs=sf::st_crs(6677)) %>%
#ラインデータを融合する為の処理
dplyr::mutate(d=1) %>%
dplyr::group_by(d) %>%
dplyr::summarise() %>%
#ポリゴン化
sf::st_polygonize()
#出力結果を返す
return(dat)
}
#被害地域ポリゴンを読み込み・投影変換・ポリゴン化した結果を融合する
#lapplyを実行して返されるlistオブジェクトを行方向に結合する
higai_po <- do.call(rbind,
#3つのkmlファイルについて、上で定義した関数を実行する
lapply(c("常総地区の推定浸水範囲_201509111000.kml",
"常総地区の推定浸水範囲_201509111300.kml",
"常総地区の推定浸水範囲_201509121530.kml"),
function(x){read_polygonize(x)})) %>%
#ポリゴンデータを融合する
dplyr::group_by(d) %>%
dplyr::summarise()ステップ5：浸水地域のダミー変数の作成
ステップ3と同じ要領で、浸水地域ダミー変数を作成し、地価ポイントに付与する。ダミー変数を地価ポイントに付与し終えたら、geometry属性を削除した上で、回帰分析用の分析データを完成させる。
#被害地域に含まれる地価ポイントの抽出
chika_ibaraki_higai <- sf::st_intersection(chika_ibaraki,higai_po) %>%
#geometry属性を削除
sf::st_set_geometry(value=NULL) %>%
#被害地域ダミーd_higaiを作成
dplyr::mutate(d_higai=1) %>%
#不要な変数を削除する
dplyr::select(FID,d_higai)
#d_higaiを地価ポイントに左結合する
chika_ibaraki_final <- chika_ibaraki %>%
dplyr::left_join(chika_ibaraki_higai,by="FID") %>%
#d_higaiが結合されなかった（=地域外）地価ポイントについて、d_higaiの値を0に置換
dplyr::mutate(d_higai=ifelse(!is.na(d_higai),d_higai,0)) %>%
#geometry属性を削除する
sf::st_set_geometry(value=NULL)ステップ6：Rの設定
ステップ0で実行済みのため省略。
ステップ7：データの読み込みと抽出
上で作成したデータセットについて下処理を行う。
#分析対象となる自治体名に関するベクトルを作成
city <- c("かすみがうら",
"つくば",
"つくばみらい",
"下妻",
"五霞",
"利根",
"取手",
"古河",
"土浦",
"坂東",
"守谷",
"常総",
"桜川",
"牛久",
"石岡",
"稲敷",
"筑西",
"結城",
"茨城八千代",
"茨城境",
"茨城河内",
"阿見",
"龍ケ崎")
#上で作成したデータセットに下処理を行う
chika <- chika_ibaraki_final %>%
#見出し番号が000（住宅地）かつ上の自治体に含まれるポイントのみ抽出
dplyr::filter(L01_001=="000"&L01_022%in%city) %>%
#必要な変数のみ残す
dplyr::select(FID,L01_085,L01_086,L01_087,L01_088,L01_089,L01_090,L01_091,L01_092,L01_093,L01_094,d_shinsui,d_higai) %>%
#変数L01_085〜L01_094について、型をcharacterからnumericに変換
dplyr::mutate(across(L01_085:L01_094,as.numeric)) %>%
#変数L01_085〜L01_094の全てが正値のポイントのみ抽出
dplyr::filter(across(L01_085:L01_094,~.>0)) %>%
#扱いやすいように変数名を変更
magrittr::set_colnames(c("FID","2012","2013","2014","2015","2016","2017","2018","2019","2020","2021","d_shinsui","d_higai"))ステップ8：パネルデータの作成・ステップ9：変数の作成
下処理を行ったデータセットをwide型からlong型に変換する。
#変数2012〜2021について、wide型からlong型に変換する。その際、FID,d_shinsui,d_higaiは同じ値を繰り返す
chika_long <- tidyr::pivot_longer(chika,cols=-c(FID,d_shinsui,d_higai),names_to="year",values_to="price") %>%
#地価の自然対数を取る
dplyr::mutate(ln_price=log(price),
#2016年以降の観測値であれば1を取るダミー変数d_afterを作成
d_after=ifelse(year>=2016,1,0),
#浸水想定区域内だが被害地域外であれば1を取るダミー変数d_nshinsuiを作成
d_nshinsui=ifelse(d_shinsui==1&d_higai==0,1,0))
#データの先頭行を表示
head(chika_long)## # A tibble: 6 × 8
## FID d_shinsui d_higai year price ln_price d_after d_nshinsui
## <int> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 119 0 0 2012 32100 10.4 0 0
## 2 119 0 0 2013 30700 10.3 0 0
## 3 119 0 0 2014 30000 10.3 0 0
## 4 119 0 0 2015 29800 10.3 0 0
## 5 119 0 0 2016 29800 10.3 1 0
## 6 119 0 0 2017 29700 10.3 1 0#時点ダミーを作成
year_d <- chika_long %>%
#データセットから時点変数yearのみを取り出し
dplyr::select(year) %>%
#yearをcharacterに型変更
dplyr::mutate_all(as.character) %>%
#yearをもとにダミー変数を作成
dummy::dummy(int=T) %>%
#ダミー変数の名前を変更
magrittr::set_colnames(gsub("year_","d_",colnames(.)))
#時点ダミーをデータセットに横結合
chika_long <- cbind(chika_long,year_d)
#データの先頭行を表示
head(chika_long)## FID d_shinsui d_higai year price ln_price d_after d_nshinsui d_2012 d_2013
## 1 119 0 0 2012 32100 10.37661 0 0 1 0
## 2 119 0 0 2013 30700 10.33202 0 0 0 1
## 3 119 0 0 2014 30000 10.30895 0 0 0 0
## 4 119 0 0 2015 29800 10.30226 0 0 0 0
## 5 119 0 0 2016 29800 10.30226 1 0 0 0
## 6 119 0 0 2017 29700 10.29890 1 0 0 0
## d_2014 d_2015 d_2016 d_2017 d_2018 d_2019 d_2020 d_2021
## 1 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0
## 3 1 0 0 0 0 0 0 0
## 4 0 1 0 0 0 0 0 0
## 5 0 0 1 0 0 0 0 0
## 6 0 0 0 1 0 0 0 0ステップ10：モデルの推定
#モデル1の推定
base_fe <- fixest::feols(ln_price~d_nshinsui:d_after+d_higai:d_after+d_after|FID,data=chika_long)
#モデル2の定式化
fmla_long_fe <- as.formula(paste0("ln_price~",
paste0("d_nshinsui:d_",c(2013:2021),collapse="+"),"+",
paste0("d_higai:d_",c(2013:2021),collapse="+"),"+",
paste0("d_",c(2013:2021),collapse="+"),"|FID"))
#モデル2の推定
long_fe <- fixest::feols(fmla_long_fe,data=chika_long)
#推定結果をまとめる
fixest::etable(base_fe,long_fe,signifCode=c("***"=0.01,"**"=0.05,"*"=0.10),digits=3,digits.stats=3)## base_fe long_fe
## Dependent Var.: ln_price ln_price
##
## d_after -0.047*** (0.003)
## d_nshinsui x d_after -0.006 (0.006)
## d_after x d_higai -0.113*** (0.007)
## d_2013 -0.035*** (0.001)
## d_2014 -0.054*** (0.003)
## d_2015 -0.065*** (0.003)
## d_2016 -0.073*** (0.004)
## d_2017 -0.079*** (0.004)
## d_2018 -0.084*** (0.005)
## d_2019 -0.088*** (0.005)
## d_2020 -0.092*** (0.005)
## d_2021 -0.097*** (0.006)
## d_nshinsui x d_2013 0.002 (0.002)
## d_nshinsui x d_2014 0.003 (0.004)
## d_nshinsui x d_2015 0.002 (0.006)
## d_nshinsui x d_2016 -0.001 (0.007)
## d_nshinsui x d_2017 -0.003 (0.007)
## d_nshinsui x d_2018 -0.004 (0.008)
## d_nshinsui x d_2019 -0.005 (0.009)
## d_nshinsui x d_2020 -0.006 (0.010)
## d_nshinsui x d_2021 -0.007 (0.010)
## d_2013 x d_higai 0.004 (0.002)
## d_2014 x d_higai 0.006. (0.004)
## d_2015 x d_higai 0.006 (0.004)
## d_2016 x d_higai -0.067*** (0.004)
## d_2017 x d_higai -0.069*** (0.005)
## d_2018 x d_higai -0.110*** (0.007)
## d_2019 x d_higai -0.128*** (0.008)
## d_2020 x d_higai -0.138*** (0.009)
## d_2021 x d_higai -0.145*** (0.010)
## Fixed-Effects: ----------------- -----------------
## FID Yes Yes
## ____________________ _________________ _________________
## S.E.: Clustered by: FID by: FID
## Observations 2,480 2,480
## R2 0.998 0.998
## Within R2 0.434 0.628
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1