##
######
######		Title: Session1Code
######		Author:	J. Kittelson
######		Date: 08.14.2013 (UW Summer Institue in Biostatistics - Group Sequential Methods)
######		Description:	Code for Session 1 (Introduction)					
######
###
#***************************
#************* Positive Predictive Value
#***************************

pwr <- c(0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.975)
pwr <- cbind(pwr,(1.96 + qnorm(pwr))^2/3.92^2 *1000)
pwr[,2] <- ceiling(pwr[,2])
pwr[,2] <-c(30,50,80,130,200,250,300,400,500,680,1000)

PPVsum <- function(pi0,a1,b1,a2,b2, N=1000000) {
	N1 <- pwr[pwr[,1] == b1,2]
	N2 <- pwr[pwr[,1] == b2,2]
	T1 <- (N/N1)/(1+(N2/N1)*(pi0*b1 + (1-pi0)*a1))
	T2 <- T1*(pi0*b1 + (1-pi0)*a1)
	TP1 <- T1*pi0*b1
	FP1 <- T1*(1-pi0)*a1
	TP2 <- TP1*b2
	FP2 <- FP1*a2
	T2 = TP1 + FP1
	rslt <- c(N1=N1,T1=T1,TP1=TP1,FP1=FP1,PPV1=TP1/(FP1+TP1),T2=T2,TP2=TP2, FP2=FP2, PPV2=TP2/(FP2+TP2))
	rslt
	}

*** scenario 3:
zz <- PPVsum(0.2,0.05,0.15,0.05,0.8)
round(zz,3)
       N1        T1       TP1       FP1      PPV1        T2       TP2       FP2      PPV2 
   50.000 11764.706   352.941   470.588     0.429   823.529   282.353    23.529     0.923 
zz[2]*50 + zz[6]*500
	
*** scenario 4:
zz <- PPVsum(0.01,0.05,0.15,0.05,0.8)
round(zz,3)
       N1        T1       TP1       FP1      PPV1        T2       TP2       FP2      PPV2 
   50.000 13245.033    19.868   655.629     0.029   675.497    15.894    32.781     0.327 
zz[2]*50 + zz[6]*500

*** scenario 5:
zz <- PPVsum(0.1,0.05,0.15,0.05,0.975)
round(zz,3)
      N1       T1      TP1      FP1     PPV1       T2      TP2      FP2     PPV2 
  50.000 9090.909  136.364  409.091    0.250  545.455  132.955   20.455    0.867 
zz[2]*50 + zz[6]*1000

*** scenario 6:
zz <- PPVsum(0.1,0.05,0.15,0.05,0.5)
round(zz,3)
       N1        T1       TP1       FP1      PPV1        T2       TP2       FP2      PPV2 
   50.000 15384.615   230.769   692.308     0.250   923.077   115.385    34.615     0.769 
zz[2]*50 + zz[6]*250


*** scenario 7:
zz <- PPVsum(0.1,0.2,0.15,0.05,0.80)
round(zz,3)
      N1       T1      TP1      FP1     PPV1       T2      TP2      FP2     PPV2 
  50.000 6779.661  101.695 1220.339    0.077 1322.034   81.356   61.017    0.571 
zz[2]*50 + zz[6]*500


*** scenario 8:
zz <- PPVsum(0.1,0.2,0.15,0.1,0.80)
round(zz,3)
     N1       T1      TP1      FP1     PPV1       T2      TP2      FP2     PPV2 
  50.000 6779.661  101.695 1220.339    0.077 1322.034   81.356  122.034    0.400 
zz[2]*50 + zz[6]*500


****** summary table:
PPVsum <- function(pi0,a1,b1,a2,b2, N=1000000) {
	N1 <- pwr[pwr[,1] == b1,2]
	N2 <- pwr[pwr[,1] == b2,2]
	T1 <- (N/N1)/(1+(N2/N1)*(pi0*b1 + (1-pi0)*a1))
	T2 <- T1*(pi0*b1 + (1-pi0)*a1)
	TP1 <- T1*pi0*b1
	FP1 <- T1*(1-pi0)*a1
	TP2 <- TP1*b2
	FP2 <- FP1*a2
	T2 = TP1 + FP1
	rslt <- c(pi0=pi0,a1=a1,b1=b1,a2=a2,b2=b2,TP2=TP2,FP2=FP2, PPV2=TP2/(FP2+TP2))
	rslt
	}


#***************************
#*****Goup sequential sample path illustration
#***************************
#***************************

Plot.Bound <- function(dsgn,col=1,newplot=T,ylim=NULL,...) {
	bnd <- as.matrix(dsgn$boundary)
	if(is.null(ylim)) ylim <- range(bnd[!is.inf(bnd)]) + c(-1,1)*2
	N <- dsgn$param$sample.size
	StopLo <- !is.inf(bnd[,1])
	StopHi <- !is.inf(bnd[,4])
	StopMid <- bnd[,2] < bnd[,3]
	if (newplot) plot(c(0,max(N)),ylim,type="n",ylab="Mean Effect",xlab="Sample Size",ylim=ylim,...,cex=2)
	bnd <- cbind(ylim[1],bnd,ylim[2])
	if(sum(StopLo) > 1) {
		lines(N,bnd[,2],lty=2,col=col)
		points(N,bnd[,2],pch=4,col=col)
		lines(rep(N[StopLo],each=3),t(cbind(bnd[StopLo,1:2],NA)),lty=4,col=col,lwd=2)
	}
	if(sum(StopHi) > 1) {
		lines(N,bnd[,5],lty=2,col=col)
		points(N,bnd[,5],pch=4,col=col)
		lines(rep(N[StopHi],each=3),t(cbind(bnd[StopHi,5:6],NA)),lty=1,col=col,lwd=2)
	}
	if(sum(StopMid) > 1) {
		lines(N[StopMid],bnd[StopMid,3],lty=2,col=col)
		points(N[StopMid],bnd[StopMid,3],pch=4,col=col)
		lines(N[StopMid],bnd[StopMid,4],lty=2,col=col)
		points(N[StopMid],bnd[StopMid,4],pch=4,col=col)
		lines(rep(N[StopMid],each=3),t(cbind(bnd[StopMid,3:4],NA)),lty=3,col=col,lwd=2)
	}
	J <- length(N)
	lines(c(N[J],N[J]),bnd[J,1:2],lty=4,col=col,lwd=2)
	lines(c(N[J],N[J]),bnd[J,3:4],lty=3,col=col,lwd=2)
	lines(c(N[J],N[J]),bnd[J,5:6],lty=1,col=col,lwd=2)
}

### simulated sample paths:
sup.D <- seqDesign(prob.model = "normal", arms = 1, null.hypothesis = 0., 
	alt.hypothesis = 3.92, variance = 1., nbr.analyses = 5, 
	sample.size = 1, test.type = "greater", power = "calculate",
	alpha = 0.025, epsilon = c(0., 1.), early.stopping = "alternative", 
	display.scale = seqScale(scaleType = "X"))
sup.A <- update(sup.D,early.stopping="null")
sup.DA <- update(sup.D,early.stopping="both")
obf.fixN <- sup.DA


dsgn <- obf.fixN
dsgn2 <- update(obf.fixN,epsilon=c(1,1),P=c(1,Inf,Inf,1))
Plot.Bound(dsgn,ylim=c(-12,12))

nrep <- 100
bnd <- matrix(as.numeric(dsgn$boundary),nc=4)
for(i in 1:nrep) {
N <- (1:100)/100
x <- rnorm(100,mean=3.92,sd=10)
x <- cumsum(x)/(1:100)
j <- (1:5)*20
StopHi <- x[j] > bnd[,4]
StopLo <- x[j] < bnd[,1]
Stop <- StopHi | StopLo
Stop <- max((1:5)[cumsum(Stop) == 0]) + 1
if(is.na(Stop)) Stop <- 1
lines(N[1:(Stop*20)],x[1:(Stop*20)],lwd=2)

}