--- title: "Rak piersi a czas przeżycia bez nawrotu choroby" output: pdf_document: includes: in_header: naglowek.tex highlight: pygments toc: false number_sections: false keep_tex: true --- ```{r, echo=FALSE} library(knitr) opts_chunk$set(comment="", message=FALSE, warning=FALSE, tidy.opts=list(keep.blank.line=TRUE, width.cutoff=120),options(width=120), cache=TRUE,fig.align='center',fig.height=6, fig.width=10,fig.path='figure/beamer-',fig.show='hold',size='footnotesize', cache=FALSE) ``` \thispagestyle{fancy} ```{r, echo=FALSE} library(foreign) dane <- read.dta("gbcs_short.dta") library(survival) ``` Analizie poddano 686 pacjentek cierpiących na raka piersi. Za głowny cel postawiono pytanie, które zmienne mają wpływ na czas przeżycia bez nawrotu choroby. Zaproponowano model \text{proporcjonalnych hazardów (PH).} \textbf{Sprawdzenie założeń - krzywe przeżycia oraz ich transformacje.} \newline Dla każdej ze zmiennych (poza ciągłymi) sprawdzono spełnienie założeń poprzez narysowanie krzywych Kaplana-Meiera (Rysunek 1) oraz wykresów transformacji \text{\textsf{log(-log)}} krzywych przeżycia \text{(Rysunek 2).} ```{r, echo=FALSE, results='hide', cache=TRUE} pdf( file = "sc_meno.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~meno, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3) # model warstwowy dev.off() pdf( file = "sc_horm.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~horm, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3) dev.off() pdf( file = "sc_prog.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~prog, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3) dev.off() pdf( file = "sc_estr.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~estr, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3) dev.off() pdf( file = "sc_grade.pdf", onefile = FALSE) plot(survfit(Surv(rectime,censrec)~grade, data = dane), col=c("orange","purple", "green"), lty=c(1:3), lwd=3) dev.off() ``` ```{r, echo=FALSE, results='hide', cache=TRUE} pdf( file = "sc_meno_log.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~meno, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3, fun=function(x) log(-log(x)), log="x", firstx=1) dev.off() pdf( file = "sc_horm_log.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~horm, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3, fun=function(x) log(-log(x)), log="x", firstx=1) dev.off() pdf( file = "sc_prog_log.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~prog, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3, fun=function(x) log(-log(x)), log="x", firstx=1) dev.off() pdf( file = "sc_estr_log.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~estr, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3, fun=function(x) log(-log(x)), log="x", firstx=1) dev.off() pdf( file = "sc_grade_log.pdf", onefile=FALSE) plot(survfit(Surv(rectime,censrec)~grade, data = dane), col=c("orange","purple", "green"), lty=c(1:3), lwd=3, fun=function(x) log(-log(x)), log="x", firstx=1) dev.off() ``` \begin{figure}[hbt!] \vspace{-10pt} \begin{center} \subfigure[\textsf{meno}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_meno.pdf}} \subfigure[\textsf{horm}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_horm.pdf}} \subfigure[\textsf{prog}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_prog.pdf}} \subfigure[\textsf{estr}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_estr.pdf}} \subfigure[\textsf{grade}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_grade.pdf}} \end{center} \vspace{-20pt} \label{fig:sc} \caption{Krzywe przeżycia Kaplana-Meiera dla zmiennych dyskretnych.} \end{figure} \begin{figure}[hbt!] \vspace{-10pt} \begin{center} \subfigure[\textsf{meno}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_meno_log.pdf}} \subfigure[\textsf{horm}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_horm_log.pdf}} \subfigure[\textsf{prog}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_prog_log.pdf}} \subfigure[\textsf{estr}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_estr_log.pdf}} \subfigure[\textsf{grade}]{ \includegraphics[width=0.18\textwidth, height=1.15in]{sc_grade_log.pdf}} \end{center} \vspace{-20pt} \label{fig:sc} \caption{Transformacje \textsf{log(-log)} krzywych przeżycia Kaplana-Meiera.} \end{figure} Krzywe przeżycia dla rożnych poziomów zmiennej \texttt{meno} przecinają się, co może oznaczać niespełnianie założeń modelu PH. Dla pozostałych zmiennych krzywe nie przecinają się oraz wykresy transformacji możemy uznać za nieodstające od równoległych. Z tej racji zaproponowano model proporcjonalnych hazardów warstwowany względem zmiennej \texttt{meno}, gdyż nie ma podstaw by nie twierdzić, że założenia modelu ph nie są spełnione. \begin{wraptable}{r}{6.8cm} \vspace{-20pt} \caption{ Wyniki testów logrank. } \begin{tabular}{lrrrrrr} \toprule% ```{r, results='asis',echo=FALSE} cat(c("\\ ", "$\\chi^2$", "St. swobody","p-wartość"), sep="&") cat(c("\\\\", "\\toprule ")) cat(c("meno","0.3","1","0.597"), sep="&") cat(c("\\\\ ")) cat(c("horm","8.6","1","0.003"), sep="&") cat(c("\\\\ ")) cat(c("prog","49.3","1","0.000"), sep="&") cat(c("\\\\ ")) cat(c("estr","14.3","1","0.000"), sep="&") cat(c("\\\\ ")) cat(c("grade","44.5","2","0.000"), sep="&") cat(c("\\\\ ","\\bottomrule")) ``` \end{tabular} \vspace{-7.5pt} \end{wraptable} W celu potwierdzenia wniosków płynących z graficznej prezentacji krzywych przeżycia, przeprowadzono formalny test logrank dla każdej zmiennej, którego podsumowanie widać w Tabeli 1. Dla zmiennej \textsf{grade} przeprowadzono test logrank dla trendu z racji na naturalne uporządkowanie poziomów. Hipoteza zerowa zakłada, że krzywe przeżycia dla zmiennych nie różnią się istotnie statystycznie. Wartości krytyczne dla zmiennych \textsf{horm}, \textsf{prog}, \textsf{estr} i \textsf{grade}, które są mniejsze od zakładanego poziomu istotności \(\alpha=0.05\) dają statystycznie istotne podstawy do odrzucenia hipotez zerowych w testach logrank oraz do przyjęcia hipotez alternatywnych o tym, że krzywe przeżycia dla tych zmiennych różnią się. Dla zmiennej \textsf{meno} p-wartość \(0.597\) nie daje podstaw do odrzucenia hipotezy zerowej (dla przyjętego poziomu istnotności \(\alpha=0.05\), mówiącej o równości krzywych przeżycia dla różnych poziomów tej zmiennej. ```{r, echo=FALSE, eval=FALSE} model <- coxph(Surv(rectime,censrec)~horm+prog+estr+as.factor(grade)+strata(meno)+log(size)+log(nodes), data = dane) cox.zph(model, transform = "identity") -> x knitr::kable( round(x$table, digits = 3), format="latex") ``` \newpage \textbf{Przekształcenia zmiennych ciągłych.} \newline W celu sprawdzenia, czy zmienne ciągłe nie powinny zostać przekształcone przed wprowadzeniem ich do modelu, wykonano wykresy reszt martyngałowych pustego modelu od każdej z tych zmiennych. Otrzymane wykresy sugerują przekształcenie zmiennych poprzez logarytm, więc dla potwierdzenia tych przypuszczeń narysowano dodatkowo wykresy reszt od logarytmów zmiennych (Rysunek 3) i zaobserwowano, że wykresy są teraz bliższe liniowym. Wprowadzono więc do modelu zmienne ciągłe przekształcone logarytmem. ```{r, echo=FALSE, results='hide', cache=TRUE} pdf( file = "log_size.pdf", onefile=FALSE) plot(log(dane$size), resid(coxph(Surv(rectime,censrec)~1, data = dane))) lines(lowess(log(dane$size), resid(coxph(Surv(rectime,censrec)~1, data = dane)), iter = 0, f = 0.6)) dev.off() pdf( file = "log_nodes.pdf", onefile=FALSE) plot(log(dane$nodes), resid(coxph(Surv(rectime,censrec)~1, data = dane))) lines(lowess(log(dane$nodes), resid(coxph(Surv(rectime,censrec)~1, data = dane)), iter = 0, f = 0.6)) dev.off() ``` \begin{figure}[hbt!] \vspace{-10pt} \begin{center} \subfigure[\textsf{nodes}]{ \includegraphics[width=0.32\textwidth, height=1.3in]{log_nodes.pdf}} \subfigure[\textsf{size}]{ \includegraphics[width=0.32\textwidth, height=1.3in]{log_size.pdf}} \end{center} \vspace{-20pt} \label{fig:sc} \caption{Wykresy reszt martyngałowych od logarytmów zmiennych ciągłych.} \end{figure} \begin{wraptable}{r}{6.5cm} \vspace{-12pt} \caption{ Wyniki testu Schoenfelda. } \begin{tabular}{lrrr} \toprule% & $\rho$ & $\chi^2$ & p-wartość\\ \toprule -horm & -0.009 & 0.024 & 0.876\\ -prog & 0.036 & 0.369 & 0.543\\ -estr & 0.071 & 1.458 & 0.227\\ -as.factor(grade)2 & -0.052 & 0.811 & 0.368\\ -as.factor(grade)3 & -0.089 & 2.347 & 0.126\\ -log(size) & 0.011 & 0.037 & 0.848\\ -log(nodes) & -0.051 & 0.879 & 0.349\\ -GLOBAL & NA & 11.322 & 0.125\\ \bottomrule \end{tabular} \vspace{-7.5pt} \end{wraptable} \textbf{Sprawdzenie założeń - formalny test Schoenfelda.} \newline W celu formalnego sprawdzenia czy współczynniki w modelu są stałe w czasie przeprowadzono test Schoenfelda, którego globalna p-wartość oraz pojedyncze p-wartości dla zmiennych dyskretnych oraz ciągłych przekształconych przez logarytm są większe \text{od zakładanego} poziomu istotności $\alpha=0.05$ co nie daje podstaw do odrzucenia hipotezy zerowej w tym teście, mówiącej o stałości współczynników w czasie. Uwagę zwrócono również na wykresy skalowanych reszt Schoenfelda dla każdej zmiennej w modelu od czasu i dopasowano do nich krzywe. Dołączone granice dla 95\% obszarów ufności sugerują, że nie ma podstaw do odrzucenia hipotez, że krzywe nie różnią się od horyzontalnych. ```{r, echo=FALSE, results='hide', eval=FALSE} pdf( file = "skal_res_shen.pdf", onefile=FALSE) par(mfrow=c(2,4)) plot(cox.zph(model, transform = "identity"), df=4,nsmo=10, se=TRUE) dev.off() ``` \begin{figure}[hbt!] \vspace{-10pt} \begin{center} \includegraphics[width=0.85\textwidth, height=3in]{skal_res_shen.pdf} \caption{Wykresy skalowanych reszt Schoenfelda.} \end{center} \end{figure} Zarówno wykresy skalowanych reszt Schoenfelda, jak i wyniki testu nie sugerują odstępstw od założenia PH. \textbf{Dopasowanie modelu proporcjonalnych hazardów.} \newline ```{r} model <- coxph(Surv(rectime,censrec)~horm+prog+estr+as.factor(grade)+strata(meno)+log(size)+log(nodes), data = dane) summary(model) ``` W stworzonym modelu zmiennymi istotnymi statystycznie na poziomie istotności \(0.05\) są: \texttt{horm}, \texttt{prog}, \texttt{as.factor(grade)2}, log(\texttt{nodes}). Współczynniki modelu przy tych zmiennych wynoszą odpowiednio: $-0.40, -0.67, 0.50, 0.47$. Oznacza to, że jeśli została użyta terapia hormonalna, to hazard zgonu lub nawrotu choroby zminiejsza się $\exp(-0.40)=0.66$ raza w stosunku do nieużycia terapii hormonalnej (gdy wszystkie inne zmienne są takie same). Jeśli wskaźnik receptorów progesteronu zmienia się z ujemnego na dodatni, to hazard zmieni się o $\exp(-0.67)=0.5$ raza, natomiast gdy logarytm liczby węzłów chłonnych z przerzutami nowotworu wzrośnie o 1, to hazard zwiększy się $\exp(0.47)=1.61$ raza. Zmiana z grade (stopień zróżnicowania komórek nowotworu (1–wysoki, 2–średni) z 1 do 2 zwiększa hazard $\exp(0.5)=1.65$ raza, gdy wszystkie inne zmienne pozostaną na tym samym poziomie. Testy Likelihood, Wald i Score dają p-wartość mniejszą od 0.05 (przyjętego poziomu istotności), co wskazuje, że odrzucamy hipotezę zerową na rzecz hipotezy alternatywnej, zatem dopasowany model jest istotnie lepszy od modelu pustego. \newpage \textbf{Sprawdzenie dopasowania modelu proporcjonalnych hazardów.} \newline W celu sprawdzenia dopasowania modelu wygenerowano wykresy reszt dewiancji/martyngałowych od liniowej kombinacji zmiennych/indeksów. Na niektórych wykresach można wskazać wyraźnie odstające reszty. Wykres reszt dewiancji od indeksów można uznać za symetryczny, pozostałe niekoniecznie, \text{co sugeruje}, że dopasowany model nie jest perfekcyjny. ```{r, echo=FALSE, results='hide', eval=FALSE} library(ggplot2) library(ggthemes) pdf( file = "res_ind_dev.pdf", onefile=FALSE) qplot(1:686,residuals(model, type="deviance"))+xlab("index")+theme_tufte(base_size=20)+geom_hline()+ geom_hline(yintercept=0, col ="orange", size = 3) dev.off() pdf( file = "res_ind_mart.pdf", onefile=FALSE) qplot(1:686,residuals(model))+xlab("index")+theme_tufte(base_size=20)+geom_hline()+ geom_hline(yintercept=0, col ="orange", size = 3) dev.off() pdf( file = "res_lp_dev.pdf", onefile=FALSE) qplot(predict(model, type="lp"),residuals(model, type="deviance" ))+theme_tufte(base_size=20)+geom_hline()+ geom_hline(yintercept=0, col ="orange", size = 3) dev.off() pdf( file = "res_lp_mart.pdf", onefile=FALSE) qplot(predict(model, type="lp"),residuals(model))+theme_tufte(base_size=20)+geom_hline()+ geom_hline(yintercept=0, col ="orange", size = 3) dev.off() ``` \begin{figure}[hbt!] \vspace{-10pt} \begin{center} \subfigure[Reszty dewiancji a indeks.]{ \includegraphics[width=0.48\textwidth, height=2in]{res_ind_dev.pdf}} \subfigure[Reszty martyngałowe a indeks.]{ \includegraphics[width=0.48\textwidth, height=2in]{res_ind_mart.pdf}} \subfigure[Reszty dewiancji a liniowe kombinacje.]{ \includegraphics[width=0.48\textwidth, height=2in]{res_lp_dev.pdf}} \subfigure[Reszty martyngałowe a liniowe kombinacje.]{ \includegraphics[width=0.48\textwidth, height=2in]{res_lp_mart.pdf}} \end{center} \vspace{-20pt} \label{fig:sc} \caption{Wykresy reszt martyngałowych i dewiancji.} \end{figure} \newpage \textbf{Kody.} \newline ```{r, eval=FALSE} library(foreign) dane <- read.dta("gbcs_short.dta") library(survival) plot(survfit(Surv(rectime,censrec)~meno, data = dane), col=c("orange","purple"), lty=c(1:2), lwd=3) #... plot(log(dane$size), resid(coxph(Surv(rectime,censrec)~1, data = dane))) lines(lowess(log(dane$size), resid(coxph(Surv(rectime,censrec)~1, data = dane)), iter = 0, f = 0.6)) #... par(mfrow=c(2,3)) plot(cox.zph(model, transform = "identity"), df=4,nsmo=10, se=TRUE) model <- coxph(Surv(rectime,censrec)~horm+prog+estr+grade+strata(meno)+log(size)+log(nodes), data = dane) summary(model) #... qplot(1:686,residuals(model, type="deviance"))+xlab("index")+theme_tufte(base_size=20)+geom_hline()+ geom_hline(yintercept=0, col ="orange", size = 3) ```