Capítulo 2 
 ============================================================================
 Script em R para obtenção de Var(tm) com correção do vício (pg 50, 51 e 54) 
 Autora: Vanessa F. Sehaber
 ============================================================================

  v.tmedio <- function(ekm){					
   su <- summary(ekm)
   temp <- su$time
   nj <- su$n.risk      
   dj <- su$n.event
   st <- su$surv         
   k <- length(temp)-1                  
   A <- sapply(1:k, function(i){
        st[i]*(temp[i+1]-temp[i])
        })
   Aj <- sapply(1:k, function(i){
        sum(A[i:k])
        })
   den <- nj*(nj-dj); den <- den[1:k]
   num <- dj[1:k]*(Aj^2)
   m<-sum(dj)
   var <- num/den; var <- sum(var)*(m)/(m-1)
   return(var)
  }
 var.tmedio <- v.tmedio(ekm); var.tmedio
 =============================================================================

 #################### Exemplo ###################################
 require(survival)
 tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 ekm<-survfit(Surv(tempos,cens)~1) 
 source("http://people.ufpr.br/~giolo/Livro/var_tm.txt")
 var.tmedio <- v.tmedio(ekm); var.tmedio
 ################################################################

			 Capítulo 3 
 =============================================================================
 Gama generalizada: ajuste usando pacote "flexsurv" do software R
 Obs: necessário instalar o pacote flexsurv
 =============================================================================
 > tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 > cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 > require(flexsurv)
 > ajust4<-flexsurvreg(Surv(tempos,cens)~1,dist='gengamma')
 > ajust4
 > ajust4$loglik
   [1] -65.69074
 =============================================================================

			 Capítulo 3         
 =============================================================================
 Além do método delta para obtenção de var(tm) no caso do modelo Lognormal ou
 outros modelos paramétricos, esta também pode ser obtida por meio do procedi-
 mento conhecido por reamostragem bootstrap como segue.  
 =============================================================================

 [A] Obtenção de var(tm) usando reamostragem bootstrap - Modelo lognormal  
 =============================================================================
 tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 data<-as.data.frame(cbind(tempos,cens))
 
 d<-5000                             # numero de reamostragens bootstrap
 n<-length(tempos)                   # n = tamanho amostral 
  tmb<-matrix(0,d,1)                 # tempos médios das reamostragens 
  for(j in 1:d){
    row<-sample(n, replace=T)        # reamostragem com reposição      
    dados<-data[row,]
    mod<-survreg(Surv(tempos,cens)~1,dist='lognorm', data=dados)
    tmb[j]<-exp(mod$coef[1]+(mod$scale^2/2))
   }
  var(tmb)			     # var(tm) bootstrap	
  hist(tmb, breaks=20)

 ### I.C.(tm)_95% utilizando a variância bootstrap ### 

  mod1<-survreg(Surv(tempos,cens)~1,dist='lognorm', data=data)
  mu<-mod1$coef[1]; sigma<-mod1$scale
  tm<-as.numeric(exp(mu+(sigma^2/2)))
  li<-tm-1.96*sqrt(var(tmb)); ls<-tm+1.96*sqrt(var(tmb))
  cbind(tm, li, ls) 

 =============================================================================
 [B] Obtenção de var(tm) usando reamostragem bootstrap - Modelo Weibull 
 =============================================================================

 tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 data<-as.data.frame(cbind(tempos,cens))
 
 d<-5000                             # numero de reamostragens bootstrap
 n<-length(tempos)                   # n = tamanho amostral 
  tmb<-matrix(0,d,1)                 # tempos médios das d reamostragens 
  for(j in 1:d){
    row<-sample(n, replace=T)        # reamostragem com reposição      
    dados<-data[row,]
    mod<-survreg(Surv(tempos,cens)~1,dist='weibull', data=dados)
    a<-exp(mod$coef[1]) 
    g<-1/mod$scale 
    tmb[j]<-a*gamma(1+1/g)
   }
  var(tmb)                           # variância bootstrap
  hist(tmb, breaks=20)

  ### I.C.(tm)_95% utilizando a variância bootstrap ### 

  mod1<-survreg(Surv(tempos,cens)~1,dist='weibull', data=data)
  a<-exp(mod$coef[1]); g<-1/mod$scale 
  tm<-as.numeric(a*gamma(1+1/g))
  li<-tm-1.96*sqrt(var(tmb)); ls<-tm+1.96*sqrt(var(tmb))
  cbind(tm, li, ls) 
 =============================================================================

   				Capítulo 3 
 =============================================================================
 Analogamente, a variância do tempo mediano também pode ser obtida por meio de 
 reamostragem bootstrap como segue.  
 =============================================================================  

 [A] Obtenção de var(tmed) usando reamostragem bootstrap - Modelo lognormal  
 =============================================================================  
 tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 data<-as.data.frame(cbind(tempos,cens)) 
 d<-5000                        # numero de reamostragens bootstrap
 n<-length(tempos)              # n = tamanho amostral 
  tmd<-matrix(0,d,1)            # tempos medianos das reamostragens 
  for(j in 1:d){
    row<-sample(n, replace=T)   # reamostragem com reposição      
    dados<-data[row,]
    mod<-survreg(Surv(tempos,cens)~1,dist='lognorm', data=dados)
    tmd[j]<-exp(mod$coef[1])
   }
  v<-var(tmd); v		# variância bootstrap
  hist(tmd, breaks=20)          

  ### I.C.(tmed)_95% utilizando a variância bootstrap ### 

  mod1<-survreg(Surv(tempos,cens)~1,dist='lognorm', data=data)
  tmed<-as.numeric(exp(mod1$coef[1]))
  li<-tmed-1.96*sqrt(v); ls<-tmed+1.96*sqrt(v)
  cbind(tmed, li, ls) 
 =============================================================================