---
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)
```

