8.3 두 군의 분율(proportion) 비교
각 군의 표본 크기(sample size, \(n_1, n_2\))가 크다면 정규분포로 근사해서 구할 수도 있을 것이나, 요즘은 컴퓨터가 발달하여 이항 분포를 사용하여 더 정확히 계산할 수 있다.
두 군의 모비율 \(\pi_1, \pi_2\)에 대한 추정치 \(p_1, p_2\)를 이용하여 다음과 같이 SE를 구하고 이를 이용하여 신뢰구간을 구할 수 있다.
\(\hat{\pi_1} = p_1 = \frac{y_1}{n_1}\)
\(\hat{\pi_2} = p_2 = \frac{y_2}{n_2}\)
\(\hat{SE} = \sqrt{\frac{p_1 ( 1 - p_1)}{n_1} + \frac{p_2 ( 1 - p_2)}{n_2}}\)
신뢰구간은 다음과 같다.
\((p_1 - p_2) \pm z_{\alpha/2} \hat{SE}\)
이런 신뢰구간을 Wald confidence interval (CI) 이라 한다. Wald CI는 소표본이거나 확률이 0이나 1에 가까울 때 성능이 나빠진다. 여기에 대한 대안으로 모든 cell에 1을 더하고 계산하거나, Score CI를 구하는 방법이 있다.
두 군의 비율이 있을 때 차이를 구할 수도 있지만, 상대위험도(Relative Risk, RR)나 오즈비(Odds Ratio, OR)을 구할 때도 있다. 각각의 정의와 표준오차 공식은 다음과 같다.
- 상대위험도(Relative Risk, RR)
\(RR = \frac{\pi_1}{\pi_2}\)
\(\hat{RR} = \frac{p_1}{p_2}\)
\(\hat{SE}(log \hat{RR}) = \sqrt{\frac{1}{y_1} - \frac{1}{n_1} + \frac{1}{y_2} - \frac{1}{n_2}}\)
- 오즈비(Odds Ratio, OR)
\(OR = \frac{\pi_1 /(1 - \pi_1)}{\pi_2 /(1 - \pi_2)}\)
\(\hat{OR} = \frac{p_1 / (1 - p_1)}{p_2 / (1 - p_2)}\)
\(\hat{SE}(log \hat{OR}) = \sqrt{\frac{1}{n_1} + \frac{1}{y_1 - n_1} + \frac{1}{n_2} + \frac{1}{y_2 - n_2}}\)
즉, RR은 분율의 비이고, OR은 비율의 비이다.
위의 3가지를 구하는 함수의 내용은 다음과 같다.
RD = function(y1, n1, y2, n2, conf.level=0.95) # Risk Difference
{
p1 = y1/n1
p2 = y2/n2
pe = p1 - p2
se = sqrt(p1*(1 - p1)/n1 + p2*(1 - p2)/n2) # SE of pe
z.crit = qnorm(0.5 + conf.level/2)
lower = pe - z.crit*se
upper = pe + z.crit*se
Res = data.frame(p1=p1, p2=p2, RD=pe, SE=se, lower=lower, upper=upper)
return(Res)
}
RR = function(y1, n1, y2, n2, conf.level=0.95) # Relative Risk
{
p1 = y1/n1
p2 = y2/n2
pe = p1/p2
selog = sqrt(1/y1 - 1/n1 + 1/y2 - 1/n2) # SE of log(pe)
z.crit = qnorm(0.5 + conf.level/2)
lower = exp(log(pe) - z.crit*selog)
upper = exp(log(pe) + z.crit*selog)
Res = data.frame(p1=p1, p2=p2, RR=pe, SElog=selog, lower=lower, upper=upper)
return(Res)
}
OR = function(y1, n1, y2, n2, conf.level=0.95) # Odds Ratio
{
o1 = y1/(n1 - y1) # odd 1
o2 = y2/(n2 - y2) # odd 2
pe = o1/o2
selog = sqrt(1/y1 + 1/(n1 - y1) + 1/y2 + 1/(n2 - y2)) # SE of log(pe)
z.crit = qnorm(0.5 + conf.level/2)
lower = exp(log(pe) - z.crit*selog)
upper = exp(log(pe) + z.crit*selog)
Res = data.frame(odd1=o1, odd2=o2, OR=pe, SElog=selog, lower=lower, upper=upper)
return(Res)
}[예제 8-3] NEJM 318:262-264 (1988) 에 출판된 논문으로 아스피린 복용과 심근경색증(myocardial infarction, MI) 발생에 대한 연구이며 결과는 다음과 같았다. RD, RR, OR을 추정하시오.
| Group | MI occured | No MI | Total |
|---|---|---|---|
| Aspirin | 104 | 10,933 | 11,037 |
| Placebo | 189 | 10,845 | 11,034 |
- Wald CI and test for Difference
RD(104, 11037, 189, 11034) # no continuity correction p1 p2 RD SE lower upper
1 0.00942285 0.01712887 -0.007706024 0.001539964 -0.0107243 -0.004687751RR(104, 11037, 189, 11034) # Active first, Placebo later p1 p2 RR SElog lower upper
1 0.00942285 0.01712887 0.550115 0.1213473 0.4336731 0.6978217OR(104, 11037, 189, 11034) # Active first, Placebo later odd1 odd2 OR SElog lower upper
1 0.009512485 0.01742739 0.5458355 0.1228416 0.429041 0.694424prop.test(c(104, 189), c(11037, 11034)) # default will use continuity correction.
2-sample test for equality of proportions with continuity correction
data: c(104, 189) out of c(11037, 11034)
X-squared = 24.429, df = 1, p-value = 7.71e-07
alternative hypothesis: two.sided
95 percent confidence interval:
-0.010814914 -0.004597134
sample estimates:
prop 1 prop 2
0.00942285 0.01712887 - Score CI and test for Difference, RR, and OR
RDmn1(104, 11037, 189, 11034) # Risk Difference p1 p2 RD lower upper
0.009422850 0.017128874 -0.007706024 -0.010788543 -0.004716840 RRmn1(104, 11037, 189, 11034) # RR (active first) p1 p2 RR lower upper
0.00942285 0.01712887 0.55011498 0.43322316 0.69847602 ORmn1(104, 11037, 189, 11034) # OR (active first) odd1 odd2 OR lower upper
0.009512485 0.017427386 0.545835457 0.429279021 0.694040668 단일군의 경우에 다음의 함수로 Score CI를 구할 수 있다.
ScoreCI = function(y, n, conf.level=0.95)
{
z = qnorm(0.5 + conf.level/2)
z2 = z^2
mp = (y + 0.5*z2)/(n + z2) # mid point but not point estimate
eb = z/(n + z2)*sqrt(y*(n - y)/n + z2/4) # error bound
Res = data.frame(Prop = y/n, lower = mp - eb, upper = mp + eb)
return(Res)
}round(ScoreCI(104, 11037)*100, 3) # in percent Prop lower upper
1 0.942 0.778 1.14위의 특이한 점은 비대칭 신뢰구간이라는 것이다.
그런데, 단일군의 경우에는 binom.test가 있으므로 굳이 위의 함수를 사용할 필요 없으나, 두 군의 비율의 차이에 대한 신뢰구간을 구하고자 할 때는 score CI가 많이 사용된다.
다음은 sasLM내의 관련 함수의 핵심내용이다.
RDmn1 = function(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)
{
p1 = y1/n1 # p of test (active) group
p2 = y2/n2 # p of control (placebo) group
RD = p1 - p2 # point estimate of risk difference (RD)
v0 = qchisq(conf.level, 1) # delta ofv for confidence interval
Obj = function(rd) { # find dp points of increased obj fx value (ofv) by v0
mLL = function(p1t) -dbinom(y1, n1, p1t, log=T) - dbinom(y2, n2, p1t - rd, log=T)
p1t = nlminb(p1, mLL, lower=max(0, RD), upper=min(1, RD + 1))$par
p2t = max(0, min(p1t - rd, 1)) # MLE p1t, p2t with fixed rd
var0 = (p1t*(1 - p1t)/n1 + p2t*(1 - p2t)/n2)*(n1 + n2)/(n1 + n2 - 1)
return(((rd - RD)^2/var0 - v0)^2) # find roots of increased ofv by v0
}
options(warn=-1)
LL = nlminb(max(eps - 1, RD - eps), Obj, lower=max(-1, RD - 1), upper=RD)$par
UL = nlminb(min(1 - eps, RD + eps), Obj, lower=RD, upper=min(RD + 1, 1))$par
options(warn=1)
return(data.frame(p1 = p1, p2 = p2, RD = RD, lower = LL, upper = UL))
}
RRmn1 = function(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)
{
p1 = y1/n1 # p of test (active) group
p2 = y2/n2 # p of control (placebo) group
RR = p1/p2 # point estimate of relative risk (RR)
v0 = qchisq(conf.level, 1) # delta ofv for confidence interval
Obj = function(rr) { # find rr points of increased obj fx value (ofv) by v0
mLL = function(p2t) -dbinom(y1, n1, p2t*rr, log=T) - dbinom(y2, n2, p2t, log=T)
p2t = nlminb(p2, mLL, lower=0, upper=1)$par
p1t = p2t*rr # MLE p1d, p2d with fixed ratio rr
var0 = (p1t*(1 - p1t)/n1 + RR*RR*p2t*(1 - p2t)/n2)*(n1 + n2)/(n1 + n2 - 1)
return(((p1t - p2t*RR)^2/var0 - v0)^2) # find the root of increased ofv by v0
}
options(warn=-1)
LL = nlminb(max(eps, RR - eps), Obj, lower=0, upper=RR)$par
UL = nlminb(RR + eps, Obj, lower=RR)$par
options(warn=1)
return(data.frame(p1 = p1, p2 = p2, RR = RR, lower = LL, upper = UL))
}
ORmn1 = function(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)
{
p1 = y1/n1
p2 = y2/n2
o1 = y1/(n1 - y1) # odd of test (active) group
o2 = y2/(n2 - y2) # odd of control (placebo) group
OR = o1/o2 # point estimate of odds ratio (OR)
v0 = qchisq(conf.level, 1) # delta ofv for confidence interval
Obj = function(or) { # find or points of increased obj fx value (ofv) by v0
mLL = function(p2t) {
p1t = p2t*or/(1 + p2t*(or - 1))
-dbinom(y1, n1, p1t, log=T) - dbinom(y2, n2, p2t, log=T)
}
p2t = nlminb(p2, mLL, lower=0, upper=1)$par
p1t = p2t*or/(1 + p2t*(or - 1)) # MLE p1d, p2d with fixed odds ratio or
div = (p1 - p1t)/(p1t*(1 - p1t)) - (p2 - p2t)/(p2t*(1 - p2t))
var0 = (1/(n1*p1t*(1 - p1t)) + 1/(n2*p2t*(1 - p2t)))*(n1 + n2)/(n1 + n2 - 1)
return((div^2/var0 - v0)^2) # find the root of increased ofv by v0
}
options(warn=-1)
LL = nlminb(max(eps, OR - eps), Obj, lower=0, upper=OR)$par
UL = nlminb(OR + eps, Obj, lower=OR)$par
options(warn=1)
return(data.frame(odd1 = o1, odd2 = o2, OR = OR, lower = LL, upper = UL))
}R의 PropCIs, PF package와 SAS에도 Miettinen와 Nurminen(이하 MN으로 표기)에 의한 score CI를 구해준다. 하지만, 이들의 방법은 risk difference(RD)의 경우에는 문제가 없지만, RR의 경우에는 3차 방정식으로 근사하여 해를 구하기 때문에 유효숫자가 2자리정도 밖에 안되며, OR의 경우에는 초기값(initial guess)에서 1.001배씩 증가 또는 감소시키면서 해를 찾기 때문에 유효숫자가 3자리 정도밖에 되지 않는다. 다른 software를 사용할 때는 이런 점에 유의한다.
NEJM에 2021-12-16 발표된 COVID-19치료제인 Molnupiravir 3상 임상시험의 자료분석에는 위의 MN방식을 사용하였고, 증상발현후 3일이내에 치료시작한 strata와 4일 또는 5일째에 치료시작한 strata로 나누어진 자료를 가중치(weighting)를 사용하여 common risk difference를 구하였다. 다음은 strata가 있는 자료를 RDmn을 사용하여 common RD를 구하는 예시이다.
d1 = data.frame(y1=c(25, 23), n1=c(339, 370), y2=c(28, 40), n2=c(335, 364)) ; d1 y1 n1 y2 n2
1 25 339 28 335
2 23 370 40 364RDmn(d1)$Strata
p1 p2 RD lower upper
1 0.07374631 0.08358209 -0.009835777 -0.05174840 0.031572839
2 0.06216216 0.10989011 -0.047727948 -0.08980705 -0.007330023
$Common
p1 p2 RD lower upper
0.067708565 0.097294030 -0.029585465 -0.058936760 -0.000855763 Percent 단위가 되려면 위의 결과에 100을 곱해야 한다.
Strata가 있는 자료에서 common RR을 구하는 경우에는 RRmn 함수를 사용하고, common OR을 구하는 경우에는 ORmn 함수를 사용하면 된다.
앞의 연구는 전향적 무작위 배정 연구(prospective randomized trial)이기 때문에 RR이 계산 가능하지만, 증례-대조군(case-control) 연구에서는 OR만 계산할 수 있고, RR은 계산할 수 없다.
발생빈도가 낮은 경우에는 OR을 RR의 근사값으로 사용할 수 있다. 만약 \(\pi_1, \pi_2\)가 0.2보다 클 때는 OR의 제곱근을 RR의 근사값으로 사용할 수 있다.
[예제 8-4] British Med J. Sep. 1950 739-748 에 출판된 흡연과 폐암의 증례 대조군 연구 결과이다. OR을 추정하시오.
| Smoker | Yes | No | Total |
|---|---|---|---|
| Case | 688 | 21 | 709 |
| Control | 650 | 59 | 709 |
OR(y1=688, n1=709, y2=650, n2=709) odd1 odd2 OR SElog lower upper
1 32.7619 11.01695 2.973773 0.2599234 1.786737 4.949427ORmn1(y1=688, n1=709, y2=650, n2=709) odd1 odd2 OR lower upper
32.761905 11.016949 2.973773 1.794878 4.925924