3-3_母子世帯の空間クラスター
Yuki Shibatsuji & Keisuke Takano
ステップ0：パッケージの読み込み
library(sf)
library(tidyverse)
library(ggplot2)
library(spdep)
library(magrittr)ステップ1：演習用データのダウンロードと確認
省略。
ステップ2：隣接関係ファイルの作成
tokyo_moran.shpを読み込んだ上で，隣接関係を生成する．
#tokyo_moranポリゴンを読み込み
tokyo_moran <- sf::st_read(dsn="tokyo_moran.shp",
#character型変数をfactor型に変換しない
stringsAsFactors=FALSE,
#読み込み結果を出力しない
quiet=TRUE,
#読み込むデータの文字コードはUTF-8と明示する
options="ENCODING=UTF-8") %>%
#読み込みの段階で「boshi%」が「boshi.」に置き換えられるので，変数名を付け直す
dplyr::rename(boshi_rate=boshi.)
#ポリゴンデータから隣接関係を生成
tokyo_moran_nb <- spdep::poly2nb(pl=tokyo_moran)
#隣接関係に関する集計情報
summary(tokyo_moran_nb)## Neighbour list object:
## Number of regions: 5206
## Number of nonzero links: 31238
## Percentage nonzero weights: 0.115259
## Average number of links: 6.000384
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 13 36 180 593 1173 1426 1012 520 147 55 18 14 2 6 2 3
## 19 20 21 25 27 37
## 1 1 1 1 1 1
## 13 least connected regions:
## 68 69 75 565 2313 3215 4084 4085 4151 4170 4171 4387 4450 with 1 link
## 1 most connected region:
## 647 with 37 links#隣接数のヒストグラム
hist(spdep::card(tokyo_moran_nb))
各種空間分析で用いる為に，隣接関係を隣接listへと変換．
#隣接関係を隣接listへと変換
tokyo_moran_nb_listw <- spdep::nb2listw(neighbours=tokyo_moran_nb)隣接関係に関するプロットを作成する．東京都全体でプロットすると潰れて見えなくなるので，国立市周辺だけに絞る．
#各小地域の重心ポイントを作成
tokyo_moran_cent <- sf::st_centroid(x=tokyo_moran)
#重心ポイント間で隣接関係を示すラインデータを作成
tokyo_moran_nb_line <- spdep::nb2lines(nb=tokyo_moran_nb,
coords=sf::st_geometry(tokyo_moran_cent),
as_sf=TRUE)
#国立市の小地域のみ取り出し
tokyo_moran_kunitachi <- tokyo_moran %>%
dplyr::filter(grepl(pattern="国立市",munic))
#国立市周辺で隣接関係をプロット
ggplot2::ggplot() +
#国立市の小地域ポリゴン
ggplot2::geom_sf(data=tokyo_moran_kunitachi) +
#重心ポイント
ggplot2::geom_sf(data=tokyo_moran_cent) +
#隣接関係ライン
ggplot2::geom_sf(data=tokyo_moran_nb_line) +
#描画範囲を国立市の小地域ポリゴンの範囲に合わせる
ggplot2::coord_sf(xlim=sf::st_bbox(tokyo_moran_kunitachi)[c("xmin","xmax")],
ylim=sf::st_bbox(tokyo_moran_kunitachi)[c("ymin","ymax")])
ステップ3：Global Moran’s Iの計算
Global Moran’s Iを計算・両側検定を行う．また，モラン散布図を作成する為に，母子世帯率boshi_rateの空間ラグ変数boshi_rate_lagを追加する．
#Global Moranを計算・両側検定を行う
spdep::moran.test(x=tokyo_moran$boshi_rate,
listw=tokyo_moran_nb_listw,
alternative="two.sided")##
## Moran I test under randomisation
##
## data: tokyo_moran$boshi_rate
## weights: tokyo_moran_nb_listw
##
## Moran I statistic standard deviate = 27.572, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2237720874 -0.0001921230 0.0000659816#母子世帯率boshi_rateの空間ラグ変数boshi_rate_lagを元データに追加
tokyo_moran <- tokyo_moran %>%
dplyr::mutate(boshi_rate_lag=spdep::lag.listw(x=tokyo_moran_nb_listw,
var=boshi_rate))x軸を母子世帯率boshi_rate，y軸を空間ラグ変数boshi_rate_lagとして，散布図をプロット（これがモラン散布図がやっていること）．
#x軸を母子世帯率boshi_rate，y軸を空間ラグ変数boshi_rate_lagとしたプロットを作る
ggplot2::ggplot(data=tokyo_moran,aes(x=boshi_rate,y=boshi_rate_lag)) +
#散布図であることを明示
ggplot2::geom_point() +
#プロットに当てはまる回帰直線
ggplot2::geom_smooth(method="lm",se=FALSE) +
#x軸を点線
ggplot2::geom_hline(yintercept=0,lty=2) +
#y軸を点線
ggplot2::geom_vline(xintercept=0,lty=2)
spdepのmoran.plot関数を使えば，実は楽に描ける（が微調整ができない）．
#モラン散布図を作成
spdep::moran.plot(x=tokyo_moran$boshi_rate,listw=tokyo_moran_nb_listw)
ステップ4：Local Moran’s Iの計算
Local Moran’s Iを計算・両側検定する．
#乱数のseedを設定
set.seed(1)
#Local Moran's Iを計算・両側検定
tokyo_moran_local <- spdep::localmoran_perm(x=tokyo_moran$boshi_rate,
nsim=999,
listw=tokyo_moran_nb_listw,
alternative="two.sided") %>%
#data.frame形式に変換
as.data.frame(stringsAsFactors=FALSE) %>%
#プロットに用いる変数のみ取り出し
dplyr::select(c("Ii","E.Ii","Var.Ii","Z.Ii","Pr(z != E(Ii)) Sim")) %>%
#計算結果の変数名を後々使いやすいように変更
magrittr::set_colnames(c("l_moran","l_moran_e","var","zval","pval"))
#元データに計算結果を結合
tokyo_moran <- cbind(tokyo_moran,tokyo_moran_local)ステップ5：Local Moran’s Iによるクラスターマップの描画
計算されたLocal Moran’s Iが統計的に有意な小地域を明示してプロット．
#有意性を表す変数signifを作成
tokyo_moran <- tokyo_moran %>%
dplyr::mutate(signif=dplyr::case_when(
pval<=0.001 ~ "p=0.001",
pval<=0.01 ~ "p=0.01",
pval<=0.05 ~ "p=0.05",
TRUE ~ "Not Significant"
))
#有意性をプロット
ggplot2::ggplot() +
#変数signifに基づき塗りつぶしを行うことを明示（その際ポリゴンの枠線を描画しない）
ggplot2::geom_sf(data=tokyo_moran,aes(fill=signif),lty=0) +
#手動で塗りつぶしの色分けを指定
ggplot2::scale_fill_manual(values=c("p=0.001"="green4",
"p=0.01"="green3",
"p=0.05"="green2",
"Not Significant"="grey"))
Local Moran’s Iの値が有意だった小地域について，モラン散布図の象限を表す変数を付与し，象限毎にプロット．
#母子世帯率boshi_rateの偏差boshi_rate_devと偏差の空間ラグ変数boshi_rate_dev_lagを追加
tokyo_moran <- tokyo_moran %>%
dplyr::mutate(boshi_rate_dev=boshi_rate-mean(boshi_rate)) %>%
dplyr::mutate(boshi_rate_dev_lag=spdep::lag.listw(x=tokyo_moran_nb_listw,
var=boshi_rate_dev)) %>%
#モラン散布図の象限を表す変数flagを追加
dplyr::mutate(flag=dplyr::case_when(
boshi_rate_dev>0&boshi_rate_dev_lag>0 ~ "HH",
boshi_rate_dev>0&boshi_rate_dev_lag<0 ~ "HL",
boshi_rate_dev<0&boshi_rate_dev_lag>0 ~ "LH",
boshi_rate_dev<0&boshi_rate_dev_lag<0 ~ "LL"
)) %>%
#有意でない小地域については，変数flagの値を「Not Significant」に変える
dplyr::mutate(flag=ifelse(signif=="Not Significant","Not Significant",flag))
#象限毎にプロット
ggplot2::ggplot() +
ggplot2::geom_sf(data=tokyo_moran,aes(fill=flag),lty=0) +
ggplot2::scale_fill_manual(values=c("HH"="red",
"HL"="pink",
"LH"="lightblue",
"LL"="blue",
"Not Significant"="grey"))
ステップ6：Local Moran’s Iの分析（クラスターの分布）
象限毎のLocal Moran’s Iのプロットを，市区町村ポリゴン・鉄道網ラインと重ねてプロット．
#市区町村ポリゴンを読み込み
tokyo_munic <- sf::st_read(dsn="tokyo_munic.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932")
#鉄道網ラインを読み込み
tokyo_railroad <- sf::st_read(dsn="tokyo_railroad.shp",
stringsAsFactors=FALSE,
quiet=TRUE,
options="ENCODING=CP932")
#市区町村ポリゴン・鉄道網ラインを重ねて象限毎にプロット
ggplot2::ggplot() +
ggplot2::geom_sf(data=tokyo_moran,aes(fill=flag),lty=0) +
ggplot2::scale_fill_manual(values=c("HH"="red",
"HL"="pink",
"LH"="lightblue",
"LL"="blue",
"Not Significant"="grey")) +
#市区町村ポリゴンを枠線のみプロット
ggplot2::geom_sf(data=tokyo_munic,fill="transparent") +
#鉄道網ラインを緑でプロット
ggplot2::geom_sf(data=tokyo_railroad,color="green")