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

library(foreign)

# change to your preferred directory
setwd("~/DocInvest/MyPapers/sqj2012/datasets/umloml")

#dataExperi1spss<-read.spss("http://www.vc.ehu.es/jiwotvim/ecdmOMLUML/data-experi4-spss.sav", to.data.frame=TRUE)
dataExperi1spss<-read.spss("data-experi4-spss.sav", to.data.frame=TRUE)

# recode the number of subjects

subjects <- as.numeric(c(1:length(dataExperi1spss$SUBJECT)))

dat <- data.frame(subj = rep(subjects,each=2),
                       group = rep(dataExperi1spss$GROUP, each=2),
                       seq = rep(dataExperi1spss$GROUP, each=2)
                    )
oldvalues <- c("Group A", "Group B")
newvalues <- c("A","B")  #
dat$seq <- newvalues[match(dat$seq, oldvalues)]

# traspose values for Score Period1 and Period2
dat$Score <- as.vector(
                t(dataExperi1spss[ ,c("NRESPER1","NRESPER2")]))
dat$Time <- as.vector(
  t(dataExperi1spss[ ,c("TSECPER1","TSECPER2")]))

dat$treat <- as.vector(
                t(dataExperi1spss[ ,c("NOTATPE1","NOTATPE2")]))
#dat$treat <- as.factor(dat$treat)

#Period
dat$Period <- ifelse(
  (dat$seq=="A" & dat$treat== "UML")|
    (dat$seq=="B" & dat$treat=="OML"),1,2)

# rename variables for reusing code 
names(dat)[names(dat)=="Time"]<- "TimeInSecs"
names(dat)[names(dat)=="treat"]<- "Formulation"

dat$seq <- ifelse((dat$seq=="A"),1,2)
names(dat)[names(dat)=="seq"]<- "Sequence"
names(dat)[names(dat)=="subj"]<- "Subject"


##################
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["OML"]
ybarSE    <- mmethod["UML"]
# 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  <- .1
# 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")
###########

####################################################################
################################  TIME  ############################
tab.mean <- aggregate(dat$TimeInSecs,list(seq=dat$Sequence,
                                          prd=dat$Period),mean)
# use ``aggregate" to get the variance
tab.var  <- aggregate(dat$TimeInSecs,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)

tapply(dat$TimeInSecs, dat[,c("Sequence","Period")], mean)
### test for carryover
Uik           <- aggregate(dat$TimeInSecs,
                           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$TimeInSecs,
                         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 TimeInSecs by Formulation (method of programming)
mmethod    <- tapply(dat$TimeInSecs, list(method=dat$Formulation), mean)
# extract the means
ybarSEF    <- mmethod["OML"]
ybarSE    <- mmethod["UML"]
# 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
# set above  alphaCI  <- .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 = -10000, upper = 10000, 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")