#USO:
# plot.basehas(formul='Surv(dat$time,dat$cens)~strata(dat$sexo)+dat$idade',dat)
"plot.basehas" <-
  function(formul,datafrm,logt=TRUE,lwd=1){
    formul<-formula(formul)
    ss <- survfit(coxph(formul,data=datafrm))
    ns <- length(unique(ss$strata))
    n0 <- 0
    a <-  log(- log(ss$surv))
    if(logt) {
      plot(c(0, 0), c(0, 0), pch = " ", xlab = "", ylab = "", xlim = 
         range(log(ss$time)), ylim = range(a),cex=1.5,lwd=lwd)
      mtext('log(Tempo)',1,l=2.3,cex=1.5)
      mtext(expression(paste("log(",Lambda,")")),2,l=2.3,cex=1.5)
      for(it in 1:ns) {
        nu <- ss$strata[it]
        j <- n0 + (1:nu)
        lines(log(ss$time[j]), a[j],lwd=lwd)
        n0 <- n0 + nu
      }
    }
    else{
      plot(c(0, 0), c(0, 0), pch = " ", xlab = "", ylab = "", xlim = 
         range(ss$time), ylim = range(a),cex=1.5,lwd=lwd)
      mtext('t',1,l=2.3,cex=2.0)      
      mtext("log A(t)",2,l=2.3,cex=2.0)
      for(it in 1:ns) {
        nu <- ss$strata[it]
        j <- n0 + (1:nu)
        lines(ss$time[j], a[j],lwd=lwd)
        n0 <- n0 + nu
      }
    }
    cat("\n")
  }

# USO:
# plot.pi('coxph object')
"plot.pi"  <-   function(fit,lwd=1,...)
{
  if(is.null(fit$x)) 
    cat("Rode a função coxph() com o argumento x=TRUE\n")
  else{
    x<-fit$x
    ss <- survfit(fit)
    cat("The solid line is the fitted model, the dashed one is the K-M\n")
    if(ncol(x)>1) xbar <- apply(x,2,"mean")
    else xbar <- mean(as.vector(x))
    
    b<-fit$coefficients
    bb<-b%*%xbar
    bx<-x%*%b
    j<-order(bx)  
    n<-fit$n
    bx1<-bx[j][round(n/3)]
    bx2<-bx[j][round(2*n/3)]
    
    i1<-bx<=bx1
    i2<-(bx>bx1)&(bx<=bx2)
    i3<-(bx>bx2)
    b1<-mean(bx[i1])
    b2<-mean(bx[i2])
    b3<-mean(bx[i3])
    s<-ss$surv
    s1<-s^{exp(b1-bb)}
    s2<-s^{exp(b2-bb)}
    s3<-s^{exp(b3-bb)}
    plot(ss$time,s1,type="s",ylim=c(0,1),xlim=c(0,max(ss$time)),lwd=lwd,...)
#    mtext("Time",1,l=2.3,cex=1.7)
#    mtext("Survival",2,l=2.3,cex=1.7)
    lines(ss$time,s2,type="s",lwd=lwd)
    lines(ss$time,s3,type="s",lwd=lwd)
    if(length(terms(fit$formula)[[2]])==4){
      start<-fit$y[,1]
      stop<-fit$y[,2]
      cens<-fit$y[,3]
      ss1<-survfit(Surv(start[i1],stop[i1],cens[i1]))
      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
      ss1<-survfit(Surv(start[i2],stop[i2],cens[i2]))
      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
      ss1<-survfit(Surv(start[i3],stop[i3],cens[i3]))
      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
    }
    else{
      time<-fit$y[,1]
      cens<-fit$y[,2]
      ss1<-survfit(Surv(time[i1],cens[i1]))
      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
      ss1<-survfit(Surv(time[i2],cens[i2]))
      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
      ss1<-survfit(Surv(time[i3],cens[i3]))
      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
    }
    cat("\n")
  }
}
  
  
  

#### USO:
#### plot.pi('coxph object')
###"plot.pi"  <-   function(fit,lwd=1,...)
###{
###  if(is.null(fit$x)) 
###    cat("Rode a função coxph() com o argumento x=TRUE\n")
###  else{
###    x<-fit$x
###    #fra <- rep(fit$frail,table(x[,ncol(x)]))
###    #x[,ncol(x)] <- fra
###        
###    ss <- survfit(fit)
###    
###    cat("The solid line is the fitted model, the dashed one is the K-M\n")
###    if(ncol(x)>1) xbar <- apply(x,2,"mean")
###    else xbar <- mean(as.vector(x))
###    
###    #b<-c(fit$coefficients,fra)
###    b<-fit$coefficients
###    
###    bb<-b%*%t(xbar)
###    bx<-x%*%t(bb)
###    j<-order(bx)  
###    n<-fit$n
###    bx1<-bx[j][round(n/3)]
###    bx2<-bx[j][round(2*n/3)]
###    
###    i1<-bx<=bx1
###    i2<-(bx>bx1)&(bx<=bx2)
###    i3<-(bx>bx2)
###    b1<-mean(bx[i1])
###    b2<-mean(bx[i2])
###    b3<-mean(bx[i3])
###    s<-ss$surv
###    s1<-s^{exp(b1-bb)}
###    s2<-s^{exp(b2-bb)}
###    s3<-s^{exp(b3-bb)}
###    plot(ss$time,s1,type="s",ylim=c(0,1),xlim=c(0,max(ss$time)),lwd=lwd,...)
####    mtext("Time",1,l=2.3,cex=1.7)
####    mtext("Survival",2,l=2.3,cex=1.7)
###    lines(ss$time,s2,type="s",lwd=lwd)
###    lines(ss$time,s3,type="s",lwd=lwd)
###    if(length(terms(fit$formula)[[2]])==4){
###      start<-fit$y[,1]
###      stop<-fit$y[,2]
###      cens<-fit$y[,3]
###      ss1<-survfit(Surv(start[i1],stop[i1],cens[i1]))
###      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
###      ss1<-survfit(Surv(start[i2],stop[i2],cens[i2]))
###      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
###      ss1<-survfit(Surv(start[i3],stop[i3],cens[i3]))
###      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
###    }
###    else{
###      time<-fit$y[,1]
###      cens<-fit$y[,2]
###      ss1<-survfit(Surv(time[i1],cens[i1]))
###      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
###      ss1<-survfit(Surv(time[i2],cens[i2]))
###      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
###      ss1<-survfit(Surv(time[i3],cens[i3]))
###      lines(ss1$time,ss1$surv,lty=2,type="s",lwd=lwd)
###    }
###    cat("\n")
###  }
###}


#change point fit
"coxph.changep"  <-   function(formul,datafrm,cptau)
{ 
  formul <- formula(formul)
  time <- datafrm[,as.character(attr(terms(formul),"variables")[[2]][[2]])]
  cens <-  datafrm[,as.character(attr(terms(formul),"variables")[[2]][[3]])]
  x <- attr(terms(formul),"term.labels")
  l <- length(x)
  v <- x[1]
  if(l>1){
    for (i in 2:l){v <- paste(v,"+",x[i])}
  }
  form <- formula(paste("Surv(time,cens)~",v))
  usedat <- data.frame(time=time,cens=cens,datafrm[,attr(terms(formul),"term.labels")])
  tt<-time
  ic<-cens
  tt1<-tt
  ic1<-ic
  i<-tt1>cptau
  tt1[i]<-cptau
  ic1[i]<-0
  tt2<-tt
  ic2<-ic
  i<-tt2<=cptau
  ic2[i]<-0
  cc<-coxph(form,data=usedat)
  l0<-cc$loglik[2]
  cat("No change\n")
  print(cc)
  cat("\n-----------------------------------------------------\n")
  usedat$time <- tt1
  usedat$cens <- ic1
  cc1<-coxph(form,data=usedat)
  l1<-cc1$loglik[2]
  cat("Before change\n")
  print(cc1)
  cat("\n-----------------------------------------------------\n")
  usedat$time <- tt2
  usedat$cens <- ic2
  cc2<-coxph(form,data=usedat)
  l2<-cc2$loglik[2]
  cat("After change\n")
  print(cc2)
  cat("\n-----------------------------------------------------\n")
  cat("Log-likelihoods\n\nNo change\n")
  print(l0)
  cat("Before change\n")
  print(l1)
  cat("After change\n")
  print(l2)
  cat("Total with change\n")
  print(l1+l2)
  
list(all=cc,before=cc1,after=cc2)  
}

#Gráfico das fragilidades e seus intervalos de confiança (Valeska)
#alterado por Marilia: 
    #fragilidades em escala log (a escala exp dificulta a visualizacao por causa da assimetria)
    #corrigidos o label q estava sem ordenaçao
    # generalizado para os dois tipo de obejto: quando tem sparse=T (default), que aparece o vetor frail
    # e sparse=F, onde a funcao cria o vetor frail
    
plot.frail<-function(unidade,model,...){
	ifelse(is.null(model$frail),{
	    vetorfra <- model$assign[length(model$assign)][[1]]
	    model$frail <- as.numeric(model$coeff[vetorfra])
	    model$fvar <- diag(model$var)[vetorfra]
	    },NA)

	if(length(model$fvar)<=1){
	       warning("Este modelo não estimou fragilidade")
	      }
	 else{
	    fragil <- data.frame(unidade=unique(unidade),fragil=model$frail,inf=(model$frail - 1.96*sqrt(model$fvar)),sup=(model$frail + 1.96*sqrt(model$fvar)))
	    ordenado <- fragil[order(fragil[,2]),]                     
	    x<-matrix(1:nrow(fragil),ncol=nrow(fragil),nrow=2,byrow=T) 
	    y<-t(ordenado[,3:4])                                       
	    matplot(x,y,type='l',col=1,lty=1,axes=F,...)               
	    box();axis(2); axis(1,at=x[1,],labels=ordenado$unidade,las=2)
	    matpoints(x[1,],y[1,],pch=24,col=1)                        
	    matpoints(x[2,],y[2,],pch=25,col=1)                        
	    points(1:ncol(x),ordenado[,2],pch=19)                      
	    abline(h=0)
	    }
}