# This file is stored as  weibull2.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
#   3    x    Treatment:
#                1=treated 0=control                
#-----------------------------------

# Enter the data into a data frame.

set1 <- read.table("c:/st565/bcancer2.dat",
      header=F, col.names=c("d", "y", "x")) 

set1

library(survival)

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

summary(weib2.out)


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

parms <- c(weib2.out$coef,weib2.out$scale) 
parms
parms[3] <- 1/parms[3] 
parms[2] <- -parms[2]*parms[3]
parms[1] <- exp(-parms[1]*parms[3]) 
parms

# Obtain a large sample estimate of the 
# covariance matrix for the parameters
# using the delta method.

V <- weib2.out$var

# Convert from log-scale to match 
# SAS output

np <- length(parms)
G <- diag(rep(1,np))
G[np,np] <- weib2.out$scale

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

# Convert to our Weibull parameterization

G <- diag(rep(-parms[np],np))
G[1,1] <-  -parms[1]*parms[np]
G[ , np] <-  -parms*parms[np]
G[1,np] <- -parms[1]*parms[3]*log(parms[1])
G[np,np] <-  -(parms[np])^(2)


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

# Obtain precentiles of the estimated
# Weibull distribution and standard
# errors for the treated population

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

tp <- ((log(1/(1-p)))/(parms[1]*exp(parms[2]*x)))^(1/parms[np])

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

# Obtain precentiles of the estimated
# Weibull distribution and standard
# errors for the treated population

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

tp <- ((log(1/(1-p)))/(parms[1]*exp(parms[2]*x)))^(1/parms[np])

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