종단 생존시간 자료의 결합 모형에 대한 고찰: R 패키지 JM을 중심으로

Joint Models for Longitudinal and Time-to-event Data with JM R Package

Article information

J Health Info Stat. 2023;48(4):291-305
Publication date (electronic) : 2023 November 30
doi : https://doi.org/10.21032/jhis.2023.48.4.291
Professor, Department of Applied Statistics, University of Suwon, Hwaseong, Korea
김진흠,orcid_icon
수원대학교 데이터과학부 교수
Corresponding author: Jinheum Kim. 17 Wauan-gil, Bongdam-eup, Hwaseong 18323, Korea
E-mail: jkimdt65@gmail.com
*This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. RS-2023-00239958).
No potential conflict of interest relevant to this article was reported.
Received 2023 August 23; Revised 2023 September 25; Accepted 2023 November 30.

Trans Abstract

Even if the value of the risk factor evolves continuously in the extended Cox model, the inference may be somewhat biased because it employs a discrete approximation method. As a result, if the risk factor's value not only fluctuates constantly but also contains measurement error, a model that substitutes the average value rather than the measured value of the risk factor might be considered. Such a model is known as a joint model. After introducing the most widely used the present-value model among joint models, this model was extended to various types of data. In addition, several residuals for model diagnosis were introduced, methods for predicting the probability of event occurrence and the value of the longitudinal risk factor by application of the estimation model were introduced, and an index that can evaluate how well the longitudinal risk factors divide patients into high-risk and low-risk groups was introduced. Finally, all the statistical inference methods introduced in this study were implemented using the JM R package and the source codes were supplied.

서 론

생존시간은 연구시작 시점부터 관심 있는 이벤트(예: 사망, 재발 등)의 발생 시점까지의 시간으로 정의한다. 생존자료분석에서는 생존시간에 영향을 미치는 위험인자를 찾기 위해 콕스(Cox) 비례위험모형을 널리 사용하고 있다[1,2]. 그러나 대부분의 종단 연구에서는 환자가 방문할 때마다 관심 있는 이벤트의 발생 여부뿐만 아니라 위험인자를 관측한다. 따라서 시간에 따라 변하는 위험인자가 있다면 위험인자의 기저 값과 함께 주기적으로 관측된 값을 모형에 포함할 수 있는데 이와 같은 모형을 확장 콕스(extended Cox) 모형이라고 한다[3]. 다만 이 모형은 위험비(hazards ratio)가 시간에 따라 변할 수 있기 때문에 더 이상 비례위험 가정을 만족하지 못할 수도 있다. 한편 확장 콕스 모형에서는 위험인자의 값이 연속적으로 변해도 이산적으로 근사하는 방법을 사용하기 때문에 다소 편향된 추론이 될 수 있다[4]. 따라서 위험인자의 값이 연속적으로 변할 뿐만 아니라 측정 오차를 포함하고 있으면 위험인자의 측정값보다는 평균값으로 대체하는 모형을 생각할 수 있는데 이와 같은 모형을 소위 결합 모형(joint model, JM)이라고 한다[5]. 이때 종단 위험인자에 대한 모형으로는 선형혼합모형이 널리 사용하고 있다.

2절에서는 본 논문에서 사용할 자료를 소개하고, 3절에서는 결합 모형 중에서 가장 널리 사용되는 현재 값 모형을 소개한다. 4절에서는 현재 값 모형을 여러 가지 유형의 자료로 확장하고, 5절에서는 모형 진단을 위한 여러 가지 잔차를 소개한다. 6절에서는 추정 모형의 응용으로 이벤트의 발생 확률과 종단 위험인자의 값을 예측하는 방법을 소개하고, 종단 위험인자가 환자들을 두 그룹(고위험과 저위험)으로 얼마나 잘 구분하는지를 평가할 수 있는 지수 등을 소개한다. 마지막으로 7절에서는 본 연구의 내용을 간략히 요약하고 본 논문의 제한점을 제시하고자 한다.

자료 소개

본 논문에서 사용하고자 하는 자료는 원발성 담낭 간경변증(primary biliary cirrhosis, PBC)을 가진 환자들을 대상으로 1974년부터 1984년까지 미국 메이오 클리닉(Mayo Clinic)에서 임상시험을 수행한 결과이다[6]. 디-페니실라민(D-penicillamine, DPCA)이 PBC 환자들의 생존시간을 향상시키는지 알아보기 위해 연구에 참여한 312명 중에서 임의추출한 158명에게는 DPCA을 처방하고 나머지 154명에게는 위약을 처방한 뒤 등록 시점부터 사망이나 간이식, 종료 시점(1986년 7월)까지의 시간(단위: 년)과 발생한 이벤트의 유형을 기록하였다. 125명이 사망했고, 19명이 간이식을 받았으며, 나머지 168명은 중도 탈락(lost to follow-up)하거나 중도절단되었다. 위험인자는 성별, 연령뿐만 아니라 혈청 내 빌리루빈(serBilir), 혈관 기형 여부(spider), 간비대 여부(hepatomegaly) 등이다. 이 자료는 JM 패키지[7]에서 얻을 수 있다.

Figure 1은 치료 방법(DPCA 대 위약)에 따라 생존함수의 카플란-마이어(Kaplan-Meier) 추정량[8]과 위험집합의 크기를 나타낸 그림이며, 이에 대한 R 코드는 다음과 같다.

Figure 1.

Kalan-Meier estimates of the survival function and the number of risk set at the selected times by drug for the PBC data. PBC, primary biliary cirrhosis.

library(“JM”)

library(“survminer”)

KM = survfit(Surv(years, status2) ~ drug, data = pbc2.id)

ggsurvplot(KM, pval = TRUE, xlab = “Time(years)”, risk.table = TRUE, legend.title = “Treatment”, legend.labs = c(“Placebo”, “DPCA”))

로그-순위검정 결과 유의수준 5%에서 치료 방법 간에 생존함수가 다르다고 할 수 없었다(p =0.99).

종단 생존시간 자료의 결합 모형

공변량은 외적 공변량과 내적 공변량으로 나누어진다. 외적 공변량은 다시 성별 등과 같이 시간에 따라 변하지 않는 공변량과 시간에 따라 변하는 공변량으로 나누어진다. 특히 후자는 관심 있는 이벤트의 발생에는 영향을 미치지만 이벤트의 발생으로 인해 영향을 받지 않으므로 이벤트의 발생 이후에도 관측이 가능한 공변량이다. 내적 공변량이란 환자가 위험집합(risk set)에 있는 동안에만 관측이 가능한 공변량을 말한다[9]. PBC 자료에서 성별, 등록 시 연령, 치료 방법 등은 시간에 따라 변하지 않으므로 외적 공변량이고, serBilir, spider, hepatomegaly 등과 같은 바이오마커(biomarkers)들은 내적 공변량이다. 외적 공변량이 이벤트의 발생에 미치는 영향을 알아보려면 확장 콕스 모형을 사용하면 된다[3]. 다만 내적 공변량은 소위 “ last value carried for-ward”라는 원칙에 따라 인접한 관측 시점들 사이 구간에서는 공변량의 값이 변하지 않는다고 가정한다. 따라서 내적 공변량을 가진 확장 콕스 모형의 통계적 추론은 편향된 결론을 내릴 수 있다[4].

본 절에서는 내적 공변량이 생존시간에 미치는 영향을 알아보기 위해 널리 사용되고 있는 결합 모형을 소개하고자 한다. 결합 모형이란 생존시간 자료에 대한 모형(이하, 생존 모형)과 내적 공변량처럼 반복 측정된 종단 자료에 대한 모형(이하, 종단 모형)을 결합한 모형을 말한다.

첨자 i(i=1,2,…,n)와j(j=1,2,…,ni) 는 각각 개체와 개체 내 관측 번호를 나타내고,T* 는 이벤트 발생 시간,C 는 중도절단 시간이라 고 하자. 본 논문에서 다룰 종단 생존시간 자료는 다음과 같다.

Dn={(Ti,δi,yi);i=1,2,,n}

여기서 Ti=min(T*i,Ci),δi = I(T*i≤Ci) 이고 yi = (yi(ti1),yi(ti2),…,Yi(tin1))' 는ni × 1 종단 자료 값 벡터이다. 또한Yi(t) 는 t 시점에서 관측한 i-번째 개체의 종단 자료 값이다.

mi(t)는 시점t 에서 i-번째 개체의 관측할 수 없는 참(true) 종단 값이라고 하고,m~i(t)={mi(s),0s<t} 는t 바로 직전 시점까지 종단 자료 값의 프로세스라고 하자. 먼저 생존 모형을 다음과 같이 정의하자.

(1) hi(t|m~i(t),wi)=h0(t)exp(γwi+αmi(t))

여기서 h0(t) 는 기저위험함수이고 Wi 는 외적 공변량 벡터이다. h0(t) 는 와이블분포나 조각지수분포, B-스플라인 등처럼 모수적으로 가정 하거나 비모수적으로 가정할 수 있다. 한편 종단 모형을 다음과 같이 정의하자.

(2) yi(t)=mi(t)+ϵi(t)=βxi(t)+bizi(t)+ϵi(t).

여기서 xi(t)는 고정 효과에 대한 디자인 벡터이고, zi(t) 는 (1,t)',(1,t,t2)' 등처럼 랜덤 효과에 대한 디자인 벡터이다. 랜덤 효과 벡터 bi ∼ N(0,D) 라고 가정하고, 측정 오차는 εi(t) ∼ N (0,σ2) 이라고 가정하자. bi와 εi(t) 는 서로 독립이라고 가정하자.

모형 (1)과 (2)에서 추정해야 할 모수 벡터는 θ = ( θ't, θ'y, θ'b)이다. 여기서 θt = (γ',α, θ'h0)' 은 생존 모형에 관한 모수 벡터이고, θy = (β'2)' 은 종단 모형에 관한 모수 벡터, θb = vec(D) 는 랜덤 효과에 관한 모수 벡터이다. 여기서 θh0 은 기저위험함수에 대응하는 모수 벡터를 의미한다. 특히 α 는 바이오마커와 이벤트 발생 위험 간의 연관성을 나타내는 모수이고 본 논문에서 관심 있는 가설은 H0:α = 0 이다.

결합 모형에 대한 우도 함수를 정의하기 위해 다음 세 가지를 가정하자.

  • (A1) P(Ti,δi,yi|bi;θ)=P(Ti,δi|bi;θ)P(y1|bi;θ).

  • (A2) P(yi|bi;θ)=j=1niP(yi(tij)|bi;θ).

  • (A3) 중도절단 시간과 관측 시점이 서로 독립이다.

특히 (A2)는 이벤트 발생 이후 결측 되는 종단 자료 값에 무관하게 추론이 가능하도록 해준다. 따라서 i-번째 개체의 (관측 자료) 우도 함수는 다음과 같이 표현된다.

(3) P(Ti,δi,yi;θ)=P(Ti,δi,yi,bi;θ)dbi=P(Ti,δi|bi;θt)P(yi|bi;θy)P(bi;θb)dbi=hi(Ti|m~i(Ti),wi;θt)δiSi(Ti|m~i(Ti),wi;θt)P(yi|bi;θy)P(bi;θb)dbi=(h0( Ti)exp(γwi+αmi(Ti)))δiexp(0Tih0( s)exp(γwi+αmi(s))ds)×j=1ni(2πσ2)12exp(12σ2(yi(tij)mi (tij))2)×(2π)r2det(D)12exp(12 biD1 bi)dbi.

θ에 대한 최대우도추정량은 EM[10] 알고리즘을 써서 로그-우도함수

logL(θ)=i=1nlogP(Ti,δi,yi;θ)

를 최대로 하는 θ로 정의된다. 그러나 EM 알고리즘은 추정량의 수렴 속도가 선형적이라는 단점이 있다[11]. (관찰 자료) 스코어 벡터는 (완전 자료) 스코어 벡터의 랜덤 효과에 대한 사후 분포의 기댓값으로 표현되는데 이를 구체적으로 표현하면 다음과 같다[5].

(4) S(θ)=i=1nθlogP(Ti,δi,yi;θ)=i=1nA(θ,bi)P(bi|Ti,δi,yi,;θ)dbi.

여기서

A(θ,bi)=θlogP(Ti,δi,yi,bi;θ)=θ{logP(Ti,δi|bi;θt)+logP(yi|bi;θy)+logP(bi;θb)}

이다. θ 의 최대우도추정량을 θ^ 라고 하면θ^ 의 공분산 행렬은 식 (4)를 써서 다음과 같이 추정된다.

(5) var(θ^)=I(θ^)1

여기서 I(θ^)=θS(θ)|(θ=θ^) 는 관찰 자료 정보 행렬이다. 한편 콕스 모형이나 확장 콕스 모형은 부분 우도(partial likelihood)로부터 회귀 모수를 추정하므로 기저위험함수를 명시하지 않아도 회귀 모수에 대한 추론이 가능하다[3]. 그러나 결합 모형에서는 부분 우도를 사용하지 못하고 완전 우도(full likelihood)를 사용하므로 누적 기저위험함수 즉H0(t)=0th0(s)ds 가 서로 다른 이벤트 발생 시점에서만 변하는 계단 함수라는 가정하에서 추정해야 한다[12,13]. 그러나 n이 크면 θh0 의 차원이 매우 커서 식 (5)의 계산이 어려우므로 θh0 의 추정량을 프로파일 추정량θ^h0 으로 대체한 후 θh0 를 제외한 나머지 모수의 추정량의 표준오차만 추정하는 방법을 사용한다[12].θ^h0 의 표준오차는 부스트랩(bootstrap)을 사용하여 추정할 수 있지만 결합 모형에서는 시간이 많이 요구되는 단점이 있다[12].

PBC 자료에서 질병 악화의 지시자(indicator)로 serBilir를 고려하고자 한다. 기저 시점, 6개월 후, 1년 후, 그 이후부터는 매년 한 번씩 serBilir를 관측하여 312명의 환자로부터 총 1,945개의 측정값을 얻었다. 환자 1인당 평균 6.2번 반복 측정한 셈이다. Figure 2는 특정 환자들의 serBilir 값(동그라미)과 평활 값(파란색: 위약, 오렌지색: DPCA)을 나타낸 그림이며, 이에 대한 R 코드는 다음과 같다.

Figure 2.

Subject-specific longitudinal profiles for sixteen patients from PBC data. PBC, primary biliary cirrhosis.

library(“lattice”)

PBC.samp = subset(pbc2, id %in% c(38, 134, 70, 82, 51, 90, 68, 93, 39, 148, 173, 200, 216, 242, 269, 290))

xyplot(log(serBilir) ~ year | id, data = PBC.samp, groups = drug, type = c(“p”, “smooth”), lwd = 2, layout = c(4, 4), as.table = TRUE, ylab = expression(log(“serBilir”)), xlab = “Time (years)”, auto.key = list(title = “Drug”, cex.title = 1, columns = 2, text = c(“placebo”,”DPCA”)))

serBilir의 평활 곡선으로 볼 때 serBilir는 시간(단위: 년)에 따라 비선형적으로 변하고 치료 방법에 따라 변하는 경향이 있다고 할 수 있다. 따라서 다음과 같은 mi(t)를 가진 종단 모형을 가정하고자 한다.

(6) mi(t)=β0+β1×drugi+β2t+β3t2+β4(t×drugi)+β5(t2×drugi)+bi0+bi1t+bi2t2.

또한 생존 모형은 다음과 같이 가정하고자 한다.

(7) hi(t)=h0(t)exp(γ×drugi+α×mi(t)).

Table 1은 기저위험함수에 따라 모형 (6)과 (7)의 결합 모형에 대한 회귀 모수와 연관성 모수의 추정값과 표준오차를 정리한 것이다. 한편 serBilir를 외적 공변량으로 정의하면 모형 (7)은 다음과 같이 표현된다.

Estimates of regression parameters and their standard errors for longitudinal biomarker and/or event parts in the joint models with three different baseline hazard functions and in the extended Cox model

(8) hi(t)=h0(t)exp(γ×drugi+α×serBiliri(t))

Table 1은 확장 콕스 모형 (8)에 대한 회귀 모수의 추정값과 표준오차도 포함하고 있으며, 이에 대한 R 코드는 다음과 같다1)

lmeFit.pbc = lme(log(serBilir) ~ drug * (year + I(year^2)), random = ~ year + I(year^2) | id, data = pbc2) #종단 모형

coxFit.pbc = coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE) #생존 모형

wb.pbc = jointModel(lmeFit.pbc, coxFit.pbc, timeVar = “year”, method = “weibull-PH-aGH”) #와이블분포

jointModel(lmeFit.pbc, coxFit.pbc, timeVar = “year”, method = “piecewise-PH-aGH”) #조각지수분포

spline.pbc = jointModel(lmeFit.pbc, coxFit.pbc, timeVar = “year”, method = “spline-PH-aGH”) #스플라인

tmp = tmerge(pbc2.id, pbc2.id, id = id, cEvent = event(years, status2))

pbc2.ec = tmerge(tmp, pbc2, id = id, serBilir.ec = tdc(year, serBilir), hepatomegaly.ec = tdc(year, hepatomegaly))

tdCox.bilir = coxph(Surv(tstart, tstop, cEvent) ~ drug + log(serBilir.ec), data = pbc2.ec) #확장 콕스 모형

c(AIC = extractAIC(tdCox.bilir)[2], BIC = extractAIC(tdCox.bilir, k = log(tdCox.bilir$nevent))[2]) #AIC & BIC

세 결합 모형에서 회귀 모수의 추론 결과는 기저위험함수에 관계없이 유사하였으며, 가장 작은 AIC [14]와 BIC [15] 값을 가진 와이블 분포 모형이 최적 모형이라고 할 수 있다. 결합 모형과 확장 콕스 모형의 생존시간에 대한 회귀 모수의 추정량을 비교하면 두드러진 차이는 없었지만 연관성 모수는 확장 콕스 모형에서 과소추정됨을 알 수 있다.

결합 모형의 확장

본 절에서는 3절에서 정의한 모형 (1)과 (2)를 확장한 모형을 소개하고자 한다.

첫째, 생존 모형 (1)은 종단 자료의 현재 값 mi(t) 에 의존하므로 모형 (1)을 소위 “현재 값” 모형이라고 부른다. 이벤트의 발생 위험이 종단 자료의 현재 값과 변화율에 의존하는 모형(“현재 값+변화율” 모형이라고 부름)으로 확장할 수 있으며, 이 모형은 모형 (1)을 변형하여 다음과 같이 정의된다.된다.

hi(t|m~i(t),wi)=h0(t)exp(γwi+α1 mi(t)+α2 mi(t))

이외에도

누적값모형:hi(t|m~i(t),wi)=h0(t)exp(γwi+α0tm(s)ds),

지연(lag)모형: hi(tm~i(t),wi)=h0(t)exp(γwi+αmi(max(tc,0)))

(단, c는 양의 상수),

상호작용모형:h1(tm¯1(t),w1)=h0(t)exp(γwi1+α1 m1(t)+α2(w12×m1(t)))

공유 랜덤 효과 모형:hi(tm~i(t),wi)=h0(t)exp(γwi+αbi) 등이 있다

PBC 자료의 경우 mi(t) 의 변화율과 누적 값은 식 (6)로부터 각각 다음과 같이 표현된다.

mi(t)=β2+β3(2t)+β4×drugi+β5(2t×d2ugi)+bi1+bi2(2t),0tmi(s)ds=β0t+β1(t×d2ugi)+β2(t22)+β3(t33)+β4(t22×drugi)+β5(t33×drugi)+bi0t+bi1(t22)+bi2(t33).

Table 2는 모형 (6)과 와이블 기저위험함수를 가진 모형 (7)의 결합 모형에서 모형 (7)을 상술한 모형으로 각각 확장하여 얻은 회귀 모수와 연관성 모수의 추정량과 표준오차를 정리한 것이며, 이에 대한 R 코드는 다음과 같다.2)

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint model under the Weibull baseline hazard function according to the association structure between the longitudinal and event outcomes

dform = list(fixed = ~ I(2 * year) + drug + I(2 * year):drug, indFixed = 3:6, random = ~ I(2 * year), indRandom = 2:3)

update(wb.pbc, parameterization = “both”, derivForm = dform) #현재 값+변화율 모형

iform = list(fixed = ~ -1 + year + I(year * (drug == “D-penicil”)) + I(year ^ 2 / 2) + I(year ^ 3 / 3) + I(year ^ 2 / 2 * (drug == “D-penicil”)) + I(year ^ 3 / 3 * (drug == “D-penicil”)), indFixed = 1:6, random = ~ -1 + year + I(year ^ 2 / 2) + I(year ^ 3 / 3), indRandom = 1:3)

update(wb.pbc, parameterization = “slope”, derivForm = iform) #누적 값 모형

update(wb.pbc, lag = 1) #지연 모형(c=1)

jointModel(lmeFit.pbc, coxFit.pbc, timeVar = “year”, method = “weibull-PH-aGH”, interFact = list(value = ~ hepatomegaly, data = pbc2.id)) #상호작용 모형3)

네 가지 모형 중에서는 “현재 값 + 변화율” 모형의 AIC (또는 BIC) 값이 가장 작았으며, 이 모형은 “현재 값” 모형보다도 AIC (또는 BIC) 값이 작았으며 가설H02 = 0 에 대한 검정 결과도 통계적으로 유의하였다(p <0.001).

둘째, 생존 모형 (1)에서 외적 공변량wi 를wi(t) 로 대체하면 식 (3)에서hi(Ti|m~i(Ti),wi;θt)hi(Ti|m~i(Ti),wi;(Ti)θi) 로 대체하고, 만약wi(t) 의 값이 구간 Ωiq(q = 1,2,…,Qi)에서 wiq 로 일정하다고 하면 Si(Ti|m~i(Ti),wi;θt)를 다음과 같이 대체한다.

Si(t|m~i(t),wi)=exp(0th0( s)exp(γwi(s)+αmi(s))ds)=exp(q=1QiΩiqh0( s)exp(γwiq+αmi(s))ds).

PBC 자료에서 시간에 따라 변하는 hepatomegaly를 모형 (7)에 hepatomegaly를 추가한 모형은 다음과 같이 정의된다.

(9) hi(t)=h0(t)exp(γ1× drug i+γ2× hepatomegaly i(t)+α×mi(t))

Table 3은 모형 (6)과 스플라인 기저위험함수를 가진 모형 (9)의 결합 모형에 대한 회귀 모수와 연관성 모수의 추정량과 표준오차를 정리한 것이며, 이에 대한 R 코드는 다음과 같다.4)

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint model with hepatomegaly (presence/absence) as an external time-dependent covariate

tdCox.hepato = coxph(Surv(tstart, tstop, cEvent) ~ drug + hepatomegaly.ec + cluster(id), data = pbc2.ec, x = TRUE, model = TRUE)

jointModel(lmeFit.pbc, tdCox.hepato, timeVar = “year”, method = “spline-PH-aGH”) #외적 공변량 모형

Hepatomegaly를 가질수록 이벤트의 발생 위험이 1.27배(95% confi-dence interval, CI: 0.81-1.99) 높지만 통계적으로 유의하지는 않았다(p =0.294).

셋째, 다기관 혹은 다국가 연구에서는 기관이나 국가마다 고유한 기저위험함수를 가정할 수 있다. 이때 기관이나 국가를 층으로 정의하고 생존 모형 (1)을 다음과 같이 확장하고자 한다.

(10) hik(t|m~i(t),wi)=hok(t)exp(γwi+αmi(t)),

혹은

(11) hik(t|m~i(t),wi)=h0k(t)exp(γkwi+αkmi(t))

여기서 h0k(t) 는 k-번째 층의 기저위험함수이고, γk와 αk 는 각각 k-번째 층의 회귀 모수와 연관성 모수이다. 모형 (10)과 (11)을 각각 공통 효과 모형과 층별 효과 모형이라고 부르자.

PBC 자료에서 층을 hepatomegaly의 기저 값(“ No”, “ Yes”)으로 정의하자. Table 4는 모형 (6)과 스플라인 기저위험함수를 가진 모형 (10), 모형 (6)과 스플라인 기저위험함수를 가진 모형 (11)의 결합 모형에 대한 회귀 모수와 연관성 모수의 추정량과 표준오차를 각각 정리한 것이며, 이에 대한 R 코드는 다음과 같다.5)

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint models stratified according to the presence or absence of hepatomegaly

stcoxFit.pbc = coxph(Surv(years, status2) ~ drug + strata(hepatomegaly), data = pbc2.id, x = TRUE)

st2coxFit.pbc = coxph(Surv(years, status2) ~ drug * hepatomegaly + strata(hepatomegaly), data = pbc2.id, x = TRUE)

st.pbc = jointModel(lmeFit.pbc, stcoxFit.pbc, timeVar = “year”, method = “spline-PH-aGH”) #공통 효과 모형

st2.pbc = jointModel(lmeFit.pbc, st2coxFit.pbc, interFact = list(value = ~ hepatomegaly, data = pbc2.id), timeVar = “year”, method = “spline-PH-aGH”) #층별 효과 모형

anova(spline.pbc, st.pbc) #우도비 검정

anova(spline.pbc, st2.pbc) #우도비 검정

anova(st.pbc, st2.pbc) #우도비 검정

Hepatomegaly 여부에 따라 기저위험함수가 통계적으로 유의하게 다르며(공통 효과 모형: p =0.002, 층별 효과 모형: p =0.001), 회귀 모수와 연관성 모수가 층별로 다른지에 대한 검정 결과는 유의수준 5%에서 통계적으로 유의하지 않았다(p =0.060). 그러나 hepatomegaly 여부에 따라 이벤트에 미치는 drug 효과는 다르지 않았지만(p =0.288), ser-Billir 효과는 유의하게 다른 것으로 나타났다(p =0.040).

넷째, 한 환자가 종단 연구 기간 동안 경험할 수 있는 이벤트가 2개 이상인 자료를 경쟁위험(competing risks) 자료라고 하는데 생존 모형 (1)을 경쟁위험 자료로 확장하고자 한다.T*k(k = 1,2,…,K) 는 k-번째 이벤트의 참 발생 시간이라고 하면 관측 가능한 자료의 형태는 다음과 같다.

{(Ti=min(T1, T2,,TK,Ci),δi{0,1,2,,K},yi),i=1,2,,n}.

여기서 δi 는 i-번째 환자의 이벤트 유형을 의미하며 중도절단된 경우에는 0으로 정의한다. k-번째 이벤트의 위험함수를 다음과 같이 정의하자. 사실 이는 층별로 기저위험함수와 회귀 모수, 연관성 모수를 고려한 모형과 동일하다.

(12) hik(t|m~i(t),wi)=h0k(t)exp(γkwi+αkmi(t)),k=1,2,,K

따라서 식 (3)에서 P(Tii|bi; θt)를 다음과 같이 대체한다.

P(Ti,δi|bi;θt,β)=k=1K(h0k(Ti)exp(γkwi+αkmi(Ti)))I(δi=k)×exp(k=1K0Tih0k(s)exp(γkwi+αkmi(s))ds).

PBC 자료에서 사망과 이식을 경쟁위험 이벤트로 정의할 수 있다. Table 5는 모형 (6)과 스플라인 기저위험함수를 가진 모형 (12)의 결합 모형에 대해 이벤트의 유형별로 회귀 모수와 연관성 모수의 추정량과 표준오차를 정리한 것이며, 이에 대한 R 코드는 다음과 같다.

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint model with competing risks (transplanted/dead)

pbc2.idCR = crLong(pbc2.id, statusVar = “status”, censLevel = “alive”, nameStrata = “CR”)

CRcoxFit.pbc = coxph(Surv(years, status2) ~ drug * CR + strata(CR), data = pbc2.idCR, x = TRUE)

jointModel(lmeFit.pbc, CRcoxFit.pbc, timeVar = “year”, method = “spline-PH-aGH”, CompRisk = TRUE, interFact = list(value = ~ CR, data = pbc2.idCR)) #경쟁위험 모형

Drug 효과는 사망이나 이식에 통계적으로 유의하지 않았으나(사망: p =0.732, 이식: p =0.267), serBilir의 현재 값이 1 유닛 증가하면 사망의 오즈는 e1.189+0.119 =3.70배(95% CI: 3.05-4.48)이고, 이식의 오즈는 e1.189 =3.28배(95% CI: 2.23-4.83) 높아지는 것으로 나타났다.

다섯째, 다중(multiple) 이벤트의 다른 형태는 같은 이벤트가 반복적으로 발생하는 자료이며 이와 같은 자료를 재발 사건(recurrent event) 자료라고 부른다[16]. 에이즈 연구에서는 실제로 환자들이 사망뿐만 아니라 종단 연구 기간 동안 감염성 질병이 반복적으로 발생할 수도 있다. 에이즈 연구에서 CD4 세포 수가 적을수록 감염성 질병의 발생 횟수는 늘어나고, 또한 CD4 세포 수의 감소와 감염성 질병의 빈번한 발생은 사망의 위험을 증가시키므로 CD4 세포 수와 감염성 질병 발생, 사망을 동시에 고려하는 결합 모형이 요구된다. 따라서 반복 자료에 대한 모형을 추가하고 모형 (1)을 다음과 같이 확장하고자 한다. Uik(k = 1,2,…,Ki)와 dik는 각각 i-번째 환자의 k-번째 재발 시간과 재발 여부라고 하자.

ri(t|mi(t),wi)=r0(t)exp(γrwi+αrmi(t)+vi)hi(t|m~i(t),wi)=h0(t)exp(γhwi+αhmi(t)+ξvi)

여기서 ri(t)는 i-번째 환자의 재발 사건의 강도(intensity)이고, r0(t)는 재발 사건의 기저강도이다. logviN(0,σv2) 은 랜덤 효과이며 프레일티(frailty)라고 부른다. αr 와αh 는 각각 종단 자료 값이 재발 강도 및 이벤트의 발생 위험과 연관된 정도를 나타내는 모수이고, 는 재발 강도가 이벤트의 발생 위험과 얼마나 강하게 연관되어 있는지를 나타내는 모수이다. 3절에서 정의한 모수 θ 를 위 모형에 맞춰 업데이트하고, 가정 (A1)을 다음과 같이 수정하자.

(Al)P(Ti,δi,Ui,di,yi|bi,vi;θ)=P(Ti,δi|bi,vi;θ)P(Ui,di|bi,vi;θ)P(yi|bi;θ).

여기서 Ui = (Ui1,Ui2,…,UiKi)'이고 di = (di1,di2,…,diKi)이다. 또한 다음 가정 (A4)를 추가하자.

 (A4) P(Ui,di|bi,vi;θ)=k=1KiP(Uik,dik|bi,vi;θ)

따라서 관찰 자료에 대한 우도 함수는 3절과 유사한 방법으로 유도할 수 있다.

여섯째, 바이오마커가 두 개 이상인 자료를 분석할 수 있도록 모형 (2)를 다음과 같이 확장하고자 한다[17]. i-번째 개체의 q-번째 (q = 1,2,…, Q)종단 변수에 대한 랜덤 효과 biq가 주어졌을 때, t 시점에서 i-번째 개체의 q-번째 종단 자료 값 yiq(t)에 대한 참 모형을 다음과 같이 가정하자.

yiq(t)=miq(t)+ϵiq(t)=βqxiq(t)+biqziq(t)+ϵiq(t).

여기서 xiq(t) 와ziq(t) 는 q-번째 종단 자료의 고정 효과와 랜덤 효과에 대한 디자인 벡터이다. 랜덤 효과 벡터biq ∼ N(0,Dq) 라고 가정하고, 측정 오차는ϵiq(t)N(0,σq2) 이라고 가정하자.biq 와 ∈iq(t) 는 서로 독립이라고 가정하자. 3절에서 정의한 모수 θ를 위 모형에 맞춰 업데이트하고, 가정 (A2)를 다음과 같이 확장하자.

 (A2) P(yi|bi;θ)=q=1QP(yiq|biq;θ)=q=1Qj=1niP(yiq(tij)|biq;θ)

여기서 yi=(yi1,yi2,...,yiQ),yiq=(yiq1,yiq2,...,yiqni),bi=(bi1,bi2,...,biQ)N(0,D) 이다. 또한 Q개의 종단 자료에 맞춰 생존 모형 (1)은 다음과 같이 표현된다.

hi(t|m~iq(t),wi)=h0(t)exp(γwi+q=1Qαqmiq(t)).

여기서 m~iq(t)={miq(s),0s<t} 이다. 따라서 관찰 자료에 대한 우도 함수는 3절과 유사한 방법으로 유도할 수 있다.

모형 진단

결합 모형에서 잔차는 종단 모형에 대한 잔차와 생존 모형에 대한 잔차로 구분된다. 종단 모형에 대한 잔차는 환자-중심(patient-specific) 잔차와 주변(marginal) 잔차로 나누어진다. 전자는 측정 오차에 대한 동질성 및 정규성 가정을 검정하는 데 사용하고, 후자는 고정 효과의 구조 및 환자 내 공분산 구조에 대한 가정을 검정하는 데 사용한다. 한편 생존 모형에 대한 잔차는 마팅게일(martingale) 잔차와 콕스-스닐(Cox-Snell) [18] 잔차가 있으며, 전자는mi(t) 의 함수 꼴에 대한 가정을 검정하는 데 사용하고, 후자는 모형의 적합성을 검정하는 데 사용한다.

Figure 3은 모형 (6)과 와이블 기저위험함수를 가진 모형 (7)의 결합 모형에 대한 환자-중심 잔차와 주변 잔차, 마팅게일 잔차, 콕스-스닐 잔차를 그린 그림이며, 이에 대한 R 코드는 다음과 같다.

Figure 3.

Subject-specific standardized residuals versus subject—specific fitted values (A), marginal standardized residuals versus marginal fitted values (B), martingale residuals versus subject-specific fitted values (C) of the log of serum bilirubin, and the Kaplan-Meier estimates of the survival function of the Cox-Snell residuals (red line: the survival function of the unit exponential distribution) (D) for the PBC data. PBC, primary biliary cirrhosis.

par(mfrow=c(2,2))

resSubY.pbc = residuals(wb.pbc, process = “Longitudinal”, type = “stand-Subject”)

fitSubY.pbc = fitted(wb.pbc, process = “Longitudinal”, type = “Subject”)

plot(fitSubY.pbc, resSubY.pbc, xlab = “Subject-specific fitted values”, ylab = “Subject-specific Residuals”) #환자-중심 잔차

lines(lowess(fitSubY.pbc, resSubY.pbc), lwd = 2, col =”red”)

abline(h = 0, lty = 3, col = “grey50”, lwd = 2)

resMarY.pbc = residuals(wb.pbc, process = “Longitudinal”, type = “stand-Marginal”)

fitMarY.pbc = fitted(wb.pbc, process = “Longitudinal”, type = “Marginal”)

plot(fitMarY.pbc, resMarY.pbc, xlab = “Marginal fitted values”, ylab = “Marginal Residuals”) #주변 잔차

lines(lowess(fitMarY.pbc, resMarY.pbc), lwd = 2, col =”red”)

abline(h = 0, lty = 3, col = “grey50”, lwd = 2)

martRes.pbc = residuals(wb.pbc, process = “Event”, type = “Martingale”)

mi.pbc = fitted(wb.pbc, process = “Longitudinal”, type = “EventTime”)

plot(mi.pbc, martRes.pbc, xlab = “Subject-specific fitted values”, ylab = “Martingale residuals”) #마팅게일 잔차

lines(lowess(mi.pbc, martRes.pbc), lwd = 2, col =”red”)

abline(h = 0, lty = 3, col = “grey50”, lwd = 2)

csRes.pbc = residuals(wb.pbc, process = “Event”, type = “CoxSnell”)

sfit.pbc = survfit(Surv(csRes.pbc, status2) ~ 1, data = pbc2.id)

plot(sfit.pbc, xlab = “Cox-Snell residuals”, ylab = “Survival probability”, conf.int = TRUE) #콕스-스닐 잔차

curve(exp(-x), from = 0, to = max(pbc2.id$years), add = TRUE, col = “red”, lwd = 2)

환자-중심 잔차나 주변 잔차에서 전반적인 추세선(평활 값: 빨간색)이 0으로부터 벗어나지 않았지만 ±2를 벗어나는 잔차는 다수 관측된다. 마팅게일 잔차의 추세선(평활 값: 빨간색)은 0에 가까웠지만 종단 자료 값이 크면 약간 벗어나는 경향이 있었고, 콕스-스닐 잔차(실선)의 생존 함수 추정량(파선: 95% CI)은 참조선(평균이 1인 지수분포: 빨간색)과 거의 같았다.

이벤트 확률 및 종단 자료 값의 예측

본 절에서는 추정된 모형을 사용하여 한 환자의 이벤트 발생 확률과 종단 자료 값을 동적으로 예측하고자 하며, 바이오마커가 환자들을 고위험 환자와 저위험 환자로 얼마나 잘 구분하는지를 평가할 수 있는 동적 ROC와 AUC 및 구분 지수(discrimination index)를 소개하고자 한다.

첫째, 새로운 환자 ℓ이t(>0) 시점까지 생존해 있다고 할 때u(> t) 시점 이후에 이벤트가 발생할 다음과 같은 확률을 예측하고자 한다.

(13) π(u|t)=P(Tu|T>t,y~(t),Dn;θ),t<u.

여기서 Y~(t)={y(s);0s<t} 이고, θ* 는 θ의 참값이다. 식 (13)은 근사 베이지안 방법[19]에 의해 다음과 같이 표현된다.

(14) π(u|t)=P(Tu|T>t,y~(t);θ)p(θ|Dn)dθ=P(Tu|T>t,b;θ)p(b|T>t,y~(t);θ)p(θ|Dn)dbbdθ=S(u|m~(u),b;θ)S(t|m~(t),b;θ)p(b|T>t,y~(t);θ)p(θ|Dn)dbdθ.

식 (14)로부터 몬테칼로(Monte Carlo) 방법을 사용하여πℓ(u|t) 를 추정하는 절차는 다음과 같다.

(Sl)θ(s)N(θ^,va^r(θ^))를 생성한다.

(S2)b(s){b|T>t,y~(t);θ(s)} 를 생성한다.

(S3)π(s)(u|t)=S(u|m~(u),b(s);θ(s))S(t|m~(t),b(s);θ(s)) 을 계산한다.

(S4) (S1)-(S3)을 S(s = 1,2,…,S) 번 반복한다.

따라서 πℓ(u|t) 를π^(u|t)=median{π(s)(u|t),s=1,2,...,S} 혹은 π^(u|t)=1Ss=1sπ(s)(u|t) 로 추정한다.

Figure 4는 모형 (6)과 와이블 기저위험함수를 가진 모형 (7)의 결합 모형에 대해 id=2인 환자의 각 t=0, 1, 4.9, 6.9년에서 이벤트(사망이나 이식) 발생 확률을 동적으로 예측한 그림이며, 이에 대한 R 코드는 다음과 같다.

Figure 4.

Dynamic survival probability at specified times of 0, 1, 5, and 7 (years) for patient id = 2 from the PBC data. PBC, primary biliary cirrhosis.

set.seed(1)

ND = pbc2[pbc2$id == 2, ]

inx = c(1, 3, 5, 7)

survPreds = vector(“list”, length(inx))

for (i in 1:length(inx)) survPreds[[i]] = survfitJM(wb.pbc, newdata = ND[1:inx[i], ])

par(mfrow = c(2, 2))

for (i in 1:length(inx)) plot(survPreds[[i]], estimator = “mean”, conf.int = TRUE, ylab = “”, ylab2 = “”, lwd = 2, xlab = “Times(years)”, include.y = TRUE, main = paste(“Follow-up time:”, round(survPreds[[i]]$last.time, 1))) #동적 이벤트 발생 확률

mtext(expression(log(“serum Bilirubin”)), side = 2, line = -2, outer = TRUE)

mtext(“Survival probability”, side = 4, line = -2, outer = TRUE)

id=2인 환자의 serBilir 값이 세 번째 관측 이후로 눈에 띄게 증가하므로 다섯 번째(Figure 4C)와 일곱 번째(Figure 4D) 관측 시점에서 구한 이벤트 발생 확률의 기울기가 급격하게 감소하고 있음을 알 수 있다.

둘째, 새로운 환자ℓ 이t(>0) 시점까지 생존해 있다고 할 때u(>t) 시점에서 종단 자료의 다음과 같은 기댓값을 예측하고자 한다.

(15) ω(u|t)=E(y(u)|T>t,y~(t),Dn;θ),t<u.

식 (15)는 근사 베이지안 방법[19]에 의해 다음과 같이 표현된다.

(16) ω(u|t)=E(y(u)|T>t,y~(t);θ)p(θ|Dn)dθ=E(y(u)|T>t,y~(t),b;θ)p(b|T>t,y~(t);θ)p(θ|Dn)dbbdθ=E(y(u)|b;θ)p(b|T>t,y~(t);θ)p(θ|Dn)dbdθ=(βx(u)+(p(b|T>t,y~(t);θ)db)z(u))p(θ|Dn)dθ.

식 (16)으로부터 몬테칼로 방법을 사용하여ω(u|t) 를 추정하는 절차는π(u|t) 를 예측하는 절차와 유사한데 (S1)과 (S2)는 동일하게 하고, (S3)만 다음과 같이 변형하면 된다.

(S3)ω(s)(u|t)=β^(s)x(u)+b(s)z(u)를 계산한다.

따라서 ω(u|t)를ω^(u|t)=median{ω(s)(u|t),s=1,2,...,S} 혹은ω^(u|t)=1Ss=1sω(s)(u|t) 로 추정한다.

Figure 5는 모형 (6)과 와이블 기저위험함수를 가진 모형 (7)의 결합 모형에 대해 id=2인 환자의 종단 자료 값의 기댓값을 각 t=0, 1, 4.9, 6.9년에서 예측한 것이며, 이에 대한 R 코드는 다음과 같다.

Figure 5.

Dynamic prediction of the log of the serum bilirubin at specified times of 0, 1, 5, and 7 (years) for patient id=2 from the PBC data. PBC, primary biliary cirrhosis.

set.seed(1)

ND = pbc2[pbc2$id == 2, ]

inx = c(1, 3, 5, 7)

longPreds = vector(“list”, length(inx))

for (i in 1:length(inx)) {

longPreds[[i]] = predict(wb.pbc, newdata = ND[1:inx[i], ], type = “Subject”, interval = “confidence”, returnData = TRUE)

longPreds[[i]]$FollowUp = round(max(ND[1:inx[i], “year”]), 1)

}

longPreds.all = do.call(rbind, longPreds)

longPreds.all$FollowUp = with(longPreds.all, factor(FollowUp, labels = paste(“Follow-up time:”, unique(FollowUp))))

xyplot(pred + low + upp ~ year | FollowUp, data = longPreds.all, lty = c(1, 2, 2), col = c(“red”, “black”, “black”), lwd = c(2, 1, 1), type = “l”, as.table = TRUE, xlab = “Time(years)”, ylab = expression(paste(“Predicted “, log(“serum Bilirubin”)))) #종단 자료의 동적 기댓값

t에 관계없이 시간이 지남에 따라 예측 구간은 폭이 넓어지고 비대칭적인 형태를 보였다. 따라서 가까운 미래 시점의 종단 자료 값을 예측하는 것은 신뢰할 수 있지만 현재 시점으로부터 멀리 떨어질수록 예측의 정확도가 떨어진다고 할 수 있다.

셋째, 새로운 환자ℓ 의y~(t) 가 주어졌을 때 구간(t,t+dt] 에서 이벤트의 발생 여부는 환자의 처치 여부 결정하는 데 매우 중요하다. 먼저 바이오마커에 대한 예측 규칙을 정해야 하는데 간단하고 가장 널리 사용되는 규칙은 시점 t에서 관측한 바이오마커의 값을 사용하여 다음과 같이 정의하는 (시간-가변) 규칙이다.

(17) "ytc"(단, c는 상수).

규칙 (17)은 “시점 t에서 바이오마커의 수준이 c 이상이면 구간(t,t+Δdt] 에서 이벤트가 발생할 가능성이 높은 환자로 분류한다.”는 의미이다. 따라서 규칙 (17)에 맞춰 동적 민감도와 특이도를 정의하면 각각 다음과 같다[20].

TPtΔt(c)=P(yi(t)c|Ti>t,Ti(t,t+Δt];θ),1FPtΔt(c)=P(yi(t)<c|Ti>t,Ti>t+Δt;θ).

민감도는 구간 (t,t+Δt]에서 이벤트가 발생할 환자를 바이오마커가 소위 시험군으로 분류할 확률을 의미하고, 특이도는(t + dt) 이후에 이벤트가 발생할 환자를 바이오마커가 소위 대조군으로 분류할 확률을 의미한다. 따라서 동적 ROC 곡선과 AUC를 정의하면 각각 다음과 같다[20].

ROCtΔt(p)=TPtΔt{[FPtΔt]1(p)},0<p<1AUCtΔt=01ROCtΔt(p)dp.

여기서 [FPtΔt]1(p)=infc{FPtΔt(c)p},cϵRy이다. 마지막으로 동적 구분 지수는 구간의 폭Δt 에만 의존하는데 정의하면 다음과 같다[21].

CΔt=0AUCtΔtu(t)dt.

여기서 u(t)는 가중값을 의미하며T* 의 주변분포로부터 다음과 같이 얻어진다.

u(t)=P(T>t)0P(T>t)dt.

몬테칼로 샘플을 사용하여 위에서 정의한 동적 민감도, 특이도, ROC 곡선, AUC 및 구분지수를 각각 추정할 수 있다[5].

Figure 6은 모형 (6)과 와이블 기저위험함수를 가진 모형 (7)의 결합 모형에 대해 id=313인 환자의 각 t=1, 3, 6, 9년에서 Δt=1, 2, 4년에 따라 동 적ROCtΔt 곡선을 그린 그림이다. Table 6은 Δt=1, 2, 4년에 따라 id=313인 환자의 각 t=1, 3, 6, 9년에서 동적AUCtΔt 와 동적 구분지수6) CΔt 를 추정한 것이다. 이에 대한 R 코드는 다음과 같다.

Figure 6.

Time-dependent ROC curves at specific times of 1, 3, 6, and 9 (years) under the simple prediction rule for patient id = 313 from the PBC data. ROC, receiver operating characteristic; PBC, primary biliary cirrhosis.

Estimated time-varying AUCs and discrimination indies of patient id=313 at four specific times with three options for Δt under the simple prediction rule

set.seed(1)

ND = data.frame(id = 313, drug = factor(“D-penicil”, levels = levels(pbc2$drug)), year = c(0, 0.01, seq(1, 14), max(pbc2$year)))

dt = c(1, 2, 4)

ROCs = vector(“list”, (nrow(ND) - 1))

for (i in seq_along(ROCs)) {

ROCs[[i]]= rocJM(wb.pbc, dt = dt, data = ND[seq_len(i+1),], directionSmaller = FALSE)

}

par(mfrow = c(2, 2))

for (i in c(2, 4, 7, 10)) {

plot(ROCs[[i]], legend = TRUE, lwd = 2) #동적 ROC 곡선

}

AUCs = sapply(ROCs, “[[“, “AUCs”)

AUCs[, c(2, 4, 7, 10)] #동적 AUC

sf = survfit(Surv(years, status2) ~ 1, data = pbc2.id)

times = ND$year[-1]

survs = summary(sf, times = times)$surv

a = head(times, -1)

b = tail(times, -1)

denom = sum((b - a) * (head(survs, -1) + tail(survs, -1)) / 2) #사다리꼴 법칙7)

Cinxs = c()

for (i in 1:length(dt)){

numer = sum((b - a) * (head(AUCs[i,] * survs,-1) + tail(AUCs[i,] * survs,-1)) / 2) #사다리꼴 법칙

Cinxs[i] = numer / denom

}

cbind(dt = dt, Cindex = Cinxs) #동적 구분 지수

t의 값에 관계없이 Δt의 값이 증가할수록 t와t + Δt 사이에서 이벤트가 발생할 환자와t + Δt 이후에 이벤트가 발생할 환자를 구분하는 바이오마커의 예측력이 증가한다고 할 수 있다. 다만 동적 ROC와 AUC 값은 t에 의존하므로Δt 의 값을 고정했을 때t 에 따른 별다른 추이는 찾아볼 수 없었다. 한편 구분지수의 값도Δt 의 값이 증가할수록 증가함을 알 수 있다.

맺음말

본 연구에서는 종단 위험인자의 값이 연속적으로 변할 뿐만 아니라 측정 오차를 포함하고 있을 때 널리 사용되고 결합 모형을 소개하고 R 패키지 JM을 사용하여 구현하는 방법을 소개하였다. 결합 모형 중에서 가장 널리 사용되는 현재 값 모형을 소개한 후, 이 모형을 여러 가지 유형의 자료로 확장하였다. 또한 모형 진단을 위한 여러 가지 잔차를 소개하였으며, 추정 모형의 응용으로 이벤트의 발생 확률 및 종단 위험인자의 값을 예측하는 방법과 종단 위험인자가 환자들을 고위험 그룹과 저위험 그룹으로 얼마나 잘 구분하는지를 평가할 수 있는 지수 등을 각각 소개하였다.

본 연구의 제한점은, 첫째, 종단 자료의 유형(범주형, 연속형 등)에 관계없이 다룰 수 있도록 종단 모형을 일반선형혼합모형으로 확장하는 문제와, 둘째, 준경쟁위험(semi-competing risks) 자료처럼 중간 이벤트가 있는 자료를 다룰 수 있도록 생존 모형을 다중상태 자료로 확장하는 문제와, 셋째, 기저위험함수와 위험인자를 승법적인 관계로 가정했는데 가법적인 관계로 가정하는 모형이나 가속수명시간 모형으로 확장하는 문제 등을 다루지 못한 점이다.

Notes

1)

R version 4.2.2까지는 기저위험함수에 대한 비모수적 추정 방법을 제공하였으나 그 이후 버전에서는 제공하고 있지 않다.

2)

JM package에는 공유 랜덤 효과 모형을 분석할 수 있는 옵션을 제공하고 있지 않아서 Table 2에서 제외하였다. SAS JMFit 매크로를 이용하면 이 모형을 분석할 수 있지만 z i(t)의 꼴에 따라 랜덤 효과가 이벤트에 미치는 효과를 해석하는 데 어려움이 따를 수 있다.

3)

Hepatomegaly는 내적 공변량이지만 상호작용 모형에서는 기저 값으로 고정하고 분석하였다.

4)

생존 모형 (1)에 시간에 따라 변하는 외적 공변량이 포함되면 JM package에서는 기저위험함수를 스플라인 방법으로만 추정할 수 있다.

5)

JM package는 층별 생존 모형에서 기저위험함수를 스플라인 방법으로만 추정할 수 있다.

6)

구분지수에서 정적분 구간의 상한은 이벤트 관측값의 최댓값인 14.1년으로 정의하였다.

7)

구분지수 계산에는 정적분이 필요한데 본 논문에서는 정적분을 사다리꼴 법칙으로 근사하였다. 즉. abf(x)dx(ba)f(a)+f(b)2

References

1. . Cox DR. Regression models and life-tables. J R Stat Soc Ser B Methodol 1972;34(2):187–220. DOI: 10.1007/978-1-4612-4380-9_37.
2. . Cox DR. Partial likelihood. Biometrika 1975;62(2):269–276. DOI: 10.2307/2335362.
3. . Andersen PK, Gill RD. Cox's regression model for counting processes: A large sample study. Ann Stat 1982;10(4):1100–1120. DOI: 10.1214/aos/1176345976.
4. . Prentice RL. Covariate measurement errors and parameter estimation in a failure time regression model. Biometrika 1982;69(2):331–342. DOI: 10.2307/2335407.
5. . Rizopoulos D. Joint models for longitudinal and time-to-event data: With applications in R Boca Raton, FL: CRC Press; 2012. DOI: 10.1201/b12208.
6. . Murtaugh PA, Dickson ER, Van Dam GM, Malinchoc M, Grambsch PM, Langworthy AL, et al. Primary biliary cirrhosis: Prediction of short-term survival based on repeated patient visits. Hepatology 1994;20(1 Pt 1):126–134. DOI: 10.1016/0270-9139(94)90144-9.
7. . Rizopoulos D. JM: An R package for the joint modelling of longitudinal and time-to-event data. J Stat Softw 2010;35(9):1–33. DOI: 10.18637/jss.v035.i09.
8. . Kaplan EL, Meier P. Nonparametric estimation from incomplete ob-servations. J Am Stat Assoc 1958;53(282):457–481. DOI: 10.2307/2281868.
9. . Kalbfleisch JD, Prentice RL. The statistical analysis of failure time data. 2nd ed.th ed. New York: John Wiley & Sons; 2002. DOI: 10.2307/1268689.
10. . Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc Ser B Methodol 1977;39(1):1–38. DOI: 10.1111/j.2517-6161.1977.tb01600.x.
11. . Rizopoulos D, Verbeke G, Lesaffre E. Fully exponential laplace approx-imations for the joint modelling of survival and longitudinal data. J R Stat Soc Ser B Methodol 2009;71(3):637–654. DOI: 10.1111/j.1467-9868.2008.00704.x.
12. . Hsieh F, Tseng YK, Wang JL. Joint modeling of survival and longitudinal data: Likelihood approach revisited. Biometrics 2006;62(4):1037–1043. DOI: 10.1111/j.1541-0420.2006.00570.x.
13. . Van der Vaart AW. Asymptotic statistics Cambridge: Cambridge University Press; 1998. DOI: 10.1017/CBO9780511802256.
14. . Akaike H. A new look at the statistical model identification. IEEE Trans-actions on Automatic Control 1974;19(6):716–723. DOI: 10.1109/TAC.1974.1100705.
15. . Schwarz G. Estimating the dimension of a model. Ann Stat 1978;6(2):461–464. DOI: 10.1214/aos/117634413.
16. . Król A, Ferrer L, Pignon JP, Proust-Lima C, Ducreux M, Bouché O, et al. Joint model for left-censored longitudinal data, recurrent events and terminal event: Predictive abilities of tumor burden for cancer evolution with application to the FFCD 2000-05 trial. Biometrics 2016;72(3):907–916. DOI: 10.1111/biom.12490.
17. . Hickey GL, Philipson P, Jorgensen A, Kolamunnage-Dona R. Joint models of longitudinal and time-to-event data with more than one event time outcome: A review. Int J Biostat 2018;14(1)DOI: 10.1515/ijb-2017-0047.
18. . Cox DR, Snell EJ. A general definition of residuals. J R Stat Soc Ser B Methodol 2018;30(2):248–265. DOI: 10.1111/j.2517-6161.1968.tb00724.x.
19. . Cox DR, Hinkley DV. Theoretical statistics Taylor & Francis; 1979. DOI: 10.1201/b14832.
20. . Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics 2005;61(1):92–105. DOI: 10.1111/j.0006-341X.2005.030814.x.
21. . Rizopoulos D. Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data. Biometrics 2011;67(3):819–829. DOI: 10.1111/j.1541-0420.2010.01546.x.

Article information Continued

Figure 1.

Kalan-Meier estimates of the survival function and the number of risk set at the selected times by drug for the PBC data. PBC, primary biliary cirrhosis.

Figure 2.

Subject-specific longitudinal profiles for sixteen patients from PBC data. PBC, primary biliary cirrhosis.

Table 1.

Estimates of regression parameters and their standard errors for longitudinal biomarker and/or event parts in the joint models with three different baseline hazard functions and in the extended Cox model

Model Longitudinal biomarker Event AIC BIC
Parameter Est SE p-value Parameter Est SE p-value
Joint Baseline hazard: Weibull
  Intercept 0.588 0.082 <0.001 DPCA 0.043 0.180 0.810 3,667.9 3,731.5
  DPCA -0.145 0.114 0.203 α 1.323 0.099 <0.001
  t 0.164 0.032 <0.001
  t2 0.002 0.004 0.498
  DPCA×t 0.005 0.044 0.908
  DPCA×t2 -0.002 0.004 0.686
Baseline hazard: piecewise exponential with 6 interior knots11
  Intercept 0.588 0.082 <0.001 DPCA 0.053 0.180 0.767 3,673.7 3,756.0
  DPCA -0.146 0.114 0.201 α 1.314 0.099 <0.001
  t 0.165 0.032 <0.001
  t2 0.002 0.004 0.516
  DPCA×t 0.004 0.044 0.920
  DPCA×t2 -0.002 0.004 0.686
Baseline hazard: spline with 5 interior knots2
  Intercept 0.586 0.082 <0.001 DPCA 0.049 0.181 0.788 3,673.6 3,763.4
  DPCA -0.143 0.114 0.209 α 1.310 0.098 <0.001
  t 0.166 0.032 <0.001
  t2 0.002 0.004 0.524
  DPCA×t 0.004 0.044 0.921
  DPCA×t2 -0.002 0.004 0.682
Separate Baseline hazard: nonparametric
DPCA 0.014 0.171 0.936 1,170.2 1,176.1
α 1.289 0.085 <0.001

SE, standard error; AIC, Akaike information criterion; BIC, Bayesian information criterion; DPCA, D-penicillamine.

1

Interior knots (yrs): 2.31, 4.12, 5.64, 6.83, 8.52, 10.31.

2

Interior knots (yrs): 2.53, 4.63, 6.30, 7.95, 9.91.

Table 2.

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint model under the Weibull baseline hazard function according to the association structure between the longitudinal and event outcomes

Association structure Longitudinal biomarker Event AIC BIC
Parameter Est SE p-value Parameter Est SE p-value
Current+ Slope Intercept 0.578 0.069 <0.001 DPCA 0.090 0.197 0.649 3,641.2 3,708.5
DPCA -0.152 0.093 0.103 α1 1.234 0.117 <0.001
t 0.184 0.031 <0.001 α2 3.250 0.682 <0.001
t2 0.002 0.003 0.629
DPCA×t -0.004 0.041 0.918
DPCA×t2 -0.001 0.004 0.826
Cumulative Intercept 0.586 0.083 <0.001 DPCA -0.109 0.173 0.527 3,824.4 3,888.1
DPCA -0.146 0.116 0.209 α 0.145 0.013 <0.001
t 0.162 0.032 <0.001
t2 0.001 0.003 0.668
DPCA×t 0.006 0.044 0.884
DPCA×t2 -0.002 0.004 0.605
Lagged Intercept 0.589 0.083 <0.001 DPCA 0.026 0.177 0.882 3,691.6 3,755.3
DPCA -0.145 0.116 0.212 α 1.268 0.094 <0.001
t 0.160 0.031 <0.001
t2 0.002 0.003 0.560
DPCA×t 0.007 0.043 0.874
DPCA×t2 -0.002 0.004 0.631
Interaction Intercept 0.588 0.082 0.000 DPCA 0.155 0.183 0.395 3,652.7 3,720.0
DPCA -0.146 0.115 0.206 α1 1.115 0.115 0.000
t 0.163 0.032 0.000 α2 0.358 0.089 0.000
t2 0.003 0.004 0.476
DPCA×t 0.004 0.044 0.933
DPCA×t2 -0.002 0.004 0.695

SE, standard error, AIC, Akaike information criterion; BIC, Bayesian information criterion; DPCA, D-penicillamine.

Table 3.

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint model with hepatomegaly (presence/absence) as an external time-dependent covariate

Longitudinal biomarker Event AIC BIC
Parameter Est SE p-value Parameter Est SE p-value
Intercept 0.589 0.082 <0.001 DPCA 0.049 0.180 0.785 3,674.5 3,768.1
DPCA -0.147 0.114 0.196 Presence 0.239 0.228 0.294
t 0.165 0.032 <0.001 α 1.272 0.105 <0.001
t2 0.002 0.004 0.518
DPCA×t 0.004 0.044 0.920
DPCA×t2 -0.002 0.004 0.678

SE, standard error; AIC, Akaike information criterion; BIC, Bayesian information criterion; DPCA, D-penicillamine.

Table 4.

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint models stratified according to the presence or absence of hepatomegaly

Longitudinal biomarker Event LRT AIC BIC
Parameter Est SE p-value Parameter Est SE p-value
Common covariate effect Intercept 0.586 0.082 <0.001 DPCA 0.113 0.183 0.535 0.002 3,665.9 3,789.4
DPCA -0.143 0.115 0.214 α 1.314 0.104 <0.001
t 0.164 0.032 <0.001
t2 0.002 0.004 0.506
DPCA×t 0.005 0.044 0.915
DPCA×t2 -0.002 0.004 0.670
Strata-specific effects Intercept 0.586 0.082 <0.001 DPCA -0.134 0.312 0.667 0.001 3,664.2 3,795.3
DPCA -0.143 0.116 0.214 DPCA×Presence 0.409 0.385 0.288
t 0.166 0.032 <0.001 α 1.103 0.140 <0.001
t2 0.002 0.004 0.543 α×Presence 0.426 0.207 0.040
DPCA×t 0.003 0.044 0.947
DPCA×t2 -0.002 0.004 0.705

SE, standard error; AIC, Akaike information criterion; BIC, Bayesian information criterion; DPCA, D-penicillamine.

Table 5.

Estimates of regression parameters and their standard errors for longitudinal biomarker and event parts in the joint model with competing risks (transplanted/dead)

Longitudinal biomarker Event AIC BIC
Parameter Est SE p-value Parameter Est SE p-value
Intercept 0.582 0.075 0.000 DPCA -0.425 0.383 0.267 3,930.8 4,061.8
DPCA -0.143 0.105 0.172 DPCA×Dead 0.486 0.415 0.242
t 0.168 0.031 0.000 α 1.189 0.197 0.000
t2 0.003 0.003 0.392 α×Dead 0.119 0.217 0.584
DPCA×t 0.002 0.043 0.960
DPCA×t2 -0.002 0.004 0.632

SE, standard error; AIC, Akaike information criterion; BIC, Bayesian information criterion; DPCA, D-penicillamine.

Figure 3.

Subject-specific standardized residuals versus subject—specific fitted values (A), marginal standardized residuals versus marginal fitted values (B), martingale residuals versus subject-specific fitted values (C) of the log of serum bilirubin, and the Kaplan-Meier estimates of the survival function of the Cox-Snell residuals (red line: the survival function of the unit exponential distribution) (D) for the PBC data. PBC, primary biliary cirrhosis.

Figure 4.

Dynamic survival probability at specified times of 0, 1, 5, and 7 (years) for patient id = 2 from the PBC data. PBC, primary biliary cirrhosis.

Figure 5.

Dynamic prediction of the log of the serum bilirubin at specified times of 0, 1, 5, and 7 (years) for patient id=2 from the PBC data. PBC, primary biliary cirrhosis.

Figure 6.

Time-dependent ROC curves at specific times of 1, 3, 6, and 9 (years) under the simple prediction rule for patient id = 313 from the PBC data. ROC, receiver operating characteristic; PBC, primary biliary cirrhosis.

Table 6.

Estimated time-varying AUCs and discrimination indies of patient id=313 at four specific times with three options for Δt under the simple prediction rule

Δt AUCtΔt CΔt
t
1 3 6 9
1 0.668 0.705 0.745 0.646 0.716
2 0.696 0.735 0.761 0.680 0.745
4 0.749 0.781 0.788 0.746 0.786

AUC, the area under the curve.