# bioequivalence from the book of 
###################################################
# Title: Chapter 10: Bioequivalence CT. Author: Dr. DG. Chen 
###################################################
###################################################
#               TRIAL 1
print("Experiment 1")
###################################################
#############   SCORE    ##########################
# use ``aggregate" to the sample size


setwd("~/DocInvest/MyPapers/sqj2012/datasets/razali")

dat <-read.csv("razaliExp1.csv")

names(dat)[names(dat)=="Treatment"]<- "Formulation"
names(dat)[names(dat)=="Efficiency"]<- "Score"

##################
tab.n    <- aggregate(dat$Score, list(seq=dat$Sequence,prd=dat$Period),length)
n1       <- tab.n[tab.n$seq==1 & tab.n$prd==1,]$x
n2       <- tab.n[tab.n$seq==2 & tab.n$prd==1,]$x
n        <- n1+n2
# use ``aggregate" to get the mean
tab.mean <- aggregate(dat$Score,list(seq=dat$Sequence, prd=dat$Period),mean)
# use ``aggregate" to get the variance
tab.var  <- aggregate(dat$Score,list(seq=dat$Sequence, prd=dat$Period),var)
# make a dataframe for the summary data
summaryTab <- data.frame(Sequence=tab.mean$seq,Period=tab.mean$prd,numSample=tab.n$x,Mean=tab.mean$x,Var=tab.var$x)
# print the summary table
round(summaryTab,2)
#attach(dat)

tapply(dat$Score, dat[,c("Sequence","Period")], mean)

### test for carryover
Uik           <- aggregate(dat$Score,
                          list(seq = dat$Sequence,sub=dat$Subject), sum) 
colnames(Uik) <- c("seq", "sub","Uik")
# sample means for each sequence
mUk  = aggregate(Uik$Uik, list(seq=Uik$seq), mean) 
colnames(mUk) <- c("seq", "mUk")
print(mUk)
# differential carryover
hatC  <-  mUk[2,2]-mUk[1,2]
hatC

dU    = merge(Uik, mUk)
sigu2 = sum((dU$Uik-dU$mUk)^2)/(n1+n2-2) 
sigu2

# the test is 
se.sigu = sqrt(sigu2*(1/n1+1/n2)) 
TC      = hatC/se.sigu
TC
# TC should be compared to t(alpha/2, n1+n2-2). if |TC < t H0 of no carryover not rejected |
pC = 2*(1-pt(abs(TC), n1+n2-2))  
pC

#######  test for direct formulation effect but use CIs instead
dik          = aggregate(dat$Score,
                         list(sub=dat$Subject,seq=dat$Sequence),diff)
dik$x        = dik$x/2  
colnames(dik)= c("sub", "seq","dik")
dik
mdk           = aggregate(dik$dik, list(seq=dik$seq), mean) 
colnames(mdk) = c("seq", "mdk")
hatF          = mdk[1,2]-mdk[2,2]
hatF
dF    = merge(dik, mdk)
sigd2 = sum((dF$dik-dF$mdk)^2)/(n1+n2-2) 
sigd2
se.sigd = sqrt(sigd2*(1/n1+1/n2)) 
TF      = hatF/se.sigd
TF
pF = 2*(1-pt(abs(TF), n1+n2-2))  
pF
#####

#########################3  Decision CIs  ###########################
# get the mean for Score by Formulation (method of programming)
mmethod    <- tapply(dat$Score, list(method=dat$Formulation), mean)
# extract the means
ybarSEF    <- mmethod["UML-B"]
ybarSE    <- mmethod["B"]
# make the decision CI. Change here the values for limits
dec2.low = theta.L = -0.2*ybarSE
dec2.up  = theta.U = 0.2*ybarSE
dec25.low <- -0.25*ybarSE
dec25.up <-  0.25*ybarSE
dec5.low <- -0.5*ybarSE
dec5.up <-  0.5*ybarSE
cat("DecisionCI.mean at 20%=(",dec2.low,",",dec2.up,")",sep="","\n")
cat("DecisionCI.mean at 25%=(",dec25.low,",",dec25.up,")",sep="","\n")
cat("DecisionCI.mean at 50%=(",dec5.low,",",dec5.up,")",sep="","\n")

### classical 90% shortest CI for the difference
# the confidence coefficient: alpha
alphaCI  <- 0.05
# the t-value
qt.alpha <- qt(1-alphaCI, n1+n2-2)
qt.alpha
# the lower and upper limits for CI1
low1 <- (ybarSEF-ybarSE)-qt.alpha*sqrt(sigd2)*sqrt(1/n1+1/n2)
up1  <- (ybarSEF-ybarSE)+qt.alpha*sqrt(sigd2)*sqrt(1/n1+1/n2)
cat("The classical CI1=(", round(low1,3),",", 
    round(up1,3),")", sep=" ","\n\n")

##### The Westlake Symmetrical CI
k12 <- 2*(ybarSE-ybarSEF)/sqrt( sigd2*(1/n1+1/n2))

k2 <- uniroot(function(k2) pt(k12-k2,n1+n2-2)- pt(k2,n1+n2-2)
             -(1-alphaCI),lower = -100, upper = 100, tol = 0.0001)$root
k1 =k12-k2
cat("The Westlake k1=",k1," and k2=",k2,sep=" ", "\n\n")
low.west = k2*sqrt(sigd2*(1/n1+1/n2))-(ybarSE-ybarSEF)
up.west  = k1*sqrt(sigd2*(1/n1+1/n2))-(ybarSE-ybarSEF)
cat("The Westlake CI for mu_SEF-mu_SE is 
(",low.west,",",up.west,")",sep=" ", "\n\n")
###########