Rによる市場アクセス分析
Keisuke Takano
はじめに
このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．予めご容赦ください．
また，いくつかのブラウザで，本文中の数式の表示が正しく行えない問題が起きています．その際には，数式のいずれかを右クリックし，「Math
Setting -> Math Renderer -> Common
HTML」という風に設定を行ってください．
目的
このRPubsでは，Rを用いた空間分析によって市場アクセス指標を作成し，回帰分析の説明変数として用いる方法を学びます．なお，この演習の内容は，An Introduction to Geographical and Urban Economicsの3・6章を主に参考にして書かれています．
お題
背景
都市経済学においては，経済活動の地理的な集中が何をもたらすか，ということが主要な関心事のひとつです．ヒト・モノの集中は，都市における集積の利益を生み出すことが知られており，その源泉としては例えば，（1）情報・知識の共有，（2）労働市場のプール，（3）中間財の共有が考えられています（Marshall,
1891）．これらを源泉とし，（1）シェアリング（大規模で不可分な施設は大都市にのみ立地），（2）マッチング（労働者・企業間での良質なマッチングを促進），（3）ラーニング（face-to-faceコンタクト）といったメカニズムを通じて，集積の利益が生み出されます（Duranton
& Puga, 2004）．
こうした集積の利益を計測する最もオーソドックスな方法が，集積を表現する何かしらの変数と，各地域の賃金（あるいは1人あたり所得）の間の関連性を見るというものです．上で記したようなメカニズムが働いているとすれば，集積と賃金の間には正の関係があります．集積を捉えるための変数として広く用いられている変数のひとつに人口密度があります．人口密度は簡便かつ分かりやすい集積の指標である一方，捉えられるのはその地域自身の経済活動のみであるという問題があります．地域間のヒト・モノの移動が盛んな現代社会では，ひとつの地域がそれ単独で経済活動を行う場面は極めて限られていると言えるでしょう．故に，隣接地域との相互作用を明示的に取り入れた変数の方が，その地域の経済活動をより正確に捉えている可能性があります．
やること
そうした問題意識の下，この演習では，人口密度に代わる変数として市場アクセス指標（market access index; MA）を用いた分析を行います．具体的には，地域\(j=1,...,n\)の時期\(t\)における経済アウトカムを\(X_{jt}\)，地域\(i\)と\(j\)との距離を\(d_{ij}\)とすると，最もシンプルな地域\(i\)のMAは以下のように定式化されます．
\[
\begin{equation}
MA_{it}=\sum_{j} \frac{X_{jt}}{d_{ij}}.
\end{equation}
\] 例えば自地域内距離\(d_{ii}=1\)と与えた場合，\(MA_{it}\)は自地域の経済アウトカムと，距離の逆数で重み付けされた周辺地域の経済アウトカムの和になります．故に，地域間の距離が近い程大きな相互作用が生じることを表現する指標として見ることができます．なお，重みは距離二乗の逆数\(1/d^2_{ij}\)や指数関数\(\exp\left(-\tau d_{ij}\right)\)（\(\tau>0\)は何かしらの定数）で与えられる場合もあります．
今回の分析では，各地域の総人口を経済アウトカム\(X_{jt}\)に用いたMAを計算し，1人あたり所得との間の関係を検証します．具体的には，2000年，2005年，2010年の3時点の総人口を網羅したパネルデータを用いて各時点のMAを計算し，1人あたり所得との関係を固定効果法（fixed
effects method; FE）で検証します．
今回用いるデータの作成方法については，「Rによる市場アクセス分析」のためのパネルデータ構築の方にまとめています．元データはいずれもオープンアクセスで利用できる無償のデータです．手動でファイルをダウンロード・解凍したり，それらをローカル環境から読み込んだりする作業を極力削減するという点を特に重視しつつ，構築方法を記しています．Rによる空間分析の他，APIやWebスクレイピングの活用方法を学びたい方におすすめです．
参考文献
『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：空間解析・可視化編に，Rを用いた空間分析や計量経済分析の基礎と，その他参考になる文献をまとめています．
事前準備
パッケージの読み込み
今回の演習では，以下のパッケージを用います．未インストールのものについては，適宜インストールを行ってください．インストールの過程で「Do you want to install from sources the packages which need compilation? (Yes/no/cancel)」というメッセージが表示された場合には，とりあえず「no」を選択することをお勧めします．
- tidyverse：各種データ操作のためのパッケージ
- openxlsx：Rでxlsx形式のファイルを読み込み・書き出しするためのパッケージ
- modelsummary：回帰モデルの推定結果を整えて表示するためのパッケージ
- estimatr：頑健標準誤差や各種発展的な回帰モデルを推定可能なパッケージ
- sf：空間データ操作のためのパッケージ
- mapview：インタラクティブな空間データの可視化を行うためのパッケージ
- units：単位付きの数値変数を扱うためのパッケージ
- fixest：固定効果推定のためのパッケージ（補足でのみ使用）
library(tidyverse)
library(openxlsx)
library(modelsummary)
library(estimatr)
library(sf)
library(mapview)
#mapviewパッケージ使用時のおまじない
mapview::mapviewOptions(fgb=FALSE)
library(units)
library(fixest)データの読み込み
所得と人口に関するパネルデータを読み込みます．このデータは，2000年，2005年，2010年の3時点の観測値を網羅しています．このデータに含まれる総人口の変数popを用いてMAを計算し，1人あたり所得との関係を検証します．今回用いる主要な変数について，定義を示します．
- 被説明変数・説明変数（の元データ）
- ln_income：納税義務者1人あたり課税対象所得（対数）
- ln_dens：総面積あたり人口密度（対数）
- pop：総人口
- ln_income：納税義務者1人あたり課税対象所得（対数）
- コントロール変数
- ln_area：総面積（対数）
- altitude：平均標高
- slope_sum：傾斜度の合計
- ln_d_coast：海岸線までの距離（対数）
- ln_area：総面積（対数）
income_dens_panel <- openxlsx::read.xlsx(xlsxFile="https://www.dropbox.com/s/yc4m71fph42t1q3/income_dens_panel.xlsx?dl=1")後々使う変数として，2000年以外の時点について時点ダミーを作成します．例えば2005年度ダミーd_05は，レコードが2005年の観測値で成っていれば1，さもなくば0を取ります．
income_dens_panel <- income_dens_panel %>%
dplyr::mutate(d_05=ifelse(year==2005,1,0),
d_10=ifelse(year==2010,1,0))続いて，市区町村ポリゴンを読み込みます．この市区町村ポリゴンでは，政令指定都市は1つの自治体として扱われます．
ab_poly <- sf::read_sf(dsn="https://www.dropbox.com/s/w81179o3vk3og5w/ab_jp_dc.geojson?dl=1")
#データの先頭行を表示
head(ab_poly)## Simple feature collection with 6 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 361034.7 ymin: 4578541 xmax: 541182.7 ymax: 4787387
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 8
## id N03_001 N03_002 N03_003 N03_004 N03_007 type geometry
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON [m]>
## 1 gci:0… 北海道 石狩振… 札幌市 <NA> 01100 desi… (((541182.7 4763255, 539…
## 2 gci:0… 北海道 渡島総… <NA> 函館市 01202 <NA> (((475880.9 4628555, 476…
## 3 gci:0… 北海道 後志総… <NA> 小樽市 01203 <NA> (((523715.7 4781405, 521…
## 4 gci:0… 北海道 胆振総… <NA> 室蘭市 01205 <NA> (((504480.3 4688873, 503…
## 5 gci:0… 北海道 渡島総… 松前郡 松前町 01331 <NA> (((399946.4 4580045, 400…
## 6 gci:0… 北海道 檜山振… 久遠郡 せたな… 01371 <NA> (((420360.8 4711939, 420…読み込んだポリゴンについて，パネルデータ内に含まれる市区町村に該当するもののみを残します．
#パネルデータから変数JCODEを取り出し，重複を削除する
JCODE_panel <- unique(income_dens_panel$JCODE)
ab_poly <- ab_poly %>%
#JCODE変数（N03_007）がパネルデータ内のJCODEと一致するレコードのみ残す
dplyr::filter(N03_007%in%JCODE_panel)データの出典
- 所得・人口パネルデータ：社会・人口統計体系
- コントロール変数
- 標高・傾斜：標高・傾斜度5次メッシュデータ（国土数値情報）
- 海岸線距離：海岸線データ（国土数値情報）
- 標高・傾斜：標高・傾斜度5次メッシュデータ（国土数値情報）
- 市区町村ポリゴン：歴史的行政区域データセット
- 日本 市区町村（政令指定都市統合版）の粗解像度（2021-01-01）を使用
市区町村間距離の算出
MAを計算するためには，市区町村間距離\(d_{ij}\)が必要です．今回は，\(d_{ij}\)を市区町村ポリゴンの重心間の直線距離とした上で，それらをまとめた距離行列を作成します．距離行列の各要素の逆数を計算することで，MAにおける重みを計算することができます．
市区町村重心ポイントの作成
ポリゴンの重心点は，sfのst_centroid関数を用いて作成できます．
#市区町村ポリゴンの重心ポイントを作成
ab_poly_cent <- sf::st_centroid(x=ab_poly)距離行列の作成
ポイント間の距離行列を作成するには，sfのst_distance関数を用います．行列の次元は「ポイント数×ポイント数」で，全て表示しようとするとエラいことになるので，行と列それぞれについて5番目までを表示します．例えば行列の(2,1)要素は，1つ目のポイントと2つ目のポイント間の直線距離[m]\(d_{12}\)を表します．
#重心ポイント間の距離を計算し行列形式で返す
ab_poly_cent_dist <- sf::st_distance(x=ab_poly_cent)
#距離行列の先頭行・列
ab_poly_cent_dist[1:5,1:5]## Units: [m]
## 1 2 3 4 5
## 1 0.00 131205.24 25560.90 72893.79 191083.54
## 2 131205.24 0.00 146695.54 59019.53 80053.77
## 3 25560.90 146695.54 0.00 87724.82 198680.44
## 4 72893.79 59019.53 87724.82 0.00 121628.53
## 5 191083.54 80053.77 198680.44 121628.53 0.00現状距離には[m]という単位が付いていますが，後々の演算では使わない（むしろ邪魔になる）ので，unitsのset_units関数を用いて単位を削除します．
#距離の単位を削除
ab_poly_cent_dist <- units::set_units(x=ab_poly_cent_dist,value=NULL)
ab_poly_cent_dist[1:5,1:5]## 1 2 3 4 5
## 1 0.00 131205.24 25560.90 72893.79 191083.54
## 2 131205.24 0.00 146695.54 59019.53 80053.77
## 3 25560.90 146695.54 0.00 87724.82 198680.44
## 4 72893.79 59019.53 87724.82 0.00 121628.53
## 5 191083.54 80053.77 198680.44 121628.53 0.00行列の各要素を1000で割ることにより，距離の単位を[m]から[km]に変換します．
#距離単位を[m]から[km]に変更
ab_poly_cent_dist <- ab_poly_cent_dist/1000
ab_poly_cent_dist[1:5,1:5]## 1 2 3 4 5
## 1 0.00000 131.20524 25.56090 72.89379 191.08354
## 2 131.20524 0.00000 146.69554 59.01953 80.05377
## 3 25.56090 146.69554 0.00000 87.72482 198.68044
## 4 72.89379 59.01953 87.72482 0.00000 121.62853
## 5 191.08354 80.05377 198.68044 121.62853 0.00000距離行列の対角要素は自地域内距離\(d_{ii}\)を表し，全て0となっています（同じ位置にあるポイント上で距離を計算しても0なので）．0の逆数をそのまま取ってしまうとInfになってしまうので，ここでは便宜上1と与えます．diag関数を用いて行列の対角要素にアクセスし，1を代入します．
#行列の対角要素を1に置き換え
diag(ab_poly_cent_dist) <- 1MAの計算
重みの作成
距離行列を用いて，距離の逆数による重み\(1/d_{ij}\)を各要素に持つ行列を作成します．行列Aの各要素について逆数を取る場合には，単純に1/Aとします．
ab_poly_cent_ddecay <- 1/ab_poly_cent_dist
ab_poly_cent_ddecay[1:5,1:5]## 1 2 3 4 5
## 1 1.000000000 0.007621647 0.039122254 0.013718590 0.005233313
## 2 0.007621647 1.000000000 0.006816840 0.016943545 0.012491603
## 3 0.039122254 0.006816840 1.000000000 0.011399282 0.005033208
## 4 0.013718590 0.016943545 0.011399282 1.000000000 0.008221755
## 5 0.005233313 0.012491603 0.005033208 0.008221755 1.000000000MAの計算
重み行列と各時点の人口から，各時点の\(MA_{it}\)を計算します．重み行列を\(\mathbf{D}\)（コード内ではab_poly_cent_ddecay），時点\(t\)における各市区町村の人口をまとめたベクトルを\(\mathbf{p}_{t}\)とすると，以下の行列演算を行うことで\(MA_{it}\)が計算できます．
\[ \begin{equation} \mathbf{D}\mathbf{p}_{t}=\begin{bmatrix} 1 & 1/d_{2,1} & \dots & 1/d_{n,1} \\ 1/d_{1,2} & 1 & \dots & 1/d_{n,2} \\ \vdots & \vdots & \vdots & \vdots \\ 1/d_{1,n} & 1/d_{2,n} & \dots & 1 \\ \end{bmatrix} \begin{bmatrix} pop_{1t} \\ pop_{2t} \\ \vdots \\ pop_{nt} \\ \end{bmatrix}=\begin{bmatrix} \sum_{j} \frac{pop_{jt}}{d_{1j}} \\ \sum_{j} \frac{pop_{jt}}{d_{2j}} \\ \vdots \\ \sum_{j} \frac{pop_{jt}}{d_{nj}} \\ \end{bmatrix}=\begin{bmatrix} MA_{1t} \\ MA_{2t} \\ \vdots \\ MA_{nt} \\ \end{bmatrix} \end{equation} \]
MAの計算は年毎に行うので，まずはパネルデータを変数yearでグループ化します．これにより，各時点の人口popを別々のベクトルとして扱えます．次に，mutate関数内で上で示したような行列・ベクトル間の掛け算を行います．行列・ベクトルの積の演算は，「%*%」という演算子で行えます．
income_dens_panel_ma <- income_dens_panel %>%
#年でグループ化
dplyr::group_by(year) %>%
#各時点のMAを計算し変数として追加
dplyr::mutate(ma=as.numeric(ab_poly_cent_ddecay%*%pop)) %>%
#グループ化を解除
dplyr::ungroup() %>%
#MAの対数値を変数として追加
dplyr::mutate(ln_ma=log(ma))補足
MAの計算はsapply関数を用いて行うことも可能です．引数XにはMAの計算対象としたい時点のベクトルを，FUNにはMAを計算するための手順をfunction関数で囲った上で指定します．sapply関数を実行すると，Xに指定したベクトルの各要素がFUNで指定した関数に渡され，結果がベクトル形式で返ってきます．
ma <- sapply(X=c(2000,2005,2010),FUN=function(x){
#パネルデータに関する処理
pop <- income_dens_panel %>%
#[x]年のレコードを抽出
dplyr::filter(year==x) %>%
#popをベクトルとして取り出し
dplyr::pull(pop)
#行列演算でMAを計算
ma_tmp <- ab_poly_cent_ddecay%*%pop
#計算結果を返す
return(ma_tmp)
})MAの可視化
計算されたMAの空間分布を地図上で可視化します．以下では例として，2000年時点のMAを可視化する手順を示します．まずは，パネルデータから2000年のレコードを抽出した上で，JCODEをキーに市区町村ポリゴンへ結合します．
ma_00 <- income_dens_panel_ma %>%
dplyr::filter(year==2000)
ab_poly_ma_00 <- ab_poly %>%
dplyr::left_join(y=ma_00,by=c("N03_007"="JCODE"))mapview関数を用いて，市区町村ポリゴンをln_maの大きさに基づいて可視化します．
mapview::mapview(x=ab_poly_ma_00,
zcol="ln_ma",
color="transparent")空間データの可視化には，ggplot2のgeom_sf関数を用いることもできます．
ggplot2::ggplot() +
ggplot2::geom_sf(data=ab_poly_ma_00,aes(fill=ln_ma),lty=0) +
ggplot2::scale_fill_gradient(low="white",high="red")
記述統計分析
変数間の相関係数
変数間の相関係数行列を計算し，変数間の関連の強さを確認しておきます．今回はベンチマークとして，人口密度ln_densを変数に用いた分析も行います．
cormat <- income_dens_panel_ma %>%
#分析に用いる変数のみ残す
dplyr::select(ln_income,ln_dens,ln_ma) %>%
#相関係数行列の算出
cor() %>%
#相関係数を小数点以下3桁で四捨五入する
round(digits=3)
cormat## ln_income ln_dens ln_ma
## ln_income 1.000 0.564 0.556
## ln_dens 0.564 1.000 0.709
## ln_ma 0.556 0.709 1.000所得とMAの散布図
1人あたり所得ln_incomeとMAの対数値ln_maの関係を，散布図を用いて把握します．まずは，全てのデータを用いたプロットです．
ggplot2::ggplot(data=income_dens_panel_ma,
aes(x=ln_ma,y=ln_income)) +
#散布図を作成
ggplot2::geom_point() +
#近似直線を重ねてプロット
ggplot2::geom_smooth()
続いて，市区部と町村部に分けて散布図をプロットします．市区部では所得とMAとの間に直線的な関係が見られる一方，町村部ではU字形の関係が見られます．
ggplot2::ggplot(data=income_dens_panel_ma,
#変数municipality_classを基準に点を色分け
aes(x=ln_ma,y=ln_income)) +
#変数municipality_classで指定したグループ毎にウィンドウを分ける
ggplot2::facet_wrap(facets="municipality_class") +
#散布図のプロット
ggplot2::geom_point() +
#近似曲線のプロット
ggplot2::geom_smooth() +
#Macで日本語を描画する際には日本語フォントの指定が必要
ggplot2::theme_gray(base_family="HiraKakuPro-W3")
FE
定式化
今回推定するベースラインのモデルを，全ての固定効果を含む場合について定式化すると，以下のようになります．
\[ \begin{equation} \ln income_{it}= \beta_{1} \ln MA_{it} + CityFE_{i} + YearFE_{t} + \varepsilon_{it}. \end{equation} \]
ベースラインの推定
MAの対数値ln_maのみを説明変数に用いたFE推定を行います．その際，人口密度の対数ln_densを説明変数に用いたFE推定をベンチマークとして行い，結果を比較します．
fe_ma <- estimatr::lm_robust(ln_income~ln_ma,
data=income_dens_panel_ma,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
fe_dens <- estimatr::lm_robust(ln_income~ln_dens,
data=income_dens_panel_ma,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
modelsummary::modelsummary(models=list("FE (MA)"=fe_ma,
"FE (Density)"=fe_dens),
stars=TRUE)| FE (MA) | FE (Density) | |
|---|---|---|
| ln_ma | 0.336*** | |
| (0.096) | ||
| ln_dens | 0.141*** | |
| (0.019) | ||
| Num.Obs. | 5220 | 5220 |
| R2 | 0.977 | 0.978 |
| R2 Adj. | 0.966 | 0.968 |
| AIC | −24078.8 | −24355.0 |
| BIC | −24065.7 | −24341.9 |
| RMSE | 0.02 | 0.02 |
| Std.Errors | by: JCODE | by: JCODE |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
市区・町村サブサンプルによる推定
市区部・町村部それぞれのサブサンプルについて分析を行います．
#市区部に関する推定
fe_ma_u <- estimatr::lm_robust(ln_income~ln_ma,
data=dplyr::filter(income_dens_panel_ma,municipality_class=="市区"),
clusters=JCODE,
fixed_effects=~JCODE+year,
se_type="stata")
#町村部に関する推定
fe_ma_r <- estimatr::lm_robust(ln_income~ln_ma,
data=dplyr::filter(income_dens_panel_ma,municipality_class=="町村"),
clusters=JCODE,
fixed_effects=~JCODE+year,
se_type="stata")
modelsummary::modelsummary(models=list("FE (Urban)"=fe_ma_u,
"FE (Rural)"=fe_ma_r),
stars=TRUE)| FE (Urban) | FE (Rural) | |
|---|---|---|
| ln_ma | 0.533*** | −0.393** |
| (0.107) | (0.144) | |
| Num.Obs. | 2445 | 2775 |
| R2 | 0.989 | 0.960 |
| R2 Adj. | 0.984 | 0.940 |
| AIC | −12944.7 | −11928.7 |
| BIC | −12933.1 | −11916.8 |
| RMSE | 0.02 | 0.03 |
| Std.Errors | by: JCODE | by: JCODE |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
コントロール変数を含む推定
総面積ln_area，平均標高altitude，傾斜度の合計slope_sum，海岸線ln_d_coastまでの距離をコントロール変数\(Geography_{i}\)として加えた以下の定式化に基づく分析を行います．これら変数は，パネルデータの期間中は時間不変なので，年ダミー（時点固定効果）との交差項を説明変数として加えます．
\[ \begin{equation} \ln income_{it}= \beta_{1} \ln MA_{it} + CityFE_{i} + YearFE_{t} + Geography_{i}×YearFE_{t} + \varepsilon_{it}. \end{equation} \]
コントロール変数を加えたモデルを推定し，ベースラインの結果と比較します．
#FE推定（含コントロール）
fe_ma_c <- estimatr::lm_robust(ln_income~ln_ma+ln_area:d_05+altitude:d_05+slope_sum:d_05+ln_d_coast:d_05+
ln_area:d_10+altitude:d_10+slope_sum:d_10+ln_d_coast:d_10,
data=income_dens_panel_ma,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
#ベースラインの結果と比較
modelsummary::modelsummary(models=list("FE"=fe_ma,
"FE"=fe_ma_c),
stars=TRUE)| FE | FE | |
|---|---|---|
| ln_ma | 0.336*** | 0.511*** |
| (0.096) | (0.107) | |
| ln_area × d_05 | 0.003* | |
| (0.001) | ||
| d_05 × altitude | 0.000** | |
| (0.000) | ||
| d_05 × slope_sum | 0.000 | |
| (0.000) | ||
| d_05 × ln_d_coast | 0.001 | |
| (0.001) | ||
| ln_area × d_10 | 0.008*** | |
| (0.002) | ||
| altitude × d_10 | 0.000*** | |
| (0.000) | ||
| slope_sum × d_10 | 0.000 | |
| (0.000) | ||
| ln_d_coast × d_10 | −0.003+ | |
| (0.002) | ||
| Num.Obs. | 5220 | 5220 |
| R2 | 0.977 | 0.978 |
| R2 Adj. | 0.966 | 0.967 |
| AIC | −24078.8 | −24303.6 |
| BIC | −24065.7 | −24238.0 |
| RMSE | 0.02 | 0.02 |
| Std.Errors | by: JCODE | by: JCODE |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
コントロール変数の推定結果を表示したくない場合，modelsummary関数の引数coef_omitで，表示したくない変数に含まれる文字列（今回の場合だと「d_」）を指定します．
#ベースラインの結果と比較
modelsummary::modelsummary(models=list("FE"=fe_ma,
"FE (w/ control)"=fe_ma_c),
stars=TRUE,
coef_omit="d_")| FE | FE (w/ control) | |
|---|---|---|
| ln_ma | 0.336*** | 0.511*** |
| (0.096) | (0.107) | |
| Num.Obs. | 5220 | 5220 |
| R2 | 0.977 | 0.978 |
| R2 Adj. | 0.966 | 0.967 |
| AIC | −24078.8 | −24303.6 |
| BIC | −24065.7 | −24238.0 |
| RMSE | 0.02 | 0.02 |
| Std.Errors | by: JCODE | by: JCODE |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
補足
異なるMAの定式化
MAの計算
本編で用いた距離に関する重みは，距離逆数\(1/d_{ij}\)でした．ここでは，距離二乗の逆数\(1/d^2_{ij}\)や指数関数\(\exp\left(-\tau d_{ij}\right)\)を重みに用いたMAの計算と，それに基づく分析を追加的に行います．距離二乗の場合，重み行列は以下のように計算できます．
#距離行列を計算->単位を削除->要素を1000で割る->要素を二乗する
ab_poly_cent_dist_sq <- (units::set_units(sf::st_distance(x=ab_poly_cent),value=NULL)/1000)^2
#距離二乗行列の逆数を計算
ab_poly_cent_ddecay_sq <- 1/ab_poly_cent_dist_sq
#対角要素を1に置き換え
diag(ab_poly_cent_ddecay_sq) <- 1
ab_poly_cent_ddecay_sq[1:5,1:5]## 1 2 3 4 5
## 1 1.000000e+00 5.808950e-05 1.530551e-03 1.881997e-04 2.738757e-05
## 2 5.808950e-05 1.000000e+00 4.646931e-05 2.870837e-04 1.560402e-04
## 3 1.530551e-03 4.646931e-05 1.000000e+00 1.299436e-04 2.533318e-05
## 4 1.881997e-04 2.870837e-04 1.299436e-04 1.000000e+00 6.759726e-05
## 5 2.738757e-05 1.560402e-04 2.533318e-05 6.759726e-05 1.000000e+00指数関数の場合，重み行列は以下のように計算できます．例えば今回は\(\tau=0.1\)と設定します，指数関数の場合は距離行列の対角成分を1に置き換える操作は行いません（\(\exp\left(0\right)=1\)なので）．
#距離行列を計算->単位を削除->要素を1000で割る->指数関数を適用する
ab_poly_cent_ddecay_exp <- exp(-0.1*(units::set_units(sf::st_distance(x=ab_poly_cent),value=NULL)/1000))
ab_poly_cent_ddecay_exp[1:5,1:5]## 1 2 3 4 5
## 1 1.000000e+00 2.003683e-06 7.760760e-02 6.827520e-04 5.027445e-09
## 2 2.003683e-06 1.000000e+00 4.256902e-07 2.734101e-03 3.336635e-04
## 3 7.760760e-02 4.256902e-07 1.000000e+00 1.549385e-04 2.351896e-09
## 4 6.827520e-04 2.734101e-03 1.549385e-04 1.000000e+00 5.220838e-06
## 5 5.027445e-09 3.336635e-04 2.351896e-09 5.220838e-06 1.000000e+00計算されたそれぞれの重み行列を用いて，MAとその対数値を計算します．距離二乗逆数の重みに基づくMAの対数値を変数ln_ma_sq，指数関数の重みに基づくMAの対数値を変数ln_ma_expとします．
income_dens_panel_ma_alt <- income_dens_panel_ma %>%
#年でグループ化
dplyr::group_by(year) %>%
#各時点のMA（距離二乗逆数）を計算し変数として追加
dplyr::mutate(ma_sq=as.numeric(ab_poly_cent_ddecay_sq%*%pop)) %>%
#各時点のMA（指数関数）を計算し変数として追加
dplyr::mutate(ma_exp=as.numeric(ab_poly_cent_ddecay_exp%*%pop)) %>%
#グループ化を解除
dplyr::ungroup() %>%
#MAの対数値を変数として追加
dplyr::mutate(ln_ma_sq=log(ma_sq),
ln_ma_exp=log(ma_exp))計算されたMAの空間分布を比較します．以下では例として，2000年の各種MAの分布を可視化します．それに先立ち，2000年の観測値を抽出した上で，市区町村ポリゴンに結合します．
ma_00_alt <- income_dens_panel_ma_alt %>%
dplyr::filter(year==2000)
ab_poly_ma_00_alt <- ab_poly %>%
dplyr::left_join(y=ma_00_alt,by=c("N03_007"="JCODE"))距離二乗逆数，指数関数の重みに基づいて計算されたMAを可視化します．
#距離二乗逆数のプロット
ggplot2::ggplot() +
ggplot2::geom_sf(data=ab_poly_ma_00_alt,aes(fill=ln_ma_sq),lty=0) +
ggplot2::scale_fill_gradient(low="white",high="red")
#指数関数のプロット
ggplot2::ggplot() +
ggplot2::geom_sf(data=ab_poly_ma_00_alt,aes(fill=ln_ma_exp),lty=0) +
ggplot2::scale_fill_gradient(low="white",high="red")
FE
距離二乗逆数，指数関数をそれぞれ重みに用いた場合のMAの対数値を説明変数として用いたFE推定を行います．
#FE推定（距離二乗逆数MA）
fe_ma_sq <- estimatr::lm_robust(ln_income~ln_ma_sq,
data=income_dens_panel_ma_alt,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
#FE推定（距離二乗逆数MA・含コントロール）
fe_ma_sq_c <- estimatr::lm_robust(ln_income~ln_ma_sq+ln_area:d_05+altitude:d_05+slope_sum:d_05+ln_d_coast:d_05+
ln_area:d_10+altitude:d_10+slope_sum:d_10+ln_d_coast:d_10,
data=income_dens_panel_ma_alt,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
#FE推定（指数関数MA）
fe_ma_exp <- estimatr::lm_robust(ln_income~ln_ma_exp,
data=income_dens_panel_ma_alt,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
#FE推定（指数関数MA・含コントロール）
fe_ma_exp_c <- estimatr::lm_robust(ln_income~ln_ma_exp+ln_area:d_05+altitude:d_05+slope_sum:d_05+ln_d_coast:d_05+
ln_area:d_10+altitude:d_10+slope_sum:d_10+ln_d_coast:d_10,
data=income_dens_panel_ma_alt,
clusters=JCODE,
#市区町村FEと時点FE両方を含める
fixed_effects=~JCODE+year,
se_type="stata")
modelsummary::modelsummary(models=list("FE (Square)"=fe_ma_sq,
"FE (Square, w/ control)"=fe_ma_sq_c,
"FE (Exponential)"=fe_ma_exp,
"FE (Exponential, w/ control)"=fe_ma_exp_c),
stars=TRUE,
coef_omit="d_")| FE (Square) | FE (Square, w/ control) | FE (Exponential) | FE (Exponential, w/ control) | |
|---|---|---|---|---|
| ln_ma_sq | 0.150*** | 0.198*** | ||
| (0.024) | (0.028) | |||
| ln_ma_exp | 0.138*** | 0.179*** | ||
| (0.030) | (0.035) | |||
| Num.Obs. | 5220 | 5220 | 5220 | 5220 |
| R2 | 0.978 | 0.979 | 0.977 | 0.978 |
| R2 Adj. | 0.967 | 0.968 | 0.966 | 0.967 |
| AIC | −24197.6 | −24426.6 | −24101.3 | −24292.2 |
| BIC | −24184.5 | −24361.0 | −24088.2 | −24226.6 |
| RMSE | 0.02 | 0.02 | 0.02 | 0.02 |
| Std.Errors | by: JCODE | by: JCODE | by: JCODE | by: JCODE |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
fixestを用いた推定
FE推定にはfixestのfeols関数を用いることも可能です．以下の推定では，回帰係数の標準誤差は1つ目の固定効果（JCODE）の単位に基づいてクラスター化されます．ちなみにfixestには，ポアソンやロジットといった非線形回帰の場合のFE推定を行うための関数も実装されています．
#FE推定（含コントロール）
fe_ma_c_fixest <- fixest::feols(fml=ln_income~ln_ma+ln_area:d_05+altitude:d_05+slope_sum:d_05+ln_d_coast:d_05+
ln_area:d_10+altitude:d_10+slope_sum:d_10+ln_d_coast:d_10|JCODE+year,
data=income_dens_panel_ma,
ssc=ssc(fixef.K="full"))
#estimatrの関数を用いた時の結果と比較
modelsummary::modelsummary(models=list("FE (lm_robust)"=fe_ma_c,
"FE (feols)"=fe_ma_c_fixest),
stars=TRUE,
coef_omit="d_")| FE (lm_robust) | FE (feols) | |
|---|---|---|
| ln_ma | 0.511*** | 0.511*** |
| (0.107) | (0.107) | |
| Num.Obs. | 5220 | 5220 |
| R2 | 0.978 | 0.978 |
| R2 Adj. | 0.967 | 0.967 |
| R2 Within | 0.058 | |
| R2 Within Adj. | 0.056 | |
| AIC | −24303.6 | −20821.6 |
| BIC | −24238.0 | −9334.6 |
| RMSE | 0.02 | 0.02 |
| Std.Errors | by: JCODE | by: JCODE |
| FE: JCODE | X | |
| FE: year | X | |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
時間可変な重み行列によるMA
今回MAの計算に用いた重み行列は，市区町村重心間の直線距離に基づく点で時間不変です．一方，交通ネットワーク上での距離や所要時間を重みとして用いる場合，交通インフラの新設・廃止等によって値が時点毎に異なります．こうした時間可変の重み行列に基づくMAは本編で示した場合ほど計算が単純ではないので，別途方法を検討する必要があります．以下では最も素朴な方法として，apply族を用いたMAの計算手順を示します．
各時点の重み行列の縦結合
各時点の重み行列を用意します．今回は，時間不変の重み行列を「あたかも時間可変として扱う」ため，それぞれ別の名前のオブジェクトに代入しています．
ab_poly_cent_ddecay_00 <- ab_poly_cent_ddecay
ab_poly_cent_ddecay_05 <- ab_poly_cent_ddecay
ab_poly_cent_ddecay_10 <- ab_poly_cent_ddecay各時点の重み行列をtibble形式（data.frameの拡張版）に変換した上で，時点を表す変数yearを追加します．
ab_poly_cent_ddecay_00 <- ab_poly_cent_ddecay_00 %>%
dplyr::as_tibble() %>%
dplyr::mutate(year=2000)
ab_poly_cent_ddecay_05 <- ab_poly_cent_ddecay_05 %>%
dplyr::as_tibble() %>%
dplyr::mutate(year=2005)
ab_poly_cent_ddecay_10 <- ab_poly_cent_ddecay_10 %>%
dplyr::as_tibble() %>%
dplyr::mutate(year=2010)rbind関数を用いて，重み行列を縦方向に結合します．
ab_poly_cent_ddecay_full <- rbind(ab_poly_cent_ddecay_00,
ab_poly_cent_ddecay_05,
ab_poly_cent_ddecay_10)MAの計算
縦結合された重み行列を用いて，例えば2000年のMAを計算する手順を以下に示します．それに先立ち，パネルデータから2000年のレコードを抽出した上で，変数popをベクトルとして取り出します．同様に，重み行列からも2000年のレコードを抽出した上で行列に変換します．
#パネルデータに関する処理
pop <- income_dens_panel %>%
#2000年のレコードを抽出
dplyr::filter(year==2000) %>%
#popをベクトルとして取り出し
dplyr::pull(pop)
#重み行列に関する処理
ddcay <- ab_poly_cent_ddecay_full %>%
#2000年のレコードを抽出
dplyr::filter(year==2000) %>%
#yearを変数から落とす
dplyr::select(-year) %>%
#行列に変換
as.matrix()本編で示した行列演算を用いて2000年のMAを計算し，計算結果をtibble形式に変換した上で，後々分かりやすいように変数yearを作成し直します．
ma_out <- ddcay%*%pop
ma_out <- ma_out %>%
dplyr::as_tibble() %>%
dplyr::mutate(year=2000) 上で示した単時点でのMAの計算を，複数時点について同時に行う際には，lapply関数を用いるのが便利です．引数XにはMAの計算対象としたい時点のベクトルを，FUNにはMAを計算するための手順をfunction関数で囲った上で指定します．lapply関数を実行すると，Xに指定したベクトルの各要素がFUNで指定した関数に渡され，結果がlist形式で返ってきます．
function関数の中身は，単時点でのMAの計算手順と同様ですが，filter関数の中で「year==2000」となっていたものが「year==x」となっている点が異なります．このxに，ベクトルXの各要素（今回のケースでは2000,2005,2010）が渡され，MAを計算する関数が実行されます．
ma <- lapply(X=c(2000,2005,2010),FUN=function(x){
#パネルデータに関する処理
pop <- income_dens_panel %>%
#[x]年のレコードを抽出
dplyr::filter(year==x) %>%
#popをベクトルとして取り出し
dplyr::pull(pop)
#重み行列に関する処理
ddcay <- ab_poly_cent_ddecay_full %>%
#[x]年のレコードを抽出
dplyr::filter(year==x) %>%
#yearを変数から落とす
dplyr::select(-year) %>%
#行列に変換
as.matrix()
#行列演算でMAを計算
ma_out <- ddcay%*%pop
#計算結果を整える
ma_out <- ma_out %>%
#tibble形式に変換
dplyr::as_tibble() %>%
#変数yearを追加
dplyr::mutate(year=x)
#計算結果を返す
return(ma_out)
})現状list内の各要素に入っている時点別のMAの計算結果を縦方向に結合し，ひとつのtibbleにまとめます．その際には，dplyrのbind_rows関数を用います．その後，縦結合を実行した結果について，変数名を分かりやすいものに変更します．
#lapply関数の実行結果を縦結合
ma <- dplyr::bind_rows(ma) %>%
#変数名をV1からmaに変更
dplyr::rename(ma=V1)
head(ma)## # A tibble: 6 × 2
## ma year
## <dbl> <dbl>
## 1 2006653. 2000
## 2 518254. 2000
## 3 395372. 2000
## 4 572537. 2000
## 5 323143. 2000
## 6 401623. 2000