####################################################################
# R codes to compute bounds and CIs for average causal effects     #
# within principal strata in three-arm trials with noncompliance   #
# Jing Cheng                                                       #
# Last revised June 20, 2006                                       #
####################################################################


# Input data

#Alcohol data

#Z<-c(rep(0,47),rep(1,50),rep(2,44))   # Z is randomization, (0,47): 47 subjects assigned to control (Z=0); (1,50): 50 to treatment A (Z=A); (2,44): 44 to treatment B (Z=B).
#D<-c(rep(0,47),rep(0,15),rep(1,35),rep(0,21),rep(2,23))   # D is treatment-received, (0,47): 47 subjects assigned to control didn't take treatment (Z=0,D=0); (0,15): 15 assigned to A didn't take A (Z=A,D=0); (1,35): 35 assigned to A took A (Z=A,D=A); (0,21): 21 assigned to B didn't take B (Z=B,D=0); (2,23): 23 assigned to B took B (Z=B,D=B).
#Y<-c(rep(0,19),rep(1,28),rep(0,8),rep(1,7),rep(0,14),rep(1,21),rep(0,6),rep(1,15),rep(0,3),rep(1,20))   # Y is outcome, (0,19): 19 subjects who assigned to control and didn't take treatment had Y equal 0 (Z=0,D=0,Y=0); (1,28): 28 subjects who assigned to control and didn't take treatment had Y equal 1 (Z=0,D=0,Y=1); (0,8): 8 subjects who assigned to A and didn't take A had Y equal 0 (Z=A,D=0,Y=0); ...

#Hypothetical data

Z<-c(rep(0,400),rep(1,400),rep(2,400))
D<-c(rep(0,400),rep(0,20),rep(1,380),rep(0,80),rep(2,320))
Y<-c(rep(0,220),rep(1,180),rep(0,16),rep(1,4),rep(0,19),rep(1,361),rep(0,60),rep(1,20),rep(0,96),rep(1,224))

truedata<-cbind(Y,D,Z)

conf<-0.95  # conf is the confidence level wanted for confidence intervals 

r<-1000     # r is the number of bootstrap resampling

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

#Functions

#Function to compute bounds

ACE<-function(data,d){

data<-data[d,]

precision<-0.01

Y<-data[,1]
D<-data[,2]
Z<-data[,3]

Paa<-length(D[((Z==1)*(D==1))==1])/length(D[Z==1])
Poa<-length(D[((Z==1)*(D==0))==1])/length(D[Z==1])
Pob<-length(D[((Z==2)*(D==0))==1])/length(D[Z==2])
Pbb<-length(D[((Z==2)*(D==2))==1])/length(D[Z==2])

Yaa<-Y[((D==1)*(Z==1))==1]
Yao<-Y[((D==0)*(Z==1))==1]
Ybb<-Y[((D==2)*(Z==2))==1]
Ybo<-Y[((D==0)*(Z==2))==1]
Yoo<-Y[((D==0)*(Z==0))==1]

Yaa1<-sum(Yaa)/length(Yaa)
Yao1<-sum(Yao)/length(Yao)
Ybb1<-sum(Ybb)/length(Ybb)
Ybo1<-sum(Ybo)/length(Ybo)
Yoo1<-sum(Yoo)/length(Yoo)

p1<-max(0,Paa-Pob)
p2<-min(Paa,Pbb)
Poabvector<-seq(p1,p2,by=precision)
k<-length(Poabvector)

AOABmina<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]

   YAOABmin <-max(0, 1-(1-Yaa1)/(Poab/Paa))
  
   YOOABmax <- min(1, max(0, 1+(Yoo1-Pob*Ybo1-Poa*Yao1+Pob-Paa)/Poab), max(0,(Yoo1-Poa*Yao1)/Poab), max(0,(Yoo1-Pob*Ybo1)/Poab))
  
   AOABmina[i]<- YAOABmin - YOOABmax 
}
ind<-order(AOABmina,decreasing=F)[1]
AOABmin<-AOABmina[ind]


AOABmaxa<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
 
   YAOABmax <-min(1, Yaa1/(Poab/Paa))

   YOOABmin <- max(0, min(1, (Yoo1-Pob*Ybo1-Poa*Yao1)/Poab), min(1, 1+(Yoo1-Poa*Yao1-Paa)/Poab), min(1, 1+(Yoo1-Pob*Ybo1-Pbb)/Poab))
   
   AOABmaxa[i]<-  YAOABmax - YOOABmin
}
ind<-order(AOABmaxa,decreasing=T)[1]
AOABmax<-AOABmaxa[ind]


BOABmina<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
   
   YBOABmin <-max(0, 1-(1-Ybb1)/(Poab/Pbb))
  
   YOOABmax <- min(1, max(0, 1+(Yoo1-Pob*Ybo1-Poa*Yao1+Pob-Paa)/Poab), max(0, (Yoo1-Poa*Yao1)/Poab), max(0,(Yoo1-Pob*Ybo1)/Poab))

   BOABmina[i]<- YBOABmin - YOOABmax 
}
ind<-order(BOABmina,decreasing=F)[1]
BOABmin<-BOABmina[ind]


BOABmaxa<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
   
   YBOABmax <-min(1, Ybb1/(Poab/Pbb))

   YOOABmin <- max(0, min(1, (Yoo1-Pob*Ybo1-Poa*Yao1)/Poab), min(1,1+(Yoo1-Poa*Yao1-Paa)/Poab), min(1,1+(Yoo1-Pob*Ybo1-Pbb)/Poab))
  
   BOABmaxa[i]<- YBOABmax - YOOABmin
}
ind<-order(BOABmaxa,decreasing=T)[1]
BOABmax<-BOABmaxa[ind]


AOA0mina<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
   
   YAOAOmin <-max(0, 1-(1-Yaa1)/(1-Poab/Paa))
  
   YOOABmin <- max(0, min(1,(Yoo1-Pob*Ybo1-Poa*Yao1)/Poab), min(1,1+(Yoo1-Poa*Yao1-Paa)/Poab), min(1,1+(Yoo1-Pob*Ybo1-Pbb)/Poab))
  
   YOOAOmax <- min(1, max(0, (Yoo1-Poa*Yao1-Poab*YOOABmin)/(Paa-Poab)))

   AOA0mina[i]<- YAOAOmin - YOOAOmax
}
ind<-order(AOA0mina,decreasing=F)[1]
AOA0min<-AOA0mina[ind]


AOA0maxa<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
  
   YAOAOmax <-min(1, Yaa1/(1-Poab/Paa))

   YOOABmax <- min(1, max(0, 1+(Yoo1-Pob*Ybo1-Poa*Yao1+Pob-Paa)/Poab), max(0,(Yoo1-Poa*Yao1)/Poab), max(0,(Yoo1-Pob*Ybo1)/Poab))

   YOOAOmin <- max(0, min(1, (Yoo1-Poa*Yao1-Poab*YOOABmax)/(Paa-Poab)))
  
   AOA0maxa[i]<- YAOAOmax - YOOAOmin
}
ind<-order(AOA0maxa,decreasing=T)[1]
AOA0max<-AOA0maxa[ind]


BO0Bmina<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
   
   YBOOBmin <-max(0, 1-(1-Ybb1)/(1-Poab/Pbb))

   YOOABmin <- max(0, min(1,(Yoo1-Pob*Ybo1-Poa*Yao1)/Poab), min(1,1+(Yoo1-Poa*Yao1-Paa)/Poab), min(1,1+(Yoo1-Pob*Ybo1-Pbb)/Poab))
  
   YOOOBmax <- min(1, max(0, (Yoo1-Pob*Ybo1-Poab*YOOABmin)/(Pbb-Poab)))

   BO0Bmina[i]<- YBOOBmin - YOOOBmax
}
ind<-order(BO0Bmina,decreasing=F)[1]
BO0Bmin<-BO0Bmina[ind]


BO0Bmaxa<-numeric(k)
for (i in 1:k){
   Poab<-Poabvector[i]
  
   YBOOBmax <-min(1, Ybb1/(1-Poab/Pbb))

   YOOABmax <- min(1, max(0, 1+(Yoo1-Pob*Ybo1-Poa*Yao1+Pob-Paa)/Poab), max(0,(Yoo1-Poa*Yao1)/Poab), max(0,(Yoo1-Pob*Ybo1)/Poab))

   YOOOBmin <- max(0, min(1, (Yoo1-Pob*Ybo1-Poab*YOOABmax)/(Pbb-Poab)))
   
   BO0Bmaxa[i]<- YBOOBmax - YOOOBmin
}
ind<-order(BO0Bmaxa,decreasing=T)[1]
BO0Bmax<-BO0Bmaxa[ind]


#Assumptions 1-5

mAOABmin<-max(0, 1-(1-Yaa1)/(Pbb/Paa)) - min(1, max(0, (1/Pbb)*(Yoo1- Pob*Ybo1)))
mAOABmax<-min(1,Yaa1/(Pbb/Paa)) - min(1, max(0, (1/Pbb)*(Yoo1- Pob*Ybo1)))

mAOAOmin<-max(0, 1-(1-Yaa1)/(1-Pbb/Paa)) - min(1, max(0, (1/(Paa-Pbb))*(Pob*Ybo1-Poa*Yao1)))
mAOAOmax<-min(1,Yaa1/(1-Pbb/Paa)) - min(1, max(0, (1/(Paa-Pbb))*(Pob*Ybo1-Poa*Yao1)))

mBOAB<-Ybb1-min(1, max(0, (1/Pbb)*(Yoo1-Pob*Ybo1)))

c(AOABmin, AOABmax, AOA0min, AOA0max, BOABmin, BOABmax, BO0Bmin, 
  BO0Bmax, mAOABmin, mAOABmax, mAOAOmin, mAOAOmax, mBOAB, mBOAB) 
 }


# Function to transform ACE results to bounds format

ToBounds<-function(ACEres){
   CIm<-array(dim=c(7,2))
   for (i in 1:7){
	CIm[i,1]<-ACEres[2*i-1]
	CIm[i,2]<-ACEres[2*i]
	}
    CIm
}


#Function to transform bootstrap CI to bounds

CI.bounds<-function(bootobj,conf,tobound){
	
     CIm<-array(dim=c(7,2))     
     n<-bootobj$R
     bootobj0<-bootobj        
     baddata<-F   
     for( i in c(1:7)){     
       
	  rind<-NULL
	  bootobj<-bootobj0

	  temp<-is.nan(bootobj$t[,2*i-1]*bootobj$t[,2*i])
	  clen<-length(temp[temp==F])
	  bootobj$t<-bootobj$t[temp==F,]	
	
	if(clen<=(n/2)){
	  baddata<-T
	  break	  
	}

	 bootobj$R<-clen
          
         if((tobound[i,1]==(-1))){
	            CIm[i,1]<-(-1)
	     }
                  else{
	          temp1<-boot.ci(bootobj,conf,type="perc",index=(2*i-1))
	          CIm[i,1]<-temp1$perc[4]
		    }
	    
          
	 if(tobound[i,2]==1){
		    CIm[i,2]<-1
                 }
                  else{
	           temp2<-boot.ci(bootobj,conf,type="perc",index=(2*i))
	           CIm[i,2]<-temp2$perc[5]

		  }
	}
     
	list(CIm=CIm,baddata<-baddata)
}

# Function to compute B-method confidence interval

Bci<-function(Ls,Us,Ln,Un,alfa){
	temp<-Ls*Us
	Ls<-Ls[is.nan(temp)==F]
	Us<-Us[is.nan(temp)==F]
	n<-length(Ls)
	Ld<-Ls-Ln
	Ud<-Un-Us
	Cd<-numeric(n)
	for (i in 1:n){
		Lcdf<-length(Ld[Ld<=Ld[i]])/n
		Ucdf<-length(Ud[Ud<=Ud[i]])/n
		Cd[i]<-max(Lcdf,Ucdf)	}
	c<-quantile(Cd,1-alfa)
	L0<-Ln-quantile(Ld,c,na.rm=T)
	U0<-Un+quantile(Ud,c,na.rm=T)
	if (L0<(-1)){L0<-(-1)}
	if (U0>1){U0<-1}	
	c(L0,U0)}


# Function to cmpute z* in H-M CI

searchZ<-function(Ls,Us,alfa,Ln,Un){
        temp<-Ls*Us
	Ls<-Ls[is.nan(temp)==F]
	Us<-Us[is.nan(temp)==F]
	z<-1/2
	z0<-0
	z1<-2
	n<-length(Ls)
	flag<-T
	if (Ln==-1||Un==1){
		z<-0
		flag<-F	}
	while(flag){
		bc<-((Ls-z)<=Ln)*((Us+z)>=Un)		
		p<-length(bc[bc==T])/n
		d<-p-1+alfa
		inte<-z1-z0
		if(abs(d)<=0.0001||inte<=0.0001){
			flag<-F		}
		
		if (d>0){			    
				  z1<-z
				  z<-(z0+z)/2				
				}else{
				  z0<-z
				  z<-(z1+z)/2			
	
                        }			
        }
	z	
}


#Bounds with the real data

n<-length(Z)
Bounds<-ToBounds(ACE(truedata,1:n))
rownames(Bounds)<-c("A_OAB","A_OAO", "B_OAB", "B_OOB", "A_OAB_monot", "A_OAO_monot", "B_OAB_monot")
colnames(Bounds)<-c("Lower Bound", "Upper Bound")
print("Bounds")
print(Bounds)



#bootstrap

btemp<-boot(truedata,ACE,R=r)
 

#bootstrap Bonferroni CI

BonfCI<-CI.bounds(btemp,conf,Bounds)[[1]]
rownames(BonfCI)<-c("A_OAB","A_OAO", "B_OAB", "B_OOB", "A_OAB_monot", "A_OAO_monot", "B_OAB_monot")
colnames(BonfCI)<-c("Lower Limit", "Upper Limit")
print("Bonferroni CI")
print(BonfCI)


# bootstrap H-M CI

HMCI<-array(dim=c(6,2))
for (s in 1:6){
Ls<-btemp$t[,2*s-1]
Us<-btemp$t[,2*s]
Ln<-btemp$t0[2*s-1]
Un<-btemp$t0[2*s]
z<-searchZ(Ls,Us,1-conf,Ln,Un)
HMCI[s,]<-c(max(-1,Ln-z),min(1,Un+z))}
rownames(HMCI)<-c("A_OAB","A_OAO", "B_OAB", "B_OOB", "A_OAB_monot", "A_OAO_monot")
colnames(HMCI)<-c("Lower Limit", "Upper Limit")
print("H-M CI")
print(HMCI)

#bootstrap B-method CI

BCI<-array(dim=c(6,2))
for (s in 1:6){
Ls<-btemp$t[,2*s-1]
Us<-btemp$t[,2*s]
Ln<-btemp$t0[2*s-1]
Un<-btemp$t0[2*s]
BCI[s,]<-Bci(Ls,Us,Ln,Un,1-conf)
}
rownames(BCI)<-c("A_OAB","A_OAO", "B_OAB", "B_OOB", "A_OAB_monot", "A_OAO_monot")
colnames(BCI)<-c("Lower Limit", "Upper Limit")
print("B-method CI")
print(BCI)