その無茶振り，（Rで）GISが解決します：フローデータ分析編
Keisuke Takano
12/18/2022
はじめに
このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．
バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．
予めご容赦ください．
やること
拙稿「その無茶振り，（Rで）GISが解決します」では，sfパッケージ等を用いた，Rによる空間分析の基礎を解説しました．このページはその応用として，空間データを用いて重力モデル（空間的相互作用モデル）推定に必要なデータを構築する方法を解説します．
重力モデルとは
重力モデルは，Figure
1に示す通り，引力の大きさは質量に比例し，距離に反比例するという物理学の万有引力を模した統計モデルです．社会科学では，質量をGDPのような経済活動の規模に，引力をヒト・モノ・コトの移動量（フローとも呼びます）に置き換えたようなモデルを推定します．
重力モデル自体は，計量地理学や国際貿易論，空間経済学のような地理空間を扱う分野ではオーソドックスな手法です．一方，初学者にとっては，データ構築の段階から様々手こずることが多いと思います．昨今では，tidyverseとその関連パッケージ群の発展により，データ構築からモデル推定までを，Rのみで完結して行えるようになってきましたので，その具体例をここでは示します．
Figure 1：万有引力の法則
ODフローデータ
重力モデルを推定する際，アウトカムはしばしばFigure
2に示すようなOD（origin-destination）表に基づいて作られます．ODデータは，有向の重み付きグラフと見ることができるので，各種ネットワーク分析の手法が適用可能です．
Rによるネットワーク分析における近年の進化として，tidygraphやggraphのような，tidyverseと組み合わせて使うことができるパッケージ群が開発されたことが挙げられます．そこで本稿では，ODフローデータからネットワークを構築し，中心性指標といったネットワーク分析の手法を実装する手順についても解説します．
Figure 2：OD表
関連リンク
なお，このページで用いられるデータの前処理・原典については，その無茶振り，（Rで）GISが解決します：フローデータ分析編（e-Stat使い倒しの巻）を参照してください．e-Stat APIを用いて政府統計を取得する方法が，合わせて学べます．加えて，このWEBページの元になっているRmd（r_de_geonw.Rmd）と演習データは，こちらのDropboxリンクに設置してあります．
今回の無茶振り
都市計画系シンクタンクに勤める私に，とある県の担当者からこのようなメールが届きました．
我が県では昨今，若者の県外流出を重要な問題として認識しています．
特に，高校を卒業した段階で他県の大学に進学し，地元に戻ってこない若者が多い状況が数十年に渡って続いています．
この問題を是正するため，どのような政策を行っていくべきかを目下検討していますが，その基礎資料として，大学生の県外進学の動向を知りたいと考えています．
県外進学に限らず，定住を伴う人の移動を分析する際には，プッシュ要因とプル要因を考えることが重要です．ここでは，山田・徳岡（2007）『地域経済学入門』の6章の内容に基づきながら，それら要因について簡単に説明します．
プッシュ要因とは，現在の場所を離れる決定に影響を与える要因です．一方プル要因とは，移転先の決定に影響を与える要因を指します．プッシュ・プル要因には，賃金（所得）や失業率のような雇用に関する要因や，家賃や物価のような家計に関する要因，公共サービスの水準や社会資本といった地域環境に関する要因が含まれます．人口移動は，これらプッシュ・プル要因の相対的な関係で決まるのです．
一方，距離の増加によって人口移動は減ることが多いです．その原因としては例えば，現在の場所で築いた人間関係の喪失，移転先に関する情報の入手困難性，食事・気候の違いによる適応の難しさが挙げられると思います．大学生の県外進学の場合はとりわけ，実家に帰るのが大変というのが理由になりそうです．今回は，これらプッシュ・プル要因及び都道府県間の距離と，県外進学者の数との関係を，重力モデルを用いて分析します．
参考文献
重力モデルに関して学びたい時には，以下の文献が参考になります．本間（2018）は都市工学，友原（2018）は国際貿易の専門家によって書かれたテキストで，重力モデルを用いる際の課題のフォーカスが，分野によって少しずつ異なることが分かります．
- 本間裕大. (2018). 空間相互作用モデル. In 貞広幸雄, 山田育穂, &
石井儀光. (編著). 空間解析入門, pp.43–50, 朝倉書店.
- 友原章典. (2018). 理論と実証から学ぶ新しい国際経済学. ミネルヴァ書房.
WEBページとしては，例えば田中鮎夢先生による連載「国際貿易と貿易政策研究メモ」が参考になります．国際貿易論全般に関する連載ですが，重力モデルの推定に関連する回は以下の通りです．
ネットワーク分析に関する最低限の知識を持っておきたい場合には，Takashi Kitanoさんによる以下のスライドが参考になります．筆者の知る限り，tidygraphやggraphに関する数少ない日本語資料のひとつです．
社会科学におけるネットワーク分析の応用に関して学びたい時には，以下の文献が参考になります．前者は，Rコードだけでなく，社会ネットワーク分析の各種手法について，比較的詳細に解説を行っているため，理論・実践両方を学びたい時に適しています．後者は，ネットワーク分析に限らず，政治科学の分野で用いられる様々なデータ分析手法も網羅したテキストですので，データ分析全般のテキストを1冊持っておきたい場合は，こちらをお勧めします．
- 鈴木努. (2017). ネットワーク分析, 第2版. 共立出版.
- Urdinez, F., & Cruz, A. (2020). R for Political Data Science: A Practical Guide. CRC Press.
なお，重力モデルを用いた県外進学フローの分析は，より発展的なモデルを用いて，以下の論文においてもなされています．大学生の県外進学というトピックをより詳しく知りたい場合は，合わせて参照することをお勧めします．
- 田村一軌. (2017). 大学進学にともなう都道府県間人口移動の定量分析. AGI Working Paper Series 2017-03. アジア成長研究所.
より発展的な分析のために
今回扱う重力モデルは極めて初歩的なものですが，より発展的な分析もR上で行うことが可能です．例えば，国際貿易の分野で扱われるような，構造的重力モデルを用いた分析を行いたい場合，gravityパッケージなどに基づく以下の文献が参考になります．
- The
gravity model of international trade: a user guide (H. Doytchinova et
al.)
- Gravity Models in R (Wolwer et al.)
一方，計量地理学や交通工学などで扱われるような制約型の重力モデルを用いたい場合，simodelsパッケージが有用です．
事前準備
必要パッケージの読み込み
今回の演習では，以下のパッケージを用います．未インストールのものについては，適宜インストールを行ってください．インストールの過程で「Do you want to install from sources the packages which need compilation? (Yes/no/cancel)」というメッセージが表示された場合には，とりあえず「no」を選択することをお勧めします．
- tidyverse：各種データ操作のためのパッケージ
- tidygraph：ネットワーク分析のためのパッケージ
- ggraph：ネットワーク可視化のためのパッケージ
- sf：空間データ（ベクタデータ）操作のためのパッケージ
- openxlsx：Rでxlsx形式のファイルを読み込み・書き出しするためのパッケージ
- estimatr：頑健標準誤差の計算を含む回帰分析を行うためのパッケージ
- modelsummary：回帰モデルの推定結果を整えて表示する為のパッケージ
library(tidyverse)
library(tidygraph)
library(ggraph)
library(sf)
#sfパッケージ使用時のおまじない
sf::sf_use_s2(FALSE)
library(openxlsx)
library(estimatr)
library(modelsummary)動作環境によっては，sfパッケージの読み込み時にエラーが出る場合があります．その原因は，sfパッケージが依存しているRcppのバージョンが古いことが多いようです．もしエラーが出る場合には，以下のコードを実行して見てください（筆者のPCではエラーが出なかったので，コメント化してあります）．
# install.packages("Rcpp")
# library(Rcpp)データの読み込み
分析を行うためには，アウトカムとして用いる県外進学者数や，進学におけるプッシュ・プル要因の情報が必要そうです．加えて，都道府県間の位置関係に関する情報も必要そうです．それぞれに対応するデータとして，以下に挙げたものを用意しました．進学者数及びプッシュ・プル要因のデータは，いずれも2010年度に観測されたものです．
- 県外進学者数
- 出典：学校基本調査
- 出典：学校基本調査
- プッシュ・プル要因
- 出典：社会・人口統計体系（統計でみる都道府県のすがた）
- 変数の定義はこちらのpdfを参照
- 出典：社会・人口統計体系（統計でみる都道府県のすがた）
- 都道府県庁ポイント
- 都道府県ポリゴン
ベクタデータ
今回用意したベクタデータは全てGeoJSON形式です．GeoJSONの読み込みにはsfのread_sf関数を用います．今回は，筆者が用意したDropboxリンクから直接ダウンロードします．
#都道府県庁ポイント
pref_office <- sf::read_sf(dsn="https://www.dropbox.com/s/7i4ck5lh829dluv/pref_office.geojson?dl=1")
head(pref_office)## Simple feature collection with 6 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 423139.2 ymin: 4232686 xmax: 528327.9 ymax: 4767974
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 2
## pref_code geometry
## <chr> <POINT [m]>
## 1 01 (528327.9 4767974)
## 2 02 (478083 4519290)
## 3 03 (513081.3 4394877)
## 4 04 (488789.6 4235648)
## 5 05 (423139.2 4396850)
## 6 06 (444280.3 4232686)#都道府県ポリゴン
pref_poly <- sf::read_sf(dsn="https://www.dropbox.com/s/r95qsekjlvu761h/pref_poly.geojson?dl=1")
head(pref_poly)## Simple feature collection with 6 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 361034.7 ymin: 4176685 xmax: 1116586 ymax: 5074233
## Projected CRS: WGS 84 / UTM zone 54N
## # A tibble: 6 × 5
## id N03_001 N03_002 N03_007 geometry
## <chr> <chr> <chr> <chr> <MULTIPOLYGON [m]>
## 1 北海道 北海道 <NA> 01 (((399115.3 4579312, 399946.4 4580045, 400928 …
## 2 青森県 青森県 <NA> 02 (((544617.2 4487185, 544110.9 4487857, 544786.…
## 3 岩手県 岩手県 <NA> 03 (((586024.8 4400680, 584783.5 4400086, 584436.…
## 4 宮城県 宮城県 <NA> 04 (((537101.5 4238584, 537539.9 4238346, 535813.…
## 5 秋田県 秋田県 <NA> 05 (((476768 4379504, 476719.3 4379278, 476587.4 …
## 6 山形県 山形県 <NA> 06 (((436844.9 4211918, 436868.9 4211840, 436928.…テキストデータ
県外進学者数データに含まれる変数の定義は以下の通りです．
- origin_name：出発地の都道府県名
- destination_name：到着地の都道府県名
- origin：出発地の都道府県コード
- destination：到着地の都道府県コード
- total：出発地・到着地ペア間の大学進学者総数
- national：出発地・到着地ペア間の国立大学進学者数
- private：出発地・到着地ペア間の私立大学進学者数
- male：出発地・到着地ペア間の大学進学者総数（男）
- female：出発地・到着地ペア間の大学進学者総数（女）
univ <- openxlsx::read.xlsx(xlsxFile="https://www.dropbox.com/s/xtnv0wcqeaq98o3/mig_2010.xlsx?dl=1")
head(univ)## origin_name destination_name origin destination total national private male
## 1 北海道 北海道 01 01 14906 3612 10499 9225
## 2 北海道 青森県 01 02 379 308 24 246
## 3 北海道 岩手県 01 03 82 52 26 49
## 4 北海道 宮城県 01 04 135 66 65 96
## 5 北海道 秋田県 01 05 49 29 3 31
## 6 北海道 山形県 01 06 74 46 27 36
## female
## 1 5681
## 2 133
## 3 33
## 4 39
## 5 18
## 6 38重力モデルの説明変数として用いる地域属性（プッシュ・プル要因）を読み込みます．変数の定義は以下の通りです．
- pref_code：都道府県ID
- pref_name：都道府県名
- income：1人あたり県民所得
- acom_univ：大学収容率
- 高校卒業の年齢になる子どもたちの人数に対し，地元の大学等の入学枠がどの程度あるのか
- 高校卒業の年齢になる子どもたちの人数に対し，地元の大学等の入学枠がどの程度あるのか
- univ_rate：全学校卒業者に占める大卒者の割合
回帰分析の推定結果を整えるため，県民所得の単位を[千円]から[万円]に変更します．
attriv <- openxlsx::read.xlsx(xlsxFile="https://www.dropbox.com/s/0w49xvdy4ouwr38/table_pref_vars.xlsx?dl=1") %>%
#県民所得の単位を[千円]から[万円]に変更
dplyr::mutate(income=income/10)
head(attriv)## pref_code pref_name income acom_univ univ_rate
## 1 00 全国 294.4 121.1 17.3
## 2 01 北海道 237.9 117.9 11.3
## 3 02 青森県 226.2 68.8 9.1
## 4 03 岩手県 227.5 59.8 9.7
## 5 04 宮城県 238.9 132.2 14.3
## 6 05 秋田県 224.2 53.7 9.0【参考】OD表をlong型に変換
フローデータはOD表（wide型）でのみ利用可能な場合も多く，その際にはデータをlong型（1行1観測値で縦に並ぶ形）に変換する作業が必要です．ここでは，進学者総数のOD表をlong型に変換する方法を示します．この方法は，後程距離行列や隣接判定行列をlong型にする場合にも用います．まず，OD表のExcelファイルを読み込みます．
univ_total_od <- openxlsx::read.xlsx(xlsxFile="https://www.dropbox.com/s/8db68n2fzeon9e9/mig_2010_total_od.xlsx?dl=1")
head(univ_total_od)## origin_name 北海道 青森県 岩手県 宮城県 秋田県 山形県 福島県 茨城県 栃木県
## 1 北海道 14906 379 82 135 49 74 38 93 30
## 2 青森県 408 2003 261 548 87 87 82 84 68
## 3 岩手県 236 289 1357 750 158 111 123 69 87
## 4 宮城県 141 82 314 6414 124 491 262 91 105
## 5 秋田県 158 210 225 536 935 106 103 68 81
## 6 山形県 115 41 74 941 59 937 192 83 75
## 群馬県 埼玉県 千葉県 東京都 神奈川県 新潟県 富山県 石川県 福井県 山梨県
## 1 86 339 302 1991 693 35 13 60 5 57
## 2 61 250 185 666 272 32 6 11 0 25
## 3 137 231 164 628 302 31 5 13 2 60
## 4 91 308 224 1175 474 70 4 26 4 37
## 5 71 199 161 685 305 165 3 20 2 24
## 6 93 316 241 844 381 269 15 25 2 35
## 長野県 岐阜県 静岡県 愛知県 三重県 滋賀県 京都府 大阪府 兵庫県 奈良県
## 1 38 29 85 127 8 82 347 214 167 38
## 2 5 13 21 12 5 3 33 35 13 4
## 3 4 8 19 20 0 7 32 12 3 4
## 4 11 9 27 31 2 26 76 32 21 5
## 5 4 4 22 21 3 2 30 4 3 5
## 6 3 11 26 19 0 5 39 18 5 5
## 和歌山県 鳥取県 島根県 岡山県 広島県 山口県 徳島県 香川県 愛媛県 高知県
## 1 4 6 5 42 51 20 10 8 5 10
## 2 1 1 1 1 1 2 0 0 0 2
## 3 0 0 0 1 1 1 0 1 4 0
## 4 0 0 1 4 6 2 0 0 0 3
## 5 0 1 1 0 1 0 1 1 3 0
## 6 0 1 0 2 0 1 1 1 0 3
## 福岡県 佐賀県 長崎県 熊本県 大分県 宮崎県 鹿児島県 沖縄県
## 1 81 2 17 12 26 9 18 33
## 2 1 0 1 3 6 2 3 0
## 3 1 0 3 1 0 0 1 2
## 4 13 1 3 4 2 2 2 8
## 5 1 0 0 0 3 1 2 1
## 6 1 0 2 4 2 1 3 1tidyrのpivot_longer関数を用いて，データをlong型に変換します．縦に並べる変数（北海道〜沖縄県）を引数colsで，横に並んでいた元の列名（到着地側の都道府県名）が格納される変数の名前を「destination_name」にすることを引数names_toで，縦に並ぶ観測値（進学者総数）が入る変数の名前を引数values_toで指定します．
univ_total_od_long <- univ_total_od %>%
#OD表をlong型に変形
tidyr::pivot_longer(cols=c(北海道:沖縄県),
names_to="destination_name",
values_to="total")
head(univ_total_od_long)## # A tibble: 6 × 3
## origin_name destination_name total
## <chr> <chr> <dbl>
## 1 北海道 北海道 14906
## 2 北海道 青森県 379
## 3 北海道 岩手県 82
## 4 北海道 宮城県 135
## 5 北海道 秋田県 49
## 6 北海道 山形県 74データ構築
県外進学データに属性を結合
地域属性データと位置関係に関するデータを，都道府県コードをキー変数として，県外進学データに結合します．その際地域属性データは，出発地・到着地それぞれに対して結合する必要があります．まず，出発地側に地域属性を結合します．
県外進学データ側では，出発地側の都道府県コードの変数名はoriginで，地域属性データ側では都道府県コードの変数名はpref_codeですので，その2つが同じであることを明示した上で，dplyrのleft_join関数を実行します．
後程，到着地側についても同じ地域属性変数を結合するので，変数名の重複を防ぐため，rename_with関数を用いて結合された変数の変数名の末尾に「_o」という添字を付け加えます．
univ_attriv <- univ %>%
#地域属性データを出発地側に結合
dplyr::left_join(y=attriv,by=c("origin"="pref_code")) %>%
#結合した変数について，末尾に「_o」という文字列を結合
dplyr::rename_with(.fn=~paste0(.,"_o"),.cols=pref_name:univ_rate)
head(univ_attriv)## origin_name destination_name origin destination total national private male
## 1 北海道 北海道 01 01 14906 3612 10499 9225
## 2 北海道 青森県 01 02 379 308 24 246
## 3 北海道 岩手県 01 03 82 52 26 49
## 4 北海道 宮城県 01 04 135 66 65 96
## 5 北海道 秋田県 01 05 49 29 3 31
## 6 北海道 山形県 01 06 74 46 27 36
## female pref_name_o income_o acom_univ_o univ_rate_o
## 1 5681 北海道 237.9 117.9 11.3
## 2 133 北海道 237.9 117.9 11.3
## 3 33 北海道 237.9 117.9 11.3
## 4 39 北海道 237.9 117.9 11.3
## 5 18 北海道 237.9 117.9 11.3
## 6 38 北海道 237.9 117.9 11.3同様にして，到着地側についても地域属性を結合します．その際，結合された変数の変数名の末尾に「_d」を付け加えます．
univ_attriv <- univ_attriv %>%
#地域属性データを到着地側に結合
dplyr::left_join(y=attriv,by=c("destination"="pref_code")) %>%
#結合した変数について，末尾に「_d」という文字列を結合
dplyr::rename_with(.fn=~paste0(.,"_d"),.cols=pref_name:univ_rate)
head(univ_attriv)## origin_name destination_name origin destination total national private male
## 1 北海道 北海道 01 01 14906 3612 10499 9225
## 2 北海道 青森県 01 02 379 308 24 246
## 3 北海道 岩手県 01 03 82 52 26 49
## 4 北海道 宮城県 01 04 135 66 65 96
## 5 北海道 秋田県 01 05 49 29 3 31
## 6 北海道 山形県 01 06 74 46 27 36
## female pref_name_o income_o acom_univ_o univ_rate_o pref_name_d income_d
## 1 5681 北海道 237.9 117.9 11.3 北海道 237.9
## 2 133 北海道 237.9 117.9 11.3 青森県 226.2
## 3 33 北海道 237.9 117.9 11.3 岩手県 227.5
## 4 39 北海道 237.9 117.9 11.3 宮城県 238.9
## 5 18 北海道 237.9 117.9 11.3 秋田県 224.2
## 6 38 北海道 237.9 117.9 11.3 山形県 242.4
## acom_univ_d univ_rate_d
## 1 117.9 11.3
## 2 68.8 9.1
## 3 59.8 9.7
## 4 132.2 14.3
## 5 53.7 9.0
## 6 59.7 10.5都道府県庁座標変数の作成
後々座標値を用いた演算で用いるため，sfのst_coordinates関数で都道府県庁ポイントの座標を取り出し，tibble形式（data.frameの拡張版）に変換したものを作成します．一旦，都道府県庁ポイントpref_office上で座標値変数を作成し，geometry属性を削除します．
pref_office_coords <- pref_office %>%
#ポイントのX座標を変数として追加
dplyr::mutate(x=sf::st_coordinates(.)[,"X"],
#ポイントのY座標を変数として追加
y=sf::st_coordinates(.)[,"Y"]) %>%
#geometry属性を削除
sf::st_drop_geometry()
head(pref_office_coords)## # A tibble: 6 × 3
## pref_code x y
## <chr> <dbl> <dbl>
## 1 01 528328. 4767974.
## 2 02 478083. 4519290.
## 3 03 513081. 4394877.
## 4 04 488790. 4235648.
## 5 05 423139. 4396850.
## 6 06 444280. 4232686.都道府県間隣接ダミーの作成
ポリゴン間の隣接判定は，sfのst_intersects関数で行えます．引数sparseをFALSEに指定することで，隣接有無をTRUE/FALSEで示された行列が返されます．距離行列の場合と同じく，tibble形式に変換した後，dplyrのmutate関数とacross関数を組み合わせることで，全ての変数の値を数値型に変換します．これにより，隣接する場合は要素の値が1，しない場合は0に置き換わります．
#隣接関係行列の取得
pref_poly_nb_mat <- sf::st_intersects(x=pref_poly,sparse=FALSE)
#1-5行目及び1-5列目を表示
pref_poly_nb_mat[1:5,1:5]## [,1] [,2] [,3] [,4] [,5]
## [1,] TRUE FALSE FALSE FALSE FALSE
## [2,] FALSE TRUE TRUE FALSE TRUE
## [3,] FALSE TRUE TRUE TRUE TRUE
## [4,] FALSE FALSE TRUE TRUE TRUE
## [5,] FALSE TRUE TRUE TRUE TRUEas_tibble関数を用いて，距離行列をtibble形式に変換します．その後，dplyrのmutate関数とacross関数を組み合わせることで，全ての変数の値を数値型に変換します．これにより，隣接する場合は要素の値が1，しない場合は0に置き換わります．
pref_poly_nb <- pref_poly_nb_mat %>%
#tibbleに変換
dplyr::as_tibble() %>%
#全ての変数の値を数値型に変換する
#これにより，TRUE/FALSEが1/0に置き換わる
dplyr::mutate(across(everything(),~as.numeric(.)))
#1-5行目及び1-5列目を表示
pref_poly_nb[1:5,1:5]## # A tibble: 5 × 5
## V1 V2 V3 V4 V5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 0 0 0
## 2 0 1 1 0 1
## 3 0 1 1 1 1
## 4 0 0 1 1 1
## 5 0 1 1 1 1都道府県コードの変数pref_code_originを追加します．都道府県庁ポイントは，都道府県コードの順に並んでいましたので，単純に行番号の先頭をゼロ埋めした変数を作ればOKです．値のゼロ埋めは，formatC関数で行えます．引数widthでゼロ埋め後の値の文字数，flagで埋めるのに使う文字列（今回は「0」）を指定します．その後，relocate関数を用いて，この変数pref_code_originを先頭列に移動します．これで，出発地側の都道府県コード変数が作成された状態です．
pref_poly_nb <- pref_poly_nb %>%
#出発地側の都道府県コード変数pref_code_originを作成
dplyr::mutate(pref_code_origin=formatC(x=dplyr::row_number(),width=2,flag="0")) %>%
#pref_code_originを先頭列に移動
dplyr::relocate(pref_code_origin)
pref_poly_nb[1:5,1:5]## # A tibble: 5 × 5
## pref_code_origin V1 V2 V3 V4
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 01 1 0 0 0
## 2 02 0 1 1 0
## 3 03 0 1 1 1
## 4 04 0 0 1 1
## 5 05 0 1 1 1現状，縦横両方に並んだ状態の隣接判定変数を，tidyrのpivot_longer関数を用いて，1行1観測値で縦に並ぶ形（long型）に変換します．縦に並べる変数（V1〜V47）を引数colsで，横に並んでいた元の列名（到着地側の都道府県コードに「V」が付いたものに相当）が格納される変数の名前を「pref_code_destination」にすることを引数names_toで，縦に並ぶ観測値（隣接判定変数）が入る変数の名前を引数values_toで指定します．
最後に，変数pref_code_destinationから文字列「V」を削除し，一旦数値型に変換した上で，先頭をゼロ埋めします．これにより出発地側と到着地側で，都道府県コードのフォーマットが統一され，県外進学データに結合可能な形式になります．
pref_poly_nb <- pref_poly_nb %>%
#行列をlong型に変形
tidyr::pivot_longer(cols=c(V1:V47),
names_to="pref_code_destination",
values_to="contig") %>%
#到着地側の都道府県コードの先頭に付いた「V」を削除
dplyr::mutate(pref_code_destination=gsub("V","",pref_code_destination)) %>%
#到着地側の都道府県コードを数値型に変換した後ゼロ埋め
dplyr::mutate(pref_code_destination=formatC(x=as.numeric(pref_code_destination),width=2,flag="0"))
head(pref_poly_nb)## # A tibble: 6 × 3
## pref_code_origin pref_code_destination contig
## <chr> <chr> <dbl>
## 1 01 01 1
## 2 01 02 0
## 3 01 03 0
## 4 01 04 0
## 5 01 05 0
## 6 01 06 0県外進学データに位置関係変数を結合
ここまでで作成した位置情報変数のうち，まずは隣接判定変数を県外進学データに結合します．left_join関数を用いる際，キー変数が出発地側と到着地側の2つある点に注意して，引数byを指定します．
univ_attriv_loc <- univ_attriv %>%
#隣接判定変数を結合
dplyr::left_join(y=pref_poly_nb,by=c("origin"="pref_code_origin","destination"="pref_code_destination"))続いて，都道府県庁座標を結合します．その際，県外進学データに属性を結合セクションと同様，出発地側と到着地側で座標変数の名前を別個にします．
univ_attriv_loc <- univ_attriv_loc %>%
#県庁座標を出発地側に結合
dplyr::left_join(y=pref_office_coords,by=c("origin"="pref_code")) %>%
#結合した変数について，末尾に「_o」という文字列を結合
dplyr::rename_with(.fn=~paste0(.,"_o"),.cols=x:y) %>%
#県庁座標を到着地側に結合
dplyr::left_join(y=pref_office_coords,by=c("destination"="pref_code")) %>%
#結合した変数について，末尾に「_d」という文字列を結合
dplyr::rename_with(.fn=~paste0(.,"_d"),.cols=x:y)出発地・到着地の座標値は投影座標系のもとで付与されているので，距離公式を用いて，地点間の直線距離を簡単に求めることができます（緯度経度の場合でも，geosphereパッケージのdistHaversine関数等を用いて簡単に計算可能です）．変数distとして，距離変数をデータに追加します．
univ_attriv_loc <- univ_attriv_loc %>%
#距離変数を追加（1/1000することで[km]単位になる）
dplyr::mutate(dist=sqrt((x_o-x_d)^2+(y_o-y_d)^2)/1000)
head(univ_attriv_loc)## origin_name destination_name origin destination total national private male
## 1 北海道 北海道 01 01 14906 3612 10499 9225
## 2 北海道 青森県 01 02 379 308 24 246
## 3 北海道 岩手県 01 03 82 52 26 49
## 4 北海道 宮城県 01 04 135 66 65 96
## 5 北海道 秋田県 01 05 49 29 3 31
## 6 北海道 山形県 01 06 74 46 27 36
## female pref_name_o income_o acom_univ_o univ_rate_o pref_name_d income_d
## 1 5681 北海道 237.9 117.9 11.3 北海道 237.9
## 2 133 北海道 237.9 117.9 11.3 青森県 226.2
## 3 33 北海道 237.9 117.9 11.3 岩手県 227.5
## 4 39 北海道 237.9 117.9 11.3 宮城県 238.9
## 5 18 北海道 237.9 117.9 11.3 秋田県 224.2
## 6 38 北海道 237.9 117.9 11.3 山形県 242.4
## acom_univ_d univ_rate_d contig x_o y_o x_d y_d dist
## 1 117.9 11.3 1 528327.9 4767974 528327.9 4767974 0.0000
## 2 68.8 9.1 0 528327.9 4767974 478083.0 4519290 253.7094
## 3 59.8 9.7 0 528327.9 4767974 513081.3 4394877 373.4085
## 4 132.2 14.3 0 528327.9 4767974 488789.6 4235648 533.7927
## 5 53.7 9.0 0 528327.9 4767974 423139.2 4396850 385.7435
## 6 59.7 10.5 0 528327.9 4767974 444280.3 4232686 541.8458最後に，同都道府県内のフロー（内々フロー）のレコードを落とします．
univ_attriv_loc_inter <- univ_attriv_loc %>%
#出発地の都道府県コードと到着地の都道府県コードが同じレコードを落とす
dplyr::filter(origin!=destination)地域属性データに位置関係変数を結合
後ほど行うネットワーク分析で用いるため，地域属性データにも，都道府県庁座標を結合します．
attriv_loc <- attriv %>%
dplyr::left_join(y=pref_office_coords,by="pref_code") %>%
#座標値が結合できたレコードのみ残す
dplyr::filter(!is.na(x))【参考】st_distance関数を用いた県庁間距離の計算
座標値を直接使って直線距離を計算する方法の他，sfのst_distance関数を用いて地物間の距離を計算する方法もあります．例として，都道府県庁ポイント間の距離行列は以下のようにして得られます．
#都道府県庁ポイント間の距離行列を計算
pref_office_dist_mat <- sf::st_distance(x=pref_office)
#1-5行目及び1-5列目を表示
pref_office_dist_mat[1:5,1:5]## Units: [m]
## 1 2 3 4 5
## 1 0.0 253709.4 373408.55 533792.7 385743.52
## 2 253709.4 0.0 129241.76 283844.0 134202.87
## 3 373408.5 129241.8 0.00 161071.5 89963.75
## 4 533792.7 283844.0 161071.48 0.0 174057.47
## 5 385743.5 134202.9 89963.75 174057.5 0.00距離行列の上側に「Units: [m]」と出ていることから，距離行列は単位付きの数値行列になっていることが分かります．ただし，単位付きだと後々の扱いが面倒なので，unitsのset_units関数で単位を落とします．
pref_office_dist_mat <- pref_office_dist_mat %>%
#単位を落とす
units::set_units(value=NULL)as_tibble関数を用いて，距離行列をtibble形式に変換します．その後，dplyrのmutate関数とacross関数を組み合わせることで，全ての変数の値を1/1000し，距離の単位をkmに変更します．
pref_office_dist <- pref_office_dist_mat %>%
#tibbleに変換
dplyr::as_tibble() %>%
#全ての変数の値を1000で割る
#これにより，距離の単位を[m]から[km]に変更
dplyr::mutate(across(everything(),~./1000))距離行列に，都道府県コードの変数pref_code_originを追加します．都道府県庁ポイントは，都道府県コードの順に並んでいましたので，単純に行番号の先頭をゼロ埋めした変数を作ればOKです．その後，relocate関数を用いて，この変数pref_code_originを先頭列に移動します．これで，出発地側の都道府県コード変数が作成された状態です．
pref_office_dist <- pref_office_dist %>%
#出発地側の都道府県コード変数pref_code_originを作成
dplyr::mutate(pref_code_origin=formatC(row_number(),width=2,flag="0")) %>%
#pref_code_originを先頭列に移動
dplyr::relocate(pref_code_origin)現状，縦横両方に並んだ状態の県庁間距離を，tidyrのpivot_longer関数を用いて，1行1観測値で縦に並ぶ形（long型）に変換します．縦に並べる変数（1〜47）を引数colsで，横に並んでいた元の列名（到着地側の都道府県コードに相当）が格納される変数の名前を「pref_code_destination」にすることを引数names_toで，縦に並ぶ観測値（県庁間距離）が入る変数の名前を引数values_toで指定します．
最後に，変数pref_code_destinationを一旦数値型に変換した上で，先頭をゼロ埋めします．これにより出発地側と到着地側で，都道府県コードのフォーマットが統一され，県外進学データに結合可能な形式になります．
pref_office_dist <- pref_office_dist %>%
#距離行列をlong型に変形
tidyr::pivot_longer(cols=c(`1`:`47`),
names_to="pref_code_destination",
values_to="dist") %>%
#到着地側の都道府県コードを数値型に変換した後ゼロ埋め
dplyr::mutate(pref_code_destination=formatC(x=as.numeric(pref_code_destination),width=2,flag="0"))
head(pref_office_dist)## # A tibble: 6 × 3
## pref_code_origin pref_code_destination dist
## <chr> <chr> <dbl>
## 1 01 01 0
## 2 01 02 254.
## 3 01 03 373.
## 4 01 04 534.
## 5 01 05 386.
## 6 01 06 542.記述的分析
重力モデルでの分析を行う前に，ggplot2の関数群を用いて，幾つかの可視化を行います．
前処理
今後の可視化で用いやすいように，各種進学者数変数について，データをlong型に変形します．
univ_attriv_loc_inter_long <- univ_attriv_loc_inter %>%
#進学者数変数と距離変数だけを残す
dplyr::select(total,national,private,male,female,dist) %>%
#進学者数変数についてデータをlong型に変形
tidyr::pivot_longer(cols=-dist)
head(univ_attriv_loc_inter_long)## # A tibble: 6 × 3
## dist name value
## <dbl> <chr> <dbl>
## 1 254. total 379
## 2 254. national 308
## 3 254. private 24
## 4 254. male 246
## 5 254. female 133
## 6 373. total 82更に，ggplot2の各種関数を適用しやすいように，進学者数変数の値だけレコードを複製します．これにより，データは擬似的に非集計な状態になります．データの複製にはtidyrのuncount関数を用います．この関数では，引数weightsに指定した変数に入った数字の数だけ，レコードを複製できます．今回は，変数value（各カテゴリの進学者数）に入った件数分，各レコードを複製します．
univ_attriv_loc_inter_long_uc <- univ_attriv_loc_inter_long %>%
#進学者数valueの値だけ各レコードを複製する
tidyr::uncount(weights=value)
head(univ_attriv_loc_inter_long_uc)## # A tibble: 6 × 2
## dist name
## <dbl> <chr>
## 1 254. total
## 2 254. total
## 3 254. total
## 4 254. total
## 5 254. total
## 6 254. total距離帯別の進学者数
密度プロットを用いて，距離帯別の進学者のシェアをプロットします．まず，非集計化されたデータから，進学者総数に該当するレコードのみを取り出します．その上で，ggplot2の関数群を用いた可視化を行います．最初にggplot関数で白地のキャンバスを作り，その上に各種可視化のための関数を「+」で重ねていくという文法は，ggplot2を用いた可視化の殆どに共通します．
geom_density関数の引数dataには可視化に用いるデータを指定します．引数mappingに指定されるaesの内部では，どの変数をx,y軸に取るか，色分けの基準に使う変数は何か等を指定します．
univ_attriv_loc_inter_long_uc_total <- univ_attriv_loc_inter_long_uc %>%
#進学者総数に該当するレコードを残す
dplyr::filter(name=="total")
ggplot2::ggplot() +
#密度プロットを作成
ggplot2::geom_density(data=univ_attriv_loc_inter_long_uc_total,mapping=aes(x=dist))
次に，国立・私立別の密度プロットを重ねて描画します．非集計化されたデータから，国立・私立別の進学者数に該当するレコードを取り出し，geom_density関数を用いて密度プロットを可視化します．その際，aesの中で「fill=name」と指定しますが，これはプロットの色分けを，変数name（今回の例ではnational/private）に基づいて行うことを意味します．
私立大学への県外進学は地理的に極めて近い都道府県間で行われる一方，国立大学への県外進学は比較的離れた地域間でも行われる傾向が見て取れます．
univ_attriv_loc_inter_long_uc_founder <- univ_attriv_loc_inter_long_uc %>%
#国立・私立別の進学者数に該当するレコードを残す
dplyr::filter(name=="national"|name=="private")
ggplot2::ggplot() +
ggplot2::geom_density(data=univ_attriv_loc_inter_long_uc_founder,
#x軸方向にdistを取り，色分けはnameに基づき行う
aes(x=dist,fill=name),
#不透明度を0.4に設定
alpha=0.4,
position="identity",
#枠線は描画しない
lty=0)
最後に，男女別の密度プロットを重ねて描画します．男女別の進学者数に該当するレコードを取り出したあとは，国立・私立別の可視化の場合と同じです．国立・私立別のプロット程の差は見られませんが，女子の県外進学の方がよりローカルである傾向が見られます．
univ_attriv_loc_inter_long_uc_gender <- univ_attriv_loc_inter_long_uc %>%
#男女別の進学者数に該当するレコードを残す
dplyr::filter(name=="male"|name=="female")
ggplot2::ggplot() +
ggplot2::geom_density(data=univ_attriv_loc_inter_long_uc_gender,
#x軸方向にdistを取り，色分けはnameに基づき行う
aes(x=dist,fill=name),
#不透明度を0.4に設定
alpha=0.4,
position="identity",
#枠線は描画しない
lty=0)
重力モデル推定
モデルの定式化
今回推定する重力モデルは，都道府県\(i\)から\(j\)への進学者数を\(flow_{i,j}\)，プッシュ要因の地域属性をまとめたベクトルを\(\mathbf{push}_{i}\)，プル要因のベクトルを\(\mathbf{pull}_{j}\)として，以下のように定式化されます．
アウトカムに対して1を足した上で対数を取る理由は，進学者数が0となる都道府県ペアが存在する一方，0の対数値は定義できないためです．
\[ \ln\left(flow_{i,j}+1\right)=\alpha+\mathbf{push}_{i}'\boldsymbol{\beta}+\mathbf{pull}_{j}'\boldsymbol{\gamma}+\tau_{1}\ln dist_{i,j}+\tau_{2}contig_{i,j}+\epsilon_{i,j}, \\ \mathbf{push}_{i}'=\begin{bmatrix} income_{i} & acomUniv_{i} & univRate_{i} \end{bmatrix}, \mathbf{pull}_{j}'=\begin{bmatrix} income_{j} & acomUniv_{j} & univRate_{j} \end{bmatrix} \]
進学者総数
進学者総数をアウトカムとする重力モデルをOLS推定します．
res_total <- estimatr::lm_robust(log(total+1)~income_o+income_d+acom_univ_o+acom_univ_d+univ_rate_o+univ_rate_d+log(dist)+contig,
data=univ_attriv_loc_inter,
se_type="stata")
modelsummary::modelsummary(models=res_total,
statistic="statistic")| Model 1 | |
|---|---|
| (Intercept) | 7.184 |
| (20.795) | |
| income_o | −0.0007 |
| (−1.112) | |
| income_d | 0.002 |
| (2.743) | |
| acom_univ_o | 0.005 |
| (7.439) | |
| acom_univ_d | 0.017 |
| (21.279) | |
| univ_rate_o | −0.027 |
| (−2.991) | |
| univ_rate_d | 0.047 |
| (5.040) | |
| log(dist) | −1.101 |
| (−28.859) | |
| contig | 1.050 |
| (12.159) | |
| Num.Obs. | 2162 |
| R2 | 0.673 |
| R2 Adj. | 0.672 |
| AIC | 34993.1 |
| BIC | 35049.9 |
| RMSE | 787.76 |
国立進学者数
同じく，国立大学進学者数をアウトカムとする重力モデルをOLS推定します．
res_national <- estimatr::lm_robust(log(national+1)~income_o+income_d+acom_univ_o+acom_univ_d+univ_rate_o+univ_rate_d+log(dist)+contig,
data=univ_attriv_loc_inter,
se_type="stata")
modelsummary::modelsummary(models=res_national,
statistic="statistic")| Model 1 | |
|---|---|
| (Intercept) | 7.278 |
| (19.731) | |
| income_o | −0.0005 |
| (−0.853) | |
| income_d | 0.001 |
| (2.267) | |
| acom_univ_o | 0.004 |
| (5.238) | |
| acom_univ_d | 0.012 |
| (15.189) | |
| univ_rate_o | −0.008 |
| (−0.864) | |
| univ_rate_d | −0.034 |
| (−3.911) | |
| log(dist) | −1.014 |
| (−25.460) | |
| contig | 0.840 |
| (8.674) | |
| Num.Obs. | 2162 |
| R2 | 0.518 |
| R2 Adj. | 0.517 |
| AIC | 24569.4 |
| BIC | 24626.2 |
| RMSE | 70.71 |
私立大学進学者数
同じく，私立大学進学者数をアウトカムとする重力モデルをOLS推定します．
res_private <- estimatr::lm_robust(log(private+1)~income_o+income_d+acom_univ_o+acom_univ_d+univ_rate_o+univ_rate_d+log(dist)+contig,
data=univ_attriv_loc_inter,
se_type="stata")
modelsummary::modelsummary(models=res_private,
statistic="statistic")| Model 1 | |
|---|---|
| (Intercept) | 3.577 |
| (9.024) | |
| income_o | −0.0003 |
| (−0.375) | |
| income_d | 0.002 |
| (3.778) | |
| acom_univ_o | 0.005 |
| (6.880) | |
| acom_univ_d | 0.018 |
| (22.618) | |
| univ_rate_o | −0.056 |
| (−5.683) | |
| univ_rate_d | 0.129 |
| (13.297) | |
| log(dist) | −0.863 |
| (−20.356) | |
| contig | 1.339 |
| (10.676) | |
| Num.Obs. | 2162 |
| R2 | 0.688 |
| R2 Adj. | 0.687 |
| AIC | 34644.8 |
| BIC | 34701.6 |
| RMSE | 726.80 |
結果をまとめて表示
modelsummary関数を用いれば，全ての結果をまとめて表示することも可能です．以下雑駁ですが，各説明変数について考察を行います．
- 1人あたり県民所得income：出発地側では負だが有意でない．到着地側では正で有意．全モデルで結果は共通．出身県の所得水準によらず，所得水準が高い都道府県の大学に進学する傾向がありそう．
- 大学収容率acom_univ：出発地・到着地両方で正に有意．全モデルで結果は共通．県外進学は多くの大学を擁する地域間で行われる傾向がありそう．
- 大卒者割合univ_rate：進学者総数（Model 1）及び私立進学者数（Model
3）では，出発地側では負で有意，到着地側では正に有意．国立進学者数（Model
2）では，到着地側でのみ負に有意．国立進学者数に関する結果は，（大学進学率が低めの）地方部の国立大学が，都市部からも入学者を集める傾向を反映していそうだが，あまり直感的ではない．
- 都道府県庁間距離dist：全モデルで共通して負で有意．距離的に近い県への進学を行う傾向がありそう．
- 隣接ダミーcontig：全モデルで共通して正で有意．隣接都道府県への進学を行う傾向がありそう．
修正済み決定係数の大きさは，私立＞総数＞国立の順です．国立大学への県外進学の傾向は，相対的に（少なくとも今回の定式化に基づく）重力モデルでは説明がしづらいようです．
modelsummary::modelsummary(models=list(res_total,res_national,res_private),
statistic="statistic")| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 7.184 | 7.278 | 3.577 |
| (20.795) | (19.731) | (9.024) | |
| income_o | −0.0007 | −0.0005 | −0.0003 |
| (−1.112) | (−0.853) | (−0.375) | |
| income_d | 0.002 | 0.001 | 0.002 |
| (2.743) | (2.267) | (3.778) | |
| acom_univ_o | 0.005 | 0.004 | 0.005 |
| (7.439) | (5.238) | (6.880) | |
| acom_univ_d | 0.017 | 0.012 | 0.018 |
| (21.279) | (15.189) | (22.618) | |
| univ_rate_o | −0.027 | −0.008 | −0.056 |
| (−2.991) | (−0.864) | (−5.683) | |
| univ_rate_d | 0.047 | −0.034 | 0.129 |
| (5.040) | (−3.911) | (13.297) | |
| log(dist) | −1.101 | −1.014 | −0.863 |
| (−28.859) | (−25.460) | (−20.356) | |
| contig | 1.050 | 0.840 | 1.339 |
| (12.159) | (8.674) | (10.676) | |
| Num.Obs. | 2162 | 2162 | 2162 |
| R2 | 0.673 | 0.518 | 0.688 |
| R2 Adj. | 0.672 | 0.517 | 0.687 |
| AIC | 34993.1 | 24569.4 | 34644.8 |
| BIC | 35049.9 | 24626.2 | 34701.6 |
| RMSE | 787.76 | 70.71 | 726.80 |
実践問題
県外進学の傾向は，男女間でも異なる可能性があります．アウトカムを進学者総数（男）及び進学者総数（女）にした場合の重力モデルを推定し，その結果を比較してください．また，男女間で結果が異なる地域属性があれば，その理由を考察してください．
【参考】フローが多数のゼロ値を含む場合の対処
重力モデルのOLS推定を上で行った際，進学者数（フロー）が0の都道府県ペアへの対処として，アウトカムに1を足して対数を取るという操作をしました．一方こうした操作は，パラメータ推定に際してバイアスをもたらすことが指摘されています（e.g.,
Mullahy & Norton
2022）．加えて，進学者数というアウトカムは非負の整数値なので，データ生成過程との整合性の観点から，カウントデータの回帰モデルを用いた方が妥当そうに思えます．
一般によく用いられるカウントデータの回帰モデルは，ポアソン回帰と負の二項回帰です．ポアソン回帰の場合，アウトカムとなる確率変数の分散が期待値より大きい際に，回帰係数の標準誤差の推定にバイアスがかかることが知られています（ポアソン分布の性質が理由です）．故に，アウトカムの過分散（分散>期待値）が生じている場合は，負の二項回帰を用いる方が望ましいです．
負の二項回帰で重力モデルを推定するに先立ち，ここまでで用いてきた各アウトカムについて0の割合を確認します．進学者総数については0の割合は少ないですが，国立・私立進学者数の場合は0の割合が1割を超えています．
univ_attriv_loc_inter %>%
#total（進学者総数）が0のレコードのみ残す
dplyr::filter(total==0) %>%
#0のレコード数を元のサンプルサイズで割る
nrow()/nrow(univ_attriv_loc_inter)## [1] 0.0453284univ_attriv_loc_inter %>%
#national（国立進学者数）が0のレコードのみ残す
dplyr::filter(national==0) %>%
nrow()/nrow(univ_attriv_loc_inter)## [1] 0.120259univ_attriv_loc_inter %>%
#private（私立進学者数）が0のレコードのみ残す
dplyr::filter(private==0) %>%
nrow()/nrow(univ_attriv_loc_inter)## [1] 0.2271045また，各アウトカムについて，平均と分散を計算すると，どのアウトカムについても分散が平均よりもはるかに大きく，やはり負の二項回帰を用いる方が妥当そうです．
univ_attriv_loc_inter %>%
#アウトカムについて平均を計算
dplyr::summarise(across(total:private,mean))## total national private
## 1 156.6707 29.59667 120.1933univ_attriv_loc_inter %>%
#アウトカムについて分散を計算
dplyr::summarise(across(total:private,var))## total national private
## 1 598380.8 4343.728 515629.5負の二項回帰に基づく重力モデルを推定します．負の二項回帰は，MASSパッケージのglm.nb関数で行えます．各進学者変数をアウトカムに用いた分析を行った後，modelsummary関数を用いて結果をまとめて表示します．今回の推定では，OLSでの推定を行った場合と傾向は変わらず，結果の頑健性が確認できました．
#進学者総数がアウトカム
res_negbin_total <- MASS::glm.nb(total~income_o+income_d+acom_univ_o+acom_univ_d+univ_rate_o+univ_rate_d+log(dist)+contig,
data=univ_attriv_loc_inter)
#国立進学者がアウトカム
res_negbin_national <- MASS::glm.nb(national~income_o+income_d+acom_univ_o+acom_univ_d+univ_rate_o+univ_rate_d+log(dist)+contig,
data=univ_attriv_loc_inter)
#私立進学者がアウトカム
res_negbin_private <- MASS::glm.nb(private~income_o+income_d+acom_univ_o+acom_univ_d+univ_rate_o+univ_rate_d+log(dist)+contig,
data=univ_attriv_loc_inter)
#結果をまとめて表示
modelsummary::modelsummary(models=list(res_negbin_total,res_negbin_national,res_negbin_private),
statistic="statistic")| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 7.646 | 8.066 | 4.899 |
| (22.724) | (21.368) | (11.098) | |
| income_o | −0.0004 | −0.0006 | −0.0004 |
| (−0.727) | (−0.915) | (−0.459) | |
| income_d | 0.0008 | 0.002 | 0.002 |
| (1.344) | (2.985) | (2.940) | |
| acom_univ_o | 0.007 | 0.006 | 0.009 |
| (10.678) | (7.159) | (10.390) | |
| acom_univ_d | 0.019 | 0.014 | 0.022 |
| (27.338) | (17.728) | (23.930) | |
| univ_rate_o | −0.055 | −0.017 | −0.084 |
| (−6.027) | (−1.622) | (−6.919) | |
| univ_rate_d | 0.043 | −0.063 | 0.121 |
| (4.783) | (−6.152) | (10.113) | |
| log(dist) | −1.060 | −1.057 | −1.013 |
| (−29.293) | (−25.937) | (−21.237) | |
| contig | 0.949 | 0.722 | 1.257 |
| (9.653) | (6.574) | (9.800) | |
| Num.Obs. | 2162 | 2162 | 2162 |
| AIC | 20198.8 | 16222.2 | 16306.5 |
| BIC | 20255.6 | 16279.0 | 16363.3 |
| Log.Lik. | −10089.400 | −8101.110 | −8143.265 |
| RMSE | 2958.48 | 273.48 | 4496.22 |
ネットワーク分析
以下のセクションでは，都道府県ペア間の県外進学数を重み付き有向グラフと見た上で，中心性のようなネットワーク分析の指標を計算したり，ネットワークを可視化する方法を示します．今回は例として，九州7県内の県外進学フローからネットワークを作成します．
ネットワークデータへの変換
九州内の県ペア間の県外進学者数を，tidygraphオブジェクトというネットワークデータの形式に変換します．ネットワークへの変換に先立ち，レコードを九州内のものに絞った上で，内々フローの観測値を除きます．その上で，tidygraphのas_tbl_graph関数を用いて，tidygraph形式に変換します．有向グラフとして扱う場合には，引数directedをTRUEに指定します．
作成されたネットワークの中身は，2つの要素で成っています．1つはNode
Dataで，ネットワークの頂点属性を示す要素で，今回の例では都道府県にあたります．もう1つはEdge
Dataで，ネットワークの辺属性を示す指標で，県外進学フローに対応します．
univ_q <- univ %>%
#出発地が九州内のレコードのみ残す
dplyr::filter(origin_name%in%c("福岡県","佐賀県","長崎県","熊本県","大分県","宮崎県","鹿児島県")) %>%
#到着地が九州内のレコードのみ残す
dplyr::filter(destination_name%in%c("福岡県","佐賀県","長崎県","熊本県","大分県","宮崎県","鹿児島県")) %>%
#内々フローを除く
dplyr::filter(origin!=destination)
#ネットワークに変換
univ_q_nw <- tidygraph::as_tbl_graph(x=univ_q,
#有向グラフ
directed=TRUE)
#ネットワークの中身
univ_q_nw## # A tbl_graph: 7 nodes and 42 edges
## #
## # A directed simple graph with 1 component
## #
## # Node Data: 7 × 1 (active)
## name
## <chr>
## 1 福岡県
## 2 佐賀県
## 3 長崎県
## 4 熊本県
## 5 大分県
## 6 宮崎県
## # … with 1 more row
## #
## # Edge Data: 42 × 9
## from to origin destination total national private male female
## <int> <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2 40 41 592 494 98 329 263
## 2 1 3 40 42 456 342 75 284 172
## 3 1 4 40 43 812 413 369 448 364
## # … with 39 more rows頂点属性の結合
上で作成した都道府県属性attriv_locを，頂点に属性として結合します．まずtidygraphのactivate関数において，引数whatを「nodes」に指定することで，頂点属性への処理を行うことを明示します．その上で，left_join関数を用いて，都道府県名をキーに都道府県属性を結合します．
univ_q_nw <- univ_q_nw %>%
#頂点への処理を行うことを明示
tidygraph::activate(what="nodes") %>%
#地域属性を都道府県名をキーに結合
tidygraph::left_join(y=attriv_loc,by=c("name"="pref_name"))中心性の計算
ネットワークから，各頂点の次数中心性（重みを考慮）を計算します．次数中心性は，多くの辺が出入りする頂点を中心性が高いと評価する中心性指標です．有向グラフの場合，頂点に入る辺の数（入次数）と頂点から出る辺の数（出次数）が異なります．
まず，activate関数で頂点への処理を行うことを明示した上で，tidygraphのcentrality_degree関数を用いて，入次数と出次数を，頂点属性に新たな変数として追加します．引数weightsには辺の重みとして使う変数を指定します．今回は，進学者総数totalを重みとして用います．引数modeで「in」を指定すると入次数，「out」を指定すると出次数が計算されます．
次数中心性の他によく用いられる指標に，PageRankと呼ばれるものがあります．PageRankは，次数中心性が高い頂点に隣接する頂点も，また中心性が高いと評価する指標です．
univ_q_nw_cent <- univ_q_nw %>%
#頂点への処理を行うことを明示
tidygraph::activate(what="nodes") %>%
#次数中心性（入次数）を計算
dplyr::mutate(deg_i=tidygraph::centrality_degree(weights=total,mode="in"),
#次数中心性（出次数）を計算
deg_o=tidygraph::centrality_degree(weights=total,mode="out")) %>%
#PageRankを計算
dplyr::mutate(PR=tidygraph::centrality_pagerank(weights=total))
univ_q_nw_cent## # A tbl_graph: 7 nodes and 42 edges
## #
## # A directed simple graph with 1 component
## #
## # Node Data: 7 × 10 (active)
## name pref_code income acom_univ univ_rate x y deg_i deg_o PR
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 福岡県 40 268. 137. 16.3 -483866. 3.77e6 7099 2722 0.350
## 2 佐賀県 41 242. 54 12.3 -499158. 3.73e6 1078 2179 0.110
## 3 長崎県 42 232. 72.5 11 -545097. 3.68e6 1289 2164 0.111
## 4 熊本県 43 225. 102. 12.4 -462730. 3.68e6 2324 2014 0.171
## 5 大分県 44 245 76.7 12.7 -376136. 3.72e6 1050 1739 0.0959
## 6 宮崎県 45 212. 59.3 10.4 -407259. 3.57e6 847 1785 0.0743
## # … with 1 more row
## #
## # Edge Data: 42 × 9
## from to origin destination total national private male female
## <int> <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2 40 41 592 494 98 329 263
## 2 1 3 40 42 456 342 75 284 172
## 3 1 4 40 43 812 413 369 448 364
## # … with 39 more rowsネットワークの可視化
シンプルな可視化
作成されたネットワークを可視化します．可視化には，ggraphの関数群を用います．ggplot2の関数を用いる場合と同じく，最初にどのネットワークを可視化するか等を指定した上で，辺や頂点を可視化する際のオプションを指定します．
辺のオプションはgeom_edge_link関数，頂点のオプションはgeom_node_point関数で指定できますが，最初は特に何も指定せずに可視化します．
#乱数の指定（可視化の度にレイアウトが変わるので）
set.seed(1)
#可視化するネットワークとレイアウトの設定
ggraph::ggraph(graph=univ_q_nw_cent,layout="fr") +
#辺を可視化
ggraph::geom_edge_link() +
#頂点を可視化
ggraph::geom_node_point()
辺や頂点の色を（ある一色に）変えたい時には，引数colorで色を指定します．
set.seed(1)
ggraph::ggraph(graph=univ_q_nw_cent,layout="fr") +
#辺を灰色で可視化
ggraph::geom_edge_link(color="gray") +
#頂点を青色で可視化
ggraph::geom_node_point(color="blue")
中心性指標に基づく可視化
例として，入次数（deg_i）の高さに応じて，頂点の大きさ・色が変わるネットワークの可視化を行います．まず，頂点にアクセスした上で，入次数が低い順に並べ替えます（こうすることで，入次数が高い頂点がより上にプロットされます）．
その上で，ネットワークを可視化します．その際，geom_node_point関数の中で，aesで囲んだ上で，頂点の大きさと色が次数中心性の値に応じて変わるような指定を行います．
univ_q_nw_cent <- univ_q_nw_cent %>%
tidygraph::activate("nodes") %>%
#次数中心性が低い順に並び替え
tidygraph::arrange(deg_i)
ggraph::ggraph(graph=univ_q_nw_cent,layout="fr") +
ggraph::geom_edge_link(color="gray") +
#頂点の大きさ・色は次数中心性に応じて変化
ggraph::geom_node_point(aes(size=deg_i,color=deg_i)) +
#背景を透明にする
ggraph::theme_graph(background="transparent")
その他細かい設定として，scale_color_viridis関数で色分けのパレットを変えたり，geom_node_text関数で頂点にラベルを表示したりできます．
ggraph::ggraph(graph=univ_q_nw_cent,layout="fr") +
ggraph::geom_edge_link(color="gray") +
ggraph::geom_node_point(aes(size=deg_i,color=deg_i)) +
#頂点の色のパレットの設定（赤→黄）
ggraph::scale_color_viridis(option="C",begin=0.5,end=1) +
#県名nameをラベルに使用
ggraph::geom_node_text(aes(label=name),
#Macの場合に文字化け防止
family="HiraKakuPro-W3",
#文字色は青
color="blue",
size=3,
check_overlap=TRUE) +
ggraph::theme_graph(background="transparent")
位置座標を組み合わせた可視化
頂点の位置座標を手動で指定したネットワークの可視化を行うことも可能です．例としてここでは，県庁の位置座標で頂点位置を指定した可視化を行います．まず，頂点属性にアクセスし，位置座標に該当する変数を取り出します．
univ_q_nw_cent_coords <- univ_q_nw_cent %>%
tidygraph::activate(nodes) %>%
#tibble形式に変換
dplyr::as_tibble() %>%
#位置座標変数x,yを取り出し
dplyr::select(x,y)まず，背景に何も表示せず，県庁座標で頂点位置を指定した可視化を行います．その際，関数ggraphの引数layoutを，上で取り出した座標で指定します．その他設定は，上の可視化と同じです．
ggraph::ggraph(graph=univ_q_nw_cent,layout=univ_q_nw_cent_coords) +
ggraph::geom_edge_link(color="gray") +
ggraph::geom_node_point(aes(size=deg_i,color=deg_i)) +
#頂点の色のパレットの設定（赤→黄）
ggraph::scale_color_viridis(option="C",begin=0.5,end=1) +
#県名nameをラベルに使用
ggraph::geom_node_text(aes(label=name),
#Macの場合に文字化け防止
family="HiraKakuPro-W3",
#文字色は青
color="blue",
size=3,
check_overlap=TRUE) +
ggraph::theme_graph(background="transparent")
次に，九州のポリゴンを背景としてネットワークを可視化します．冒頭で読み込んだ都道府県ポリゴンから，九州内の県のポリゴンのみを残します．ggplot2のgeom_sf関数を用いて，最初に九州ポリゴンを描画し，その上に先ほどのネットワークを描画します．進学者の流入数の面では福岡県が最も優勢で，その次に熊本県が続くという構造が見られます．
#九州の都道府県ポリゴンのみ残す
poly_q <- pref_poly %>%
dplyr::filter(N03_001%in%c("福岡県","佐賀県","長崎県","熊本県","大分県","宮崎県","鹿児島県"))
ggraph::ggraph(graph=univ_q_nw_cent,layout=univ_q_nw_cent_coords) +
#都道府県ポリゴンを枠線無し・黄緑で可視化
ggplot2::geom_sf(data=poly_q,fill="greenyellow",lty=0) +
ggraph::geom_edge_link(color="gray") +
ggraph::geom_node_point(aes(size=deg_i,color=deg_i)) +
#頂点の色のパレットの設定（赤→黄）
ggraph::scale_color_viridis(option="C",begin=0.5,end=1) +
#県名nameをラベルに使用
ggraph::geom_node_text(aes(label=name),
#Macの場合に文字化け防止
family="HiraKakuPro-W3",
#文字色は青
color="blue",
size=3,
check_overlap=TRUE) +
ggraph::theme_graph(background="transparent")
同じ手続きで，今度は出次数が大きい頂点を大きく・明るい色で可視化します．入次数の際には福岡県のみが際立って大きかったのが，出次数の場合は北部九州の県と鹿児島県も頂点が大きくなっています．
ggraph::ggraph(graph=univ_q_nw_cent,layout=univ_q_nw_cent_coords) +
#都道府県ポリゴンを枠線無し・黄緑で可視化
ggplot2::geom_sf(data=poly_q,fill="greenyellow",lty=0) +
ggraph::geom_edge_link(color="gray") +
ggraph::geom_node_point(aes(size=deg_o,color=deg_o)) +
#頂点の色のパレットの設定（赤→黄）
ggraph::scale_color_viridis(option="C",begin=0.5,end=1) +
#県名nameをラベルに使用
ggraph::geom_node_text(aes(label=name),
#Macの場合に文字化け防止
family="HiraKakuPro-W3",
#文字色は青
color="blue",
size=3,
check_overlap=TRUE) +
ggraph::theme_graph(background="transparent")