輝瑞的疫苗功效研究設計中使用了哪種統計模型?
我知道這裡有一個類似的問題:
但目前還沒有答案。另外,我的問題不同:另一個問題詢問如何使用 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)
此外,通過查看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 給出了給定特定 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 類錯誤將是偶數較少的)