R：疫苗接種後死亡人數 (1)

🌞 快要開學了，讓我們來複習一下統計學 …

高端疫苗接種後一天有4人死亡，所以有人說要暫停施打，我們用統計學的方法來分析看看好了。

⚠ 請注意：這是一份商業數據分析的教材，以下所引用的數字並不精確，也沒有對死亡率做出任何具體結論，我們只是希望利用時事提高學生們的興趣，用來說明數據分析與詮釋技巧。

1. 敘述統計

🌻 先做簡單的敘述統計

p = 17.8/2300/365    # 用109年的死亡人數估計每人每天死亡機率
n = 280000           # 接種人數
eVal = n * p         # 期望值
c(p, eVal)

[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人的機率不大於：

pbinom(4,n,p) 

[1] 0.2936

p > 0.05，所以不能拒絕「打疫苗會增加死亡機率」的假設

2.3 使用臨界值做檢定

依據假設檢定的邏輯

• 假設：如果打疫苗不會增加死亡機率，我們可以推論 …
• 檢定臨界值：在95%信心水準之下，要拒絕這個假設的臨界值是：

qbinom(0.95,n,p)

[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())

data.frame(alpha=q, criValue=qbinom(q,n,p))

  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)