##############################################################
##  DRAFT OF CODE: Current version available at			##
##  http://faculty.washington.edu/acook/software.html		##
##############################################################

##########################################################################
## CODE FOR THE GS EE COVARIATE ADJUSTED METHOD FOR				## 
## GROUP SEQUENTIALLY MONITORING POSTMARKETING SURVEILLANCE DATA		##
## By: Andrea J Cook									##
##########################################################################

################## Function to run the GS EE Approach ####################
##input:	ATime: a vector of analysis times observed and expected 	## 
##			to be observed							##
##		CurTime: which look current analysis is at 			## 
##		PrevBounds: give vector of previous boundaries if in		## 
##				middle of study						## 
##		Nend: specify total sample size at end of study 		## 
##			(leave blank if at end of study				##
##		propinfo: proportion of statistical information up to each	##
##			observed and expected analysis times (0,1] 		##
##			(example: proportion of total sample size at each look) ##	
##		Y: a vector of observed outcomes					##
##		X: a vector of indicators of being exposed or unexposed 	##
##		Z: a vector of the confounders 					##	
##		S: a vector of the start times			 		##	
##		E: a vector of the exposure durations observed up to		## 
##			CurTime for each individual    				##
##			(only needed for chronic exposure setting)		##	
##        family: specify "binomial" for single exposure binary outcome ##
##			or "poisson" for chronic exposure				##
##		alpha: total type I error desired to spend across all looks	##
##		delta: shape parameter for the boundary 				##
##			(delta=.5 is Pocock) (delta=0 is Fleming)			##
##        nsim: Number of simulations to calculate the boundary		##
##			(use at least 1000, but higher the better			##	
##output: maxstatY: a list that returns alpha, delta, signal (T or F),	##
##		signal time (time of signal or end of study), and 		##
##		test a dataset (ATimet, LLRt, boundt, sigt) 			##
##		at each time point t							## 
##########################################################################
seqEECovAdjUnif<-function(ATime,CurTime=max(ATime),PrevBounds=NULL, 
Nend=length(Y),propinfo,Y,X,Z,S,E=rep(1,length(X)),family="binomial",
	delta=0.5,alpha=0.05,nsim=10000)
{
	Y<-as.matrix(Y[order(S)])
	X<-X[order(S)]
	Z<-as.matrix(Z)
	Z<-Z[order(S),]
	Z<-as.matrix(Z)
	E<-E[order(S)]
	S<-S[order(S)]
	Zint<-cbind(1,Z)

	# Create a vector with amount of sample size observed up to time t
	Ncum<-NULL
	for(T in 1:length(ATime[ATime<=CurTime]))
	{
		Ncum<-c(Ncum,sum(S<=ATime[T]))
	}

	## BINOMIAL MODEL USING LOGISTIC REGRESSION LINK ##
	if(family=="binomial")
	{
# If you are in the middle of the study and need to simulate what 
# will happen at future looks of the dist of Y, X, and E #
		if(length(Y)<Nend)
		{
			numA<-length(propinfo)
			NcumAll<-c(Ncum, round(Nend*propinfo[(length(Ncum)+1):numA],0))

			Nt<-Ncum[1]
			for(T in 2:numA)
			{
				Nt<-c(Nt,NcumAll[T]-NcumAll[T-1])
			}
	
			datcur<-gen.dataFO.bin(Nt=Nt,Yt=Y, Xt=X, Zt=Z, Acur=length(Ncum))
			Y<-datcur$Y
			X<-datcur$X
			Z<-datcur$Z
			Ncum<-NcumAll
		}

		betaZ<-NULL
		for(T in 1:length(ATime))
		{
		betaZ<-cbind(betaZ,glm(Y[1:Ncum[T]]~Z[1:Ncum[T],],family=binomial)$coef)
		}
		#Skip looks that you could not estimate beta
		bfound<-apply(is.na(betaZ),2,sum)
		betaZ<-betaZ[,bfound==0]
		ATime<-ATime[bfound==0]
		propinfo<-propinfo[bfound==0]
		Ncum<-Ncum[bfound==0]

		bd<-seqEECovAdjUnif.boundary.bin(PrevBound=PrevBounds,propinfo=propinfo,Y=Y,X=X,Z=Z,S=S,Ncum=Ncum,delta=delta,betaZ=betaZ,alpha=alpha,nsim=nsim)	 		
		sig<-0
	 	sigtime<-max(ATime)
	 	ScT<-NULL

	 	boundScTt<-NULL
	 	boundESt<-NULL
	 	ATimet<-NULL
	 	sigt<-NULL
	 	T<-1
	 	while(T<=length(ATime) & ATime[T]<=CurTime)
	 	{
			ScTcur<-ScTt.bin(Y=Y[1:Ncum[T]],X=X[1:Ncum[T]], Z=Z[1:Ncum[T],],betaZ=betaZ[,T])
				ScT<-c(ScT,ScTcur)
				boundScTt<-c(boundScTt,bd[T])
				ATimet<-c(ATimet,ATime[T])
				sigt<-c(sigt,as.numeric(ScTcur>boundScTt[T]))
			if(ScTcur>boundScTt[T])
			{
				sig<-1
				sigtime<-ATime[T]
				#T<-length(ATime)+1
			}
			T<-T+1
	 	}
	 	test<-data.frame(cbind(ATimet,ScT,boundScTt,sigt))
	 	return(list(alpha=alpha,delta=delta,signal=(sig==1), sigtime=sigtime,test=test,bd=bd))
	}
	## POISSON MODEL USING LOG LINK ##
	if(family=="poisson")
	{
# If you are in the middle of the study and need to simulate what # will happen at future looks of the dist of Y, X, and E #
		if(length(Y)<Nend)
		{
			numA<-length(propinfo)
			NcumAll<-c(Ncum, round(Nend*propinfo[(length(Ncum)+1):numA],0))

			Nt<-Ncum[1]
			for(T in 2:numA)
			{
				Nt<-c(Nt,NcumAll[T]-NcumAll[T-1])
			}
			datcur<-gen.dataFO.pois(Nt,ATime,Yt=Y,Xt=X,Zt=Z,Et=E,St=S,
 				Acur=length(Ncum))
			Y<-datcur$Y
			X<-datcur$X
			Z<-datcur$Z
			E<-datcur$E
			S<-datcur$S
			Ncum<-NcumAll
		}

		Etime<-NULL
		Ytime<-NULL
		for(T in 1:length(ATime))
		{
			Ecur<-E
			Ecur[(ATime[T]-S+1)<E]<-(ATime[T]-S+1)[(ATime[T]-S+1)<E]
			Ecur[S>ATime[T]]<-0
			Ycur<-Y
			Ycur[(ATime[T]-S+1)<E]<-0
			Etime<-cbind(Etime,Ecur)
			Ytime<-cbind(Ytime,Ycur)
		}
		#Skip looks if there are not at least 5 cases ##
		sumY<-apply(Ytime,2,sum)
		Etime<-Etime[,sumY>=5]
		Ytime<-Ytime[,sumY>=5]
		ATime<-ATime[sumY>=5]
		propinfo<-propinfo[sumY>=5]
		Ncum<-Ncum[sumY>=5]

		betaZ<-NULL

		for(T in 1:length(ATime))
		{
			betaZcur<- glm(Ytime[1:Ncum[T],T]~Z[1:Ncum[T],]+ offset(log(Etime[1:Ncum[T],T])),family="poisson")$coef
			betaZ<-cbind(betaZ,betaZcur)
		}
		#Skip looks that you could not estimate beta
		bfound<-apply(is.na(betaZ),2,sum)
		betaZ<-betaZ[,bfound==0]
		Etime<-Etime[,bfound==0]
		Ytime<-Ytime[,bfound==0]
		ATime<-ATime[bfound==0]
		propinfo<-propinfo[bfound==0]
		Ncum<-Ncum[bfound==0]
	
		bd<-seqEECovAdjUnif.boundary.pois(PrevBound=PrevBounds,propinfo=propinfo,Ytime=Ytime,X=X,Z=Z,S=S,Etime=Etime,
			Ncum=Ncum,delta=delta,betaZ=betaZ,alpha=alpha,nsim=nsim)			
		sig<-0
		sigtime<-max(ATime)
		ScT<-NULL

		boundScTt<-NULL
		boundESt<-NULL
		ATimet<-NULL
		sigt<-NULL
		T<-1
		while(T<=length(ATime) & ATime[T]<=CurTime)
	 	{
			ScTcur<-ScTt.pois(Y=Ytime[1:Ncum[T],T],X=X[1:Ncum[T]],
 				Z=Z[1:Ncum[T],],E=Etime[1:Ncum[T],T], betaZ=betaZ[,T])
			ScT<-c(ScT,ScTcur)
			boundScTt<-c(boundScTt,bd[T])
			ATimet<-c(ATimet,ATime[T])
			sigt<-c(sigt,as.numeric(ScTcur>boundScTt[T]))
			if(ScTcur>boundScTt[T])
			{
				sig<-1
				sigtime<-ATime[T]
				#T<-length(ATime)+1
			}
			T<-T+1
		}
		test<-data.frame(cbind(ATimet,ScT,boundScTt,sigt))
	 	return(list(alpha=alpha,delta=delta,signal=(sig==1), 
sigtime=sigtime,test=test))
	}
}

### REST OF CODE JUST FOR INTERNAL FUNCTION CALLS ###
###############Generate Data for future Observations #####################
##Observations as Unit                                                  ##
##input:   Nt: Vector of number of people you plan to observe/have	##
##			observed at each look						##
##         Yt: Vector of number of drugs you have observed up to time t	##
##         Xt: Vector of indicator of exposed/unexposed you have 		##
##			observed up to time t						##
##         Zt: Matrix of confounders that you have observed up to time t##	
##         Et: Vector of exposure times observed up to time t		##
##         Acur: Current look that we are on at time t			##	
##output:  datcur: a list of outcome data Y, X, Z, E, S			##
##########################################################################
gen.dataFO.bin <- function(Nt, Yt, Xt, Zt, Acur)
{
	Y<-Yt
	X<-Xt
	Z<-as.matrix(Zt)

	## Assumes that current distribution of Y, X, and Z will stay constant
 	## throughout the rest of the study 
	## sample with replacement from current Y, X and Z 
	id<-c(1:length(Xt))
	for(i in (Acur+1):length(Nt))
	{
		lvals<-sample(id,Nt[i],replace=T)
		Y<-c(Y,Yt[lvals])
		X<-c(X,Xt[lvals])
		Z<-rbind(Z,as.matrix(Zt[lvals,]))
	}
	return(list(Y=Y,X=X,Z=Z))
}

gen.dataFO.pois <- function(Nt, ATime, Yt, Xt, Zt, Et, St, Acur)
{
	Y<-Yt
	X<-Xt
	Z<-as.matrix(Zt)
	E<-Et
	S<-St
## Assumes that current distribution of X, Z, and E will stay constant ## throughout the rest of the study
	## sample with replacement from current X, E, and Z 
	id<-c(1:length(Xt))
	for(i in (Acur+1):length(Nt))
	{
		lvals<-sample(id,Nt[i],replace=T)
		Y<-c(Y,Yt[lvals])
		X<-c(X,Xt[lvals])
		Z<-rbind(Z,as.matrix(Zt[lvals,]))
		E<-c(E,Et[lvals])
		S<-c(S,runif(Nt[i],ATime[i-1]+1,ATime[i]))
	}
	return(list(Y=Y,X=X,Z=Z,E=E,S=S))
}

##############Calculate stat#############################################
#input: St: a vector of score test statistics each look       		#
#	  propinfo: vector of total sample size up to look t			#
#       delta: shape parameter (0.5=Pocock, 0=O'Brien Fleming		#
#output: stat: a vector of stat each look (same dimension as score)	#
#########################################################################
stat <- function(St, propinfo, delta)
{
	stat <- (propinfo)^(1-2*delta)*St
	return(stat)
}

#####################   Calculate Score Statistic  #######################
##input:    X: a vector indicating if on exposure of interest		##
##          Y: a vector or matrix of outcomes for look t			##
##          Z: a matrix of confounders for X at look t			##
##          E: a vector of exposure durations at look t			##	 
##output: ScT: the value of the Score Statistic at each look statistic	##
##########################################################################
## Case 1: Single Exposure Time
ScTt.bin <- function(Y, X, Z,betaZ=NULL)
{
	Y<-as.matrix(Y)
	Zint<-cbind(1,Z)

	# Makes sure that there are at least 5 cases
	if((sum(X*Y)==0)|(sum((1-X)*Y)==0)|sum(Y)<5)
	{
		return(0)
	}
	else
	{
		n<-length(Y)
		mui<-as.vector(exp(Zint%*%betaZ)/(1+exp(Zint%*%betaZ)))
		ei<-as.vector(Y-mui)
		XZ<-cbind(X,Zint)
		Ui<-XZ*ei
		U<-t(rep(1,length(Y)))%*%Ui
		W<-ginv(t(XZ)%*%(XZ*mui*(1-mui)))
		V<-W%*%(t(Ui)%*%Ui)%*%W
		ScT<-U%*%W[,1]%*%W[1,]%*%t(U)/V[1,1]
		ScT[((U[1,1]<0)|(is.na(t(X)%*%ei)))]<-0
		return(ScT)
	}
}

## Case 2: Chronic Exposure Time
ScTt.pois <- function(Y, X, Z, E,betaZ=NULL)
{
	Y<-as.matrix(Y)
	Zint<-cbind(1,Z)

	# Makes sure that there are at least 5 cases
	if((sum(X*Y)==0)|(sum((1-X)*Y)==0)|sum(Y)<5)
	{
		return(0)
	}
	else
	{
		n<-length(Y)
		
		mui<-as.vector(exp(Zint%*%betaZ)*E)
		ei<-as.vector(Y-mui)

		XZ<-cbind(X,Zint)

		Ui<-XZ*ei
		U<-t(rep(1,length(Y)))%*%Ui
		W<-ginv(t(XZ)%*%(XZ*mui))
		V<-W%*%(t(Ui)%*%Ui)%*%W
		ScT<-U%*%W[,1]%*%W[1,]%*%t(U)/V[1,1]
		ScT[((U[1,1]<0)|(is.na(ScT)))]<-0
	
		return(ScT)
	}
}


### Find Boundary on the score statistic given unifying boundary function ###
##input:  Ncum: a vector of the cumulative number of people to be		## 
##			observed at each look						##
##	propinfo: vector of total sample size up to look t			##
##		X: a vector of indicator of exposure data observed at each  ##
##			look up to current look 					##	
##		Z: a vector of confounder data observed at each look up to	##
## 			current look							##
##		S: a vector of the start times observed up to current look	##
##		Etime: a matrix of exposure times observed at each look	## 
##				up to current look					##
##      delta: shape of the boundary (delta=0.5 is Pocock; 0 is Fleming)##	
##output: maxstatY: a vector of the critical values at each look t     	## 
##########################################################################
do.Cr.bin <- function(Ncum,PrevBounds,propinfo,Y,X,Z,S,delta,betaZ)
{
	T<-length(Ncum)
	N<-Ncum[T]
	Ncum0<-c(1,Ncum)
	
	ScTY<-NULL
	Xcur<-NULL
	for(i in 1:T)
	{
		Xcur<-c(Xcur,sample(X[Ncum0[i]:Ncum0[i+1]],Ncum0[i+1]-
Ncum0[i]+1,replace=F))
		ScTY<-c(ScTY,ScTt.bin(Y=Y[1:Ncum[i]], X=Xcur[1:Ncum[i]],
Z=Z[1:Ncum[i],],betaZ=betaZ[,i]))
	}
	
	if(!is.null(PrevBounds))
	{
		CurTime<-length(PrevBounds)
		PrevBoundsComp<-c(PrevBounds,rep(10000,T-CurTime))
		sigPL<-as.numeric(sum(ScTY>PrevBoundsComp)>0)
		maxstatY <- max(stat(ScTY[CurTime+1:T], propinfo, delta),na.rm=T)
		return(cbind(sigPL,maxstatY))
	}
	else 
	{
		maxstatY <- max(stat(ScTY, propinfo, delta),na.rm=T)
		return(cbind(sigPL=0,maxstatY))
	}
}

do.Cr.pois <- function(Ncum,PrevBounds,propinfo,Ytime,X,Z,S,Etime, 
delta,betaZ)
{
	T<-length(Ncum)
	N<-Ncum[T]
	Ncum0<-c(1,Ncum)
	
	ScTY<-NULL
	Xcur<-NULL
	for(i in 1:T)
	{
		Xcur<-c(Xcur,sample(X[Ncum0[i]:Ncum0[i+1]],Ncum0[i+1]-
Ncum0[i]+1,replace=F))
		ScTY<-c(ScTY,ScTt.pois(Y=Ytime[1:Ncum[i],i], X=Xcur[1:Ncum[i]],
Z=Z[1:Ncum[i],], E=Etime[1:Ncum[i],i],betaZ=betaZ[,i]))
	}
	if(!is.null(PrevBounds))
	{
		CurTime<-length(PrevBounds)
		PrevBoundsComp<-c(PrevBounds,rep(10000,T-CurTime))
		sigPL<-as.numeric(sum(ScTY>PrevBoundsComp)>0)
		maxstatY <- max(stat(ScTY[CurTime+1:T], propinfo, delta),na.rm=T)
		return(cbind(sigPL,maxstatY))
	}
	else 
	{
		maxstatY <- max(stat(ScTY, propinfo, delta),na.rm=T)
		return(cbind(sigPL=0,maxstatY))
	}
}

################  Function to obtain Critical Boundary Values	########
##input:	PrevBounds: a vector of previous boundaries if not 		##
##			at first look							##
##		propInfo: Proportion of statistical information observed at ##
##			each time look							##
##		X: a vector of the current indicator of exposure data 	##
##			observed up to time t 						##
##		Z: a vector of the confounders observed up to time t 		##
##		S: a vector of the start times observed up to time t 		##
##		E: a vector of the exposure durations observed up to time t ##
##		Bo: Estimated baseline from model without X			##
##		Bz: Estimated confounder relationship without X			##
##         delta: shape parameter of boundary 					##
##			(delta=.5 is Pocock; 0 is Fleming)				##
##         nsim: Number of simulations to calculate the boundary		##
##output: maxstatY: a vector of the critical values at each look t     	## 
##########################################################################
seqEECovAdjUnif.boundary.bin<-function(PrevBounds,propinfo, 
Y,X,Z,S,Ncum,delta,betaZ,alpha=0.05,nsim=1000)
{
	temp<-replicate(nsim,do.Cr.bin(Ncum=Ncum,PrevBounds=PrevBounds, 
propinfo=propinfo,Y=Y,X=X,Z=Z,S=S,delta,betaZ=betaZ))
	# Previous Looks alpha spent #
	alphaP<-mean(temp[1,])
	print(alphaP)
	qtail<-quantile(temp[2,temp[1,]==0],1-(alpha-alphaP),na.rm=T,type=9)
	return(CV=c(PrevBounds, 
(propinfo^(2*delta-1)*qtail)[length(PrevBounds):length(propinfo)]))
}

seqEECovAdjUnif.boundary.pois<-function(PrevBounds,propinfo, 
Ytime,X,Z,S,Etime,Ncum,delta,betaZ,alpha=0.05,nsim=1000)
{ 	
temp<-replicate(nsim,do.Cr.pois(Ncum=Ncum,PrevBounds=PrevBounds,propinfo=propinfo, 
 Ytime=Ytime,X=X,Z=Z,S=S,Etime=Etime,delta,betaZ=betaZ))
	# Previous Looks alpha spent #
	alphaP<-mean(temp[1,])	
	print(alphaP)
	qtail<-quantile(temp[2,temp[1,]==0],1-(alpha-alphaP),na.rm=T,type=9)
	return(CV=c(PrevBounds,(propinfo^(2*delta-1)*qtail)[length(PrevBounds):length(propinfo)]))
	
}