新型コロナウイルス SARS-CoV-2 とそれによる感染症 COVID-19 については、生物学的にはまだわかってないことが多い。さらに、感染力もよく分かってないし、流行のメカニズムもよく分かってない。しかし感染者と死者は世界中で毎日、数えられている。そのカウントデータに数理モデルを当てはめること自体は容易である。そしてそれで、収束の予測ができる。予測は予測でしかないが、簡単なので以下に解析例を示す。
データ
日本国内のデータを外国からもらう
データは欧州疾病予防対策センター でまとめているのを使う (使いやすいし、自分でまとめるより楽だから)。ここから CSV 形式のファイルをもらってきて、GNU R で読み込む。日付は文字列から日付に変更する。tidyverse はいつも使うから読み込む。また plot() で日本語を使いたいのでフォントを指定しておく。パソコンは mac を使う (楽だから)。NA を取り除くのは、まぁしてもしなくてもいい。気分の問題。
library(tidyverse);
library(ggfortify);
par(family = "HiraKakuProN-W3");
theme_set(theme_bw(base_family = "HiraKakuProN-W3"));
data <- read_csv("COVID-19-geographic-disbtribution-worldwide-2020-03-26.csv");
data <- mutate(data, dateRep = as.Date(dateRep, "%d-%m-%Y"));
data <- data[!is.na(data$dateRep), ];
感染者数の推移の様子を見てみる
読み込んだら、まずはプロット。日本のところだけ取り出して、 ggplot2 で一日あたりの発生感染者数のデータそのままと、その対数値を見る。これにはクルーズ船は入ってない。2019/12/31 から 2020/03/26 までの 87 日のデータである。
jpn <- data[data$countriesAndTerritories == "Japan", ];
ggplot(jpn, aes(x = dateRep, y = cases)) + geom_point();
ggplot(jpn, aes(x = dateRep, y = log(cases))) + geom_point();
感染者数は指数的／加速度的に増えている。ただ上下のバラつきは大きい。対数値は直線っぽくも見えるが、曲線と言われればそうかもしれない。いずれにせよ、直線っぽいこととバラつきの幅の狭さからして、対数の方が扱いやすそうに見える。
しかし対数のモデルはデータに０があると困る (log(0) は負の無限大に発散するから)。そんな時に役に立つのが一般化線型モデル (GLM) である。これなら０を含む非負のカウントデータをうまく扱える、ってこの本にも書いてある (突然の宣伝。いい本だから)。
一般化線形モデル
と言うわけで GLM でフィッティングしてプロットして、モデリングの良し悪しを見てみる。
GJM <- glm(cases ~ dateRep, data = jpn, family = poisson);
autoplot(GJM);
GLM で得られるモデルの診断プロットを autoplot で見ると、あんまりよくない。残差の分散は不均一だし、ＱＱプロットも両端のハズレ具合が大きい。でもまぁせっかくモデルができたので、一応、データと重ねて見てみる。
DAT <- data.frame(jpn, predict(GJM, se.fit = TRUE));
SMT <- mutate(DAT, cases = exp(fit), lwr = exp(fit - 1.96*se.fit), upr = exp(fit + 1.96*se.fit));
ggplot(DAT, aes(x = dateRep)) +
geom_point(aes(y = cases)) +
geom_line(aes(y = exp(fit))) +
geom_smooth(data = SMT, aes(y = cases, ymin = lwr, ymax = upr), stat = 'identity');
モデルはなんとなくうまくデータに乗っているように見えるが、これはいいモデルではないことは診断プロットが示している。仮にいい性質のモデルだったとしても、これは単調増加のモデルなので、今後についてはどこまでも感染者が増えていくという予測しかできない。もともとダメなモデルだったのだ。
(参考書はいい本なんですけどね...)。
対象に合ったモデルを使う
ベルカーブを対数変換する
そこで、データはもともと何だったのか、どんな推移が想定されるのかを思い出す。これは感染者数の推移なので、いわゆる流行曲線はベルカーブが想定されることが多い。指数的に増加し、頭を打って、指数的に減少する。ベルカーブといえば正規分布の曲線が代表的である。正規分布の曲線は A * exp(-(x-m)^2/s) という形式をしていて、対数を取ると、ただの二次多項式になる。二次式ならフィッティングのしやすさはお墨付きなので、喜び勇んで当てはめてみる。GNU R では線形回帰でやる方法もある (し、それも楽だ) が、ここでは非線形回帰の関数を使ってやってみる (好きだから)。感染者数が０の日はモデルに入れられないので、そこは無視するデータを作って回帰する。
y <- log10(jpn$cases);
t <- length(y):1;
t <- t[is.finite(y)];
y <- y[is.finite(y)];
M <- nls(y ~ a*(t-b)^2+c, start = c(a = -10, b = y[length(y)]/2, c = 10))
Y <- predict(M);
summary(M);
summary(M) で表示されるのは、得られたモデルの係数などの情報である。こんな感じ。
Formula: y ~ a * (t - b)^2 + c
Parameters:
Estimate Std. Error t value Pr(>|t|)
a -2.285e-04 1.142e-04 -2.001 0.051070 .
b 1.150e+02 3.028e+01 3.797 0.000412 ***
c 1.918e+00 3.875e-01 4.950 9.57e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2819 on 48 degrees of freedom
Number of iterations to convergence: 3
Achieved convergence tolerance: 5.621e-06
収束は早い (繰り返し計算３回)。そして係数のｐ値はかなり小さく追い込めていて、悪くない感じである。データ値 (感染者数の対数) とモデル、元のデータ値とモデルの指数値をそれぞれプロットするとこんな感じになる。
plot(t, log10(jpn$cases[is.finite(log10(jpn$cases))]),
xlab = "2019/12/31 からの日数",
ylab = "log10(感染者数)");
lines(t, log10(jpn$cases[is.finite(log10(jpn$cases))]));
lines(t, Y, col = "gray", lwd = 7);
plot(t, jpn$cases[is.finite(log10(jpn$cases))],
xlab = "2019/12/31 からの日数",
ylab = "感染者数",
ylim = c(0, 80));
lines(t, jpn$cases[is.finite(log10(jpn$cases))]);
lines(t, 10^Y, col = "gray", lwd = 7);
対数値の方 (上のプロット、モデリングした方) では、データとモデルの残差は外れ値があるにしても均一と言えば均一、「均一ではない！」と言い切ることはできないと思う。悪くなさそうだ。データを元の空間に戻した方 (下のプロット) も、データに乗っていると言っていいだろう。
最大値とその日を予測してみる
ここからはちょっと度胸のいるところである。
モデルの係数値では、b = 115 であった。これは、モデルの値が最大になるのは 2019/12/31 からの日数が 115 日のときということである。流行曲線がベルカーブを描くなら、その頂点は 115 日目であり、そのあとは感染者数は減少することをモデルは示している。そのあたりをプロットしてみるとこんな感じである。
x <- 1:150;
F <- predict(M, list(t = x));
max(10^F);
plot(x, 10^F, type = "l",xlab = "2019/12/31 からの日数", ylab = "感染者数");
points(t, jpn$cases[is.finite(log10(jpn$cases))]);
max(10^F) で表示される値は 82.8 である。つまりこのモデルでは、一日の感染者数は 115 日目に 82.8 人になり、そこが最大である。まぁもしモデル通りになったとしても、上下のバラつきがかなり大きいから、倍くらいになっても驚きはしない。
モデルはまともか
プロットを見ると「モデルの値に比べて感染者数のバラつきが広いな...」そして「右に行くほどバラつきが大きいな...」という気がする。しかもデータのある範囲はおおよそモデル曲線の変曲点のあたりまでなので、正直、少なくともデータの上昇が鈍るのがはっきり見えるところまでデータがないと、係数値がいくらデータにぴったりハマると言っても、その値を真に受けるのは危険だと思う。一方で、モデルはデータのその辺りに変曲点がある、と捉えているとも言えるが、これは変曲点を当てはめるとしたらそこしかない、ということでもある。
なんにせよ、感染者の発生数は４月下旬がピークで、そこから減少に転じることがモデルから予測できる。
このモデルが当たる条件
このモデルによる未来の予測が当たるための条件は、以下のようになろう。
- 感染者の発生と減少のメカニズムが途中で変わらないこと。
- 感染者数の対数が間違いなく二次関数に従うこと。
- モデルがデータにちゃんと乗っていて、残差の分布が全体的に一定なこと。
- データの観測方法が一貫していること。
- これは観測方法が正当かどうか、科学的な意味を正しく反映しているかという意味ではなく、途中で感染者の判定方法や集計方法が変わっていないか、ということ。
まぁ率直に言って、確実に当たるモデルとは言えないと思う。
しかし今後の推移が読めず社会全体に暗い雰囲気が立ち込めていることを考えると、当たってほしいなぁ...と思う。
追記
データが一日分増えたので、3/27 分までのデータで計算し直してみたところ、
ピーク日数：122.6 日 (大型連休突入)
ピーク高さ：106.5 人/日
になった。3/27 時点では強い増加が続いている。
3/28 のデータを使ったら、二次曲線が下に凸になっちゃったね... 微生物培養で言えば対数増殖期、直線を当てはめるべき状況かな。