####################################
##### loading the dataset ##########
####################################


load("c://DataFOS.RData")

library(DoseFinding)
dim(datafos)

##########################
## Data for compound A  ##
##########################
datafos1 <- datafos[datafos$compound=="CompA",]
head(datafos1 )

##########################
## Defining the models ###
##########################

models1 <- list(linear = NULL, linlog = NULL, emax = 1.288150, sigEmax=c( 1.3762727, 0.8783642 ), 
	logistic= c(-8.7298, 2.461841))

#############################################
### Model Averaging using MCPMod function ###
#############################################

dfeA1 <- MCPMod(resp ~ dose,datafos1 , models1, alpha = 0.05, doseEstPar = 0.5, pVal = TRUE,
           selModel = "aveAIC", doseEst = "ED", off = 1,critV = TRUE, optimizer = "bndnls", lenDose=1000)

summary(dfeA1 )

##############################################
## Obtaining Bootstrap CI for the estimates ##
## Note: it may take a while for processing ##
##############################################

resA1 <- bootMCPMod(dfeA1 , nSim = 1000, seed=123)
resA1 

## !!! Note: load first all the required      ##
## functions provided at the end of this file ##



#################
## Fig 13.4(a)###
#################

PlotModel (data=datafos1, toplot="each")
#PlotModel (data=datafos1, toplot="both")

#################
## Fig 13.4(b)###
#################

PlotModel (data=datafos1, toplot="MA")
doseSeq <-  seq(min(datafos1$dose),max(datafos1$dose),0.001)
bound <- getBnds(mD = (max(datafos1$dose)))
preddat <- data.frame(dose = doseSeq  )
Emax <- fitDRModel(resp ~ dose,datafos1, "emax",  optimizer = "bndnls",bnd=bound$emax)
predHEmax<- predict(Emax , type = "fullModel",newdata=preddat )
lines(doseSeq ,predHEmax,lwd=2, lty=2)
legend("bottomright",c("Averaged Model", "Hyp Emax"), col=c(2,1), lty=c(1,2),lwd=2, merge = TRUE )


#################
## Fig 13.5(a)###
#################

########################################
### Perform Bootstrap for FOS Comp A  ###
########################################

## Bootstrap Index  ##
set.seed(123)
indArray <- indGen(nSim=1000, datafos1$dose)

#################################
## Bootstrap, may take a while ##
#################################

boot <- t(apply(indArray, 1, CalcEDeach , data = datafos1 , model=models1 ,doseEstPar =0.5, doseEst ="ED" ))

#####################################
## Estimate from the original data ##
#####################################

dfeA1 <- MCPMod(resp ~ dose,datafos1 , models1, alpha = 0.05, doseEstPar = 0.5, pVal = TRUE,
           selModel = "aveAIC", doseEst = "ED", off = 1,critV = TRUE, optimizer = "bndnls", lenDose=1000)
estEdMod <- data.frame(attr(dfeA1 $estDose,"tdModels"))

MA.res    <- data.frame(lo=quantile(boot[,1] , c(0.025), na.rm=TRUE) , ED50= dfeA1$estDose , hi=quantile(boot[,1] , c(0.975), na.rm=TRUE))
Lin.res   <-  data.frame(lo=quantile(boot[,3] , c(0.025), na.rm=TRUE) , ED50= estEdMod $linear, hi=quantile(boot[,3] , c(0.975), na.rm=TRUE))
LinLog.res<-  data.frame(lo=quantile(boot[,4] , c(0.025), na.rm=TRUE) , ED50= estEdMod $linlog, hi=quantile(boot[,4] , c(0.975), na.rm=TRUE))
HEmax.res <-  data.frame(lo=quantile(boot[,2] , c(0.025), na.rm=TRUE) , ED50= estEdMod $emax, hi=quantile(boot[,2] , c(0.975), na.rm=TRUE))
SEmax.res <-  data.frame(lo=quantile(boot[,6] , c(0.025), na.rm=TRUE) , ED50= estEdMod $sigEmax,  hi=quantile(boot[,6] , c(0.975), na.rm=TRUE))
Logis.res <-  data.frame(lo=quantile(boot[,5] , c(0.025), na.rm=TRUE) , ED50= estEdMod $logistic, hi=quantile(boot[,5] , c(0.975), na.rm=TRUE))

	
Result.Fos.compA <- rbind(Lin.res, LinLog.res,HEmax.res, SEmax.res, Logis.res,MA.res   )
rownames(Result.Fos.compA ) <- c("Lin","LinLog","HEmax","SEmax","Logis","MA")

MAPlot2 (Result.Fos.compA ,title="",xlab="ED50")


#################
## Fig 13.5(b)###
#################


#################################
### ED50 for Model average ######
#################################

## calculate MA Prediction ##
 doseSeq <-  seq(min(datafos1$dose),max(datafos1$dose),0.001)
 PredFos <- FitModels (datafos1, doseSeq=doseSeq, model=models1)
 FosMA <- cbind(doseSeq ,PredFos $PredMA)


## Predict the Y value for ED50 and the low and high bound ###
loED50 = quantile(boot[,1] , c(0.025), na.rm=TRUE)
YloED50 <- min(FosMA [which (FosMA[,1] >= loED50  ),2])

ED50 <- dfeA1 $estDose
YED50 <- min(FosMA [which (FosMA[,1] >= ED50),2])

hiED50 = quantile(boot[,1] , c(0.975), na.rm=TRUE)
YhiED50 <- min(FosMA [which (FosMA[,1] >= hiED50),2])

dd1 <- data.frame (c( loED50 , loED50 ,0),c(0,YloED50 ,YloED50 ))
dd2 <- data.frame (c( ED50, ED50,0),c(0,YED50 ,YED50 ))
dd3 <- data.frame (c( hiED50 , hiED50 ,0),c(0,YhiED50 ,YhiED50 ))


plot(datafos1$dose, datafos1$resp, xlab = "dose", ylab = "gene expression",lwd=2 , main="Model Average")

## The predicted MA ##
lines(FosMA ,lwd=2, lty=1)
lines(dd1, lty=2)
lines(dd2)
lines(dd3, lty=2)
points(ED50,YED50 ,col="red",pch=19,lwd=2)
abline(h=0, v=0)


########################################
## Comparison Between DR-curve for a  ##
## single gene with several compounds ##
########################################


library(nlme)
miny <- min(datafos$resp)
rng <- range(datafos$resp)
dif <- diff(rng )
#############  

#######################################
## Model 0:  Without Compound Effect ##
#######################################

fos.Com0 <- gnls(resp  ~ E0 + ( (dose *Emax)/ (xmid+dose)),data = datafos,
  start = c(Emax = dif,E0 = miny ,xmid = 0.63))

summary(fos.Com0)

###################################
##  Model 1: With the common E0  ##
###################################

fos.Com1 <- gnls(resp  ~ E0 + ( (dose *Emax)/ (xmid+dose)),data = datafos , 
  params= list(Emax+xmid~compound,E0~1),
  start = c(Emax = dif,E0 = miny, rep(0,5), xmid = 0.63, rep(0,5)))

summary(fos.Com1)

#######################################
##  Model 2: With common E0 and Emax ##
#######################################

fos.Com2 <- gnls(resp  ~ E0 + ( (dose*Emax)/ (xmid+dose)),data = datafos , 
  params= list(Emax+E0~1,xmid~compound),
  start = c(Emax = dif,E0 = miny, xmid = 0.63,rep(0,5)))

summary(fos.Com2)


######################################################################
## Model 3: Full Model with compound spesific on E0, Emax, and xmid ##
######################################################################


fos.Com3 <- gnls(resp  ~ E0 + ( (dose *Emax)/ (xmid+dose)),data = datafos, 
  params=E0 + Emax+xmid~compound, start = c(Emax = dif, rep(0,5),E0 = miny, 
  rep(0,5),xmid = 0.63,rep(0,5)))

summary(fos.Com3 )

################################
## anova test for all models  ##
################################

anova(fos.Com0, fos.Com1 , fos.Com2, fos.Com3)

### plot the predicted values ###
## Model with common E0 ##
## Fig 13.7 (a) Extrapolating the dose ###
plotextrap (datafos,fos.Com1 )

## Fig 13.7 (c) within the dose range ###
plotextrap (datafos,fos.Com1,extrap=F )

## Model with common E0 and Emax ##

## Fig 13.7 (b) Extrapolating the dose ###

plotextrap (datafos,fos.Com2  )

## Fig 13.7 (d) Extrapolating the dose ###

plotextrap (datafos,fos.Com2,extrap=F  )


######################
### MCPmod vs gnls ###
######################

###  With E0 * Compound using gnls() function ###

fos.Com4 <- gnls(resp  ~ E0 + ( (dose *Emax)/ (xmid+dose)),data = datafos , 
params= list(E0~compound,Emax+xmid~1),
start = c(Emax = 2.3851944,rep(0,5), E0 = 9.5487984,xmid = 0.63))

summary(fos.Com4)

###  With E0 * Compound  using MCPMod() function ###

modelEmax <- list( emax = 1.288150)

dfeB <- MCPMod(resp ~ dose,datafos,modelEmax ,
                 addCovars=~compound,alpha=0.05,
                 doseEstPar = 0.5,pVal=TRUE,
                 selModel = "aveAIC",
                 doseEst = "ED",critV = TRUE,
                 optimizer = "bndnls", lenDose=1000)
summary(dfeB )

###  With E0 * Compound  using fitDRModel() function ###

bound <- getBnds(mD = 40)
fit <- fitDRModel(resp ~ dose,datafos, model ="emax",
                    optimizer="bndnls",
                    bnd=bound$emax,
                    addCovars=~as.factor(compound))

fit 




####################################
### Additional Functions needed ####
####################################



################################################
### Function for plotting the fitted models ####
################################################

PlotModel <- function (data, toplot=c("each", "MA", "both")) {
      toplot<- match.arg(toplot)
      doseSeq <-  seq(min(data$dose),max(data$dose),0.001)
      preddat <- data.frame(dose = doseSeq  )

      linear <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "linear"),error= function (x) x="error")
  	if (linear !="error"){
      predLin <- predict(linear , type = "fullModel",newdata=preddat )
  	}
      ## obtain the bound for emax, sigmoid emax and logistic model ##
      bound <- getBnds(mD = 40)
      
      ## Hyperbolic Emax model ##

	Emax <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "emax",  optimizer = "bndnls",bnd=bound$emax), error= function (x) x="error")

	if (Emax !="error")  predHEmax <- predict(Emax , type = "fullModel",newdata=preddat )

      ## Sigmoidal Emax model ##

	EmaxSig <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "sigEmax", optimizer = "bndnls",bnd=bound$sigEmax), error= function (x) x="error")

	if (EmaxSig !="error")	predSEmax<- predict(fit , type = "fullModel",newdata=preddat )

	fpl2 <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "logistic", optimizer = "bndnls",  bnd=bound$logistic ) ,error= function (err) err="error")

	if (fpl2 !="error")	predLogis<- predict(fit , type = "fullModel",newdata=preddat )
      
      loglinear <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "linlog",off = 1),error= function (x) x="error")
      if (loglinear !="error")       predLinLog <- predict(loglinear , type = "fullModel",newdata=preddat )
   	
      if (toplot=="MA" | toplot=="both") {

        models1 <- list(linear = NULL, linlog = NULL, emax = 1.288150, sigEmax=c( 1.3762727, 0.8783642 ), logistic= c(-8.7298, 2.461841))

        dfeA1 <- MCPMod(resp ~ dose,data, models1, alpha = 0.05, doseEstPar = 0.5, selModel = "aveAIC", doseEst = "ED", 
          off = 1,critV = TRUE, optimizer = "bndnls", lenDose=1000)

        ## Obtaining the weights ##
        weight <- attr(dfeA1$fm,"weights")

        ## Obtaining the Model averaged predicted  values ##
        ## Predicted values of all models are multiplied by the weights ##

        predict_all <- data.frame(predSEmax, predHEmax, predLogis, predLinLog,predLin )
        PredMA  <- as.matrix(predict_all )%*% as.vector(weight )
       }


    if (toplot=="MA") {
       plot(data$dose, data$resp, xlab="doses", ylab="gene expression")
       lines(preddat[,1] ,PredMA ,col="red" ,lwd=2)
      }

    if (toplot=="each"|toplot=="both") {
    	plot(data$dose, data$resp, xlab = "dose", ylab = "gene expression",lwd=2)
    	try(lines(doseSeq ,predLin ,lwd=2, lty=1),silent =T)
      try(lines(doseSeq ,predHEmax,lwd=2, lty=3),silent =T)
	try(lines(doseSeq ,predSEmax,lwd=2, lty=4),silent =T)
	try(lines(doseSeq ,predLogis,lwd=2, lty=5),silent =T)
      try(lines(doseSeq ,predLinLog , lwd=2, lty=2),silent =T)
      }

   if (toplot=="each")      legend("bottomright",c("Linear","Linear Log","Hyp Emax","Sig Emax","Logistic"), lty=c(1:5),lwd=2, merge = TRUE )
  
   if (toplot=="both")    { 
  
     lines(preddat[,1] ,PredMA ,col="red" ,lwd=2)
     legend("bottomright",c("Linear","Linear Log","Hyp Emax","Sig Emax","Logistic","Averaged Model"), col=c(1,1,1,1,1,2), lty=c(1:5),lwd=2, merge = TRUE )
    }
}



MAPlot2  <- function (result, title=NULL, ylab=NULL,xlab) 

  {
if (is.null(ylab)) ylab="Model"
maxV <- max(result[,3] )

plot(result[,2],1:nrow(result),axes = FALSE, xlab=xlab, ylab=ylab, xlim= c(0,maxV ), main=title)
 axis(2, 1:nrow(result), rownames(result))
        axis(1)
        box()

    for (i in 1:nrow(result)) {
      arrows( result[i,1], i,  result[i,3],i, length = 0.1, 
             code = 3, col = par("fg"), lty = NULL, 
            lwd = par("lwd"), xpd = NULL)
      }

}
########################################
## obtaining permutations matrix      ##
## extracted from DoseFinding Package ##
########################################
indGen <- function (nSim, doseVec) 
{
    nObs <- length(doseVec)
    output <- matrix(nrow = nSim, ncol = nObs)
    inds <- sort(unique(doseVec))
    for (is in inds) {
        gp <- (1:nObs)[doseVec == is]
        sampl <- sample(gp, nSim * length(gp), replace = TRUE)
        output[, gp] <- matrix(sampl, nrow = nSim)
    }
    output
}



#### Calculating IC50/ED50 ####

ICdose  <- function (pred,doseSeq  ) {
 pEff <- 0.5 * (max(pred) - pred[1])
 ind <- pred >= pEff + pred[1]
  val <- getDose(doseSeq, ind)
  return(val)
   }


########################################
## Get the dose in grid search        ##
## extracted from DoseFinding Package ##
########################################


getDose <- function (dose, ind) 
{
    aa <- !is.na(ind)
    if (!all(aa)) {
        ind <- ind[aa]
        dose <- dose[aa]
    }
    if (any(ind)) 
        min(dose[ind])
    else NA
}



##############################
## Calculate ED50 for each  ##
## model for bootstrap      ##
##############################

CalcEDeach <- function(index, data,model,doseEstPar, doseEst  ){

 data <- data [index,]
 dfeA1 <- MCPMod(resp ~ dose,data, model, alpha = 0.05, doseEstPar = doseEstPar , selModel = "aveAIC", doseEst = doseEst , 
  off = 1,critV = TRUE, optimizer = "bndnls", lenDose=1000)

estEdMod <- attr(dfeA1 $estDose,"tdModels")
ResED50 <-  c(dfeA1 $estDose, estEdMod [,order(colnames(estEdMod ) )] )
ResED50 
}

### Fit the model and predict the MA fitted values ###



FitModels <- function (data, doseSeq, model ) {
      preddat <- data.frame(dose = doseSeq  )

      linear <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "linear"),error= function (x) x="error")
  	if (linear !="error")      predLin <- predict(linear , type = "fullModel",newdata=preddat )
  	
      ## obtain the bound for emax, sigmoid emax and logistic model ##
      bound <- getBnds(mD = (max(data$dose)))
      
      ## Hyperbolic Emax model ##

	Emax <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "emax",  optimizer = "bndnls",bnd=bound$emax), error= function (x) x="error")

	if (Emax !="error")  predHEmax <- predict(Emax , type = "fullModel",newdata=preddat )

      ## Sigmoidal Emax model ##

	EmaxSig <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "sigEmax", optimizer = "bndnls",bnd=bound$sigEmax), error= function (x) x="error")

	if (EmaxSig !="error")	predSEmax<- predict(fit , type = "fullModel",newdata=preddat )

	fpl2 <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "logistic", optimizer = "bndnls",  bnd=bound$logistic ) ,error= function (err) err="error")

	if (fpl2 !="error")	predLogis<- predict(fit , type = "fullModel",newdata=preddat )
      
      loglinear <- tryCatch(fit <- fitDRModel(resp ~ dose,data, "linlog",off = 1),error= function (x) x="error")
      if (loglinear !="error")       predLinLog <- predict(loglinear , type = "fullModel",newdata=preddat )
   	
      dfeA1 <- MCPMod(resp ~ dose,data, model, alpha = 0.05, doseEstPar = 0.5, selModel = "aveAIC", doseEst = "ED", 
          off = 1,critV = TRUE, optimizer = "bndnls", lenDose=1000)

        ## Obtaining the weights ##

        weight <- attr(dfeA1$fm,"weights")

        ## Obtaining the Model averaged predicted  values ##
        ## Predicted values of all models are multiplied by the weights ##

        predict_all <- data.frame(predSEmax, predHEmax, predLogis, predLinLog,predLin )
        PredMA  <- as.matrix(predict_all )%*% as.vector(weight )
      
   Result <- data.frame(PredMA  ,predict_all)
   Result
}


########################################
## Function to plot the fitted values ## 
## of model with compound effects     ## 
########################################
compName <- c("CompA","JnJa","CompB","CompC","JnJb","JnJc")

plotextrap  <- function (dat,fit, extrap=TRUE, legend =TRUE) {

  dose2 <- dat$dose 
  dose2 [dose2 == 0] <- 0.04

  dose= sort(unique(dat$dose))
  doses <- sort(unique(log10(dose2 )))

   if (extrap)  plot(log10(dose2 ),dat$resp, xlim=c(-2,4), ylim=c(8,13),xlab="doses", ylab="Gene expression",axes = FALSE,frame.plot=T)
   else {
     plot(log10(dose2 ),dat$resp,,xlab="doses", ylab="Gene expression",axes = FALSE,frame.plot=T)
     } 
  axis(1, doses  , as.character(dose))
  axis(2)
   
  x1 <- seq (10,10000,10)
  x1a <-seq ( 0.04,10,0.1) 
  
  x1b <- c(x1a, x1)
  x1b [x1b ==  0.04] <- 0
  nsamp <- length(x1b)
  x2 <- data.frame(dose=rep(x1b ,6 ),compound=rep(compName , each=nsamp ) )
  x3 <- predict(fit ,x2)
  x11 <- x1b 
  x11 [x1b == 0] <- 0.04
  lines(log10(x11 ),x3 [1:nsamp ], lty=1, lwd=2)
  lines(log10(x11 ),x3 [(nsamp +1): (2*nsamp) ], lty=2,lwd=2)
  lines(log10(x11 ),x3 [(2*nsamp +1): (3*nsamp) ], lty=3,lwd=2)
  lines(log10(x11 ),x3 [(3*nsamp +1): (4*nsamp) ], lty=4,lwd=2)
  lines(log10(x11 ),x3 [(4*nsamp +1): (5*nsamp) ], lty=5,lwd=2)
  lines(log10(x11 ),x3 [(5*nsamp +1): (6*nsamp) ], lty=6,lwd=2)

  if (legend){
      if (extrap)  legend(2,9.5,compName , lty=c(1:6),lwd=2, merge = TRUE )

     else {
        legend(0.7,9.7,compName , lty=c(1:6),lwd=2, merge = TRUE ) 
        }

    }

}