

library(shiny)
#library(shinyjs)
#library(shinyAce)
library(shinydashboard)
#library(shinycssloaders)
library(purrr)
library(gslnls)
library(tidyverse)
library(ggplot2)
library(reshape2)
library(openxlsx)
library(DT)
library(ggpubr)
library(gridExtra)
library(drc)
library(twopartm)
library(car)
library(dplyr)


Dat <- reactiveValues()
REP <- reactiveValues()

dilFUN2 <- function(cs_,dils,Faktor) {
  av <- cs_
  dils_av <- dils_av
  dils_avsc <- dils_av*Faktor
  dils2 <- dils_avsc+av
  dilfactors <- 1/exp(dils2-lag(dils2))
  return(dilfactors)
}

plot_f <- function(dat, sigmoid,det_sig) {
  CORdat <- cor(dat[,1],dat[,ncol(dat)])
  #browser() 
  all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout")
  isRef <- rep(c(1,0),1,each=nrow(all_l)/2) 
  isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
  all_l2 <- cbind(all_l, isRef, isSample)
  #browser()  
  if(is.null(det_sig)) {
    if (CORdat<0) {
      startlist <- list(a=sigmoid[1], b=-sigmoid[5],cs=sigmoid[7],
                        d=sigmoid[3],r=sigmoid[8])
    } else {
      
      startlist <- list(a=sigmoid[1],b=sigmoid[5],cs=sigmoid[7],
                        d=sigmoid[3],r=sigmoid[8])
    }
  } else {
    startlist <- list(a=det_sig[5], b=det_sig[1],cs=det_sig[7],
                      d=det_sig[3],r=det_sig[7] - det_sig[8])
  }
  
  mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(log_dose-(cs-r*isSample)))),
                data=all_l,
                start=startlist,
                control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  s_mr <- tryCatch({
    s_mr <- summary(mr)
  },
  error = function(err) {
    s_mr <- NULL
  })
  
  a <- s_mr$coefficients[1,1]
  b <- s_mr$coefficients[2,1]
  cs <- s_mr$coefficients[3,1]
  d <- s_mr$coefficients[4,1]
  r <- s_mr$coefficients[5,1]
  
  log_dose <- unique(all_l$log_dose)
  seq_x <- seq(min(log_dose),max(log_dose),0.1)
  SAMPLE <- a+(d-a)/(1+exp(b*(seq_x-(cs-r))))
  REF <- a+(d-a)/(1+exp(b*(seq_x-(cs))))
  
  if (is.null(det_sig)) {
    SAMPLEtrue <- sigmoid[2] + (sigmoid[4] -sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-sigmoid[8]-seq_x))))
    REFtrue <- sigmoid[1] + (sigmoid[3] -sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x))))
  } else {
    SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(det_sig[2]*(det_sig[8]-seq_x)))
    REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-seq_x)))
  }
  #browser()
  pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue)
  all_l2$readout[all_l2$readout < 0] <- 0.01
  all_l2$readouttrans <- log(all_l2$readout)
  slopeEC50 <- b*(a-d)/4
  
  Xbendl3 <- cs-(1.31696/b)
  Xbendu3 <- cs+(1.31696/b)
  XbendlT <- cs-r-(1.31696/b)
  XbenduT <- cs-r+(1.31696/b)
  bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3),
                  bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,3))
  Dat$bendpoints <- bendpoints
  Dat$cfordils <- cs
  
  p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
    geom_point(shape=factor(isRef), alpha=0.8) +
    labs(title = paste("restricted 4pl; bendp:", round(Xbendl3,3),round(Xbendu3,3),round(XbendlT,3),round(XbenduT,3)),
         color="product") +
    scale_color_manual(labels=c("test","reference"), values=c("red","blue")) +
    scale_shape_manual(labels=c("test","reference")) +
    theme_bw() +
    theme(axis.text = element_text(size=14))
  
  p2 <- p + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLE), color="red",
                      inherit.aes = F) +
    geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="blue",
              inherit.aes = F) +
    geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="red", linetype=2, alpha=0.4,
              inherit.aes = F) +
    geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="blue", linetype=2, alpha=0.4,
              inherit.aes = F) +
    geom_vline(xintercept=c(Xbendl3, Xbendu3), col="blue",linetype=2) +
    geom_vline(xintercept=c(XbendlT, XbenduT), col="red",linetype=2) +
    annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) +
    theme(legend.position="none")
  Dat$p2 <- p2
  
  # transformed plots
  p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) +
    geom_point(shape=factor(isRef), alpha=0.8) +
    labs(title = paste("restricted transformed 4pl"), color="product") +
    scale_color_manual(labels=c("test","reference"), values=c("red","blue")) +
    theme_bw()
  
  mrt <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(log_dose-(cs-r*isSample)))),
                 data=all_l2,
                 start=startlist,
                 control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  s_mrt <- summary(mrt)
  a_trans <- s_mrt$coefficients[1,1]
  b_trans <- s_mrt$coefficients[2,1]
  cs_trans <- s_mrt$coefficients[3,1]
  d_trans <- s_mrt$coefficients[4,1]
  r_trans <- s_mrt$coefficients[5,1]
  
  XbendlTrans <- cs_trans-(1.31696/b_trans)
  XbenduTrans <- cs_trans+(1.31696/b_trans)
  XbendlTransT <- cs_trans-r_trans-(1.31696/b_trans)
  XbenduTransT <- cs_trans-r_trans+(1.31696/b_trans)
  bendpointsTRANS <- c(bendREF_lower = round(XbendlTrans,3), bendREF_upper=round(XbenduTrans,3),
                       bendSAMPLE_lower = round(XbendlTransT,3), bendSAMPLE_upper=round(XbenduTransT,3))
  Dat$bendpointsTRANS <- bendpointsTRANS
  SAMPLEtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*(seq_x-(cs_trans-r_trans))))
  REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*(seq_x-(cs_trans))))
  
  pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans)
  p_rt2 <- p_rt + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=SAMPLEtrans), color="red",
                            inherit.aes = F) +
    geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="blue",
              inherit.aes = F) +
    geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="blue",linetype=2) +
    geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="red",linetype=2) +
    theme(legend.position = "none", axis.text=element_text(size=14))
  
  if (is.null(det_sig)) {
    unrestr <- drm(readout ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(),
                   pmodels=data.frame(isSample, isSample,isSample,isSample))
    Sum_u <- summary(unrestr)
    ast <- Sum_u$coefficients[3,1]
    ate <- Sum_u$coefficients[4,1]
    bst <- Sum_u$coefficients[1,1]
    bte <- Sum_u$coefficients[2,1]
    cst <- log(Sum_u$coefficients[7,1])
    cte <- log(Sum_u$coefficients[8,1])
    dst <- Sum_u$coefficients[5,1]
    dte <- Sum_u$coefficients[6,1]
  } else {
    ast <- det_sig[5]
    ate <- det_sig[6]
    bst <- det_sig[1]
    bte <- det_sig[2]
    cst <- det_sig[7]
    cte <- det_sig[8]
    dst <- det_sig[3]
    dte <- det_sig[4]
  }
  REFu <- ast + (dst-ast)/(1+exp(bst*(seq_x-cst)))
  SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(seq_x-cte)))
  pl_df2 <- cbind(seq_x, SAMPLEu, REFu)
  #browser()  
  pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
    geom_point() +
    labs(title="unrestricted 4_pl-Model", color="product") +
    scale_color_manual(labels = c("test","reference"), values=c("red","blue")) +
    theme_bw()
  pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu),
                        color="red", inherit.aes = F) +
    geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu),
              color="blue", inherit.aes = F,
              show.legend = F) 
  pu2_ <- pu2 +
    theme(legend.position = "none", axis.text = element_text(size=14))
  putrans <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) +
    geom_point() +
    labs(title="unrestricted transformed 4_pl-Model", color="product") +
    scale_color_manual(labels = c("test","reference"), values=c("red","blue")) +
    theme_bw()
  
  unrestr_trans <- drm(readouttrans ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(),
                       pmodels=data.frame(isSample, isSample,isSample,isSample))
  Sum_ut <- summary(unrestr_trans)
  ast_t <- Sum_ut$coefficients[3,1]
  ate_t <- Sum_ut$coefficients[4,1]
  bst_t <- Sum_ut$coefficients[1,1]
  bte_t <- Sum_ut$coefficients[2,1]
  cst_t <- log(Sum_ut$coefficients[7,1])
  cte_t <- log(Sum_ut$coefficients[8,1])
  dst_t <- Sum_ut$coefficients[5,1]
  dte_t <- Sum_ut$coefficients[6,1]
  
  REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(seq_x-cst_t)))
  SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(seq_x-cte_t)))
  pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans)
  
  pu2_t <- putrans + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=SAMPLEu_trans),
                               color="red", inherit.aes = F) +
    geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans),
              color="blue", inherit.aes = F,
              show.legend = F) 
  pu3_t <- pu2_t
  grid.arrange(p2,p_rt2,pu2_,pu3_t, nrow=2)
}

Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1,  dt=10,r=0.0001,ct=cs-r,
                        sd_fac=0.1, gt=1, gs=1, log_conc,
                        heteroNoise=FALSE, noDilSeries, noDils) {
  yAxfac <- (ds-as)
  log_dose <- log_conc
  isRef <- rep(c(1,0),1,each=length(log_conc)*noDilSeries)
  isSample <- rep(c(0,1),1,each=length(log_conc)*noDilSeries)
  #browser()  
  av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
    (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose)))
  if (heteroNoise) {
    # heterosc noise
    ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac/100))), nrow=noDils, ncol=noDilSeries*2)
  } else {
    # homosc noise
    ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac*yAxfac/100))), nrow=noDils, ncol=noDilSeries*2)
  }
  ro_jit <- abs(ro_jit)
  
  ro_jit2 <- cbind(ro_jit, log_dose)
  if (noDilSeries==3) {
    colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose")
  } else {
    colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose")
  }
  return(ro_jit2)
}

LinPotTab <- function(circles, Lim, PureErrFlag) {
  circ_ABl <- circles
  circ_Al <- circ_ABl[circ_ABl$isSample ==1,]
  circ_Al <- circ_ABl[circ_ABl$isSample ==0,]
  # restr CSSI model
  modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
  coeffs <- modAB$coefficients
  SU_modAB <- tryCatch({
    SU_modAB <- summary(modAB)
  }, error = function(msg) {
    return(NA)
  })
  # Intercept diff/slope modAB
  linPot <- exp(modAB$coefficients[3]/modAB$coefficients[2])
  
  if(PureErrFlag) {
    FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
    meanPureErr <- FitAnova[4,3]
    DFsPure <- FitAnova[4,1]
    VCOV <- vcov(modAB)
    V_V <- VCOV/SU_modAB$sigma^2
    VCOVpure <- V_V*meanPureErr
    SEsPure <- sqrt(diag(V_V)*meanPureErr)
  }
  
  log_pot_delta <- deltaMethod(modAB, "isSample/log_dose")
  if (PureErrFlag) {
    V_ <- log_pot_delta$SE^2/SU_modAB$sigma^2
    V_p <- V_*meanPureErr
    potDeltaPureSE <- sqrt(V_p)
    CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure)*potDeltaPureSE
    CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure)*potDeltaPureSE
  } else {
    CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE
    CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE
  }
  #browser()
  ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))
  if (ExpLinPot[2]*100>Lim[[9]] & ExpLinPot[3]*100<Lim[[10]]) test_potCI <- 0 else test_potCI <- 1
  # Rel pot CI
  relLinpotCI <- ExpLinPot/linPot*100
  pottab <- cbind(round(linPot*100,3), round(ExpLinPot[2]*100,3), round(ExpLinPot[3]*100,3),
                  round(test_potCI,3), round(relLinpotCI[2],3),round(relLinpotCI[3],3))
  colnames(pottab) <- c("Potency","lower 95%CI", "upper 95%CI", "test_result", "lowerRel95%CI","upperRel95%CI")
  return(pottab)
}

ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1,0),1,each=nrow(all_l)/2)
  isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_lA <- all_l[all_l$isSample == 1,] # TEST
  all_lB <- all_l[all_l$isSample == 0,] # REF
  
  circ_ABl <- circles
  circ_Al <- circ_ABl[circ_ABl$isSample ==1,]
  circ_Bl <- circ_ABl[circ_ABl$isSample ==0,]
  # restr CSSI model
  modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
  # unrestr with interact term SSSI
  modABu <- lm(readout ~ log_dose + isSample + log_dose*isSample, circ_ABl)
  modA <- lm(readout ~ log_dose + isSample, circ_Al)
  modB <- lm(readout ~ log_dose + isSample, circ_Bl)
  
  log_pot_delta <- deltaMethod(modAB, "isSample/log_dose")
  CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE
  CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE
  ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))
  if (ExpLinPot[2]*100>Lim[9] & ExpLinPot[3]*100>Lim[10]) test_potCI <- 0 else test_potCI <- 1
  
  su_mod <- summary(modAB)$coefficients
  su_mod2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.")), su_mod)
  su_modU <- summary(modABu)$coefficients
  su_modU2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.","slope difference")), su_modU)
  
  uCI_SloDiff <- su_modU[4,1] + qt(0.975,8)*su_modU[4,2]
  lCI_SloDiff <- su_modU[4,1] - qt(0.975,8)*su_modU[4,2]
  SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiff,uCI_SloDiff)
  
  lenCirc <- nrow(circ_ABl)
  dfTreat <- 3
  dfPrep <- 1
  dfReg <- 1
  dfnonP <- 1
  dfRMSE <- c(lenCirc-3-1)
  dfTotal <- lenCirc-1
  dfPureE <- lenCirc-6
  dfNonLin <- dfRMSE-dfPureE
  
  RSS <- sum(resid(lm(readout ~ log_dose*isSample, circ_ABl))^2)
  MSE <-  RSS/dfRMSE
  SSE <-  sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)
  MSpure <- SSE/dfPureE
  
  if (PureErrFlag) {
    FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
    meanPureErr <- FitAnova[4,3]
    SU_modAB <- tryCatch({
      SU_modAB <- summary(modAB)
    }, error = function(msg) {
      return(NA)
    })
    if (length(SU_modAB)>1) s_modABcoeffs <- summary(modAB)$coefficients
    
    DFsPure <- FitAnova[4,1]
    VCOV <- vcov(modAB)
    V_V <- VCOV/SU_modAB$sigma^2
    VCOVpure <- V_V*meanPureErr
    SEsPure <- sqrt(diag(V_V)*meanPureErr)
    su_mod2[,3] <-  SEsPure
    su_mod2[,4] <-  su_mod2[,2]/su_mod2[,3]
    su_mod2[,5] <-  2*(1-pt(abs(su_mod[,4]), FitAnova[4,1]))
    
    s_mu <- summary(modABu)$coefficients
    SU_modABu <- summary(modABu)
    VCOVu <- vcov(modABu)
    V_Vu <- VCOVu/SU_modABu$sigma^2
    SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
    
    su_modU2[,3] <-  SEsPureU
    su_modU2[,4] <-  su_modU2[,2]/su_modU2[,3]
    su_modU2[,5] <-  2*(1-pt(abs(su_modU2[,4]), FitAnova[4,1]))
    
    uCI_SloDiffP <- su_modU[4,1] + qt(0.975,8)*SEsPureU[4]
    lCI_SloDiffP <- su_modU[4,1] - qt(0.975,8)*SEsPureU[4]
    SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiffP,uCI_SloDiffP)
    
    SSRes <- SSE
    dfRes <- dfPureE
    
  } else {
    SSRes <- RSS
    dfRes <- dfRMSE
  }
  
  # treatment
  SStreat <- print(sum((predict(lm(readout ~ factor(log_dose)*isSample, circ_ABl))-mean(circ_ABl$readout))^2))
  F_treat <- (SStreat/dfTreat)/(SSRes/dfRes)
  # Preparation
  SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl))-mean(circ_ABl$readout))^2))
  F_prep <- (SSprep/dfTreat)/(SSRes/dfRes)
  # Regression
  # ANOVA tape II SS of regression
  SSreg <- Anova(lm(readout ~log_dose + isSample, circ_ABl))[1,1]
  # Non-parallelism
  # diff of RSS of restricted and unrestricted model
  SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2)
  F_nonpar <- SSnonpar/(sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)/(lenCirc-4))
  
  # non-linearity
  SSnonlin <- sum((predict(modABu)-predict(lm(readout ~ as.factor(log_dose)*isSample, circ_ABl)))^2)
  # = RSS-SSE
  # Total SS
  SStot <- sum((circ_ABl$readout-mean(circ_ABl$readout))^2)
  # Significance of R^2 F-ratio
  # MSR/MSE
  # sample A
  F_R2_A <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2)/
    (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - circ_Al$readout)^2)/(nrow(circ_Al)-3))
  pFR2_A <- round(pf(F_R2_A,1,6),4)
  # sample B
  F_R2_B <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2)/
    (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2)/(nrow(circ_Bl)-3))
  pFR2_B <- round(pf(F_R2_B,1,6),4)
  # sign of non-lin with pure error: MSSnonlin/MSSE
  F_nonlin <- (SSnonlin/2)/(SSE/dfPureE)
  
  # sign of slope
  F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2)/(sum((circ_Bl$readout - predict(modB))^2)/(nrow(circ_Bl)-2))
  F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2)/(sum((circ_Al$readout - predict(modA))^2)/(nrow(circ_Al)-2))
  # F-test on regression: MSSreg/MSSE
  if (is.na(F_nonlin)) F_nonlin <- 0
  if (F_nonlin > 0) {
    p_F_nonlin <- round(pf(F_nonlin,2,dfPureE, lower.tail = F),5)
  } else { p_F_nonlin <- "SSnonlin neg or 0"; }
  
  # significances
  F_regr <- (SSreg/1)/(SSRes/dfRes)
  p_F_regr <- round(pf(F_regr,1,dfRes, lower.tail = F),3)
  p_F_treat <- round(pf(F_treat,3,dfRes, lower.tail = F),3)
  p_F_prep <- round(pf(F_prep,1,dfRes, lower.tail = F),3)
  p_F_slope_A <- round(pf(F_slope_A,1,(nrow(circ_Al)-2), lower.tail = F),3)
  p_F_slope_B <- round(pf(F_slope_B,1,(nrow(circ_Bl)-2), lower.tail = F),3)
  p_F_nonp <- round(pf(F_nonpar,1,dfRes, lower.tail = F),3)
  p_F_LoF <- p_F_nonlin
  
  res_tab_lin <- data.frame(test = c("F-test on sign. of regression", "F_test on non-lin",
                                     "F-test on R^2 A","F_test on R^2 B",
                                     "F-test on slope A","F-test on slope B",
                                     "F-test on non-parallelism","F-test on preparation"),
                            test_results = c(ifelse(p_F_regr<0.05,0,1),ifelse(p_F_nonlin<0.05,1,0),
                                             ifelse(pFR2_A<0.05,1,0),ifelse(pFR2_B<0.05,1,0),
                                             ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1),
                                             ifelse(p_F_nonp<0.05,1,0),ifelse(p_F_prep<0.05,0,1)),
                            estimate = c(p_F_regr, p_F_nonlin,pFR2_A,pFR2_B,p_F_slope_A,
                                         p_F_slope_B,p_F_nonp,p_F_prep),
                            Source = c("Treatment","Preparation","Regression","Non-parallelism",
                                       "Resid Error","Non-linearity","Pure error", "Total"),
                            df = c(dfTreat,1,1,1,dfRMSE,2,dfPureE,lenCirc-1),
                            SumSquares = c(round(SStreat,5),round(SSprep,5),round(SSreg,5),
                                           round(SSnonpar,5),round(RSS,5),round(SSnonlin,5),
                                           round(SSE,5),round(SStot,5)),
                            MS = c(round(SStreat/dfTreat,5),round(SSprep,5),round(SSreg,5),
                                   round(SSnonpar,5),round(RSS/dfRMSE,5),round(SSnonlin/2,5),
                                   round(SSE/dfPureE,5),round(SStot/dfTotal,5)),
                            "F-value" = c(round(F_treat,5), round(F_prep,5),round(F_regr,5),
                                          round(F_nonpar,5),"",round(F_nonlin,5),"",""),
                            "p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "",""))
  RET  <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2)
  return(RET)
}

pot4plFUNC <- function(ro_new, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
  isRef <- rep(c(1,0),1,each=nrow(all_l)/2) 
  isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01
  all_l$readouttrans <- log(all_l$readout)
  
  CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)])
  if (CORdat<0) SLOPE <- -1 else SLOPE <- 1
  
  startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0)
  tryCatch({
    mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                  data=all_l,
                  start=startlist,
                  control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose),
                      at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0)
  
  tryCatch({
    mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
                    (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))),
                  data=all_l,
                  start=startlistmu,
                  control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  if (!PureErrFlag) {
    pot_est <- exp(confintd(mr, "r", method="asymptotic"))
    potU_est <- exp(confintd(mu, "r", method="asymptotic"))
  } else {
    FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
    meanPureErr <- FitAnova[4,3]
    SU_mr <- tryCatch({
      SU_mr <- summary(mr)
    }, error = function(msg) {
      return()
    })
    
    #browser() 
    if (length(SU_mr)>1) {
      s_mr <- SU_mr$coefficients
    } else { SU_mr <- rep(NA,5) }
    
    VCOV <- vcov(mr)
    V_V <- VCOV/SU_mr$sigma^2
    SEsPure <- sqrt(diag(V_V)*meanPureErr)
    pot_est <- c(exp(s_mr[5,1]),exp(s_mr[5,1]-qt(0.975,nrow(all_l)-5)*SEsPure[5]),
                 exp(s_mr[5,1]+qt(0.975,nrow(all_l)-5)*SEsPure[5]))
    # unrestricted
    s_mu <- summary(mu)$coefficients
    SU_mu <- summary(mu)
    VCOVu <- vcov(mu)
    V_Vu <- VCOVu/SU_mu$sigma^2
    SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
    potU_est <- c(exp(s_mu[7,1]),exp(s_mu[7,1]-qt(0.975,nrow(all_l)-8)*SEsPureU[7]),
                  exp(s_mu[7,1]+qt(0.975,nrow(all_l)-8)*SEsPureU[7]))
  } # PureErrFlag
  
  startlistmr_log <- list(a=max(all_l$readouttrans), b=SLOPE, d=min(all_l$readouttrans), cs=mean(all_l$log_dose),r=0)
  
  tryCatch({
    mr_log <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                      data=all_l,
                      start=startlistmr_log,
                      control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  
  startlistmu_log <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose),
                          at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0)
  
  tryCatch({
    mu_log <- gsl_nls(fn = readouttrans ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
                        (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))),
                      data=all_l,
                      start=startlistmu_log,
                      control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic"))
  potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic"))
  colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI")
  colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI")
  #browser()   
  su_mr_log <- summary(mr_log)
  Dat$RMSE_Rlog <- su_mr_log$sigma
  su_mu_log <- summary(mu_log)
  Dat$RMSE_Ulog <- su_mu_log$sigma
  Dat$up_lowAslog <- su_mu_log$coefficients[1,1] - su_mu_log$coefficients[4,1]
  potALL <- rbind(pot_est, potU_est, pot_est_log, potU_est_log)
  
  potALL2 <- cbind(c("restricted","unrestricted","transformed restr","untransf restr"), potALL)
  return(potALL2)
}

ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVarlog,DFs, Conf=0.975) {
  log_xs <- log(abs(xs))
  log_xt <- log(abs(xt))
  var_log_xs <- (se_xs/xs)^2 # approximate variance of log(bs)
  var_log_xt <- (se_xt/xt)^2
  se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVarlog)
  
  lower_log_ratio <- log_xt-log_xs - qt(Conf,DFs)*se_log_ratio
  upper_log_ratio <- log_xt-log_xs + qt(Conf,DFs)*se_log_ratio
  ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio))
  return(ci_ratio)
}

tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
  isRef <- rep(c(1,0),1,each=nrow(all_l)/2) 
  isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01
  
  pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
             pmodels=data.frame(1,1,1,isSample))
  compParm(pot, "c",display=T)
  ED50 <- ED(pot,c(50), interval="delta")
  PotEst <- ED50[1,1]/ED50[2,1]
  potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE)
  potAll2 <- potAll[1:3]
  
  CORro <- cor(ro_new[,1], ro_new[,ncol(ro_new)])
  
  if (CORro<0) SLOPE <- -1 else SLOPE <- 1
  startlist <- list(a=max(ro_new[,2]), b=SLOPE, d=min(ro_new[,2]), cs=mean(all_l$log_dose),r=0)
  tryCatch({
    mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                  data=all_l,
                  start=startlist,
                  control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose),
                      at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0)
  tryCatch({
    mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
                    (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))),
                  data=all_l,
                  start=startlistmu,
                  control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  smu <- tryCatch({ summary(mu) },
                  error=function(err) {
                    return(0)
                  })
  
  POTr_CI <- potAll2[2:3]
  
  FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
  # pure error
  pureSS <- FitAnova[4,2]
  pureSS_df <- FitAnova[4,1]
  meanPureErr <- FitAnova[4,3]
  vcovMU <- vcov(mu)
  V_V <- vcovMU/smu$sigma^2
  SEsPure <- sqrt(diag(V_V)*meanPureErr)
  VCOVpure <- V_V*meanPureErr
  DFsPure <- FitAnova[4,1]
  
  
  
  testPOTr <- logical()
  if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100<Lim[[10]] ) testPOTr <- 0 else testPOTr <- 1
  
  potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
              pmodels=data.frame(isSample, isSample,isSample,isSample))
  potAllU <- EDcomp(potU, percVec=c(50,50), interval="delta", display=FALSE)
  potAllU2 <- potAllU[1:3]
  sum_potU <- summary(potU)
  coeffs <- potU$coefficients
  coeffs[1] <- ifelse(CORro<0, -coeffs[1], coeffs[1])
  coeffs[2] <- ifelse(CORro<0, -coeffs[2], coeffs[2])
  names(coeffs) <- c("bs","bt","ds","dt","as","at","cs","ct")
  
  e_c_ref <- coeffs[7]
  e_c_test <- coeffs[8]
  coeffs[7:8] <- log(coeffs[7:8])
  test_c <- logical()
  if((potAllU2[2] >Lim[[9]]/100 & potAllU2[3] <Lim[[10]]/100)) test_c <- 0 else test_c <- 1
  
  #### ANOVA ----
  noConc <- length(unique(all_l$Conc))
  nofitted <- noConc
  AnovaDFs <- c(nofitted-1,1,3,nofitted-4-1,nrow(all_l)-nofitted, nofitted,nrow(all_l)-2*nofitted,nrow(all_l)-1)
  SStreat <- round(sum((predict(potU)-mean(all_l$readout))^2),5)
  SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5)
  # non-parallelism
  SSnonparall <- round(sum(resid(pot)^2)-sum(resid(potU)^2),5)
  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l))-mean(all_l$readout))^2),5)
  
  RSS <- round(sum(potU$predres[,2]^2),5)
  RSS_df <- AnovaDFs[5]
  MSEunr <- RSS/RSS_df
  RMSEunr <- sqrt(RSS/RSS_df)
  # Pure Err
  FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
  SSE <- sum(resid(lm(readout ~ factor(Conc)*isSample, all_l))^2) # =FitAnova[4,2]
  SSE_df <- FitAnova[4,1]
  PureMSE <- SSE/SSE_df
  RMSE_pure <- sqrt(PureMSE)
  ## non-lin = LoF
  if (PureErrFlag) { ERR <- PureMSE; ERR_df <- SSE_df } else { ERR <- MSEunr; ERR_df <- RSS_df }
  SSnonlin <- sum((predict(lm(readout ~ factor(Conc)*isSample, all_l))-predict(potU))^2)
  LoF_df <- FitAnova[1,1]+FitAnova[2,1]
  F_regr <- (SSregr/AnovaDFs[3])/ERR
  p_F_regr <- round(pf(F_regr, AnovaDFs[3], ERR_df, lower.tail = F),5)
  if (ncol(ro_new)<4) F_nonlin <- 0 else F_nonlin <- (SSnonlin/AnovaDFs[6])/ERR
  if (F_nonlin > 0) {
    p_F_nonlin <- round(pf(F_nonlin, AnovaDFs[6], ERR_df, lower.tail = F),5)
  } else { p_F_nonlin <- "SSnonlin neg or single dilutions" }
  
  test_a <- test_b <- test_d <- test_ad <- logical()
  
  RSS_r <- round(sum(pot$predres[,2]^2),5)
  MSE_r <- RSS_r/(nrow(all_l)-5)
  RMSE_r <- round(sqrt(MSE_r),6)
  Dat$RMSE_r <- RMSE_r
  Dat$RMSE_pure <- RMSE_pure
  Dat$RMSE_unr <- round(RMSEunr,6)
  #browser()  
  ## EQ test on lower As diff
  ds <- coeffs["ds"]
  dt <- coeffs["dt"]
  lAs_diff <- (dt-ds)
  uCI_laDiff <- lAs_diff+qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,2]^2)
  lCI_laDiff <- lAs_diff-qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,2]^2)
  if (uCI_laDiff < Lim[[2]] & lCI_laDiff > Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1
  
  #### EQ test on upper asymptote ratio ----
  as <- coeffs["as"]
  at <- coeffs["at"]
  uAsRatio <- compParm(potU, "a","/",display=F)
  uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2])
  #browser()  
  ds <- smu$coefficients["ds",1]
  dt <- smu$coefficients["dt",1]
  if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds","ds"]) else se_ds <- smu$coefficients["ds",2]
  if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt","dt"]) else se_dt <- smu$coefficients["dt",2]
  if (PureErrFlag) CoVarlog_d <- VCOVpure["dt","ds"] else CoVarlog_d <- vcovMU["dt","ds"]
  if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-noConc
  uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.9975)
  if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1
  estUppA <- round(at/as,5)
  
  Dat$uAsCI <- uAsCI2
  
  #### EQ test on slope ratio ----
  bs <- coeffs["bs"]
  bt <- coeffs["bt"]
  slopeRatio <- compParm(potU, "b","/",display=F)
  slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2])
  
  bs <- smu$coefficients["bs",1]
  bt <- smu$coefficients["bt",1]
  if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs","bs"]) else se_bs <- smu$coefficients["bs",2]
  if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt","bt"]) else se_bt <- smu$coefficients["bt",2]
  if (PureErrFlag) CoVarlog_b <- VCOVpure["bt","bs"] else CoVarlog_b <- vcovMU["bt","bs"]
  slopeCI2 <- ParamCI_F(bt,bs,se_bt, se_bs,CoVarlog_b, DFs, Conf=0.975)
  if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1
  estUppA <- round(at/as,5)
  
  Dat$slopeRatioCI <- slopeCI
  
  #### EQ test on lower As ratio ----
  
  lAsRatio <- compParm(potU, "d","/",display=F)
  slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2])
  
  as <- smu$coefficients["as",1]
  at <- smu$coefficients["at",1]
  if (PureErrFlag) se_as <- sqrt(VCOVpure["as","as"]) else se_as <- smu$coefficients["as",2]
  if (PureErrFlag) se_at <- sqrt(VCOVpure["at","at"]) else se_at <- smu$coefficients["at",2]
  if (PureErrFlag) CoVarlog_a <- VCOVpure["at","as"] else CoVarlog_a <- vcovMU["at","as"]
  lAsCI2 <- ParamCI_F(at,as,se_at, se_as,CoVarlog_a, DFs, Conf=0.975)
  if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1
  estLowA <- round(at/as,5)
  
  Dat$lAsCI <- lAsCI2
  
  #### EQtest on ratio of As difference ----
  AsDiffRatio <- (at-dt)/(as-ds)
  
  at_dt <- (at-dt)
  as_ds <- (as-ds)
  se_ds_asPure <- sqrt(VCOVpure["as","as"]+VCOVpure["ds","ds"]-2*VCOVpure["as","ds"])
  se_dt_atPure <- sqrt(VCOVpure["at","at"]+VCOVpure["dt","dt"]-2*VCOVpure["at","dt"])
  se_ds_asRMSE <- sqrt(vcovMU["as","as"]+vcovMU["ds","ds"]-2*vcovMU["as","ds"])
  se_dt_atRMSE <- sqrt(vcovMU["at","at"]+vcovMU["dt","dt"]-2*vcovMU["at","dt"])
  if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE
  if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE
  
  AsDiffCI2 <- ParamCI_F(at_dt,as_ds,se_dt_at, se_ds_as,CoVarlog=0, DFs, Conf=0.975)
  if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1
  estLowA <- round(at/as,5)
  
  Dat$up_lowAs <- (as-ds)
  
  lowerCIlowerA <- lAsCI2[1]; lowerCIupperA <- uAsCI2[1]; upperCIlowerA <- lAsCI2[2]; upperCIupperA <- uAsCI2[2]
  test_lowA <- test_d; test_uppA <- test_a
  #browser() 
  res_tab <- data.frame(test= c("F-test on sign. of regression*",
                                "EQ test on lower asymptotes difference",
                                "EQ test ratio of lower asymptotes",
                                "EQ test ratio of Hill slopes",
                                "EQ test ratio of upper asymptotes",
                                "F-test on non-linearity*",
                                "EQ test ratio of asymptote difference",
                                "geom. rel. CI restr. model",
                                "geom. rel. CI unrestr. model"),
                        test_results = c(ifelse(p_F_regr<0.05,0,1), test_la_diff, test_lowA, test_b, test_uppA,
                                         ifelse(p_F_nonlin>1,1, ifelse(p_F_nonlin<0.05,1,0)), test_ad,
                                         testPOTr, test_c),
                        estimate = c(round(p_F_regr, 3), round(lAs_diff, 5),
                                     estLowA, round(bs/bt,5), estUppA, p_F_nonlin, 
                                     round(at_dt/as_ds, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)),
                        lower_limit = c("-",Lim[[1]],Lim[[3]],Lim[[5]],Lim[[7]],"-",Lim[[11]],Lim[[9]],Lim[[9]]),
                        upper_limit = c("-",Lim[[2]],Lim[[4]],Lim[[6]],Lim[[8]],"-",Lim[[12]],Lim[[10]],Lim[[10]]),
                        lower_CI = c(RMSE_r, round(lCI_laDiff,3), round(lAsCI2[1],5), round(slopeCI2[1],5),
                                     round(uAsCI2[1],5), "-", round(AsDiffCI2[1],5), round(potAll2[2],2), round(potAllU2[2],2)),
                        upper_CI = c(RMSE_pure, round(uCI_laDiff,3), round(lAsCI2[2],5), round(slopeCI2[2],5),
                                     round(uAsCI2[2],5), "-", round(AsDiffCI2[2],5), round(potAll2[3],2), round(potAllU2[3],2))
  )
  return(res_tab)
}

ANOVA4plUnresfunc <- function(ro_new, sigmoid) {
  all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
  all_len <- nrow(all_l)
  isRef <- rep(c(1,0),1,each=all_len/2) 
  isSample <- rep(c(0,1),1,each=all_len/2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01
  
  pot <- drm(readout ~ Conc, isSample, data=all_l2, fct=LL.4(names=c("b","d","a","c")),
             pmodels=data.frame(1,1,1,isSample))
  potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
              pmodels=data.frame(isSample, isSample,isSample,isSample))
  SStreat <- round(sum((potU$predres[,1] - mean(all_l$readout))^2),5)
  SStreat_df <- length(unique(all_l$log_dose))-1
  SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5)
  ## Non-parallel
  SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5)
  ## Preparation
  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5)
  ## Resid Err
  RSS <- round(sum(potU$predres[,2]^2),5)
  RSS_df <- nrow(all_l)-SStreat_df-1
  FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
  # PureErr
  SSE <- FitAnova[4,3]
  SSE_df <- FitAnova[4,1]
  # Non-Linearity
  SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^2),4)
  LoF_df <- FitAnova[1,1]+FitAnova[2,1]
  ## Total
  SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5)
  MSE <- RSS/RSS_df
  noConc <- length(unique(all_l$Conc))
  AnovaDFs <- c(noConc-1, 1,3,noConc-4-1, nrow(all_l)-noConc, noConc, nrow(all_l)-noConc-noConc, nrow(all_l)-1)
  p_SStreat <- round(pf((SStreat/AnovaDFs[1])/MSE, AnovaDFs[1],RSS_df, lower.tail = F),3)
  p_SSprep <- round(pf((SSprep/AnovaDFs[2])/MSE, AnovaDFs[2],RSS_df, lower.tail = F),3)
  p_SSregr <- round(pf((SSregr/AnovaDFs[3])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3)
  p_SSnonp <- round(pf((SSnonparallel/AnovaDFs[4])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3)
  p_SSLoF <- round(pf((SSnonlin/LoF_df)/(SSE/SSE_df), LoF_df,SSE_df, lower.tail = F),5)
  
  ANOVAtab <- data.frame(Source = c("Treatment","Preparation","Regression",
                                    "Non-Parallelism","Residual Error","Non-linearity",
                                    "Pure Error","Total"),
                         DF = AnovaDFs,
                         SumSquares = c(SStreat, SSprep,SSregr, SSnonparallel,
                                        RSS, SSnonlin,SSE, SStot),
                         MeanSquares = c(round(SStreat/AnovaDFs[1],3), SSprep, round(SStreat/AnovaDFs[3],3),round(SSnonparallel/AnovaDFs[4],3),
                                         round(MSE,5), round(SSnonlin/LoF_df,5), round(SSE/SSE_df,5),""),
                         "F-value" = c(round((SStreat/AnovaDFs[1])/MSE,5), round((SSprep/AnovaDFs[2])/MSE,5),
                                       round((SSregr/AnovaDFs[3])/MSE,5),round((SSnonparallel/AnovaDFs[4])/MSE,5),
                                       "",round((SSnonlin/LoF_df)/(SSE/SSE_df),5),"",""),
                         "p_value" = c(round(p_SStreat,3), p_SSprep, round(p_SSregr,3), p_SSnonp,"",p_SSLoF,"","")
  )
  
  return(ANOVAtab)
}

perConcTab <- function(ro_new, noDilSeries) {
  Reftab <- ro_new[,c(1:noDilSeries)]
  Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))]
  tReftab <- t(Reftab)
  colnames(tReftab) <- round(ro_new[,ncol(ro_new)],5)
  
  avs <- apply(tReftab,2,mean)
  sds <- apply(tReftab,2,sd)
  cv <- sds/avs*100
  tReftab2 <- rbind(tReftab, avs,sds,cv)
  
  tTesttab <- t(Testtab)
  colnames(tTesttab) <- round(ro_new[,ncol(ro_new)],5)
  
  avs_test <- apply(tTesttab,2,mean)
  sds_test <- apply(tTesttab,2,sd)
  cv_test <- sds_test/avs_test*100
  tTesttab2 <- rbind(tTesttab, avs_test,sds_test,cv_test)
  concTab <- rbind(tReftab2, tTesttab2)
  return(concTab)
  
}

divFUN <- function(x,Div,N,res,noDil) {
  N <- N+1
  y <- x/Div
  res <- c(res,y)
  if (N==noDil)  { return(res) }
  divFUN(y,Div,N,res,noDil)
}





server = function(input,output,session) {
  ReportParS <- reactiveValues()
  IPReportParS <- reactiveValues()


  output$homePage <- renderUI({
    navbarPage("Home",
               tabPanel("Introduction",
                        tags$img(src="logo.png", class="adv_logo"),
                        h4("Introduction to the bioassay software"),
                        tags$mark("linear regression"), br(),

                        ),
               tabPanel("Documentation",
                        h4("Introduction "),
                        h4("Information on dilution setting"),
                        "(for details see:", a(href="ADONIS.pdf","Whitepaper", download=NA, target="_blank"),")",tags$br(),
                        "Bend points are calculated according to following formula:",
                        withMathJax(" $$bp_{1/2} = \\pm\\frac{1.31696}{Hill's slope}$$")),
               tabPanel("Configuration",
                        verbatimTextOutput("sessioninfo"))
               )
  })

  output$Dataupload <- renderUI({
    navbarPage(title="Information",
               tabPanel(title = "Upload all worksheets",
                        box(title = "Upload", status="warning",solidHeader = T, width=4, "Please upload your EXCEL file here",
                            fileInput("XLfile",'',accept=".xlsx")),
                        "Two possibilities of data structures...",
                        img(src="Screenshot.png", width=400),
                        tags$head(tags$style(HTML("label {font-size:80%;margin-bottom: 3px;margin-top: 3px;}"))),
                        column(2,
                               tags$image(src="logo.png", class="adv_logo"),
                               h4("upload EXCEL"),
                               fileInput(inputId = "iFile", label = "", accept="ms-excel"),
                               uiOutput(outputId = "sheetName"),
                               div(checkboxInput("PureErr", "Should pure error be used for calculation of CIs?", FALSE),
                                   style = "font-size: 24px !important;color: red"),
                               verbatimTextOutput("stats"),
                               h5("\n\n\n Author: Franz Innerbichler, InnerAnalytics")),
                        div(id="parameter",
                            column(1,style = "background: lightgrey",
                                   #actionButton("StartCalc", "Click, when calculations to be started"),
                                   h4("curve settings"),
                                   numericInput("lowAsymptREF", "lower asymptote REF",10,step=1,min=0),
                                   numericInput("lowAsymptTEST", "lower asymptote TEST",10,step=1,min=0),
                                   numericInput("uppAsymptREF", "upper asymptote REF",110,step=1,min=0),
                                   numericInput("uppAsymptTEST", "upper asymptote TEST",110,step=1,min=0)
                            ),
                            column(1,style = "background: lightgrey",
                                   numericInput("slopeREF", "slope REF",1,step=0.1,min=-10),
                                   numericInput("slopeTEST", "slope TEST",1,step=0.1,min=-10),
                                   numericInput("EC50", "EC50 REF",-3.5),
                                   numericInput("potDiff", "potency difference",0)
                            ),
                            column(1,style = "background: lightgreen",
                                   h4("dilutions"),
                                   numericInput("CONC1", "highest concentration",0.3, min=-3.5),
                                   numericInput("CONC2", "2nd concentration",0.15),
                                   numericInput("CONC3", "3rd concentration",0.075),
                                   numericInput("CONC4", "4th concentration",0.0375),
                                   numericInput("CONC5", "5th concentration",0.01875),

                                   numericInput("CONC6", "6th concentration",0.00938)
                            ),
                            column(1,style = "background: lightgreen",

                                   numericInput("CONC7", "7th concentration",0.00469),
                                   numericInput("CONC8", "8thd concentration",0.00235),
                                   numericInput("CONC9", "9thd concentration",NA),
                                   numericInput("CONC10", "10th concentration",NA),
                                   numericInput("CONC11", "11th concentration",NA),

                                   numericInput("CONC12", "lowest concentration",NA)
                            ),
                            column(1,style = "background: lightblue",
                                   h4("geometric dilution scheme"),
                                   numericInput("ConcStart", "starting concentration",NA),
                                   numericInput("dilutionFac", "dilution factor",NA),
                                   numericInput("NoDil", "no. of dilutions",NA),
                                   numericInput("NoDilSer", "no. of dil. series",NA),
                                   verbatimTextOutput("dilutions")
                            ),
                            column(4,
                                   h4("log-dilution scheme"),
                                   verbatimTextOutput("logdil"),
                                   h4("User help"),
                                   h5("If new dilutions are entered, please adjust EC50 to avoid calculation errors"))
                        )

                        ))
  })

  output$fourPL <- renderUI({
    navbarPage(title="4PL",
               tabPanel("Analysis and Plots",
                        sidebarLayout(
                          sidebarPanel(
                            width=3,
                          fluidRow(
                            column(6,
                                   numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")),
                                                bsPopover(id="q4", title="", content="The calculated limits ...")))
                          )),
                          mainPanel(
                            tabsetPanel(id="tabs",
                                        tabPanel("4pl",
                                                 box(title="ANOVA table", status="primary",solidHeader = T, width=12,
                                                     tableOutput("Anovatab")),
                                                 column(2,
                                                        h5("Diagnostics only shown if EXCEL is uploaded"),
                                                        htmlOutput("PureErrW2"),
                                                        tags$head(tags$style("#PureErrW2{color: red;
                                                             font-size: 16px;
                                                             font_style: italic;}")),
                                                        tableOutput("FileSAmpl"),
                                                        downloadButton("downloadXLReport", label="Download PDF report", class="butt"),
                                                        tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"),
                                                        plotOutput("relpotTestPlot", width="300px", height="150px"),
                                                        h4("Akaike Information Criterion"),
                                                        tableOutput("AIC"),
                                                        h5("First row: restricted model; 2nd row: unrestricted model"),
                                                        h5("Smaller values of AIC indicate better fit to the data"),
                                                        tableOutput("VarDiagn")
                                                 ),
                                                 column(4,
                                                        plotOutput("XLplot"),
                                                        column(6,
                                                               br(),
                                                               "Regression results restricted",
                                                               tableOutput("coeffs_r"),
                                                               "Bend points restricted",
                                                               tableOutput("coeffs_r2")),
                                                        column(6,
                                                               br(),
                                                               "Regression results unrestricted",
                                                               tableOutput("coeffs_unr"))),
                                                 column(4,
                                                        plotOutput("diagnplot"),
                                                        column(6,
                                                               tableOutput("logcoeffs_r"),
                                                               tableOutput("coeffs_unr2")),
                                                        column(6,
                                                               tableOutput("logcoeffs_unr"))),
                                                 column(2,
                                                        tableOutput("ANOVAXLS"),
                                                        DT::renderDataTable("EQtests"))
                                        ),
                                        tabPanel("Report",
                                                 h4("Settings for report"),
                                                 ))
                          )
                        )))
  })

  output$pla <- renderUI({
    navbarPage(title="pla",
               tabPanel("Analysis and Plots",
                        sidebarLayout(
                          sidebarPanel(
                            width=3,
                            fluidRow(
                              column(6,
                                     numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")),
                                                  bsPopover(id="q4", title="", content="The calculated limits ...")))
                            )),
                          mainPanel(
                            tabsetPanel(id="tabs",
                                        tabPanel("4pl",
                                                 box(title="ANOVA table", status="primary",solidHeader = T, width=12,
                                                     tableOutput("Anovatab")),
                                                 column(6,
                                                        htmlOutput("PureErrW3"),
                                                        tags$head(tags$style("#PureErrW2{color: red;
                                                             font-size: 16px;
                                                             font_style: italic;}")),
                                                        plotOutput("plotLin"),
                                                        "Delta method is used for potency CIs",
                                                        DT::dataTableOutput("pottab"),
                                                        h4("Unrestricted linear model (SSSI):"),
                                                        tableOutput("SummaryModABu"),
                                                        h4("Restricted linear model (CSSI):"),
                                                        tableOutput("SummaryModAB")),
                                                 column(3,
                                                        h3("Tests for linear PLA):"),
                                                        DT::dataTableOutput("TESTSlin"),
                                                        h5("The estimate is the p-value of the test"),
                                                        h5("F-tests on regression, significance of slopes, and preparation need to have a p-value <0.05 to pass"),
                                                        h5("All other tests pass if p-value > 0.05"),
                                                        "SST CI for difference of slopes:",
                                                        tableOutput("SlopeDiffCI")),
                                                 column(3,
                                                        h3("ANOVA for parallel line assay"),
                                                        DT::dataTableOutput("ANOVAlin"))),
                                        tabPanel("Report",
                                                 h4("Settings for report")
                                        ))
                          )
                )))
  })

  output$wizard <- renderUI({
    navbarPage(title="Dilution setting",
               tabPanel("Plots",
                        sidebarLayout(
                          sidebarPanel(
                            width=3,
                            fluidRow(
                              column(6,
                                     numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")),
                                                  bsPopover(id="q4", title="", content="The calculated limits ...")))
                            )),
                          mainPanel(
                            tabsetPanel(id="tabs",
                                        tabPanel("4pl",
                                                 box(title="ANOVA table", status="primary",solidHeader = T, width=12,
                                                     tableOutput("Anovatab")),
                                                 column(4,
                                                        h3("Confidence intervals"),
                                                        tableOutput("CIs"),
                                                        "The confidence intrval table is interaactive for changes in: variability slider (%SD), potency of test-slider,
                                       and 'Adjust the dilutions'-slider",
                                                        tableOutput("optimalDils"),
                                                        selectInput(inputId="scenario", label= "Select an 'optimal' scenario:", choices = c("scenario 2","scenario 3","scenario 6","steep slope"))),
                                                 column(5,
                                                        plotOutput("plotfordilutions"),
                                                        h4("in grey: most extreme bend point lines of theoretical samples with 50% and 200% potency"),
                                                        sliderInput("dilslider", "Adjust the dilutions(+-change in %)", min = -100,max=100, value=0, step=1, round=0),
                                                        checkboxInput("fixupper","Fix highest concentration (if unticked, the center is fixed)",FALSE),
                                                        h5("Dilution factors"),
                                                        tableOutput("adjlogdil"),
                                                        "Short guidance: wider dilution ranges increase the CIs of rel. potency, and decrease the CIs of upper and lower asymptote ratios, as well as Hill's slope ratios", br(),
                                                        "Narrower dilution ranges decrease the CIs of rel. potency, and increase the CIs of upper and lower asymptote ratios, ands Hill's slope ratios",),
                                                 column(3,
                                                        h3("Bend points"),
                                                        tableOutput("bps"),
                                                        tableOutput("extremebps"),
                                                        h4("Explanation of the plots")
                                                 )),
                                        tabPanel("Report",
                                                 h4("Settings for report")
                                        ))
                          )
                        )))
  })

  
}
  

