##install.packages(c("lattice", "vcd","survival","colorspace"));

##library them in order to use directly;

  library(lattice)

  library(grid, splines,vcd)

  library(survival)

  library(colorspace)

##Set the seed to make the result produciple;  

set.seed(126)

##Assign the time point;

times <- abs(c(rnorm(200, 5, 10), rnorm(200, 9, 10), rnorm(200, 13, 10), rnorm(200, 20, 10),

rnorm(200, 7, 10), rnorm(200, 11, 10), rnorm(200, 15, 10), rnorm(200, 24, 10),

rnorm(200, 13, 10), rnorm(200, 18, 10), rnorm(200, 23, 10), rnorm(200, 28, 10),

rnorm(200, 16, 10), rnorm(200, 21, 10), rnorm(200, 26, 10), rnorm(200, 31, 10)))

times[times < 0] <- 3

treat <- rep(rep(c('High','Mid','Low','Placebo'),each=200),4)

term <- rep(c('Headache','Dizziness','Nausea','Abdominal Pain'), each=800)

bin1h <- c(rep(0, 40),rep(1,60))

bin1m <- c(rep(0, 45),rep(1,55))

bin1l <- c(rep(0, 50),rep(1,50))

bin1p <- c(rep(0, 65),rep(1,35))

bin2h <- c(rep(0, 65),rep(1,35))

bin2m <- c(rep(0, 70),rep(1,30))

bin2l <- c(rep(0, 80),rep(1,20))

bin2p <- c(rep(0, 90),rep(1,10))

cens <- c(sample(bin1h, 200, replace=TRUE), sample(bin1m, 200, replace=TRUE),

sample(bin1l, 200, replace=TRUE), sample(bin1p, 200, replace=TRUE),

sample(bin2h, 200, replace=TRUE), sample(bin2m, 200, replace=TRUE),

sample(bin2l, 200, replace=TRUE), sample(bin2p, 200, replace=TRUE),

sample(bin1h, 200, replace=TRUE), sample(bin1m, 200, replace=TRUE),

sample(bin1l, 200, replace=TRUE), sample(bin1p, 200, replace=TRUE),

sample(bin2h, 200, replace=TRUE), sample(bin2m, 200, replace=TRUE),

sample(bin2l, 200, replace=TRUE), sample(bin2p, 200, replace=TRUE))

dat <- data.frame(times, treat, term, cens)

aes <- c('Headache','Dizziness','Nausea','Abdominal Pain')

survlist <- NULL

timelist <- NULL

trtlist <- NULL

for(j in 1:4){

              aedat <- dat[dat$term%in%aes[1],]

              nd <- length(unique(aedat$id))

              aest <- NULL

              # Run a Cox model

              fit <- survfit(Surv(times,cens)~treat, data=aedat)

              #Store output for plot

              trtlistj <- list(rep(names(fit$strata),fit$strata))

              survlistj <- list(fit$surv)

              timelistj <- list(fit$time)

             }

# Create one big data set for plotting

nae <- c(800, 800, 800, 800)

aelabel <- c('Headache','Dizziness','Nausea','Abdominal Pain')

aelist <- NULL

for(k in 1:4){

aelistk <- rep(aelabel[k],times=nae[k])

}

panarg <- unlist(aelist)

trtarg <- unlist(trtlist)

survarg <- unlist(survlist)

timearg <- unlist(timelist)

library(vcd)

hclcolors <- rainbow_hcl(4, l=50, start=50, c=75)[c(3,2,4,1)]

colorpalette <- c(hclcolors[1:3], 'grey45')

new.back <- trellis.par.get("background")

new.back$col <- "white"

newcol <- trellis.par.get("superpose.symbol")

newcol$col <- colorpalette[c(3,1,2,4)]

newcol$pch <- c(16,1,4,8)

new.line <- trellis.par.get('superpose.line')

new.line$col <- colorpalette[c(3,1,2,4)]

new.pan <- trellis.par.get("strip.background")                                          

new.pan$col <- c('gray90','white')

trellis.par.set("background", new.back)

trellis.par.set("superpose.symbol", newcol)

trellis.par.set("superpose.line", new.line)

trellis.par.set("strip.background",new.pan)

lines_style=c(1,2,3,6)

lines_widths=c(6,4,3,1)

xyplot(

       (1-survarg)~timearg|panarg, groups=trtarg, type='s',

       scales=list( x=list(alternating=c(1, 1)), y=list(alternating=c(1, 1))

                  ) ,

                   lwd=lines_widths,

                   lty = lines_style,

                   xlab='Time (days)',

                   ylab='Proportion of Subjects Reporting Event',

                   col=trellis.par.get("superpose.symbol")$col[1:4],

                   key=list(lines=list(lty=lines_style, lwd=lines_widths,

                   #the numbers for N are just dummy numbers. In the actual program they would need to be calculated

                   col=trellis.par.get("superpose.symbol")$col[1:4]),

                   text=list(

                             c('High (N=100)','Mid (N=110)','Low (N=94)','Placebo (N=130)'),

                   col=trellis.par.get("superpose.symbol")$col[1:4]

                             ),

                   columns=4,

                   title='Treatment',

                   between=.6)

      )