Introduction to empirical Bayes estimation (with R code)
延續上一篇關於Bayesian Statistics的研究, 本篇Post將進一步闡述貝氏統計的應用。以下的內容與程式碼主要翻譯自此post: http://varianceexplained.org/r/empirical_bayes_baseball/ . 這篇文章主要為研讀之後的整理和其他發想。
下列兩個比率數字, 哪個比較大?
裝有10個球的盒子中有4個紅球和6個白球 è 紅球的比例?
裝有1000個球的袋子裡有300個紅球和700個白球 è 紅球的比例?
很明顯, 當然是
4/10 = 0.4 大於 300/1000=0.3
但是, 假設今天你是球隊老闆,
正在評估兩位潛力球員。你以下列兩位球員的打擊成績做為評估標準:
球員A上場10次,擊出4支安打
球員B上場1000次,共擊出300支安打
雖然球員A有較高的打擊率,
但是僅僅10次的打擊紀錄, 並無法提供足夠的可信度。一般職棒選手的打擊率大約為0.27, 球員A高達四成的打擊率,期中運氣的成分居多;
反倒球員B 1000次的打擊紀錄, 較能證明他是一個優於平均的打擊者。
Empirical Bayes estimation
這篇文章同樣將使用棒球比賽的例子, 來說明一個十分有效用以估計資料比率的統計技術,
來幫助我們分析類似下列的數據資料:
表一
Success
|
Total
|
11
|
104
|
82
|
1351
|
2
|
26
|
0
|
40
|
1203
|
7592
|
5
|
166
|
一般我們可能取得表一這種success(成功)/total(總計)成對型式的數據,
然後用以估計特定事件的成功比率。每筆資料可能代表著:
- 廣告點擊率:公司投放了許多不同的廣告, 你想知道哪一個有比較高的Clickthrough rates?
- 網站使用者類型: 你想知道到訪公司網站的使用者, 有多少會點擊閱讀一篇文章, 或是點擊某個商品後決定購買?
當處理表一這類數字時, 我們需特別注意分母total數字的大小,
對於分析判斷的影響。1/2和50/100不會表達出同樣的意義;
0/1和0/1000也是一樣! 處理這種issue的方法之一是於分析前先訂立一個最小值,然後將total數小於最小值的資料排除。但這種方法並無法應用至所有的狀況, 因為有些時候我們可能將重要的資訊在尚未分析前就先被我們篩選掉了。
這篇post將介紹empirical Bayes estimation來處理上述的問題。我們將透過Beta distribution來進行預測分析。Empirical Bayes
estimation的一個特點是, 如果分析的資料數據夠多, 我們可以不需要prior expectations. 我們使用Lahman package, Batting
這個 Dataset來進行empirical Bayes estimation分析, 預測的目標是提升對於每位球員的打擊率預測。在進行模型建立前, 先進行Data prepare和clean的工作
library(dplyr)
library(tidyr)
library(Lahman)
career <- Batting %>%
filter(AB > 0) %>%
anti_join(Pitching, by = "playerID") %>%
group_by(playerID) %>%
summarize(H = sum(H), AB = sum(AB)) %>%
mutate(average = H / AB)
# 將nameFirst 和nameLast 組合成names欄位後, 使用player IDs進行join
career <- Master %>%
tbl_df() %>%
select(playerID, nameFirst, nameLast) %>%
unite(name, nameFirst, nameLast, sep = " ") %>%
inner_join(career, by = "playerID") %>%
select(-playerID)
- 欄位AB表示[At Bats-打數], 基本上就簡單地當作是上場打擊的次數; 欄位H表示[Hits-安打數]
- 排除投手(filter(AB > 0))
- 加總球員生涯的安打數和打數, 計算出他們的打擊率
先看一下Data prepare和clean後的dataset:
首先,我們想知道誰是最厲害的打者?? 回答這個問題, 我們直覺會依照打擊率(average)欄位進行排序:
career %>%
arrange(desc(average)) %>%
head(5) %>%
kable()
這幾位看起來並不是我們真的要找的強打球員。他們僅僅是曾經上場打擊個一、兩次,然後成功的擊出安打幸運兒。那麼最差的球員情況?
career %>%
arrange((average)) %>%
head(5) %>%
kable()
資料看起來還是不對。很明顯[打擊率(average)]這個數字不是一個好的預估值。
我們需要建立一個可信的分析預測模型。在進入模型建立前, 我另外用data.table建立了同樣的Career資料, 並用data.table的指令,快速的瀏覽了一次資料內容如下:
Estimate a prior from all your data
首先, 看看所有球員的打擊率(average)分佈
career %>%
filter(AB >= 500) %>%
ggplot(aes(average)) +
geom_histogram(binwidth = .005)
 |
| 圖一 |
Empirical Bayes estimation的第一步是使用已有的資料去預估beta prior. 使用已有的資料去預估beta prior 不是典型(typical)的Bayesian方法 – 一般是在事前就應該先決定預估模型的prior. 雖然, 關於何時或是甚麼情況下適合使用Empirical Bayes方法仍有許多的辯論和探討, 但基本上問題的關鍵, 是在我們實際上已經握有多少資料。如果資料數量夠大, 我們其實可以依此做出一個良好的預估值, 而不需仰賴過於單一來源。在此我使用原文來闡述這樣個觀點:
Empirical Bayes is an approximation to more exact Bayesian methods- and with the amount of data we have, it’s a very good approximation.
就圖一來看, 此數據分配是一個不錯的beta prior選擇。(那甚麼是不好的? 例如直方圖(Histogram)有兩個以上的峰柱(Peak), 此時我們需要更複雜的模型。)。
因此就可以建立下列方程式:
X ~ Beta(α0,β0)
與上一篇Post情況不一樣的地方在於, 上篇Post我們握有的資訊僅為對於打擊率的Expected value以及此數值的偏差範圍; 在這一篇Post, 我們擁有一群數據, 我們要做的是找出Beta distribution α 和 β的數值可以match這一群數據的分佈。R提供許多方法可計算出α 和 β的數值, 在這裡, 我們使用MASS package中的fitdistr function.
為了提高Prior distribution的預測精準度, 原文作者僅抓取AB(打數) >=500的球員進行Beta distribution的模型建置。計算的結果α0=78.611 和 β0=224.875; 我們來看看此Beta distribution是否fit全體球員的數據分佈:
ggplot(career_filtered) +
geom_histogram(aes(average, y = ..density..), binwidth = .005) +
stat_function(fun = function(x) dbeta(x, alpha0, beta0), color = "red",size = 1) +
xlab("Batting average")
還算不錯,雖然不是100% fit, 但已經是一個可供我們進行後續分析的模型!
Using the Beta distribution as a prior to update each individual Estimate
我們將使用上述建立的Beta distribution, 來更新我們對於每位球員的打擊率預估值。更新的方法就如同上篇Post描述, 將α0加上球員的安打數, 以及將α0+β0加上球員的打數。以文章一開頭的例子來說:
球員A上場10次,擊出4支安打
球員B上場1000次,共擊出300支安打
球員A的打擊率預估值=
(4+α0) / (10+α0+β0) = (4+78.7)/ (10+78.7+224.9) = 0.264
球員B的打擊率預估值=
(300+α0) / (1000+α0+β0) = (300+78.7)/ (1000+78.7+224.9)
= 0.29
因此, 即使4/10>300/1000, 此時我們對於球員B的打擊率預估值仍將是會高於球員A的打擊率預估值。
我們可以使用下列的程式碼更新所有球員的打擊率預估值:
career_eb <- career %>%
mutate(eb_estimate = (H + alpha0) / (AB + alpha0 + beta0))
現在,我們可以來回答: 誰是最強的打者??
options(digits = 3)
career_eb %>%
arrange(desc(eb_estimate)) %>%
head(5) %>%
kable()
倒數5位的球員是??
options(digits = 3)
career_eb %>%
arrange(eb_estimate) %>%
head(5) %>%
kable()
由於empirical Bayes estimates將球員群體的長期數據做為分析的依據, 我們將不會再受到0/1或是1/1這類數據影響而破壞了我們的分析判斷。總體來說, 我們可以利用下圖看出empirical Bayes estimates如何改變球員的打擊率預估。
ggplot(career_eb, aes(average, eb_estimate, color = AB)) +
geom_hline(yintercept = alpha0 / (alpha0 + beta0), color = "red", lty = 2) +
geom_point() +
geom_abline(color = "red") +
scale_colour_gradient(trans = "log", breaks = 10 ^ (1:5)) +
xlab("Batting average") +
ylab("Empirical Bayes batting average")
中間的那條水平虛線, 座落在Y軸=0.259的數值(此數值=α0/(α0+β0)), 這個數字為在沒有任何球員的實際打擊數據下, 我們對於每一球員的打擊率預測值。原文對於這水平虛線上下的數據點描述如下:
Notice that
points above that line tend to move down towards it, while points below it move
up.
另外一條略微傾斜的垂直紅色實線, 為X=Y的數據位置。資料靠進這條紅色實線的數據,
為那些球員的平均打擊率(即使用球員各自的安打數/打數),
較小被empirical Bayes縮減(shrink)的數據點; 可以注意到這些數據點有個特性就是打數較高(數據顏色偏向淺亮的藍色), 表示我們會傾向去相信使用球員的平均打擊率做為預測數值。我們稱此方法為shrinkage: 我們將所有的預測數值朝平均數方向移動,
移動的大小取決於我們握有多少的個別實際數據; 如果某一個別實際數據很少(例如10打數4安打),
則移動的數值就比較大; 反之, 若實際數據很多(例如1000打數300支安打),
則移動的數值就很小。概括說來, shrinkage就是:
Extraordinary outliers require extraordinary
evidence.Implementation codes to the production system
總結empirical Bayes estimation包含了兩個步驟:
1. 預估你所有資料,產生分佈圖(the
overall distribution of your data)
2. 使用上述分佈(distribution)作為prior去預測每一個獨立目標的平均
第一個步驟可以”離線”執行 – 分析所有的資料然後產生分佈(distribution); 第二步則是在每一個獨立目標或事件發生時,
產生/更新該目標的平均值; 我們可以將此方法運用在如預估某一個廣告的成功率, 或分析單一使用者在其各種行為模式下, 哪些類型的行為模式會產生特定的選擇。
此外,由於我們使用的是beta和binomial
distribution, 第二步驟可以簡單的將此方法寫入我們的系統Code中; 將一個數字(α0)加上success,
將另外一個數字(α0+β0)加上total:
float estimate = (successes + α0) / (total + α0 + β0);
by J.D.
