コロナウイルスのデータで遊ぶ
久しぶりの更新。 最近流行のコロナウイルスのデータで少し遊んでみる。 なお、コロナウイルス関連の情報は錯綜しており、本記事も必ずしも正しいことを言っているわけではないことには留意されたい。 また、分析ミス等があるかもしれない。
Rのパッケージ
こんなパッケージを発見した。
コロナウイルスの累積感染者数のデータを簡単に取得できる。今回はこれで少し遊んでみた。
データの取得
上記のパッケージを利用すると簡単に国別の累積感染者数のデータをリアルタイムで取得できる。 コードは下記の通り。多分。
# 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日)で感染者数が特に多い国、日本、中国、イタリア、韓国、イラン、フランスの累積感染者数を可視化してみる。 なお、日本のデータにはクルーズ船乗客も含まれているようだ。ただし、以下では、このクルーズ船の影響は無視して議論する。
これをみると、イタリア、韓国、イラン、フランスで指数関数的に累積感染者数が増加していることがわかる。中国は収束に向かいつつあるようだ。 日本は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)
次に、その国の人口で累積感染者数を割ってみた。すると、上記の図は以下のようになった。
これをみると日本がその人口のわりに累積感染者数が少ないことがわかる。一方で、イタリアや韓国はその割合が非常に高い。 これらの国に比べると、日本の感染対策はうまくいっているように思える。 コードは下記の通り。
# 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)
分析
最後に、日本と中国のデータを使って、簡単に分析をしてみる。 今回のような累積感染者数はどのような関数にしたがっているのだろうか。 よく用いられるのはロジスティック関数らしい。 それは以下の通り。
今回の文脈では、が感染者数、が最大感染者数、が増加率、が時間、が初期感染者数である。多分。 この関数のパラメータを日本と中国のデータに当てはめてみる。 推定すべきパラメータは と である。 はデータから経験的に分かるので推定する必要がない。今回はそれを2とした。
使用するRの関数はnls()である。これはformulaを自分で好きに決めることができ、そのformulaのパラメータを(おそらく)最小二乗法で推定してくれる。
これを用いてパラメータを推定し、その関数と実際のデータを図示してみた。 まずは、中国。
そして日本。
中国は、予想されるよりも遅く感染者が増加している。あえて解釈すると、当初、中国当局は感染者が増加していることを隠したかったのかもしれない。 また、予想されるより早く感染者が減少している。中国当局はもしかすると早めにコロナウイルス収束宣言を出したいのかもしれない。
また、日本は予想されるよりも収束していない。既に言及したように、最近検査体制の変化があったから感染者数が予想されるよりも増加したのかもしれない。
色んな解釈が成り立ちそうである。もっとも、累積感染者数がロジスティック関数に従うという前提のもとでだが。
コードは下記の通り。 なお、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))