library(splines)
library(survival)

bmt <- read.table("c:/Teaching/bmt.dat", header=F)
d <-as.matrix(bmt)

names(bmt) <-c("g", "T1", "T2", "d1", "d2", "d3", "TA", "A", "TC", "C", "TP", "P", "Z1",
     "Z2", "Z3", "Z4", "Z5", "Z6", "Z7", "Z8", "Z9", "Z10")

# The following is to prepare handing ties between TA and T2, and between TC and T2.
# In a tied case, the trick it to set TA (or TC) a little bit earlier than T2
# To do that, we need to fine out the minimum gap between all the event times

allt <-sort(unique(c(bmt$T2),bmt$TA, bmt$TC))
mintdif <- min(c(allt,max(allt)+1)-c(-0.01,allt))
epsi=mintdif/2

n<-nrow(bmt)
nt <-rep(0,n)

for (i in 1:n)
	{

if (bmt$T2[i]< bmt$TA[i]) 
  {bmt$TA[i] <- bmt$T2[i]
   bmt$A[i] <- 0
  } # We should not use acute GVHD data after T2[i], the relapse time

if (bmt$T2[i]< bmt$TC[i]) 
  {bmt$TC[i] <- bmt$T2[i]
   bmt$C[i] <- 0
  } # We should not use data chronic GVHD data after T2[i], the relapse time


if (bmt$T2[i]==bmt$TA[i] && bmt$A==1) bmt$TA[i]=bmt$TA[i]-epsi
if (bmt$T2[i]==bmt$TC[i] && bmt$C==1) bmt$TC[i]=bmt$TC[i]-epsi



nt[i] <- length(unique(c(bmt$T2[i],bmt$TA[i], bmt$TC[i])))

	}

nnew=sum(nt)

x <- matrix(0,nnew,ncol(bmt)+6)

k <- 0

for (i in 1:n)
	{

	ts <- sort(unique(c(0, bmt$T2[i],bmt$TA[i], bmt$TC[i])))

	for (j in 1:nt[i])
		{
		k <- k+1
		
		start <- ts[j]
		end <- ts[j+1]
		acute <- ifelse( (bmt$TA[i]<=start & bmt$A[i]==1),1,0)
		chronic <- ifelse( (bmt$TC[i]<=start & bmt$C[i]==1),1,0)
		status <- 0  # set all relapse indicator=0 at this moment

	    # write to the new data set

		x[k,1] <- i   # use this as an ID for easy data checking 
		x[k,2] <- start
		x[k,3] <- end
		x[k,4] <- acute
		x[k,5] <- chronic
		x[k,6] <- status  # set all relapse indicator=0 at this moment
		x[k,7:28] <- d[i,1:22] # copy the data to the new data set x
		}
	      x[k,6]=d[i,5] # for the last interval, set the relapse indicator = the original relapse indicator

	}

y <- data.frame(x)


names(y)=c("ID","start", "end","acute","chronic","status","g", "T1", "T2", "d1", "d2", "d3", "TA", "A", "TC", "C", "TP", "P", "Z1",
     "Z2", "Z3", "Z4", "Z5", "Z6", "Z7", "Z8", "Z9", "Z10")




result <- summary(coxph(Surv(start,end,status)~acute+chronic+acute*chronic+factor(g)+Z8, data=y))