This paper describes survival analysis methods and shows how the methods can be applied to analyses of ophthalmologic data using the open source statistical program R. The R codes are presented in the Appendix. We discuss some papers published in the ophthalmology area that used survival analysis methods. We hope that researchers interested in applying survival analysis to studies of ophthalmology develop easier to understand survival analysis methods and be able to apply the appropriate methods for their research purposes.
생존분석은 관심 사건(event of interest)이 발생하기까지 걸린 시간을 분석하는 통계적 방법론으로, 중도탈락을 포함한 모든 연구대상자에서 생존시간과 관심 사건 발생 여부를 종속변수로 하여 생존함수를 추정하거나 또는 생존함수에 유의한 영향을 미치는 독립변수를 파악하고 그 효과를 추정하는 데 사용 된다. 예를 들어, 혼탁이 생긴 각막 대신 투명한 각막을 이식하여 시력을 개선하는 각막이식술을 살펴보자. 이식편이 투명성을 유지하지 못한 이식 실패가 관심사건으로, 특정 시점에서 이식 실패 여부를 종속변수로 하여 로지스틱 회귀분석으로 분석할 수 있다. 특정 시점 전 중도탈락한 연구대상자들은 특정 시점에서의 이식 실패 여부를 모르므로 로지스틱 회귀분석에서 제외하게 된다. 이 경우 중도탈락할 때까지 투명성을 유지했다는 정보를 활용하지 못하며, 선택 편향(selection bias)이 일어 난다. 한편, 수술 이후 이식 실패까지의 시간(투명성을 유지한 시간, 이식편의 생존시간)을 연구하고자 한다면, 이식편의 생존시간은 연속형 변수이므로 회귀분석으로 분석할 수 있다. 그러나 양의 값을 가지는 시간 자료는 극단으로 치우친 값이 있을 경우 정규분포를 따르지 않아 회귀분석을 적용하기 어렵다. 또한 연구기간 동안 이식 실패가 발생하지 않아 생존시간을 알 수 없는 경우도 회귀분석 시 문제가 된다. 이처럼 관심사건 발생까지의 시간을 고려하며 중도탈락한 자료가 포함되는 연구는 생존분석을 적용하는 것이 적절하다.
본 논문에서는 생존자료의 특성을 설명하고 생존분석기법을 소개한다. 또한 안과학 논문에서 생존분석이 어떻게 적용되었는지 예제를 통해 살펴보고자 한다. 이와 관련된 R 프로그램 코드는 부록에 수록하였다.
생존자료는 생존시간의 시작점과 관심 사건 발생 시점을 모두 관찰한 완결 자료(complete data)와 두 시점 중 한 시점이라도 정확히 관찰하지 못하여 생존시간을 불완전하게 반영하는 절단자료(censored data)를 포함한다. 관심 사건이 발생하지 않은 채 연구가 종료되거나 예상하지 못한 이유(이사, 전원, 다른 질병이나 사고, 치료 거부 등)로 관찰이 중단될 때 절단(censoring)이 일어난다. 이런 절단자료는 절단 직전까지 관심 사건이 발생하지 않았다는 부분 정보를 제공한다.
시간에 따라 이식편이 투명성을 유지할 확률, 즉 이식편 생존율을 구하기 위해 각막이식술을 시행한 10명을 10년 동안 관찰하여 연구종료시점에 6명의 이식 실패, 1명의 이식 성공, 3명의 중도탈락 자료를 얻었다고 하자(Fig. 1). 시작점과 관심 사건 발생 시점을 모두 정확히 관찰한 4, 5, 7, 9번 환자는 완결자료이다. 1번 환자는 연구 이전 각막이식술을 시행하여 정확한 시작점을 모르는 좌중도절단자료이며, 8번 환자는 관심 사건 발생 시점이 7-9년 사이라는 정보만 아는 구간중도절단자료이다. 그리고 이식 성공으로 분류된 3번 환자와 관찰이 중단된 2, 6, 10번 환자의 자료는 우중도절단자료이다. 절단자료의 종류 (
Example of censored data. The symbol + of table and X of graph mean the censored data and the occurrence of the event, respectively. ID = identification.
본 논문은 임상연구에서 가장 흔한 우중도절단자료를 다룬다. 좌중도절단자료는 연구 개시일을 생존시간의 시작점으로 가정할 수 있다. 구간중도절단자료는 구간 중 어느 시점을 사건 발생 시점으로 정하느냐에 따라 결과에 선택 편향이 발생할 수 있다. 구간중도절단자료가 포함된 자료를 연구할 때 Yun and Kim [1]과 Kim et al. [2]을 참고할 수 있다.
발생률은 특정 집단에서 특정 시간 구간 내 사건 발생 빈도를 정량화한 지표이다. 사건 발생까지의 시간 정보 수집 유무에 따라 누적발생률(proportion cumulative incidence) 또는 발생률 (incidence rate, incidence density rate)로 나뉜다. 누적발생률은 특정 집단의 수를 분모로, 특정 구간 내 관심 사건 발생 빈도를 분자로 계산한 관심 사건 발생분율이다. 만약 관심 사건이 발생한 시간 정보를 알고 있다면 이 집단에서 관찰 기간 동안 사건이 얼마나 빠르게 발생했는지 확인할 수 있다. 발생률은 이런 관심 사건의 총 생존시간을 분모로, 관심 사건 발생 빈도를 분자로 계산한 관심 사건 발생 비율이다. 예를 들어, A와 B 두 각막이식술 후 10년간 이식 실패 발생률을 비교하는 연구를 진행한다고 하자(Table 1). 두 군에서 각각 10안 중 8안에서 이식 실패가 발생하였다면, 10년 동안 이식 실패가 일어난 안의 비율인 누적발생률은 0.8 [관심 사건 발생 수/집단의 수]로 동일하다. 여기에 시간 정보를 추가하여 각 집단에서 이식 실패가 얼마나 빨리 일어나는지 발생률로 비교할 수 있다. 이 예에서 발생률은 각막이식술을 받은 안의 생존시간의 합에 대해 10년 동안 발생한 이식 실패 건수로 계산된다. 두 군의 발생률은 각각 0.21 [관심 사건 발생 수/총 생존시간]과 0.10 [관심 사건 발생 수/총 생존시간]으로, 초반에 사건이 더 많이 발생한 A군이 더 높다.
한편, 생존율은 특정 시간 이후 생존할 확률이다. 생존분석은 시간의 함수인 생존율이 관심변수로, 생존곡선(survival curve)과 위험함수(hazard function)로 시간에 따른 생존율의 양상을 비교한다. 생존곡선은 시간에 따른 생존율을 나타낸(year)그래프이며 위험함수는 특정 시점까지 생존한 직후 사망 할 확률함수이다.
생존분석의 이해를 위해 당뇨병 환자의 시력 손실을 늦추는 레이저치료법의 효과를 확인하려 한 ‘retinopathy’ (R 패키지 ‘survival’에서 발췌) 자료를 살펴보자. 총 197명의 당뇨병 환자에서 한쪽 안에는 레이저치료를 진행하고 다른 안에는 치료를 진행하지 않은 채 각 안의 시력상실까지의 시간을 관찰하였다 (Table 2). Appendix에서 이 자료를 생존분석으로 분석하는 방법을 R프로그램 기반으로 설명하였다.
생존곡선은 시간 0에서 100% 생존율로 시작하여 시간이 지나면서 생존율이 감소하는 생존함수를 나타내는 곡선이다. 특정 시점의 생존율이 아니라 모든 시점의 생존율을 비교하게 된다. 즉 생존곡선은 모든 시간의 누적생존율 곡선이다. 누적생존율은 각 구간의 생존율을 곱하여 계산하며, 각 구간마다 유지되거나 감소하는 계단형태로 표현된다. 구간을 정의하는 방법에 따라 연구자가 임의로 구간을 정하는 생명표법과 사건이 일어난 시점을 기준으로 구간을 정하는 카플란마이어 분석법[3]이 있다. 생명표법은 충분히 큰 표본(군당 50 이상)에 사용되며, 카플란마이어 분석법은 작은 표본에도 적용할 수 있고 절단자료의 처리가 쉽다. 카플란마이어 누적생존율을 계산하는 방법은
간혹 생존율의 반대 개념으로 시간에 따른 사망률(1-생존율), 즉 사망 곡선(mortality curve)을 제시하고 싶을 때가 있다. 시력상실 발생 비율이 적어 생존곡선으로 생존율 감소 양상이 명확하지 않을 때 시력상실곡선으로 나타낼 수 있다. 또는 질병 회복 같은 긍정적 성격을 갖는 관심 사건은 시간에 따라 회복률이 높아지는 양상을 나타내기 위해 병의 회복곡선으로 제시할 수 있다.
군별 생존곡선이 차이가 있는지 통계적으로 검정하기 위해 로그순위 분석법을 적용한다(
로그순위 검정통계량과 유의확률은 ‘survdiff’ 함수로 구하며 유의확률이 유의수준보다 작다면 귀무가설을 기각한다. 즉 레이저 치료를 한 군의 생존곡선이 레이저 치료를 하지 않은 군의 생존곡선보다 느리게 감소하며, 로그순위 분석법을 통해 레이저 치료가 당뇨병 환자의 시력 손실을 늦추는 데 효과가 있음을 통계적으로 주장할 수 있다(Fig. 2, Appendix). 이를 논문에 표기할 때는 로그순위 검정통계량과 유의확률을 같이 제시한다.
Survival curves with log-rank
로그순위 분석법은 모든 연구 시점에 같은 가중치를 주고 두 집단 또는 여러 집단 간의 생존함수를 비교한다. 하지만 연구 목적에 따라 연구 초기 시점 또는 후기 시점에 가중치를 주고 집단 간의 생존함수를 비교할 수 있는데, 이를 가중로그순위 분석법이라 한다. 가중로그순위 분석법은 부여하는 가중치 종류에 따라 다양하며 연구 초반의 생존함수 차이를 부각시키기 위해 초기 시점에 가중치를 부여하는 윌콕슨 분석법(Wilcoxon method) 등이 있다[4,5].
로그순위 분석법으로 레이저 치료 여부에 따른 두 군의 생존 함수를 비교하였다. 하지만 로그순위 분석법은 단변수 분석법 (univariable analysis)으로 두 개 이상의 독립변수가 생존함수 또는 위험함수에 미치는 영향을 파악하고자 할 때는 사용할 수 없다. 여러 독립변수들이 생존함수에 미치는 영향을 조사하는 다변수 분석법(multivariable analysis)으로 콕스의 비례위험모형[6]이 있다. 로그순위 분석법은 비례위험모형 적용 전 각각의 변수가 생존함수 또는 위험함수에 미치는 개별적인 영향을 보기 위한 기초적인 분석 자료로 활용할 수 있다.
비례위험모형은 위험함수비율(hazard ratio [HR])을 통해 독립변수가 위험률(hazard rate)에 미치는 영향을 보인다. 위험률이란 시점까지 생존한 사람이 그 직후 사망할 확률, 즉 순간사망률을 의미하며, HR은 이 위험률의 비로 정의된다. 예를 들어, 레이저 치료군과 비치료군을 각각 실험군과 비교(통제)군으로 설정했다면, HR은 레이저 비치료군에 대한 레이저 치료군의 위험률 비로 계산한다. 일반적으로 HR > 1일 때 비교군보다 치료군이 더 위험함을 의미하며 HR < 1일 때 덜 위험함을 의미한다. 하지만 관심 사건이 긍정 사건이면 HR의 해석을 조심해야 한다. 예를 들어, 관심 사건이 회복(remission)이라면 HR > 1일 때 비교군보다 치료군이 회복 가능성이 높음을 의미한다.
비례위험모형은 용어 정의에서 드러나듯 HR이 시간에 관계없이 일정하다는 가정을 만족해야 한다. 위 예제처럼 두 집단의 비교라면 각 집단별 log (-log s(t)) 곡선이 서로 평행하면 이는 비례위험가정을 만족함을 의미한다. 만약 비례위험가정이 만족되지 않으면 변형된 방법인 층화된(stratified) 비례위험모형을 사용하거나 시간 가변 독립변수(time-dependent covariate)를 포함시켜 분석할 수 있다. 이 경우는 본 논문의 범위를 벗어 나므로 자세한 설명은 생략한다.
비례위험모형을 만들기 위해 수집한 여러 독립변수 중 통계적으로 유의한 변수(관심 있는 사건의 위험함수에 유의한 영향을 미치는 변수)를 찾아야 한다. 찾아진 변수들을 이용하여 비례위험모형을 만든다. 이때 선택된 변수들의 조합으로 만들어지는 비례위험모형 중 가장 적합한 모형을 선택해야 한다. 설명력과 모형의 복잡도를 고려하여 최종 모형을 선택한다.
모형에 포함될 변수를 선택하는 전통적 방법으로 전진선택법(forward selection), 후진제거법(backward selection), 단계별 선택법(stepwise selection)이 있다. 하지만 이러한 방법들은 중요한 변수가 제거될 가능성이 있고 변수들이 많은 경우 계산이 많아진다. 다양한 변수선택법이 연구되고 있으므로 연구 주제에 따라 적절하게 적용해야 한다. 연구자의 판단 혹은 기존 결과에서 중요하다고 확인된 변수는 통계적으로 유의하지 않더라도 포함할 수 있으며 중요한 교호작용이 있으면 선택할 수 있다. 선택된 변수는 마팅게일잔차(Martingale residual)와 선형 관계가 되도록 모형에 적합하여야 하는데 이를 확인하기 위해 각 독립변수들과 위험함수와의 로그선형(log-linear) 관계를 확인하여 적절한 형태를 선택한다(
‘Retinopathy’ 자료로 ‘당뇨병 진단 연령에 따라 레이저 치료 여부에 대한 효과가 다르다’라는 가설을 검정한다고 하자. 이를 위해 레이저 치료 여부에 따라 효과가 다르고, 당뇨병 진단 연령에 따라 추가적으로 효과가 달라지는 상호작용을 고려한 모형을 세워야 한다. 이 경우에는 나이와 치료 여부(trt 변수)를 고려해야 한다. 본 예제에서는 진단 당시 나이를 기준으로 20세 미만(trt = juvenile 군)과 20세 이상(trt = adult 군)으로 구별 해 놓은 type 변수를 사용하여 콕스의 비례위험모형을 적용한다. 로그선형관계를 고려하여 모든 변수를 선형 형태로 표현하였다.
모형에 포함할 변수를 선택하는데, R 프로그램에서는 전진선택법과 후진제거법, 단계별선택법을 제공한다. 이 예제에 단계별 선택법을 적용하여 모형에 포함된 변수는 레이저 치료 여부(trt)와 당뇨병 종류(type) 그리고 그 교호작용이며, 최종모형으로부터 추정된 위험함수는
이다(Appendix). 위의 위험함수는 20세 이상(type = adult)에서는 레이저 치료군(trt=1)이 비치료군(trt=0)에 비해 0.40배 (=exp(0.341-0.425-0.846))의 위험률을 가지며 20세 미만(type = juvenile)에서는 레이저 치료군이 비치료군에 비해 0.65배(=exp(- 0.425))의 위험률을 가짐을 의미한다. 즉 20세 이상과 20세 미만 모두 치료 효과가 있지만 20세 이상에서 더 큼을 알 수 있다.
최종모형에 포함된 독립변수들이 비례위험가정을 만족하는지 검토하기 위하여 ‘cox.zph’ 함수를 이용한다. 위 예제에서 모형에 포함된 모든 독립변수들에 대하여 유의수준 0.05에서 ‘비례위험가정을 만족한다’는 귀무가설을 기각할 수 없고, 따라서 모형에 포함된 모든 변수는 비례위험가정을 만족한다고 주장할 수 있다(Appendix).
만일 반복 측정된 자료나 다기관 임상시험자료 같이 군집 내에 연관성이 있어 생존시간의 독립을 가정할 수 없으면 비례위험모형의 확장 모형인 프레일티모형(frailty model)을 적용할 수 있다. R에서 ‘survival’ 패키지 또는 프레일티 모형의 변수 선택 및 모형 적합을 위해 개발된 ‘frailtyHL’ 패키지 등에서 적용할 수 있다[7,8].
생존분석을 적용한 안과학 분야 논문들을 소개하여 안과학 분야에서 생존분석이 어떻게 사용되고 있는지 보여주고자 한다.
Kim et al. [9]은 국내공여각막과 해외공여각막의 유효성을 확인할 목적으로 2004년 11월부터 2011년 8월까지 전층각막이식술을 시행받은 총 211명의 236안에서 국내공여각막과 해외공여각막의 전층각막이식술 후 임상 결과를 비교하였다. 관심 사건을 원인과 상관 없이 이식편이 혼탁해지는 것으로 정의하고, 관심 사건이 관찰된 이식편의 생존시간은 수술한 날부터 실패의 징후가 나타난 날까지, 관심 사건이 관찰되지 않은 이식편의 생존시간은 수술한 날부터 마지막 내원일까지로 정의하였다. 각 군에 대한 카플란마이어 생존곡선을 그렸고 두 군의 생존곡선을 비교한 로그순위 분석법의 유의확률이 유의수준 0.05보다 커 통계적으로 다르다고 주장할 수 없었다. 즉 국내공여각막을 이식받은 환자군과 해외공여각막을 이식받은 환자군 사이에 유의한 임상적 차이가 있지 않다고 주장할 수 있다.
성별, 나이, 전층각막이식술 전 녹내장 수술 병력 및 수술 전 후의 수정체 상태를 나타내는 변수들이 생존곡선에 미치는 영향을 로그순위 분석법으로 살펴본 결과 각 변수들이 생존곡선에 미치는 영향은 통계적으로 유의하지 않았다. 그리고 공여편 내피세포 밀도에 영향을 줄 수 있는 백내장수술을 병행 또는 추후 시행한 군에 대하여 추가 분석하였는데, 생존율을 비교 한 로그순위 분석법의 유의확률이 유의수준 0.05보다 커 통계적으로 유의한 차이가 있다고 할 수 없었다.
Park et al. [10]은 망막모세포종 환자의 생존율 연구를 계획하고 국립암센터 자료를 활용하여 1993년부터 2010년까지 총 600 명을 대상으로 진단 연도별, 성별, 진단 나이별 발생률 및 생존율을 비교하였다. 관심 사건은 사망으로 정의하였고, 생존시간은 첫 진단일부터 사망일, 관찰이 중단된 날 또는 관찰 종료 시점으로 정의하였다.
600명의 생존곡선을 제시하고, 진단 시기별, 성별, 진단 나이 별로 각 군에 대한 생존곡선과 로그순위 분석법을 시행한 결과를 그래프에 함께 제시하였다. 위에서 고려한 세 가지 변수가 생존율에 미치는 영향을 함께 살펴보려면 콕스의 비례위험 모형을 적용할 수 있다. 하지만 성별, 진단 나이별 생존곡선이 교차하고 있어 비례위험가정을 만족하지 않을 수 있으므로 적용에 조심해야 한다.
Jeon et al. [11]은 괴사성공막염의 치료를 위한 면역억제제의 효능과 안전성을 알아보기 위해 다기관 연구를 수행하였다. 11개의 3차의료기관에서 2002년부터 2012년까지의 괴사성공막염 환자 50명의 의무기록으로부터 cyclophosphamide를 사용한 그룹과 그 외 면역억제제를 사용한 그룹으로 나누어 회복률, 재발률, 시력상실률 등을 비교하였다.
이 연구는 총 50명의 환자의 두 군별 완전회복(Complete Remission) 및 부분 회복(Partial Remission) 여부와 시간을 관찰한 것으로 각 군별 회복 여부를 Fisher’s exact test로 검정하고 Mantel-Haenszel test로 비교하였다. 그리고 카플란마이어 생존곡선을 제시하여 각 군의 회복률을 제시하였다.
이 자료는 회복 여부와 시간을 종합하여 생존분석으로 분석할 수 있다. 각 군별 회복률의 차이를 확인하기 위해 카플란마이어 생존곡선을 그리고 로그순위 분석법으로 그 차이를 비교 한다. 이때 이 연구의 관심 사건이 긍정적 의미를 지니고 있는 회복이므로 해석에 주의해야 한다. 이 논문에서 생존율은 회복되지 않은 상태(병증이 있는 상태)로 남아있는 비율을 나타내며 생존곡선이 빠르게 감소할수록 환자가 빠르게 회복되었다는 것을 의미한다. 이 경우 발생(회복)곡선을 제시하여 시간에 따른 회복률의 증가 양상을 제시할 수 있다(Fig. 3).
Remission rate curve of 50 patients with necrotizing scleritis treated with immunosuppressive agents in CYC and OISA groups. No significant difference was seen between CYC group and OISA group in either complete remission (
본 논문의 목적은 안과영역 연구자들이 생존분석 자료 분석 시 적절한 생존분석방법을 적용할 수 있도록 돕는 것이다. 이를 위해 생존분석방법론 중 로그순위 분석법과 콕스의 비례위험모형에 대하여 설명하고 안과 영역의 예제 자료를 통계프로그램인 R로 분석할 수 있게 하였다. 실제 출판된 논문들을 통해 현재 연구 영역에서 생존분석을 어떻게 적용하였는지 확인하여 봄으로써 연구자들이 각자의 연구를 설계할 때 생존분석 방법을 이해하고 적절한 방법을 선택하여 분석할 수 있기를 바란다.