Co to jest cenzurowanie?

Występują zatem sytuacje, dla których nie obserwujemy dokładnego czasu zdarzenia \(T\).

Mamy wtedy do czynienia z cenzurowaniem.

Możemy obserwować różne typy cenzurowania, takie jak:

• prawostronne,
• lewostronne,
• przedziałowe.

Zapis formalny

• \(\tilde{T} \geq 0\) - zmienna losowa opisująca czas do wystąpienia zdarzenia
• \(C \geq 0\) - zmienna losowa opisująca czas cenzurowania

Dla cenzurowania prawostronnego obserwujemy czas do wystąpenia zdarzenia: \[ T = min(\tilde{T}, C) \]

W praktyce obserwujemy pary \((T,\delta)\), gdzie \(\delta=\mathbb{I}(\tilde{T}, C)\) jest wskaźnikiem cenzurowania.

Zapis formalny

Funkcja przeżycia: \(S(t) = \mathbb{P}(\tilde{T} \geq t)\). Uwaga! \(S(0)=1\)

Funkcja hazardu: \(\lambda(t) = \lim\limits_{h \rightarrow 0^+} \frac{\mathbb{P}(t \leq \tilde{T} < t+h | \tilde{T} \geq t)}{h}\).

Zależność występująca pomiędzypowyższymi funkcjami: \[ S(t) = e^{-\int\limits_0^t \lambda(u)du} \]

Model Coxa

Załóżmy, że mamy wektor zmiennych niezależnych: \[ \textbf{Z}(t)'=(Z_1(t),Z_2(t),\ldots,Z_p(t)) \] Model Coxa zakłada następującą postać funkcji hazardu: \[ \lambda(t,\textbf{Z}(t))=\lambda_0(t)\cdot e^{\beta'\textbf{Z}(t)} \]

Dla \(\textbf{Z}(t)=\textbf{Z}\) otrzymujemy model proporcjonalnych hazardów.

Model proporcjonalnych hazardów

Interpretacja założenia: Niech \(\textbf{Z}\) będzie zmienną opisującą wielkość guza:

• \(\textbf{Z} = 0\) ozn. mały guz (I),
• \(\textbf{Z} = 1\) ozn. średni guz (II),
• \(\textbf{Z} = 2\) ozn. duży guz (III).

Wówczas funkcja hazardu dla poszczególnych rozmiarów guza wynosi:

• mały guz: \(\lambda_I(t) = \lambda_0(t) \cdot e^{0 \cdot \beta} = \lambda_0(t)\)
• średni guz: \(\lambda_{II}(t) = \lambda_0(t) \cdot e^{1 \cdot \beta} = \lambda_0(t) \cdot e^{\beta}\)
• duży guz: \(\lambda_{III}(t) = \lambda_0(t) \cdot e^{2 \cdot \beta} = \lambda_0(t) \cdot e^{2\beta}\)

Model zakłada, że wzrost \(\textbf{Z}\) o 1 zwiększa hazard o \(e^\beta\) razy.

Jak sprawdzić założenie modelu?

Wykres skalowanych reszt Schoenfelda dla pewnej funkcji czasu \(g(t)\).

Dla czasów zdarzeń \(t_{(k)}\): \[ \hat{r}_{Sk} = \sum\limits_{i: t_{(i)}=t_{(k)}} \big\{ Z_{(k)}(t_{(k)}) - \overline{Z}(\hat{\beta},t_{(k)}) \big\}^T \] \[ \hat{r}_{Sk}^* = \hat{r}_{Sk} Var^{-1}(\hat{\beta},t_{(k)}) \]

model_PH <- coxph(Surv(futime,fustat)~., data=ovarian) 
model_PH_test<-cox.zph(model_PH,transform = "identity")
plot(model_PH_test, df=3, nsmo=10, se=TRUE, var=1)
abline(model_PH$coeff[1], 0, lty=5)

Postać funkcyjna zmiennych niezależnych

Reszty martyngałowe vis ciągła zmienna niezależna.

\[ M(t) = N(t) - \int\limits_0^t Y(u) \lambda_0 (u) e^{\beta^t Z(u)} du \] gdzie \(N(t)\) to procez zliczający zdarzenia, \(Y(t)\) to proces zbioru ekspozycji na ryzyko zdarzenia.

Postać funkcyjna zmiennych niezależnych

Jak znaleźć odpowiednie przekształcenie?

• dopasowujemy model “pusty”
• rysujemy wykres reszt martyngałowych dla danej zmiennej
• wygładzamy wykres dla ustalenia funkcji

model_PH <- coxph(Surv(futime,fustat)~1, data=ovarian) 
res_mart <- resid(model_PH)
plot(ovarian$age, res_mart)
lines(lowess(ovarian$age, res_mart, iter=0, f=0.6))

Dopasowanie modelu

Wykres reszt martyngałowych (lub dewiancji) vis liniowa kombinacja zmiennych.

Reszty oparte na dewiancji dla zmiennych niezależnych stałych w czasie: \[ \hat{r}_{Di} = \text{sign}(\hat{M}_i) \sqrt{-\hat{M}_i - \delta_i \log(\delta_i - \hat{M}_i)} \]

Uwaga! Z reguły reszty są symetrczne i mieszczą się w przedziale \((-2,2)\).

model_PH <- coxph(Surv(survtime, survind) ~ mutation + expression, data=nsclc) 
res_dev <- residuals(model_PH,type="deviance") 
fit_val <- predict(model_PH,type="lp") 
plot(fit_val, res_dev) 

Obserwacje wpływowe

Reszty score vis zmienna niezależna.

Reszta “score” dla obserwacji i i zmiennej j: \[ \hat{r}_{U_{ij}} = U_{ij} (\hat{\beta},\infty), \qquad U_{ijk}(\beta) = \int\limits_{t_{(k-1)}}^{t_k} \big\{ Z_{ij}(u) - \overline{Z}_j(\beta,u) \big\} dM_i(u) \]

model_PH <- coxph(Surv(futime,fustat)~., data=ovarian) 
score <- resid(model_PH, type="score")
plot(ovarian$age, score[,1], ylab="Score Residuals", xlab="Age")

Konkurencyjne ryzyka a reszty - pakiet survival

Postać danych

Dane przekształcone

Model PH dal funkcji hazardu

coxph_aids <- coxph(Surv(time, status == 1) ~ ccr5, data = aidssi)
coxph_si <- coxph(Surv(time, status == 2) ~ ccr5, data = aidssi)

Reszty Schoenfelda

model_PH_test<-cox.zph(coxph_aids,transform = "identity")
plot(model_PH_test, df=3, nsmo=10, se=TRUE, var=1)
abline(coxph_aids$coeff[1], 0, lty=5)

Reszty dewiancji

res_dev <- residuals(coxph_aids,type="deviance") 
fit_val <- predict(coxph_aids,type="lp") 
plot(fit_val, res_dev)

Konkurencyjne ryzyka a reszty - pakiet mstate

library(mstate)
tmat <- transMat(list(c(2,3,5,6), c(4,5,6), c(4,5,6), c(5,6), c(), c()),
  names=c("Tx", "Rec", "AE", "Rec+AE", "Rel", "Death"))
msebmt <- msprep(
  data = ebmt4, trans = tmat, 
  time = c(NA, "rec", "ae", "recae", "rel", "srv"),
status = c(NA, "rec.s", "ae.s", "recae.s", "rel.s", "srv.s"),
  keep = c("match", "proph", "year", "agecl"))
msebmt2 <- expand.covs(msebmt, longnames = FALSE, c("match", "proph", "year", "agecl"))
msebmt2[, c("Tstart", "Tstop", "time")] <- msebmt2[, c("Tstart", "Tstop", "time")]/365.25
knitr::kable(msebmt2[msebmt2$id == 1, (1:9)])

model <- coxph(Surv(Tstart, Tstop, status)~strata(trans), data = msebmt2, method = "breslow")
model_fit <- msfit(object = model, vartype = "greenwood", trans = tmat)

Podsumowanie

• Wszystkie funkcje dzięki, którym wyliczamy reszty wymagają jako argumenty obiektu coxph.
• Dla reszt martyngałowych dodatkowo wymagany jest model pusty, co w przypadku modelu wielostanowego dla funkcji coxph jest niemożliwe.