#=========================================
#     ANÁLISE DE SOBREVIVÊNCIA
#     Prof.: Enrico A. Colosimo
#=========================================
#
#=========================================
#  Códigos em R utilizados nas análises
#=========================================
#
#
# Dados sem Censura
# Leucemia (Feigl e Zelen, 1965)
# Livro: Cox e Snell (1981), p. 148
# Y:tempo do diagnóstico até a morte (em semanas)
# X: log10(contagem de células brancas no diagnóstico)
# n=17
# Objetivo: descrever a relação entre Y e X
#=====================================================
#
 y<-c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
 x<-c(3.36,2.88,3.63,3.41,3.78,4.02,4,4.23,3.73,3.85,
 3.97,4.51,4.54,5,5,4.72,5)
 plot(x,y,xlab="log10 leucócitos", ylab="tempo até a morte")
 summary(y)
 summary(x)
# lines(lowess(x,y),col="red",lwd=2)  
# Dicotomizando x: teste t (4 é o limite superior de normalidade)
 x1<-ifelse(x>4,0,1)
 x1<-factor(x1)
 levels(x1)<-c(">4","<4")
 boxplot(y~x1,ylab="tempo até a morte")
 bartlett.test(y,x1) # teste de Bartlett
 var.test(y~x1) # teste F
 t.test(y~x1,var.equal=T) # teste t
 t.test(y~x1) # teste t de Welch
 out<-aov(y~x1)
 summary(out)
 qqnorm(out$res,ylab="Residuos", main=NULL)
 qqline(out$res)
 title("Gráfico Normal de Probabilidade dos Resíduos")
#
# O teste de Shapiro-Wilks verifica a normalidade.
#
 shapiro.test(out$res)
#
# X continuo
#
  plot(x,y,xlab="log10 leucócitos", ylab="tempo até a morte")
#
 a<-glm(y~x)
 plot(a) 
 summary(a)
# Teste de Normalidade
 shapiro.test(a$resid)
#
#  Um teste simples para verificar variância não-constante
#   Este teste não é completamente correto pois algum peso
#   deve ser utilizado e os graus de liberdade ajustados
#   Veja Faraway (2004)
#
 summary(lm(abs(residuals(a))~fitted(a)))
#  valor-p=0,357 para a hipótese de homocedasticidade
#
 plot(x,y,xlab="log10 leucócitos", ylab="tempo até a morte")
 abline(coef(a),lty=5)
 lines(lowess(x,y),col="red",lwd=2)
#
#============================================================
#
#  Modelagem Alternativa
#
#   Componente Sistemático: log E(Y) = modelo linear
#   Implica em função de ligação logaritmo
#
  a1<-glm(y~x,family=gaussian(link="log"))
  summary(a1)
  plot(a1)
  x1<-seq(2.5,5.5,0.1)
  yp1<-exp(a2$coefficients[1]+(a2$coefficients[2]*x1))
  plot(x,y,xlab="log10 leucócitos", ylab="tempo até a morte")
  abline(coef(a),lty=5,col=2)
  lines(x1,yp1,col=4)
#
#==============================================================
#
#  Ajuste do Livro de Cox e Snell (1981)
#
#   Y = exp(modelo linear)*erro
#       erro: distribuição exponencial
#
  y1<-log(y)
  plot(x,y1,ylab="log Y",xlab="log10 leucócitos")
  lines(lowess(x,y1),col="red",lwd=2)
#
  require(survival)
  cens<-rep(1,17)
  a2<-survreg(Surv(y,cens)~x,dist='exponential')
  summary(a2)
#
  yp2<-exp(a2$coefficients[1]+(a2$coefficients[2]*x1))
  plot(x,y,xlab="log10 leucócitos", ylab="tempo até a morte")
  abline(coef(a),col=1)
  lines(x1,yp1,col=2)
  lines(x1,yp2,col=4)
  legend(4.2, 150, legend = c("M1: linear-normal","M2: log-linear-normal", "M3: exponencial"), 
  lty=1,col = c(1,2,4), bty = "n", cex = 1)
#
# AIC: 1- Linear-normal: 178.3
       2- Log-normal: 178.4
#      3- Exponencial: -2*loglik + 2*2 = 2*83.9 + 4 = 171.8
#      4- Log-normal: -2*loglik + 2*3 = 2*83.8 + 6 = 173.6
#
# BIC: 1- (linear-normal): -2* (-86.15) + 3*log(17) = 180.8
#      2- (log-normal): -2* (-86.2) + 3*log(17) = 180.9
#      3- (exponencial): -2*loglik + plog(n) = -2*(-83.9) + 2*log(17) = 173,47
#      2- (log-normal): -2* (-83.8) + 3*log(17) = 176,10
#      3- (gama): + 167,97 + 3*log(17) = 176,47
#     ***  A favor do modelo exponencial ****
#
# Refazendo análise de Cox-Snell (1981, p.148)
#
  x1<-x-mean(x)
  a3<-survreg(Surv(y,cens)~x1,dist='exponential')
  summary(a3)
  exp(a3$coefficients)
#
# Respondendo às perguntas
#  1- Intervalo de Confiança para E[Y/x]
#      log E^[Y/x] = \beta_0^ + \beta_1^ x
   x1<-seq(2.5,5.5,0.1)
   log_yme<- a2$coefficients[1] + (x1*a2$coefficients[2])
   var_log_yme<- a2$var[1,1] + ((x1^2)*a2$var[2,2]) + (2*x1*a2$var[1,2])
   ICymed_sup<- exp(log_yme + 1.96*var_log_yme)
   ICymed_inf<- exp(log_yme - 1.96*var_log_yme)
   plot(x1,exp(log_yme), ylim=c(0,500),lty=2, type="l")
   lines(x1,ICymed_sup, col=2)
   lines(x1,ICymed_inf, col=2)
#  2- Intervalo de Predição para E[Y/x_new]
#
   var_pred<-var_log_yme + ((pi^2)/6)
   ICpred_sup<- exp(log_yme + 1.96*var_pred)
   ICpred_inf<- exp(log_yme - 1.96*var_pred)
   plot(x1,exp(log_yme), ylim=c(0,3000),lty=2, type="l")
   lines(x1,ICymed_sup, col=2)
   lines(x1,ICymed_inf, col=2)
   lines(x1,ICpred_sup, col=4)
   lines(x1,ICpred_inf, col=4)
#
#=====================================================================
#
#  Ajuste do Livro do Farrell (2001)
#
#   mesmo modelo do Cox mas com splines no x
#
  require(rms)
  f<-psm(Surv(y,cens)~rcs(x,3),dist='exponential')
  f
#
#  aparentemente o spline não é necessário
#
#=========================================================================