###############################################

## Joshua M. Tebbs

## Date: 22 May 2010

## Update: 9 August 2011

## STAT 513 course notes: R Code

## Chapter 13

###############################################

######################################################

## IMPORTANT ##

## The user (YOU) must install the splinesurv package in R.

## The user (YOU) must then load the package.

######################################################

# Example 13.1.

# CSP+MTX vs MTX only clinical trial data

# Page 133-135

# Input CSP+MTX data

csp.mtx<-c(3,8,10,12,16,17,22,64,65,77,82,98,155,189,199,247,

324,356,378,408,411,420,449,490,528,547,691,769,1111,1173,1213,1357)

# Create censoring indicator for CSP+MTX group (0=censored)

cens.csp.mtx<-c(0,1,1,0,1,1,1,0,1,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1)

# Input MTX only data

mtx<-c(9,11,12,20,20,22,25,25,25,28,28,31,35,35,46,49,

104,106,156,218,230,231,316,393,395,428,469,602,681,690,1112,1180)

# Create censoring indicator for MTX only group (0=censored)

cens.mtx<-c(0,1,1,0,1,1,1,0,1,1,1,1,0,0,1,0,0,0,1,1,0,0,0,1,0,0,1,1,0,1,0,1)

# Combine time and censoring indicator data into one vector

agvhd.times<-c(csp.mtx,mtx)

cens.agvhd.times<-c(cens.csp.mtx,cens.mtx)

# Treatment indicator (1=CSP+MTX; 2=MTX only)

# "rep(a,b)" creates a vector of length b, consisting entirely of a's

treat<-c(rep(1,32),rep(2,32))

# Fit KM estimates for each treatment

fit<-survfit(Surv(agvhd.times,cens.agvhd.times)~treat)

# Plot KM estimates

plot(fit,xlab="Time to AGVHD (in days)",ylab="Survival Probability",lty=c(1,4))

abline(h=0)

# Create legend; click where you want legend to go

names<-c("CSP+MTX","MTX only")

legend(locator(1),legend=names,lty=c(1,4),bty="n")

# Example 13.2.

# Exponential survivor function

# Page 136-137

time<-seq(0,5,0.001)

surv<-exp(-time)

plot(time,surv,type='l',ylab="Surival function, S(t)",xlab="Time (t)")

abline(h=0)

# Example 13.3

# Plots of Weibull hazard functions

# Page 139-140

t<-seq(0,4,0.001)

haz.1<-3*t^2

haz.2<-1.5*t^(0.5)

haz.3<-1+t*0

haz.4<-0.5*t^(-0.5)

par(mfrow=c(2,2))

plot(t,haz.1,type='l',xlab="Time (t)",ylab="Hazard function")

plot(t,haz.2,type='l',xlab="Time (t)",ylab="Hazard function")

plot(t,haz.3,type='l',xlab="Time (t)",ylab="Hazard function")

plot(t,haz.4,type='l',xlab="Time (t)",ylab="Hazard function")

# Example 13.5.

# KM estimate with fictitious data

# Page 146-148

# Input the failure/censoring times

survtime<-c(4.5,7.5,8.5,11.5,13.5,15.5,16.7,17.5,19.5,21.5)

# Input the failure/censoring indicators

# 1 = failure; 0 = censored

status<-c(1,1,0,1,0,1,1,0,1,0)

# Use internal 'survfit' function

# Use conf.type="plain" for confidence bands

# Use '~1' because there is only one survival curve

fit<-survfit(Surv(survtime,status)~1,conf.type="plain")

summary(fit)

plot(fit,xlab="Patient time (in years)",ylab="Survival probability")

# Example 13.6.

# Simulation exercise

# Page 150-151

# Generate the data

# Failure mean = 5; Censoring mean = 10; n=100 patients

survtime<-rexp(100,1/5) # T_i's

censtime<-rexp(100,1/10) # C_i's

# Creates indicator variable (Delta),

# 1 (TRUE), if surv < cens; 0 (FALSE), otherwise

delta<-(survtime <= censtime)

# obstime are the observe data we would see (the X_i's)

obstime<-survtime*delta+censtime*(1-delta)

# Compute KM estimate

# Use '~1' because there is only one survival curve

fit<-survfit(Surv(obstime,delta)~1,conf.type="plain")

summary(fit)

plot(fit,xlab="Patient time (in years)",ylab="Survival probability")

# Place horizontal line at S(t)=0.5

abline(h=0.5,lty=1)

# Display first 5 simulated patients

# These will be different than those in the notes!

results<-data.frame(survtime[1:5],censtime[1:5],obstime[1:5])

results<-round(results,3)

results

# Example 13.7.

# Tongue cancer data

# Page 152-153

# Load data in R

tongue<-read.table("C:\\texfiles\\Classes\\USC\\stat513\\f10\\data\\tongue.txt",header=TRUE)

cancer.type<-tongue[,1]

time<-tongue[,2]

delta<-tongue[,3]

# Compute KM estimate

# Use '~1' because there is only one survival curve

fit<-survfit(Surv(time,delta)~1,conf.type="plain")

summary(fit)

plot(fit,xlab="Patient time (in weeks)",ylab="Survival probability")

abline(h=0.5,lty=1)

# Gives 95 percent CI for median survival

fit

# Example 13.8.

# HAART AIDS data

# Page 159-162

# Input treatment 1 data

# Put censored observations at the end

survtime.1<-c(14,17,128,129,164,228,333,444,558,568,677,702,706,909,

1213, 1420,1527,1730,1834,2246,2565,3004,1216,2244)

# Input treatment 2 data

# Put censored observations at the end

survtime.2<-c(64,178,478,533,742,756,863,998,1205,1232,1232,1433,

1873,1993,1999,2140,2361,2380,2680,2696,2896,3223,2204,3344)

# Combines data into one vector

survtime<-c(survtime.1,survtime.2)

# Create death/censoring indicator

# rep(a,b) creates a vector of length b, consisting entirely of a's

delta<-c(rep(1,22),0,0,rep(1,22),0,0)

# Treatment indicator: 1's and 2's

treat<-c(rep(1,24),rep(2,24))

# fit KM estimates for each treatment

fit.1<-survfit(Surv(survtime,delta) ~ treat)

# Logrank test for difference in survivor distributions

fit.2<-survdiff(Surv(survtime,delta) ~ treat)

# Plot KM estimates

plot(fit.1,xlab="Patient time (in days)", ylab="Survival probability",

lty=c(1,4))

# Create legend; click where you want legend to go

names<-c("Trt 1 (monotherapy)","Trt 2 (HAART)")

legend(locator(1),legend=names,lty=c(1,4),bty="n")

# Returns CIs and logrank test results

fit.1

fit.2

# Example 13.10.

# Design study example

# Page 168-169

# Constant accrual rate (no. patients per unit of time)

a=100

# Hazard rates (constant, since exponential assumption)

lambda.1=0.173

lambda.2=0.116

# Number of events (deaths)

d=256

# root.func is a function of x

# A=L (assumed to be equal), both equal to "x"

root.func<-function(x)

{(a/2)*(x-(exp(-lambda.1*x)/lambda.1)*(exp(lambda.1*x)-1)) +

(a/2)*(x-(exp(-lambda.2*x)/lambda.2)*(exp(lambda.2*x)-1)) - d}

## find the value of x which solves root.func

## option $root returns the root (solution)

uniroot(root.func, c(0,50), tol=0.00001)$root

# Example 13.11.

# Breast cancer example

# Page 172-173

# Input treatment 1 data

survtime.1<-c(501,1721,4280,3350,3142,4167,3266,894,4454,2360,4610,

665,3660,2067,3260,653,3684,4197,2320,3905)

# Input treatment 2 data

survtime.2<-c(1959,354,1157,95,2729,2385,625,1716,2574,1169,357,1666,

1464,3146,76,1199,3750,1206,391,1847)

# Input treatment 3 data

survtime.3<-c(4067,3494,1323,1992,1482,1305,3230,71,4117,4002,974,4205,

2734,3634,3302,3436,4372,1853,989,1712)

# Combines data into one vector

survtime<-c(survtime.1,survtime.2,survtime.3)

# Create death/censoring indicator

delta<-c(rep(1,16),0,0,0,0,rep(1,18),0,0,rep(1,17),0,0,0)

# Treatment indicator: 1's, 2's, and 3's

treat<-c(rep(1,20),rep(2,20),rep(3,20))

# Fit KM estimates for each treatment

fit.1<-survfit(Surv(survtime,delta) ~ treat, conf.type="plain")

# Logrank test for difference in survivor distributions

fit.2<-survdiff(Surv(survtime,delta) ~ treat)

# Plot KM estimates

# Create legend; click where you want legend to go

plot(fit.1,xlab="Patient time (in days)", ylab="Survival probability",lty=c(1,2,4))

names<-c("Intensive CAF","Low CAF","Standard CAF")

legend(locator(1),legend=names,lty=c(1,2,4),bty="n")

# Returns CIs and logrank test results

fit.1

fit.2