12: 重力モデル

作者

田村一軌

公開

2024年7月3日

更新日

2024年8月7日

はじめに

今回の演習では、重力モデルについて学びます。 重力モデルは、人々の移動の発生頻度と距離との関係を分析する、最もプリミティブなモデルの1つです。

今回は、佐賀市と久留米市への通勤者数のデータを使って、通勤者数と距離と関係を分析することで、両市の特徴について議論します。

重力モデル

地域\(i\)から地域\(j\)への通勤者数は、

  1. 地域\(i\)に就業者（仕事に就いている人）が多いほど、
  2. 地域\(j\)の従業者（そこで働いている人）が多いほど、
  3. 地域\(i\)と地域\(j\)が近いほど

多くなりそうです。

1と2は地域の大きさ、3は地域間の距離が、それぞれ影響していることを示唆しています。

ところで、大きさと距離でその間の関係性が決まる物理法則に、万有引力の法則があります。

万有引力の法則

質量\(M\)の物体と質量\(m\)の物体の間に働く万有引力の大きさ\(F\)は、物体間の距離が\(d\)のとき、 \[ F=G\frac{Mm}{d^2} \] です（\(G\)は万有引力定数）。

質量≒大きさと距離が重要な要素になっています。 これからヒントを得て、地域間の社会移動を説明する重力モデルが着想されます。

（古典的）重力モデル

地域\(i\)のから地域\(j\)への移動量\(t_{ij}\)は、それぞれの地域の魅力度を\(R_i\)、\(S_j\)、地域間の距離を\(d_{ij}\)としたとき、 \[ t_{ij}=k\frac{R_iS_j}{d_{ij}^2} \] と推定されます（\(k\)は比例定数）。

（古典的）重力モデルでは、地域間の移動量は、地域の魅力\(R_i\)と\(S_j\)に正比例し、距離\(d_{ij}\)の2乗に反比例しています。 しかし、正比例・2乗に反比例というのは、万有引力の法則に倣ったに過ぎません。 実際の移動に当てはめるには、地域の魅力や距離が移動量に与える影響は、移動の種類などによって異なると考えた方が良さそうです。 そこで、(古典的）重力モデルを拡張した（修正）重力モデルが提案されました。

（修正）重力モデル

地域\(i\)のから地域\(j\)への移動量\(t_{ij}\)は、それぞれの地域の魅力度を\(R_i\)、\(S_j\)、地域間の距離を\(d_{ij}\)としたとき、 \[ t_{ij}=k\frac{R_i^\alpha S_j^\beta}{d_{ij}^\gamma} \] と推定されます（\(k\)は比例定数、\(\alpha\)、\(\beta\)、\(\gamma\)はパラメータ）。

現在では「重力モデル」というと、（修正）重力モデルのことを指すことが多いです。

重力モデルは、万有引力の法則に倣ったとても単純なモデルですが、現実の移動によく当てはまることが、多くの実証研究によって示されています。

距離・移動モデル

今回の演習では、重力モデルを簡略化した以下のモデルを考えます。

地域\(i\)から地域\(j\)への通勤者数\(t_{ij}\)は、地域\(i\)の就業者数を\(R_i\)、地域\(j\)の従業者数を\(S_j\)、地域間の距離を\(d_{ij}\)とするとき、 \[ t_{ij}=k\frac{R_iS_j}{d_{ij}^\alpha} \] に従うものとします。さらに、次のように変形します。 \[ \begin{align} \frac{t_{ij}}{R_i}&=k\frac{S_j}{d_{ij}^\alpha}\\ f_{ij}&=\frac{\beta}{d_{ij}^\alpha} \end{align} \] ここで、\(f_{ij}\)は地域\(i\)の就業者のうち地域\(j\)への通勤者の比率を表しています。 このようなモデルを移動・距離モデルと呼ぶことがあります。

両辺の対数をとると、 \[ \log(f_{ij})=\log(\beta)-\alpha\log(d_{ij}) \] となります。\(\log(f_{ij})=Y\)、\(\log(\beta)=b\)、\(-\alpha=a\)、\(\log(d_{ij})=X\)とおくと、上式は \[ Y=aX+b \] という一次式であることがわかります。つまり、線形（単）回帰分析によってパラメータ\(a\)、\(b\)を求めることができます。元のパラメータは、\(\alpha=-a\)、\(\beta=\exp(b)\)として求めることができます。

準備

パッケージ

今回の演習では、以下のパッケージを使用します。 library関数を使ってロードしてください。

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(jpstat)
「このサービスは、政府統計総合窓口(e-Stat)のAPI機能を使用していますが、サービスの内容は国によって保証されたものではありません。」

データの入手

空間データの入手

国土数値情報から、 福岡県、佐賀県、長崎県、熊本県、大分県の「行政区域データ」と「市町村役場等及び公的集会施設データ」 をダウンロードしてください。

合計10個のZIPファイルのダウンロードが完了したら、 ダウンロードフォルダで展開し、できた10個のフォルダを全てプロジェクトのdataフォルダに移動してください。

空間データの読み込み

空間データを読み込みます。

まず、福岡県、佐賀県、長崎県、熊本県、大分県の行政区域データを読み込みます。 これは、塗り分け地図を作成するために使います。

5県のデータを、

fukuoka <- st_read("data/N03-20240101_40_GML/N03-20240101_40.geojson")

などとして、これを5回繰り返してもいいですが、今日は少し違った方法でファイルを読み込んでみましょう。

まず、読み込むGeoJSONファイルのリストを作ります。

reg_files <- str_c("data/N03-20240101_", seq(40, 44), 
                   "_GML/N03-20240101_", seq(40, 44), ".geojson")
reg_files
[1] "data/N03-20240101_40_GML/N03-20240101_40.geojson"
[2] "data/N03-20240101_41_GML/N03-20240101_41.geojson"
[3] "data/N03-20240101_42_GML/N03-20240101_42.geojson"
[4] "data/N03-20240101_43_GML/N03-20240101_43.geojson"
[5] "data/N03-20240101_44_GML/N03-20240101_44.geojson"

stringr::str_c関数は、複数の文字列を1つの文字列に結合する関数です。 seq関数は連続する数字のベクトルを作る関数です。 ここでは、40から44までの整数のベクトルを作成します。 つまり、seq(40, 44)は、c(40, 41, 42, 43, 44)と同じ結果を返します。

この2つの関数を組み合わせて、GeoJSONファイルのリストを作成します。

このリストに対してpurrr::mapを使うことで、リストにあるファイルを繰り返しsf::st_read関数で読み込みます。

purrr::map関数は、ベクトルの要素に対して、1つずつ関数を適用するための関数です。 ここでは、reg_filesの要素、つまり読み込むべきGeoJSONのファイル名を1つずつsf::st_read関数に渡して、それぞれのファイルを読み込みます。

region <- map(reg_files, st_read) |> bind_rows()
Reading layer `N03-20240101_40' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/N03-20240101_40_GML/N03-20240101_40.geojson' 
  using driver `GeoJSON'
Simple feature collection with 2380 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 129.9814 ymin: 33.00002 xmax: 131.1906 ymax: 34.25024
Geodetic CRS:  WGS 84
Reading layer `N03-20240101_41' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/N03-20240101_41_GML/N03-20240101_41.geojson' 
  using driver `GeoJSON'
Simple feature collection with 316 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 129.7368 ymin: 32.95048 xmax: 130.5421 ymax: 33.61899
Geodetic CRS:  WGS 84
Reading layer `N03-20240101_42' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/N03-20240101_42_GML/N03-20240101_42.geojson' 
  using driver `GeoJSON'
Simple feature collection with 9208 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 128.1045 ymin: 31.96711 xmax: 130.3903 ymax: 34.7286
Geodetic CRS:  WGS 84
Reading layer `N03-20240101_43' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/N03-20240101_43_GML/N03-20240101_43.geojson' 
  using driver `GeoJSON'
Simple feature collection with 2045 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 129.9388 ymin: 32.09492 xmax: 131.3295 ymax: 33.19518
Geodetic CRS:  WGS 84
Reading layer `N03-20240101_44' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/N03-20240101_44_GML/N03-20240101_44.geojson' 
  using driver `GeoJSON'
Simple feature collection with 2367 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 130.8247 ymin: 32.71444 xmax: 132.1774 ymax: 33.74028
Geodetic CRS:  WGS 84

読み込んだ結果、5つのsfオブフェクト（データフレームでもある）のリストが出来上がるのですが、それに対して、dplyr::bind_rows関数を適用することで、1つの大きなsfオブジェクト（データフレーム）に結合しています。

次に、読み込んだ地図を、自治体ごとにまとめます。 この行政区域データは、1つの自治体が複数のポリゴンから成り立っています。 これを1つのデータに集約するために、dplyr::group_by関数とdplyr::summarise関数を使った空間データに対する集計を実行します（少し時間がかかります）。 また、北九州空港の近くに「羽島」という小さな島があるのですが、これが所属未定地なので、このデータを取り除いています。

region <- region |> filter(N03_004 != "所属未定地") |>
  st_make_valid() |> 
  group_by(N03_001, N03_004, N03_007) |> 
  summarise(.groups = "drop") |> ungroup()

北部九州5県には、180の自治体があることがわかります。 国土数値情報の行政区域データは、塗り分け地図を作成するという用途には細かすぎるので、 sf::st_simplify関数を使って単純化しておきます（データ量も減らせます）。 そして、平面直角座標系II系に投影変換しておきます。

region <- region |> 
  st_simplify(dTolerance = 100) |> st_transform(6670) |> 
  select(origin_code = N03_007, geometry)

次に、自治体間の距離を図るために、市区町村役場のポイントデータを読み込みます。 行政区域データでやったのと同じように、GeoJSONファイル名のリストを作成して、 purrr::map関数で順番に読み込みます。

loc_files <- 
  str_c("data/P05-22_", seq(40, 44), "_GML/P05-22_", seq(40, 44), ".geojson")
location <- map(loc_files, st_read) |> bind_rows()
Reading layer `P05-22_40' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P05-22_40_GML/P05-22_40.geojson' 
  using driver `GeoJSON'
Simple feature collection with 3021 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 130.0375 ymin: 33.0057 xmax: 131.1882 ymax: 33.99226
Geodetic CRS:  JGD2011
Reading layer `P05-22_41' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P05-22_41_GML/P05-22_41.geojson' 
  using driver `GeoJSON'
Simple feature collection with 636 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 129.7736 ymin: 32.97361 xmax: 130.5416 ymax: 33.54343
Geodetic CRS:  JGD2011
Reading layer `P05-22_42' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P05-22_42_GML/P05-22_42.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1688 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 128.6016 ymin: 32.56244 xmax: 130.3792 ymax: 34.69298
Geodetic CRS:  JGD2011
Reading layer `P05-22_43' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P05-22_43_GML/P05-22_43.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1660 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 129.9844 ymin: 32.14 xmax: 131.2616 ymax: 33.181
Geodetic CRS:  JGD2011
Reading layer `P05-22_44' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P05-22_44_GML/P05-22_44.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1526 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 130.8488 ymin: 32.72246 xmax: 132.075 ymax: 33.73561
Geodetic CRS:  JGD2011

読み込んだ市区町村役場のデータから、まず、支所・出張所・連絡所のデータを削除します。 自治体コードと自治体名のデータをそれぞれorigin_codeとnameにdplyr::renameし、 stringr::str_remove関数を使って、自治体名から「役所」「役場」の文字列を削除します。

location <- location |> filter(P05_002 == 1) |> 
  select(origin_code = P05_001, name = P05_003, geometry) |> 
  mutate(name = str_remove(name, "役所|役場"))

政令指定都市については、区（区役所）のデータを使うことにして、指令指定都市の市役所のデータは取り除いておきます。 また、福岡県那珂川町は2018年に市政施行し那珂川市）に変わっていますが、那珂川町役場のデータが残っていますので、これも併せて削除しておきます。

location <- location |> 
  filter(name != "福岡市", name != "北九州市", 
         name != "熊本市", name != "那珂川町") |>
  st_transform(6670)

読み込んだ境界データと役所データを地図にしてみましょう。

ggplot() + 
  geom_sf(data = region) + 
  geom_sf(data = location) +
  theme_minimal()

他市区町村就業者数と通勤者数

e-statから、国勢調査のデータを入手します。

e-stat API利用のためのアプリケーションIDの取得

以下では、e-Stat APIを利用したデータ取得の演習を行います。 そのために、以下のページの案内にしたがって、各自でe-Statのユーザー登録を行い、アプリケーションIDを入手してください。

e-Statのユーザー登録とアプリケーションIDの入手

アプリケーションIDの管理

アプリケーションIDは、e-Statのマイページに行けば確認することができますが、いちいち確認するためにe-Statにログインするのも手間です。 アプリケーションIDは暗記するには長すぎますし、かといって、RスクリプトファイルにアプリケーションIDを平文で書き込んでしまうのも不安です。 このような場合には、OSが用意しているパスワード管理システム1を利用すると便利です。 ここでは、keyringパッケージを利用したアプリケーションIDの管理方法を紹介します。

アプリケーションIDの登録

RStudioのコンソール（Console）ペインに、以下のコマンドを入力します。

keyring::key_set('e-stat')

すると、パスワードを入力するウィンドウが開きますので、e-StatのページにあるアプリケーションIDをコピーして、入力欄にペーストしましょう。 そして［OK］ボタンをクリックすれば、登録完了です。 この情報は、自分でパスワード管理システムを操作して情報を削除しない限り、PC内部に暗号化された形で永続的に保持されます。

さて、jpstatパッケージを使って、e-Statのデータを取得する準備をします。 jpstatパッケージでは、アプリケーションIDを環境変数に登録する必要があります。 ここでは、アプリケーションIDを直接入力するのではなく、keyring::key_get関数を使って、パスワード管理システムに登録したアプリケーションIDを呼び出しています（このとき、もしかしたら、パスワードの入力を求められるかもしれませんが、その場合は適切なパスワードを入力してください）。

Sys.setenv(ESTAT_API_KEY = keyring::key_get("e-stat"))

これで、e-Stat APIを使う準備が整いましたので、早速、データを取得してみましょう。

他市区町村就業者数

jpstatパッケージを利用してe-Statからデータを取得するには、取得したいデータの統計表表示IDを把握する必要があります。 今回の他市区町村就業者数のデータは、男女，年齢（5歳階級），常住地又は従業地・通学地別就業者数－全国，都道府県，市区町村という統計表から取得します。 この統計表の統計表表示IDは”0003454500”ですので、jpstat::estat関数のstatsDataId引数にこれを与えます。 これで、統計表のメタデータを取得することができます。

worker <- estat(statsDataId = '0003454500')
worker
# ☐ tab:   表章事項                 [1] <code, name, level, unit>
# ☐ cat01: 男女                     [3] <code, name, level>
# ☐ cat02: 年齢                     [25] <code, name, level>
# ☐ cat03: 常住地又は従業地・通学地 [19] <code, name, level>
# ☐ area:  全国，都道府県，市区町村 [1,965] <code, name, level, parentCode>
# ☐ time:  時間軸（年次）           [1] <code, name, level>
# 
# Please `activate()`.

e-Stat APIを使ったデータ抽出では、表彰事項、カテゴリー（年齢、性別など、統計表によって異なる）、地域、時間軸ごとにデータが格納されているので、データが必要なカテゴリー項目を設定する必要があります。 そのために、設定項目ごとに、navigatr::activate（必要であればnavigatr::rekeyも）とdplyr::filterを繰り返し適用して、最後にdplyr::collectでデータを取得します。

worker <- worker |> 
  activate(tab) |> 
  select() |> 
  activate(cat01) |> 
  rekey('gender') |> 
  select(name) |> 
  activate(cat02) |> 
  filter(name == '総数') |> 
  select() |> 
  activate(cat03) |> 
  filter(name == '他市区町村で従業・通学') |> 
  select() |> 
  activate(area) |> 
  rekey('origin') |> 
  select(code) |> 
  activate(time) |> 
  select() |> 
  collect(n = "worker") |> 
  mutate(worker = parse_number(worker, na = "-"))
The total number of data is 5895.
e-Statデータ取得の処理内容（就業者データ）
rekey filter select
tab
cat01 gender name
cat02 “総数”のみ
cat03 “他市区町村で従業・通学”のみ
area origin code
time

ここでの処理内容を表に整理しています。

  • 表彰事項（tab）には処理はしていません。select関数の中身も空白なので、最終的なデータからも落としています。
  • カテゴリー1（cat01）の項目名をgenderに設定しています。データとして、name（日本語文字列：“総数”、“男”、“女”）を取得します。
  • カテゴリー2（cat02）の年齢については、“総数”データだけを取得します。
  • カテゴリー3（cat03）は、“他市区町村で従業・通学”しているデータだけを取り出します。
  • 地域（area）は、名前をoriginにしたうえで、全てのデータのコード（自治体コード）を取得します。
  • 時間軸（time）も、tabと同様にデータから落とします。

最後に、readr::parse_number関数を使って就業者数を文字列から数値にしています。

dplyr::glimpse関数で、取得したデータの内容を確認してみましょう。

glimpse(worker)
Rows: 5,895
Columns: 3
$ gender_name <chr> "総数", "総数", "総数", "総数", "総数", "総数", "総数", "…
$ origin_code <chr> "00000", "01000", "01100", "01101", "01102", "01103", "011…
$ worker      <dbl> 25015142, 694338, 417629, 28085, 61207, 53219, 47849, 5929…

通勤者数

他市区町村就業者数のデータは、男女，就業・通学，常住地（全国，都道府県，市区町村）別就業者・通学者数－全国［総数］，都道府県，市区町村（従業地・通学地）という統計表から取得します。 この統計表の統計表表示IDは”0003454525”ですので、jpstat::estat関数のstatsDataId引数にこれを与えます。 これで、統計表のメタデータを取得することができます。

commuter <- estat(statsDataId = '0003454525')
commuter
# ☐ tab:   表章事項                                           [1] <code, name, level, unit>
# ☐ cat01: 男女                                               [3] <code, name, level>
# ☐ cat02: 常住地                                             [11] <code, name, level>
# ☐ cat03: 全国，都道府県，市区町村（常住地）                 [1,965] <code, name, level, parentCode>
# ☐ cat04: 就業・通学                                         [5] <code, name, level>
# ☐ area:  全国［総数］，都道府県，市区町村（従業地・通学地） [1,965] <code, name, level, parentCode>
# ☐ time:  時間軸（年次）                                     [1] <code, name, level>
# 
# Please `activate()`.
commuter <- commuter |> 
  activate(tab) |> 
  select() |> 
  activate(cat01) |> 
  rekey('gender') |> 
  select(name) |> 
  activate(cat02) |>
  filter(name == '他市区町村に常住') |>
  select() |>
  activate(cat03) |>
  rekey('origin') |> 
  select(code) |>
  activate(cat04) |>
  filter(name == '15歳以上就業者') |>
  select() |>
  activate(area) |>
  rekey('destination') |> 
  filter(name %in% c('佐賀市', '久留米市')) |> 
  select(name) |>
  activate(time) |>
  select() |>
  collect(n = "commuter") |> 
  mutate(commuter = parse_number(commuter, na = "-"))
The total number of data is 11790.
e-Statデータ取得の処理内容（通勤者データ）
rekey filter select
tab
cat01 gender name
cat02 “他市区町村に常住”のみ
cat03 origin code
cat04 “15歳以上就業者”のみ
area destination “佐賀市”または”久留米市” name
time

ここでの処理内容を表に整理しています。

glimpse(commuter)
Rows: 11,790
Columns: 4
$ gender_name      <chr> "総数", "総数", "総数", "総数", "総数", "総数", "総数…
$ origin_code      <chr> "00000", "00000", "01000", "01000", "01100", "01100",…
$ destination_name <chr> "久留米市", "佐賀市", "久留米市", "佐賀市", "久留米市…
$ commuter         <dbl> 36433, 31163, NA, 1, NA, NA, NA, NA, NA, NA, NA, NA, …

データの可視化

性別の他市区町村就業者数の塗り分け地図を作ってみましょう。

region_w <- region |> 
  left_join(worker, by = join_by(origin_code))
glimpse(region_w)
Rows: 540
Columns: 4
$ origin_code <chr> "41346", "41346", "41346", "41345", "41345", "41345", "412…
$ geometry    <MULTIPOLYGON [m]> MULTIPOLYGON (((-55156.34 3..., MULTIPOLYGON …
$ gender_name <chr> "総数", "男", "女", "総数", "男", "女", "総数", "男", "女"…
$ worker      <dbl> 7050, 3915, 3135, 3051, 1621, 1430, 5685, 3276, 2409, 2065…

行数が540になっているのは、市町村数（180）の3倍です。 性別の種類が3（“総数”、“男”、“女”）なので、left_joinしたことで、自動的に540行のデータになっています。

ggplotで地図にしてみましょう。 ggplot2::scale_fill_gradient関数のtrans引数を"log10"にすることで、 塗り分けの軸を常用対数にしています。 また、ggplot2::facet_wrap関数で、性別の種類（“総数”、“男”、“女”）ごとの地図を並べて表示しています。

ggplot() + geom_sf(aes(fill = worker), data = region_w) +
  scale_fill_gradient(low = 'white', high = 'blue', trans = 'log10') +
  facet_wrap(vars(gender_name)) +
  theme_bw()

次に、佐賀市・久留米市への通勤者数による塗り分け地図を作ってみましょう。

region_c <- region |> 
  left_join(commuter, by = join_by(origin_code))
glimpse(region_c)
Rows: 1,080
Columns: 5
$ origin_code      <chr> "41346", "41346", "41346", "41346", "41346", "41346",…
$ geometry         <MULTIPOLYGON [m]> MULTIPOLYGON (((-55156.34 3..., MULTIPOL…
$ gender_name      <chr> "総数", "総数", "男", "男", "女", "女", "総数", "総数…
$ destination_name <chr> "久留米市", "佐賀市", "久留米市", "佐賀市", "久留米市…
$ commuter         <dbl> 1537, 874, 748, 499, 789, 375, 290, 455, 143, 236, 14…

left_joinの結果、今度はデータの列数が1,080になっていますね。 これは、性別の3区分に加えて、目的地の2区分（“佐賀市”、“久留米市”）が加わったために、3×2で6倍になるためです。

ggplot() + geom_sf(aes(fill = commuter), data = region_c) +
  scale_fill_gradient(low = 'white', high = 'red', trans = 'log10') +
  facet_grid(cols = vars(gender_name), rows = vars(destination_name)) +
  theme_bw()

佐賀市・久留米市までの距離

各市区町村から佐賀市・久留米市までの通勤者数と距離との関係を分析するために、 各市区町村から佐賀市・久留米市までの距離のデータが必要です。 ここでは、sf::st_distanceを使って計算しています。

sf::st_drop_geometry関数を使って、locationデータのgeometryを落として、自治体コードと自治体名だけのデータフレームを作っています。 そのデータフレームに対して、dplyr::mutate関数を使って、 各自治体から佐賀市までの距離のデータ列（佐賀市）と、 各自治体から久留米市までの距離のデータ列（久留米市）を追加しています。

saga   <- location |> filter(name == "佐賀市")
kurume <- location |> filter(name == "久留米市")
distance <- location |> st_drop_geometry() |> 
  mutate(佐賀市 = st_distance(location, saga) |> as.vector(),
         久留米市 = st_distance(location, kurume) |> as.vector())
glimpse(distance)
Rows: 180
Columns: 4
$ origin_code <chr> "40101", "40103", "40105", "40106", "40107", "40108", "401…
$ name        <chr> "門司区", "若松区", "戸畑区", "小倉北区", "小倉南区", "八…
$ 佐賀市      <dbl> 97041.164, 85711.545, 85576.524, 86843.958, 84538.224, 819…
$ 久留米市    <dbl> 80686.67, 70828.86, 70314.36, 70902.12, 68120.66, 66640.05…

これを「tidy」なデータ形式にしておきます。 tidyr::pivot_longer関数を使って、佐賀市と久留米市の2つのデータ列を、 cityというデータ列（佐賀市または久留米市のどちらかの文字列がデータとして入る）と distanceというデータ列（佐賀市または久留米市までの距離がデータとして入る）に再編しています。 その後、距離を1,000で割ることで、距離の単位をmからkmに変換しています。

distance <- distance |> 
  pivot_longer(佐賀市:久留米市, names_to = 'destination_name', values_to = 'distance') |> 
  mutate(distance = distance / 1000)
glimpse(distance)
Rows: 360
Columns: 4
$ origin_code      <chr> "40101", "40101", "40103", "40103", "40105", "40105",…
$ name             <chr> "門司区", "門司区", "若松区", "若松区", "戸畑区", "戸…
$ destination_name <chr> "佐賀市", "久留米市", "佐賀市", "久留米市", "佐賀市",…
$ distance         <dbl> 97.04116, 80.68667, 85.71155, 70.82886, 85.57652, 70.…
tidyr::pivot_longer

tidyr::pivot_longerは、「tidy」なデータを作るためには不可欠で便利な関数ですが、 最初はどのような働きをするのかイメージしにくいと思います。 tidyrのCheatsheet にある図がわかりやすいかもしれません。

逆にlong型のデータをwide型に変換する関数はtidyr::pivot_widerです。

距離と通勤率との関係

距離と通勤率との関係を見てみます。

dplyr::left_join関数を使って、行政区域データに、就業者数、通勤者数、距離のデータを結合します。 さらに、dplyr::mutate関数を使って、通勤者比率（＝通勤者÷就業者）のデータ列を追加します。

region_d <- region |> 
  left_join(worker, by = join_by(origin_code)) |> 
  left_join(commuter, by = join_by(origin_code, gender_name)) |> 
  left_join(distance, by = join_by(origin_code, destination_name)) |> 
  mutate(rate = commuter / worker)
glimpse(region_d)
Rows: 1,080
Columns: 9
$ origin_code      <chr> "41346", "41346", "41346", "41346", "41346", "41346",…
$ geometry         <MULTIPOLYGON [m]> MULTIPOLYGON (((-55156.34 3..., MULTIPOL…
$ gender_name      <chr> "総数", "総数", "男", "男", "女", "女", "総数", "総数…
$ worker           <dbl> 7050, 7050, 3915, 3915, 3135, 3135, 3051, 3051, 1621,…
$ destination_name <chr> "久留米市", "佐賀市", "久留米市", "佐賀市", "久留米市…
$ commuter         <dbl> 1537, 874, 748, 499, 789, 375, 290, 455, 143, 236, 14…
$ name             <chr> "みやき町", "みやき町", "みやき町", "みやき町", "みや…
$ distance         <dbl> 5.044407, 15.849471, 5.044407, 15.849471, 5.044407, 1…
$ rate             <dbl> 0.218014184, 0.123971631, 0.191060026, 0.127458493, 0…

ggplot2::ggplotを使って、距離と通勤者比率の散布図を描きましょう。 ggplot2::facet_grid関数を使って、都市（佐賀市と久留米市）と性別ごとの散布図にします。

ggplot(region_d, aes(x = distance, y = rate)) + 
  geom_point() +
  facet_grid(rows = vars(destination_name), cols = vars(gender_name)) +
  theme_bw()
Warning: Removed 335 rows containing missing values or values outside the scale range
(`geom_point()`).

距離が遠くなるほど、通勤者比率が減少しているように見えますが、両者にどのような関係があるのかよくわかりません。 そこで、x軸とy軸を対数軸にしてみます。

ggplot(region_d, aes(x = distance, y = rate)) + 
  geom_point() +
  scale_y_log10() + scale_x_log10() + 
  facet_grid(rows = vars(destination_name), cols = vars(gender_name)) +
  theme_bw()
Warning in scale_x_log10(): log-10 transformation introduced infinite values.
Warning: Removed 335 rows containing missing values or values outside the scale range
(`geom_point()`).

両対数をとると、散布図は右下がりの線分に沿っているように見えます。 そこで、ggplot2::geom_smooth関数を使って、回帰直線を表示してみましょう。

ggplot(region_d, aes(x = distance, y = rate)) + 
  geom_point() +
  scale_y_log10() + scale_x_log10() + 
  geom_smooth(method = lm, se = FALSE) +
  facet_grid(rows = vars(destination_name), cols = vars(gender_name)) +
  theme_bw() 
Warning in scale_x_log10(): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 335 rows containing non-finite outside the scale range
(`stat_smooth()`).
Warning: Removed 335 rows containing missing values or values outside the scale range
(`geom_point()`).

回帰係数を確認してみましょう。 都市と性別でデータをグループ化し、tidyr::nest関数を使って、グループ化したデータフレームをリスト型の列の要素にしてしまいます。 そのネストしたデータフレームに対して、順番に回帰分析（lm関数）を実行し、その結果をlm列に格納しています（purrr::mapを使って処理しています）。 さらに回帰分析の結果を、broom::tidy関数でtidyな形式に変換したものをcoef列に追加し、 最後にtidyr::unnestでcoefをアンネスト（ネストを解除）しています。

region_d |> 
  group_by(destination_name, gender_name) |> 
  nest() |> 
  mutate(lm = map(data, function(x){lm(log(rate) ~ log(distance), data = x)})) |> 
  mutate(coef = map(lm, broom::tidy)) |> 
  unnest(coef)
# A tibble: 12 × 9
# Groups:   destination_name, gender_name [6]
   gender_name destination_name data   lm     term  estimate std.error statistic
   <chr>       <chr>            <list> <list> <chr>    <dbl>     <dbl>     <dbl>
 1 総数        久留米市         <sf>   <lm>   (Int…     3.80     0.492      7.72
 2 総数        久留米市         <sf>   <lm>   log(…    -2.52     0.127    -19.8 
 3 総数        佐賀市           <sf>   <lm>   (Int…     5.17     0.783      6.61
 4 総数        佐賀市           <sf>   <lm>   log(…    -2.92     0.200    -14.6 
 5 男          久留米市         <sf>   <lm>   (Int…     3.59     0.455      7.89
 6 男          久留米市         <sf>   <lm>   log(…    -2.39     0.118    -20.3 
 7 男          佐賀市           <sf>   <lm>   (Int…     4.70     0.770      6.10
 8 男          佐賀市           <sf>   <lm>   log(…    -2.71     0.198    -13.7 
 9 女          久留米市         <sf>   <lm>   (Int…     4.40     0.630      6.98
10 女          久留米市         <sf>   <lm>   log(…    -2.78     0.170    -16.3 
11 女          佐賀市           <sf>   <lm>   (Int…     6.61     0.803      8.23
12 女          佐賀市           <sf>   <lm>   log(…    -3.42     0.215    -15.9 
# ℹ 1 more variable: p.value <dbl>
tidyr::nest

tidyr::nestとtidyr::unnestについても、Cheatsheet の図を貼り付けておきます。 参考にしてください。

n

線形回帰の結果を表にまとめると、以下のようになります。 \(\beta\)については、回帰の係数（\(b\)）ではなく\(\exp(b)\)の値であることに注意してください。

都市 性別 \(\beta\) \(\alpha\)
久留米市 女 81.34652 2.780654
久留米市 男 36.05946 2.390588
久留米市 総数 44.55006 2.520781
佐賀市 女 744.36591 3.416006
佐賀市 男 109.67832 2.712030
佐賀市 総数 176.53639 2.918836

\(\beta\)はその都市が通勤者を惹きつける力（魅力）のようなものを表しています。 \(\alpha\)は距離の影響に関するパラメータで、値が大きいほど距離が効く（距離が遠くなると通勤者比率が大きく減る）ことを意味しています。

この結果表を見ると、佐賀市においても、久留米市においても、女性の方が男性よりも\(\alpha\)も\(\beta\)も大きいことがわかります。 女性は男性に比べ遠くに通勤しないが、近場であれば他市区町村への通勤を厭わない、と解釈できるかもしれません。

性別による回帰直線を重ねて表示すると、次の図のようになります。

ggplot(region_d, aes(x = distance, y = rate, color = gender_name)) + 
  geom_point() +
  scale_y_log10() + scale_x_log10() + 
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(vars(destination_name)) +
  theme_bw()
Warning in scale_x_log10(): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 335 rows containing non-finite outside the scale range
(`stat_smooth()`).
Warning: Removed 335 rows containing missing values or values outside the scale range
(`geom_point()`).

また、いずれの性別においても、久留米市に比べて佐賀市の\(\alpha\)及び\(\beta\)が大きくなっています。 佐賀市の方が久留米市よりも市外からの通勤者を惹きつけるが、遠くからの通勤者比率は少なくなる、という結果になりました。

通勤先都市による回帰直線を重ねて表示すると、次の図のようになります。

ggplot(region_d, aes(x = distance, y = rate, color = destination_name)) + 
  geom_point() +
  scale_y_log10() + scale_x_log10() + 
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(vars(gender_name)) +
  theme_bw() 
Warning in scale_x_log10(): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 335 rows containing non-finite outside the scale range
(`stat_smooth()`).
Warning: Removed 335 rows containing missing values or values outside the scale range
(`geom_point()`).

おわりに

今回は、重力モデルについて実際にデータを操作しながら、その雰囲気を味わってもらいました。 重力モデルといいつつ、実際にはそれを簡略化した形（距離・移動モデル）ですので、距離と移動との関係を体感していただけたのではないかと思います。

脚注

  1. Windowsでは資格情報マネージャー、MacOSではキーチェーンアクセスが該当します。↩︎