한동안 휴지기를 가졌던 R 통계 포스팅을 다시 진행하겠습니다. 앞서 Cox 비례위험모형의 적합도를 판정하기 위해 Schoenfeld 잔차를 분석했습니다. 하지만 모형 적합도를 판정하는 잔차는 여러 가지가 있습니다. 마팅게일 잔차 (martingale residuals) 역시 그 중 하나로 만약 모델이 적합할 경우 잔차들의 합이 0에 근접한 모습을 보여줍니다.
여기서는 새로운 데이터인 lung 데이터를 기준으로 알아보겠습니다. 이 데이터는 survival 패키지의 기본 데이터 중 하나입니다. North Central Cancer Treatment Group의 폐암 환자의 생존 데이터를 공개한 것입니다. 이벤트는 status에 기록되어 있습니다.
inst: Institution code
time: Survival time in days
status: censoring status 1=censored, 2=dead
age: Age in years
sex: Male=1 Female=2
ph.ecog: ECOG performance score (0=good 5=dead)
ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
pat.karno: Karnofsky performance score as rated by patient
meal.cal: Calories consumed at meals
wt.loss: Weight loss in last six months
이 데이터의 구조는 str(lung), 설명은 ?lung에서 볼 수 있습니다. 여기서는 survival 외에 survminer 패키지를 사용할 것입니다.
library("survival")
library("survminer")
data(lung)
str(lung)
?lung
우선 나이와 생존율의 관계를 알아보겠습니다. 상식적으로 생각할 때 나이가 많은 폐암 환자일수록 사망 가능성이 높고 생존 기간이 짧을 것입니다. 여기서는 간단한 Cox 비례위험모형 결과와 마팅게일 잔차 그래프를 같이 보겠습니다.
lung2<-coxph data=" " lung="" span="" status="=2)~age," ties="efron" time="" urv="">
summary(lung2)
lung$resid_mart <- lung2="" residuals="" span="" type="martingale">
ggplot(data = lung, mapping = aes(x = age, y = resid_mart)) +
geom_point() +
geom_smooth() +
labs(title = "age") +
theme_bw() + theme(legend.key = element_blank())
> lung2<-coxph data=" " lung="" span="" status="=2)~age," ties="efron" time="" urv="">
> summary(lung2)
Call:
coxph(formula = Surv(time, status == 2) ~ age, data = lung, ties = "efron")
n= 228, number of events= 165
coef exp(coef) se(coef) z Pr(>|z|)
age 0.018720 1.018897 0.009199 2.035 0.0419 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.019 0.9815 1.001 1.037
Concordance= 0.55 (se = 0.025 )
Rsquare= 0.018 (max possible= 0.999 )
Likelihood ratio test= 4.24 on 1 df, p=0.04
Wald test = 4.14 on 1 df, p=0.04
Score (logrank) test = 4.15 on 1 df, p=0.04
결과는 HR 1.019 (1.001 - 1.037)로 나이가 증가함에 따라 사망 가능성이 증가했습니다.
마팅게일 잔차의 평균은 거의 0에 수렴하고 있어 모델은 적합한 것으로 평가할 수 있습니다. 만약 적합하지 않은 경우 변수 변환을 시도할 수 있는데, 로그 함수를 적용하거나 제곱을 하는 등 몇 가지 전통적인 변형법을 사용할 수 있습니다. ggcoxfunctional 함수를 통해 한번에 그래프를 그려 비교할 수 있습니다.
ggcoxfunctional(Surv(time, status) ~ age + log(age) + sqrt(age), data = lung)
여기서는 변수 변환은 필요하지 않은 것으로 보입니다.
또 다른 잔차 확인법이 콕스 - 스넬 잔차 (Cox-Snell residuals)로 일종의 표준화 잔차입니다. 이 잔차가 적합하다면 관측치와 예측치가 같은 경향을 보여 점들이 45도 각도로 배열됩니다. 여기서 많이 벗어나면 문제가 있는 것이죠. 이를 구현하기 위한 코드가 좀 복잡합니다.
#Cox-Snell residual
lung2<-coxph data=" " lung="" span="" status="=2)~age," ties="efron" time="" urv="">
summary(lung2)
lung$resid_mart <- lung2="" residuals="" span="" type="martingale">
## Cox-Snell residuals
lung$resid_coxsnell <- -="" 1="" lung="" resid_mart="" span="" status="">
## Fit model on Cox-Snell residuals (Approximately Expo(1) distributed under correct model)
fit_coxsnell <- 1="" coxph="" formula="Surv(resid_coxsnell," span="" status="">
data = lung,
ties = c("efron","breslow","exact")[1])
## Nelson-Aalen estimator for baseline hazard (all covariates zero)
df_base_haz <- basehaz="" centered="FALSE)</span" fit_coxsnell="">
## Plot
ggplot(data = df_base_haz, mapping = aes(x = time, y = hazard)) +
geom_point() +
scale_x_continuous(limit = c(0,3.5)) +
scale_y_continuous(limit = c(0,3.5)) +
labs(x = "Cox-Snell residuals as pseudo observed times",
y = "Estimated cumulative hazard at pseudo observed times") +
theme_bw() + theme(legend.key = element_blank())
그래프는 거의 45도 각도로 올라가는데 이벤트가 2로 코딩되어 있어 0부터 시작하지 않는 모습을 보입니다. 아무튼 큰 문제 없이 해석은 가능합니다.
참고
댓글
댓글 쓰기