第 82 章 共軛先驗概率 Conjugate priors

本章節我們重溫一下最早在貝葉斯統計學入門部分(Chapter 40)介紹過的一些基本原則。特別是關於共軛先驗概率的概念,並提供一些使用BUGS模型的例子來展示如何運算這些模型。

82.1 貝葉斯推斷的基礎

在一個臨牀試驗中,作爲一名貝葉斯統計學者,必須清晰明瞭地闡述如下幾個問題:

  1. 合理地描述目前爲止,在瞭解本次試驗數據的結果之前,類似研究曾經給出過的療效差異的報告結果,可能的取值範圍 (the prior distribution);
  2. 本次試驗數據得到的結果,支持怎樣的療效差異 (the likelihood);
    之後需要將上述兩個資源通過合理的數學模型結合在一起,用於產生
  3. 最終療效是否存在的意見和建議,證據的總結 (the posterior distribution)。

最後第三步將先驗概率和似然相結合的過程,用到的是貝葉斯定理(Bayes Theorem)。通過貝葉斯定理把目前位置的經驗作爲先驗概率統合現階段試驗數據給出的似然的過程,其實是一個自我學習不斷更新知識的過程。經過貝葉斯推斷,給出事後概率分布之時,我們可以拿它來做什麼呢?

  • 估計和評價治療效果,治療差異。
  • 估計和評價模型中參數的不確定性。
  • 計算你感興趣的那個變量(可以是療效差異,可以是模型中的參數)達到或者超過某個特定目標值的概率。
  • 預測你感興趣的那個變量可能存在的範圍。
  • 作爲未來要進行的試驗設計階段的先驗概率分布。
  • 給決策者提供證據。

你是否還記得貝葉斯定理的公式: 如果\(A, B\)分別標記兩個事件,那麼有

\[ p(A|B) = \frac{p(B|A)p(A)}{p(B)} \]

如果\(A_i\)是一系列互斥不相交事件,也就是\(p(\cup_iA_i) = \sum_ip(A_i) = 1\),貝葉斯定理可以被改寫成爲:

\[ p(A_i|B) = \frac{p(B|A_i)p(A_i)}{\sum_jp(B|A_j)p(A_j)} \]

貝葉斯統計學推斷從根本上的特點在於嚴格區分:

  1. 觀測數據 \(y\),也就是試驗獲得的數據。
  2. 未知參數 \(\theta\),這裏,\(\theta\)可以用統計學工具來描述,它可以是統計學參數,可以是缺失值,可以是測量誤差數據等等。
  • 這裏貝葉斯統計學把未知參數當做一個可以變化的隨機變量(parameters are treated as random variables);
  • 在貝葉斯統計學框架下,我們對參數的不確定性進行描述(we make probability statements about model parameters)。
  • 在概率論統計學框架下,統計學參數是未知,但是確實不變的。使用概率論統計學進行推斷時,我們只對數據進行不確定性的描述(parameters are fixed non-random quantities and the probability statements concern the data.)

貝葉斯統計學推斷中,我們仍然需要建立模型用來描述 \(p(y|\theta)\),這個也就是概率論統計學中常見的似然(likelihood)。似然是把各個變量關聯起來的完整的概率模型(full probability model)

從貝葉斯統計學推斷的觀點出發,

  • 在實施試驗,收集數據之前,參數(\(\theta\))是未知的,所以它需要由一個概率分布(probability distribution)來反應它的不確定性,也就是說,我們需要先對參數可能來自的分布進行描述,指定一個先驗概率(prior distribution)\(p(\theta)\)
  • 試驗進行完了,數據整理分析之時,我們知道了\(y\),這就是我們來和先驗概率結合的似然,使用貝葉斯定理,從而獲得給定了觀測數據之後(conditional on),服從先驗概率的參數現在服從的概率分布,這被叫做事後概率分布(posterior distribution)

\[ p(\theta|y) = \frac{p(\theta)p(y|\theta)}{\int p(\theta)p(y|\theta)d\theta} \propto p(\theta)p(y|\theta) \]

總結一下就是,

  1. 先驗概率(prior distribution),\(p(\theta)\)描述的是在收集數據之前參數的不確定性。
  2. 事後概率(posterior distribution),\(p(\theta | y)\) 描述的是在收集數據之後參數的不確定性。

82.2 二項分布(似然)數據的共軛先驗概率

沿用前一個章節新藥試驗的例子。我們在實施試驗之前對藥物的認知是,我們認爲它的藥效概率可能在 0.2-0.6 之間。我們把這個試驗前對藥物療效的估計認知翻譯成一個服從均值爲0.4,標準差爲0.1的Beta分布。這個Beta分布使用的參數(參數的參數被叫做超參數, hyperparameter),是9.2, 13.8,寫作\(\text{Beta}(9.2, 13.8)\)。那麼現在我們假設試驗結束,收集的20名患者中15名療效顯著。接下來貝葉斯統計學家要回答的問題是,這個試驗結果對先驗概率分布產生了多大的影響(How should this trial change our opinion about the positive response rate)?

在這個例子中,我們現在來詳細給出先驗概率和似然。

  • 似然 likelihood (distribution of the data):
    如果患者可以被認爲是相互獨立的,他們來自一個相同的總體,在這個總體中有一個未知的對藥物療效有效反應(positive response)的概率 \(\theta\),這樣的數據可以用一個二項分布似然來描述 binomial likelihood:

\[ p(y | n, \theta) = \binom{n}{y}\theta^y(1-\theta)^{n-y} \propto \theta^y(1-\theta)^{n-y} \]

  • 描述試驗前我們對\(\theta\)的認知的先驗概率 prior distribution,這是一個連續型先驗概率分布。
    對於百分比,我們用Beta分布來描述:

\[ \begin{aligned} \theta & \sim \text{Beta}(a,b) \\ p(\theta) & = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \theta^{a-1}(1-\theta)^{b-1}\\ \end{aligned} \]

根據貝葉斯定理,我們來把先驗概率分布和似然相結合(相乘),來獲取事後概率分布:

\[ \begin{aligned} p(\theta | y, n) & \propto p(y|\theta, n)p(\theta) \\ & \propto \theta^y(1-\theta)^{n-y}\theta^{a-1}(1-\theta)^{b-1} \\ & = \theta^{y + a -1}(1-\theta)^{n - y + b -1} \end{aligned} \] 眼尖的人立刻能看出來,這個事後概率分布本身也是一個Beta分布的概率方程,只是它的超參數和先驗概率相比發生了變化(更新):

\[ p(\theta | y,n) = \text{Beta}(a + y, b + n -y) \]

像這樣,先驗概率和事後概率兩個概率分布都來自相同家族的情況,先驗概率又被叫做和似然成共軛關系的先驗概率(共軛先驗概率, conjugate prior)。

par(mfrow=c(3,1))
# Plot exact prior probability density 
# values of the hyperparameters
a <- 9.2 
b <- 13.8

# prior function
prior <- Vectorize(function(theta) dbeta(theta, a, b))
# Plot 
curve(prior, 0, 1, type = "l", main = "Prior for "~theta, ylim = c(0, 4.5), frame = F,
      xlab = "Probability of positive response", ylab = "Density", lwd = 2, cex.axis = 1.5, cex.lab = 1.5)

# binomial likelihood function (likelihood)

Likelihood <- Vectorize(function(theta) dbinom(15, 20, theta))

# Plot
curve(Likelihood, 0, 1, type = "l", main = "Likelihood for the data", 
      frame = FALSE, xlab = "Probability of positive response", ylab = "Density", 
      lwd = 2, cex.axis =  1.5, cex.lab = 1.5)
# n <- 0; r <- 0; a <- 9.2; b <- 13.8; np <- 20
# plot(0:20, BetaBinom(0:20), type = "b", xlab = "r*", ylab = "P(R = r*)", 
#      main = "Prior predictive: a = 9.2, b = 13.8")

# Posterior function 

posterior <- Vectorize(function(theta) dbeta(theta, a+15, b+20-15))

# Plot

curve(posterior, 0, 1, type = "l", main = "Posterior for "~theta, 
      frame = FALSE, xlab = "Probability of positive response", ylab = "Density", 
      lwd = 2, cex.axis = 1.5, cex.lab = 1.5)
Prior, likelihood, and posterior for Drug example

圖 82.1: Prior, likelihood, and posterior for Drug example

本次試驗的模型,它的三個部分(先驗概率,似然,事後概率),分別從上到下繪制在圖 82.1 中。由於我們使用了共軛先驗概率,所以我們也可以通過數學的計算(甚至不需要計算機的輔助)也能算出事後概率分布。可是,並不是所有的模型都有共軛先驗概率分布供我們選擇,這時候,蒙特卡羅模擬試驗的算法就提供了強有力的工具。在BUGS語言中,我們可以用蒙特卡羅算法,忽視掉那些無法在數學上推導出事後概率分布方程的模型。BUGS本身很厲害,它可以自動識別出我們給它的先驗概率分布是否和似然之間是共軛的,如果是,那麼它會計算出共軛的事後概率分布方程,然後從事後概率分布方程中選取蒙特卡羅樣本。這個新藥試驗的BUGS模型可以寫作:

#  Monte Carlo model for Drug example

model{
    theta   ~ dbeta(9.2,13.8)          # prior distribution
    y       ~ dbin(theta,20)           # sampling distribution (likelihood)
    y       <- 15                      # data
}

你可以看到這個模型和我們在前一章做預測的模型只有第三行指令發生了變化。當時我們是打算要來做試驗結果的預測。此時,我們試驗完畢,觀察到15名患者的疼痛症狀得到了改善,所以試驗數據是15。BUGS本身會自動識別出我們是否給似然增加了觀察數據。當它識別到我們不是用這個模型做結果預測時,它會自動明白我們現在要來做事後概率分布的計算了。這個在似然裏面的數據,是它要拿來放到模型中做條件的(observed values, i.e. data needs to be conditioned on)。

82.2.1 事後概率分布預測

假如這個新藥的效果仍然無法讓人覺得信服,我們考慮再做一次試驗徵集更多的患者,如果在這個試驗中,40名患者中有25名或者更多的患者症狀得到緩解,可以考慮把該藥物加入下一次發展計劃當中。這時候,又一次回到了預測概率的問題上來,我們想知道,“再做40人的試驗時,有25名或者更多的患者的症狀可以得到緩解”這件事可能發生的概率。這時候的模型可以被擴展如下:

\[ \begin{split} \theta & \sim \text{Beta}(a,b) & \text{ Prior distribution} \\ y & \sim \text{Binomial}(\theta, n) & \text{ Sampling distribution} \\ y_{\text{pred}} & \sim \text{Binomial}(\theta, m) & \text{ Predictive distribution} \\ P_{\text{crit}} & \sim P(y_{\text{pred}} \geqslant m_{\text{crit}}) & \text{ Probability of exceeding critical threshold} \end{split} \]

這段模型翻譯成BUGS語言可以描述爲:

model{
  theta     ~ dbeta(a, b)                  # prior distribution
  y         ~ dbin(theta, n)               # sampling distribution
  y.pred    ~ dbin(theta, m)               # predictive distribution
  P.crit   <- step(y.pred - mcrit + 0.5)   # = 1 if y.pred >= mcrit, 0 otherwise
}

我們可以把數據寫在另一個txt文件裏面:

list(a = 9.2,             # parameters of prior distribution 
     b = 13.8, 
     y = 15,              # number of successes in completed trial
     n = 20,              # number of patients in completed trial
     m = 40,              # number of patients in future trial 
     mcrit = 25)          # critical value of future successes

當然這是一個很簡單的例子,你完全可以把數據和模型寫在一起:

model{
  theta     ~ dbeta(9.2, 13.8)                  # prior distribution
  y         ~ dbin(theta, 20)                   # sampling distribution
  y.pred    ~ dbin(theta, 40)                   # predictive distribution
  P.crit   <- step(y.pred - 24.5)               # = 1 if y.pred >= mcrit, 0 otherwise
  y        <- 15
}

```{r R-OpenBUGS08, cache=TRUE, fig.width=7, fig.height=3, fig.cap=‘Posterior and predictive distributions for Drug example’, fig.align=‘center’, out.width=‘80%’, message=TRUE, warning=FALSE}