Stratified Cox 비례위험모형

개념 : 비례위험 가정을 만족하지 않는 변수의 경우, 층화를 하여서 분석하여야함

비례 위험 가정이란 ? => Hazard ratio가 시간에 따라 모두 같아야한다는 것

기본적인 Cox 비례위험모형의 식은 다음과 같이 쓸 수 있다. 두 가지 변수 dose와 Sex 가 있다고 해보자. (각각 연속형, 범주형 변수이다.)

$$ h(t, X) = h_0(t)exp(\beta_1 dose + \beta_2 Sex) $$

예를 들어, 성별 (Sex)에 대해서 비례위험 가정이 만족하지 않는다고 해보자. 

$$ h_{1} (t, X) = h_{1}(t)exp(\beta_1 dose) $$

$$ h_{2} (t, X) = h_{2}(t)exp(\beta_1 dose) $$

그러면 위와 같이 Sex에 따라 따로 따로 모형을 fitting 시키는 것이 Stratified Cox regression의 기본 개념이라고 할 수 있다. 또한 두 식의 dose에 대한 beta 값이 같다는 것을 확인할 수 있다. Sex에 대해 층화를 한 후, 공통 beta 값을 추정하는 것이라고 할 수 있다. baseline hazard 인 h(t) 값만 두 식에서 다르게 추정한다. 이는 Sex의 hazard ratio를 추정하는 것이 아니라 비모수적 방법으로 각각의 Sex의 h(t)를 추정하는 것으로 이해할 수 있다. 

R에서의 구현 (survival 라이브러리)

coxph(Surv(time, event) ~ dose + strata(sex), data=data)

Stratified cox 모형의 한 가지 이슈는 층화한 변수가 다른 변수와 교호작용이 있는 경우이다. 이 경우에는 두식에서 공통된 beta를 추정하는 것이 아니라, 각각의 beta를 추정하는 것이 좋다. 이를 Interaction model이라고 한다. 

$$ h_{1} (t, X) = h_{1}(t)exp(\beta_{1} dose) $$

$$ h_{2} (t, X) = h_{2}(t)exp(\beta_{2} dose) $$

coxph(Surv(time, event) ~ dose, data=data[data$sex== 0,])
coxph(Surv(time, event) ~ dose, data=data[data$sex== 1,])

또는 아래와 같은 방식으로 표현할 수 있다.

$$ h_g (t, X) = h_g (t)exp(\beta_{1} dose + \beta_2(Sex*dose)) $$

위 식에서 (Sex=0 or 1) 이고, Sex=1일 때, beta2 값이 beta1 값과 더해져 dose의 beta가 된다. 이렇게 식을 쓰면 beta2를 이용하여 교호작용이 유의한지에 대한 검정을 할 수 있다. 교호작용을 검정하는 방법에는 두 가지 방법이 있다.

1. Wald test : 교호작용에 대한 beta값이 통계적으로 유의한지 검정

2. Likelihood ratio test : 두 모델의 likelihood ratio를 이용해 likelihood의 증가값이 통계적으로 유의한지 검정

결론적으로 Stratified cox regression 의 장점과 단점은 아래와 같이 요약해볼 수 있다.

장점 : PH 가정이 맞지 않는 변수에 대하여 층화 (Stratification)를 수행하여 분석하기 때문에, 다른 계수 추정을 더욱 신뢰성 있게 할 수 있다. 

단점 : 층화를 한 변수에 대해서는 HR을 추정할 수 없다. 층화한 변수와 다른 변수와의 교호작용을 고려하여야한다. 교호작용이 있는 경우 추정된 계수값은 바이어스가 생긴다. 

Competing Risk를 고려한 생존 분석

Competing risk 분석은 생존분석에서 관심 event 외의 다른 event가 있을 때, 분석하는 방법입니다. 예를 들어, 특정 약이 폐암으로 인한 사망을 줄이는가에 대해 관심을 갖고 연구하려고 하는데 심혈관 질환, 뇌혈관 질환 등으로 사망한 경우 Competing risk 분석을 통해 이 약이 폐암으로 인한 사망은 줄여도, 다른 사망에 대해서는 반대의 효과를 줄 수 있기 때문에 이를 확인할 필요가 있습니다.

Competing risk는 기본적인 생존분석 모형의 확장된 버전입니다. 데이터의 형태는 기본 생존 분석의 [time to event, event( censoring, event 를 나타내는 변수)]를 따르지만, vent가 censoring, event 0/1 로 나누어진 것이 아니라, 3개 이상의 범주형 변수를 갖게 됩니다. 

R의 경우 timereg, cmprsk 를 이용하여 competing risk 분석을 수행할 수 있습니다. 

먼저 데이터를 읽어옵니다. 

fol <- read.csv("data.csv", header=TRUE)

이 데이터에는 1,2 두 가지의 사망원인이 있습니다. 

cause
  0   1   2 
193 272  76 

1. 관심 사건 외 다른 이벤트들을 censoring으로 취급하여 분석

library(data.table)
fol <- setDT(fol)

fol$y1 = fol$cause
fol$y2 = fol$cause
fol[cause==2, y1:=0]
fol[cause==1, y2:=0]
fol[y2==2 , y2:=1]

원래는 0이 censoring 인데, 관심요인이 1일 때는 2를 censoring, 관심요인이 2일 때는 1을 censoring으로 바꾸어주는 코드입니다. 이 때, data.table으로 dataframe으로 변환한 후, := 키워드를 이용하였습니다.

① Cause = 1에 대한 Model (Cause=2 를 Censoring 취급)

=> age, hgb, clinstg를 보정해서 blktxcat의 effect 파악

library(survival)

cause1 <- Surv(fol$dftime, fol$y1)
cox1 <- coxph(cause1 ~ age+hgb+blktxcat+clinstg, data=fol)

Call:
coxph(formula = y0 ~ age + hgb + blktxcat + clinstg, data = fol)

  n= 541, number of events= 272 

             coef exp(coef) se(coef)     z Pr(>|z|)    
age      0.023784  1.024069 0.004733 5.025 5.02e-07 ***
hgb      0.003078  1.003083 0.004075 0.755 0.449974    
blktxcat 0.207383  1.230454 0.057329 3.617 0.000298 ***
clinstg  0.329831  1.390733 0.138171 2.387 0.016981 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
age          1.024     0.9765    1.0146     1.034
hgb          1.003     0.9969    0.9951     1.011
blktxcat     1.230     0.8127    1.0997     1.377
clinstg      1.391     0.7190    1.0608     1.823

Concordance= 0.623  (se = 0.019 )
Rsquare= 0.082   (max possible= 0.997 )
Likelihood ratio test= 46.23  on 4 df,   p=2e-09
Wald test            = 45.16  on 4 df,   p=4e-09
Score (logrank) test = 45.84  on 4 df,   p=3e-09

② Cause = 2에 대한 Model (Cause=1 를 Censoring 취급)

=> age, hgb, clinstg를 보정해서 blktxcat의 effect 파악

cause2 <- Surv(fol$dftime, fol$y2)

cox2 <- coxph(cause2 ~ age+hgb+blktxcat+clinstg, data=fol)

summary(cox2)

Call:
coxph(formula = cause2 ~ age + hgb + blktxcat + clinstg, data = fol)

  n= 541, number of events= 76 

             coef exp(coef) se(coef)     z Pr(>|z|)    
age      0.087125  1.091033 0.011440 7.616 2.62e-14 ***
hgb      0.005945  1.005962 0.008452 0.703    0.482    
blktxcat 0.164626  1.178952 0.109133 1.508    0.131    
clinstg  0.353410  1.423915 0.285801 1.237    0.216    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
age          1.091     0.9166    1.0668     1.116
hgb          1.006     0.9941    0.9894     1.023
blktxcat     1.179     0.8482    0.9519     1.460
clinstg      1.424     0.7023    0.8132     2.493

Concordance= 0.764  (se = 0.039 )
Rsquare= 0.131   (max possible= 0.759 )
Likelihood ratio test= 75.68  on 4 df,   p=1e-15
Wald test            = 62.84  on 4 df,   p=7e-13
Score (logrank) test = 67.68  on 4 df,   p=7e-14

2. Cumulative incidence curve를 통한 분석

Cumulative incidence curve의 정의 

개념 : Cumulative incidence curve란 competing risk를 고려한 1-생존함수 의 추정값 이다. 

$$ CIC_c(t) = P(T_c \lt t) $$

위 식에서 CIC란 어떤 사람의 c라는 event의 발생 시각인 Tc가 t보다 작을 확률입니다. 만약 competing risk가 없다면 이것은  간단하게 계산됩니다. 일반적으로 competring risk 가 없는 상황에서는 이것은 1 - (Kaplen-Meier 생존함수)로 추정하게 됩니다. 하지만 competing risk가 있는 경우, 다른 원인으로도 사람이 죽을 수 있기 때문에 Kaplen-Meier 방법을 그대로 사용하면 특정 시점 t 이전에 event 가 발생했을 확률이 과대추정됩니다. 

CIC 방법에서는 proportion of dying을 아래와 같은 방식으로 추정하게 됩니다. 

$$ CIC_c(t) = \int_{0}^{t} h_c(u)S(u) du $$

이 방식을 직관적으로 이해하면 Incidence를 모든 t에 대해서 더해준다고 이해할 수 있는데 incidence는 특정 시점에서의 cause-specific hazard와 생존함수를 곱해준 값입니다. 이 때, 생존함수는 인구집단에서의 생존함수입니다. 만약 이 생존 함수를 해당 cause에 대한 생존함수를 사용한다면 1-KM과 같은 값이 나오게 됩니다. 하지만 인구집단에서의 생존함수를 사용하여 competing risk를 고려하는 것으로 이해할 수 있습니다. 

Cumulative incidence curve는 cmprsk 라는 라이브러리를 통해 분석할 수 있습니다. Cumulative incidence는 cause-specific hazard function으로 계산되며, proportion of dying을 계산해줍니다. 이 때, 관심 요인 blktxcat을 지정해줍니다. blktxcat을 약의 종류라고 생각을하고 0,1,2,3 총 4개의 약의 종류가 있는 것입니다. 이 4개의 약이 1,2 두 가지 사망 원인에 어떻게 작용하는지를 아래 표를 통해 확인할 수 있습니다=> 각각의 strata 별로 proportion of dying의 예측값을 보여주고 있습니다. 

예를 들어, 25달이 지난 후, 사망원인 1을 보면, 0의 약을 쓴 사람은 52%가 죽은 반면, 3 약을 쓴 사람은 66%가 죽었습니다.  

library(cmprsk)
CIC <- cuminc(fol$dftime, fol$cause, fol$blktxcat)
CIC

Tests:
       stat           pv df
1 17.189363 0.0006461085  3
2  4.596607 0.2038333549  3
Estimates and Variances:
$est
             5         10         15        20        25        30
0 1 0.28456914 0.40301775 0.47654623 0.5235252 0.5235252 0.5235252
1 1 0.36620667 0.44808409 0.47719815 0.5096313 0.5541112        NA
2 1 0.42528724 0.54412907 0.59678622 0.5967862 0.5967862        NA
3 1 0.50604626 0.63733974 0.66670849 0.6667085 0.6667085        NA
0 2 0.02540834 0.08048914 0.15465458 0.1546546 0.1546546 0.2619280
1 2 0.11604642 0.16329519 0.18920323 0.2309031 0.2679696        NA
2 2 0.05245607 0.07941283 0.13288099 0.2003667 0.2409361        NA
3 2 0.03960396 0.06897272 0.09834147 0.1374998 0.2353957        NA

$var
               5           10          15          20          25          30
0 1 0.0010438691 0.0013684404 0.001665324 0.002070841 0.002070841 0.002070841
1 1 0.0022693095 0.0025757378 0.002723652 0.003459139 0.004840443          NA
2 1 0.0018165618 0.0021166396 0.002288543 0.002288543 0.002288543          NA
3 1 0.0025253539 0.0027149325 0.002701947 0.002701947 0.002701947          NA
0 2 0.0001266488 0.0004372677 0.001007219 0.001007219 0.001007219 0.006742183
1 2 0.0010075372 0.0014431850 0.001683726 0.002346285 0.003481490          NA
2 2 0.0003774854 0.0005948500 0.001226766 0.002385591 0.003610060          NA
3 2 0.0003819491 0.0007880136 0.001162941 0.001791598 0.011245328          NA

plot(CIC, xlab="month", ylab="proportion dying")

위와 같이 그림도 그릴 수 있습니다. 위쪽에 위치한 4개의 그래프가 사망원인 1에 대한 4가지 약의 효과 그래프이고, 아래쪽에 있는 4개의 그래프가 사망원인 2에 대한 것입니다. 전반적으로 사망원인 1로 죽는 사람이 많다는 것을 확인할 수 있습니다. 또한 약의 경우 사망원인 1에 대해 어느 정도 연관성이 있는 것으로 보입니다. 각 4개의 약에 대해 proportion of dying에 차이를 보여줍니다. 하지만 사망원인 2에 대해서는 약과 별로 상관 없는 것을 확인할 수 있습니다. 

3. Conditional probability function

library(Cprob)

CPC <- cpf(Hist(dftime, cause)~blktxcat, data=fol, failcode=1)
CPC
plot(CPC, xlab="month", ylab="proportion dying")

Conditional probability function도 마찬가지로 그릴 수 있습니다. 이 때, failcode로 관심 event를 지정해줍니다. 

4. Proportional-odds Model for the Conditional Probability Function

CPC_reg=cpfpo(Hist(dftime, cause)~blktxcat+age, data=fol, failcode=1, tis=seq(0, 70, 5))
CPC_reg

Warning message:

“glm.fit: algorithm did not converge”

Test for non-significant effect 

                  stat       pvalue

(Intercept) -26.769122 0.000000e+00

blktxcat      3.952040 7.748788e-05

age           5.041975 4.607511e-07

 Test for constant fit 

                   coef     exp.coef     SE.coef       stat     pvalue

(Intercept) -9.51827834 7.349609e-05 0.355569317 -24.964365 0.00000000

blktxcat     0.27674405 1.318829e+00 0.070025626   1.790873 0.07331363

age          0.02991248 1.030364e+00 0.005932691  10.309364 0.00000000

Conditional probability function의 proportional odds 모델을 이용하면 regression을 통해 각 factor들의 odds ratio를 구할 수 있습니다. 

data.csv

데이터 출처

Thomas H. Scheike, Mei-Jie Zhang

Title: Analyzing Competing Risk Data Using the R timereg Package

Journal of statistical software

https://www.jstatsoft.org/article/view/v038i02

h(t) = f(t)/S(t) 로 나타내어 지는데, 아래 식을 통해 확인해 보자.

생존분석 관련 함수의 정리

$$ S(t) = P(T> t) $$

$$ F(t) = 1-S(t) = P(T \leq t) $$

$$ f(t) = F'(t) $$

$$ h(t) = \lim_{\Delta{t}\to\ 0}P(t\leq T<t+\Delta{t} | T > t) = \frac{f(t)}{S(t)} $$

why? 직접 f(t)/S(t) 를 계산해보면 조건부 확률의 계산 공식을 통해 양변이 같다는 것을 알 수 있다.

$$ f(t) = \lim_{\Delta{t}\to\ 0}P(t<T<t+\Delta{t}) $$

$$ \frac{f(t)}{S(t)} =  \frac{\lim_{\Delta{t}\to\ 0}P(t<T<t+\Delta{t})}{P(T>t)} = \lim_{\Delta{t}\to\ 0}P(t\leq T<t+\Delta{t} | T > t) $$

$$ S(t) = e^{-\int_{0}^{t}h(u)du} $$

why? h(t) 라는 것은 1-S(t)를 미분한 것을 다시 S(t) 로 나눈 것이다. 이를 만족하는 S(t)는 위의 S(t) 밖에 없다. (if and only if 이다.)

카플란마이어 estimation

Censoring이 있는 데이터에서 생존함수를 추정하는 비모수적인 방법이다 만약 censoring이 아예 없다면, 생존함수는 그 시점에서 살아있는 사람을 보면 된다.  하지만 right-censoring이 있는 경우 해당시점에서 살아있는 사람은 censoring 된 사람을 제외한 사람일 것이고, 이 경우에 살아있는 사람만 계산하게 되면  생존 함수가 잘못 추정되게 된다.  따라서 censoring이 있을 때, 그 사람이 t시점까지 살았다는것을 활용하여 각 시점에서 survival rate을 구하여 계속 곱하면서 survival function을 추정한다.

http://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/BS704_Survival/BS704_Survival_print.html

case/control 디자인에서 샘플 수를 계산하는 방법 

샘플 수를 계산한다는 것의 의미  

실제로 노출 여부가 case/control 여부에 영향을 줄 때, (association이 있을 때) 해당 샘플 수에서 이를 충분히 detection 할 수 있는가?

이는 다음과 같은 귀무가설/대립가설 하에서 가설 검정할 수 있습니다. 

즉, Control에서 expose된 사람이랑 case에서 expose된 사람의 비율이 같은지를 검정하는 것이죠.

이를 검정하기 위해서는 p0, p1, OR을 알아야 합니다. 근데 p0, p1, OR 중에 2개를 알면 1개를 아래와 같은 식으로 구할 수 있습니다. 

따라서 일반적인 샘플 수 구하는 공식처럼, 유의수준, 검정력, p0, p1 을 주고 샘플 수를 구할 수가 있습니다. 한가지 더 주어야할 것은 case/control 디자인에서는 보통 case/control의 비를 1 이상으로 맞추기 때문에 r이라고 하는 case/control ratio를 주어야합니다.  

이 때, Kelsey의 공식에 따르면 필요한 case의 숫자는 아래와 같이 계산됩니다. 

Kelsey의 방법 

이 때, x = (p1+rp0) / (r+1) 라고 놓고, p=x(1-x) 로 계산됨 

참고문헌

Kelsey J.L., Whittemore A.S., Evans A.S.,and Thompson W.D. Methods in Observational Epidemiology. Oxford University Press, 1996. Print.

http://www.openepi.com/PDFDocs/SSCCDoc.pdf

벅슨의 역설 (Berkson's Paradox)

벅슨의 역설은 병원 기반 데이터로 case-control 연구를 할 때 발생하는 현상으로 1946년 Joseph Berkson이 보고한 역설입니다. 

벅슨의 역설이 생기는 이유

A와 B 둘다 있는 사람이나, A와 B 둘다 없는 사람이 연구 대상에서 제외되었을 때, A-B 간의 가짜 관계가 생깁니다. 벅슨의 역설을 병원 기반 데이터로 보여줄 때는 주로, A, B 둘다 없는 사람이 연구 대상에서 제외되었기 때문에 발생합니다. 

간단한 예

이 예를 통해 전혀 상관 없어보이는 두 가지 질병, 인플루엔자와 맹장염의 가짜 연관성(spurious association)이 어떻게 생기는지를 보입니다. (이미지 출처)

우선, 병원에 입원한 환자가 100명이라고 했을 때 다음과 같이 가정합니다.

1. 전체 인구집단에서 인플루엔자의 유병률은 10% 이다.

2. 병원에 온 사람 중 인플루엔자로 입원한 사람은 30%이다.

3. 병원에 온 사람 중 맹장염으로 입원한 사람은 10%이다. 

4. 인플루엔자와 맹장염은 독립적이다. 즉, 맹장염으로 입원한 사람중 10%가 인플루엔자를 가졌다고 가정. (맹장염의 유병률이 매우 낮다고 가정) 

(파란색=인플루엔자, 빨간색=맹장염, 검은색=기타, 빨간색+파란색=인플루엔자, 맹장염을 동시에 가진 사람)

이 상황에서 인플루엔자 따른 맹장염의 발생 양상을 파악하기 위해, 인플루엔자를 가진사람과 안 가진 사람으로 나뉘어 맹장염에 걸린 사람을 나눠봅니다. 총 30명의 인플루엔자를 가진 사람 중  1명이 맹장염에 걸렸습니다. 총 70명의 인플루엔자를 안 가진 사람중 9명이 맹장염에 걸렸습니다. 

위 노란박스는 70명의 인플루엔자를 안 가진 사람 중 9명이 맹장염에 걸린 상황을 나타냅니다.

이 때, 인플루엔자를 기준으로 한 맹장염의 상대위험도는 (1/30)/(9/70) = 0.004 로 인플루엔자가 맹장염에 보호적인 효과가 있는 것으로 나옵니다. 이 때, 이런 가짜 연관성이 생긴 이유는 병원 데이터의 특성 때문입니다. 즉, 인플루엔자와 맹장염이 둘 다 없는 사람이 일반 인구집단과 비교하여 더 적기 때문에, 인플루엔자가 없는 사람 중 맹장염에 걸린 사람이 9/70 = 12% 로 과추정되었습니다. 

참고

http://www.statisticshowto.com/berksons-paradox-definition/

로지스틱 회귀분석

로지스틱 회귀분석은 반응변수가 1 또는 0인 이진형 변수에서 쓰이는 회귀분석 방법입니다. 종속변수에 로짓변환을 실시하기 때문에 로지스틱 회귀분석이라고 불립니다. 로지스틱 회귀분석의 좋은점은 우선 계수가 Log Odds ratio가 되기 때문에 해석이 매우 편리하고, case-control과 같이 반응 변수에 따라 샘플링된 데이터에 대해서 편의(bias)가 없는 타당한 계수 추정치를 계산할 수 있다는 것입니다. 이 부분에 초점을 맞추어 로지스틱 회귀분석에 대해 간단히 정리해보려고합니다. 

1. 로지스틱 회귀분석의 원리

반응변수(y)가 1/0인 상황에서 마치 선형 회귀분석처럼 아래의 식에서 계수를 추정한다고 해봅시다. 

$$ y = \beta_0 + \beta X $$

로지스틱 회귀분석에서는 보통 우리가 예측하고자 하는 y의 "예측값"이 P[Y=1] 이라는 믿음을 갖는데, 그러면, y의 예측값은 0-1사이로 고정이 되어야합니다. 하지만 위 식의 경우 특정 X 값에서는 y의 예측값이 0-1 사이로 고정되는 것이 아니라 0보다 작아질 수도 있고, 1보다 커질 수도 있습니다. 

즉, 좌변의 범위는 0-1이며 우변의 범위는 마이너스 무한대에서 플러스 무한대의 값을 갖게됩니다. 이를 해결하기 위해 로지스틱 회귀분석에서는 반응 변수 대한 로짓 변환을 실시합니다. 로짓변환은 y를 log(y/1-y)로 만드는 함수적 변환을 말합니다.

우변을 0-1 사이로 변환하기

$$ P(y=1) = \frac {1} {(1+e^{-(\beta_0 + \beta X)})} $$

이 함수를 시그모이드 함수라고 부릅니다. y의 값이 0-1 사이로 고정됩니다. 이것이 왜 로짓변환이냐고 할 수 있지만 아래의 식을 보면 이해할 수 있습니다.

좌변을 마이너스 무한대에서 플러스 무한대로 변환하기

$$ \log(\frac{P(y=1)}{1-P(y=1)}) = \log(\frac{P(y=1)}{P(y=0)}) = \beta_0 + \beta X $$

이것이 y에 로짓변환을 한 식입니다. 위 두식은 모양만 다르지 결국 풀어서 보면 같은 식입니다. 보통 머신러닝 쪽에서는 위의 식을, 통계학에서는 아래식을 더 좋아하는 것 같아 보입니다. 위의 식이 더 이해하기 쉽지만, 아래 식은 해석이 좋은 특징이 있습니다.

왜냐하면 좌변이 X가 주어진 상황에서 로그 오즈(log odds)가 되기 때문입니다.  즉, 로지스틱 회귀분석은의 계수 추정은 독립변수 X에 대한 선형 회귀식을 X가 주어졌을 때의 반응 변수 Y의 로그 오즈에 적합시킨다." 라고 할 수 있습니다. 

2. 계수의 해석

로지스틱 회귀분석에서 계수에 exponential을 취하면, 해당 변수가 한 단위 증가했을 때, OR 이 됩니다.

$$ \log(\frac{P(y=1| x_1)}{P(y=0 | x_1)}) = \beta_0 + \beta_1 x_1 $$

$$ \log(\frac{P(y=1| x_1 + 1)}{P(y=0 | x_1 + 1)}) = \beta_0 + \beta_1 (x_1 +1) $$

이 때, 아래 식에서 위 식을 빼면, $$ \beta_1 $$ 만 남고, $$ \beta_1 = \log(\frac{P(y=1| x_1 + 1)}{P(y=0 | x_1 + 1)}) -\log(\frac{P(y=1| x_1)}{P(y=0 | x_1)}) $$ 가 됩니다. 

양변에 exponential을 취해주고, 로그함수의 성질에 의해 exponential된 계수는 X가 X+1이 됐을 때, 오즈의 비율 즉, 오즈 Ratio가 됩니다 (OR). 

3. Case-Control 연구에서 로지스틱 회귀분석의 계수 추정은 편의가 없다.

case-control 연구는 인구집단에서 case와 control을 샘플링해서 이를 대상으로 질병의 원인이 되거나 연관성이 있는 요인을 원인을 찾는 연구입니다. 예를 들어 전체 인구집단 중에 유병률이 0.001 밖에 되지 않는 질병이 있다면, 전체 인구 집단을 대상으로 연구를 수행하면 질병에 걸린 사람(case)를 연구에 포함하기가 매우힘들어지겠죠. 따라서 case-control study에서는 질병이 있는 사람들을 더 많이 샘플링을 해서 case와 control의 비율을 최대 1:1 까지 맞춥니다. 그렇다면 이 연구용 집단은 전체 모집단과는 특성이 다르게 됩니다. 따라서 많은 수치가 잘못된 계산값을 내놓게 됩니다. 예를 들어 상대위험도(Relative Risk)의 경우 참값과 다른 값이 나오게 됩니다.

이상적으로 case-control 연구는 case, control 여부에 따라서만 샘플링이 결정되어야 합니다. 로지스틱 회귀분석이 좋은 한 가지 이유는 이 가정이 맞은 경우, 편의가 없는 계수 추정치를 추정하게 됩니다. 편의가 없다는 말은 연구 집단에서 구한 계수의 기댓값이 전체 인구집단에서 구할 수 있는 계수의 "참값" 이라는 것입니다. 

왜 이렇게 되는 것인지 살펴보겠습니다. 샘플링된 데이터에서 X가 주어졌을 때, Y의 오즈는 아래와 같이 나타낼 수 있습니다. P의 아래첨자는 각각 population, sample을 나타냅니다. 이 때, 파이(y=1|x)는 case 샘플에 포함될 확률입니다. 파이(y=0|x)도 마찬가지입니다. 

$$ \frac{P_s(y=1| x)}{P_s(y=0 | x)} = \frac{\pi(y=1|x)}{\pi(y=0|x)} \frac{P_p(y=1| x)}{P_p(y=0| x)} $$

$$ \frac{P_s(y=1| x)}{P_s(y=0 | x)} = K \frac{P_p(y=1| x)}{P_p(y=0| x)} $$

$$ K \frac{P_p(y=1| x)}{P_p(y=0| x)} = K \exp(\beta_0 + \beta_1 x) = \exp(\beta_0' + \beta_1 x)$$

따라서, 샘플링 결과에 따라 계수의 추정량에는 영향이 없고 오직 절편만 바뀐다는 것을 알 수 있습니다.

Matching 데이터의 예

Matching이란 case의 control에서 혼란변수의 분포를 맞추어주는 방법으로 데이터셋을 구성할 때 자주 쓰이는 방식입니다. 이는 confounding을 방지하는 방법으로 알려져 있습니다. 이번엔 Matching 데이터를 분석할 때 매우 자주 쓰이는 Conditional Logistic Regression에 대해 아주 간단히 알아보겠습니다. R에서 먼저 사용할 데이터셋을 임포트합니다. 종속변수는 d이며, 이진변수입니다. 그래서 logistic regression으로 연관성 분석을 할 수 있는데요.

library(Epi)

data(bdendo)

데이터셋의 대한 설명은 help(bdendo)를 입력하면 나오고, 아래를 참고 바랍니다.

Format

This data frame contains the following columns:

set: Case-control set: a numeric vector

d: Case or control: a numeric vector (1=case, 0=control)

gall: Gall bladder disease: a factor with levels No Yes.

hyp: Hypertension: a factor with levels No Yes.

ob: Obesity: a factor with levels No Yes.

est: A factor with levels No Yes.

dur: Duration of conjugated oestrogen therapy: an ordered factor with levels 0 < 1 < 2 < 3 < 4.

non: Use of non oestrogen drugs: a factor with levels No Yes.

duration: Months of oestrogen therapy: a numeric vector.

age: A numeric vector.

cest: Conjugated oestrogen dose: an ordered factor with levels 0 < 1 < 2 < 3.

agegrp: A factor with levels 55-59 60-64 65-69 70-74 75-79 80-84

age3: a factor with levels <64 65-74 75+

하지만 Matching 데이터에 일반적인 Unconditional Logistic regression을 쓰면 bias가 생깁니다. 왜냐하면 Matching을 하면서 가져온 데이터에는 데이터 고유의 특성이 있기 때문입니다. 예를 들어서 time matching을 한 경우에, 비슷한 시기의 데이터를 샘플링을해서 하나의 strata를 만들고 데이터셋을 구성하게 됩니다. 그러면 이 시기에 의한 효과가 추정량에 영향을 주게 됩니다. 따라서 conditional logistic regression 이라는 조금 더 개선된 logistic regression 방법을 사용하여야합니다. (수학적인 설명은 생략하겠습니다..)

## Analysis

res.clogistic <- clogistic(d ~ cest + dur, strata = set, data = bdendo)

R로 conditional logistic regression(clr)을 하는 방법은 간단한 데 Epi 패키지의 clogistic을 활용하면 됩니다.  분석하고자 하는 공변량을 ~ 뒤에 넣고, strata에 matching pair를 나타내는 변수값을 입력합니다. 이 데이터의 경우, set 이라는 변수가 strata의 정보를 갖고 있습니다. 총 5개의 row가 같은 strata임을 알 수 있습니다. 

결과를 돌리면 다음과 같이 나오는 것을 확인할 수 있습니다.

res.clogistic

Call: 

clogistic(formula = d ~ cest + dur, strata = set, data = bdendo)

         coef exp(coef) se(coef)      z    p

cest.L  0.240     1.271    2.276  0.105 0.92

cest.Q  0.890     2.435    1.812  0.491 0.62

cest.C  0.113     1.120    0.891  0.127 0.90

dur.L   1.965     7.134    2.222  0.884 0.38

dur.Q  -0.716     0.489    1.858 -0.385 0.70

dur.C   0.136     1.146    1.168  0.117 0.91

dur^4      NA        NA    0.000     NA   NA

Likelihood ratio test=35.3  on 6 df, p=3.8e-06, n=254

association을 알아보기 위한 변수 cest, dur의 추정량을 볼 수 있습니다. 

쉽게 이해하는 민감도, 특이도, False Positive, False Negative, 양성예측도

민감도 (Sensitivity), 특이도(Specificity), 양성 예측도(Positive Predictive Value, PPV) 는 바로 무언가를 예측하는 상황에서 쓰이는 통계적인 지표입니다. 이 지표들이 가장 많이 쓰이는 곳이 바로 검진(Diagnostic)인데요.

우선 검진의 정의에 대해서 명확히 짚고 넘어가보겠습니다. 검진이란 "건강 상태와 질병의 유무를 알아보기 위하여 증상이나 상태를 살피는 일" 입니다. 이처럼 검진은 그 사람의 병의 유무를 알기 위해 일종의 검사를 하여 병의 유무를 예측하는 것입니다. 중요한 건, 검진은 실제와 다를 수 있다는 것입니다.  

예를 들어, 폐암 검진을 위해 흉부를 X-ray로 촬영한다고 해봅시다. X-ray는 저비용으로 폐암 여부를 알아보는 효과적인 방법입니다. 하지만 문제는 X-ray 결과로 그 사람의 폐암 여부가 결정이 나는 것이 아니라는 것입니다. X-ray에서 폐암이 의심이 가는 경우, MRI, CT 등의 더 정밀한 검사를 통해 폐암의 실제 여부를 결정하게 됩니다. 

이 때, X-ray가 "폐암 여부를 정확하게 가려내는 정도" 를 수치화할 수 있겠죠. 이렇게 수치화하는 지표 중에 유명한 것이 바로 민감도, 특이도, 양성 예측도입니다. 이를 이해하기 쉽게 설명하기 위해 다음과 같은 표로 많이 설명합니다.

폐암 X-ray 검진 예제

총 1000명이 X-ray 검사를 받았고, 이 중에 폐암환자는 10%인 100명, 비환자는 90% 입니다. 폐암환자 100명 중 90명이 X-ray 양성 판정을 받았고, 비환자 900명중 800명이 X-ray 음성 판정을 받았습니다. X-ray가 꽤 잘 맞춘 것이죠.

이 때, 폐암환자 중 양성 판정을 받은 사람의 비율을 민감도,

비환자 중 음성 판정을 받은 사람의 비율을 특이도라고 합니다. 

이 예에서, 민감도는 90/100 = 0.9, 특이도는 800/900 = 0.89 입니다. 이 때 1-0.9=0.1 = False Negative Rate, 1-0.89 = False Positive Rate 라고 합니다. 

즉, False Negative는 환자인데 검진에서 가려내지 못한 사람의 수, False Positive는 정상인인데 검진에서 환자로 판단한 사람의 수입니다. Rate가 붙으면 분모로 나누어 0과 1 사이의 값으로 나타내게 되는 것입니다. 민감도, 특이도, False Negative, False Positive는 이렇게 세트로 묶어서 이해할 수 있습니다.

그렇다면 양성 예측도, 음성예측도는 무엇일까요?

민감도와 특이도는 검진을 받은 "사람"의 관점에서 그 기기의 검진의 정확도를 판단한 것입니다. 반면, 양성 예측도, 혹은 음성 예측도의 경우 "기기의 관점"에서 검진의 정확도를 판단하게 됩니다. X-ray가 양성이라고 판단했을 때, 실제 폐암 환자일 확률은 90/190 = 0.47 입니다. 바로 이것을 양성 예측도라고 합니다. 반대로, X-ray가 음성이라고 판단했을 때, 진짜 정상인일 확률은 음성 예측도라고 하며, 800/810=0.98 입니다. 

문제는 양성 예측도가 왜이렇게 작게 나왔는가? 하는 것입니다. 양성 예측도는 민감도(0.9)와 특이도(0.89)처럼 그 값이 크지 않다는 것을 확인할 수 있습니다. 관찰력이 좋으신 분들은 직감으로 이걸 이해할 수 있으실 것입니다. 문제는 바로 "데이터의 불균형" 때문인데요. 비환자 900명중에 단지 100명만이 X-ray 양성으로 잘못 판단되었지만 (False Positive) 상대적으로 비환자가 1:9로 많기 때문에 그 값이 상대적으로 크게 잡혔다는 것을 볼 수 있습니다. 

이를 식으로 살펴보면, 양성 예측도 = 90 / 190 = 90 / 90+100 으로 풀어쓸 수 있겠죠. 이 때, 비환자가 환자와 마찬가지로 100명이었다면 90/90+100/9 = 약 0.9가 되었어야합니다. 하지만, 비환자가 환자에 비해 상대적으로 많기 때문에 False Positive가 많아져서 양성 예측도가 적게 나왔다는 것을 이해할 수 있습니다.

이를 두고, 양성 예측도는 유병률이 작을 수록 작다. 라고 표현합니다. 이 예제에서 폐암 유병률은 0.1이겠죠. 이 값이 작으면 작을 수록 False Positive에 의해 양성 예측도, PPV 는 작아지게 됩니다. 이를 보통 아래 식을 통해서 이해를 하게 되는데요. 

이 식은 위 테이블에서 변수를 놓고 방정식을 풀면 나오게 됩니다. (시간이 되시는 분들은 직접 손으로 풀어보시면 나옵니다.) 이 식을 보면, 분모는 검진의 양성 판단 횟수이며, 분자는 환자를 양성으로 판단한 횟수임을 알 수 있습니다. 

근데, sensitivity*prevalance는 분모, 분자에 같이 있으므로, (1-specificity)*(1-prevalence)가 실제 PPV 값에 영향을 주게 되겠죠. specificity가 고정되어있을 경우, prevalence, 유병률이 작아질 수록 분모가 커져 PPV 값은 작아진다는 것을 쉽게 알 수 있습니다. 

민감도, 특이도, 양성 예측도 중 어떤 지표로 검진 기기를 판단해야 할까요?

명확하게 무엇이 절대적인 기준이라는 것은 없습니다. 다만, 민감도의 경우 환자의 생명과 직접적인 관련이 있을 수 있습니다. 민감도가 낮은 검진의 경우, False Negative(환자인데 정상으로 판단)가 많아지게 되겠죠. 환자인데 정상으로 판단하여 적절한 치료를 받지 못하면, 환자의 생명에 지장이 있을 수도 있습니다. 이것이 많은 사람들이 민감도를 중요시하는 이유입니다. 반면, 특이도의 경우 불필요한 비용의 낭비와 관련이 있습니다. 특이도가 낮으면, False Positive(정상인데 환자로 판단)가 많아지게 되겠죠. 정상인이 환자로 판단되면, 환자가 불안감에 시달리거나, 불필요한 정밀한 추가검사를 받아 환자가 아니라는 것을 밝혀야합니다. 따라서 불필요한 손실이 일어나게 됩니다. 하지만 민감도와 특이도는 현실적인 관점에서 trade off의 관계에 있습니다. 민감도를 높이면 특이도는 보통 낮아지게 되죠. 따라서 적절한 수준에서 검진의 민감도와 특이도를 설정할 필요가 있습니다. 양성 예측도의 경우 유병률이 매우 낮은 질병에서는 높은 것을 기대하기가 힘듭니다. 위 식에서 봤듯이 유병률이 매우 낮으면, 데이터의 불균형으로 인해 False positive가 많아지고 양성 예측도는 낮아지게 됩니다. 이를 방지하기 위해서는 검진의 특이도를 매우 높여야하는데, 현실적으로 민감도를 포기할 수 없기 때문에 양성 예측도는 유병률이 낮은 질병에서 보통 낮은 경향이 있습니다. 

Nested case-control study

우선 nested case-control study란, 코호트 데이터에서 질병 상태에 따라 case-control을 뽑아내서 case-control study를 하는 것을 말한다. 예를 들어 어떤 코호트 스터디에서  10년간 추적관찰을 하였는데 10만명의 사람중 2천명이 유방암에 걸렸다고 하자. 이 때, case와 control을 적절하게 추출해 gene expression의 차이를 보고 싶다고 하자. 이것이 nested case-control study의 예라고 할 수 있다. 만약 cohort 초기에 gene expression 데이터셋을 구축하고자 한다면, 10만명의 사람들의 피를 추출해 gene expression DB를 구축해야할 것이다. 하지만 nested case-control study 방법으로 한다면, 보통, case의 1-4배 정도 되는 control을 추출해 gene expression을 분석하면 되기 때문에 시간과 비용을 아낄 수 있다. 실제로 nested case-control study는 biomarker를 이용하는 연구에서 많이 사용하는 방법이다. 

Retrospective cohort study

retrospective cohort study(후향적 코호트 연구)란 cohort study의 크게 두 가지 범주(전향적, 후향적 코호트 연구) 중 하나로, 연구 시점이 코호트 데이터를 모은 '후' 가 되는 연구를 말한다. 이는 nested case-control study와 비슷한데, 마찬가지로 연구 시점이 코호트 모집 전이 아니라 후이기 때문이다.  하지만 후향적 코호트 연구의 경우, 기본적으로 코호트 연구이며, case-control처럼 질병 상태에 따라 연구 집단을 뽑지 않는다는 차이가 있다. 예를 들어, 광산에서 일하는 광부들을 10년간 추적 관찰했다고 하자. 추후에 이 광부들을 대상으로 lung cancer으로 인한 사망률을 일반 인구집단과 비교해 알아보고 싶은 경우, 이것이 후향적 코호트 연구의 예가 될 수 있다. 이 광부 코호트로 nested cohort study를 하는 경우는 lung cancer에 따라 case-control을 추출 한 후, 이 case-control을 대상으로 exposure을 비교하는 경우이다.

참고 

https://en.wikipedia.org/wiki/Nested_case-control_study

http://epidemiologyclass.blogspot.kr/2004/10/difference-between-nested-case-control.html