はじめに
久保研介研究会のサブゼミの資料です！今回から本格的にメッシュ統計を扱っていきます。グループワーク多めなので、Rの操作を復習しつつ頑張っていきましょう〜
ちなみに、今日のゴールはいつも通り2つです。
- メッシュ統計データの活用の仕方がわかった！
- Rの操作をなんとなく思い出した！
やること
1. メッシュ統計データについて知る
2. 緯度経度をメッシュコードに変換する
3. メッシュ統計データを組み合わせて分析してみる
1. メッシュ統計データとは？
- メッシュ統計データとは、緯度・経度に基づき地域を隙間なく網の目（メッシュ）の区域に分けて，それぞれの区域に関する統計データです。
- ちなみに、メッシュコードとは網の目（メッシュ）の区域の番号になります。
- **これ**を読んでおけば、大体わかります。
- **より詳しくは、こちら**の資料を読んでください。
- メッシュ統計データは**e-stat**からダウンロード出来ます。
2. 緯度経度をメッシュコードに変換する
メッシュコードは緯度経度から算出することが出来ます。今回は2次メッシュ（10km四方）、3次メッシュ（1km四方）、4次メッシュ（500m四方）のメッシュコードを算出してみます。jpmeshというライブラリにあるcoords_to_meshというコマンドで簡単に出来ます。
# メッシュコード変換に使用するライブラリ
install.packages("jpmesh")
library(jpmesh)
# （）内では、経度、緯度の順番に座標を指定します。
gas_naha$meshcode3 <- coords_to_mesh(gas_naha$longitude, gas_naha$latitude)
# 同じ様にして、2次メッシュを算出してみる。mesh_sizeでメッシュの大きさを指定します。
gas_naha$meshcode2 <- coords_to_mesh(gas_naha$longitude, gas_naha$latitude, mesh_size = 10)
# つまり、500m四方の4次メッシュだと…
gas_naha$meshcode4 <- coords_to_mesh(gas_naha$longitude, gas_naha$latitude, mesh_size = 0.5)
##Groupwork①
3次メッシュコードごとにガソリンスタンドの数を集計してみてください。集計結果はgas_mesh3としてください。
※ヒント：gas_mesh3 <- gas_naha %>%...
解答
gas_mesh3 <- gas_naha %>%
group_by(meshcode3) %>%
summarise(number_store = n())
3. メッシュ統計データを組み合わせて分析してみる
3.0. テーマ設定
ここで、ガソリンスタンドの出店行動を規定する要因を分析してみようと思います。
まずは、**人口が多い地域ほど出店の確率は高まるのか?**という仮説を検証してみましょう。データの観察単位(データの行に相当)は3次メッシュ（1km四方）とします。
分析の全体像
- 那覇市の3次メッシュコードを集める
- その3次メッシュコードに対してガソリンスタンドが出店している地域なら1をとるダミー変数を作成する。（被説明変数）
- 3次メッシュごとの人口データを追加する（説明変数）
3.1.那覇市の3次メッシュコードを集める
市区町村別のメッシュコードは**このサイト**からゲット出来ます。各自で沖縄県のデータを取得してください。
##Groupwork②
取得したデータから那覇市のデータだけ抽出してください。また、メッシュコード以外の列（2つ）はいらないので、削除してください。抽出したデータはnaha_mesh3としましょう。
※ヒント：subsetコマンドでデータを抽出することが出来ます。
解答
## 回答
okinawa_mesh3 <- read.csv('47.csv', fileEncoding = "shift-jis")
# データの抽出にはsubsetというコマンドを使います。
naha_mesh3 <- subset(okinawa_mesh3, 市区町村名 == "那覇市", 基準メッシュコード)
# 列名をmeshcode3とします。
colnames(naha_mesh3) <-"meshcode3"
## 別解
naha_mesh3 <- okinawa_mesh3 %>%
filter(市区町村名 == "那覇市") %>%
select(基準メッシュコード) %>%
rename(meshcode3 = 基準メッシュコード)
これで**「3.1.那覇市の3次メッシュコードを集める」**は達成ですね。
3.2.ガソリンスタンドが出店している地域なら1をとるダミー変数を作成する。（被説明変数）
つまりは、那覇市の3次メッシュコードのデータ**(naha_mesh3)に**、メッシュコードごとのガソリンスタンドの数のデータ**(gas_mesh3)**を結合すればOKですね。
##Groupwork③
まず、naha_mesh3とgas_mesh3を結合して、那覇市の3次メッシュコードに対してガソリンスタンド出店数の情報を加えてください。結合したデータはgas_naha_3とします。次に、出店数のデータをもとに出店している地域なら1をとるダミー変数を作成してください。ダミー変数の名前はentryとしておきましょう。
※ヒント：mergeでデータを結合しましたね。結合する時に注意することはなんでしたっけ？？ダミー変数を作る時はifelseを使いましたね。
解答
# naha_mesh3のデータは残しておきたいことに注意
gas_naha_3 <- merge(naha_mesh3, gas_mesh3, by= "meshcode3", all.x = TRUE)
# 欠損値は出店していない地域を示すので、0を代入
gas_naha_3$number_store[is.na(gas_naha_3$number_store)] <- 0
# number_store > 0の時に1をとればよい
gas_naha_3$entry <- ifelse(gas_naha_3$number_store > 0, 1, 0)
これで**「3.2.ガソリンスタンドが出店している地域なら1をとるダミー変数を作成する。（被説明変数）」**は達成です。
3.3.人口データを追加する（説明変数）
ついに最後です！メッシュ統計データを使っていきましょう。メッシュ統計データは政府統計サイトe-statからとってこれます。下の画像にもあるようにデータの種類は豊富です。（総務省統計局：第１章 地域メッシュ統計の特質・沿革より）
今回は人口のデータを使用するので、**こちら**からダウンロードしてください。那覇市の3次メッシュコードは3927...となっているので、M3972というデータをダウンロードしましょう。
早速データをR上に読み込んでみます。データの形が汚いので、ちょっとめんどくさいですが…
# sepの指定に注意する
pop <- read.csv("tblT000846S3927.txt", sep = ",", fileEncoding = "shift-jis")
# 必要なデータだけに絞る。必要な行は1行目以外。列は1列目と5列目だけ。
pop <- pop[-1,c(1, 5)]
# 変数名を変える
colnames(pop) <- c("meshcode3", "pop")
##Groupwork④
まず、gas_naha_3とpopを結合して、地域人口のデータを加えてみましょう。結合したデータ名はgas_naha_3とします。地域人口を示す変数名はpopとします。また、popに欠損値を含むデータは人が住めない地域として削除してください。
※ヒント：mergeでデータを結合しましたね。結合する時に注意することはなんでしたっけ？？
解答
推定
プロビットモデルで分析してみました。有意な影響は観察出来ませんでしたね、、もっと他の変数を入れて分析する必要がありそうです。
# 推定
out <- glm(entry ~ pop, family = binomial(link = "probit"), data = gas_naha_3)
install.packages("margins")
library(margins)
summary(margins(out))
さいごに
お疲れ様でした！公的統計を組み合わせて分析してみると、分析してる感が増しますよね〜。