R：疫苗接種後死亡人數 (1)
卓雍然, 中山大學 管理學術研究中心
2021-08-30 13:37:51
🌞 快要開學了，讓我們來複習一下統計學 …
高端疫苗接種後一天有4人死亡，所以有人說要暫停施打，我們用統計學的方法來分析看看好了。
⚠ 請注意：這是一份商業數據分析的教材，以下所引用的數字並不精確，也沒有對死亡率做出任何具體結論，我們只是希望利用時事提高學生們的興趣，用來說明數據分析與詮釋技巧。
1. 敘述統計
🌻 先做簡單的敘述統計
[1] 0.0000212 5.9368672由於接種人數很多(n=280,000)，所以雖然每人每天死亡機率(p=0.0000212)不大， 背景期望值(5.937)還是大於實際死亡人數(4)，這樣看起來好像還好， 不過，可能有人會說這個說法很「不科學」。
2. 推論統計
那我們就用比較「嚴謹」的推論統計方法分析看看。
2.1 背景死亡人數分布
為了避免不必要的複雜度，我們假設這28萬人是從2300萬人中隨機抽出，而且他們的死亡率與族群平均值之間沒有太大的差異；依據這一個假設，在沒有打疫苗的狀況下，這28萬人每一天的死亡人數應該呈現二項分佈：Binom[n=280000, p=0.0000212]
par(cex=0.7, mar=c(3,3,2,1)) # ploting options
k = 1:12
dbinom(k, n, p) %>%
barplot(names=k, border="white", main="PDF, Binom[n=280000, p=0.0000212]")
2.2 使用p值做檢定
依據二項分佈的累計密度函數
- 假設：如果打疫苗會增加死亡機率，我們可以推論 …
p值：看見一天死亡人數少於(含)4人的機率不大於：
[1] 0.2936p > 0.05，所以不能拒絕「打疫苗會增加死亡機率」的假設
2.3 使用臨界值做檢定
依據假設檢定的邏輯
- 假設：如果打疫苗不會增加死亡機率，我們可以推論 …
- 檢定臨界值：在95%信心水準之下，要拒絕這個假設的臨界值是：
[1] 10實際值(4)小於臨界值(10)，所以也不能拒絕「打疫苗不會增加死亡機率」的假設
🌷 所以，根據資料分析，我們的結論是 …
- 不能拒絕「打疫苗會增加死亡機率」
- 也不能拒絕「打疫苗不會增加死亡機率」
🙄 WTF … 你老闆會說： 這什麼鬼啊？ 這種結論根本沒有用啊！
3. 有圖有真相
📊 人家都說有圖有真相 …
par(cex=0.75, mar=c(3,3,2,1)) # ploting options
pbinom(1:12,n,p) %>% barplot( # plot cdf
names=1:12, ylim=c(0,1.15), border='white', family="BL", cex.axis=0.7,
main="二項分佈累計密度函數, n = 接種人數, p = 每天死亡機率")
abline(h=c(0.05,0.95), col='magenta',lty=2)
abline(h=seq(0,1,0.2),col='lightgray')
🙄 Hmm… 那請問你從上面這張圖都看到了什麼？
p=0.05和0.95那兩條水平的補助線是做什麼的呢？
🗿： 如果你是CDC或相關單位的幕僚，
- 你認為像上面這樣、看起來沒有結論的分析結果，對決策會有幫助嗎？
- 你要如何報告沒有結論的分析結果，老闆才會覺得你(的分析)是有用的呢？
- 換句話說，你要如何報告以上的分析結果，才不會被老闆打槍呢？
4. 不顯著的檢定
當檢定的結果不顯著，我們做報告的時候，最好把各種假設、各種信心水準的臨界值值都交代清楚，這樣看起來會比較專業，好比說：在這28萬人是由2,300萬人之中隨機抽出，並且族群之中死亡機率沒有太大變異的前提之下 …
🌻 28萬人打完疫苗之後1天之內的死亡人數，
- 要小於2個人 或
- 大於10個人
我們才能以95%的信心水準，分別推論打疫苗不會(或者會)增加死亡率。
這個講法會比“不能拒絕會增加死亡機率，也不能拒絕不會增加死亡機率” 容易懂，看起來也會比較專業。
🌷 我們還需要注意：
- 假設檢定是一種帶有不確定性的推論方法
- 以「95%的信心水準」推論打疫苗會(或者不會)「顯著」增加死亡率，是很容易造成誤會的說法，它不代表打疫苗會(或者不會)增加死亡率的機會高於95%，也不能估計死亡率增加或減少的幅度
- 一般人，即使是學過機率統計的同學，都很容易誤解假設檢定的結果，所以在社群網路之中傳播某些事件、或者某一天的統計數字，是很不恰當的做法
- 政府權責機構或是值得信賴的領域專家，才是比較可靠的消息來源。
5. 錯誤推論的機率與改進的方法
95%的信心水準，代表有5%的機率會推論錯誤，如果我們連續兩個禮拜，每天都用當天的資料做這個檢定，就會有大於50%的機會會做出至少一次錯誤的推論。
k = 1:20
par(cex=0.65, mar=c(3,3,2,1))
pgeom(k, 0.05) %>%
barplot(names=k,main="CDF,Geom[p=0.05]",border="white")
abline(h=0.5, col="magenta")
所以如果我們有很多天的資料，比較合理的做法是用所有的資料來做檢定，這樣做可以提高信心水準，降低錯誤判斷的機會。
🌻 當接種人數增加到60萬，觀察期間擴大到3天時 …
n = 600000
p = 3*17.8/2300/365
k = 20:60
q = c(0.005,0.01,0.05,0.95,0.99,0.995) # level of conf.
data.frame(k, prob=dbinom(k,n,p)) %>%
ggplot(aes(k, prob)) + geom_col(fill="gray",alpha=0.75) +
geom_vline(xintercept=qbinom(q,n,p),col='red',linetype='dashed') +
theme_bw() + ggtitle("PDF, dbinom[n=600K, p=6.4e-5]") +
theme(axis.title.x=element_blank())
alpha criValue
1 0.005 23
2 0.010 25
3 0.050 28
4 0.950 49
5 0.990 53
6 0.995 55🌻 在相同的前提下，60萬人打完三天之內，觀察到的死亡人數：
- 大於49/53/55人時
- 小於28/25/23人時
分別在95/99/99.5%的信心水準之下，我們可以推論打疫苗會(不會)增加死亡機率。
這個說法會比顯著、不顯著，或者能不能拒絕容易懂，對決策者來說，也比較有用對嗎？
PS：有同學說，這算小機率事件啊，要不用Poison做看看；其實dbinom[n.p]和dpois[n*p]是一樣的，不信我畫給你看 …
p = 17.8/2300/365; n = 280000; k = 1:12
data.frame(
k = k, binom = dbinom(k, n, p), ppois = dpois(k, n*p)
) %>% gather('distribution','density',-k) %>%
ggplot(aes(k, density, fill=distribution)) +
geom_col(position="dodge",alpha=0.5) +
scale_x_continuous(breaks=k)