Tensorflow GPU 버전을 설치



설치버전

Cuda 9.0v

Cudnn 7.0.5v

tensorflow-gpu 1.9.0


GPU 

1080TI


운영체제

Ubuntu 16.04


설치방법

아래포스트를 기본으로 오류를 수정해나가면서 설치함

https://medium.com/@zhanwenchen/install-cuda-and-cudnn-for-tensorflow-gpu-on-ubuntu-79306e4ac04e


설치 팁

- Cuda, Cudnn을 먼저 설치하고, (이 때, Cuda Cudnn의 버전이 tensorflow 에서 지원하는지를 확인해야 함. 최신버전의 경우 tensorflow가 지원을 안하는 경우가 있을 수도 있음.) 이에 맞는 tensorflow-gpu를 설치하는 것이 좋다.

- Cuda, Cudnn, tensorflow-gpu는 서로 의존적이므로, 버전에 특히 신경 써야하고, 함부로 업그레이드 하지 않아야하고 현재 버전을 기록해 두어야한다. 

- 굳이 하나하나 Tool 별로 Document를 찾아가면서 설치하는 것보다, 먼저 설치한 사람의 설치 방법을 보고 설치하는 것이 훨씬 빠르고 안전하다. 

- 위 포스트에서 cuda와 cudnn이 설치 되었다는 것을 확인하는 방법으로 sample 코드를 실행하는 방법을 사용하였는데, 확실한 설치 여부를 알 수 있으므로 좋다. 


설치할 때 발생한 오류

1. libstdc++6 의 문제 

해결 : https://askubuntu.com/questions/671791/lib32stdc6-package-depends-on-gcc-base-but-my-installed-version-is-newer


2. tensorflow-gpu 버전 문제

pip install tensorflow-gpu를 통해 설치하니 tensorflow-gpu 1.12.0 버전이 설치됨. 근데, 이 버전에서는 구버전의 cudnn을 쓰지 않는 문제가 발생. 따라서 tensorflow-gpu를 1.9.0 버전으로 다운그레이드 함. 


설치 완료 후, Multi GPU 환경에서 default로 0번 gpu만 사용하는 문제


tensorflow를 import하기 전에 아래 문구를 삽입함으로써, 특정 gpu를 사용하도록 강제할 수 있다.


import os

os.environ["CUDA_VISIBLE_DEVICES"]= "2"   



또한 CPU를 사용하도록 강제하는 방법은 tensorflow를 import하기 전 아래 문구를 삽입하면 된다.

import os

os.environ["CUDA_VISIBLE_DEVICES"]= ""  



tensorflow의 이전 방법에서는 아래 문구를 통해 특정 gpu를 사용할 수 있었는데 1.9.0 버전에서는 아래 문구가 안 먹힌다. 


import keras.backend.tensorflow_backend as K
with K.tf.device('/gpu:2'):


R igraph 설치 오류 해결


문제 : libgfortran.so.4 오류로 인해 igraph 패키지가 설치되지 않음


해결 : 

sudo apt-get install gcc-4.8-base libgfortran4


library(devtools)

install_github("igraph/rigraph")

R Default library path 바꾸기 (공용 라이브러리 경로 설정)


.libPaths() 를 하면 현재 사용하고 있는 라이브러리 경로들이 나온다.


[1] "C:/Users/사용자이름/Documents/R/win-library/3.3"

[2] "C:/Program Files/R/R-3.3.2/library"


R은 여기서 첫 번째 element를 Default로 쓰는데, 이것은 보통 개인 폴더로 잡혀 있다. 

하지만 다중 사용자 환경에서는 내가 설치한 라이브러리를 남이 사용하도록 해야할 경우도 있는데 이때는 아래와 같이 하면 된다. 


1. R_LIBS 환경변수를 공용 라이브러리 경로로 설정함 

2. R에서 Sys.getenv("R_LIBS")를 통해 환경변수가 잘 잡혔는지를 볼 수 있다.

3. .libPaths()를 다시 입력하면 공용 라이브러리 경로가 1번으로 잡혀있는 것을 확인할 수 있다. 

 

[1] "C:/Program Files/R/R-3.3.2/library" 

[2] "C:/Users/사용자이름/Documents/R/win-library/3.3"


이후에는 패키지를 설치하면 모든 사용자가 이용할 수 있게 된다. 


Python PDPbox 패키지 설치시 문제 해결


문제


 PDPbox 패키지를 설치하던 도중 아래 에러가 발생하면서 설치가 완료되지 않음


/usr/bin/ld: cannot find -lpython3.5m

 collect2: error: ld returned 1 exit status

 error: command 'gcc' failed with exit status 1


해결


sudo apt install python3.5-dev



구글링해도 안나와서 올립니다.


도움을 얻은 페이지

https://github.com/giampaolo/psutil/issues/1143

안녕하세요. 하이엔드 기계식 키보드 레오폴드 fc750r 텐키리스 갈축 구매를 하였는데 간단한 리뷰를 해보려고 합니다! 





구성품은 단촐합니다. USB 및 PS2 연결잭, 사용 설명서, 키보드 덮개 이 정도만 들어있네요! 처음 키보드를 꺼내보고 느낀점은 '생각보다 무겁다' 입니다. 타이핑을 해보니 키보드 자체에 무게감이 있다 보니까 굉장히 안정감이 느껴집니다. 어차피 키보드를 들고 움직일 일이 별로 상관없는 것 같습니다. 



레오폴드 fc750r의 특징은 PBT 이충사출 키캡을 사용하였다는 점과 체리 MX 스위치, 그리고 흡음패드를 사용하였다는 점입니다. 이중사출 키캡은 두 개의 키캡을 합쳐서 키캡을 만드는 방식입니다. 겉쪽 플라스틱에 글자를 파고 안쪽에 흰색 플라스틱을 덧대어 글자 부분을 매꾸는 방식으로 키캡을 만듭니다. 그렇기 때문에 오래 사용하더라도 글자가 지워질 염려가 없습니다. 



외관이 고급스럽고 이쁩니다. 확실히 레오폴드 명성에 맞게 잘 만들어진 키보드라는 느낌이 듭니다. 




위와 같이 키보드 높이 조절 하는 것이 있어서 개인의 선호에 따라 각도를 올려서 사용할 수도 있습니다. 



fc750r을 보고 느낀점은 기본에 충실하다는 느낌이 들었습니다. 화려한 불빛 등으로 치장된 키보드들도 많은데 레오폴드 fc750r의 경우, 별다른 치장을 하지 않은 클래식한 멋이 있는 기계식 키보드라는 생각이드네요. 또 키보드 자체가 굉장히 단단하고 고급스러운 느낌이 납니다. 화려한 걸 별로 좋아하지 않는 분, 또는 튀고 싶지 않은 직장인 분들이 사용하면 좋을 키보드라는 생각이 듭니다!

 

또한 흡음패드를 사용하였다고 하는데, 기계식 키보드를 처음 사용하는 입장에서 타이핑할 때 나는 소리가 그렇게 작다는 생각은 들지 않습니다. 저소음을 원하는 분들은 구매할 때 고려하시면 좋을 것 같네요. ^^ 타건 영상 위에 올려놓았습니다. 


이상 레오폴드 fd750r 텐키리스 갈축 리뷰였습니다.



생존분석 - 모수적 모형



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가 다르다는 알 수 있습니다.




의생명 분야에서의 신경망 모델 (Deep learning in biomedicine)

본 포스팅에서는 Nation biotechnology에서 publish된 논문을 기반으로 하여 의생명 분야에서 적용되고 있는 신경망 모델 (Deep learning) 에 대해서 살펴보려고합니다. (https://www.nature.com/articles/nbt.4233)

자율 주행차, 게임, 음성 인식, 텍스트 인식 등 Deep learning은 인공지능 연구자들과 일반 대중들의 마음을 사로잡고 있습니다. 최근에는 Deep learning은 임상 의사들의 관심도 끌고 있습니다. 지금까지의 많은 분야에서의 AI를 통한 모델링의 목표는 인간 레벨의 인공지능 (human-level AI)였습니다. 이미 사람이 잘 하는 것들을 높은 정확도로 모방하는 것입니다. 

이러한 분야들과 의생명 분야 (Biomedical) 데이터의 차이점은 의생명 분야의 데이터는 사람도 이해하기 힘들다는 것입니다. 예를 들어, genome 데이터는 사람이 맨눈으로 보고 바로 이해하기 어렵습니다. Biomedical 분야에 AI를 적용하는 것의 한가지 목표는 이렇게 사람이 쉽게 하기 힘든 분야에서 AI가 추론을 통해 사람에게 통찰력을 전해줄 수 있을까 하는 것입니다. 

딥러닝과 Biomedical 분야의 간략한 역사


실제로 이러한 것들을 실현 가능하도록하는 기회가 만들어지고 있습니다. 큰 규모의 Biomedical dataset들이 수집되고 있기 때문입니다. 예를 들어, NGS 데이터를 들 수 있습니다. Biomedical 분야에 Deep learning을 적용한다는 것은, 다양한 분야에서 각광 받고 있는 AI의 방법론인 Deep learning 라는 도구를 이용해서 이 데이터들로부터 유용한 가치를 이끌어내고, 과학적인 발견을 하고자 하는 것입니다. 


NGS (Next-generation sequencing) 데이터의 증가




머신러닝과 Deep learning의 핵심 아이디어를 간단하게 이야기해보겠습니다. 간단한 모델, 예를 들어, Linear regression과 다르게 Deep learning은 조금 더 복잡하고 유연한 모델링을 할 수 있습니다. Label을 Input feature로 바로 mapping 하는 것이 아니라 intermediate variable을 만들고, 이 intermediate variable의 function으로 output을 예측하는 모델을 만드는 것입니다. Deep learning의 강점은 이러한 Feature (독립변수)와 Label(종속변수) 사이의 어떠한 복잡한 함수관계도 유연하게 모델링할 수 있다는 것이며, 이론적으로도 이것이 가능하다는 것입니다 (Universal approximation theorem). 1980-90년대쯤 이러한 복잡한 모델링이 대한 이론적이 이미 정립이 되어있었지만, 이를 실현 가능한 하드웨어 기술이 부족해 실현할 수 없었습니다. 현 시대에서는 이러한 복잡한 모델에서 계수를 찾는 것이 계산적으로 가능합니다. 이는 이러한 Deep learning model의 수많은 계수 추정을 효율적으로 할 수 있는 다양한 알고리즘 및 라이브러리들 (Pytorch, tensorflow) 의 등장과 GPU 등의 발전으로 인한 컴퓨팅 파워의 발전의 결과입니다.  


Shallow model vs. Deep model.
 


우선 Shallow model과 Deep model의 차이에 대해 간단하게 설명을 하고 넘어가려고 합니다. Linear regression의 경우, Input feature와 Label 사이의 Linear 한 관계를 가정합니다. 그래서 예측값은 Input feature와 Weight의 linear combination으로 만들어 집니다. (여기에 bias를 더해준 값이 됩니다.) 

하지만 딥러닝의 한 종류인 Multi-layer perceptron의 경우, Layer라는 개념을 도입해서 output을 바로 예측하지 않고, Intermediate variable들을 만듭니다. Deep learning 이란 이러한 Layer 구조를 도입하여 Input feature와 Label 사이의 복잡한 함수관계를 모델링 하는 것을 말합니다. 이러한 구조를 도입하면, X와 Y의 관계가 linear한 관계가 아닌 경우, 혹은 다양한 형태의 Interaction이 존재하는 경우에, Y를 Linear regression 보다 더 잘 예측할 수 있게 됩니다. 즉, Deep model은 bias를 줄인다고 할 수 있습니다. 반면, Overfitting 등의 문제로 variance가 커질 수 있다는 것이 Deep model의 단점이기도 합니다. 


Genome sequencing 데이터에 대한 신경망 모델의 응용


다음으로는 Biomedical 데이터의 한 종류라고 할 수 있는 Genome sequencing 데이터에 어떻게 Deep learning이 어떻게 적용되는지를 알아보겠습니다. 예를 들어, motif detection (transcription factor binding site detection) 같은 분야를 예로 들어보면, 이 분야에는 Bioinformatics 분야에서 전통적으로 자주 사용되었던 Position weight matrix라는 방법이 있습니다. 어떤 문제에 대해 이러한 DNA Sequence들이 알려져 있을 때, 이것에 기반하여 아래와 같은 matrix를 만듭니다. 

Position weight matrix (position probability matrix)


이 Matrix는 해당 위치에서 어떤 sequence가 발견된 확률을 나타내는 matrix입니다. 이것을 기반으로 새로운 sequence가 들어왔을 때, 어떤 score를 내주고 이 score를 기반으로 sequence에 존재하는 어떤 pattern을 detection 할 수 있습니다. 이것은 Sequence와 어떤 pattern을 직접적으로 mapping 시킨 것으로 볼 수 있습니다. 반면 Deep learning에서의 방법은 layer를 더 만들어, Sequence와 어떤 pattern 사이에 존재할 수 있는 복잡한 관계를 모델링할 수 있습니다. 


이 그림은 CAACTT 라고 하는 sequence pattern을 찾는 CNN Model을 나타낸 것입니다. 물론 실제 상황에서는 CAACTT라고 하는 Pattern을 우리가 알 수는 없습니다. Deep learning에서는 수많은 데이터를 주고, 그 속에서 CAACTT라고 하는 Pattern을 딥러닝이 직접 학습하도록 하는 것입니다. 그리고 딥러닝의 강점이 바로 이런 것입니다. CNN에 익숙하신 분들은 잘 아시겠지만, 이 것은 1-D Convolution에 Same padding을 적용한 것으로 볼 수 있습니다. 그림에 나와있듯, Filter size는 3으로 잡고, 총 Input sequence로는 18의 길이를 갖는 sequence를 넣어준 것입니다. 이것은 일렬로 죽 늘어선 1차원 공간상의 이미지로 볼 수 있습니다. 2차원 이미지에는 W-H-(RGB) 3개의 차원이 있다면, Genome data의 경우 W-(ATCG) 2개의 차원이 있는 것입니다. 

딥러닝을 통해 Sequencing 데이터에 무슨 일을 할 수 있는가?

다음으로는 이러한 Deep learning이 Biomedical 분야에 어떻게 적용되고 있는지를 설명하였습니다. genetic data에 어떻게 되고 있는지를 먼저 예로 들었습니다. 현재, Genotype-Phenotype 관계를 규명하는데 GWAS (Genome-wise association study) 라고 하는 도구가 사용됩니다. 수많은 샘플을 보아 variant랑 phenotype의 association을 통계적으로 보는 것입니다. GWAS를 통해 variant를 찾는 것에 추가적으로 variant의 function을 연구하는 것도 한 가지 주제입니다. 왜냐하면, GWAS를 통해 찾아낸 variant (SNPs)는 그것이 질병과 인과관계를 갖는다고 보기 어렵기 때문입니다. GWAS의 경우 Mendelian disease와는 다르게 LD, 작은 effect size, regulatory network의 복잡한 구조 등으로 인해 causal variant를 찾기가 힘듭니다. 몇몇 coding 지역에 위치한 causal variant는 코돈을 통해서 그 변이의 effect를 예측할 수 있지만, non-coding variant 같은 경우, 해석이 매우 어렵습니다. Deep learning은 이 분야에서 적합합니다. 바로 transcription, splicing, regulation 등의 Molecular phenotype과 genetic variant의 관계를 보는 것입니다. 

신경망을 통한 Molecular phenotype 예측 소프트웨어, 논문

1. SPIDEX
    · DNA sequence → percent-spliced-in of cassette exons across tissues
2. DeepBind
    · DNA and RNA sequence → transcription factor and RNA-binding protein binding
3. Basset
    · DNA sequence → DNase hypersensitivity
4. DeepSEA and DanQ
    · DNA sequence → transcription factor binding
5. TITER
    · DNA sequence → translation initiation sites

이것들이 현재 논문으로 나와 있거나 소프트웨어로서 구현된 variant를 통해 molecular phenotype을 예측하는 구현체들입니다. 이러한 non-coding variant로부터 molecular phenotype을 예측하는 일에는 주로 CHIP-seq이나 DNase-seq 데이터를 트레이닝 데이터로 이용하며, DNA sequence로부터, transcription factor binding이나, DNase hypersensitivity (이를 chromatin feature 라고도 합니다.) 같은 것들을 예측합니다. 


DeepSEA (Deep learning based sequence analyzer)

DeepSEA는 Deep learning based sequence analyzer의 약자인데, Genome Sequence를 input으로 받아들여, variant의 effect를 chromatin에 미치는 영향을 바탕으로 예측합니다. ENCODE, Roadmap epigenomics 와 같은 genome의 function을 찾아내려는 목적을 갖고 수집된 데이터들, Chip-seq 또는 DNase-seq 데이터를 학습해, variant가 chromatin feature 에 미치는 영향을 학습한 후, 최종적으로 학습된 모델에 wild-type과 variant가 있는 두 개의 input을 주고, output을 뽑아내서 이 output의 비율로 그 variant의 effect를 예측하는 방식입니다. 


DeepSEA의 데이터 구축 (Ground truth) 및 딥러닝 모델 구조

DeepSEA의 데이터 구축 과정 및 모델 아키텍쳐입니다. 데이터 구축 과정에서는 총 919개의 chromatin feature를 학습하도록 Ground truth가 마련되습니다. 시퀀스를 통해 이것이 chromatin feature (예를 들어, Dnase hypersensitive site) 인가? 하는 것이 바로  Y, 즉 예측하고자 하는 것입니다. X로는 whole genome을 200bp 로 나눈 후, 이 시퀀스 중 절반 이상이 919개의 chromatin feature의 peak region에 포함되면 1 아니면 0으로 코딩되었습니다. 모델로는 Convolutional neural network를 사용하였고, Regularization을 위해 L2, L1 regularization, Dropout을 사용하였습니다. 



모델 구축


DeepSEA의 목적은 총 2개로 나눌 수 있습니다.


1. 해당 Variant가 Chromatin feature에 미치는 영향 파악, 이를 통해 variant의 각각의 chromatin feature (919개) 에 대한 기능을 파악할 수 있습니다.

2. 해당 Variant의 overall한 functional prediction 


1을 위해서 variant가 없는 sequence (1000bp 단위) 와 variant가 있는 sequence를 모델에 넣어서 log(P(reference)/P(alternative))를 통해 해당 chromatin에 variant가 미치는 영향을 파악합니다. 2를 위해서는 총919개의 chromatin feature에 대한 예측값과 함께 Evolutionary conservation score를 이용합니다. 저는 이 부분이 DeepSEA 가 대단한 부분이라고 생각합니다. 기존에 알려진 과학적 지식을 Deep learning 모델에 통합하여 더욱 잘 functional score를 예측하는 것이죠. 일종의 앙상블 모델이라고 볼 수 있습니다. 이를 통해 딥러닝 모델로부터 발생할 수 있는 오버피팅을 방지하고 더욱 robust한 모델이 될 수 있습니다. 



모델 테스트


위 링크에 방문하면 실제로 웹으로 구현된 DeepSEA 를 이용해볼 수 있습니다. variant를 나타내는 VCF file을 인풋으로 넣어주면, VCF file에 있는 각각의 variant에 대해, 919개의 chromatin feature에 대한 영향과, overall한 functional prediction score를 구합니다. 

신경망을 통한 Phenotype 을 예측하는 것이 가능한가?

Sequence를 통해 바로 phenotype을 예측할 수 있지 않겠냐 라는 의문이 들 수 있습니다. 물론, 최종적으로 나아가야할 목표는 sequence를 통해 phenotype 보는 것입니다. 하지만 아직까지 딥 러닝이 그 정도의 수준은 아닌듯 합니다. Molecular phenotype은, genotype에서 phenotype으로 이어지는 복잡한 메커니즘 중 최하단에 있다고 말할 수 있는 것이고, 아직까지 그것조차 제대로 해결되어지고 있지 않습니다. 그 다음으로 해결할 과제는 network-level gene interaction, physiological process 등 여러가지 해결해야할 과제들이 남아있습니다. 하지만 본 논문에서는 궁극적인 목표는 딥러닝에 genotype-phenotype 데이터를 학습시키고, 여기에 여러가지 생물학적 지식, 실험 데이터를 합쳐서 바로 phenotype을 예측하는 것이라고 주장하고 있습니다. 저는 DeepSEA에서의 사례처럼 Evolutionary conservation score 처럼, 학습된 딥러닝 모델에 이러한 explicit한 "지식" 들을 어떻게 통합하냐가 하나의 해결해야할 과제로 보입니다.


Medical Image에서의 신경망의 응용

다음으로는 medical image 분야에서의 딥러닝의 응용입니다. 이 분야는 딥러닝이 가장 직접적으로 응용이되고 있는 분야이고, 실제로 임상에서 활용이 되고 있는 분야이기도 합니다. Medical image의 특징은, Multi-modal, 즉 MRI, X-ray, CT 등 다양한 방법으로 이미지가 얻어지고, 같은 MRI 데이터라도 세팅값, 기기 종류에 따라 intensity가 모두 다릅니다. 또한, CT 같은 경우 3D 이미지 이기 때문에 다루기 까다롭습니다. 하지만 그럼에도 불구하고 딥러닝이 가장 성공을 한 분야이기도 합니다.

신경망을 Medical Image에 적용할 때 발생하는 이슈들

1. 이미지에 대한 높은 수준의 해석은 Automation이 힘듭니다. 이는 사람에 있어서도 사람마다 주관적인 기준이 있기 때문에 Intra-class variation이 크기 때문입니다. 

2. 또한 딥러닝 자체가 블랙박스의 성격을 띈다는 것입니다. 특히, 이미지 진단의 경우에는 이해 관계가 매우 크기 때문에, 딥 러닝의 결과를 어떻게 설명하느냐가 매우 중요합니다. 이 논문에서는 딥러닝이 그렇게 예측한 이유, 그 지역을 highlighting을 해주는 등의 전략이 필요하다고 언급하였습니다.

3. 세 번째는 사람의 예측값을 Ground truth로 해서 트레이닝 했을 때, 사람의 성능을 뛰어넘기 힘들다는 것입니다. 물론, 딥러닝은 계속 같은 값을 내기 때문에 Reliability는 좋습니다. 다만, 그 트레이닝 데이터가 한 사람의 기준에만 맞춘 경우, 다른 데이터를 대상으로 했을 때, bias가 생길 수 있다는 것입니다. 그러므로 Multi-expert consensus가 중요하다고 할 수 있습니다. 

국내 Medical Image 응용

국내에서는 2018년도에 Vuno에서 최초로 식약처에서 의료기기 사용허가를 받았습니다. 성장기 소아에서 X-ray 영상에서 인공지능을 통해 골연령을 자동으로 측정해주는 Vunomed-Boneage 소프트 웨어를 통해서입니다. Vuno에서는 이 방법의 효율성을 보이기 위해 임상시험 논문을 냈습니다. 결과에 따르면 AI가 2명의 의사들보다 consensus와의 concordance가 더 좋았고, 특히 의사가 AI 를 보조적으로 활용했을 때, 그렇지 않았을 때보다 concordance가 증가하는 결과를 보였습니다. (이 부분에 대해 궁금하신 분들은 최윤섭님의 유투브를 참고하시기 바랍니다. https://www.youtube.com/watch?v=wqXzmChH3N0&t=349s)




딥러닝과 실제 세계의 괴리

다음으로는 Deep learning이 실제 현실에 사용되었을 때 발생할 수 있는 문제점에 대한 것입니다. 

1. 가장 중요한 것은 반드시 Deep learning의 성능이 보장 되어야한다는 것입니다. 
    • · 이를 위해서는 Performance 측정을 해야합니다. C.V 나 hold out validation 같은 방법을 통해 충분한 validation이 이루어져야 합니다.
    • · 두 번째는 Deep learning을 Overfitting이 큰 문제이기 때문에, 모델의 불확실성이 어느정도인지를 보여주어야 한다는 것입니다. 통계적인 모델처럼 딥러닝은 결과의 신뢰구간을 통계적인 방법으로 얻어낼 수 없기 때문에 부트스트랩이나 베이지안 방법 등을 통해 모델의 신뢰구간을 추정할 수 있습니다. 

2. 또한 딥러닝의 문제점은 딥러닝이 목표하는 바와 실제 목표하는 바가 다를 수 있다는 것입니다. 이것이 Target mismatchloss function mismatch 입니다. 예를 들어, 임상에서는 종양의 크기가 일정한 임계치를 넘는 것이 중요한데, 딥러닝 모델을 트레이닝 할 때는 Intersect over union을 최소화 하기 위해 트레이닝을 보통 이용하곤 합니다. 이 경우에 딥 러닝 모델은 실제로 원하는 결과에서 bias가 생길 수 있습니다.

3. 다음으로는 오로지 딥러닝은 현재 이용가능한 데이터에 기반하여 예측 모형을 만드는 것이기 때문에 selection biasconfounding 이 생길 수 있고, 이 경우에 Causality 는 추론하는것이 매우 힘들다는 것입니다. 예를 들어, 여자라는 단어는 인문학, 남자라는 단어는 이공계와 연관시키거나, 흑인을 백인보다 더 위험하다고 학습할 수 있습니다. 


신경망 모델이 신뢰를 얻기 위해서는?

1. Performance를 보장해야합니다.
    • · Stakeholder가 원하는 메트릭을 제공해야 한다.
    • · Performance를 보장하기 위해서는 데이터 가공, 모델 선택, 오버피팅, 아웃라이어 제거, 혼란 변수등을 잘 해야한다는 것입니다.
    • · 원하는 메트릭이 Stakeholder 마다 다르기 때문에 여러가지 메트릭에 대해 성능을 테스트해서 로버스트 한지를 봐야합니다.


2. Stakeholder가 그것을 사용할만한 Rationale이 있어야합니다.

  • · 이들은 Small-test를 한다거나 직관, 사고실험 등으로 신뢰할만한지를 판단하기 때문에 이런것들에 도움을 주어야 합니다.
  • · 가장 좋은 Rationale은 인과적인 설명입니다. 그래서 causal relationship 에 대한 설명이 반드시 이루어저야한다고 말을 하고 있습니다. 예를 들어, Breast cacner를 높인다고 알려진 Mammographic texture를 뽑아내는 deep learning system을 개발 한 경우, 이것의 Breast cancer에 대한 causal effect size를 명확하게 제시를 해주어야합니다.

결론

1. 의생명 분야는 복잡하기 때문에 정확하게 이해하기 힘들고, 그렇기 때문에 AI의 서포트가 필요합니다.

2. 딥러닝은 수많은 데이터셋을 포함하여 복잡한 모델링을 하는데 유망한 방법이고, 그렇기 때문에 Deep learning 은 의생명 분야에서 중요한 역할을 할 것입니다. 



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에서 보이는 생존률에 대한 경향성이 시간이 지날 수록 강해지고 있는 것이다. 이는 위의 카플란 마이서 생존곡선에서도 확인할 수 있다. 



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/


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을 추정할 수 없다. 층화한 변수와 다른 변수와의 교호작용을 고려하여야한다. 교호작용이 있는 경우 추정된 계수값은 바이어스가 생긴다. 







티스토리 툴바