Follow @data_no_memo

メモ

個人的なメモです。他者にわかりやすく書くよりも未来の自分にわかりやすく書いています。なお、記事内容の正確さは保証できません。勉強中の身ですので、間違い等ご指摘頂けたら幸いです。

コロナウイルスのデータで遊ぶ

久しぶりの更新。 最近流行のコロナウイルスのデータで少し遊んでみる。 なお、コロナウイルス関連の情報は錯綜しており、本記事も必ずしも正しいことを言っているわけではないことには留意されたい。 また、分析ミス等があるかもしれない。

Rのパッケージ

こんなパッケージを発見した。

www.medrxiv.org

コロナウイルスの累積感染者数のデータを簡単に取得できる。今回はこれで少し遊んでみた。

データの取得

上記のパッケージを利用すると簡単に国別の累積感染者数のデータをリアルタイムで取得できる。 コードは下記の通り。多分。

# package of corona
remotes::install_github("GuangchuangYu/nCov2019")
library(nCov2019)

# probably list 3 is world data
last_update <- nCov2019::load_nCov2019(lang='en')
data <- nCov2019::load_nCov2019(lang='en')[[3]]

データの可視化

まずはこの記事の執筆段階(2020年3月8日)で感染者数が特に多い国、日本、中国、イタリア、韓国、イラン、フランスの累積感染者数を可視化してみる。 なお、日本のデータにはクルーズ船乗客も含まれているようだ。ただし、以下では、このクルーズ船の影響は無視して議論する。

f:id:abcxyzonetwothree:20200308191849p:plain

これをみると、イタリア、韓国、イラン、フランスで指数関数的に累積感染者数が増加していることがわかる。中国は収束に向かいつつあるようだ。 日本は2月中旬から一気に増加し、2月の後半でその増加率は減少しているものの、最近になってまた指数関数的に増加しているような傾向にある。 いろんな要因が考えられるが、1つの解釈は検査体制の変化が挙げられる。3月6日からコロナ検査の保険適用が決まり、検査が受けやすくなっているようだ。

コードは下記の通り。

# for plot and tidy data
library(patchwork)
library(tidyverse)

# plot Japan data
g_japan <- data %>% 
  filter(country == 'Japan') %>% 
  ggplot(aes(x = time, y = cum_confirm)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number') + ggtitle('Japan') 

# plot Italy data
g_italy <- data %>% 
  filter(country == 'Italy') %>% 
  ggplot(aes(x = time, y = cum_confirm)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number') + ggtitle('Italy') 

# plot Korea data
g_korea <- data %>% 
  filter(country == 'South Korea') %>%
  ggplot(aes(x = time, y = cum_confirm)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number') + ggtitle('South Korea')

# plot China data
g_china <- data %>% 
  filter(country == 'China') %>% 
  ggplot(aes(x = time, y = cum_confirm)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number') + ggtitle('China')

# plot Iran data
g_iran<- data %>% 
  filter(country == 'Iran') %>% 
  ggplot(aes(x = time, y = cum_confirm)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number') + ggtitle('Iran')

# plot France data
g_france<- data %>% 
  filter(country == 'France') %>% 
  ggplot(aes(x = time, y = cum_confirm)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number') + ggtitle('France')

(g_japan | g_china | g_italy) / (g_korea | g_iran | g_france)

次に、その国の人口で累積感染者数を割ってみた。すると、上記の図は以下のようになった。

f:id:abcxyzonetwothree:20200308192644p:plain

これをみると日本がその人口のわりに累積感染者数が少ないことがわかる。一方で、イタリアや韓国はその割合が非常に高い。 これらの国に比べると、日本の感染対策はうまくいっているように思える。 コードは下記の通り。

# for plot and tidy data
library(patchwork)
library(tidyverse)

# divide by population
pop_japan <- 126800000
pop_italy <- 60480000
pop_korea <- 51470000
pop_china <- 1386000000
pop_iran <- 81160000
pop_france <- 66990000

# plot Japan data
g_japan <- data %>% 
  filter(country == 'Japan') %>% 
  mutate(cum_per_pop = cum_confirm/pop_japan) %>% 
  ggplot(aes(x = time, y = cum_per_pop)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number (Percentage)') + ggtitle('Japan') +
  ylim(0, 0.0001)

# plot Italy data
g_italy <- data %>% 
  filter(country == 'Italy') %>% 
  mutate(cum_per_pop = cum_confirm/pop_italy) %>% 
  ggplot(aes(x = time, y = cum_per_pop)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number (Percentage)') + ggtitle('Italy') +
  ylim(0, 0.0001)

# plot Korea data
g_korea <- data %>% 
  filter(country == 'South Korea') %>%
  mutate(cum_per_pop = cum_confirm/pop_korea) %>% 
  ggplot(aes(x = time, y = cum_per_pop)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number (Percentage)') + ggtitle('South Korea')  +
  ylim(0, 0.0001)

# plot China data
g_china <- data %>% 
  filter(country == 'China') %>% 
  mutate(cum_per_pop = cum_confirm/pop_china) %>% 
  ggplot(aes(x = time, y = cum_per_pop)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number (Percentage)') + ggtitle('China')  +
  ylim(0, 0.0001)

# plot Iran data
g_iran<- data %>% 
  filter(country == 'Iran') %>% 
  mutate(cum_per_pop = cum_confirm/pop_iran) %>% 
  ggplot(aes(x = time, y = cum_per_pop)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number (Percentage)') + ggtitle('Iran') +
  ylim(0, 0.0001)

# plot France data
g_france<- data %>% 
  filter(country == 'France') %>% 
  mutate(cum_per_pop = cum_confirm/pop_france) %>% 
  ggplot(aes(x = time, y = cum_per_pop)) +
  geom_bar(stat = "identity") +
  xlab('Time') + ylab('Cumulative Number (Percentage)') + ggtitle('France') +
  ylim(0, 0.0001)

(g_japan | g_china | g_italy) / (g_korea | g_iran | g_france)

分析

最後に、日本と中国のデータを使って、簡単に分析をしてみる。 今回のような累積感染者数はどのような関数にしたがっているのだろうか。 よく用いられるのはロジスティック関数らしい。 それは以下の通り。

N = \frac{K}{1 + (K / N_0 - 1)exp(-rt)}

今回の文脈では、 Nが感染者数、 Kが最大感染者数、 rが増加率、 tが時間、 N_0が初期感染者数である。多分。 この関数のパラメータを日本と中国のデータに当てはめてみる。 推定すべきパラメータは K r である。 N_0 はデータから経験的に分かるので推定する必要がない。今回はそれを2とした。

使用するRの関数はnls()である。これはformulaを自分で好きに決めることができ、そのformulaのパラメータを(おそらく)最小二乗法で推定してくれる。

これを用いてパラメータを推定し、その関数と実際のデータを図示してみた。 まずは、中国。

f:id:abcxyzonetwothree:20200308194426p:plain

そして日本。

f:id:abcxyzonetwothree:20200308194448p:plain

中国は、予想されるよりも遅く感染者が増加している。あえて解釈すると、当初、中国当局は感染者が増加していることを隠したかったのかもしれない。 また、予想されるより早く感染者が減少している。中国当局はもしかすると早めにコロナウイルス収束宣言を出したいのかもしれない。

また、日本は予想されるよりも収束していない。既に言及したように、最近検査体制の変化があったから感染者数が予想されるよりも増加したのかもしれない。

色んな解釈が成り立ちそうである。もっとも、累積感染者数がロジスティック関数に従うという前提のもとでだが。

コードは下記の通り。 なお、nls()では推定すべきパラメータの初期値をうまく与えてやらないと推定値を求めてくれないので、試行錯誤で決める必要がある。 ただし、いろいろ初期値をいじってみたが、推定されるパラメータはほとんどその初期値に依存せず、初期値に対してロバストであった。

# china
temp_data_china <- data %>% 
  filter(country == 'China')
data_china <- data %>% 
  filter(country == 'China') %>% 
  mutate(past_time = 1:nrow(temp_data_china))

reg_china <- nls(cum_confirm ~ K / (1 + (K / 20 - 1) * exp( - r*past_time)), 
           start = c(K = 100000,r = 2), 
           data = data_china)

plot(data_china$past_time, data_china$cum_confirm, 
     xlab = 'Past Time', ylab = 'Cumulative Number',
     xlim = c(0, 100), ylim = c(0, 90000), title('China')) 
lines(data_china$past_time, fitted(reg_china))

# japan
temp_data_japan <- data %>% 
  filter(country == 'Japan')
data_japan <- data %>% 
  filter(country == 'Japan') %>% 
  mutate(past_time = 1:nrow(temp_data_japan))

reg_japan <- nls(cum_confirm ~ K / (1 + (K / 2 - 1) * exp( - r*past_time)), 
           start = c(K = 10000,r = 2), 
           data = data_japan)

plot(data_japan$past_time, data_japan$cum_confirm, 
     xlab = 'Past Time', ylab = 'Cumulative Number',
     xlim = c(0, 50), ylim = c(0, 1200), title('Japan')) 
lines(data_japan$past_time, fitted(reg_japan))