Rでe-Statの統計表から都道府県別SMRを描く
このテキストでは、政府統計e-StatのAPIを使って、Rに直接統計表を読み込み、それを処理して特定の疾患（ICD-10に基づいた死因簡単分類で指定）の都道府県別の標準化死亡比(Standard Mortality Ratio, SMR)を算出し、その値によって都道府県を塗り分けた地図をかいてみます。すべてRの中だけで完結し、またプログラムの長さも短くて済んでいることに注目してください。
# appId はe-Statのサイトでユーザー登録をして取得したものを記載してください。
#appId <- "0000000000000000000000000000000000000000"
# 必要なパッケージを読み込みます。estatapi, tidyverse, NipponMap を使います
library(estatapi)## このサービスは、政府統計総合窓口(e-Stat)のAPI機能を使用していますが、サービスの内容は国によって保証されたものではありません。library(tidyverse)## -- Attaching packages ------------------------------------------------------ tidyverse 1.2.1 --## √ ggplot2 3.1.0 √ purrr 0.2.5
## √ tibble 1.4.2 √ dplyr 0.7.7
## √ tidyr 0.8.2 √ stringr 1.3.1
## √ readr 1.1.1 √ forcats 0.3.0## -- Conflicts --------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(NipponMap)標準化死亡比は、年齢調整死亡率（間接法）の算出過程で求められる値で、死亡数の観測数と期待数の比で表現されます。
\(\frac{ある地域で実際に観察された死亡数}{\sum{基準集団の年齢別死亡率\times 地域の年齢別人口}}\)
です。今回は基準集団としては日本全体を用います。
まずは都道府県別に年齢別人口が必要です。2015年国勢調査の値を使います。
# まず年齢階級別・都道府県別人口（期待死亡数の算出につかう）
# 平成２７年度国勢調査「出生の月(4区分)，年齢(5歳階級)，男女別人口(総数及び日本人) 」すなわちstatsDataId "0003149862"
d <- estat_getMetaInfo(appId, "0003149862")
# cdTab="020": "人口"
# cdCat01=全域
# cdCat02=最初の「総数」および末尾の「不詳」をのぞいたすべての年齢階級
# cdCat03=総数（男女別)
# cdCat04=総数（国籍）
# cdCat05=総数（出生の月）
# area=すべての都道府県（県コード01-47）
data_pop <- estat_getStatsData(appId, statsDataId = "0003149862", cdTab="020", cdCat01="00710", cdCat02=d$cat02$`@code`[c(-1,-23)], cdCat03="0000", cdCat04="0000", cdCat05="0000", cdArea=sprintf("%02d000",1:47))## Fetching record 1-987... (total: 987 records)str(data_pop)## Classes 'tbl_df', 'tbl' and 'data.frame': 987 obs. of 18 variables:
## $ tab_code : chr "020" "020" "020" "020" ...
## $ 表章項目 : chr "人口" "人口" "人口" "人口" ...
## $ cat01_code : chr "00710" "00710" "00710" "00710" ...
## $ 全域・人口集中地区（2015）: chr "全域" "全域" "全域" "全域" ...
## $ cat02_code : chr "1180" "1180" "1180" "1180" ...
## $ 年齢_2015 : chr "0～4歳" "0～4歳" "0～4歳" "0～4歳" ...
## $ cat03_code : chr "0000" "0000" "0000" "0000" ...
## $ 男女別_2015 : chr "総数（男女別）" "総数（男女別）" "総数（男女別）" "総数（男女別）" ...
## $ cat04_code : chr "0000" "0000" "0000" "0000" ...
## $ 国籍_2015 : chr "総数（国籍）" "総数（国籍）" "総数（国籍）" "総数（国籍）" ...
## $ cat05_code : chr "0000" "0000" "0000" "0000" ...
## $ 出生の月_2015 : chr "総数（出生の月）" "総数（出生の月）" "総数（出生の月）" "総数（出生の月）" ...
## $ area_code : chr "01000" "02000" "03000" "04000" ...
## $ 地域（2015） : chr "北海道" "青森県" "岩手県" "宮城県" ...
## $ time_code : chr "2015000000" "2015000000" "2015000000" "2015000000" ...
## $ 時間軸(年次) : chr "2015年" "2015年" "2015年" "2015年" ...
## $ unit : chr "人" "人" "人" "人" ...
## $ value : num 186010 42943 44460 88787 30148 ...次に、都道府県別の特定の疾患の死亡数を取得します。今回は疾患の例として、死因簡単分類の「結腸の悪性新生物」を用いてみます。このデータは人口動態統計のうち、報告書には掲載されていない「保管表」にありますが、2016年以降はこの表もe-statで取得できるようになっています。本来は人口と同じ年度の2015年度データを使いたいところですが、この年度の保管表は簡単に手に入らないので今回はちょっとずれますが2016年を使ってやってみます。
# そして死因別・都道府県別死亡数(観測死亡数実データってこと)
## 2015年人口動態データは不完全のようなので、2016年で。
# 2016年
# 人口動態統計 確定数 保管統計表 都道府県編（報告書非掲載表） 死因 すなわち statsDataId "0003215053"
d <- estat_getMetaInfo(appId, "0003215053")
# 02104_結腸の悪性新生物 についてやろう
# cdTab=死亡数
# cdCat01=総数（男女）
# cdCat02=02104_結腸の悪性新生物("02104")
# area=すべての都道府県（県コード01-47）
data_dd <- estat_getStatsData(appId, statsDataId = "0003215053", cdTab="10100", cdCat01="00100", cdCat02="02104", cdArea=sprintf("%02d000",1:47))## Fetching record 1-47... (total: 47 records)str(data_dd)## Classes 'tbl_df', 'tbl' and 'data.frame': 47 obs. of 12 variables:
## $ tab_code : chr "10100" "10100" "10100" "10100" ...
## $ 表章項目 : chr "死亡数" "死亡数" "死亡数" "死亡数" ...
## $ cat01_code : chr "00100" "00100" "00100" "00100" ...
## $ 性別 : chr "総数" "総数" "総数" "総数" ...
## $ cat02_code : chr "02104" "02104" "02104" "02104" ...
## $ 死因簡単分類 : chr "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" ...
## $ area_code : chr "01000" "02000" "03000" "04000" ...
## $ 都道府県・保健所: chr "北海道" "青森県" "岩手県" "宮城県" ...
## $ time_code : chr "2016000000" "2016000000" "2016000000" "2016000000" ...
## $ 時間軸(年次) : chr "2016年" "2016年" "2016年" "2016年" ...
## $ unit : chr "人" "人" "人" "人" ...
## $ value : num 1816 537 448 637 418 ...では、最後に基準集団すなわち日本全体での疾患別・年齢階級別死亡率を取得します。こちらは人口動態統計にそのものの表があるので使ってみます。ただ、この死亡率は有効数字があまり大きくないので誤差が大きいから、本当は死亡実数の表と人口から死亡率を算出したほうがいいですが、ここでは簡単のために人口動態統計に掲載されている死亡率をそのまま用いてみます。
# 最後に死因別・年齢階級別の全国での死亡率(期待死亡数に使う)
# 人口動態調査 / 人口動態統計 確定数 死亡 すなわち statsDataId "0003214724"
d <- estat_getMetaInfo(appId, "0003214724")
# cdTab=死亡率
# cdCat01=最初の「総数」「０歳」「１歳」「２歳」「３歳」「４歳」「１－４歳」および末尾の「（再掲）」たちを除いたすべて
# cdCat02=総数（男女）
# cdCat03=02104_結腸の悪性新生物("02104")
data_dr <- estat_getStatsData(appId, statsDataId = "0003214724", cdTab="10110", cdCat01=d$cat01$`@code`[c(-1:-7,-29:-32)], cdCat02="00100", cdCat03="02104")## Fetching record 1-21... (total: 21 records)str(data_dr)## Classes 'tbl_df', 'tbl' and 'data.frame': 21 obs. of 12 variables:
## $ tab_code : chr "10110" "10110" "10110" "10110" ...
## $ 表章項目 : chr "死亡率" "死亡率" "死亡率" "死亡率" ...
## $ cat01_code : chr "00170" "00180" "00200" "00220" ...
## $ 年齢(5歳階級): chr "0～4歳" "5～9歳" "10～14歳" "15～19歳" ...
## $ cat02_code : chr "00100" "00100" "00100" "00100" ...
## $ 性別 : chr "総数" "総数" "総数" "総数" ...
## $ cat03_code : chr "02104" "02104" "02104" "02104" ...
## $ 死因簡単分類 : chr "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" ...
## $ time_code : chr "2016000000" "2016000000" "2016000000" "2016000000" ...
## $ 時間軸(年次) : chr "2016年" "2016年" "2016年" "2016年" ...
## $ unit : chr "人口10万対" "人口10万対" "人口10万対" "人口10万対" ...
## $ value : num NA NA NA 0.1 0.1 0.2 0.7 1.5 2.4 4.5 ...
## ..- attr(*, "problems")=Classes 'tbl_df', 'tbl' and 'data.frame': 3 obs. of 4 variables:
## .. ..$ row : int 1 2 3
## .. ..$ col : int NA NA NA
## .. ..$ expected: chr "a number" "a number" "a number"
## .. ..$ actual : chr "-" "-" "-"それでは、これで準備が整いました。SMRを都道府県別に算出してみましょう。
data_pr <- data_pop %>% # 都道府県別・年齢別人口表から始めます
# 基準集団の年齢別死亡率表を、年齢階級をキーにして結合します
left_join(data_dr, by=c("年齢_2015"="年齢(5歳階級)"), suffix=c(".pop", ".ref")) %>%
# 都道府県の年齢階級別人口と全国の年齢階級別死亡率の積をとって期待死亡数とします
mutate(exp_d = value.pop * value.ref / 100000.0) %>%
# 都道府県別にグループ化することを宣言
group_by( area_code ) %>%
# ０－５歳などで期待死亡数がゼロのところは省きます
filter(!is.na(exp_d)) %>%
# グループ（都道府県）ごとに期待死亡数の全年齢層の和をとります
summarise(exp_dd = sum(exp_d)) %>%
# 都道府県別の実死亡数の表を結合します
left_join(data_dd, by="area_code") %>%
# 実死亡数/期待死亡数を計算してSMRとします
mutate(SMR= value / exp_dd) %>%
# 都道府県コードにより並べ替えておきます
arrange(area_code)
str(data_pr)## Classes 'tbl_df', 'tbl' and 'data.frame': 47 obs. of 14 variables:
## $ area_code : chr "01000" "02000" "03000" "04000" ...
## $ exp_dd : num 1575 397 407 607 362 ...
## ..- attr(*, "problems")=Classes 'tbl_df', 'tbl' and 'data.frame': 3 obs. of 4 variables:
## .. ..$ row : int 1 2 3
## .. ..$ col : int NA NA NA
## .. ..$ expected: chr "a number" "a number" "a number"
## .. ..$ actual : chr "-" "-" "-"
## $ tab_code : chr "10100" "10100" "10100" "10100" ...
## $ 表章項目 : chr "死亡数" "死亡数" "死亡数" "死亡数" ...
## $ cat01_code : chr "00100" "00100" "00100" "00100" ...
## $ 性別 : chr "総数" "総数" "総数" "総数" ...
## $ cat02_code : chr "02104" "02104" "02104" "02104" ...
## $ 死因簡単分類 : chr "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" "02104_結腸の悪性新生物" ...
## $ 都道府県・保健所: chr "北海道" "青森県" "岩手県" "宮城県" ...
## $ time_code : chr "2016000000" "2016000000" "2016000000" "2016000000" ...
## $ 時間軸(年次) : chr "2016年" "2016年" "2016年" "2016年" ...
## $ unit : chr "人" "人" "人" "人" ...
## $ value : num 1816 537 448 637 418 ...
## $ SMR : num 1.15 1.35 1.1 1.05 1.16 ...
## ..- attr(*, "problems")=Classes 'tbl_df', 'tbl' and 'data.frame': 3 obs. of 4 variables:
## .. ..$ row : int 1 2 3
## .. ..$ col : int NA NA NA
## .. ..$ expected: chr "a number" "a number" "a number"
## .. ..$ actual : chr "-" "-" "-"これで都道府県別の結腸の悪性新生物のSMRが出ましたから、日本地図をSMRの値で塗り分けてみましょう。
JapanPrefMap(col=heat.colors(nrow(data_pr))[order(data_pr$SMR,decreasing=TRUE)])
このように、データの取得から加工、計算、視覚化まですべてRだけで完了させることができました。