dplyrメモ
・関数
mutate_cond <- function(.data, condition, ..., envir = parent.frame()) {
condition <- eval(substitute(condition), .data, envir)
.data[condition, ] <- .data[condition, ] %>% mutate(...)
.data
}
・例
mutate_cond(行名 == "*", 行名 = 0)
資料_2022/05/27
このサイトから「都道府県ごとの感染状況」のオープンデータを入手できる。
https://www3.nhk.or.jp/n-data/opendata/coronavirus/nhk_news_covid19_prefectures_daily_data.csv
これをもとに、現在も進行中と思われる新型コロナ流行第6波（2021/12/15から2022/05/26まで）の累積死者数を各都道府県ごとにプロットすると Fig. 1 のようになる。
また、感染者数と死亡者数のそれぞれについて、各都道府県の人口10万人あたりの値を散布図としてプロットすると、
一見して、大都市を抱える自治体で感染者、死者が多いのは明らかであるが、しかしたとえば東京都と大阪府では感染割合に比べて死亡率の差が大きい。
首都圏の埼玉県と千葉県の比較でも同じような傾向が認められる。
以下では両県のこの差異について、過去に遡ってもう少し仔細に検討してみる。
⬛千葉県と埼玉県における感染者数および死者数の推移
Fig. 3 は2020/03/01から2022/05/26までの感染者数と死者数の推移である（どちらも7日間移動平均）。
これを見ると、感染者数推移については、新型コロナ流行の初期から現在まで、両県の差はかなり小さい。
他方、死者数推移に関しては、第1波から第4波までは感染者数と同様に顕著な差は見られないが、第5波、第6波では千葉県の死者数が埼玉県のそれを上回る結果となっている。
第5波が始まった当時、千葉県の死者数が目立って増加したという指摘はSNS上でも散見された。
＞第5波だけで作成したら、千葉県が一番なはず。
— nagaya (@nagaya2013) September 4, 2021
過去30日の100万人あたりの死亡者累積です。たしかに千葉県が一番です。死亡は遅行指数ですので、あと10日くらいは見る必要がありますが、千葉は埼玉の3倍くらいです。知事の差でしょうか。こうしてみると広島は優秀ですね。https://t.co/iZDnJ8GB5R pic.twitter.com/LLGrHZ1yLL
Fig. 3 の死者数推移を見ると、この千葉県と埼玉県の差は、第6波になってさらに拡大しているように思われる。
そこで次に、もう少し定量的に両県のデータを比較検討してみよう。
⬛「波」の定義
コロナ流行の「波」の定義は多様であり一概に定められないが、厚労省データ（後述）の集計方式を参考に、まず死亡の波をそれぞれ以下のような期間とした。
- 第1〜4波： 〜2021/07/27
- 第5波： 2021/07/28〜2022/01/04
- 第6波： 2022/01/05〜2022/05/24
感染の波は死亡の波よりも先行するのが普通である。
その度合いは波ごとに、また地域によって異なると考えられるが、ここでは機械的に感染の波が１ヶ月先行するとした。すなわち、
- 第1〜4波： 〜2021/06/27
- 第5波： 2021/06/28〜2021/12/04
- 第6波： 2021/12/05〜2022/05/24
⬛千葉県と埼玉県における感染者数の比較
NHKのオープンデータから求めた各波の累積感染者数を下表に示す。単位は人口１０万人あたりの感染者数である。
| 第1〜4波 | 第5波 | 第6波 | |
| 千葉県 | 636 | 965 | 5485 |
| 埼玉県 | 629 | 950 | 5828 |
ちなみに、千葉県の人口は約630万人、埼玉県の人口は約740万人である。
両県の累積感染者と、この人口の比率の差は、第5波と第6波で有意である。ただし大小関係は逆転している。
⬛千葉県と埼玉県における年齢調整死亡率の比較
年齢調整死亡率は年齢構造の異なる集団や時点などの死亡率を比較するために用いられる。
計算式は、$$ \frac{\rm{（観察集団の年齢階級別死亡率×基準人口集団の年齢階級別人口）の総和}}{\rm{基準人口集団の総人口}} $$
ここで、観察集団の年齢階級別死亡率は、実際に報告された年齢階級別死亡数と実際の年齢階級別人口から算出する。
埼玉県の場合、年齢階級別死亡数のデータが、厚生労働省のウェブサイト：データからわかる－新型コロナウイルス感染症情報－より入手できる。（*1）
参照：「性別・年代別死亡者数（累積）」https://covid19.mhlw.go.jp/public/opendata/deaths_detail_cumulative_weekly.csv
このデータファイルには、全国およびいくつかの県を除く各都道府県の週ごとの集計が記載されている。
千葉県のデータは上記ファイルには記載されていないので、県が日々公表している「死亡者の発生状況」の情報をもとに、独自に作成したデータファイルを用いた。
参照：https://github.com/igproj-fusion/covid19_data/blob/main/covid_chiba_deaths.csv
各県の年齢階級別人口は、令和3年住民基本台帳年齢階級別人口（都道府県別）による。
参照：https://www.soumu.go.jp/main_content/000762442.xlsx
また、基準人口集団の年齢階級別人口とは、いわゆる“モデル人口”であり、従来は昭和60年モデル人口が用いられてきたが、ここでは平成27年モデル人口とした。
参照：https://www.mhlw.go.jp/content/12601000/000638712.pdf
以上の諸データにより算出した埼玉、千葉、および全国の年齢調整死亡率を下表に示す。なお、単位は人口10万人あたりの死亡数である。
| 第1〜4波 | 第5波 | 第6波 | |
| 埼玉県 | 11.23 | 3.51 | 7.83 |
| 千葉県 | 10.20 | 4.65 | 11.47 |
| 全 国 | 9.03 | 3.59 | 7.81 |
各「波」における全国の死亡率を1.00として、千葉県と埼玉県の推移を図示するとFig. 4のようになる。
コロナの流行が続くなか、このように死亡率の県間格差は拡大している。
全国との比較では、埼玉県では改善、千葉県は悪化傾向といえる。
感染者数にも差はあるが、死亡率この差を単に感染割合の違いだけで説明するのは難しい。いったい他のどんな要因がこうした差を生んでいるのだろうか。
⬛考察
オミクロン株の第6波では、全国で過去最大の死者数を記録しているが、国立感染研の最近のレポートには、
侵襲性の高い治療を希望されない場合や基礎疾患の悪化等の影響で重症の定義を満たさずに死亡する方など、新型コロナウイルス感染症が直接の死因でない事例も少なくないことが報告されており．．．
「新型コロナウイルス感染症の直近の感染状況等（2022年5月11日現在） 」
https://www.niid.go.jp/niid/ja/2019-ncov/11135-covid19-ab83th.html
とある。
たとえば、とりわけ死者数の多い大阪府の分析では、以下のように指摘している。
「感染した場合のリスクは、高齢者の方が圧倒的に高い。府の集計では、40～50代の死亡率は0.02％にとどまるが、60代以上では1.26％に跳ね上がる。コロナ自体は軽症でも、持病が悪化して亡くなるケースも相次ぐ。死者のうち70歳以上が93％を占める。」
「第６波」死者、大阪がなぜ全国で突出するのか…カギ握る高齢者：読売新聞オンライン
また、東京都の分析では、死亡した高齢者の多くが施設内感染によるとされている。
第6波の死者、施設での「看取り」相次ぐ 都内コロナ死536人分析：朝日新聞デジタル https://www.asahi.com/articles/ASQ336T90Q33UTIL03S.html
こうしたコロナ感染症の変化する実態への自治体の対処能力の違いが、死亡率の格差として各所で表面化したのではないか。
ちなみに、Fig. 2の東京都と神奈川県を見ると、感染者数は東京都よりも神奈川県は少ないものの、死亡率は同等ないし僅かに多い。これは感染者がほぼ同じで死亡率が異なる千葉県と埼玉県の関係とは対照的だが、格差を生み出す原因は同じではないかと考えられる。
先ごろ神奈川県の医療危機対策本部室・阿南英明統括官が表明した対策は、
「「早く」ということがキー。陽性者が出たら早く介入して治療につなげる。 治療までを３日以内に終わらせようと目標にしている」
次の感染流行に備えて 高齢者の感染拡大抑止策とは 神奈川
大きな損失を目の当たりにして、やっとこうした改善策を提示する専門家がいるというのは信じがたいことで、たとえば埼玉県のある医師は、第6波の始まる前に、その被害を予測して早期発見早期治療の重要さを強調していた。
開業医は自分の普段診てる患者さんがコロナ疑いでも検査もできないなら、患者さんの早期発見早期治療はできないから、悪化してから救急車→受け入れ病院なくて死亡、ですよ。夏の第5波より酷いことになるぞ。
— 職場換気担当köttur-lover22🐱 (@kottur_lover22) December 8, 2021
人命や健康の損失を極小化するには、感染を制御するのが第一義的であるのは言うまでもないが、いったん“防壁”が突破されて感染拡大が起こったときには、医療分野と介護・福祉分野の迅速な連携を図るのが肝要で、それには行政の適切な介入が欠かせない。
おそらくそうした行政能力の優劣が自治体間の死亡率格差に反映されている。
事後的にではあれ、神奈川県は改善策を示したが、さらに大きな損失を被った千葉県は、今後どのように対応していくのだろうか。
第5波、第6波で何が起こっていたのか。先ずはその精査が必要だが、埼玉県との死者数の差を指摘されて、公表基準が違うからなどと虚偽説明をするような知事では大変心もとない。
そのグラフの元は各都道府県の公表データで、千葉県は陽性の死亡者は原則公表、埼玉県は遺族の同意が前提で公表、基準が違うので当然差が出るということを申し上げました。
— 熊谷俊人（千葉県知事） (@kumagai_chiba) February 26, 2022
そして超過死亡数では千葉と埼玉に有意な差は無く、人口当たりの陽性者数と超過死亡には一定の相関性があるとも申し上げました https://t.co/XZaDapbkAP
*1:埼玉県が公表している死亡データについて、詳細が分かるのは全体の半分にも満たないらしい（2022/3時点）。なので、年齢調整に利用可能な情報は、厚労省の「性別・年代別死亡者数（累積）」に限られる。
参照：
パソコン教室・キュリオステーション志木店【公式】 on Twitter: "このデータから、埼玉県の死亡日別死者数のグラフを作りました。
総数の半分以下しか反映していない点にご注意ください。おおまかな波のありようは分かると思います。
https://t.co/oij0o35ZI4… "
千葉コロナマップ
千葉日報の 市町村別の感染者数（一覧表） | 千葉日報オンライン が毎日更新されており、そのデータをCSVファイルとして Gist に保存してある。
ファイル名：covid_chiba.csv
これの最終日のデータをマップとしてプロットする。
使用する千葉県のシェープファイルは geojson 形式で、筆者の Dropbox に置いてある。
URL：https://www.dropbox.com/s/4310z8hktjijzol/chiba_admin.geojson?dl=0
R スクリプト：
library(sf)
library(dplyr)
library(ggplot2)
library(ggthemes)
chiba <- st_read("https://www.dropbox.com/s/4310z8hktjijzol/chiba_admin.geojson?dl=1")
df <- read.csv("df <- read.csv("https://gist.githubusercontent.com/fusion0202/fdc1bfb603603b9dd8b42930483bbdac/raw/58e650e9087f2d121d19f1088335bd6c82fb867f/covid_chiba.csv", check.names = FALSE)
Date <- tail(colnames(df), 1)
df %>% rename(cases = Date) %>%
select(sichoson, cases) %>%
mutate(cases = replace(cases, cases == 0, NA)) -> dtleft_join(chiba, dt, by = c("SIKUCHOSON" = "sichoson")) %>%
ggplot() +
geom_sf(aes(fill = cases), alpha = 0.8, colour = 'grey5', size = 0.1) +
scale_fill_gradient(low = "#fef9f9", high = "#cd0505", na.value = "white", name = "No. of Cases") +
geom_sf_text(aes(label = cases)) +
theme_map() +
theme(legend.position = c(0.80, 0.05)) +
labs(title = "COVID-19 Cases in Chiba",
subtitle = paste("Data as of ", Date),
caption = "Data Source: https://www.chibanippo.co.jp/news/national/681627")
プロット例：
２色コロプレスマップ
先ごろ実施された第２５回参議院選挙の千葉県選挙区における３位当選者と４位落選者の得票について、５９の各市区町村での優劣をコロプレスマップとして描画する。
選挙の投開票結果については千葉県選挙管理委員会の〈投開票速報ページ（第25回参議院議員通常選挙）〉を参照（将来、リンクが変更される可能性がある）。
千葉県行政区分のシェープファイルは、過去記事「千葉県市区町村の高齢化率」と同じ手順で生成した。
各候補者の得票状況はCSVファイルとしてGitHubにアップロードしてある。
変数名：
- sikuchoson 市区町村名
- ishii 石井準一氏（自由民主党、得票率１位、当選）
- nagahama 長浜博行氏（立憲民主党、得票率２位、当選）
- toyoda 豊田俊郎氏（自由民主党、得票率３位、当選）
- kadota（省略）
- asano 浅野史子氏（日本共産党、得票率４位、落選）
- hiratuka（省略）
- total 総投票数
以下に、３位当選の豊田氏と４位落選の浅野氏について、得票数で競り勝った候補者によって地図上で塗り分けるためのRスクリプトを示す。
library(dplyr)
library(ggplot2)
library(ggthemes)
jp <- sf::st_read("japan_ver81.shp", stringsAsFactors = FALSE)
jp %>% filter(KEN %in% "千葉県") -> chiba
chiba <- chiba[-60,]
cand <- read.csv("sangiin2019vote.csv", stringsAsFactors = FALSE)
cand$who <- cand$toyoda > cand$asano
chiba2 <- left_join(chiba, cand, by = c("SIKUCHOSON" = "sikuchoson"))
ggplot() +
geom_sf(data = chiba2, aes(fill = who),
alpha = 0.8, colour = 'grey20', size = 0.1) +
scale_fill_manual(labels = c("浅野", "豊田"),
values = c("#FF6060", "#6060FF")) +
theme_bw() +
theme(axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.border = element_blank()) +
theme(legend.position = c(0.95, 0.1),
legend.title = element_blank())
実行結果：
LuaTeX-ja に乗り換えようか
長いこと、platex＋dvipdfmx でソースからPDFを作っていたが、texlive をインストールしてから設定までの手順が煩雑だし、世の中の流れは luaTeX に変わりつつあるようなので、さきほど texlive 2019 をインストールしたのを契機に luaTeX-ja を試してみた。
なお、フロントエンドには KDE アプリの kile を使っているので、その設定の仕方を以下に示しておく。
kile のメニュー（Settings > Configure Kile ...）から設定ダイアログを出して、まず LuaLaTeX の設定。
texlive はインストールしただけで実行ファイルへのパスを通していないので、ここでそれを指定している（/usr/local/texlive/2019/bin/x86_64-linux/lualatex）。
つぎに、ViewPDF の設定というか確認。kile に搭載されている Document Viewer が選択されていればよい。
最後に、QuickBuild。当然、LuaLaTeX+ViewPDF を選ぶ。
TeX のソースのプリアンブル部は至って簡素で、普通の文書なら、
\documentclass{ltjsarticle}
また、PDF に埋め込む和文フォントにヒラギノを使うときは、
\usepackage[hiragino-pron]{luatexja-preset}
これだけでよい。
dvipdfmx で苦労したときのことを思うと隔世の感がある。
なお、ヒラギノ・フォントファイルは /usr/local/texlive/texmf-local/fonts の配下に適当なサブディレクトリを作って放り込み、mktexlsr を実行しておく。
LuaLaTeX コマンドの実行時に２つ Warning が出るが、特に問題はなさそうである。
platex＋dvipdfmx と比べると処理時間は長くなるが、ストレスになるほどではない。
当面は、これで様子を見てみようと思う。
過去７日間の地震マップ
この記事を参考に、ごく基本的な地図作成のスクリプト。
データソースはUSGSの地震データで、過去１時間から１日、７日、３０日、また、マグニチュードの下限についても、いくつか選択できる。
ここでは、Past 7 Days > All Earthquakes を用いる。
ダウンロードした csv ファイルから、日付のフォーマットを修正し、続いて太平洋中心の地図（map_data の world2）にプロットするために、通常の経度：-180°〜180° を 0°〜360° に変換する処理などに dplyr の関数を利用しているのがポイント。
library(lubridate)
library(dplyr)
library(ggplot2)
library(ggthemes)
url <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.csv"
download.file(url = url, destfile = "all_week.csv")
df <- read.csv("all_week.csv", stringsAsFactors = FALSE)
df$time <- ymd_hms(df$time)
df <- arrange(df, -row_number())
quakes <- df %>% filter(type == "earthquake") %>%
mutate(longitude = ifelse(longitude < 0, longitude + 360, longitude))quakes$mag <- factor(round(quakes$mag))
map <- map_data('world2') %>% filter(region != "Antarctica")
title <- paste("Earthquakes map from",
as.Date(quakes$time[1]), "to", as.Date(quakes$time[nrow(df)]))
p <- ggplot()
p <- p + geom_map(data = map, map = map,
aes(x = long, y = lat, group = group, map_id = region),
fill = "grey60", colour = "grey90", size = 0.1)
p <- p + geom_point(data = quakes,
aes(x = longitude, y = latitude, colour = mag))
p <- p + scale_color_brewer(palette = 8)
p <- p + theme_map() + coord_equal() + ggtitle(title)
p <- p + theme(legend.position = "bottom")
p <- p + guides(colour = guide_legend(nrow = 1))
print(p)
estatapi パッケージを使って…
政府統計 e-Stat では、統計データがエクセル形式などのファイルのほかにデータベースとして提供されているものもある。後者のデータにAPI機能でアクセスするときに、estatapi パッケージを使うとRで比較的簡単に必要とするデータを抽出・取得できる。
e-Stat の API 機能を利用するには、①ユーザー登録を行い、② appID を取得する必要がある。②に関して、多少細かい手順を以下に示す。
- https://www.e-stat.go.jp/mypage/login からログイン
- 「マイページ」へ移行
- 「▷API機能（アプリケーションID発行）」へ
- 「*名称」には適当な文字列を、「*URL」は http://localhost で可
- 「発行」ボタンをクリック
この結果、appid 欄に表示される文字列が「アプリケーションID；appID」である。
以下では、都道府県の潜在扶養率のコロプレス・マップを描いてみる。
ちなみに潜在扶養率とは、65歳以上の人口に対する25～64歳の人口の比率のことで、国連経済社会局が毎年、国際比較を公表している。最新の結果は2019年版報告の中で言及されている。
同様の指標を、都道府県毎に計算してみようというわけである。
e-Stat の次のページにアクセスすると対象データを閲覧できる。
Rで、まず初めにこのデータのメタデータを取得する。
library(estatapi)appID <- "****************************************"dataID <- "0003312321"info <- estat_getMetaInfo(appID, dataID)
オブジェクト：info の内容は、RStudio の Environment でチェックするのがよい。
今回は「総人口（cat03/"001"）」の「男女計（cat01/"000"）」を抽出して、（全国分を除外し）47都道府県の潜在扶養率（ratio）を算出する。
dt <- estat_getStatsData(appID, dataID, cdCat01 = "000", cdCat03 = "001")
pop <- matrix(dt$value, nrow = 48, ncol = 19, byrow = F)
pop <- pop[-1,]
ratio <- apply(pop[,7:14], 1, sum) / apply(pop[,15:19], 1, sum)
以下は、この値をもとに、コロプレス・マップの描画。
都道府県のシェープファイルはウェブから何種類か入手できるが、ここでは三重大・奥村晴彦先生が公開されている GeoJSON形式のものを利用する（参照：シェープファイルを読む）。
library(sf)
library(dplyr)
library(ggplot2)
library(ggthemes)
jp = st_read("https://oku.edu.mie-u.ac.jp/~okumura/stat/data/japan.geojson",
stringsAsFactors = FALSE)jp$code <- as.integer(jp$code)
pref <- 1:47
dat <- data.frame(pref, ratio)
jp2 <- left_join(jp, dat, by = c("code" = "pref"))
ggplot() +
geom_sf(data = jp2,
aes(fill = cut_interval(ratio, length = 0.2)),
colour = 'grey20', size = 0.1) +
xlim(120, 150) +
scale_fill_brewer(palette = "YlOrRd",
name = "Potential\nSupport Ratio") +
theme_bw() +
theme(legend.position = c(0.85, 0.22))
