3-5_高速道路開業と雇用水準の変化
Keisuke Takano
【内容】
高速道路インターチェンジ（IC）開業前後での、IC周辺の従業者数の変化を観測します。具体的には、圏央道（首都圏中央連絡自動車道）のICのうち、2011年〜2015年に開業したICの周辺地域における従業者数変化を、500mメッシュ単位で定量的に計測します。
本演習ではまず、ICのポイントデータと路線網のラインデータを国土数値情報から、500mメッシュ単位の経済センサスの従業者数データを「政府統計の総合窓口」（e-Stat）の「地図で見る統計（統計GIS）」からダウンロードします。次に、ダウンロードしたデータをR上に読み込んだ上で、ICからの距離帯毎に従業者数データを抽出します。その上で、IC開業前後における従業者数増減を地図上に可視化し、その空間的傾向を把握します。また、ICに隣接した地域とそうでない地域の間での、増減傾向の比較もあわせて行います。
なお、この章では、データの読み込み、加工、可視化、集計までを、一貫してRを用いて行う方法を解説します。
【学ぶポイント】
- Rのsfパッケージとdplyrパッケージを組み合わせて、空間データを加工する方法
- Rのggplot2パッケージを用いて、空間データを可視化する方法
【R Markdownの実行方法】
このRmdファイルの中身のうち、Rコードとして実行されるのは「```」という記号で上下挟まれた部分で、以後「チャンク」と呼びます。チャンク内のコードは、実行したい行を選択した状態で、Sourceペイン右上の[Run]ボタンをクリックすれば実行できます。また、Sourceペイン左上の[Knit]ボタンをクリックすると、Rmd内のテキスト、チャンク内のコード、コードの実行結果をHTML形式で出力することができます。その出力結果が「analysis_road_main.html」です。なお、チャンク内のコードにエラーがある場合、Knitは実行途中で止まってしまうので注意しましょう。
【使用データ】
今回分析で用いるRコードは全て、「analysis_road_main.Rmd」というR
Markdown（Rmd）ファイルにまとめています。コードを実行する際には、このRmdファイルをRStudioで開いた上で実行してください。また、以下の分析データは全てRmdファイルと同じフォルダーに格納してください。
高速道路に関する空間データは、国土数値情報ダウンロードにある、「高速道路時系列データ」です。「高速道路時系列データ」は、路線網のラインデータと、ICのポイントデータの2つを含んでいます。
- N06-20_HighwaySection.shp：令和2年の全国の路線網のラインデータ。
- N06-20_Joint.shp：令和2年の全国のICを含む各種接合部（ICの他に、ジャンクションを含む）のポイントデータ。
経済センサスに関する空間データは、e-Statの地図で見る統計にある、500ｍメッシュのポリゴンデータと、統計データです。紙幅の都合上、統計データについては、ダウンロードされた2009年・2016年のデータそのままではなく、事前にクレンジングを行ったものを使用します。クレンジングの詳細な工程は、「analysis_road_app.Rmd」を参照してください。
- MESH05339.shp：圏央道開業区間をカバーする500ｍメッシュの境界データ。
- ec_2009_2016.csv：500mメッシュ単位で観測された、2009年・2016年の従業者数を収録したcsvファイル。
【演習】
ステップ1：各種データのダウンロード
注：このステップは，演習データを用いることで省略できます．
高速道路時系列データを、以下の要領でダウンロードします。国土数値情報ダウンロードを開きます。「高速道路時系列（ライン）（ポイント）」をクリックします。高速道路時系列データのダウンロードページが開いたら、下にスクロールし、「ダウンロードするデータの選択」で「全国世界測地系 令和2年」をクリックし、データをダウンロードします。
fig1
ダウンロードしたzipファイル（N06-20_GML.zip）を展開し、全てのファイルを「analysis_road_main.Rmd」を設置したフォルダに移動します（その際、ファイルをサブフォルダから取り出し、「analysis_road_main.Rmd」と同じ階層に設置してください）。
続いて、経済センサスのデータをダウンロードします。e-Statの地図で見る統計（統計GIS）を開きます。下にスクロールすると、「統計データダウンロード」と「境界データダウンロード」の2つがありますが、メッシュの境界データは「境界データダウンロード」からダウンロードします。
fig2
「境界データダウンロード」をクリックします。「境界一覧」の「4次メッシュ（500mメッシュ）」をクリックします。
fig3
「データ形式一覧」ページから、「世界測地系緯度経度・Shapefile」をクリックします。
fig4
「データダウンロード」ページには、ダウンロード可能なメッシュデータが一覧表示されます。「M」の後に続く4桁の番号は第1次地域区画の地域メッシュコードです。
fig5
「データダウンロード」ページ右上の「1次メッシュ枠情報」をクリックし、圏央道開業区間に対応するメッシュ番号を確認します（今回のケースでは、5339です）。
fig6
ページを移動して（5ページ目あたりにあります）「M5339」が見つかったら、対応する「世界測地系緯度経度・Shapefile」をクリックし、zipファイルをダウンロードします。zipを展開し、全ファイルを「analysis_road_main.Rmd」を設置したフォルダに移動します。
fig7
ダウンロードした「境界データ」は、経済センサスの各種変数を含まないので、別途「統計データ」をダウンロードする必要があります。再びe-Statの「地図で見る統計（統計GIS）」を開き、「統計データダウンロード」をクリックします。
fig8
「政府統計名」から、2009年・2016年の統計データをダウンロードします。2009年は「経済センサス－基礎調査」、2016年は「経済センサス－活動調査」からダウンロードしますが、各々手順は同じですので、ここでは2009年データについてのみ手順を示します。
fig9
「経済センサス－基礎調査」➡「4次メッシュ（500mメッシュ）」➡「全産業事業所数及び全産業従業者数」とクリックしていくと、境界データの場合と同様、メッシュコード別にダウンロード可能なデータの一覧が出ます。境界データと同じく「M5339」の「csv」ボタンをクリックし、zipファイルをダウンロードします。zipファイルを展開し、展開されたファイルを「analysis_road_main.Rmd」を設置したフォルダに移動します。
fig10
2016年についても同様の手順でデータをダウンロードします。
「analysis_road_app.Rmd」コードを実行すると、整形済みの経済センサスデータ「ec_2009_2016.csv」が作成されます。
ステップ2：Rを用いたGISデータの読み込みと抽出
パッケージの読み込み
Rで行うデータ加工・分析に必要なパッケージを，install.packages関数を用いてインストールします（既にインストール済みの場合はスキップしてください）。筆者は既に、これらパッケージをインストール済みなので、実行コードの先頭に「#」記号を付けることで，その行が実行されない（すなわち単なるコメントとして扱う）ようにしています。
sfはシェープファイルを始めとするベクタデータを用いた空間解析のためのパッケージ、dplyrは空間データを含む様々なデータをハンドリング・集計するためのパッケージ、ggplot2は可視化するためのパッケージです。
# install.packages("sf")
# install.packages("dplyr")
# install.packages("ggplot2")インストールしたパッケージを呼び出します。インストールは1回行えばよいですが、呼び出しはRを起動する度に行います。
library(sf)
library(dplyr)
library(ggplot2)シェープファイルの読み込み
分析に用いるシェープファイル一式を読み込みます。シェープファイルの読み込みには、sfのst_read関数を用います。引数dsnには読み込むシェープファイルのパス（Rmdファイルと同じフォルダにシェープファイルがある場合にはファイル名だけで可）を指定します。引数stringsAsFactorsは、シェープファイル内の変数（主に文字列）をfactor型という特殊な形式として扱うかを指定する変数ですが、今回は文字列を文字列として扱うので、FALSEを指定します。引数optionsには、シェープファイルに含まれる変数の文字コードを明示的に指定し、文字化けを防ぎます。今回読み込むシェープファイルの文字コードは、CP932という、Windows環境下で作成された多くのシェープファイルが当てはまるものです。
#高速道路の路線データ（ライン）
hw_line <- sf::st_read(dsn="N06-20_HighwaySection.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932")以下では、全ての空間データの地理座標系をJGD2011、投影座標系を平面直角座標第9系に統一した上で分析を行うので、必要に応じて投影変換を行います。路線データと接合部データでは、既に地理座標系がJGD2011に定義されているので、投影座標系への変換のみを行います。投影変換には、sfのst_transform関数を用います。引数crsには、EPSGコードを元に生成される座標系の情報を指定します。平面直角座標系第9系（JGD2011）のEPSGコードは「6677」ですが、コードから座標系の情報を呼び出す際には、sfのst_crs関数を用います。
#高速道路の路線データ（ライン）
hw_line <- sf::st_read(dsn="N06-20_HighwaySection.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932") %>%
#地理座標系（JGD2011）から投影座標系（平面直角座標第9系）へ投影変換
sf::st_transform(crs=sf::st_crs(6677))st_crs関数の実行結果な中身は以下のようになっていて、緯度経度換算の原点や座標系がカバーしているエリアに関する情報を含むことが分かります。
sf::st_crs(6677)## Coordinate Reference System:
## User input: EPSG:6677
## wkt:
## PROJCRS["JGD2011 / Japan Plane Rectangular CS IX",
## BASEGEOGCRS["JGD2011",
## DATUM["Japanese Geodetic Datum 2011",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",6668]],
## CONVERSION["Japan Plane Rectangular CS zone IX",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",36,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",139.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9999,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (X)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (Y)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
## AREA["Japan - onshore - Honshu - Tokyo-to. (Excludes offshore island areas of Tokyo-to covered by Japan Plane Rectangular Coordinate System zones XIV, XVIII and XIX)."],
## BBOX[29.31,138.4,37.98,141.11]],
## ID["EPSG",6677]]接続部ポイントと500mメッシュのシェープファイルも同様に読み込みます。接続部ポイントでは、路線データと同様、既に地理座標系が定義されていますので、平面直角座標系第9系への投影変換のみ行います。一方、500mメッシュの元の地理座標系はJGD2000になっているので、その他データと合わせるために、まずは地理座標系をJGD2011（EPSG：6668）へ投影変換し、さらに投影座標系の平面直角座標系第9系へ投影変換します。
#高速道路のICデータ（ポイント）
hw_point <- sf::st_read(dsn="N06-20_Joint.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932") %>%
sf::st_transform(crs=sf::st_crs(6677))
#500mメッシュデータ（ポリゴン）
grid <- sf::st_read(dsn="MESH05339.shp",
stringsAsFactors=FALSE,
quiet=TRUE) %>%
#地理座標系（JGD2000）から地理座標系（JGD2011）へ投影変換
sf::st_transform(crs=sf::st_crs(6668)) %>%
#地理座標系（JGD2011）から投影座標系（平面直角座標第9系）へ投影変換
sf::st_transform(crs=sf::st_crs(6677))読み込んだshapefileの先頭行を、head関数を用いて表示させます。「Projected CRS」を見て、全てのデータで座標系が揃っているかを確認しましょう。
#各データの先頭6行を表示
head(hw_line)## Simple feature collection with 6 features and 11 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -848952.7 ymin: -339532.6 xmax: -80542.18 ymax: -74364.29
## Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
## N06_001 N06_002 N06_003 N06_004 N06_005 N06_006 N06_007
## 1 2014 2014 9999 EA02_136002 <NA> <NA> 伊豆縦貫自動車道
## 2 2014 2014 9999 EA02_043021 <NA> <NA> 東九州自動車道
## 3 2014 2014 9999 EA02_310003 <NA> <NA> 高知東部自動車道
## 4 2014 2014 9999 EA02_329001 <NA> <NA> 仁摩温泉津道路
## 5 2014 2014 9999 EA02_043023 <NA> <NA> 東九州自動車道
## 6 2014 2014 9999 EA02_330001 <NA> <NA> 九州横断自動車道延岡線
## N06_008 N06_009 N06_010 N06_011 geometry
## 1 3 2 2 <NA> LINESTRING (-81120.36 -9606...
## 2 1 2 2 <NA> LINESTRING (-752626.6 -3362...
## 3 3 2 2 <NA> LINESTRING (-566842.1 -2548...
## 4 2 2 2 <NA> LINESTRING (-680464 -74364....
## 5 1 2 2 <NA> LINESTRING (-822270.8 -2091...
## 6 1 2 2 <NA> LINESTRING (-848952.7 -3235...head(hw_point)## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -809450.8 ymin: -236696.5 xmax: 95796.89 ymax: 190895.7
## Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
## N06_012 N06_013 N06_014 N06_015 N06_016 N06_017 N06_018 N06_019
## 1 2015 2015 9999 EA03_513042 <NA> <NA> 南相馬鹿島SIC 2
## 2 2015 2015 9999 EA03_524062 <NA> <NA> 高岡砺波SIC 2
## 3 2015 2015 9999 EA03_521021 <NA> <NA> 南砺SIC 2
## 4 2015 2015 9999 EA03_517015 <NA> <NA> 府中SIC 2
## 5 2015 2015 9999 EA03_543043 <NA> <NA> 豊前 1
## 6 2015 2015 9999 EA03_543044 <NA> <NA> 上毛SIC 2
## geometry
## 1 POINT (95796.89 190895.7)
## 2 POINT (-253274.3 76467.12)
## 3 POINT (-262271.8 70589.06)
## 4 POINT (-30653.41 -37664.76)
## 5 POINT (-809450.8 -233404.3)
## 6 POINT (-805557.9 -236696.5)head(grid)## Simple feature collection with 6 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -75757 ymin: -73661.37 xmax: -73480.47 ymax: -72718.04
## Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
## KEY_CODE MESH1_ID MESH2_ID MESH3_ID MESH4_ID OBJ_ID
## 1 533900001 5339 00 00 1 1
## 2 533900002 5339 00 00 2 2
## 3 533900003 5339 00 00 3 3
## 4 533900004 5339 00 00 4 4
## 5 533900011 5339 00 01 1 5
## 6 533900012 5339 00 01 2 6
## geometry
## 1 POLYGON ((-75188.81 -73647....
## 2 POLYGON ((-74620.62 -73652....
## 3 POLYGON ((-75184.95 -73185....
## 4 POLYGON ((-74616.79 -73189....
## 5 POLYGON ((-74052.43 -73656....
## 6 POLYGON ((-73484.24 -73661....分析対象となる圏央道の区間上にあるICポイントのみ抽出
接続部ポイントから、開業前後の変化を見る対象となるICのポイントのみを抽出します。まず、2009-2016年に開業した圏央道ICの名前を格納したベクトルを作成します。次に、dplyrのfilter関数を用いて、変数「N06_018」（接続部の名称）の値がベクトルの要素に一致するポイントだけを抽出します。
もちろんfilter関数は、「filter(N06_018==“寒川南”)」というように、特定のひとつの条件に一致するポイントを抽出するためにも用いることができますが、今回の例のように、ベクトルの形式で条件をまとめて渡せる点に特長があります。
#2009-2016年の間に開業した圏央道のICの内，分析対象とするICの名前をまとめた文字列ベクトル
ic <- c("寒川南","寒川北","海老名","圏央厚木","相模原愛川","相模原","高尾山")
#IC・JCTポイントから，ベクトルicに含まれるものだけを抽出する
hw_point_keno <- hw_point %>%
dplyr::filter(N06_018%in%ic)
#抽出したICポイントの中身を表示（先頭10行を表示）
head(hw_point_keno,10)## Simple feature collection with 7 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -51624.03 ymin: -71758.7 xmax: -41193 ymax: -41576.02
## Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
## N06_012 N06_013 N06_014 N06_015 N06_016 N06_017 N06_018 N06_019
## 1 2015 2015 9999 EA03_623051 <NA> <NA> 相模原 1
## 2 2013 2013 9999 EA03_623041 <NA> <NA> 寒川北 1
## 3 2010 2013 9999 EA03_520043 <NA> <NA> 海老名 1
## 4 2013 2013 9999 EA03_520044 <NA> <NA> 圏央厚木 1
## 5 2013 2013 9999 EA03_623040 <NA> <NA> 寒川南 1
## 6 2012 2012 9999 EA03_623038 <NA> <NA> 高尾山 1
## 7 2013 2013 9999 EA03_520045 <NA> <NA> 相模原愛川 1
## geometry
## 1 POINT (-48929.15 -46310.83)
## 2 POINT (-41316.61 -68181.66)
## 3 POINT (-41475.35 -62145.54)
## 4 POINT (-41790.72 -57600.47)
## 5 POINT (-41193 -71758.7)
## 6 POINT (-51624.03 -41576.02)
## 7 POINT (-43487.36 -52801.34)ステップ3：ICから0m-500mと500m-1000mの距離帯にあるメッシュの抽出
ICポイントから500m・1kmのバッファを生成
開業したICから0m-500mの距離帯にあるメッシュを開業の影響を受ける地域、500m-1000mの距離帯にあるメッシュを開業の影響を受けない地域と仮定した上で、開業前後での従業者数の変化を比較していきます。そのためにはまず、0m-500mの距離帯と500m-1000mの距離帯それぞれに含まれるメッシュを抽出する必要があります。抽出の方法はいくつかありますが、例として、バッファを用いた抽出方法を示します。
まず、開業したICから500mと1000mのバッファを生成します。バッファの生成には、sfのst_buffer関数を用います。引数distには、m単位のバッファ距離を指定します。
#ICポイントから500mのバッファを発生させる
buffer_0500m <- hw_point_keno %>%
sf::st_buffer(dist=500)
#ICポイントから1kmのバッファを発生させる
buffer_1000m <- hw_point_keno %>%
sf::st_buffer(dist=1000)生成した500mバッファの形状をプロットしてみましょう。sfパッケージ群で扱われる空間データをプロットする際には、ggplot2のgeom_sf関数を用います。引数dataには、プロットしたい空間データを指定します。
ggplot2::ggplot() +
ggplot2::geom_sf(data=buffer_0500m)
ggplot2では、2つ以上の空間データを重ねてプロットできます。その際には、geom_sf関数によるプロットのコマンドを、「+」記号（「%>%」記号を使わないよう注意）でつなげます。例えば、1000mバッファの上に500mバッファを重ねてプロットしてみます。
ggplot2::ggplot() +
ggplot2::geom_sf(data=buffer_1000m) +
ggplot2::geom_sf(data=buffer_0500m)
各バッファとメッシュの共通部分を取り出す
ここで、各バッファと500mメッシュとの共通部分を取り出します。2つの空間データ間で共通部分を取る時には、sfのst_intersection関数を用います。引数x,yには、共通部分を取りたい2つの空間データを指定します。
#500mバッファとメッシュの共通部分を取る（ポリゴンが生成される）
grid_buffer_0500m <- sf::st_intersection(x=grid,y=buffer_0500m)
#1kmバッファとメッシュの共通部分を取る（ポリゴンが生成される）
grid_buffer_1000m <- sf::st_intersection(x=grid,y=buffer_1000m)抽出された共通部分の中の、10桁の数字で表される変数KEY_CODE（元々メッシュデータ内にある変数です）が、各距離帯に含まれるメッシュを特定する際の手がかりとなるキー変数で、以下ではメッシュコードとも呼称します。head関数を用いて、データの先頭6行を表示し、KEY_CODEの形式を見てみましょう。
#共通部分ポリゴンの先頭行を表示
head(grid_buffer_0500m)## Simple feature collection with 6 features and 14 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -49429.15 ymin: -46810.83 xmax: -48429.15 ymax: -45810.83
## Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
## KEY_CODE MESH1_ID MESH2_ID MESH3_ID MESH4_ID OBJ_ID N06_012 N06_013
## 7573 533922931 5339 22 93 1 7573 2015 2015
## 7574 533922932 5339 22 93 2 7574 2015 2015
## 7575 533922933 5339 22 93 3 7575 2015 2015
## 7576 533922934 5339 22 93 4 7576 2015 2015
## 10013 533932031 5339 32 03 1 10013 2015 2015
## 10014 533932032 5339 32 03 2 10014 2015 2015
## N06_014 N06_015 N06_016 N06_017 N06_018 N06_019
## 7573 9999 EA03_623051 <NA> <NA> 相模原 1
## 7574 9999 EA03_623051 <NA> <NA> 相模原 1
## 7575 9999 EA03_623051 <NA> <NA> 相模原 1
## 7576 9999 EA03_623051 <NA> <NA> 相模原 1
## 10013 9999 EA03_623051 <NA> <NA> 相模原 1
## 10014 9999 EA03_623051 <NA> <NA> 相模原 1
## geometry
## 7573 POLYGON ((-48903.33 -46554....
## 7574 POLYGON ((-48903.33 -46554....
## 7575 POLYGON ((-49366.73 -46552....
## 7576 POLYGON ((-48903.33 -46554....
## 10013 POLYGON ((-49377.62 -46090....
## 10014 POLYGON ((-48900.8 -46092.7...head(grid_buffer_1000m)## Simple feature collection with 6 features and 14 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -49929.15 ymin: -47310.83 xmax: -48225.17 ymax: -46087.14
## Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
## KEY_CODE MESH1_ID MESH2_ID MESH3_ID MESH4_ID OBJ_ID N06_012 N06_013
## 7532 533922824 5339 22 82 4 7532 2015 2015
## 7535 533922833 5339 22 83 3 7535 2015 2015
## 7536 533922834 5339 22 83 4 7536 2015 2015
## 7539 533922843 5339 22 84 3 7539 2015 2015
## 7570 533922922 5339 22 92 2 7570 2015 2015
## 7572 533922924 5339 22 92 4 7572 2015 2015
## N06_014 N06_015 N06_016 N06_017 N06_018 N06_019
## 7532 9999 EA03_623051 <NA> <NA> 相模原 1
## 7535 9999 EA03_623051 <NA> <NA> 相模原 1
## 7536 9999 EA03_623051 <NA> <NA> 相模原 1
## 7539 9999 EA03_623051 <NA> <NA> 相模原 1
## 7570 9999 EA03_623051 <NA> <NA> 相模原 1
## 7572 9999 EA03_623051 <NA> <NA> 相模原 1
## geometry
## 7532 POLYGON ((-49472.35 -47014....
## 7535 POLYGON ((-49472.35 -47014....
## 7536 POLYGON ((-48905.86 -47017....
## 7539 POLYGON ((-48339.38 -47020....
## 7570 POLYGON ((-49640.79 -47013....
## 7572 POLYGON ((-49899.92 -46549....ICポイントから0m-500m圏内にあるメッシュ、500m-1km圏内にあるメッシュそれぞれのコードを特定する
各距離帯の共通部分ポリゴンから、変数KEY_CODEの値を取り出しベクトル化します。これにより、ICから500m圏内にあるメッシュコードと、1000m圏内にあるメッシュコードに関するリストを作成できます。
#500mバッファとメッシュの共通部分ポリゴンから変数KEY_CODEを選択し，ベクトル化する
grid_id_0000_0500 <- grid_buffer_0500m$KEY_CODE
#1kmバッファとメッシュの共通部分ポリゴンから変数KEY_CODEを選択し，ベクトル化する
grid_id_0000_1000 <- grid_buffer_1000m$KEY_CODE2つのメッシュコードベクトルを用いて、500m-1000mの距離帯にあるメッシュコードを特定します。そのためには、1000m圏内メッシュコードから、500m圏内メッシュコードを除いたIDベクトルを作成する必要があります。その操作は、dplyrのsetdiff関数を用いて行えます。下では、1000m圏内メッシュコードベクトルには含まれるが、500m圏内メッシュコードベクトルには含まれない要素を抽出しています。
#1kmバッファ内のKEY_CODEベクトルと，500mバッファ内のKEY_CODEベクトルの差分を取る
#これにより，ICポイントから500m-1km圏内に位置するメッシュ（そのKEY_CODE）が取得できる
grid_id_0500_1000 <- dplyr::setdiff(grid_id_0000_1000,grid_id_0000_0500)特定したコードに対応するメッシュを抽出する
生成したメッシュコードベクトルを用いて、0m-500mの距離帯と500m-1000mの距離帯に含まれるメッシュを抽出します。ここでは、filter関数を用いて、ベクトルに基づく抽出を行います。
#500m圏内に含まれるメッシュを抽出する
grid_0000_0500 <- grid %>%
dplyr::filter(KEY_CODE%in%grid_id_0000_0500)
#500m-1km圏内に含まれるメッシュを抽出する
grid_0500_1000 <- grid %>%
dplyr::filter(KEY_CODE%in%grid_id_0500_1000)ステップ4：地物のプロット
メッシュ・路線・ICを1枚の地図にプロットする
抽出されたメッシュを、ICポイントや路線網ラインと重ねてプロットする方法を、段階を追って示します。ここでは先ほどとは異なり、ポリゴンの色や枠線に条件を指定した上でプロットします。geom_sfで、色は引数fillで指定できます。下の例では、0m-500m圏内にあるメッシュを赤、500m-1000m圏内にあるメッシュを青でプロットしています。枠線無しでプロットしたい場合には、引数ltyを0に指定します。
#2つのオブジェクトを重ねてプロット
ggplot2::ggplot() +
#抽出されたICから0m-500m圏内にあるメッシュを，赤＆枠線無でプロット
ggplot2::geom_sf(data=grid_0000_0500,fill="red",lty=0) +
#抽出されたICから500m-1km圏内にあるメッシュを，青＆枠線無でプロット
ggplot2::geom_sf(data=grid_0500_1000,fill="blue",lty=0)
メッシュとICポイント・路線網ラインを重ねてプロットします。ポイントやラインの場合、色の指定を行う引数は、fillではなくcolorです。例えば、路線を緑で、ICを黒でプロットすると以下のようになります。
#全ての地物をプロットする為の設定
ggplot2::ggplot() +
#抽出されたICから0m-500m圏内にあるメッシュを，赤＆枠線無でプロット
ggplot2::geom_sf(data=grid_0000_0500,fill="red",lty=0) +
#抽出されたICから500m-1km圏内にあるメッシュを，青＆枠線無でプロット
ggplot2::geom_sf(data=grid_0500_1000,fill="blue",lty=0) +
#路線を緑でプロット
ggplot2::geom_sf(data=hw_line,color="green") +
#ICを黒でプロット
ggplot2::geom_sf(data=hw_point_keno,color="black")
ICについては、2009-2016年に開業した圏央道ICだけを抽出済みですが、路線については抽出処理をしていないため、全国の路線網がそのままプロットされています。
ここでは例えば、プロットの描画領域を、抽出されたメッシュの空間範囲だけに限定する方法を示します。空間範囲の取得には、sfのst_bbox関数を用います。例えば、抽出されたメッシュの空間範囲は以下のようになります。x（東西）,y（南北）軸の各々について、最大値・最小値が格納されたベクトルが取得できます。
#coord_sfの行で使った描画領域の中身
sf::st_bbox(obj=grid_0500_1000)## xmin ymin xmax ymax
## -52843.85 -72948.03 -39908.08 -40522.96ちなみに、このベクトルからx軸方向に関する要素だけを取り出したい場合は以下のように書きます。
sf::st_bbox(obj=grid_0500_1000)[c("xmin","xmax")]## xmin xmax
## -52843.85 -39908.08st_bbox関数を用いてプロットの空間範囲を限定する際には、下のように、ggplot2のcoord_sf関数を用います。引数xlimには、最大値・最小値ベクトルの要素のうちx軸方向に関するもの、ylimにはy軸方向に関するものを指定します。
#全ての地物をプロットする為の設定
ggplot2::ggplot() +
#抽出されたICから0m-500m圏内にあるメッシュを，赤＆枠線無でプロット
ggplot2::geom_sf(data=grid_0000_0500,fill="red",lty=0) +
#抽出されたICから500m-1km圏内にあるメッシュを，青＆枠線無でプロット
ggplot2::geom_sf(data=grid_0500_1000,fill="blue",lty=0) +
#路線を緑でプロット
ggplot2::geom_sf(data=hw_line,color="green") +
#ICを黒でプロット
ggplot2::geom_sf(data=hw_point_keno,color="black") +
#500m-1km圏内にあるメッシュが全てプロットできるように描画領域を指定
ggplot2::coord_sf(xlim=sf::st_bbox(grid_0500_1000)[c("xmin","xmax")],ylim=sf::st_bbox(grid_0500_1000)[c("ymin","ymax")])
仕上げに、ICの名前を重ねてプロットしてみます。空間データの属性をラベルとして地図上にプロットする際には、ggplot2のgeom_sf_text関数を使います。引数dataには、ラベルにしたい変数が入っている空間データ、aes()で囲まれた引数labelには、ラベルにしたい変数を指定します。引数familyに「HiraKakuPro-W3」を指定していますが、これは日本語を文字化けなくプロットする為のものです。文字のサイズは引数sizeで、文字色はcolorで指定できます。以降、フォントのWarningが出る場合がありますが、無視して先に進んで構いません。
#全ての地物をプロットする為の設定
ggplot2::ggplot() +
#抽出されたICから0m-500m圏内にあるメッシュを，赤＆枠線無でプロット
ggplot2::geom_sf(data=grid_0000_0500,fill="red",lty=0) +
#抽出されたICから500m-1km圏内にあるメッシュを，青＆枠線無でプロット
ggplot2::geom_sf(data=grid_0500_1000,fill="blue",lty=0) +
#路線を緑でプロット
ggplot2::geom_sf(data=hw_line,color="green") +
#ICを黒でプロット
ggplot2::geom_sf(data=hw_point_keno,color="black") +
#500m-1km圏内にあるメッシュが全てプロットできるように描画領域を指定
ggplot2::coord_sf(xlim=sf::st_bbox(grid_0500_1000)[c("xmin","xmax")],ylim=sf::st_bbox(grid_0500_1000)[c("ymin","ymax")]) +
#IC名（変数N06_018）を白でプロット
#Windowsで実行の際は引数familyの指定は不要）
ggplot2::geom_sf_text(data=hw_point_keno,aes(label=N06_018),family="HiraKakuPro-W3",size=3,color="white")
ステップ5：経済センサスに関する前処理
2009・2016年経済センサスのcsvファイルを読み込み
整形済みの2009・2016年経済センサスのcsvファイルを、read.csv関数を用いて読み込みます。csvの文字コードはCP932ですので、引数fileEncodingでそのことを明示します。
ec_2009_2016 <- read.csv("ec_2009_2016.csv",fileEncoding="CP932")従業者数（対数値）の変化率を計算
従業者数について対数をとった上で、変化率を計算する方法を、順を追って示します。従業員数の対数を取る理由は、その分布の形にあります。hist関数を用いて、従業者数のヒストグラムをプロットすると、分布がかなり左に偏っていることが分かります。
#プロット時の日本語の文字化け防止のための設定
par(family="HiraKakuPro-W3")
#2009年の従業者数のヒストグラムをプロット
hist(ec_2009_2016$従業者数,breaks=100)
分布の偏りを減らすために、従業者数の対数を取ります。データに新たな変数を追加する際には、dplyrのmutate関数を用います。下のコードでは2つの新たな変数を追加しています。1つは、従業者数に1を足した上で対数を取った変数ln_emp（0の対数は定義されないため、1を足しています）、もう1つは、変数KEY_CODEを数値型から文字列に変換して元の変数を上書きしたものです。この型変換は、メッシュデータの変数KEY_CODEが文字列型になっており、それに合わせるために行います。
#各時点の従業者数に1を足し，対数を取った値を新たに変数として加える
ec_2009_2016 <- ec_2009_2016 %>%
dplyr::mutate(ln_emp=log(従業者数+1)) %>%
#メッシュとの結合の為にKEY_CODEを文字列型に変換
dplyr::mutate(KEY_CODE=as.character(KEY_CODE))
head(ec_2009_2016)## KEY_CODE year 従業者数 ln_emp
## 1 533900043 2009 6 1.945910
## 2 533900043 2016 7 2.079442
## 3 533900051 2009 10 2.397895
## 4 533900051 2016 9 2.302585
## 5 533900053 2009 27 3.332205
## 6 533900053 2016 5 1.791759対数変換後の従業者数ln_empについてヒストグラムを作成すると、未だ0付近に値が集中しているものの、全体としての分布の偏りは幾分減ったことが分かります。
#従業者数（対数）のヒストグラムをプロット
hist(ec_2009_2016$ln_emp,breaks=100)
各KEY_CODE（すなわちメッシュ）について、従業者数（対数）の変化率を計算します。その為にまず、dplyrのgroup_by関数を用いてデータを変数KEY_CODEによってグループ化した上で、dplyrのlag関数を用いて、変数ln_emp_lagを時間方向に1行ずらした（「1時点のラグを取る」とも言います）変数ln_emp_lagを追加します。
#従業者数（対数）について変化率を計算
ec_2009_2016_rc <- ec_2009_2016 %>%
dplyr::group_by(KEY_CODE) %>%
dplyr::mutate(ln_emp_lag=dplyr::lag(ln_emp))
head(ec_2009_2016_rc)## # A tibble: 6 × 5
## # Groups: KEY_CODE [3]
## KEY_CODE year 従業者数 ln_emp ln_emp_lag
## <chr> <int> <int> <dbl> <dbl>
## 1 533900043 2009 6 1.95 NA
## 2 533900043 2016 7 2.08 1.95
## 3 533900051 2009 10 2.40 NA
## 4 533900051 2016 9 2.30 2.40
## 5 533900053 2009 27 3.33 NA
## 6 533900053 2016 5 1.79 3.33作成したラグ変数を用いて、従業者数（対数）変化率の変数ln_emp_rcを計算します。上の処理とひとまとめにすると、以下のようになります。
#従業者数（対数）について変化率を計算
ec_2009_2016_rc <- ec_2009_2016 %>%
dplyr::group_by(KEY_CODE) %>%
dplyr::mutate(ln_emp_lag=dplyr::lag(ln_emp)) %>%
dplyr::mutate(ln_emp_rc=(ln_emp-ln_emp_lag)/ln_emp_lag)
head(ec_2009_2016_rc)## # A tibble: 6 × 6
## # Groups: KEY_CODE [3]
## KEY_CODE year 従業者数 ln_emp ln_emp_lag ln_emp_rc
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 533900043 2009 6 1.95 NA NA
## 2 533900043 2016 7 2.08 1.95 0.0686
## 3 533900051 2009 10 2.40 NA NA
## 4 533900051 2016 9 2.30 2.40 -0.0397
## 5 533900053 2009 27 3.33 NA NA
## 6 533900053 2016 5 1.79 3.33 -0.462is.finite関数とfilter関数を使って、2009年時点の従業者数（対数）が0のために変化率が無限大になったレコードを除きます。is.finite関数は、入力された値が無限大ならばTRUEを返す関数ですが、それに否定を表す演算子「!」を付けることにより、無限大「でない」場合にTRUEを返す関数として使えます。下では、変数ln_emp_rcが無限大「でない」レコードだけを抽出しています。また、変数ln_emp_rcがNAでないレコードをのみを更に抽出します。
#従業者数（対数）について変化率を計算
ec_2009_2016_rc <- ec_2009_2016_rc %>%
#変化率が計算不能だった（2009年時点の従業者数が0）のレコードを除く
dplyr::filter(!is.infinite(ln_emp_rc)) %>%
#変化率が入った行だけ残す
dplyr::filter(!is.na(ln_emp_rc))ステップ6：従業者数（対数値）の変化率をプロット
プロットの前準備
上で作成した、ICから1000m圏内に含まれるメッシュコードベクトルを用いて、ICから1000m圏内のメッシュを抽出します。次に、抽出したメッシュに変化率の計算結果を結合します。データの結合にはdplyrのleft_join関数を用います。引数byには、結合の基準となるキー変数の名前を指定します。最後に、変化率が計算されなかったメッシュを、filter関数を用いて除きます。
#1km圏内にある全てのメッシュを抽出
grid_0000_1000 <- grid %>%
dplyr::filter(KEY_CODE%in%grid_id_0000_1000)
#抽出したメッシュに，変化率の計算結果を結合
grid_rc <- grid_0000_1000 %>%
dplyr::left_join(ec_2009_2016_rc,by="KEY_CODE") %>%
#変化率が計算されなかったメッシュを除く
dplyr::filter(!is.na(ln_emp_rc))変化率のプロット
500mメッシュを、変化率の値に応じて色分けしてプロットする方法を、段階を踏んで示します。まず、色分けに関する指定を特にせずに、変化率の変数ln_emp_rcで色分けを行った場合を示します。先ほどと同様、geom_sf関数を使ってメッシュをプロットしますが、空間データ内のある特定の変数に基づいて色分けを行いたい場合には、aes()で囲まれた引数fillに、色分けの基準となる変数を指定します。
作図すると以下のようになりますが、正負両方の値に同系統の色が使われている等、プロットとしては実用に耐えるものとは言えません。
ggplot2::ggplot() +
#変化率が結合された500mメッシュをプロット
ggplot2::geom_sf(data=grid_rc,aes(fill=ln_emp_rc),lty=0)
作図して意味のあるプロットを作るためには幾分工夫が必要です。ここでは、「0を境に負の値を寒色（青系）、正の値を暖色（赤系）で塗る」、「色の区切りを自分で指定する」という2つの条件の下で変化率をプロットする方法を示します。そのためにまず、引数fillに色を区切る境界値を要素として持つベクトルを指定します。その際に用いるのがcut関数で、引数xにプロットしたい変数、breaksに境界値ベクトルを指定します。
ここでは例として、「ln_emp_rcの最小値,-0.5,0,0.5,1.0,ln_emp_rcの最大値」というように、0.5刻みで階級を区切ったプロットを作成する指定をします。
その上で、階級の色分けに使うパレットをggplot2のscale_fill_brewer関数で指定します。下のコードで指定している「Spectral」というパレットは元々、値が大きくなるにつれて暖色→寒色と色が変わるパレットですが、引数directionを-1に指定することで、寒色→暖色と順序を変えることができます。
ggplot2::ggplot() +
#変化率が結合された500mメッシュをプロット（変化率の色分けの階級は，最小値から最大値までを0.5刻み）
ggplot2::geom_sf(data=grid_rc,
aes(fill=cut(x=ln_emp_rc,breaks=c(min(ln_emp_rc),-0.5,0,0.5,1.0,max(ln_emp_rc)))),
lty=0) +
#階級の色分けには「Spectral」というパレットを使う
ggplot2::scale_fill_brewer(palette="Spectral",direction=-1)
これで、色分けの問題は解決できましたが、今度は凡例名に色分け条件がそのまま使われてしまいます。凡例名を自分で変えたい時は、ggplot2のguides関数とguide_legend関数を組み合わせて用います。guide_legend関数の引数titleに、使いたい凡例名を指定します。
ggplot2::ggplot() +
#変化率が結合された500mメッシュをプロット（変化率の色分けの階級は，最小値から最大値までを0.5刻み）
ggplot2::geom_sf(data=grid_rc,
aes(fill=cut(ln_emp_rc,breaks=c(min(ln_emp_rc),-0.5,0,0.5,1.0,max(ln_emp_rc)))),
lty=0) +
#階級の色分けには「Spectral」というパレットを使う
ggplot2::scale_fill_brewer(palette="Spectral",direction=-1) +
#図の凡例のタイトルを手動で指定する
ggplot2::guides(fill=ggplot2::guide_legend(title="RC in employment"))
ここで、各メッシュが0m-500mの距離帯の中か外かを枠線で区別できるようにします。そのためにまず、0m-500mの距離帯にあるメッシュ群を1つのポリゴンに融合します。回りくどい方法ですが、
- 500m圏内に含まれる全てのメッシュで同じ値を取る変数treatを作成（変数名は何でも構いません）
- 変数treatでメッシュをグループ化
- dplyrのsummarise関数を用いてまとめる
という手順を踏むことでまとめられます。
#プロット用に，500m圏内に含まれるメッシュをひとまとめにする
grid_0000_0500_geom <- grid_0000_0500 %>%
#変数をひとまとめにする為，500m圏内に含まれる全てのメッシュで同じ値を取る変数treatを作成
dplyr::mutate(treat=1) %>%
#treatでメッシュをグループ化
dplyr::group_by(treat) %>%
#グループをひとまとめにする
dplyr::summarise()上で作成したポリゴンを、色分けのプロットの上から塗りつぶしなしで重ねてプロットすることで、0m-500mの距離帯の内側と外側で色の傾向を比較することができます。先ほどのプロットのコードの最後に、塗りつぶしなしのポリゴンをプロットするコードを追加します。塗りつぶしを透明にするためには、引数fillに「transparent」を指定します。
ggplot2::ggplot() +
#変化率が結合された500mメッシュをプロット（変化率の色分けの階級は，最小値から最大値までを0.5刻み）
ggplot2::geom_sf(data=grid_rc,
aes(fill=cut(ln_emp_rc,breaks=c(min(ln_emp_rc),-0.5,0,0.5,1.0,max(ln_emp_rc)))),
lty=0) +
#階級の色分けには「Spectral」というパレットを使う
ggplot2::scale_fill_brewer(palette="Spectral",direction=-1) +
#図の凡例のタイトルを手動で指定する
ggplot2::guides(fill=guide_legend(title="RC in employment")) +
#500m圏内メッシュ（をひとまとめにしたもの）を，塗りつぶし無しでプロット
ggplot2::geom_sf(data=grid_0000_0500_geom,size=0.5,fill="transparent")
先ほどと同様、プロットの仕上げとして、ICポイントと路線網ラインを重ねてプロットします。あまりはっきりとは見て取れませんが、枠線の内側では黄色やオレンジのような暖色のメッシュが多い印象を受けます。
ggplot2::ggplot() +
#変化率が結合された500mメッシュをプロット（変化率の色分けの階級は，最小値から最大値までを0.5刻み）
ggplot2::geom_sf(data=grid_rc,
aes(fill=cut(ln_emp_rc,breaks=c(min(ln_emp_rc),-0.5,0,0.5,1.0,max(ln_emp_rc)))),
lty=0) +
#階級の色分けには「Spectral」というパレットを使う
ggplot2::scale_fill_brewer(palette="Spectral",direction=-1) +
#図の凡例のタイトルを手動で指定する
ggplot2::guides(fill=guide_legend(title="RC in employment")) +
#500m圏内メッシュ（をひとまとめにしたもの）を，塗りつぶし無しでプロット
ggplot2::geom_sf(data=grid_0000_0500_geom,size=0.5,fill="transparent") +
#路線ラインを緑でプロット
ggplot2::geom_sf(data=hw_line,color="green") +
#ICポイントを黒でプロット
ggplot2::geom_sf(data=hw_point_keno,color="black",size=1) +
#描画領域を，変化率が結合された500mメッシュがカバーしている範囲に絞る
ggplot2::coord_sf(xlim=sf::st_bbox(grid_rc)[c("xmin","xmax")],ylim=sf::st_bbox(grid_rc)[c("ymin","ymax")])
ステップ7：従業員数変化の平均的傾向の比較
経済センサスのデータに500m圏内フラグと2016年フラグを変数として追加
従業者数（対数）変化率プロットから、ICから500m圏内のメッシュでは従業者数が増えることが多い傾向が見て取れましたが、このことを別の方法を用いて確かめることも可能です。まず、読み込んだ経済センサスのデータに、
- after：観測年の変数yearが2016なら1、2009なら0
- treat：0m-500mの距離帯にメッシュが含まれれば1、500m-1000mの距離帯にメッシュが含まれれば0、それ以外ならNA
という値をとるフラグ変数を追加します。
条件式に基づいて変数を作成する方法は主に2通りあります。1つはifelse関数を用いた方法で、変数afterを追加する際に用いています。もう1つはdplyrのcase_when関数を用いた方法で、変数treatを追加する際に用いています。case_when関数は、
case_when(
[条件式1]が真 ~ [値A],
[条件式2]が真 ~ [値B],
[条件式3]が真 ~ [値C],
…
)
というように、条件分岐が3つ以上ある場合に、可読性が高いコードを書くことができます。
#500m圏内外を表すフラグ変数treatを作成
ec_2009_2016 <- ec_2009_2016 %>%
#時点変数yearが2016なら1，2009なら0となる変数afterを作成
dplyr::mutate(after=ifelse(year==2016,1,0)) %>%
dplyr::mutate(treat=dplyr::case_when(
#そのKEY_CODEが500m圏内のメッシュのものであれば1
KEY_CODE%in%grid_id_0000_0500 ~ 1,
#そのKEY_CODEが500m-1km圏内のメッシュのものであれば0
KEY_CODE%in%grid_id_0500_1000 ~ 0
)) %>%
#変数treatがNAになる（ICから1km以上離れている）レコードを除く
dplyr::filter(!is.na(treat))500m圏内・圏外で従業者数を前後比較
先ほど作成したフラグ変数の組み合わせ毎に、従業者数（対数）ln_empの平均値を計算します。本事例では、変数treatが0か1か（500m圏内のメッシュか否か）と変数afterが0か1か（IC開業後か否か）に関する2×2=4通りで平均値を計算します。このような4つのグループ間での平均値の比較は、3-2で扱われたDID分析と概念的に同じものです。
最初に、group_by関数を用いて、変数treat・afterを基準にデータをグループ化します。次に、summarise関数を用いて、グループ毎に変数ln_empの平均値を計算します。
集計結果を見ると、500m圏内のメッシュ（treat=1）では開通前後でln_empが上昇した一方で、500m圏外のメッシュ（treat=0）では下落した傾向が見て取れます。
#従業者数（対数）をクロス集計する
ec_2009_2016_sum <- ec_2009_2016 %>%
#変数treat（500m圏内のメッシュか）とafter（IC開業後か）についてデータをグループ化する
dplyr::group_by(treat,after) %>%
#各グループ（2×2=4グループ）について従業者数（対数）の平均値をクロス集計する
dplyr::summarise(avg=mean(ln_emp))
head(ec_2009_2016_sum)## # A tibble: 4 × 3
## # Groups: treat [2]
## treat after avg
## <dbl> <dbl> <dbl>
## 1 0 0 5.37
## 2 0 1 5.16
## 3 1 0 5.07
## 4 1 1 5.15上の集計結果から、500m圏内メッシュでの、500m圏外メッシュに比した「相対的な」ln_empの変化を計算することができます。具体的には、500m圏内・圏外のメッシュ群それぞれについて、開通前後でのln_empの増減幅を計算し、その差を取ることで計算されます。
#500m圏内・圏外で従業者数を前後比較
(5.148044-5.074181)-(5.155492-5.374023)## [1] 0.292394プロットの結果と同様、500m圏内では相対的に従業者数（対数）が増加したことが分かります。
【発展】DID分析で前後比較
紙幅の関係で書籍の本編からは割愛しましたが、500m圏内での相対的な従業者数（対数）の増加が、統計的に意味があるかを調べる為の仮説検定を行うことができます。上記のような500m圏内外での前後比較の手続きは、実は3-2でも用いられた差分の差分（difference-in-differences;
DID）分析と同じものです。
今回のDID分析の回帰モデルを用いた最も素朴な定式化は、以下のようになります。\(\ln emp_{it}\)は時点\(t\)におけるメッシュ\(i\)での従業者数（対数）、\(treat_{i}\)はメッシュ\(i\)がICから500m圏内にあれば1を取るフラグ変数（ダミー変数とも呼びます）\(after_{t}\)はIC開業後であれば1を取るフラグ変数、\(\epsilon_{it}\)は誤差項です。500m圏内での相対的な従業者数（対数）の増減の程度は、\(treat_{i} \times after_{t}\)の回帰係数\(\beta\)で捉えられます。
\[ \ln emp_{it} = \beta_{0} + treat_{i} \beta_{1} + after_{t} \beta_{2} + treat_{i} \times after_{t} \beta_{3} + \epsilon_{it} \tag{1} \]
この定式化での\(treat_{i}\)は、ICから500m圏内のメッシュ全てに関する、不時変要因の平均的傾向を捉える変数です。その代わりに、メッシュレベルの個別効果をモデルに導入することによって、メッシュ毎の不時変要因の平均的傾向を捉えることができます。\(after_{t}\)も同様に、開業後の全ての時点・全てのメッシュに共通する時変要因の平均的傾向を捉える変数ですが、その代わりに時点毎に変動する効果を導入することで、開業前後の各時点で異なる時変要因の平均的傾向を捉えることができます。
これら個別効果・時点効果をモデルに導入した回帰モデルは、二方向の「固定効果」を制御するという意味で、two-ways
fixed effect
（TW-FE）の回帰分析と呼ばれます。TW-FEとDID分析を組み合わせることで、より精緻な前後比較が可能となります。個別効果を\(\kappa_{i}\)、時点効果を\(\rho_{t}\)としてモデルに導入した定式化は、以下のようになります。
\[ \ln emp_{it} = \kappa_{i} + \rho_{t} + treat_{i} \times after_{t} \beta + \epsilon_{it} \tag{2} \]
上のような、TW-FEの回帰モデルをRで推定する方法はいくつかありますが、ここでは例えば、fixestパッケージを用いた推定法を説明します。まずパッケージをインストールした上で、Rに読み込みます。
#fixest：固定効果推定の為のパッケージ
# install.packages("fixest")
library(fixest) DID分析を行う際には「政策介入（今回の文脈ではIC開業）の前後で、アウトカム（今回の文脈では従業者数）のトレンドに差がない」という条件（平行トレンド）が満たされる必要があります。今回は素朴な方法として、IC開業前における、500m圏内外の従業者数の平均値についての差の検定を行います。本来ならば、介入前の時点を更に遡った上で、「トレンドを」比較する必要がありますが、現状利用可能な最古のデータは2009年のものなので、次善の策として差の検定を行う点に注意してください。
dplyrのfilter関数を用いて、500m圏内メッシュと500m-1km圏内メッシュの各々について、2009年時点の観測値を抽出します。
#500m圏内メッシュの2009年時点のレコードを抽出
ec_2009_0000_0500 <- ec_2009_2016 %>%
dplyr::filter(after==0&treat==1)
#500m-1km圏内メッシュの2009年時点のレコードを抽出
ec_2009_0500_1000 <- ec_2009_2016 %>%
dplyr::filter(after==0&treat==0)抽出されたレコードの内、従業員数（対数）の列を各々ベクトルとして取り出し、t検定を行います。帰無仮説を「500m圏内外で従業者数の平均値に差がない」とした検定の結果、t値は-1.0086となりました。故に、500m圏内外では従業者数の平均値に差があるとは必ずしも言えない、という結果が得られました。
#500m圏内外で、2009年時点で従業者数（対数）に差があるかを検定
t.test(ec_2009_0000_0500$ln_emp,ec_2009_0500_1000$ln_emp,var.equal=F)##
## Welch Two Sample t-test
##
## data: ec_2009_0000_0500$ln_emp and ec_2009_0500_1000$ln_emp
## t = -1.0086, df = 103.63, p-value = 0.3155
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.8893902 0.2897050
## sample estimates:
## mean of x mean of y
## 5.074181 5.374023 (2)式の定式化に基づき、500m圏内での相対的な従業者数（対数）の増加があったかを検証します。線形回帰の定式化に基づき、固定効果モデルを推定する際には、fixestのfeols関数を用います。1つ目の引数でモデルの定式化を指定しますが、下のコードでの引数の意味は、「被説明変数をln_emp、説明変数をtreatとafterの交差項とした上で、KEY_CODEに基づく個別効果とyearに基づく時間効果を更に導入する」です。2番目の引数にはデータを指定します。
推定結果を表示する際には、fixestのetable関数を用います。引数signifCodeでは回帰係数の有意水準とそれに対応する記号を、引数digitsでは回帰係数・標準誤差の有効数字を、引数digitsではその他の統計量（例えば決定係数）の有効数字を指定します。回帰係数\(\beta\)の推定値の符号は正であり、10%水準で統計的に有意であることが分かります。故に、開業したICから500m圏内のメッシュでは、相対的な従業者数（対数）の増加があったことが示唆されます。
#DID推定で、前後比較で計算された差の統計的有意性を検証
res_fe <- fixest::feols(ln_emp~treat:after|KEY_CODE+year,ec_2009_2016)
#推定結果をまとめる
fixest::etable(res_fe,signifCode=c("***"=0.01,"**"=0.05,"*"=0.10),digits=3,digits.stats=3)## res_fe
## Dependent Var.: ln_emp
##
## treat x after 0.292. (0.167)
## Fixed-Effects: --------------
## KEY_CODE Yes
## year Yes
## _______________ ______________
## S.E.: Clustered by: KEY_CODE
## Observations 246
## R2 0.917
## Within R2 0.022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ここで、今回の分析では対処できていない幾つかの問題について注記しておきます。
- 今回の分析では直線距離に基づきICへのアクセス性を数値化しました。一方で、IC周辺の道路ネットワークや、走行速度に関するデータが利用可能な場合、道路距離や所要時間など、より精緻なアクセス性の変数に基づく分析を行うことが望ましいです。
- DID分析の前提条件のひとつである、「介入（IC開業）の影響は、介入を受けなかった集団（500m圏外のメッシュ）には及ばない」というStable
Unit Treatment Value
Assumption（SUTVA）が満たされているかは、今回の分析デザインでは疑わしいです。IC開業によって元々500m圏外にあった事業所が500m圏内に移転してきたり、ICの開業効果が500mという範囲を超えてスピルオーバーしたりする場合、SUTVAは破れてしまう為、それに対処した分析（例えば、比較対象のメッシュを変更する等）を追加的に行う必要が出てきます。
- 例えば、地価公示や物流拠点数に着目して、圏央道開業前後における地域経済の変化を分析した高速道路調査会による2021年のレポート「高速道路を利用した物流に関する調査研究」では、開業前後で最寄りICまでの所要時間が減少した地域を、開業の影響を受けた地域として定義しています。
- 例えば、地価公示や物流拠点数に着目して、圏央道開業前後における地域経済の変化を分析した高速道路調査会による2021年のレポート「高速道路を利用した物流に関する調査研究」では、開業前後で最寄りICまでの所要時間が減少した地域を、開業の影響を受けた地域として定義しています。
- 今回は500mという境界値を決め打ちで設定しましたが、IC開業の効果は500mよりも広い（もしくは狭い）範囲で生じる可能性があります。その為、厳密に境界値を設定する上では、追加的な情報・データを踏まえた判断が必要となります。
- 今回のTW-FEのDID分析では、メッシュ毎に異なる時間可変の地域要因については制御を行っていません。故に、それら地域要因に関する変数をコントロール変数として回帰モデルに追加しても結果が変わらないか、更に分析を行う必要があります。
- ICの場所や高速道路のルートが、各メッシュが持つ経済的属性（例えば、潜在的な開発需要）に応じて決定されている可能性があります。その場合、処置変数\(treat_{i} \times after_{t}\)と誤差項\(\epsilon_{it}\)の間に相関が生じてしまい、推定される介入効果が過大（もしくは過小）になります。これについても、例えば操作変数法と組み合わせたDID分析を行う等して、追加的な対処が必要です。