Bayesian

輝瑞的疫苗功效研究設計中使用了哪種統計模型?

  • November 17, 2020

我知道這裡有一個類似的問題:

如何計算90%效力疫苗的95% CI?

但目前還沒有答案。另外,我的問題不同:另一個問題詢問如何使用 R 包中的函數計算 VE。我想知道為什麼疫苗功效的定義如本頁底部所示:

$$ \text{VE} = 1 - \text{IRR} $$

在哪裡

$$ \text{IRR} = \frac{\text{illness rate in vaccine group}}{\text{illness rate in placebo group}} $$

這就是它背後的統計模型。

我的嘗試:我認為這些研究適合邏輯回歸模型和單個二元預測變量 $ X $ ,識別接種疫苗的受試者( $ X=1 $ ) 或不 ( $ X=0 $ ):

$ p(Y|X) = \frac{1}{1+\exp{-(\beta_0 +\beta_1 X)}} $

然而,情況顯然不是這樣,因為對於 Moderna 疫苗,我們知道疫苗組有 5 例,安慰劑組有 90 例,這對應於 $ \text{VE} $ 的 $ 94.\bar{4}% $ . 僅這些數據就足以確定 $ \text{VE} $ ,但肯定它們不足以擬合 LR 模型,因此確定 $ \beta_1 $ .

此外,通過查看Pfizer 文檔的第111-113頁,似乎執行了不同的(貝葉斯?)分析。同樣,點估計似乎是 $ \text{VE} = 1 - \text{IRR} $ ,但提到了測試的力量,並提供了兩個表 7 和 8,它們顯示了成功和失敗的概率。你能告訴我如何在這些表格中獲得結果嗎?

效率與疾病風險比的關係

我想知道為什麼疫苗功效的定義如本頁底部所示:

$$ \text{VE} = 1 - \text{IRR} $$

在哪裡

$$ \text{IRR} = \frac{\text{illness rate in vaccine group}}{\text{illness rate in placebo group}} $$

這只是一個定義。可能以下表達式可以幫助您對它有不同的直覺

$$ \begin{array}{} VE &=& \text{relative illness rate reduction}\ &=& \frac{\text{change (reduction) in illness rate}}{\text{illness rate}}\ &=& \frac{\text{illness rate in placebo group} -\text{illness rate in vaccine group}}{\text{illness rate in placebo group}}\ &=& 1-IRR \end{array} $$

使用邏輯回歸建模

僅這些數據就足以確定 $ \text{VE} $ ,但肯定它們不足以擬合 LR 模型,因此確定 $ \beta_1 $ .

注意

$$ \text{logit}(p(Y|X)) = \log \left( \frac{p(Y|X)}{1-p(Y|X)} \right) = \beta_0 + \beta_1 X $$

並給出了兩個觀察結果 $ \text{logit}(p(Y|X=0)) $ 和 $ \text{logit}(p(Y|X=1)) $ 兩個參數 $ \beta_0 $ 和 $ \beta_1 $ 可以計算

R代碼示例:

請注意以下代碼cbind在 glm 函數中的使用。有關輸入此內容的更多信息,請參見此處的答案

vaccindata <- data.frame(sick    = c(5,90), 
                        healthy = c(15000-5,15000-90),
                        X       = c(1,0) 
                       )
mod <- glm(cbind(sick,healthy) ~ X, family = binomial, data = vaccindata)
summary(mod)

這給出了結果:

Call:
glm(formula = cbind(sick, healthy) ~ X, family = binomial, data = vaccindata)

Deviance Residuals: 
[1]  0  0

Coefficients:
           Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -5.1100     0.1057 -48.332  < 2e-16 ***
X            -2.8961     0.4596  -6.301 2.96e-10 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

   Null deviance: 9.2763e+01  on 1  degrees of freedom
Residual deviance: 2.3825e-12  on 0  degrees of freedom
AIC: 13.814

Number of Fisher Scoring iterations: 3

所以參數 $ \beta_1 $ 估計為 $ -2.8961 $ 有標準差 $ 0.4596 $

由此,您可以計算(估計)機率、效率及其​​置信區間。另請參閱:Moderna 和 Pfizer 疫苗試驗中的“有效性”究竟是如何估計的?

貝葉斯模型(表 6)

表 6

此外,通過查看Pfizer 文檔的第111-113頁,似乎執行了不同的(貝葉斯?)分析。同樣,點估計似乎是 $ \text{VE} = 1 - \text{IRR} $ ,但提到了測試的力量,並提供了兩個表 7 和 8,它們顯示了成功和失敗的概率。你能告訴我如何在這些表格中獲得結果嗎?

這些分析是在早期階段進行的,以根據結果驗證疫苗是否有效。表格給出了假設觀察結果,它們將達到臨界點,宣布失敗(成功的後驗概率<5%)或巨大成功(VE>30% 的概率大於 0.995)。

這些臨界點的百分比實際上是基於控制 I 類錯誤(更多內容見下文)。他們控制了總體的I 類錯誤,但不清楚這是如何在多個通過/不通過點之間分佈的。

考慮的結果是所有感染者中接種疫苗者的比率/數量。以感染者總數為條件,該比率服從二項分佈*。有關在這種情況下計算後驗的更多詳細信息,請參閱:beta 先驗如何影響二項似然下的後驗

*這裡可能對此有疑問;我仍然需要為此找到一個鏈接;但是您可以根據兩組近似泊松分佈(更準確地說它們是二項式分佈)和觀察特定案例組合的概率來得出這一點 $ k $ 和 $ n-k $ 以達到為條件 $ n $ 病例總數為$$ \frac{\lambda_1^k e^{-\lambda_1}/k! \cdot \lambda_2^{n-k}e^{-\lambda_2}/(n-k)! }{\lambda_2^ne^{-(\lambda_1\lambda_2)}/n! } = {n \choose k} \left(\frac{\lambda_1}{\lambda_1+\lambda_2}\right)^k \left(1- \frac{\lambda_1}{\lambda_1+\lambda_2}\right)^{n-l} $$

下圖顯示了此類計算的輸出圖

示例圖

  • 成功邊界 這是由值的後驗分佈計算的$$ \begin{array}{}\theta &=& (1-VE)/(2-VE)\ &=& RR/(1-RR) \&=& \text{vaccinated among infected}\end{array} $$ 例如,在前 32 名感染者中有 6 名接種疫苗和 26 名安慰劑的情況下,後驗是 Beta 分佈,參數為 0.7+6 和 1+26,累積分佈為 $ \theta < (1-0.3)/(2-0.3) $ 將會 $ \approx 0.996476 $ 對於 7 個接種疫苗和 25 個安慰劑,它將是 0.989,低於該水平。在 R 中,您可以將這些數字計算為pbeta(7/17,0.700102+6,1+26)
  • 無效邊界 為此,他們計算成功的概率,這是測試的力量。假設對於給定的假設,測試標準可以是在前 164 個病例中觀察疫苗組中的 53 個或更少的病例。然後,作為真實 VE 的函數,您可以估計通過測試的可能性。

在表 6 中,他們不是將其計算為單個 VE 的函數,而是作為 VE 的後驗分佈的積分或 $ \theta $ (還有這個 $ \theta $ 是β分佈,測試結果將是β-二項分佈)。似乎他們使用了以下內容:

### predict the probability of success (observing 53 or less in 164 cases at the end)
### k is the number of infections from vaccine
### n is the total number of infections
### based on k and n the posterior distribution can be computed
### based on the posterior distribution (which is a beta distribution)
### we can compute the success probability

predictedPOS <- function(k,n) {
  #### posterior alpha and beta
  alpha = 0.7+k
  beta = 1+n-k
  ### dispersion and mean
  s = alpha + beta
  m = alpha/(alpha+beta)
  ### probability to observe 53 or less out of 164 in final test
  ### given we allread have observed k out of n (so 53-k to go for the next 164-n infections)
  POS <- rmutil::pbetabinom(53-k,164-n,m,s)
  return(POS)
}

# 0.03114652
predictedPOS(15,32)
# 0.02486854
predictedPOS(26,62)
# 0.04704588
predictedPOS(35,92)

# 0.07194807
predictedPOS(14,32)
# 0.07194807
predictedPOS(25,62)
# 0.05228662
predictedPOS(34,92)

值 14、25、34 是後驗 POS 仍高於 0.05 的最高值。對於值 15、26、35,它在下面。

控制 I 類錯誤(表 7 和表 8)

表 7 和表 8

表 7 和表 8 給出了給定特定 VE 的成功概率分析(它們顯示為 30、50、60、70、80%)。它給出了在中期分析或最終分析期間分析通過成功標準的概率。

第一列很容易計算。它是二項式分佈的。例如,第一列中的概率 0.006, 0.054, 0.150, 0.368, 0.722 是在以下情況下出現 6 個或更少案例的概率 $ p=(100-VE)/(200-VE) $ 和 $ n = 32 $ .

其他列不是類似的二項分佈。如果在早期分析期間沒有成功,它們代表達到成功標準的概率。我不確定他們是如何計算的(他們指的是統計分析計劃 SAP,但不清楚在哪裡可以找到它以及它是否是開放訪問的)。但是,我們可以用一些 R 代碼來模擬它

### function to simulate succes for vaccine efficiency analysis
sim <- function(true_p = 0.3) {
 p <- (1-true_p)/(2-true_p)
 numbers <- c(32,62,92,120,164)
 success <- c(6,15,25,35,53)
 failure <- c(15,26,35)
 n <- c()
 ### simulate whether the infection cases are from vaccine or placebo group
 n[1] <- rbinom(1,numbers[1],p)
 n[2] <- rbinom(1,numbers[2]-numbers[1],p)
 n[3] <- rbinom(1,numbers[3]-numbers[2],p)
 n[4] <- rbinom(1,numbers[4]-numbers[3],p)
 n[5] <- rbinom(1,numbers[5]-numbers[4],p)

 ### days with succes or failure
 s <- cumsum(n) <= success
 f <- cumsum(n)[1:3] >= failure
 
 ### earliest day with success or failure
 min_s <- min(which(s==TRUE),7)
 min_f <- min(which(f==TRUE),6)
 
 ### check whether success occured before failure
 ### if no success occured then it has value 7 and will be highest
 ### if no failure occured then it will be 6 and be highest unless no success occured either
 result <- (min_s<min_f)
 
 return(result)
}

### compute power (probability of success)
### for different efficienc<y of vaccine
set.seed(1)
nt <- 10^5
x <- c(sum(replicate(nt,sim(0.3)))/nt,
      sum(replicate(nt,sim(0.5)))/nt,
      sum(replicate(nt,sim(0.6)))/nt,
      sum(replicate(nt,sim(0.7)))/nt,
      sum(replicate(nt,sim(0.8)))/nt)
x

這給出了 0.02073 0.43670 0.86610 0.99465 0.99992 接近最後一列中成功的總體概率。

雖然他們使用貝葉斯分析來計算表 6 中的值。他們已經選擇了邊界,根據控制 I 類錯誤(我認為他們使用在 VE = 0.3 的情況下取得成功的概率)來執行貝葉斯分析, p=0.021, 作為 I 類錯誤的基礎。這意味著如果真正的 VE = 0.3,那麼他們可能仍然錯誤地以 0.021 的概率宣布成功,如果真正的 VE<0.3,這種 I 類錯誤將是偶數較少的)

引用自:https://stats.stackexchange.com/questions/496774

comments powered by Disqus