# This file is stored as  weibull1.R

# This is code for analyzing time to death
# or censoring for women with breast cancer. 
# Data are posted as bcancer.dat

# The variables are as follows:

#Number  Name      Description      
#----------------------------------
#   1    d    Failure indicator:
#               1=death 0=censored
#   2    y    Observed time
#----------------------------------

# Enter the data into a data frame.

set1 <- read.table("c:/stat565/bcancer1.dat",
      header=F, col.names=c("d", "y")) 
set1$z <- rep(1, length(set1$y))

set1

#  Attach the survival analysis library

library(survival)


weib.out <-survreg( Surv(y,d) ~ z , data=set1, 
                    dist=c("weibull"))

summary(weib.out)

names(weib.out)

# Obtain estimates of the parameters
# in our parameterization of the Weibull
# distribution

parms <- weib.out$coef 
parms
parms[2] <- 1/(weib.out$scale) 
parms[1] <- exp(-parms[1]/weib.out$scale) 
parms

# Obtain a large sample estimate of the 
# covariance matrix for the parameters
# using the delta method. First delete the 
# row and column of the covariance matrix
# corresponding to the zero coefficient for z

V <- weib.out$var[-2, -2]

# Convert from log-scale to match 
# SAS output

G <- matrix(0, nrow=2, ncol=2)
G[1,1] <- 1
G[2,2] <- weib.out$scale

V <- G%*%V%*%t(G)
V

# Convert to our Weibull parameterization

G <- matrix(0, nrow=2, ncol=2)
G[1,1] <-  -parms[1]*parms[2]
G[1,2] <-  -parms[1]*log(parms[1])*parms[2]
G[2,2] <-  -(parms[2])^(2)


covparms <- G%*%V%*%t(G)
stdparms <- sqrt(diag(covparms))
parms
stdparms
covparms

# Obtain precentiles of the estimated
# Weibull distribution and standard
# errors

p <- c(.1, .25, .5, .75, .9)

tp <- ((log(1/(1-p)))/parms[1])^(1/parms[2])

G2 <- matrix(0, nrow=length(tp), ncol=2)
G2[ ,2] <-  -tp*log(tp)/parms[2]
G2[ ,1] <-  -tp/(parms[1]*parms[2])
stdtp <- sqrt(diag(G2%*%covparms%*%t(G2)))
tp
stdtp