28 GIS の基本
28.1 概要
データに空間情報があることで、アウトブレイクの状況について多くの洞察が得られ、次のような質問に答えることができます。
- 現在の疾病のホットスポットはどこか？
- ホットスポットは時系列でどのように変化しているか？
- 医療施設へのアクセスはどうなっているか？改善が必要か？
この GIS の章では、疫学業務担当者がアウトブレイク対応において知りたいニーズに応えることを目的としています。ここでは、tmap と ggplot2 パッケージを使って、基本的な空間データの可視化方法を学びます。また、sf パッケージを使用して、基本的な空間データの管理とクエリの方法についても説明します。最後に、spdep パッケージを用いて、空間的関係、空間的自己相関、空間回帰など、空間統計の概念に簡単に触れます。
28.2 重要な用語
ここでは、いくつかの重要な用語を紹介します。GIS と空間分析について詳しく知りたい方は、「リソース」節に掲載されている長いチュートリアルやコースをご覧になることをお勧めします。
GIS（Geographic Information System, 地理情報システム） - GIS とは、空間データを収集、管理、分析、視覚化するためのフレームワークまたは環境のことです。
GIS ソフトウェア
一般的な GIS ソフトウェアの中には、マウス操作で地図の作成や空間分析ができるものがあります。これらのツールには、コードを覚える必要がない、アイコンや機能を手動で選択して地図上に配置するのが簡単であるなどの利点があります。ここでは、人気のある2つのソフトウェアを紹介します。
ArcGIS - ESRI 社が開発した商用の GIS ソフトウェアで、非常に人気がありますますが、かなり高価です。
QGIS - フリーのオープンソース GIS ソフトウェアで、ArcGIS でできることもおおむねできます。ここでQGIS をダウンロードすることができます。
R を GIS として使用することは、「マウス操作」ではなく「コマンドライン・インターフェース」（目的の結果を得るためにはコーディングが必要）であるため、最初は敷居が高く感じられるかもしれません。しかし、繰り返し地図を作成したり、再現性のある分析を行う必要がある場合には、これは大きな利点となります。
空間データ
GIS で使用される空間データには、主にベクトルデータとラスターデータの 2 種類があります。
ベクトルデータ - GIS で使用される空間データの最も一般的な形式です。ベクトルデータは、頂点とパスの幾何学的なフィーチャ (feature; 「地物」とも訳される) で構成されています。ベクトルの空間データは、さらに3つの広く使用されるタイプに分類されます。
点 - 点は、座標系における特定の位置を表す座標ペア(x,y)で構成されます。点は空間データの最も基本的な形態であり、地図上の症例（例：患者の家）や場所（例：病院）を表すのに使用されます。
線 - 線は、2つの接続された点で構成されます。線には長さがあり、道路や川などを表すのに使われます。
ポリゴン - ポリゴンは、点で接続された少なくとも3つの線分で構成されています。ポリゴンのフィーチャは、面積だけでなく、長さ（すなわち、面の外周の長さ）を持っています。ポリゴンは、地域（例：村）や構造（例：病院の実際の面積）を示すときに使用されます。
ラスターデータ - 空間データの代替フォーマットであるラスターデータは、セル（ピクセルなど）のマトリックスで、各セルには高さ、温度、傾斜、森林被覆などの情報が含まれています。これは、航空写真や衛星画像などでよく見られます。ラスターデータは、ベクトルデータの下の「基図」としても使用されます。
空間データの可視化
GIS ソフトウェアでは、空間データを地図上に視覚的に表現するために、それぞれのフィーチャがどこにあるのか、お互いの関係について十分な情報を提供する必要があります。ベクトルデータを使用している場合、ほとんどのユース症例では、この情報は通常シェープファイルに保存されます。
シェープファイル - シェープファイルは、線、点、ポリゴンで構成される「ベクトル」空間データを保存するための一般的なデータフォーマットです。1つのシェープファイルは、実際には少なくとも3つのファイル（.shp、.shx、.dbf）の集合体です。シェープファイルを読み取るためには、これらのサブコンポーネントファイルがすべて同一のディレクトリ（フォルダ）に存在する必要があります。これらの関連ファイルは、ZIPフォルダに圧縮してメールで送信したり、ウェブサイトからダウンロードすることができます。
シェープファイルには、フィーチャ自体の情報だけでなく、地球の表面のどこにあるかという情報も含まれています。空間データをどのように「平面化」するかは、地図の見た目や解釈に大きな影響を与えるため、重要なポイントとなります。
座標参照系 (Coordinate Reference Systems, CRS) - CRS とは、地球上の地理的フィーチャを特定するために使用される座標ベースのシステムです。CRSにはいくつかの重要な要素があります。
座標系 (Coordinate System) - 非常にたくさんの座標系があるため、自分が使用する座標がどの座標系を用いているのかを必ず確認してください。緯度/経度の度数を用いたものは一般的ですが、UTM 座標が利用されていることもあります。
単位 (Units) - 座標系の単位を知ることができます（例：小数点で表された度、メートル）。
測地系 (Datum) - 地球上の特定のモデルの、特定のバージョンを指します。測地系は長年にわたって改訂されてきたので、地図として表示されるレイヤー（マップレイヤー）が同じ測地系を使用していることを確認してください。
投影 (Projection) - 実際には丸い地球を平らな表面（地図）に投影するために使用された数式への参照。
以下に示すマッピングツールを使用しなくても、空間データを要約することができることも覚えておいてください。地域別（例：地区、国など）のシンプルな表があれば十分なこともあります。
訳注: 日本で主に利用されている CRS は次のような ものがあります。
- 平面直角座標系（世界測地系）I〜XIII EPSG: 2443-2455
- JGD2000 GRS80楕円体 EPSG:4612
- JGD2011 GRS80楕円体 EPSG:6668 （東日本大震災以降の大規模地殻変動を反映）
28.3 はじめての GIS
地図を作るためには、以下の通り、いくつかの重要な要素について、考え、また入手しなければありません。
データセット - これは空間データ形式（上述のようにシェープファイルなど）の場合もあれば、空間形式ではない場合もあります（例えば、単なる csv など）。
-
データセットが空間フォーマットでない場合は、参照データも必要です。参照データは、データの空間表現と関連する属性で構成されており、特定のフィーチャの位置情報やアドレス情報を含む資料などがあります。
事前に定義された地理的境界線（例えば、行政区域）を扱う場合、参照シェープファイルは、政府機関やデータ共有組織から自由にダウンロードできることが多いです。迷ったときは、「[地域名] shapefile」とグーグルで検索してみるとよいでしょう。
住所情報はあっても緯度・経度がないような場合は、 ジオコーディング・エンジン（geocoding engine）を使ってデータの空間参照データを取得する必要があるかもしれません。
データセットの情報を対象者に、どのように見せたいかというアイデア。地図には様々な種類がありますが、どの種類の地図が自分のニーズに最も適しているかを考えることが重要です。
訳注：日本の主なGISデータ： * 国勢調査 * 国土数値地図 * 国土地理院基盤地図 * G空間情報センター
データを可視化するための地図のタイプ
色分け（Choropleth）地図 - 主題図の一種で、色、濃淡、またはパターンを使って、ある属性の値に関連して地理的な地域を表現するもの。例えば、大きな値は小さな値よりも濃い色で示されます。このタイプの地図は、ある変数と、それが定義された地域や地政学的エリアでどのように変化するかを視覚化する場合に特に便利です。
症例密度ヒートマップ - 主題図の一種で、値の強さを色で表しますが、データをグループ化するために定義された地域や地政学的な境界線は使用しません。このタイプの地図は、一般的に「ホットスポット」や、点が高密度または集中しているエリアを示すために使用されます。
点密度マップ - は、データの属性値をドットで表現する主題図タイプです。このタイプのマップは、データの散らばりを視覚化し、クラスターを視覚的にスキャンするのに適しています。
等級シンボルマップ（段階のあるシンボルマップ） - 色分け地図に似ている主題図ですが、属性の値を色で表すのではなく、その値に対応するシンボル（通常は円）を使用します。例えば、大きな値は小さな値よりも大きなシンボルで示されます。このタイプの地図は、地域ごとのデータの大きさや量を視覚化したい場合に最適です。
また、複数の異なるタイプの視覚化を組み合わせて、複雑な地理的パターンを示すこともできます。例えば、下の地図の症例（点）は、最も近い医療施設に応じて色分けされています（凡例参照）。大きな赤い円は、一定の半径を持つ医療施設のキャッチメントエリアを示し、赤色の症例点は、キャッチメントの範囲外にある症例を示しています。
注：この GIS の章の主な焦点は、現場でのアウトブレイク対応の状況に基づいています。そのため、この章の内容は、基本的な空間データの操作、視覚化、および分析をカバーしています。
28.4 準備
パッケージを読み込む
以下のコードを実行すると、分析に必要なパッケージが読み込まれます。このハンドブックでは、パッケージを読み込むために、pacman パッケージの p_load() を主に使用しています。p_load() は、必要に応じてパッケージをインストールし、現在の R セッションで使用するためにパッケージを読み込む関数です。また、すでにインストールされたパッケージは、R の基本パッケージである base （以下、base R）の library() を使用して読み込むこともできます。R のパッケージに関する詳細はR の基礎の章をご覧ください。
pacman::p_load(
rio, # データをインポート
here, # ファイルの位置を探す
tidyverse, # データを処理、プロット (ggplot2 パッケージを含む)
sf, # Simple Featureフォーマットによる空間データの管理
tmap, # シンプルな地図を作成、インタラクティブな地図と静的な地図の両方に対応
janitor, # 列名のクリーニング
OpenStreetMap, # ggplot map に OSM 基図を追加
spdep # 空間統計
) CRAN “Spatial Task View”では、空間データを扱うすべての R パッケージの概要を見ることができます。
症例データサンプル
ここでは分かりやすくするため、シミュレーションされたエボラ出血熱の linelist データフレームから1000件のランダムなサンプルを使って作業します（計算上、少ない件数で作業すると、このハンドブックでの表示が容易になります）。続いて、「前処理された」ラインリスト（linelist） をクリックしてダウンロードできます（.rdsファイル）。
ここではランダムにサンプルを取っているので、実際にコードを実行してみると、ここに表示されている結果とは若干異なる結果になるかもしれません。
データを rio パッケージの import() を使ってインポートします（rio パッケージは、.xlsx、.csv、.rds など様々な種類のファイルを扱えるパッケージです。詳細はデータのインポート・エクスポートの章を参照してください）。
# import clean case linelist
linelist <- import("linelist_cleaned.rds") 次に、base R の sample() を使って 1000 行のランダムなサンプルを抽出します。
# ラインリストの行数から1000個のランダムな行番号を生成
sample_rows <- sample(nrow(linelist), 1000)
# サンプル行の全ての列だけを保持するラインリストにサブセットする
linelist <- linelist[sample_rows,]そして、データフレーム型である linelist を “sf”（空間フィーチャ）型のオブジェクトに変換していきます。この ラインリストには、各症例の居住地の経度と緯度を表す “lon” と “lat” の2つの列があるので、これは簡単です。
ここでは、sf（空間的フィーチャ）パッケージとその関数、 st_as_sf() を使用して、linelist_sf という新しいオブジェクトを作成します。この新しいオブジェクトは、基本的には linelist と同じですが、列 lon と lat を座標列として指定し、点を表示する際の座標参照系 (CRS) を割り当てています。4326では、GPS座標の標準である World Geodetic System 1984 (WGS84) に基づいた座標が指定されています。
これがオリジナルの linelist データフレームの中身です。このデモでは、date_onset と geometry（上記の経度と緯度のフィールドから作成されたもので、データフレームの最後の列）の列のみを使用します。
行政境界シェープファイル
シエラレオネ (Sierra Leone): 行政境界シェープファイル
事前に、Humanitarian Data Exchange (HDX) ウェブサイトからシエラレオネのすべての行政（admin）境界データをダウンロードしておきました。あるいは、ハンドブックとデータのダウンロードの章で説明しているように、私たちの R パッケージを介して、このデータやその他のサンプルデータをダウンロードすることもできます。
では、いよいよ Admin Level 3 のシェープファイルを R に保存するために、以下の作業を行います。
- シェープファイルをインポートする
- 列名を加工する
- 興味のある地域だけを残すために行を抽出する
シェープファイルをインポートするには、sf パッケージの read_sf() を使います。ファイルのパスは here() で与えられます。- 今回の場合、ファイルはRプロジェクトの “data”、“gis”、“shp” というサブフォルダに “sle_adm3.shp” として格納されています（詳細は データのインポート・エクスポート と R プロジェクトの設定 の章を参照してください）。独自のファイルパスを設定してください。
次に、janitor パッケージの clean_names() を使って、シェープファイルの列名を標準化します。また、filter() を使って、admin2name が “Western Area Urban” または “Western Area Rural” の行だけを残します。
# ADM3 level clean
sle_adm3 <- sle_adm3_raw %>%
janitor::clean_names() %>% # 列名の標準化
dplyr::filter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # 特定の地域を残してフィルター下の表は、インポートし、クリーニングした後のシェープファイルの様子です。右側にスクロールして、admin レベル 0（国）、admin レベル 1、admin レベル 2、そして最後に admin レベル 3 の列があることがわかります。それぞれのレベルには、キャラクター名と固有の識別子 “pcode” があります。たとえば、SL（Sierra Leone） → 04（Western） → L0410（Western Area Rural）→ L0040101（Koya Rural）のように、pcode は admin レベルが上がるごとに拡張されます。
人口データ
シエラレオネ: ADM3 ごとの人口
これらのデータは、HDX からダウンロードすることもできますし（リンクはこちら）、この章ハンドブックとデータのダウンロードで説明しているように、epirhandbook パッケージを介してダウンロードすることもできます。.csv ファイルの読み込みには、import() を使用します。また、インポートしたファイルを clean_names() に渡して列名の構文を標準化します。
# ADM3 ごとの人口
sle_adm3_pop <- import(here("data", "gis", "population", "sle_admpop_adm3_2020.csv")) %>%
janitor::clean_names()人口ファイルはこのようになっています。右にスクロールすると、各行政区域に male、female、total の列があり、さらに年齢層別の列に人口の内訳が表示されています。
医療施設
シエラレオネ: OpenStreetMap からの医療施設
今回も HDX から医療施設の位置情報をダウンロードしました。こちら または ハンドブックとデータのダウンロードの章の指示に従ってください。
施設の位置を示す点のシェープファイルを read_sf() でインポートし、列名をクリーニングしてから、病院（“hospital”）、診療所（“clinic”）、医師（“doctor”） のいずれかのタグが付けられた点だけを残すように抽出しました。
# OSM 医療施設の shapefile
sle_hf <- sf::read_sf(here("data", "gis", "shp", "sle_hf.shp")) %>%
janitor::clean_names() %>%
dplyr::filter(amenity %in% c("hospital", "clinic", "doctors"))右にスクロールすると、施設名と geometry 座標が表示されます。
28.5 座標のプロット
X-Y座標（経度・緯度で表される点）をプロットする最も簡単な方法は、今回の場合、準備編で作成した linelist_sf オブジェクトから直接、点として描画することです。
tmap パッケージ は、わずか数行のコードで、静的（“plot” モード）とインタラクティブ（“view” モード）の両方に対応したシンプルなマッピング機能を提供します。tmap パッケージ の構文は ggplot2 パッケージ と似ていて、コマンド同士を + で追加していきます。詳しくはこの vignette をご覧ください。
- tmap モードを設定。ここでは plot モードを使い、静的なアウトプットを作成します。
tmap_mode("plot") # "view" または "plot" を選択以下では、点を単独でプロットしています。tm_shape() は、linelist_sf オブジェクトとともに提供されます。次に、tm_dots() でサイズと色を指定して点を追加します。linelist_sf は sf オブジェクトなので、緯度・経度座標と座標系（CRS）を格納する2つの列はすでに指定されています。
# 症例（点）のみ
tm_shape(linelist_sf) + tm_dots(size=0.08, col='blue')
点だけでは、あまり意味がありません。そこで、行政境界もマッピングする必要があります。
ここでも tm_shape() を使いますが（documentation 参照）、症例の点のシェープファイルを提供する代わりに、行政境界のシェープファイル（ポリゴン）を提供します。
bbox = の引数（bbox は “bounding box” の略）で、座標の境界を指定することができます。まず bbox なしの地図表示を行い、次に bbox ありの地図表示を行います。
# 行政境界（ポリゴン）のみ
tm_shape(sle_adm3) + # 行政境界のシェープファイル
tm_polygons(col = "#F7F7F7")+ # ポリゴンを薄い灰色で表示
tm_borders(col = "#000000", # 境界を色と線の太さで表示
lwd = 2) +
tm_text("admin3name") # 各ポリゴンについて表示する列テキスト
# 上と同様、ただしバウンディングボックス引数から縮尺を指定
tm_shape(sle_adm3,
bbox = c(-13.3, 8.43, # 角
-13.2, 8.5)) + # 角
tm_polygons(col = "#F7F7F7") +
tm_borders(col = "#000000", lwd = 2) +
tm_text("admin3name")
そして、今度は点とポリゴンの両方を一緒に表示します。
# すべてまとめる
tm_shape(sle_adm3, bbox = c(-13.3, 8.43, -13.2, 8.5)) + #
tm_polygons(col = "#F7F7F7") +
tm_borders(col = "#000000", lwd = 2) +
tm_text("admin3name")+
tm_shape(linelist_sf) +
tm_dots(size=0.08, col='blue', alpha = 0.5) +
tm_layout(title = "Distribution of Ebola cases") # 地図にタイトルを付ける
R でのマッピングオプションの比較については、こちらの ブログ記事 をご覧ください。
28.6 空間結合
データセットと別のデータセットを 結合 することについては知っているでしょう。いくつかの方法は、このハンドブックのデータの結合の章で説明しています。空間結合は、同様の目的を持ちますが、空間的な関係を利用します。観測値を正しく一致させるために列の共通の値に頼るのではなく、あるフィーチャが他のフィーチャの中に入っているとか、他のフィーチャに最近傍であるとか、他のフィーチャから一定の半径のバッファ内にあるなど、空間的な関係を利用することができます。
sf パッケージには、空間結合のための様々なメソッドが用意されています。st_join メソッドや空間結合の種類については、こちらの参考資料に詳しい説明があります。
ポリゴン内の点
症例に行政単位を空間割当て
ここで興味深い問題があります。症例のラインリストには、症例の admin 単位に関する情報が含まれていません。このような情報は、最初のデータ収集段階で収集するのが理想的ですが、空間的な関係（つまり、点がポリゴンと交差する）に基づいて、個々の症例に admin 単位を割り当てることもできます。
以下では、症例の位置（点）と ADM3 の境界（ポリゴン）を空間的に交差させます。
- linelist（点）から始める
- 結合のタイプを “st_intersects” に設定して、境界線に空間的に結合する
-
select()を使用して、新しい行政境界列のうち特定のものだけを残す
linelist_adm <- linelist_sf %>%
# 空間交差に基づいて行政境界を linelist に結合
sf::st_join(sle_adm3, join = st_intersects)sle_adms のすべての列がラインリストに追加されました！各症例には、該当する行政区域レベルの詳細を示す列が追加されました。この例では、新しい列のうち 2 つ（admin レベル3）だけを残したいので、古い列名を select() して、追加したい2つだけを選択します。
linelist_adm <- linelist_sf %>%
# 空間交差に基づいて、行政境界ファイルをラインリストに結合
sf::st_join(sle_adm3, join = st_intersects) %>%
# 古い列名はそのままに、新たに2つの admin の列名を追加
select(names(linelist_sf), admin3name, admin3pcod)下の図は、最初の10件と、ポリゴン図形と空間的に交差する点に基づいて付けられた admin レベル3（ADM3）の行政境界を表示しています。
## Simple feature collection with 1000 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -13.27122 ymin: 8.447961 xmax: -13.20576 ymax: 8.490898
## Geodetic CRS: WGS 84
## First 10 features:
## case_id admin3name admin3pcod geometry
## 3546 7cfdb0 Mountain Rural SL040102 POINT (-13.211 8.466934)
## 237 ee1b06 West III SL040208 POINT (-13.25373 8.462149)
## 2334 5da423 West I SL040206 POINT (-13.2454 8.484667)
## 2164 4360db East I SL040203 POINT (-13.21523 8.485395)
## 636 6c6bc5 East III SL040205 POINT (-13.20848 8.464067)
## 2602 2759d7 Central I SL040201 POINT (-13.23106 8.483106)
## 4054 cfd58c East II SL040204 POINT (-13.21283 8.476574)
## 4091 86f641 West III SL040208 POINT (-13.2701 8.451366)
## 1684 61751f West II SL040207 POINT (-13.23935 8.474191)
## 4686 576d23 East I SL040203 POINT (-13.21303 8.486304)空間結合の前にはできなかった、行政単位での症例の表現が可能になりました。
# 行政単位ごとに症例の数を含むデータフレームを新規作成
case_adm3 <- linelist_adm %>% # 新しい admin 列を含むラインリストで始める
as_tibble() %>% # 見やすいように tibble に変換
group_by(admin3pcod, admin3name) %>% # 名前と pcode で行政単位をグループ化
summarise(cases = n()) %>% # 要約と行数のカウント
arrange(desc(cases)) # 下り順に並べ替え
case_adm3## # A tibble: 10 × 3
## # Groups: admin3pcod [10]
## admin3pcod admin3name cases
## <chr> <chr> <int>
## 1 SL040102 Mountain Rural 285
## 2 SL040208 West III 228
## 3 SL040207 West II 184
## 4 SL040204 East II 111
## 5 SL040203 East I 58
## 6 SL040201 Central I 52
## 7 SL040206 West I 39
## 8 SL040205 East III 26
## 9 SL040202 Central II 15
## 10 <NA> <NA> 2また、行政単位ごとの症例数の棒グラフを作成することもできます。
この例では、ggplot() を linelist_adm で始めているので、頻度で棒グラフを並べる fct_infreq() などの因子型のための関数を適用することができます（ヒントは 因子（ファクタ）型データ の章を参照してください）。
ggplot(
data = linelist_adm, # admin unit 情報を含むラインリストで始める
mapping = aes(
x = fct_rev(fct_infreq(admin3name))))+ # x軸は行政単位、件数で並べ替え
geom_bar()+ # 帽を作成、高さは行数
coord_flip()+ # adm を読みやすくするためにXとYを入れ替える
theme_classic()+ # 背景を単純に
labs( # タイトルとラベル
x = "Admin level 3",
y = "Number of cases",
title = "Number of cases, by adminstative unit",
caption = "As determined by a spatial join, from 1000 randomly sampled cases from linelist"
)
最近傍
最寄の医療施設／集客エリアを探す
病気のホットスポットに関連して、医療施設がどこにあるかを知っておくと便利かもしれません。
ここでは、st_join() (sf パッケージ)の st_nearest_feature 結合メソッドを使って、個々の症例に最も近い医療施設を可視化することができます。
- まず、シェープファイルのラインリスト
linelist_sfを用意します。 - 保健施設や診療所の位置（点）である
sle_hfと空間的に結合します。
# 各症例に最も近い医療施設
linelist_sf_hf <- linelist_sf %>% # linelist shapefile で始める
st_join(sle_hf, join = st_nearest_feature) %>% # 症例データから最も近い診療所からのデータ
select(case_id, osm_id, name, amenity) %>% # 残しておくべき列、例えば id, name, type, と医療施設の位置情報
rename("nearest_clinic" = "name") # 分かりやすいように名前を変更以下（最初の50行）のように、それぞれの症例には最寄りの診療所や病院のデータがあることがわかります。
約30％の症例で “Den Clinic” が最も近い医療機関であることがわかります。
# 医療施設ごとに症例を数える
hf_catchment <- linelist_sf_hf %>% # 最寄りの診療所データを含むラインリストで始める
as.data.frame() %>% # shapefile をデータフレームに変換
count(nearest_clinic, # （診療所の）"name" で行を数える
name = "case_n") %>% # 数えたデータの列に "case_n" と命名
arrange(desc(case_n)) # 下がり順で並べ替え
hf_catchment # console に表示## nearest_clinic case_n
## 1 Den Clinic 361
## 2 Shriners Hospitals for Children 321
## 3 GINER HALL COMMUNITY HOSPITAL 191
## 4 panasonic 43
## 5 Princess Christian Maternity Hospital 33
## 6 ARAB EGYPT CLINIC 23
## 7 MABELL HEALTH CENTER 14
## 8 <NA> 14結果を可視化するために、tmap パッケージ を使用することができます。
tmap_mode("view") # tmap モードをインタラクティブに設定
# 症例と診療所の点をプロット
tm_shape(linelist_sf_hf) + # 症例をプロット
tm_dots(size=0.08, # 最も近い診療所で症例を色分け
col='nearest_clinic') +
tm_shape(sle_hf) + # 診療所を大きい黒い点でプロット
tm_dots(size=0.3, col='black', alpha = 0.4) +
tm_text("name") + # 施設名でオーバーレイ
tm_view(set.view = c(-13.2284, 8.4699, 13), # 縮尺を調整 (中心座標, zoom)
set.zoom.limits = c(13,14))+
tm_layout(title = "Cases, colored by nearest clinic")バッファ
また、最も近い医療施設から徒歩2.5km（約30分）以内にある症例がどれくらいあるかを調べることもできます。
注: より正確な距離を計算するためには、UTM（地球を平面に投影したもの）などの各地域の地図投影系に sf オブジェクトを再投影するのがよいでしょう。この例では、簡単にするために、世界測地系（WGS84）の地理的座標系（地球は球状/円形の表面で表現され、そのため単位は10進法の度数）にしておきます。ここでは、一般的な換算として、1度＝約111kmとします。
地図投影法と座標系については、こちらのesriの記事をご覧ください。こちらのブログ では、様々なタイプの地図投影について、また、興味のある分野や地図・分析の文脈に応じて、どのように適切な投影を選ぶことができるかについて述べられています。
まず、各医療施設の周囲に半径約2.5kmの円形バッファを作成します。これは tmap パッケージ の st_buffer() で行います。地図の単位は緯度・経度の10進法であるため、“0.02” は度数として解釈されます。地図の座標系がメートル単位の場合は、数値もメートル単位で指定する必要があります。
sle_hf_2k <- sle_hf %>%
st_buffer(dist=0.02) # 約2.5kmに度数を変換以下のように、バッファゾーン自体を作ります：
tmap_mode("plot")
# 円バッファを作成
tm_shape(sle_hf_2k) +
tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) + # 診療所を大きい赤丸でプロット
tm_dots(size=0.3, col='black') 
次に、st_join() の st_intersects 結合タイプを使って、これらのバッファと症例（点）を交差させます。これによって、バッファからのデータが、交差している点に結合されます。
# バッファで症例を交差
linelist_sf_hf_2k <- linelist_sf_hf %>%
st_join(sle_hf_2k, join = st_intersects, left = TRUE) %>%
filter(osm_id.x==osm_id.y | is.na(osm_id.y)) %>%
select(case_id, osm_id.x, nearest_clinic, amenity.x, osm_id.y)これによって結果を数えることができます。1000件のうち1000件は、どのバッファとも交差していない（値が欠けている）ので、最寄りの医療施設から徒歩30分以上の場所に住んでいることになります。
## [1] 1000その結果、どのバッファとも交わらなかった症例が赤で表示されるように可視化されます。
tmap_mode("view")
# まず症例を点で表示
tm_shape(linelist_sf_hf) +
tm_dots(size=0.08, col='nearest_clinic') +
# 診療所を大きい黒点でプロット
tm_shape(sle_hf) +
tm_dots(size=0.3, col='black')+
# 医療施設バッファをポリラインで重ねる
tm_shape(sle_hf_2k) +
tm_borders(col = "black", lwd = 2) +
# どの医療施設バッファにもない症例を赤点で強調
tm_shape(linelist_sf_hf_2k %>% filter(is.na(osm_id.y))) +
tm_dots(size=0.1, col='red') +
tm_view(set.view = c(-13.2284,8.4699, 13), set.zoom.limits = c(13,14))+
# タイトルを追加
tm_layout(title = "Cases by clinic catchment area")他の空間結合
joinの引数は、このほかに以下の値をとります。 (ドキュメントより)
- st_contains_properly
- st_contains
- st_covered_by
- st_covers
- st_crosses
- st_disjoint
- st_equals_exact
- st_equals
- st_is_within_distance
- st_nearest_feature
- st_overlaps
- st_touches
- st_within
28.7 色分け（Choropleth）地図
色分け地図は、あらかじめ定義されたエリア（通常は行政単位や保健エリア）ごとにデータを視覚化するのに役立ちます。例えば、感染発生対応では、発生率の高い特定の地域にリソースを割り当てる際に役立ちます。
すべての症例に行政単位名が割り当てられたので（上記の「空間結合」のセクションを参照）、症例数をエリア別にマッピングすることができます（色分け地図）。
ADM3 ごとの人口データもあるので、この情報を先に作成した case_adm3 テーブルに追加することができます。
まず、前のステップで作成したデータフレーム case_adm3 から始めます。これは、各行政単位とその症例数の要約表です。
- 人口データ
sle_adm3_popは、dplyr パッケージ のleft_join()を用いて、case_adm3データフレームのadmin3pcod列とsle_adm3_popデータフレームのadm_pcode列に共通する値に基づいて結合されています。データの結合の章を参照）。 -
select()を新しいデータフレームに適用して、有用な列だけを残す -totalは総人口です。 -
mutate()を用いて、人口 10,000 人あたりの症例数を新しい列として計算します。
# 人口データを追加し、１万人当たりの症例数を計算
case_adm3 <- case_adm3 %>%
left_join(sle_adm3_pop, # 人口データから列を追加
by = c("admin3pcod" = "adm3_pcode")) %>% # 二つの列の共通の値に基づく結合
select(names(case_adm3), total) %>% # 重要な列のみ保持、総人口など
mutate(case_10kpop = round(cases/total * 10000, 3)) # 10000あたりの症例数の列を作成、小数点以下３桁で四捨五入
case_adm3 # console に表示## # A tibble: 10 × 5
## # Groups: admin3pcod [10]
## admin3pcod admin3name cases total case_10kpop
## <chr> <chr> <int> <int> <dbl>
## 1 SL040102 Mountain Rural 285 33993 83.8
## 2 SL040208 West III 228 210252 10.8
## 3 SL040207 West II 184 145109 12.7
## 4 SL040204 East II 111 99821 11.1
## 5 SL040203 East I 58 68284 8.49
## 6 SL040201 Central I 52 69683 7.46
## 7 SL040206 West I 39 60186 6.48
## 8 SL040205 East III 26 500134 0.52
## 9 SL040202 Central II 15 23874 6.28
## 10 <NA> <NA> 2 NA NAこのテーブルを ADM3 のポリゴンシェープファイル と結合してマッピングします。
case_adm3_sf <- case_adm3 %>% # 行政単位での症例と感染率で始める
left_join(sle_adm3, by="admin3pcod") %>% # shapefile データと共通列で結合
select(objectid, admin3pcod, # 指定した列のみ保持
admin3name = admin3name.x, # 一つの列名をきれいにする
admin2name, admin1name,
cases, total, case_10kpop,
geometry) %>% # プロットできるように座標を保持
drop_na(objectid) %>% # NA 行を取り除く
st_as_sf() # shapefile に変換結果をマッピング。
case_adm3_sf <- na.omit(case_adm3_sf)
# tmap mode
tmap_mode("plot") # 静的地図を表示
# ポリゴンをプロット
tm_shape(case_adm3_sf) +
tm_polygons("cases") + # 症例数の列で色分け
tm_text("admin3name") # 表示に名前を付ける
インシデンス率もまた地図にします。
# 1万人あたりの症例
tmap_mode("plot") # 静的ビューモード
# プロット
tm_shape(case_adm3_sf) + # ポリゴンをプロット
tm_polygons("case_10kpop", # 症例率を含む列で色分け
breaks=c(0, 10, 50, 100), # 色分けの値を定義
palette = "Purples" # 紫色のカラーパレットを使用
) +
tm_text("admin3name") # テキストを表示
28.8 ggplot2 で地図作成
すでに ggplot2 パッケージ の使用に慣れている場合は、データの静的なマップを作成する代わりに、ggplot2 パッケージを使用することができます。geom_sf() は、データに含まれる特徴 (点、線、多角形) に応じて異なるオブジェクトを描画します。例えば、ポリゴンを含む sf データを使った ggplot() の中で geom_sf() を使うと、色分け地図を作ることができます。
これがどのように機能するかを説明するために、先ほど使用した ADM3 ポリゴンのシェープファイルから始めます。このポリゴンは、シエラレオネの Admin Level 3 の地域であることを思い出してください。
sle_adm3## Simple feature collection with 12 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13.29894 ymin: 8.094272 xmax: -12.91333 ymax: 8.499809
## Geodetic CRS: WGS 84
## # A tibble: 12 × 20
## objectid admin3name admin3pcod admin3ref_n admin2name admin2pcod admin1name
## * <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 155 Koya Rural SL040101 Koya Rural Western A… SL0401 Western
## 2 156 Mountain Ru… SL040102 Mountain R… Western A… SL0401 Western
## 3 157 Waterloo Ru… SL040103 Waterloo R… Western A… SL0401 Western
## 4 158 York Rural SL040104 York Rural Western A… SL0401 Western
## 5 159 Central I SL040201 Central I Western A… SL0402 Western
## 6 160 East I SL040203 East I Western A… SL0402 Western
## 7 161 East II SL040204 East II Western A… SL0402 Western
## 8 162 Central II SL040202 Central II Western A… SL0402 Western
## 9 163 West III SL040208 West III Western A… SL0402 Western
## 10 164 West I SL040206 West I Western A… SL0402 Western
## 11 165 West II SL040207 West II Western A… SL0402 Western
## 12 167 East III SL040205 East III Western A… SL0402 Western
## # ℹ 13 more variables: admin1pcod <chr>, admin0name <chr>, admin0pcod <chr>,
## # date <date>, valid_on <date>, valid_to <date>, shape_leng <dbl>,
## # shape_area <dbl>, rowcacode0 <chr>, rowcacode1 <chr>, rowcacode2 <chr>,
## # rowcacode3 <chr>, geometry <MULTIPOLYGON [°]>dplyr パッケージ の left_join() を使って、マッピングしたいデータをシェープファイルオブジェクトに追加することができます。この例では、先ほど作成した case_adm3 データフレームを使って、行政区ごとの症例数をまとめていますが、データフレームに格納されているどのようなデータでも、同様の方法でマッピングすることができます。
sle_adm3_dat <- sle_adm3 %>%
inner_join(case_adm3, by = "admin3pcod") # inner join = どちらのデータオブジェクトにもある場合にのみ残す
select(sle_adm3_dat, admin3name.x, cases) # console に選択した変数を表示## Simple feature collection with 9 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13.29894 ymin: 8.384533 xmax: -13.12612 ymax: 8.499809
## Geodetic CRS: WGS 84
## # A tibble: 9 × 3
## admin3name.x cases geometry
## <chr> <int> <MULTIPOLYGON [°]>
## 1 Mountain Rural 285 (((-13.21496 8.474341, -13.21479 8.474289, -13.21465 8.4…
## 2 Central I 52 (((-13.22646 8.489716, -13.22648 8.48955, -13.22644 8.48…
## 3 East I 58 (((-13.2129 8.494033, -13.21076 8.494026, -13.21013 8.49…
## 4 East II 111 (((-13.22653 8.491883, -13.22647 8.491853, -13.22642 8.4…
## 5 Central II 15 (((-13.23154 8.491768, -13.23141 8.491566, -13.23144 8.4…
## 6 West III 228 (((-13.28529 8.497354, -13.28456 8.496497, -13.28403 8.4…
## 7 West I 39 (((-13.24677 8.493453, -13.24669 8.493285, -13.2464 8.49…
## 8 West II 184 (((-13.25698 8.485518, -13.25685 8.485501, -13.25668 8.4…
## 9 East III 26 (((-13.20465 8.485758, -13.20461 8.485698, -13.20449 8.4…ggplot2 パッケージ を使って、地域別の症例数の棒グラフを作成するには、次のように geom_col() を呼び出します。
ggplot(data=sle_adm3_dat) +
geom_col(aes(x=fct_reorder(admin3name.x, cases, .desc=T), # 'cases' 下がり順でx軸を並べ替え
y=cases)) + # y軸は地域ごとの症例数
theme_bw() +
labs( # 図のテキストを設定
title="Number of cases, by administrative unit",
x="Admin level 3",
y="Number of cases"
) +
guides(x=guide_axis(angle=45)) # 見やすいようにx軸を45度傾ける
もし、ggplot2 パッケージ を使って症例数の色分け地図を作りたい場合は、同様の構文で geom_sf() を呼び出すことができます。
そして、例えば、ggplot2 パッケージ で統一されている文法を使って、地図の外観をカスタマイズすることができます。
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=cases)) +
scale_fill_continuous(high="#54278f", low="#f2f0f7") + # 色の段階を変更
theme_bw() +
labs(title = "Number of cases, by administrative unit", # 図のテキストを設定
subtitle = "Admin level 3"
)
ggplot2 パッケージ に慣れているRユーザーにとって、geom_sf()は、実装がシンプルで直接的であり、基本的な地図の可視化に適しています。詳しくは、geom_sf() ビニエット（vignette） または ggplot2 book をご覧ください。
28.9 基図
OpenStreetMap
以下では、OpenStreetMap の機能を使って ggplot2 パッケージを利用したマップの基図を作成する方法を説明します。他の方法としては、Google への無料登録が必要な ggmap パッケージ を使用する方法があります(詳細)。
OpenStreetMap パッケージは、自由に編集可能な世界地図を作成する共同プロジェクトです。基本となるジオロケーション・データ（都市、道路、自然、空港、学校、病院、道路などの位置情報）は、このプロジェクトの主要なアウトプットと考えられています。
まず、OpenStreetMap パッケージをロードして、そこから基図を取得します。
次に、オブジェクト map を作成します。これは、OpenStreetMap パッケージの openmap() を使用して定義します (documentation)。以下のものを用意します。
-
upperLeftとlowerRight基図タイルの境界を指定する2つの座標ペア。- ここでは、ラインリストの行の最大値と最小値を入れているので、地図は動的にデータに反応します。
-
zoom =(null の場合は自動的に決定されます)
-
type =基図の種類 - ここではいくつかの可能性を挙げていますが、コードは現在、最初のもの（[1]）である “osm” を使用しています。 -
mergeTiles =TRUEを選択したので、ベースタイルはすべて1つに統合されます。
# パッケージをロード
pacman::p_load(OpenStreetMap)
# 緯度経度の範囲に基図を切り取る。タイル種別を選択。
map <- OpenStreetMap::openmap(
upperLeft = c(max(linelist$lat, na.rm=T), max(linelist$lon, na.rm=T)), # 基図のタイルの制限
lowerRight = c(min(linelist$lat, na.rm=T), min(linelist$lon, na.rm=T)),
zoom = NULL,
type = c("osm", "stamen-toner", "stamen-terrain", "stamen-watercolor", "esri","esri-topo")[1])今、この基図を OpenStreetMap パッケージの autoplot.OpenStreetMap() を使ってプロットしてみると、軸の単位が緯度・経度座標ではない違う座標系を使っていることがわかります。緯度/経度で保存されている）症例の位置を正しく表示するには、これを変更する必要があります。
autoplot.OpenStreetMap(map)そこで、 OpenStreetMap パッケージの openproj() を使って、地図を緯度/経度に変換したいと思います。基図 map と、必要な座標系（CRS）を指定します。ここでは、WGS 1984投影の “proj.4”文字列を指定していますが、他の方法でCRSを指定することもできます。(proj.4文字列が何であるかを理解するには、このページを参照してください)
# 座標系 WGS84
map_latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")プロットを作成すると、軸に沿って緯度と経度の座標が表示されています。座標系が変換されたことになります。これで、症例を重ね合わせても正しくプロットできるようになりました。
# 地図をプロット。ggplot を使うためには、"autoplot" が必須。
autoplot.OpenStreetMap(map_latlon)
28.10 等高線密度ヒートマップ
以下では、基図の上に症例の等高線付き密度ヒートマップを作成する方法を説明します。まず、ラインリストから始めます（1症例につき1行）。
- 上述のように OpenStreetMap から基図タイルを作成する
-
linelistの症例を、緯度と経度の列を使ってプロットする。 -
ggplot2 パッケージの
stat_density_2d()を使って、点を密度ヒートマップに変換する。
緯度・経度座標を持つ基図があれば、その上に事例の居住地の緯度・経度座標を使ってプロットすることができます。
基図を作成するための関数 autoplot.OpenStreetMap() をベースにして、ggplot2 パッケージ の関数を使えば、以下の geom_point() のように簡単に上に追加することができます。
# 地図をプロット。ggplot と使うためには autoplot が必須。
autoplot.OpenStreetMap(map_latlon)+ # 基図から始める
geom_point( # ラインリストの経度と緯度の列から xy の点を追加
data = linelist,
aes(x = lon, y = lat),
size = 1,
alpha = 0.5,
show.legend = FALSE) + # 凡例を完全に削除
labs(x = "Longitude", # タイトルとラベル
y = "Latitude",
title = "Cumulative cases")
上のマップは、特に点が重なっているところは、解釈が難しいかもしれません。そこで、 ggplot2 パッケージの stat_density_2d() を用いて 2 次元密度マップを作成することができます。ラインリストの緯度/軽度座標を使用していますが、2D カーネル密度推定を行い、その結果を等高線で表示します（地形図のように）。完全なドキュメントはこちらをご覧ください。
# 基図から始める
autoplot.OpenStreetMap(map_latlon)+
# 密度プロットを追加
ggplot2::stat_density_2d(
data = linelist,
aes(
x = lon,
y = lat,
fill = ..level..,
alpha = ..level..),
bins = 10,
geom = "polygon",
contour_var = "count",
show.legend = F) +
# カラースケールを指定
scale_fill_gradient(low = "black", high = "red")+
# ラベル
labs(x = "Longitude",
y = "Latitude",
title = "Distribution of cumulative cases")
時系列ヒートマップ
上の密度ヒートマップは、累積症例を示しています。このヒートマップをラインリストから得られた症状の発症月に基づいて分割 （facet） することで、時空間的なアウトブレイクを調べることができます。
まず、linelistで、発症した年と月の新しい列を作成します。base R の format() は、日付の表示方法を変更します。ここでは “YYYY-MM”とします。
# 発症した月を抽出
linelist <- linelist %>%
mutate(date_onset_ym = format(date_onset, "%Y-%m"))
# 値を調査
table(linelist$date_onset_ym, useNA = "always")##
## 2014-04 2014-05 2014-06 2014-07 2014-08 2014-09 2014-10 2014-11 2014-12 2015-01
## 2 10 13 36 73 204 199 123 91 71
## 2015-02 2015-03 2015-04 <NA>
## 55 44 34 45ここでは、密度ヒートマップに ggplot2 パッケージによる分割列を導入するだけです。facet_wrap() が適用され、新しい列が行として使用されます。わかりやすくするために、分割列の数を 3 に設定しました。
# パッケージ
pacman::p_load(OpenStreetMap, tidyverse)
# 基図から始める
autoplot.OpenStreetMap(map_latlon)+
# 密度プロットを追加
ggplot2::stat_density_2d(
data = linelist,
aes(
x = lon,
y = lat,
fill = ..level..,
alpha = ..level..),
bins = 10,
geom = "polygon",
contour_var = "count",
show.legend = F) +
# カラースケールを指定
scale_fill_gradient(low = "black", high = "red")+
# ラベル
labs(x = "Longitude",
y = "Latitude",
title = "Distribution of cumulative cases over time")+
# 発症した month-year でプロットを分割
facet_wrap(~ date_onset_ym, ncol = 4)
28.11 空間統計
これまでの議論のほとんどは、空間データの視覚化に焦点を当ててきました。しかし、空間統計を使用し、データ内の属性の空間的な関係を定量化したいような場合もあるでしょう。ここでは、空間統計の主要な概念を簡単に説明し、より包括的な空間分析を行いたい場合に参考となるリソースを紹介します。
空間的な関係
空間統計を計算する前に、データのフィーチャ間の関係を特定する必要があります。空間的な関係を概念化する方法は数多くありますが、シンプルで一般的に適用できるモデルは「隣接関係」です。具体的には、境界を共有している、または「隣り合っている」エリア間に地理的な関係があると考えます。
spdep パッケージで使用している sle_adm3 データの行政区域ポリゴン間の隣接関係を定量化することができます。ここでは、共点隣接 (queen contiguity) を指定します。こちらは、地域がその境界に沿って少なくとも1つの点を共有していれば、隣り合っていることを意味します。この他には、共辺隣接 (rook) を指定することもできます。こちらは、地域がその境界に沿って境界縁を共有していれば、隣り合っていることを意味します。この違いは些細なことですが、この事例は不規則なポリゴンですので、場合によっては queen と rook の選択が影響を与えることもあります。
sle_nb <- spdep::poly2nb(sle_adm3_dat, queen=T) # 隣接を作成
sle_adjmat <- spdep::nb2mat(sle_nb) # 隣接関係をまとめた行列の作成
sle_listw <- spdep::nb2listw(sle_nb) # listw (重みのリスト) オブジェクトを作成 -- 後で使います
sle_nb## Neighbour list object:
## Number of regions: 9
## Number of nonzero links: 30
## Percentage nonzero weights: 37.03704
## Average number of links: 3.333333
round(sle_adjmat, digits = 2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## 1 0.00 0.20 0.00 0.20 0.00 0.2 0.00 0.20 0.20
## 2 0.25 0.00 0.00 0.25 0.25 0.0 0.00 0.25 0.00
## 3 0.00 0.00 0.00 0.50 0.00 0.0 0.00 0.00 0.50
## 4 0.25 0.25 0.25 0.00 0.00 0.0 0.00 0.00 0.25
## 5 0.00 0.33 0.00 0.00 0.00 0.0 0.33 0.33 0.00
## 6 0.50 0.00 0.00 0.00 0.00 0.0 0.00 0.50 0.00
## 7 0.00 0.00 0.00 0.00 0.50 0.0 0.00 0.50 0.00
## 8 0.20 0.20 0.00 0.00 0.20 0.2 0.20 0.00 0.00
## 9 0.33 0.00 0.33 0.33 0.00 0.0 0.00 0.00 0.00
## attr(,"call")
## spdep::nb2mat(neighbours = sle_nb)上の図は、sle_adm3 データの 9 つのリージョン間の関係を示しています。スコア 0 は2つの地域が隣り合っていないことを示し、0以外の値は隣り合っていることを示しています。行列の値は、各地域の行の重みの合計が 1 になるようにスケーリングされています。
これらの隣接関係を視覚化するより良い方法は、プロットすることです。
plot(sle_adm3_dat$geometry) + # 地域境界をプロット
spdep::plot.nb(sle_nb,as(sle_adm3_dat, 'Spatial'), col='grey', add=T) # 隣接関係を追加
ここでは、隣接するポリゴンを特定するために、隣接アプローチを用いました。特定した隣接は、隣接ベース近傍 (contiguity-based neighbors) と呼ばれることもあります。しかし、これは、どの地域が地理的な関係を持っていると予想されるかを選択する一つの方法に過ぎません。地理的関係を特定するため別の方法で、最も一般的な方法は、以下のような距離ベース近傍を生成することです。
空間的自己相関
よく引用されるトブラーの地理学の第一法則では、「すべてのものは他のすべてのものに関係しているが、近くのものは遠くのものより関係している」とされています。疫学においては、ある地域における特定の健康結果のリスクは、遠くの地域よりも近隣の地域に類似していることを意味していることが多いです。この概念は、「空間的自己相関」として正式に定義されています。これは、類似した値を持つ地理的フィーチャが空間的に集まっているという統計的特性です。空間的自己相関の統計的測定は、データにおける空間的クラスタリングの程度を定量化し、クラスタリングが発生する場所を特定し、データ内の異なる変数間の空間的自己相関の共通パターンを特定するために使用することができます。このセクションでは、空間的自己相関の一般的な測定方法とRでの計算方法について説明します。
Moran’s I - これは、ある地域での変数の値と、近隣の地域での同じ変数の値との間の相関のグローバルな要約統計です。通常、Moran’s I 統計量は、-1 から 1 の範囲をとります。0 は、空間的な相関のパターンがないことを示し、1 や -1 に近い値は、それぞれ空間的な自己相関（似たような値が近くにある）や空間的なばらつき（似ていない値が近くにある）が強いことを示しています。
例として、先ほどマッピングしたエボラ出血熱の症例の空間的自己相関を定量化するために、Moran’s I 統計量を計算してみます（これは、シミュレーションされた linelist データフレームからの症例のサブセットであることを覚えておいてください）。spdep パッケージには、この計算を行うための、moran.test() があります。
moran_i <-spdep::moran.test(sle_adm3_dat$cases, # 指定の変数で数値ベクトル
listw=sle_listw) # 隣接関係を要約したlistwオブジェクト
moran_i # Moran's I 検定の結果を表示##
## Moran I test under randomisation
##
## data: sle_adm3_dat$cases
## weights: sle_listw
##
## Moran I statistic standard deviate = 1.6045, p-value = 0.0543
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.21087821 -0.12500000 0.04382128moran.test() の出力を見ると、Moran’s I 統計量は round(moran_i$estimate[1],2) となっています。これは、データに空間的な自己相関が存在することを示しています。具体的には、エボラ出血熱の患者数が同程度の地域は、近くにある可能性が高いということです。moran.test() で得られるp値は、空間的自己相関がないという帰無仮説の下での期待値との比較によって生成され、正式な仮説検定の結果を報告する必要がある場合に使用できます。
Local Moran’s I - 上で計算した（グローバルな）Moran’s I 統計を分解することで、局所的な空間的な自己相関を識別することができます。これを利用すると、データ内の特定のクラスターの識別を行うことができます。この統計は、LISA（Local Indicator of Spatial Association）統計量と呼ばれることもあり、個々の地域における空間的自己相関の程度を要約したものです。この値は、地図上の「ホット」と「コールド」な場所を見つけるのに役立ちます。
例として、上で使用したエボラ出血熱の症例数に対して、spdep パッケージの local_moran() を用いて、Local Moran’s I を計算し、地図上に描画することができます。
# local Moran's I を計算
local_moran <- spdep::localmoran(
sle_adm3_dat$cases, # 指定の変数
listw=sle_listw # 隣接の重み付けをした listw オブジェクト
)
# 結果を sf データに結合
sle_adm3_dat<- cbind(sle_adm3_dat, local_moran)
# 地図をプロット
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=Ii)) +
theme_bw() +
scale_fill_gradient2(low="#2c7bb6", mid="#ffffbf", high="#d7191c",
name="Local Moran's I") +
labs(title="Local Moran's I statistic for Ebola cases",
subtitle="Admin level 3 regions, Sierra Leone")
Getis-Ord Gi* - これもホットスポットの分析によく使われる統計です。この統計が人気である理由は、ArcGIS のホットスポット分析ツールで使われていることに関係しています。この統計は、通常、近隣地域間の変数の値の差は正規分布に従うはずだという仮定に基づいています。このツールは、z-score を利用したアプローチを用いて、指定された変数の値が近隣の地域に比べて統計的に特定の変数の値が有意に高い（ホットスポット）または低い（コールドスポット）地域を特定します。
spdep パッケージの localG() を使って Gi* 統計量を計算し、マッピングすることができます。
# local G 分析を実行
getis_ord <- spdep::localG(
sle_adm3_dat$cases,
sle_listw
)
# 結果を sf データに結合
sle_adm3_dat$getis_ord <- as.numeric(getis_ord)
# 地図をプロット
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=getis_ord)) +
theme_bw() +
scale_fill_gradient2(low="#2c7bb6", mid="#ffffbf", high="#d7191c",
name="Gi*") +
labs(title="Getis-Ord Gi* statistic for Ebola cases",
subtitle="Admin level 3 regions, Sierra Leone")
上のとおり、Getis-Ord Gi* のマップは、先に作成した Local Moran’s の地図とは若干異なっています。これは、これらの2つの統計量を計算するために使用される方法がわずかに異なることを反映しています。どちらを使用すべきかは、特定の事例とあなたが関心を持っているリサーチクエスチョンに依存します。
Lee’s L test - これは、二変量の空間相関に関する統計的な検定です。与えられた変数 x の空間パターンが、x に空間的に関連すると仮定された別の変数 y の空間パターンに類似しているかどうかを検定することができます。
例として、シミュレーションしたエボラ出血熱の患者数の空間パターンが、人口の空間パターンと相関しているかどうかを検証してみましょう。まず始めに、sle_adm3 データに population 変数が必要です。先ほど読み込んだ sle_adm3_pop データフレームの total 変数を使うことができます。
2つの変数の空間パターンを並べて視覚化することで、似ているかどうかをすぐに確認することができます。
tmap_mode("plot")
cases_map <- tm_shape(sle_adm3_dat) + tm_polygons("cases") +
tm_layout(main.title="Cases")
pop_map <- tm_shape(sle_adm3_dat) +
tm_polygons("population") +
tm_layout(main.title="Population")
tmap_arrange(cases_map, pop_map, ncol=2) # arrange into 2x1 facets
視覚的には、この 2 つのパターンは似ていないように見えます。spdep パッケージの lee.test() を使って、2 つの変数の空間的自己相関のパターンが関連しているかどうかを統計的に検定することができます。L 統計量は、パターン間に相関がなければ 0 に近く、強い正の相関があれば（パターンが似ている場合は）1 に近く、強い負の相関があれば（パターンが逆の場合は）-1 に近くなります。
lee_test <- spdep::lee.test(
x=sle_adm3_dat$cases, # 比較する変数１
y=sle_adm3_dat$population, # 比較する変数２
listw=sle_listw # 隣接重みづけのある listw オブジェクト
)
lee_test##
## Lee's L statistic randomisation
##
## data: sle_adm3_dat$cases , sle_adm3_dat$population
## weights: sle_listw
##
## Lee's L statistic standard deviate = -0.94983, p-value = 0.8289
## alternative hypothesis: greater
## sample estimates:
## Lee's L statistic Expectation Variance
## -0.13999151 -0.03451068 0.01233270上記の出力では、2つの変数の Lee’s L 統計量が round(lee_test$estimate[1],2) となり、弱い負の相関があることを示しています。これは、症例のパターンと人口が互いに関連していないという私たちの視覚的な評価を裏付けるものであり、症例の空間的なパターンは、厳密には高リスク地域の人口密度の結果ではないという証拠となります。
Lee’s L 統計は、空間的に分布する変数間の関係について、このような推論を行うのに役立ちます。しかし、2 つの変数間の関係の性質をより詳細に記述したり、交絡を調整したりするには、空間回帰技術が必要になります。以下では、これらについて簡単に説明します。
空間回帰
空間データの変数間の関係を統計的に推論したい場合があります。このような場合には、空間回帰の手法、つまり、データ内のユニットの空間的な構成を明示的に考慮した回帰のアプローチを検討することが有用です。GLM のような標準的な回帰モデルではなく、空間回帰モデルを検討する必要がある理由は以下の通りです。
標準的な回帰モデルは、残差が互いに独立していることを前提としています。標準回帰モデルは、残差が互いに独立であると仮定していますが、強い空間的自己相関がある場合、標準回帰モデルの残差は空間的にも自己相関している可能性が高く、この仮定に反することになります。これは、モデルの結果を解釈する際に問題となる可能性があり、そのような場合には、空間モデルを使用することが望ましいでしょう。
また、回帰モデルでは、変数 x の効果がすべての観測点で一定であることを仮定します。空間的に異質である場合、推定したい効果は空間によって異なる可能性があり、その違いを定量化することに興味があるかもしれません。このような場合、空間回帰モデルは効果の推定と解釈をより柔軟に行うことができます。
空間回帰アプローチの詳細については、このハンドブックの範囲外です。このセクションでは、最も一般的な空間回帰モデルとその用途の概要を説明し、この分野をさらに探求したい場合に役立つ参考文献を紹介します。
空間誤差モデル - これらのモデルは、空間ユニット間の誤差項が相関していると仮定していますが、その場合、データは標準的な最小二乗法(OLS) モデルの仮定に反することになります。空間誤差モデルは、「同時自己回帰（SAR）モデル」と呼ばれることもあります。空間誤差モデルは、 spatialreg パッケージの errorsarlm()（空間回帰関数、以前は spdep パッケージの一部だった）を使ってフィットさせることができます。
空間ラグモデル - このモデルは、ある地域 i の従属変数が、i の独立変数の値だけでなく、i に隣接する地域のそれらの変数の値にも影響されると仮定しています。空間誤差モデルと同様に、空間ラグモデルもまた、同時自己回帰（SAR）モデルと呼ばれることがあります。 空間ラグモデルは、spatialreg パッケージの lagsarlm() を使ってフィットさせることができます。
spdep パッケージには、標準的な最小二乗法 (OLS) モデル、空間ラグモデル、空間誤差モデルを決定するための便利な診断テストがいくつか含まれています。これらのテストは、ラグランジュ乗数診断と呼ばれ、データの空間依存性のタイプを特定し、どのモデルが最も適切かを選択するために使用することができます。lm.LMtests() を用いると、すべてのラグランジュ乗数検定を計算することができます。Anselin (1988) は、ラグランジュ乗数検定の結果に基づいて、どの空間回帰モデルを使用するかを決定するための便利なフローチャートツールも提供しています。
ベイズ階層モデル - ベイズを用いた手法は、いくつかの応用例において空間分析で一般的に使用されています。最も一般的なのは 疾患マッピング です。症例データの分布がまばらであったり（例えば、稀な結果の場合）、統計的に「ノイズが多い」場合には、潜在的な空間プロセスを考慮して疾病リスクの「平滑化された」推定値を生成するために使用することができるので、ベイズを用いた手法の方が好まれます。これにより、推定値の質が向上する可能性があります。また、データに存在する可能性のある複雑な空間的相関パターンを（選択することで）事前に指定することができ、独立変数と従属変数の両方における空間依存および非依存の変動を説明することができます。Rでは、CARbayes パッケージ（ ビニエット（vignette） 参照）や R-INLA（website や textbook 参照）を用いてベイズ階層モデルをフィットさせることができます。また、R は JAGS や WinBUGS のようなベイズ推定を行う外部ソフトウェアを呼び出すのにも使うことができます。
28.12 参考資料
R Simple Features and sf package vignette
R tmap package vignette
Spatial Data in R (EarthLab course)
Applied Spatial Data Analysis in R textbook
SpatialEpiApp - a Shiny app that is downloadable as an R package, allowing you to provide your own data and conduct mapping, cluster analysis, and spatial statistics.
An Introduction to Spatial Econometrics in R workshop