Wald Confidence Interval for a Difference of Binomial Proportions Based on Paired Data

2x2 table에서 민감도와 특이도는 본래 McNemar’s test를 통해 $z=\frac{(f_{21}-f_{12})}{\sqrt{f_{12}+f_{21}}}$로 정의된 z값의 p-value를 확인하는 방법으로 진행되었습니다. 하지만 위 방식으로는 민감도와 특이도의 신뢰 구간에 대한 정보를 알 수 없기에, 신뢰 구간까지 정할 수 있는 방법인 Wald confidence interval을 적용함으로써 더 많은 정보를 얻을 수 있습니다.

개요

이 문서에서는 진단 테스트 1과 테스트 2의 민감도 차이를 기반으로 신뢰 구간(CI)을 계산하는 방법을 다룹니다. 이를 위해 두 테스트의 결과를 기반으로 2x2 테이블을 구성하고, 각 민감도를 계산합니다. 최종적으로 두 민감도의 차이에 대한 신뢰 구간을 계산하는 수식을 도출합니다.

1. 기본 개념

각 테스트의 민감도 $p_1$과 $p_2$는 아래와 같이 정의됩니다:

• $p_1$은 테스트 1이 진정한 양성을 탐지한 비율입니다.
• $p_2$은 테스트 2가 진정한 양성을 탐지한 비율입니다.

2x2 테이블은 아래와 같이 구성될 수 있습니다:

$$\begin{array}{cccc}& \text{Test 2 Positive}& \text{Test 2 Negative}& \text{Total}\\ \text{Test 1 Positive}& a& b& a+b\\ \text{Test 1 Negative}& c& d& c+d\\ \text{Total}& a+c& b+d& n\end{array}$$

여기서 각 테스트의 민감도는 다음과 같이 계산됩니다.(특이도 역시 모든 데이터들이 양성인 것을 음성으로 바꾸어 생각하면 아래 과정이 동일해집니다)

$$p_1=\frac{a+b}{n},p_2=\frac{a+c}{n}$$

2. 신뢰 구간의 도출

두 민감도의 차이에 대한 신뢰 구간을 계산하기 위해 Wald 방법을 사용할 수 있습니다. Wald 신뢰 구간은 아래와 같이 정의됩니다:

2.1. 분산의 기본 식

두 확률 변수 $p_1$과 $p_2$의 차이의 분산은 다음과 같이 계산됩니다:

$$\text{Var}(p_2-p_1)=\text{Var}(p_2)+\text{Var}(p_1)-2\text{Cov}(p_2,p_1)$$

2.2. 각 분산의 계산

우리는 $p_1$과 $p_2$를 비율로 간주할 수 있습니다. 이 비율들의 분산은 다음과 같이 계산됩니다:

$$\text{Var}(p_1)=\frac{p_1(1-p_1)}{n}=\frac{p_{12}+p_{11}}{n}\cdot (1-\frac{p_{12}+p_{11}}{n})$$

$$\text{Var}(p_2)=\frac{p_2(1-p_2)}{n}=\frac{p_{21}+p_{11}}{n}\cdot (1-\frac{p_{21}+p_{11}}{n})$$

2.3. 공분산 $\text{Cov}(p_2,p_1)$ 의 계산

$p_2$ 와 $p_1$ 의 공분산은 다음과 같이 정의됩니다:

$$\text{Cov}(p_2,p_1)=\text{Cov}(p_{11}+p_{21},p_{11}+p_{12})$$

이를 확장하면 다음과 같은 네 가지 항으로 분리할 수 있습니다:

$$\text{Cov}(p_2,p_1)=\text{Cov}(p_{11},p_{11})+\text{Cov}(p_{11},p_{12})+\text{Cov}(p_{21},p_{11})+\text{Cov}(p_{21},p_{12})$$

공분산은 다음과 같이 정리할 수 있습니다.

$$\text{Var}(p_{11})=\frac{p_{11}(1-p_{11})}{n}$$

$$\text{Cov}(p_{11},p_{12})=-\frac{p_{11}{p}_{12}}{n}$$

$$\text{Cov}(p_{11},p_{21})=-\frac{p_{11}{p}_{21}}{n}$$

$$\text{Cov}(p_{21},p_{12})=-\frac{p_{21}{p}_{12}}{n}$$

이를 바탕으로 $\text{Cov}(p_2,p_1)$를 계산하면:

$$\text{Cov}(p_2,p_1)=\frac{p_{11}(1-p_{11}-p_{12}-p_{21})-p_{21}{p}_{12}}{n}$$

여기서 $p_{11}+p_{12}+p_{21}+p_{22}=1$임을 이용하여, 공분산을 다음과 같이 단순화할 수 있습니다:

$$\text{Cov}(p_2,p_1)=\frac{p_{11}{p}_{22}-p_{21}{p}_{12}}{n}$$

2.4. 최종적으로 $\text{Var}(p_2-p_1)$ 의 유도

이제 분산 공식을 대입하여 최종적으로 $p_2-p_1$의 분산을 계산합니다:

$$\text{Var}(p_2-p_1)=\text{Var}(p_2)+\text{Var}(p_1)-2\text{Cov}(p_2,p_1)$$

따라서 $p_2-p_1$의 분산은 다음과 같습니다:

$$\text{Var}(p_2-p_1)=\frac{(p_{12}+p_{21})-(p_{21}-p_{12}{)}^2/n}{n}$$

$$\text{Var}(p_2-p_1)=\frac{(b+c)-(b-c{)}^2/n}{n}$$

이를 적용한 민감도의 차이에 대한 신뢰 구간은 다음과 같이 계산됩니다. $$C{I}_{1-\alpha /2}(\hat{\theta })=[\hat{\theta }\pm z_{1-\alpha /2}\cdot \frac{1}{n}\cdot \sqrt{b+c-\frac{(b-c{)}^2}{n}}]$$

  if ( (ci.method == "wald") & (cont.corr == FALSE) ) {
    # sensitivity
    b <- tab$diseased[1,2]; c <- tab$diseased[2,1]; n <- tab$diseased[3,3]
    sens.diff.se <- sqrt((b+c) - ((b-c)**2) / n) / n
    sens.diff.cl <- sens.diff + c(-1,1) * qnorm(1-alpha/2) * sens.diff.se}

만약 이항분포로 얻어진 값을 정규분포에 맞도록 조정하기 위해 continuity correction을 진행한다면 wald에서는 확률 $p$ 하나당 $\frac{1}{2n}$만큼 분산을 늘려, 아래와 같은 값을 갖게 된다.

$$C{I}_{1-\alpha /2}(\hat{\theta })=[\hat{\theta }\pm (z_{1-\alpha /2}\cdot \frac{1}{n}\sqrt{b+c-\frac{(b-c{)}^2}{n}}+\frac{1}{n})].$$

  if ( (ci.method == "wald") & (cont.corr == TRUE) ) {
    # sensitivity
    b <- tab$diseased[1,2]; c <- tab$diseased[2,1]; n <- tab$diseased[3,3]
    sens.diff.se <- (sqrt((b+c) - ((b-c)**2) / n) / n) + 1/n
    sens.diff.cl <- sens.diff + c(-1,1) * qnorm(1-alpha/2) * sens.diff.se}

3. Wald 이외 신뢰구간 계산 방법

3.1. Agresti 신뢰 구간(Agresti CI)

Wald 신뢰 구간을 수정한 방법 중 하나로, 샘플 크기에 특정 상수를 더한 후 위 wald와 같은 방식으로 전체 샘플 크기 𝑛에 1에서 4까지의 상수를 더한 후 이 방법들을 비교했습니다. 그 결과, n+2를 사용하는 것이 표준 Wald 신뢰 구간과 비교했을 때 표본 증가의 효과로 신뢰 구간의 성능(coverage probability)을 향상시킴을 확인했습니다. $$C{I}_{1-\alpha /2}(\hat{\theta })=[\hat{\theta }\pm z_{1-\alpha /2}\cdot \frac{1}{n+2}\cdot \sqrt{(b+0.5)+(c+0.5)-\frac{(b-c{)}^2}{n+2}}]$$

  if (ci.method == "agresti-min") {
    k <- 0.5
    # sensitivity    
    b <- tab$diseased[1,2]+k; c <- tab$diseased[2,1]+k; n <- tab$diseased[3,3]+4*k
    sens.diff.se <- (sqrt((b+c) - ((b-c)**2) / n) / n) 
    sens.diff.cl <- sens.diff + c(-1,1) * qnorm(1-alpha/2) * sens.diff.se}

3.2. Tango

Tango는 귀무가설과 함수는 두 비율 간의 차이와, 주어진 정보 행렬($b$, $c$, $n$)을 기반으로 likelihood를 계산해 신뢰구간을 계산하는 방법입니다. R에서는 scoreci.mp라는 별도의 함수를 사용해 진행합니다.

  if (ci.method == "tango") {
    # sensitivity    
    b <- tab$diseased[1,2]; c <- tab$diseased[2,1]; n <- tab$diseased[3,3]
    tango <- scoreci.mp(b, c, n, conf.level=1-alpha)    
    sens.diff.se <- NA    
    sens.diff.cl <- sort(c(tango$conf.int[1], tango$conf.int[2]))
    if ( (tango$conf.int[1] > sens.diff) | (tango$conf.int[2] < sens.diff))
      sens.diff.cl <- sort(-1*sens.diff.cl)}

4. 예시 코드

4.1. 일반적인 민감도, 특이도 검정

library(DTComPair)

Loading required package: PropCIs

t1 <- read.tab.paired(18, 14, 0, 18,
                      18, 12, 2, 18)
t1

Two binary diagnostic tests (paired design)

Test1: 'Noname 1'
Test2: 'Noname 2'

Diseased:
           Test1 pos. Test1 neg. Total
Test2 pos.         18         14    32
Test2 neg.          0         18    18
Total              18         32    50

Non-diseased:
           Test1 pos. Test1 neg. Total
Test2 pos.         18         12    30
Test2 neg.          2         18    20
Total              20         30    50

sesp.diff.ci(t1, ci.method="wald", cont.corr=FALSE)

$sensitivity
     test1      test2       diff    diff.se   diff.lcl   diff.ucl 
0.36000000 0.64000000 0.28000000 0.06349803 0.15554615 0.40445385 

$specificity
      test1       test2        diff     diff.se    diff.lcl    diff.ucl 
 0.60000000  0.40000000 -0.20000000  0.06928203 -0.33579029 -0.06420971 

$ci.method
[1] "wald"

$alpha
[1] 0.05

$cont.corr
[1] FALSE

sesp.diff.ci(t1, ci.method="wald", cont.corr=TRUE)

$sensitivity
     test1      test2       diff    diff.se   diff.lcl   diff.ucl 
0.36000000 0.64000000 0.28000000 0.08349803 0.11634687 0.44365313 

$specificity
      test1       test2        diff     diff.se    diff.lcl    diff.ucl 
 0.60000000  0.40000000 -0.20000000  0.08928203 -0.37498957 -0.02501043 

$ci.method
[1] "wald"

$alpha
[1] 0.05

$cont.corr
[1] TRUE

sesp.diff.ci(t1, ci.method="agresti-min")

$sensitivity
     test1      test2       diff    diff.se   diff.lcl   diff.ucl 
0.36000000 0.64000000 0.28000000 0.06444681 0.15368658 0.40631342 

$specificity
      test1       test2        diff     diff.se    diff.lcl    diff.ucl 
 0.60000000  0.40000000 -0.20000000  0.06954236 -0.33630053 -0.06369947 

$ci.method
[1] "agresti-min"

$alpha
[1] 0.05

$cont.corr
[1] FALSE

sesp.diff.ci(t1, ci.method="tango")

$sensitivity
    test1     test2      diff   diff.se  diff.lcl  diff.ucl 
0.3600000 0.6400000 0.2800000        NA 0.1747417 0.4166512 

$specificity
      test1       test2        diff     diff.se    diff.lcl    diff.ucl 
 0.60000000  0.40000000 -0.20000000          NA -0.34470882 -0.06111243 

$ci.method
[1] "tango"

$alpha
[1] 0.05

$cont.corr
[1] FALSE

아래와 같은 방법으로 answer, test1, test2를 열로 가지고, 각 결과 데이터가 1,0(1은 diseased sample, 2는 non-diseased sample)로 표시되어 있는 경우에는 바로 paired table을 제작할 수 있습니다.

library(DTComPair)
tb <- tab.paired(answer, test1, test2, data = na.omit(sample_data))

4.2. 민감도, 특이도의 비열등성 검정

어떤 테스트의 민감도와 특이도가 대조 테스트에 비해 떨어지지 않는다는 것을 확인하기 위해서는 비열등성 검정을 사용해야 합니다. 이 때는 sesp.diff.ci 함수를 통해 구한 standard error 값과 두 민감도/특이도의 차이를 통해 p value를 계산합니다. 아래 예시에서는 민감도에 대한 비열등성 마진(sens_margin)을 5%, 특이도에 대한 비열등성 마진(spec_margin)을 10%로 설정했습니다.
민감도를 예로 들자면 귀무가설은 아래와 같고, 아래 수식에 따라 p value를 계산합니다.

$$H_0:\text{test1 sensitivity}-\text{test2 sensitivity}\le -\text{ sensitivity margin}$$ $$p-value=1-\Phi \left(\frac{\text{test1 sensitivity}-\text{test2 sensitivity}+\text{sensitivity margin}}{\text{sensitivity diff.SE}}\right)$$

library(DTComPair)
t1 <- read.tab.paired(18, 14, 0, 18,
                      18, 12, 2, 18)
t1.wald <- sesp.diff.ci(t1, ci.method="wald", cont.corr=FALSE)

sens_margin <-  0.05
spec_margin <-  0.1

p_value_sensitivity <- pnorm((t1.wald$sensitivity['diff'] + sens_margin) / t1.wald$sensitivity['diff.se'], lower.tail = FALSE)
p_value_specificity <- pnorm((t1.wald$specificity['diff'] + spec_margin) / t1.wald$specificity['diff.se'], lower.tail = FALSE)
p_value_sensitivity

        diff 
1.012589e-07 

p_value_specificity

     diff 
0.9255427 

따라서 민감도에 대해서는 위 귀무가설이 기각 되었음으로 test1의 민감도가 test2보다 열등하지 않다고 말할 수 있지만, 특이도에 대해서는 귀무가설이 기각되지 않았기 때문에 test1의 특이도가 test2에 비해 열등하지 않다는 결론을 내릴 수 없습니다.