백신 효능에 대한 화이자 연구 설계에서 사용되는 통계 모델은 무엇입니까?
여기에 비슷한 질문이 있다는 것을 알고 있습니다.
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$.
또한 화이자 문서 111-113 페이지를 보면 다른 (Bayesian?) 분석이 수행되는 것처럼 보입니다. 다시 말하지만, 포인트 추정치는$ \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 코드 예 :
아래 코드는 cbindglm 함수에서 사용 합니다. 입력에 대한 자세한 내용은 여기에서이 답변을 참조 하십시오 .
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)
또한 화이자 문서 111-113 페이지를 보면 다른 (Bayesian?) 분석이 수행되는 것처럼 보입니다. 다시 말하지만, 포인트 추정치는$ \text{VE} = 1 - \text{IRR}$, 그러나 테스트의 힘이 언급되고 성공 및 실패 확률을 보여주는 두 개의 표 7과 8이 제시됩니다. 이러한 테이블에서 결과를 얻는 방법을 보여줄 수 있습니까?
이러한 분석은 결과가 주어지면 백신이 효과적인지 확인하기 위해 초기 단계에서 수행됩니다. 표는 실패 (사후 성공 확률 <5 %) 또는 대성공 (VE> 30 %가 0.995보다 클 확률)을 선언하기 위해 티핑 포인트에 도달 할 수있는 가설 적 관측치를 제공합니다.
티핑 포인트에 대한 이러한 백분율은 실제로 제 1 종 오류 제어를 기반으로합니다 (자세한 내용은 아래 참조). 전체적인 제 1 종 오류를 제어 하지만 이것이 여러 진행 / 불가 점에 어떻게 분산되는지는 명확하지 않습니다 .
고려되는 결과는 모든 감염자 중 예방 접종을받은 사람의 비율 / 수입니다. 총 감염자에 대한 조건부로이 비율은 이항 분포 *를 따릅니다. 이 경우 사후 계산에 대한 자세한 내용은 이항 우도에서 베타 사전 이 사후에 어떻게 영향을 미치는가를 참조하십시오.
* 여기에 질문이있을 것입니다. 여전히 이것에 대한 링크를 찾아야합니다. 그러나 두 그룹이 대략 포아송 분포 (더 정확하게는 이항 분포)라는 아이디어와 특정 케이스 조합을 관찰 할 확률을 바탕으로이를 도출 할 수 있습니다.$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 명의 위약의 경우, 후부는 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 또는 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의 값을 계산하지만, 제 1 종 오류를 제어하여 베이지안 분석을 수행 한 경계를 선택했습니다 (VE = 0.3 일 때 성공할 확률을 사용한다고 생각합니다. , p = 0.021, 제 1 종 오류의 기초가됩니다. 즉, 참 VE = 0.3이면 오류가 여전히 0.021 확률로 성공을 선언 할 수 있으며 참 VE <0.3이면이 제 1 종 오류가 짝수입니다. 적게)
이 모든 결과는 fisher의 정확한 테스트의 기본 R 구현에서 구현 된 조건부 최대 가능성 추정을 사용하는 것과 일치합니다.
splits <- matrix(c(6,26,15,47,25,67,35,85,53,111), ncol = 2, byrow = T)
total <- 43000
for(interim in 1:nrow(splits)) {
positive_vax <- splits[interim, 1]
positive_pla <- splits[interim, 2]
negative_vax <- (total / 2 ) - positive_vax
negative_pla <- (total / 2 ) - positive_pla
cont_tab <- matrix(c(positive_vax, positive_pla, negative_vax, negative_pla), nrow = 2)
test <- fisher.test(cont_tab)
VE <- 1 - test$estimate
print(paste(VE, "% (", positive_vax, ":", positive_pla, ")"))
}
결과:
[1] "0.769425572629548 % ( 6 : 26 )"
[1] "0.681342630733629 % ( 15 : 47 )"
[1] "0.627606975573189 % ( 25 : 67 )"
[1] "0.589208653283242 % ( 35 : 85 )"
[1] "0.523803347975998 % ( 53 : 111 )"