俄羅斯地區的每日 COVID-19 病例圖表在我看來令人懷疑——從統計的角度來看是這樣嗎？

May 21, 2020

下面是俄羅斯克拉斯諾達爾邊疆區從 4 月 29 日到 5 月 19日新發現的 COVID 感染的每日圖表。該地區的人口為 550 萬人。

我讀到它並想知道 - 從統計的角度來看，這（新病例的相對平穩的動態）看起來好嗎？或者這看起來很可疑？疫情期間曲線可以這麼平坦，而不需要該地區當局對數據進行任何修改嗎？例如，在我的家鄉斯維爾德洛夫斯克州，圖表要混亂得多。

我是統計學的業餘愛好者，所以也許我錯了，這張圖表並沒有什麼特別之處。

根據 2020 年 5 月 18 日的新聞報導，自流行期開始到當天，該地區共進行了 136695 次 COVID-19 檢測。

截至 2020 年 5 月 21 日，該地區共記錄了 2974 例感染病例。

PS 這是我發現的一個鏈接，該鏈接指向具有更好看的統計數據的頁面，並且涵蓋了更長的時間，特別是針對克拉斯諾達爾邊疆區。在該頁面上，您可以將光標懸停在圖表上以獲取當天的具體數字。（標題使用“每日引出”病例數，條形標題“每日確診”病例數）：

這絕對是與眾不同的。

原因是像這樣的計數往往具有泊松分佈。這意味著它們的固有方差等於計數。對於附近的計數 $ 100, $ 的差異 $ 100 $ 意味著標準差幾乎 $ 10. $ 除非結果存在極端的序列相關性（這在生物學或醫學上不合理），否則這意味著大多數個體值應該隨機偏離基本假設的“真實”率，最多可達 $ 10 $ （上下），並且在相當多的情況下（大約三分之一）應該偏離更多。

這很難以真正穩健的方式進行測試，但一種方法是過度擬合數據，嘗試非常準確地描述它們，並查看殘差有多大。例如，這裡有兩個這樣的擬合，一個低平滑和一個過擬合 Poisson GLM：

此廣義線性模型 (GLM) 擬合（在 logit 尺度上）的殘差方差僅為 $ 0.07. $ 對於具有（視覺上）緊密擬合的其他模型，方差往往來自 $ 0.05 $ 到 $ 0.10. $ 這太小了。

你怎麼知道？ 引導它。 我選擇了一個參數引導程序，其中數據被從分佈中提取的獨立泊松值替換，其參數等於預測值。這是一個這樣的自舉數據集：

您可以看到單個值比以前波動了多少，以及波動了多少。

這樣做 $ 2000 $ 產生的時間 $ 2001 $ 方差（在兩到三秒的計算中）。這是他們的直方圖：

垂直的紅線標記數據的方差值。

（在一個擬合良好的模型中，這個直方圖的平均值應該接近於 $ 1. $ 平均值是 $ 0.75, $ 略少於 $ 1, $ 給出過擬合程度的指示。）

此測試的 p 值是這些測試的分數 $ 2001 $ 方差等於或小於觀察到的方差。由於每個自舉方差都較大，因此 p 值僅 $ 1/2001, $ 基本上為零。

我對其他模型重複了這個計算。在下面的R代碼中，模型根據樣條的節數k和度數而變化d。在每種情況下，p 值都保持在 $ 1/2001. $

這證實了數據的可疑外觀。事實上，如果你沒有說這些是案件的數量，我會猜到它們是某種東西的百分比。對於接近的百分比 $ 100 $ 這種變化將比這個泊松模型中的要小得多，而且數據看起來也不會那麼可疑。

這是產生第一個和第三個數字的代碼。（一個輕微的變體產生了第二個，在開始時替換X為X0。）
y <- c(63, 66, 66, 79, 82, 96, 97, 97, 99, 99, 98, 99, 98, 
      99, 95, 97, 99, 92, 95, 94, 93)
X <- data.frame(x=seq_along(y), y=y)

library(splines)
k <- 6
d <- 4
form <- y ~ bs(x, knots=k, degree=d)
fit <- glm(form, data=X, family="poisson")
X$y.hat <- predict(fit, type="response")

library(ggplot2)
ggplot(X, aes(x,y)) + 
 geom_point() + 
 geom_smooth(span=0.4) + 
 geom_line(aes(x, y.hat), size=1.25) + 
 xlab("Day") + ylab("Count") + 
 ggtitle("Data with Smooth (Blue) and GLM Fit (Black)",
         paste(k, "knots of degree", d))

stat <- function(fit) var(residuals(fit))
X0 <- X
set.seed(17)
sim <- replicate(2e3, {
 X0$y <- rpois(nrow(X0), X0$y.hat)
 stat(glm(form, data=X0, family="poisson"))
})

z <- stat(fit)
p <- mean(c(1, sim <= z))
hist(c(z, sim), breaks=25, col="#f0f0f0",
    xlab = "Residual Variance", 
    main=paste("Bootstrapped variances; p =", round(p, log10(length(sim)))))
abline(v = z, col='Red', lwd=2)