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

Paired data에서 wald 신뢰 구간을 구하는 방법을 통해 민감도와 특이도의 신뢰 구간을 구하는 방법을 소개합니다.

R
statistics
Author

Hyojong Myung

Published

August 26, 2024

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}{c|c|c|c} & \text{Test 2 Positive} & \text{Test 2 Negative} & \text{Total} \\ \hline \text{Test 1 Positive} & a & b & a + b \\ \hline \text{Test 1 Negative} & c & d & c + d \\ \hline \text{Total} & a + c & b + d & n \\ \end{array}\]

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

\[p_1 = \frac{a + b}{n}, \quad 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 \left(1 - \frac{p_{12} + p_{11}}{n}\right) \]

\[ \text{Var}(p_2) = \frac{p_2(1 - p_2)}{n} = \frac{p_{21} + p_{11}}{n} \cdot \left(1 - \frac{p_{21} + p_{11}}{n}\right) \]

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} \]

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

  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}\)만큼 분산을 늘려, 아래와 같은 값을 갖게 된다.

\[ CI_{1-\alpha/2}(\hat{\theta}) = \left[ \hat{\theta} \pm \left( z_{1-\alpha/2} \cdot \frac{1}{n} \sqrt{b + c - \frac{(b - c)^2}{n} }+ \frac{1}{n} \right) \right]. \]

  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)을 향상시킴을 확인했습니다. \[ \ CI_{1-\alpha/2}(\hat{\theta}) = \left[ \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}} \right] \ \]

  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. 일반적인 민감도, 특이도 검정

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 \colon \text{test1 sensitivity} - \text{test2 sensitivity} \leq -\text{ sensitivity margin}\] \[p-value = 1 - Φ\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에 비해 열등하지 않다는 결론을 내릴 수 없습니다.

Reuse

Citation

BibTeX citation:
@online{myung2024,
  author = {Myung, Hyojong},
  title = {Wald {Confidence} {Interval} for a {Difference} of {Binomial}
    {Proportions} {Based} on {Paired} {Data}},
  date = {2024-08-26},
  url = {https://blog.zarathu.com/posts/2024-08-26-sensspec-waldinterval/},
  langid = {en}
}
For attribution, please cite this work as:
Myung, Hyojong. 2024. “Wald Confidence Interval for a Difference of Binomial Proportions Based on Paired Data.” August 26, 2024. https://blog.zarathu.com/posts/2024-08-26-sensspec-waldinterval/.