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

##########################################################################
## DRAFT OF CODE FOR THE CONDITIONAL SEQUENTIAL SAMPLING PROCEDURE 	## 
## (CSSP) FOR GROUP SEQUENTIAL MONITORING OF POSTMARKETING SURVEILLANCE ##
## DATA. 											##
## Programming: Robert D. Wellman							##
## CSSP method developed by LingLing Li (see references)			##
##												##
## Current version available at							##
## http://faculty.washington.edu/acook/software.html				##
##########################################################################

##########################################################################
## DESCRIPTION
## cssp is used to analyze postmarketing surveillance data sequentially
## using stratification for confounding control.
##########################################################################

##########################################################################
## USAGE
## 
## cssp(data,path=NULL,this.look,tstat.sim=10000,alpha=0.05,total.size,
##      current.size,boundary.type=1,rho=1) 
##########################################################################


##########################################################################
## ARGUMENTS
##
## data
## a data frame or list (or object coercible by as.data.frame to a data
## frame) containing variables needed to use method (see details below  
## for data specification).
##  
## path
## File path indicator folder where interim files are stored. If NULL 
## then uses working directory (getwd()). 
##
## this.look
## a integer giving the current look in terms of a sequence of looks
## starting from 1. Defaults to 1.
##
## tstat.sim
## an integer giving the number of simulations used to obtain the 
## distribution of the test statistic under the null in each strata.
##
## alpha
## overall alpha level for the study, e.g., 0.05 is commonly used and
## equates to a type I error of 5%.
##
## total.size
## total sample size to be accumulated by at the end of the study
## (last analysis).
##
## current.size
## sample size at the time of this.look.
##
## boundary.type
## an integer giving the desired boundary shape: 1=Pocock, 
## 2=O'Brien-Fleming, 3=Power function. If the power function is used,
## an additional parameter, rho, must be specified. Defaults to 1.
##
## rho
## a numeric constant > 0 indicating the desired exponent for 
## power function type boundary.
##########################################################################

##########################################################################
## DETAILS 
## 
## Data must be formatted as given in the included CSSPexample.Rdata file.
## The data must already be summarized by stratum and include the
## following variables named exactly as specified here:
##
## look	integer giving look number in a sequence beginning at 1 with
##		no maximum
## 
## strata	integer indexing strata valued from 1 to number of strata
##
## subj_cmp	integer giving the number of subjects in the comparison
##		group within each look and each stratum
##
## subj_doi	integer giving the number of subjects in the group on the
##		drug of interest within each look and each stratum
##
## exp_cmp	numeric giving the aggregated exposure time for the comparison
##		group within each look and each stratum
##
## exp_doi	numeric giving the aggregated exposure time for the group on
##		the drug of interest (DOI) with each look and each stratum
##
## evt_cmp	integer giving the aggregated number of events in the comparison
##		group within each look and each stratum
##
## evt_doi	integer giving the aggregated number of events in the group on
##		the drug of interest within each look and each stratum
##########################################################################

##########################################################################
## VALUE
##
## cssp returns an object of class "CSSP".
##
## A summary of the analysis by look is printed. The summary function will
## extract the same summary from the CSSP object. After each analysis the
## resulting CSSP object is stored in a user-specified directly. If no
## directory is passed into the cssp function then the current working
## directory is used. The file is saved as "t_sim.Rdata". After the first
## look, cssp will load the previously saved file from the directory
## specified by the user. The augmented object will then be saved over
## the previous version but will contain all information through the 
## current look. Directories and file names can be manipulated in
## the code if they need to be modified for a particular situation.
##########################################################################

##########################################################################
## EXAMPLES
## source("H:/CSSPfunctions.R")
## options(digits=10)
## ex1<-read.csv("H:/example1.csv", header=T)
## look1<-cssp(ex1,this.look=1,tstat.sim=10000,alpha=0.05,total.size=10000,current.size=1250,boundary.type=1,rho=1)
## look2<-cssp(ex1,this.look=2,tstat.sim=10000,alpha=0.05,total.size=10000,current.size=2500,boundary.type=2,rho=1)
## summary(look2)
##
##        Events P-value Look Alpha Cum. Alpha Signal Num. Sims
## Look 1      1  0.5539     0.0063     0.0063      0     10000
## Look 2     10  0.2331     0.0063     0.0126      0     10000
##########################################################################


##########################################################################
## REFERENCES
##
## Li, L. A conditional sequential sampling procedure for drug safety
## surveillance. Statistics in Medicine 2009; 28:3124-3138. 
##########################################################################


##########################################################################
## FUNCTIONS (Load before using cssp)

binom.gen <- function(data, H, R) rbinom(1,data[H],data[R])

comp.r<-function(x,c1,c2) x[,c1]/(x[,c1]+x[,c2])

sim.T<-function(data, H,R,S) sum(rbinom(max(data[,S]),data[,H],data[,R]))
	
cssp.current<-function(data,nsim,H,R,S){  
  res.sim<-replicate(nsim,sim.T(data, H,R,S))
  return(res.sim)
}

#Error Spending Function (esf)
esf<-function(type=1,t,T,alpha,rho=NULL){
  if(type==1) {alpha*log(1 + (exp(1)-1)*(t/T))}
    else if(type==2) {2*(1-pnorm(abs(qnorm(alpha/2))/(sqrt(t/T))))}
    else if(type==3) {alpha*(t/T)^(rho)}
  }

## Class constructor function

CSSP<-function(input.list){
  CSSP<-input.list
  class(CSSP)<-"CSSP"
  return(CSSP)
}

summary.CSSP<-function(x){
  print(x[['Summary']])
}

cssp<-function(data,path=NULL,this.look,tstat.sim=10000,alpha=0.05,total.size,current.size,boundary.type=1,rho=1){

  if(is.null(path)){path<-getwd()}
  dset<-as.data.frame(data)
  l<-this.look
  if(l==1){
    dset<-subset(dset,look==1)
    alpha.1<-esf(boundary.type,current.size,total.size,alpha,rho)
    cut.rank<-ceiling(tstat.sim*(1-alpha.1))
    dset$h<-apply(dset[,c("evt_doi","evt_cmp")], 1, sum)
    dset$r<-comp.r(dset,"exp_doi","exp_cmp")
    res<-cssp.current(subset(dset,look==l),tstat.sim,H="h",R="r",S="strata")
    t.sim.raw<-list('Look 1'=res)
    save(t.sim.raw,file=paste(path,"/tsimraw.Rdata",sep=""))
    res.sort<-rank(res,ties.method='first')
    a.spent.t<-sum(res.sort>=cut.rank)/tstat.sim
    obs.t<-sum(dset$evt_doi)
    order.obs.t<-sum(res>=obs.t)
    p.val.t<-order.obs.t/tstat.sim
    signal<-p.val.t<a.spent.t
    res1<-ifelse(res.sort>cut.rank, 9999, res)
    output<-list('Summary'=matrix(c(obs.t,p.val.t,a.spent.t,a.spent.t,signal,tstat.sim),nrow=1,
      dimnames=list(c('Look 1'),c('Events','P-value','Look Alpha','Cum. Alpha','Signal',
      'Num. Sims'))),'Look1'=res1)
    save(output,file=paste(path,"/t_sim.Rdata",sep=""))
    #print(output[["Summary"]])
    class(output)<-"CSSP"
    output
    }
  else if(l>1){
    load(file=paste(path,"/t_sim.Rdata",sep=""))
    load(file=paste(path,"/tsimraw.Rdata",sep=""))
    if(l > dim(output$Summary)[1] + 1){stop('Looks out of sequence.')}
    dset.this<-subset(data,look==l)
    dset.all<-data
    alpha.l<-esf(boundary.type,current.size,total.size,alpha,rho)
    cut.rank<-ceiling(tstat.sim*(1-alpha.l))
    dset.this$h<-apply(dset.this[,c("evt_doi","evt_cmp")], 1, sum)
    dset.this$r<-comp.r(dset.this,"exp_doi","exp_cmp")
    res<-cssp.current(subset(dset.this,look==l),tstat.sim,H="h",R="r",S="strata")
    t.sim.raw[[paste('Look',l)]]<-res
    save(t.sim.raw,file=paste(path,"/tsimraw.Rdata",sep="")) 
    res.current<-output[[l]] + res
    res.sort<-rank(res.current, ties.method='first')
    a.spent.t<-sum(res.sort>=cut.rank)/tstat.sim
    prev.a.spent.t<-a.spent.t-output[['Summary']][l-1,'Cum. Alpha']
    obs.t<-sum(subset(dset.all,look<=l)$evt_doi)
    order.obs.t<-sum(res.current>=obs.t)
    p.val.t<-order.obs.t/tstat.sim
    signal<-p.val.t<a.spent.t
    res.current<-ifelse(res.sort>cut.rank, 9999, res.current)
    output[['Summary']]<-rbind(output[['Summary']],
      c('Events'=obs.t,'P-value'=p.val.t,'Look Alpha'=prev.a.spent.t,
      'Cum. Alpha'=a.spent.t,'Signal'=signal,'Num. Sims'=tstat.sim))
    output[[paste("Look",l,sep='')]]<-res.current
    rownames(output[['Summary']])[l]<-c(paste('Look',l))
    save(output,file=paste(path,"/t_sim.Rdata",sep=""))
    #print(output[["Summary"]])
    class(output)<-"CSSP"
    output
    }
}

##########################################################################