Hard skills/Statistic (43)

Histogram Equalization 이란 무엇인가? 


IDEA



Histogram equalization (HE) 은 이미지를 전처리하는 방법이다. 위 그림만 보고도 HE 의 기본적인 아이디어를 알 수 있는데, HE 는 이미지의 contrast 가 적을 때 매우 유용하게 사용할 수 있는 방법이다. "Before Histogram Equalization" 이미지는 contrast 가 매우 떨어진다는 볼 수 있다. 오른쪽 histogram 을 보면 이미지 픽셀이 0~255에 고르게 퍼져있는 게 아니라, 일부분에 몰려있다는 것을 확인할 수 있다. "After Histogram Equalization" 을 보면, 이미지의 contrast 가 증가한 것을 확인할 수 있다. 이렇게 이미지의 밝기를 조정하는 것이 HE 로, HE 이후에는 밝기가 0~255 사이에 고르게 분포함을 알 수 있다. 이는 Cumulative density function (cdf) 을 통해 볼 수 있는데, 0~255 사이에서 cdf 가 일직선 형태를 나타내게 된다. 



Formula

L = 266 (일반적으로 0~255 사이에서 pixel intensity 를 표현하기 때문에)    


Example

 

HE 를 실제로 하는 방법은 해당 pixel intensity 의 cdf function 값을 구한 후, 0~255 사이의 uniform 분포에 맞게 값을 변환을 시켜서 넣어주면 된다. 예를 들어, 총 픽셀 수가 64 인 이미지가 있다고 하자 이 이미지에서 가장 낮은 pixel intensity 는 52 이다. 이를 어떻게 변환하는지를 보면 cdf(v) = 1, cdf_min = 1 MxN=64, cdf_min = 1, L-1 = 255 이고 분자가 0이므로 값은 0이 된다. pixel intensity 가 55 일 때, cdf(v) = 4, cdf_min = 1, MxN = 64 이므로, h(v) = 12 이다. cdf_min 은 cdf(v) 중 가장 작은 값으로 이 값을 분자 분모에서 빼주는 이유는 h(v) 를 0부터 시작하도록 하기 위함이다. 


v, Pixel Intensitycdf(v)h(v), Equalized v
5210
55412
58620
59932
601036
611453
621557
631765
641973
652285
662493
672597
6830117
6933130
7037146
7139154
7240158
7342166
7543170
7644174
7745178
7846182
7948190
8349194
8551202
8752206
8853210
9054215
9455219
10457227
10658231
10959235
11360239
12261243
12662247
14463251
15464255


참고

https://en.wikipedia.org/wiki/Histogram_equalization



Linearly weighted kappa 란 무엇인가?



Formula


$$ k_w = (P_0 - P_e) / (1 - P_e) $$



사용하는 이유


두 명의 측정자 (rater) 의 일치도를 평가 (inter-rater agreement) 할 때 쓰이며, 측정치가 순서가 있는 정성적 지표 (ordinal categorical variable) 일 때 쓰인다.  그냥 kappa coefficient 대신에 linearly weighted kappa 를 사용하는 이유는, 이러한 순서형 척도 변수에서 근접하게 측정된 값을 그냥 버리기 아깝기 때문에 이를 고려해서 일치도를 평가하기 위함이다. 


방법


P_0 =  proportion weighted observed agreement

P_e =  proportion weighted chance agreement


위의 두 값을 구하면 된다. 이를 위해 우선 아래와 같은 weight matrix 를 생성한다.  weight matrix는 측정치가 5개의 범주로 나누어지는 경우 아래와 같이 구할 수 있다. 



두 명의 측정자의 측정치가 아래와 같다고 한다. 예를 들어, 아래 테이블은 측정자 A와 측정자 B가 동시에 1이라고 판단한 경우는 8건이고 이는 전체의 0.07 이며, 측정자 A가 1, 측정자 B 가 2라고 판단한 경우는 2건으로 전체의 0.02 임을 의미한다.  





이 때, $$ P_0 = 0.07 * 1 + 0.02 * 0.75 + 0.009 * 0.5 + 0.026 * 0.75 + 0.09 * 1 ... + 0.02 * 1 = 0.93 $$ 로 구할 수 있다.


또한, $$ P_e =1 *(0.094 * 0.094) + 0.75 * (0.094 * 0.1709) + 0.50 * (0.094 * 0.530) ... + 0.02 * (0.0171*0.0171) = 0.78 $$ 로 구할 수 있다. 


결국 $$ k_w = (P_0 - P_e) / (1 - P_e) = (0.93 - 0.78) / (1 - 0.78) = 0.68 $$




Reference

http://www.anestesiarianimazione.com/DWLDocuments/The%20Linearly%20Weighted%20Kappa.pdf




로지스틱 회귀분석의 최대우도추정법 직관적 이해



Logistic regression의 Random component 



$$ Y \sim B(n, \pi) $$ 


Logistic regression은 반응변수 y가 위와 같은 이항 분포를 따른다고 가정한다.  



Model 


$$ \pi = \frac{1}{1+e^{-{(X\beta+\alpha)}}} $$ 


이 때, X를 통해 모수 pi 를 위와 같이 추정하는 것이다. 또는 일반적으로 아래와같이 표현되기도 한다. 


$$ logit(\pi) = X\beta + \alpha $$ 


Logistic regression에서 하고자하는건 결국, Random component에 주어진 pi를 추정하는 것인데, 이것이 0-1 사이의 값을 갖기 때문에 logit 변환을 수행하는 것이다. logit link function은 0~1사이의 값을 음의 무한대에서 양의 무한대로 바꾸어준다. 


Random component가 이항 분포이므로 likelihood는 아래와 같이 계산된다.



이곳에 Model을 집어넣음으로써, beta와 alpha 를 maximum likelihood 추정법으로 계산할 수 있게 된다. 근데 여기서 n과 y는 n번 시도했을 때 성공횟수 y를 나타낸다. 근데, 실제 각각의 데이터를 보면 시행횟수 n=1이며, y=0 or 1이다. 따라서 위 likelihood 는 베르누이 likelihood 와 같게된다. 



결국, 위의 likelihood를 최대화하면 계수 추정값을 얻을 수 있는데, log likelihood를 최대화하는 방식으로 쉽게 계수를 추정할 수 있다. 로지스틱 회귀분석에서는 보통 iteratively reweighted least squares (IRLS) 라는 방법을 통해 계수를 추정한다. 


$$ \sum_{i=1}^{N} y_i log(\pi_i) + (1-y_i) log (1-\pi_i) $$


Logistic regression의 log likelihood이다. 이는 negative binary cross entropy 이다. 결국, log likelihood를 최대화 하는 것은 negative binary cross entropy를 최대화 하는 것과 같고 이는 binary cross entropy를 최소화화는 것과 같다. 



고급통계 - Net reclassification improvements


Net reclassification index란 예측 모형에서 새로운 모델이 과거 모델보다 얼마나 얼마나 더 나아졌는지를 측정하는 지표이다. AUC 값을 기준으로 이를 판단하는게 과거에는 일반적이었으나, 그 한계점이 지적되외서 NRI가 등장하게 되었다. Framingham Heart Study에서 HDL 관련 연구를 수행할 때, AUC의 한계를 볼 수 있었는데, HDL을 예측 모형에 추가했을 때, AUC값이 크게 증가하지 않는 것이었다. 하지만 AUC 측면에서가 아니라, 모형의 계수를 통해 HDL의 effect를 파악할 때, HDL은 심장질환에 중요한 예측 인자였다. 


AUC의 한계는 다음과 같이 제시되었다. 1) 임상적인 의미와 연결시키기 힘들다. 2) 작은 변화는 잡아내기가 힘들다.


하지만 NRI는 이전 모델과 비교해 새로운 모델에서 얼마나 많은 사람들이 제대로 재분류되었는가를 파악함으로써 모델의 개선을 측정한다. 

NRI example table
EventTest 1Total, splitTotal
Non-eventAbnormalNormal
Test 2Abnormal1842228
246
Normal26872
85664
Total, split201030
106070
Total3070100

출처 - https://en.wikipedia.org/wiki/Net_reclassification_improvement


만약 환자 30명, 비환자 70명에 대해서 환자-비환자 여부를 판단하는 모델 Test1, Test2를 만들었다고 하자. Test 2가 Test1에 비해 얼마나 우월한지를 판단하고자 한다. 위 테이블에서 여기서 검은색 글씨는 Test 1, Test 2 둘 다 맞춘거고, 하얀색 글씨는 Test1, Test 2 둘 다 틀린 경우, 초록색 글씨는 Test 2에서만 맞추고, Test 1에서는 틀린 것, 빨간색 글씨는 Test 1에서 맞추고, Test 2 에서 틑린 것이다. 


이 때, NRI의 계산은 매우 간단한데, Test 2를 통해서 얻은 이득이 얼마나 더 많은지를 계산해주면 된다. 이것은 Test 2를 도입했을 때, 추가적으로 맞춘것-오히려 틀린것 을 해줌으로써 계산되며, Event, Non-event 두 그룹의 대상자에 대해서 계산한 다음에 더함으로써 최종 NRI를 계산하게 된다. 


즉, Event가 발생한 사람들에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_e = (4-2)/30 = 0.067 $$ 

다음으로 Non-event에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_{ne} = (8-4)/70 = 0.057 $$ 

따라서 NRI는 아래와 같이 계산된다.

$$ NRI = NRI_e + NRI_{ne}= 0.124 $$ 


생존분석 - 모수적 모형



Cox의 비례위험 모형


일반적으로 생존분석에 사용되는 모델인 Cox의 비례위험모형의 경우 준모수적 (Semi-parametric) 방법입니다. 이 방법은 생존함수 (Survival function)을 추정하기 위해 Hazard function에 대한 모델링을 합니다. 시간에 따라 변하는 Baseline hazard가 존재하고, Covariate 들의 effect에 따라 hazard가 달라진다고 모델링합니다. 그리고 hazard의 effect인 Hazard ratio가 시간에 따라 동일하다는 비례위험을 가정하여 Hazard function과 Survival function 을 추정하게 됩니다. 이 방법의 장점은 Survival function에 분포에 대해서 가정을하지 않기 때문에 강건하며, Covariate를 보정할 수 있다는 것입니다. 또한 Baseline hazard function이 어떤 함수인지를 알 필요가 없이 계수 추정을 통해 Hazard ratio만 구할 수 있다면 Covariate의 effect를 파악할 수 있습니다. 


모수적 모형


반면 모수적 모형은 Hazard function과 Survival function의 모양에 대해 Cox 보다 더 강한 '가정'을 합니다. 이 가정은 Hazard function과 Survival function이 어떤 분포를 따를 것이라는 믿음에서 출발합니다. 가장 많이 사용되는 모수적 모델은 Weibull 모형입니다. 


Weibull 모형의 Hazard function과 Survival function 은 다음과 같이 정의됩니다.


$$ h(t) = \lambda p t^{p-1} $$

$$ S(t) = exp(-\lambda t^p) $$ 


h(t)와 S(t)에 정의에 따라 h(t)가 저렇게 정의되면 S(t)가 저렇게 계산되는 것을 확인할 수 있습니다.


$$ \lambda = exp(\beta_0 + \beta_1 TRT) $$ 


만약 Treatment 변수에 대해 Weibull model을 만들 경우, 이 때 위와 같이 식을 놓고, beta0, beta1, p 3개의 계수를 추정하게 됩니다. 


또는 Hazard를 기준으로 계수를 추정하는 것이 아니라 Survival time을 기준으로 추정할 수 있는데, 이 경우가 모수적 방법에서 많이 사용되는 방법입니다. 이를 Acceleration Failure Time Model이라고 합니다. 


$$ S(t) = exp(-\lambda t^p) $$

$$ t = [-lnS(t)]^{1/p} \times \frac{1}{\lambda^{\frac{1}{p}}} $$

$$ \frac{1}{\lambda^{\frac{1}{p}}} = exp(\alpha_0 + \alpha_1 TRT) $$


AFT 모델에서는 위와 같이 놓고 alpha0, alpha1, p를 추정하게 됩니다. 


이렇게 두 가지 방법으로 구한 계수 beta와 alpha의 의미는 다른데, 만약 p가 고정되있다고 한다면, beta는 exponential을 취하면 Hazard ratio를 뜻하며, alpha는 exponential을 취하면 Accerleration factor가 됩니다. Acceleration factor에 대해서는 본 포스팅에서는 다루지 않을 예정이므로 자료를 참고 바랍니다. 


Strata 별로 p가 동일하다고 가정한 경우, Weibull model은 Proportional Hazard 가정과 Acceleration Failure Time 가정을 만족합니다. (https://en.wikipedia.org/wiki/Accelerated_failure_time_model)


PBC 데이터를 통한 Weibull model의 구축


pbc.dat

library(survival)

data <- read.table("pbc.dat")

# Column rename
colnames(data) <- c('id', 'futime', 'status', 'drug', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'edema', 'bili', 'chol', 'albumin', 'copper', 'alk_phos', 'sgot', 'trig', 
                   'platelet', 'protime', 'stage')

head(data)
data <- data[data$drug != '.', ]
data$drug <- as.character(data$drug)
data$drug[(data$drug == 2)] = 0
data$drug <- factor(data$drug)


Drug = 0 플라시보, Drug = 1 페니실린에 대해 Death에 대한 Time to event 분석을 수행합니다.


Weibull Model Fitting

library(survival)
library(eha)

# AFT  
swei = survreg(s~drug, dist = "weibull",data = data) 

# PH  
sweiph = phreg(s~drug, dist = "weibull",data = data)

summary(swei) 
summary(sweiph)
Call:
survreg(formula = s ~ drug, data = data, dist = "weibull")
              Value Std. Error     z      p
(Intercept)  8.3167     0.1094 75.99 <2e-16
drug1       -0.0438     0.1431 -0.31  0.759
Log(scale)  -0.1532     0.0729 -2.10  0.036

Scale= 0.858 

Weibull distribution
Loglik(model)= -1348.2   Loglik(intercept only)= -1348.3
	Chisq= 0.09 on 1 degrees of freedom, p= 0.76 
Number of Newton-Raphson Iterations: 5 
n= 312 
Call:
phreg(formula = s ~ drug, data = data, dist = "weibull")

Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
drug 
               0    0.491     0         1           (reference)
               1    0.509     0.051     1.052     0.167     0.759 

log(scale)                    8.317               0.109     0.000 
log(shape)                    0.153               0.073     0.036 

Events                    144 
Total time at risk        625985 
Max. log. likelihood      -1348.2 
LR test statistic         0.09 
Degrees of freedom        1 
Overall p-value           0.759387


Median Survival Time Comparison

predict(swei,newdata=data.frame(drug='0'),type=c("quantile"),se.fit=FALSE,p=0.5) predict(swei,newdata=data.frame(drug='1'),type=c("quantile"),se.fit=FALSE,p=0.5)

1: 2987.74956079141
1: 2859.63735172644



Weibull Model을 통해 구한 Survival Function과 Kaplen-Meier estimation의 비교

data$event_all=ifelse(data$status==0,0,1)
data$futime = ifelse(data$futime==0,0.1,data$futime)
s = Surv(data$futime,data$event_all)
sKM <- survfit(s ~ drug, data=data)
plot(sKM)
lines(predict(swei,newdata=list(drug='0'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')
lines(predict(swei,newdata=list(drug='1'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')

위 그림에서 검은색 직선이 Kaplen-Meier Survival function estimation이며, 붉은색과 파란색이 Weibull 모델로 추정한 Survival function 입니다. Weibull 모델의 S(t)의 모양에 대한 가정이 맞다면 두 직선은 가까울 것입니다. 그리고 위 직선들은 p가 같다는 가정하에 구해진 직선들입니다. 즉, PH, AFT 가정을 한 후에 그려진 직선입니다. 그렇기 때문에 Weibull 모형의 Survival function에 대한 가정, 그리고, PH, AFT 가정을 해도 될지에 대해서 체크를 해보아야합니다.



Weibull 가정의 체크 


(https://stat.ethz.ch/education/semesters/ss2011/seminar/contents/handout_9.pdf)


Weibull 모형에서 가정을 체크하는 방법은 X 축을 log(t), Y 축을 log(-log(S(t)) 로 그래프를 그려보는 것입니다. 이 때, S(t)는 Kaplen-Meier 방법으로 구한 Survival function 입니다. 이 때, 그룹간의 그래프의 모양을 보고 Weibull 생존모형 가정이 맞는지, 또한 두 그룹에서 PH, AFT 가정이 맞는지를 체크할 수 있습니다. 


먼저, 두 그룹의 직선이 평행하게 나오면 Weibull 가정을 만족합니다. 그렇지 않으면 Weibull 가정이 만족되지 못하며 AFT 가정도 만족되지 못합니다. 또한 두 그룹의 기울기의 비가 일정하면 PH 가정을 만족합니다. 


sKM <- survfit(s~ drug, data=data)
kmfit.TRT1 = summary(sKM)
kmfit.TRT2 = data.frame(kmfit.TRT1$strata, kmfit.TRT1$time, kmfit.TRT1$surv)
names(kmfit.TRT2) = c("TRT","time","survival")
TRT0=kmfit.TRT2[kmfit.TRT2$TRT=="drug=0", ]
TRT1=kmfit.TRT2[kmfit.TRT2$TRT=="drug=1", ]
plot(log(TRT0$time), log(-log(TRT0$survival)), pch=16, xlab="time in days using log-arithmic scale", ylab="log-log survival",type="p", xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2, col="Black",main="log-log curves by treatment")
par(new=T)
plot(log(TRT1$time), log(-log(TRT1$survival)), pch=3, col="grey50", ylab="",xlab="",  xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2)


이 그래프의 경우 우선 그룹에서 거의 직선이라는 것을 확인할 수 있지만, 두 그래프가 평행하지는 않습니다. (교차가 일어납니다.) 그렇기 때문에 Weibull 가정을 만족하지만 PH 가정과 AFT 가정을 만족하지 않으며, 두 그룹에서 p가 다르다는 알 수 있습니다.




Time-dependent variable의 정의 


Definition : 생존 분석에 있어서 시간에 따라 변수의 값이 바뀌는 변수를 time-dependent variable이라고한다. 


예 : 흡연 여부, 직장에서의 지위 등


Extended cox model : Time-dependent variable을 생존분석에서 사용하는 방법 


개념 : 만약 시간에 따라 변수의 값이 변화하는 변수가 있다면 extended cox model을 수행하여야한다. 


- 기존 데이터의 형태


ID
clinic
event
time
prison
dose
111428050
211275155
311262055
411183030
511259165
611714055

- Extended cox model을 위한 데이터의 형태


ID
clinic
prison
dose
ID2
start
time
event
110501070
1105017130
11050113170
11050117190
11050119260
11050126290

모든 time에 대해 데이터를 쪼개서 start와 time으로 정의한다. 하나의 row를 시간에 따라 쪼개므로 한 사람당 여러개의 row가 생기게 된다. 만약 해당 시간 내에서 event가 발생한 경우, event는 1이된다.


library(data.table)
library(survival)
data <- read.csv("addict.csv")

Extended cox model을 위한 데이터의 형태 만드는 법

data.cp <- survSplit(data, cut=data$time[data$event==1], end="time", event="event", start="start", id="ID2")

만약 아래와 같이 logdose라는 time-dependent variable이 있다고 하면 

data.cp$logdose <- data.cp$dose * log(data.cp$time)

Extended cox model 은 아래와 같이 수행할 수 있다.

Y = Surv(data.cp$start, data.cp$time, data.cp$event) coxph(Y ~ clinic + logdose +prison+ cluster(ID2), data=data.cp)


Time-dependent variable의 활용 1 : PH (proportional hazard) Assumption check


개념 : PH Assumption에 위반이 되는 것이 의심되는 변수가 있을 때, 이 변수에 대한 time-dependent variable을 정의함으로서 PH 가정을 체크할 수 있다. 


ID
clinic
event
time
prison
dose
111428050
211275155
311262055
411183030
511259165
6117140

55

만약에 위와 같은 데이터에 대해서 dose에 대해 PH Assumption 가정이 위반된다고 해보자. 즉, 이말은 시간에 따라 dose의 hazard ratio가 달라진다는 것이다. 이를 테스트하기위해 dose*g(t)라는 새로운 변수를 정의하고 이 변수의 계수가 유의한지 확인하는 방법을 통해 해당 변수에 대한 PH Assumption을 체크할 수가 있다. 


Y = Surv(data.cp$start, data.cp$time, data.cp$event)
coxph(Y ~ clinic + dose + logdose +prison+ cluster(ID2), data=data.cp)
Call:
coxph(formula = Y ~ clinic + dose + logdose + prison + cluster(ID2), 
    data = data.cp)

            coef exp(coef) se(coef) robust se     z       p
clinic  -1.01988   0.36064  0.21542   0.23637 -4.31 1.6e-05
dose    -0.08262   0.92070  0.03598   0.02960 -2.79  0.0053
logdose  0.00862   1.00865  0.00645   0.00525  1.64  0.1007
prison   0.34063   1.40584  0.16747   0.15972  2.13  0.0329

Likelihood ratio test=66.35  on 4 df, p=1e-13
n= 18708, number of events= 150 
   (1379 observations deleted due to missingness)


위 결과에서 logdose의 p-value > 0.05이므로, dose가 PH Assumption을 만족한다고 결론내릴 수 있다. 위 방법은 Wald test를 한 것이며, 다른 방법으로는 Likelihood ratio test를 통해 p-value를 구할 수도 있다. 


Likelihood ratio test

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
mod1 <- coxph(Y ~ clinic + dose + logdose +prison+ cluster(ID2), data=data.cp)
mod2 <- coxph(Y ~ clinic + dose +prison+ cluster(ID2), data=data.cp)

print(mod1$loglik[2])
print(mod2$loglik[2])

LRT <- (-2) * (mod2$loglik[2]-mod1$loglik[2])
pvalue <- 1 - pchisq(LRT, 2) 

이는 logdose를 포함한 full model과 logdose를 포함하지 않은 restricted model의 likelihood ratio를 구하고 이것이 카이제곱 분포를 따른다는 것을 이용한 검정이다. p-value를 계산한 결과 0.4 정도로 위의 wald test와 마찬가지로 logdose는 의 계수는 유의하지 않았다. 


따라서 dose에 대해 PH Assumption이 어긋나지 않는다고 결론 내릴 수 있다.



Time-dependent variable의 활용 2 : PH Assumption 만족 안하는 경우, 시간에 따라 달라지는 Hazard ratio 모델링


개념 : 어떤 변수에 대하여 heavy side function을 이용한 time-dependent variable을 정의함으로서 시간에 따라 달라지는 hazard ratio를 모델링할 수 있다. 



heavy side function은 어떤 값을 기준으로 그 이상일 때는 1, 미만일 때는 0을 돌려주는 함수를 말한다. 


예제 : 만약 clinic에 대해 시간에 따라 달라지는 HR을 모델링 하고 싶다고 해보자. 


fit <- survfit(Surv(data$time, data$event)~data$clinic)
plot(fit)

이것이 clinic=1과 clinic=2일 때 각각 그린 카플란 마이어 생존함수 추정값이다. 위 그림에서 hazard ratio가 시간에 따라 바뀐다는 것이 의심스러울 수 있다. 이 때, 어떤 시점을 기준으로 heavy side function을 이용한 time-dependent variable을 정의함으로써 이를 모델링할 수 있다. 


365를 기준으로 데이터를 쪼갠다. 

data.cp365 <- survSplit(data, cut=365, end="time", event="event", start="start", id="ID2")
ID
clinic
prison
dose
ID2
start
time
event
11050103650
1105013654281
21155202751
31055302621
41030401831
51165502591
data.cp365$hv1 <- data.cp365$clinic * (data.cp365$start < 365)
data.cp365$hv2 <- data.cp365$clinic * (data.cp365$start >= 365)

위와 같이 두 개의 heavy side function을 원래 clinic의 값에 곱해준 것을 두 개의 변수로 정의한다.  


ID
clinic
prison
dose
ID2
start
time
event
hv1
hv2
1105010365010
110501365428101
2115520275110
3105530262110
4103040183110
5116550259110
6105560365010
610556365714101

Extended cox model

Y365 <- Surv(data.cp365$start, data.cp365$time, data.cp365$event)
coxph(Y365 ~ dose + hv2 + hv1 + prison + cluster(ID2), data=data.cp365)


Call:
coxph(formula = Y365 ~ dose + hv2 + hv1 + prison + cluster(ID2), 
    data = data.cp365)

           coef exp(coef) se(coef) robust se     z       p
dose   -0.03548   0.96514  0.00643   0.00652 -5.44 5.3e-08
hv2    -1.83052   0.16033  0.38595   0.39838 -4.59 4.3e-06
hv1    -0.45937   0.63168  0.25529   0.25998 -1.77   0.077
prison  0.37795   1.45929  0.16841   0.16765  2.25   0.024

Likelihood ratio test=74.25  on 4 df, p=3e-15
n= 360, number of events= 150 


이 때 결과를 보면 hv1의 HR은 0.63, hv2의 HR은 0.16이다. 이 말은 clinic 1과 clinic2의 hazard ratio가 365미만에서는 0.63이었는데, 365 이상에서는 0.16이라고 해석할 수 있다. 즉, 두 clinic에서 보이는 생존률에 대한 경향성이 시간이 지날 수록 강해지고 있는 것이다. 이는 위의 카플란 마이서 생존곡선에서도 확인할 수 있다. 


  • 용코리 2021.04.07 16:09

    안녕하세요, 질문이 있습니다.

    제가 재발사건생존분석(recurrent survival analysis)를 시행 중에 있습니다.

    이 분석은 패널데이터로 구축 후 추정하게 되있습니다.

    그러나, 패널데이터로 구축시 독립변수들은 시간에 따라 변화하는 변수값으로 측정이 되어있기 마련입니다.

    1) 이와 같은 상황에서 time-varing covariate PH를 바로 추정하는건지,

    아니면 비례위험가정을 체크 후, 위반되야지 하는건지가 궁금합니다.

    또한,

    2) 재발사건분석에서 독립변수를 시간에 따라서 변화하는 변수(각 년도별 데이터)를 수집 안하고,

    한 년도의 데이터만 수집해도 괜찮은건지 궁금합니다. 종속변수의 경우 재발사건분석이기에 시간에 따라 변화하는 기간이 들어가구요.

    답변 부탁드리겠습니다.


Cox 비례위험 모형의 가정 체크


개념 : Cox 비례위험 모형은 "비례위험 가정" 을한다. 


비례위험가정은 시간에 따라 어떤 변수에서의 Hazard ratio가 일정해야한다는 것이다.

비례위험가정이 어긋나면, Stratified cox regression 또는 time-dependant variable 정의를 통한 extended cox regression을 수행할 수 있다. 


비례위험 가정의 체크는 아래방식으로 할 수 있다. 


1. Log-log plot

2. Observed-expected plot

3. Goodness of fit test


본 포스팅에서는 R의 survival, survminer 패키지를 통한 실습을 통해 cox 비례 위험 모형의 가정을 체크하는 방법을 알아본다. 


addict.csv


library(data.table)
library(survival)
data <- read.csv("addict.csv")

head(data)
ID
clinic
event
time
prison
dose
111428050
211275155
311262055
411183030
511259165

먼저 데이터의 형태는 위와 같이 생겼다. 


관심변수 : clinic 에 대하여 Kaplan-Meier 플롯 그리기

# Treatment 별로 그룹지어서 KM Estimation
fit = survfit(y~clinic, data=data)
par(mai=c(1,1,1,1))
plot(fit,main="KM Curve", xlab="Time(Week)", ylab=expression(paste(hat(S),"(t)")), col=c("blue", "red"))


우리는 이것이 비례위험 가정을 만족하는지를 확인하고 싶다. 비례위험 가정을 만족하려면 모든 timepoints에서 Hazard ratio가 같아야한다. hazard란 f(t)/S(t) (참고 : https://3months.tistory.com/349?category=743476) 로 위 그림만 봐서는 가정을 만족하는지 한 눈에 알 수 없다.  


로그로그 플롯


개념 : 로그로그 플롯을 그리고 그룹간에 평행하게 나오면 비례위험가정을 만족한다. 

plot(fit, fun="cloglog", lty=1:2, col=c("Black", "Grey50"), lwd=2, font=2, ylim=c(-3,2), font.lab=2, main="log-log KM curves by Rx", ylab="log-log survival", xlab="time in weeks using logarithmic scale")
legend("bottomright",lty=1:2,legend=c("Group1","Group2"),bty="n", text.font=2, lwd=2, col=c("Black", "Grey50"))


관찰값-기댓갓 플롯


개념 : Kaplen-Meier 생존 함수와 Cox regression을 통해 구한 생존함수를 동시에 그려본다.


이 때, Expected plot은 Cox 비례위험 모형에서 나오는 생존함수로 관찰한 생존함수와 모양의 차이가 많이나면 비례위험 가정을 만족하지 않는 것이다. 

kmfit <- survfit(y~clinic, data=data)
#observed
plot(kmfit, lty="dashed", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2, main="Observed Versus Expected Plots by Clinic", ylab="Survival probability", xlab="Survival time")
par(new=T)
#expected
exp <- coxph(y~clinic, data=data)
new_df <- data.frame(clinic=c(1,2))
kmfit.exp <- survfit(exp, newdata=new_df)
plot(kmfit.exp, lty="solid", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2)



Goodness of fit test


개념 : Goodness of fit test(GoF Test) 는 관찰값-기댓값의 차이를 통계적으로 검정하는 방법이다. (카이스퀘어 분포 이용)

# Gooodness of fit test
mod1 <- coxph(y~clinic, data=data)
cox.zph(mod1, transform=rank)

plot(cox.zph(mod1, transform=rank), se=F, var="clinic")
abline(h=0, lty=3)
plot(cox.zph(mod1, transform=rank), se=F, var="prison")
abline(h=0, lty=3)


        rho chisq       p
clinic -0.255  10.5 0.00118

p-value가 0.001이므로 clinic 변수는 비례위험 가정을 만족하지 않는다. 이는 그림을 통해서도 확인할 수 있다. 위 그림에서 직선이 평행하지를 확인하면되는데, 직선은 추정된 beta값을 의미하는데 직선이 평행하지 않고 어떠한 경향성이 보이면 beta 값이 시간에 따라 바뀐다는 의미로 비례위험 가정에 어긋난다. 




참고 : https://www.r-bloggers.com/cox-model-assumptions/

  • 2020.08.10 12:29

    비밀댓글입니다


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) $$


Interaction model은 R에서 다음과 같이 각각의 strata 별로 모델을 따로 만들어주면 된다. 
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)

age
path1
hgb
ldh
clinstg
blktxcat
relsite
ch
rt
survtime
stat
dftime
dfcens
resp
stnum
56NH04140NA21BY0.698152010.2381930181CR1
36NH02130NA21DY14.5023956112.4188911701CR2
39NH02140NA23YY4.914442210.0027378511NR3
37NH03140NA11Y15.6851472115.6851471591CR4
61NH04110NA22Y0.235455210.0027378511NR5
69NH02120NA11Y8.418891218.4188911701CR6

이 데이터에는 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

  • 2020.05.01 01:23

    비밀댓글입니다

  • 2020.08.28 10:43

    비밀댓글입니다


생존분석은 왜 하는가?

생존분석은 Censoring을 고려하여 Time to event에 대해 분석하기 위해 수행한다. 어떤 약의 효과를 판단한다고 해보자. 만약에 암에 걸린 사람을 대상으로 A라는 약을 처방하였는데, 이 약의 효과를 판단하기 위해서 여러가지 디자인을 해볼 수 있다. 

생존분석의 접근법은 약을 처방 받은 그룹 X, 처방받지 않은 그룹 Y에 대해서 생존 시간을 비교하고, 이 생존시간의 차이가 유의한지를 확인한다. 그럼 T-test로 할 수 있다고 생각할 수 있지만, 생존분석에서는 Censoring을 고려한다. 이것이 무슨말이냐하면 중도탈락한 데이터라도 그 데이터가 있었던 시점까지의 정보는 활용한다는 것이다. 생존분석은 중도절단된 자료의 부분적 정보를 최대한 이용한다.

즉, 기존 통계 모형을 쓰지 않고 생존분석을 쓰는 이유는 다음과 같이 두 가지로 정리할 수 있다.

1. Time to event를 알고 싶고,

2. Censoring 데이터를 고려하고 싶을 때


이는 linear regression, t-test, logistic regression 등의 다른 통계적 방법으로는 해결할 수 없다.

또한 생존분석의 목적은 다음과 같이 정리할 수 있다.

1. 어떤 사람의 Time to event를 예측 하고 싶을 때
2. 둘 이상의 그룹 간의 Time to event (생존예후) 를 비교하고 싶을 때
3. 변수들이 Event (생존)에 미치는 영향 파악 및 비교

Censoring이란 무엇인가?

Censoring은 생존분석에서 중요한 개념으로 두 종류로 나누어볼 수 있다.

1. left censoring : 관찰 기간보다 event가 발생 시각빠른 경우를 말한다.  

● 만약, 나무늘보가 언제 나무에서 내려오는지를 관찰하려고 한다. 근데 나무늘보는 항상 사람이 잠드는 새벽시간에만 땅으로 내려온다고 해보자 (예를 들어, 5am 정도). 그리고 사람은 아침 9시에 일어나서 나무늘보를 관찰한다. 그렇다면 사람은 평생 나무늘보다 땅으로 내려오는 것을 관찰할 수 없을지도 모른다. 이는 관찰기간 (9am~) 보다 event의 발생시각(5am) 이 전에 있기 때문이다. 
● 언제 담배를 처음 피었는가?
● 치매의 발생 (언제 정확히 발생했는지 알기 힘들다.)

2. right censoring :   관찰기간보다 event 발생 시각느린 경우를 말한다. 

● 가장 일반적인 censoring의 의미이다. censoring하면 보통 right censoring을 의미하는 경우가 많다. 
● 연구자는 어떤 연구 참여자를 평생 관찰할 수 없다. 예를 들어, 흡연자와 비흡연자의 폐암 발생까지의 시간을 비교하는 연구를 해본다고 해보자. 흡연자가 생각보다 건강해서 100살 넘도록 살수도 있을 것이다. 하지만 100년 넘도록 폐암의 발생여부를 주기적으로 관찰하는 것은 매우 어렵다. 따라서 적절한 시점에서 연구 종료를 선언하고, 지금까지 수집된 정보를 바탕으로 흡연자와 비흡연자를 비교하는 것이 적절할 것이다. 또는 중간에 폐암이 아닌 다른 원인으로 사망한다거나, 더 이상의 연구 참여를 하지 않겠다고 통보할 수도 있다. 이러한 경우를 right censoring 이라 한다. 
● 즉, right censoring이 발생하는 경우는 다음과 같은 상황을 예로 들어볼 수 있다. -> 스터디의 종료, 관심 event가 아닌 이유로 사망, 중도탈락
● 실제 분석 시, right censoring 이 있는 경우, 그 사람의 정보가 어떤 특정시점 t까지는 available 한 것을 이용하게 된다. 이러한 right censoring을 활용하는 것이 생존분석의 강점이라고 할 수 있다.

생존분석 관련 함수

생존분석 함수는 개인의 생존시간 T가 확률변수 (random variable) 라고 생각했을 때 이와 관련된 함수를 의미한다.

Survival function
S(t) = P(T > t) : 특정 시점 t에서 살아 있을 확률을 나타내는 함수이다. 즉, 이것은 event time T가 t보다 클 확률이다. 

F(t) : 특정시점 t까지 event가 발생했을 확률을 말한다. 이는 1-S(t)이다. f(t)의 cdf (cumulative density function)이다.
f(t) : 특정 시점 t에서 event가 발생할 확률을 나타내는 함수이다. (이것은 probability density function 이다.) 그리고 f(t)는 F(t)의 t에 대한 미분이다.

Hazard function
h(t) : t까지 살았을 때, 직후에 바로 event가 일어날 조건부 확률을 나타낸다.

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