##
######
######		Title:			Session5Code
######		Author:			D. Gillen
######		Date:			08.14.2013 (UW Summer Institue in Biostatistics - Group Sequential Methods)
######		Description:	Code for Session 5 (Monitoring and Constrained Boundaries)					
######
##
library( RCTdesign )
##
#####
#####	Impact of changing the number and timing of analyses
#####
##
dsn <- seqDesign( 	prob.model="normal", arms=1, null.hypothesis=0, 
		alt.hypothesis=0.5, test.type="greater", variance=4,
		power=0.975, P=0.5, nbr.analyses=4, early.stopping="both" ) 
dsn

##
#####	Timing of analyses
##
dsn.late.power <- update(dsn, sample.size=c(.4,.6,.8,1) )
dsn.late.power

dsn.late.n <- update(dsn, sample.size=c(.4,.6,.8,1)*max(dsn$parameters$sample.size), alt.hypothesis="calculate" )
dsn.late.n

##
#####	Number of analyses
##
dsn.5.power <- update(dsn, sample.size=c(.2,.4,.6,.8,1) )
dsn.5.power

dsn.5.n <- update(dsn, sample.size=c(.2,.4,.6,.8,1)*max(dsn$parameters$sample.size), alt.hypothesis="calculate" )
dsn.5.n


##
#####
#####	Modify the OBF design at the first analysis using exact constraint	
#####
##
obf <- seqDesign( 	prob.model="normal", arms=1, null.hypothesis=0, 
					alt.hypothesis=0.5, test.type="greater", variance=4,
					power=0.975, P=1, nbr.analyses=4, early.stopping="both" ) 
obf
seqBoundary(obf, scale="P")

bnd.const <- as.seqBoundary( cbind(matrix(NA,nrow=4,ncol=3),c(.0005,rep(NA,3))), scale="P" )
obf.const <- update( obf, exact.constraint=bnd.const )
obf.const
seqBoundary(obf.const, scale="P")

##
#####	Compare operating characteristics
##
plot( obf, obf.const, dsnLbls=c("OBF", "Constrained OBF") )
seqPlotPower( obf, obf.const, lwd=2, dsnLbls=c("OBF", "Constrained OBF"), fixed=FALSE  )
seqPlotPower( obf.const, lwd=2, reference=obf, fixed=FALSE  )
seqPlotASN( obf, obf.const, lwd=2, dsnLbls=c("OBF", "Constrained OBF") )


##
#####
#####	Example of Error Spending Functions
#####
##
#####	Pre-planned timing of analyses (information time)	
##
sepsis.fix <- seqDesign(prob.model="proportions", arms=2, size=.025, power="calculate",
null.hypothesis= c(.30, .30), alt.hypothesis=c(0.25,0.30), sample.size=1700,
test.type="less") 

#****** pre-trial monitoring plan
sepsis.obf <- update(sepsis.fix,nbr.analyses=4,P=1)
sepsis.obf

V <- .25*(.75) + .3*.7
I <-  (1700/2)/V
cbind(V, I)
(c(425,850,1275,1700)/2)/.3975

##
#####	Hypothetical data at first analysis
##
Y.1 <- c(rep(1,52),rep(0,263-52),rep(1,65),rep(0,257-65))
tx.1 <- c(rep(1,263),rep(0,257))
y1 <- sum(Y.1[tx.1==1])
n1 <- length(Y.1[tx.1==1])
theta1 <- y1/n1 

y0 <- sum(Y.1[tx.1==0])
n0 <- length(Y.1[tx.1==0])
theta0 <- y0/n0
cbind(n1,y1,theta1,n0,y0,theta0)

s.1 <- ( theta1*(1-theta1)/n1 ) + ( theta0*(1-theta0)/n0 )
I.1 <- 1/s.1
Pi.1 <- I.1 / I
cbind(s.1, I.1, Pi.1)


## Pre-planned error spending realizations
seqOC(sepsis.obf, theta=0)
cbind( 	t(seqOC(sepsis.obf, theta=0)$lower.stop.prob),
cumsum( seqOC(sepsis.obf, theta=0)$lower.stop.prob ),
cumsum( seqOC(sepsis.obf, theta=0)$lower.stop.prob ) / .025 )
## Compare with scale="E"		
seqBoundary( sepsis.obf, scale="E" )

##
####	Use linear interpolation to estimate cumulative error to be spent
##
alpha <- cumsum( seqOC(sepsis.obf, theta=0)$lower.stop.prob )
cum.alpha.1 <- alpha[1] + (alpha[2]-alpha[1])*((Pi.1-.25)/(.5-.25))  
cum.alpha.1	## Compare with 0.00091872 in notes --> rounding error

##	New stopping rule on Z-stat scale
qnorm(cum.alpha.1)
qnorm(0.00091872)

## Stopping boundary on sample mean scale (w/ and w/o roundoff error)
qnorm(cum.alpha.1) * sqrt(s.1)
qnorm(0.00091872) * sqrt(0.0013384)

##	Compare with observed difference in sample proportions
theta1-theta0


##
#####
#####	Constrained boundaries example
#####
##
##	Find design parameter (G)
sepsis.obf$parameters
G <- 2.0032
Pi <- 1:4/4
aj <- -G*(1/Pi)*sqrt(V/850)
dj <- (-2*G + G*(1/Pi))*sqrt(V/850)
cbind(Pi,Pi*1700,aj,dj)


sepsis.IA1 <- update(sepsis.obf,sample.size=c(520,850,1275,1700),null.hypothesis=c(65/257,65/257), alt.hypothesis=c(52/263,65/257))
sepsis.IA1
sepsis.IA1$parameters
G.1 <- 2.0036
V.1 <- theta1*(1-theta1) + theta0*(1-theta0)
new.I <- I <-  (1700/2)/V.1
new.Pi <- c(520,850,1275,1700)/1700
new.aj <- -G.1*(1/new.Pi)*sqrt(V.1/850)
new.dj <- (-2*G.1 + G.1*(1/new.Pi))*sqrt(V.1/850)
cbind(new.Pi,new.Pi*1700,new.aj,new.dj)


##
#####	Alternatively, use seqMonitor()
##
IA1 <- seqMonitor(sepsis.obf,response=Y.1,treatment=tx.1,future.analyses=c(850,1275,1700))
IA1





##
#####
#####  Monitoring the Hodgkins trial with constrained boundaries
#####
##
##
#####	Specification of designs
##
Fixed <- seqDesign( prob.model = "hazard", arms = 2, null.hypothesis = 1., 
							alt.hypothesis = 0.67, ratio = c(1., 1.), nbr.analyses = 1, 
							test.type = "less", power = 0.80, alpha = 0.025 )
Eff11.Fut8 <- update( Fixed, nbr.analyses = 4, P=c(1.1,.8), sample.size=196, power="calculate" )

##
#####	Table HR, Zstat, and Fixed P boundaries for Eff11.Fut8 design
##
bndStat <- round( cbind(	seqBoundary( Eff11.Fut8, scale="X" )[,1],
							seqBoundary( Eff11.Fut8, scale="Z" )[,1],
							1-seqBoundary( Eff11.Fut8, scale="P" )[,1],
							seqBoundary( Eff11.Fut8, scale="X" )[,4],
							seqBoundary( Eff11.Fut8, scale="Z" )[,4],
							1-seqBoundary( Eff11.Fut8, scale="P" )[,4] ), 3 )
colnames( bndStat ) <- cbind( "lo.hr", "lo.ztat", "lo.pval", "up.hr", "up.zstat", "up.pval" )
bndStat

##
#####	Simulation of (full) data for monitoring example
set.seed( 123456 )
n <- 120
grp1 <- rexp( n, rate=log(2)/1.0 )
grp2 <- rexp( n, rate=(log(2)/1.0)*.70 ) 
trueSurv <- c( grp1, grp2 )
entry <- runif( 2*n, 0, 3 )
grp <- rep( 0:1, each=n ) 

##
#####	Original analyses scheduled at 49, 98, 147, and 196 events
##
seqPHSubjects( Eff11.Fut8, controlMedian=0.75, accrualTime=3, followupTime=1 )

##
#####
#####	First DSMB meeting takes place at 1.5 years after study start 
#####
##
analysisTime <- 1.5
obsSurv <- ifelse( trueSurv + entry <= analysisTime, trueSurv, analysisTime-entry )
event <- ifelse( obsSurv == trueSurv, 1, 0 )
hodgData <- as.data.frame( cbind( grp, entry, obsSurv, event ) )
hodgData <- hodgData[ hodgData$obsSurv > 0, ]
dim( hodgData )
sum( hodgData$event )

##
#####	First analysis takes place after 35 events have occured.
#####		Constrain boundaries and schedule following analyses at 98, 147, 196 events.
##
resp <- Surv( hodgData$obsSurv, hodgData$event )
interim1 <- seqMonitor( Eff11.Fut8, response=resp, treatment=hodgData$grp, future.analyses=c(98,147,196) )
interim1
plot( interim1, lwd=2 )

##
#####	Estimate overall hazard with POOLED data for planning future analyses/accrual rates
##
expFit <- survreg( Surv(obsSurv, event)~1, dist="exponential", data=hodgData )
estHaz <- exp(-expFit$coef)
seqPHSubjects( Eff11.Fut8, controlMedian=log(2)/estHaz, accrualTime=3, followupTime=1 )
seqPHSubjects( Eff11.Fut8, controlMedian=log(2)/estHaz, rate=80, accrualTime=3, followupTime=NULL )


##
#####
#####	Second DSMB meeting takes place at 2.75 years after study start 
#####
##
analysisTime <- 2.75
obsSurv <- ifelse( trueSurv + entry <= analysisTime, trueSurv, analysisTime-entry )
event <- ifelse( obsSurv == trueSurv, 1, 0 )
hodgData <- as.data.frame( cbind( grp, entry, obsSurv, event ) )
hodgData <- hodgData[ hodgData$obsSurv > 0, ]

##
#####	Constrain boundaries and schedule following analyses at 147, 196 events.
##
resp <- Surv( hodgData$obsSurv, hodgData$event )
interim2 <- seqMonitor( interim1, response=resp, treatment=hodgData$grp, future.analyses=c(147,196) )
interim2
plot( interim2, lwd=2 )


##
#####	Estimate overall hazard with POOLED data for planning future analyses/accrual rates
##
expFit <- survreg( Surv(obsSurv, event)~1, dist="exponential", data=hodgData )
estHaz <- exp(-expFit$coef)
seqPHSubjects( interim1, controlMedian=log(2)/estHaz, accrualTime=3, followupTime=1 )
seqPHSubjects( interim1, controlMedian=log(2)/estHaz, rate=80, accrualTime=3, followupTime=NULL )

##
#####
#####	Third DSMB meeting takes place at 3.5 years after study start 
#####
##
analysisTime <- 3.5
obsSurv <- ifelse( trueSurv + entry <= analysisTime, trueSurv, analysisTime-entry )
event <- ifelse( obsSurv == trueSurv, 1, 0 )
hodgData <- as.data.frame( cbind( grp, entry, obsSurv, event ) )
hodgData <- hodgData[ hodgData$obsSurv > 0, ]

##
#####	Constrain boundaries and schedule final analysis at 196 events.
##
resp <- Surv( hodgData$obsSurv, hodgData$event )
interim3 <- seqMonitor( interim2, response=resp, treatment=hodgData$grp, future.analyses=c(196) )
interim3
plot( interim3, lwd=2 )


##
#####	Estimate overall hazard with POOLED data for planning future analyses/accrual rates
##
expFit <- survreg( Surv(obsSurv, event)~1, dist="exponential", data=hodgData )
estHaz <- exp(-expFit$coef)
seqPHSubjects( interim2, controlMedian=log(2)/estHaz, accrualTime=3, followupTime=1 )
seqPHSubjects( interim2, controlMedian=log(2)/estHaz, rate=80, accrualTime=3, followupTime=NULL )


##
#####
#####	Final DSMB meeting takes place at 5 years after study start 
#####
##
analysisTime <- 5
obsSurv <- ifelse( trueSurv + entry <= analysisTime, trueSurv, analysisTime-entry )
event <- ifelse( obsSurv == trueSurv, 1, 0 )
hodgData <- as.data.frame( cbind( grp, entry, obsSurv, event ) )
hodgData <- hodgData[ hodgData$obsSurv > 0, ]
resp <- Surv( hodgData$obsSurv, hodgData$event )
interim4 <- seqMonitor( interim3, response=resp, treatment=hodgData$grp )
interim4
plot( interim4, lwd=2 )