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

setwd("~/DocInvest/MyPapers/sqj2012/datasets/razali")
dat <-read.csv("razaliExp2.csv")

# Recode N/A as  NA
NAs <- dat == "N/A"
dat[NAs] <- NA #   for overall comprehension delete all rows with N/As
# find subjects with NA
subjects.with.na <- as.vector(dat$SUBJECTS[is.na(dat$Question.1)])
# select only the subjects without NA
dat <- dat[which(!(dat$SUBJECTS %in%  subjects.with.na)),]

#-> -> Compute Overall Comprehension (and then Retention and Transfer) Efficiency- and Modification 
# as the sum of the SIX questions

dat$Comprehension <- 
  as.numeric(as.vector(dat$Question.1))

names(dat)[names(dat)=="Treatment"]<- "Formulation"
names(dat)[names(dat)=="Comprehension"]<- "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["Event-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")
###########
######################   Figures EFFICIENCY
# actual values. Add this script to the bioeq....R
ci1score <- c(low1, up1)


## from other R scripts
lowdelta20 <- dec2.low
lowdelta25 <- dec25.low
lowdelta50 <- dec5.low
updelta20 <-  dec2.up
updelta25 <-  dec25.up
updelta50 <-  dec5.up


limxscore <- 0.7
namethisfile <- paste(getwd(),"/razaliFigQuestion1-2Exp",alphaCI,".pdf",sep="")
pdf(file=namethisfile, height=1.5, width=4.5) #, paper="a4",pagecentre=TRUE) #, pagecentre=TRUE)

par(mar=c(3,0,0,0))
plot(0,0,type="n", xlim=c(lowdelta20,limxscore), ylim=c(0,1.5),axes=FALSE, ann=FALSE)

equivintervals <- c(limxscore,lowdelta20,0,updelta20,updelta50)
equivintervals <- sort(equivintervals)

labels <- c(
  #           expression(-Delta[paste(50,'%',sep="")]),
  expression(-Delta[paste(20,'%',sep="")]),
  0,
  expression(Delta[paste(20,'%',sep="")]),
  expression(Delta[paste(50,'%',sep="")]),
  limxscore)
axis(1, at=equivintervals, labels = labels, pos=0)
mtext(expression(mu[T]-mu[R]),at=0,side=1,line=2)
mtext("More Efficiency in UML-B",at=updelta50,side=1,line=2)

#abline(v=equivintervals,lty="dotted")
segments(c(lowdelta50,updelta50),c(0,0),c(lowdelta50,updelta50),c(1.5,1.5), lty=3,lwd=1)
segments(c(lowdelta20,updelta20),c(0,0),c(lowdelta20,updelta20),c(1,1), lty=3,lwd=1)

segments(ci1score[1],0.2,ci1score[2],lty=1,lwd=5, lend=2)
text(ci1score[1],0.2,"|",offset=0.5,cex=1.1,adj=1,font=2)
text(ci1score[2],0.2,"|",offset=0.5,cex=1.1,adj=0,font=2)


text(0.2,0.2,
     "Exp.2 Question1",
     offset=0.5,cex=0.9,adj=1)
text(0.4,0.2,
     paste("(",
           round(ci1score[1],2),",",
           round(ci1score[2],2),
           ")"),
     pos=3,offset=0.5,cex=0.8,adj=1)

dev.off()  # for pdf et al