経営のための統計学II
第9回 空間データ分析1
高橋 慎
2021/11/19
空間データ
空間データ (spatial data) は、地図 (map) を使ってビジュアル化することで最もうまく分析できる。
ここでは、以下の2種類の空間データを扱う。
空間ポイント・データ (spatial point data)：地図上に点の集合としてプロットできる。
空間ポリゴン・データ (spatial polygon data)：群、区、県などの、特定の地域の境界に対応する地図上の一連の連結された点の連続体として表される。
また、ある期間記録された空間ポイントやポリゴン・データの集合で、時間の経過による空間パターンの変化を表す時空間データ (spatial-temporal data) についても検討する。
ロンドンにおける1854年のコレラ発生
イギリスの内科医ジョン・スノーは、著書『コレラの伝染形態』(On the Mode of Communication of Cholera) において、コレラ感染による死亡事例の空間的な分布をビジュアル化し、地図の効果的な使い方を示した。
スノーは、1854年のコレラ発生時にロンドンのソーホー近郊における感染による死亡事例についての空間ポイント・データを収集し、その情報を地図上にプロットした。
ジョン・スノーによるロンドンにおけるコレラ感染による死亡事例の地図を再現したもの
コレラ感染による死亡事例を示す黒い長方形の領域は、地図の中央に位置しているブロード通り (Broad Street) の揚水ポンプの周りに集中している。
- 地図上にはすべての揚水ポンプが「PUMP」というラベル付きの黒丸で示されている。
この地図からスノーは、コレラ感染による死亡事例はブロード通り周辺に集中していることを発見した。
彼は、コレラが下水で汚染された水によって広がったのではないかと推測した。
水道水の詳しい検査や地域住民へのインタビューを含む大規模な調査の後、スノーはブロード通りとケンブリッジ通りとの交差点にある揚水ポンプがコレラ発生源であると結論づけた。
また、スノーはサザーク&ボクソール社 (Southwark and Vauxhall Company) の水道水供給にロンドンのコレラ拡大の責任があることを示すために「壮大な自然実験 (grand natural experiment)」を採用した。
自然実験 (natural experiment)とは、研究者の介入はないが実験と似た現実世界の状況を意味する。
以下の図は、スノーが自然実験を行った地域をビジュアル化するために用いた空間ポリゴン地図を再現したものである。
ジョン・スノーによる自然実験の地図を再現したもの
ランベス社 (Lambeth Company) は南部地域（赤色で示されている）へ、サザーク&ボクソール社はテムズ川沿岸地域（青色で示されている）へ水道水を供給していた。
スノーは、ロンドンのテムズ川周辺の地域で、一部の人々はランベス社から（清潔な）水道水供給を受け、別の人々はサザーク&ボクソール社から（汚染された）水道水供給を受けていたことに注目した。
十分な調査の後、スノーはいかなる交絡因子 (confounder) もこの自然実験には影響を与えなかったと結論づけた。
- この文脈における交絡因子とは、水道会社とある地域でのコレラ発生率の両方と関連をもつ変数を意味する。
つまり、コレラ感染率における違いは水道会社の選択に起因する。
コレラにより死亡した人の住所と、彼らに水道水を供給していた会社とを照合することで、スノーは大多数の死亡事例はサザーク&ボクソール社によって水道水を供給されていた家庭で起きたことを明らかにした。
- スノーの著書は「地図を用いて空間データをビジュアル化することで研究者は過去に知ることのできなかったパターンを発見し、それらを納得のいく形で説明できる」という空間データ分析の威力を示している。
Rにおける空間データ
- 以下の図は、2008年のアメリカ大統領選挙の結果を効果的にビジュアル化している。
2008年アメリカ大統領選挙での選挙人団の投票先マップ
これは、空間ポリゴン・データの例で、各州がいくつかの点をつないで境界を作るポリゴン（多角形）を描く。
こうしておくと、バラク・オバマ（ジョン・マケイン）がある州で相対多数の票を得れば、そのポリゴンすなわち州を青（赤）に色づけすることができる。
Rでは、mapsパッケージが様々な地図作成ツールと多くの空間データベースを提供しており、世界の様々な都市の空間データベースが含まれている。例えば、
us.citiesと呼ばれるアメリカ合衆国の都市のデータフレームなどがある。パッケージに収録されているどのデータフレームも、
data()関数を使って読み込むことができる。
## name country.etc pop lat long capital
## 1 Abilene TX TX 113888 32.45 -99.74 0
## 2 Akron OH OH 206634 41.08 -81.52 0
## 3 Alameda CA CA 70069 37.77 -122.26 0
## 4 Albany GA GA 75510 31.58 -84.18 0
## 5 Albany NY NY 93576 42.67 -73.80 2
## 6 Albany OR OR 45535 44.62 -123.09 0このデータセットには、以下の変数が含まれている。
name：名前country.etc：州pop：人口lat：緯度long：経度capital：首都であればcapital = 1、州都であればcapital = 2、どちらでもなければcapital = 0
map()関数を用いれば、空間データベースにアクセスし、その中のデータをビジュアル化することができる。例えば、アメリカ合衆国を図示するには、引数
databaseを"usa"に設定する。空間上の点は、経度と緯度の情報を
x座標とy座標の入力値としてpoints()関数を使って簡単に地図上へ加えることができる。各州都は実線の円で描かれ、その大きさを
cex = capitals$pop / 500000のように州都の人口に比例させることができる。地図が描かれた後で
title関数を用いてタイトルを加えることもできる。
map(database = "usa")
capitals <- subset(us.cities, capital == 2) # 州都を部分集合化
## 経度と緯度を用いて人口に比例した点を追加する
points(x = capitals$long, y = capitals$lat, col = "blue", cex = capitals$pop / 500000, pch = 19)
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
title("アメリカ合衆国の州都") # タイトルを追加する
別の例として、カリフォルニア州をプロットしてみる。
アメリカ合衆国の空間ポリゴン・データを含む
stateデータベースを用い、引数regionsを"California"に指定する。
このカリフォルニア州の地図に、州内の人口の多い7大都市を加える。
これらの都市をデータから抽出するには、順番が付された指標ベクトルを返す
order関数を用いる。order()関数では、引数decreasingをTRUEやFALSEに指定することで降順や昇順で並べ替え（ソート）することができる。sort()関数は、順序づけられたベクトルそのものを返す。
## name country.etc pop lat long capital
## 3 Alameda CA CA 70069 37.77 -122.26 0
## 10 Alhambra CA CA 88857 34.08 -118.13 0
## 11 Aliso Viejo CA CA 41975 33.57 -117.73 0
## 15 Altadena CA CA 43280 34.19 -118.13 0
## 20 Anaheim CA CA 334909 33.84 -117.87 0
## 27 Antioch CA CA 109485 37.99 -121.80 0sind <- order(cal.cities$pop, decreasing = TRUE) # 人口を降順に並べる
top7 <- sind[1:7] # 人口の多い上位7都市
cal.cities[top7, ] # 人口の多い上位7都市のデータを表示## name country.etc pop lat long capital
## 521 Los Angeles CA CA 3911500 34.11 -118.41 0
## 801 San Diego CA CA 1299352 32.81 -117.14 0
## 804 San Jose CA CA 897883 37.30 -121.85 0
## 802 San Francisco CA CA 723724 37.77 -122.45 0
## 517 Long Beach CA CA 486571 33.79 -118.16 0
## 778 Sacramento CA CA 480392 38.57 -121.47 2
## 337 Fresno CA CA 472517 36.78 -119.79 0## [1] 3911500 1299352 897883 723724 486571 480392 472517points()関数を用いてこれらの都市を地図に加え、それらの名前もtext()関数を用いて追加する。cex = 0.6として文字サイズを小さく（0.6倍に）している。
map(database = "state", regions = "California")
points(x = cal.cities$long[top7], y = cal.cities$lat[top7], pch = 19, col = "blue")
## 円と重なることを避けるため、経度に定数を加える
text(x = cal.cities$long[top7] + 2, y = cal.cities$lat[top7], label = cal.cities$name[top7],
cex = 0.6, col = "blue")
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
title("カリフォルニア州の大都市")
ここで、
Rで空間ポリゴン・データがどのようになっているのか検討する。プロット作成が行われないようにするため、
map関数の引数plotをFALSEに設定すると、x（\(x\)座標すなわち経度）とy（\(y\)座標すなわち緯度）として保存された座標の数列を伴うリスト・オブジェクトを返す。このリストの中では名前が
namesとして保存されている異なるポリゴンはNAによって分割されている。
## [1] "x" "y" "range" "names"ベクトル
xの長さを計算して、アメリカ合衆国の地図を作成するために用いられた座標の数を確認することができる。cbind()関数を使って\(x\)座標と\(y\)座標を行列にまとめた後、最初の数件の観察値を表示する。
## [1] 7252## [,1] [,2]
## [1,] -101.4078 29.74224
## [2,] -101.3906 29.74224
## [3,] -101.3620 29.65056
## [4,] -101.3505 29.63911
## [5,] -101.3219 29.63338
## [6,] -101.3047 29.64484アメリカ合衆国の地図は、7252件の座標ペアで構成されていることがわかる。
map関数はこれらの点を連結し、地図を作成する。
空間データは、空間に広がるパターンについての情報を含み、地図を用いてビジュアル化できる。
空間ポイント・データでは事象の地点を地図上の点として表すのに対し、空間ポリゴン・データは地理的領域を地図上の点を連結することで表す。
日本の空間データ
- 日本の情報を確認したい場合は
mapdataパッケージを用いる。
緯度と経度の範囲を
xlimとylimで指定して拡大することもできる。オプション
interior = FALSEで都道府県の境界を非表示にして、boundary =以降で境界を波線で表示するように指定する。
map("japan", interior = FALSE, xlim = c(138,141), ylim = c(34,36))
map("japan", boundary = FALSE, lty = 2, add = TRUE)
map.axes()
法政大学（BT）の緯度と経度を調べて表示する。
map("japan", interior = FALSE, xlim = c(138.5,140.5), ylim = c(35,36))
map("japan", boundary = FALSE, lty = 2, add = TRUE)
points(139.740631, 35.695374, col = "blue", pch=20)
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
text(139.740631, 35.695374 + .05, "法政大学", col = "blue")
map.axes()
参考文献
今回の内容は下巻の第5章の内容に基づく。
原著 Quantitative Social Science (Princeton University Press) のウェブ・ページでデータなどをダウンロードできる。