lexibel wenn fehlende Werte

R/Global.R

@@ -23,6 +23,26 @@ library(car)
 library(dplyr)
 library(scales)
 
+#' Estimate correlations
+#'
+#' returns the correlation of 2 vectors
+#'
+#' @param vec1 The 1st vector.
+#' @param vec2 The 2nd vector. 
+#' @returns A float as correlatioin estimate
+#' @export
+#' @examples
+#' suppressMessages(source("../../dev/setup.R"))
+#' vector1 <- c(1,2,3,4,5)
+#' vector2 <- c(5.1,4.3,NA,1.9,1.2)
+#' te <- COR_FUNC(vector1,vector2)
+#' print(te)
+COR_FUNC <- function(vec1, vec2) {
+  df <- data.frame(v1 = vec1, v2 = vec2)
+  df2 <- df[complete.cases(df),]
+  #browser()
+  return(cor(df2[,1],df2[,2]))
+}
 
 #' Levenberg Marquard fit of 4 pl
 #'
@@ -46,8 +66,9 @@ library(scales)
 #' Dat <- list()
 #' te <- Fitting_FUNC(dat, TransF)
 #' print(te)
-Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
+Fitting_FUNC <- function(ro_new, TransFlag = FALSE, nameWS ="") {
-  CORro <- cor(ro_new[, 1], ro_new[, ncol(ro_new)])
+  #browser()
+  CORro <- COR_FUNC(ro_new[, 1], ro_new[, ncol(ro_new)])
   # browser()
   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)
@@ -69,6 +90,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
       },
       warning = function(e) {
         mr <<- "In nlsModel singular gradient matrix"
+        
       }
     )
     # Stop if singular gradient matrix
@@ -82,6 +104,12 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
       },
       error = function(err) {
         s_mr <- NULL
+        showModal(modalDialog(
+          title = " fit",
+          paste("fit not possible: EC50 outside concentration range for dataset", nameWS),
+          easyClose = TRUE,
+          footer = NULL
+        ))
       }
     )
   } else {
@@ -98,7 +126,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
   if (!TransFlag) {
     startlistmu <- list(
       as = min(ro_new[, 2]), bs = SLOPE, ds = max(ro_new[, 2]), cs = mean(all_l$log_dose),
-      at = min(ro_new[, 2]), bt = SLOPE, dt = max(ro_new[, 2]), r = 0
+      at = min(ro_new[, 4]), bt = SLOPE, dt = max(ro_new[, 4]), r = 0
     )
     tryCatch(
       {
@@ -120,13 +148,19 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
         summary(mu)
       },
       error = function(msg) {
+        showModal(modalDialog(
+          title = "4PL fit",
+          paste("fit not possible: EC50 outside concentration range for dataset", nameWS),
+          easyClose = TRUE,
+          footer = NULL
+        ))
         return(0)
       }
     )
   } else {
     startlistmu <- list(
       as = log(min(ro_new[, 2])), bs = SLOPE, ds = log(max(ro_new[, 2])), cs = mean(all_l$log_dose),
-      at = log(min(ro_new[, 2])), bt = SLOPE, dt = log(max(ro_new[, 2])), r = 0
+      at = log(min(ro_new[, 4])), bt = SLOPE, dt = log(max(ro_new[, 4])), r = 0
     )
     tryCatch(
       {
@@ -152,16 +186,27 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
       }
     )
   }
+  
+#browser()
   if (!TransFlag) {
-    pot_est <- exp(confintd(mr, "r", method = "asymptotic"))
+#browser()
-    potU_est <- exp(confintd(mu, "r", method = "asymptotic"))
+    if (length(s_mr) ==1 | length(Sum_u) ==1) { 
-    PRED <- predict(mr)
+      return("failed")
-    PREDu <- predict(mu)
+      } else {
+        pot_est <- exp(confintd(mr, "r", method = "asymptotic"))
+        potU_est <- exp(confintd(mu, "r", method = "asymptotic"))
+        PRED <- predict(mr)
+        PREDu <- predict(mu)
+      }
   } else {
-    pot_est <- exp(confintd(mrT, "r", method = "asymptotic"))
+    if (length(s_mr) ==1 | length(Sum_u) ==1) { 
-    potU_est <- exp(confintd(muT, "r", method = "asymptotic"))
+      return("failed")
-    PRED <- predict(mrT)
+      }else {
-    PREDu <- predict(muT)
+        pot_est <- exp(confintd(mrT, "r", method = "asymptotic"))
+        potU_est <- exp(confintd(muT, "r", method = "asymptotic"))
+        PRED <- predict(mrT)
+        PREDu <- predict(muT)
+      }
   }
   return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu))
 }
@@ -212,7 +257,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
 #' p <- plotSingularity(dat)
 #' print(p)
 plotSingularity <- function(dat) { # sigmoid,det_sig,
-  CORdat <- cor(dat[, 1], dat[, ncol(dat)])
+  CORdat <- COR_FUNC(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)
@@ -272,7 +317,7 @@ plotSingularity <- function(dat) { # sigmoid,det_sig,
 #' p <- plot_f(dat, TransFlag)
 #' print(p)
 plot_f <- function(dat, TransFlag = FALSE) { # sigmoid,det_sig,
-  CORdat <- cor(dat[, 1], dat[, ncol(dat)])
+  CORdat <- COR_FUNC(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)
@@ -666,9 +711,10 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
   all_l$isRef <- isRef
   all_l$isSample <- isSample
   all_l$Conc <- exp(all_l$log_dose)
+  all_l <- all_l[complete.cases(all_l),]
   all_lA <- all_l[all_l$isSample == 1, ] # TEST
   all_lB <- all_l[all_l$isSample == 0, ] # REF
-  # browser()
+#browser()
   circ_ABl <- circles
   circ_Al <- circ_ABl[circ_ABl$isSample == 1, ]
   circ_Bl <- circ_ABl[circ_ABl$isSample == 0, ]
@@ -753,40 +799,40 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
   }
 
   # treatment
-  SStreat <- print(sum((predict(lm(readout ~ factor(log_dose) * isSample, circ_ABl)) - mean(circ_ABl$readout))^2))
+  SStreat <- print(sum((predict(lm(readout ~ factor(log_dose) * isSample, circ_ABl)) - mean(circ_ABl$readout, na.rm = T))^2, na.rm = T))
   F_treat <- (SStreat / dfTreat) / (SSRes / dfRes)
   # Preparation
-  SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl)) - mean(circ_ABl$readout))^2))
+  SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl)) - mean(circ_ABl$readout, na.rm = T))^2, na.rm = T))
   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)
+  SSnonpar <- sum(resid(modAB)^2, na.rm = T) - sum(resid(modABu)^2, na.rm = T)
-  F_nonpar <- SSnonpar / (sum(resid(lm(readout ~ factor(log_dose) * isSample, circ_ABl))^2) / (lenCirc - 4))
+  F_nonpar <- SSnonpar / (sum(resid(lm(readout ~ factor(log_dose) * isSample, circ_ABl))^2, na.rm = T) / (lenCirc - 4))
 
   # non-linearity
-  SSnonlin <- sum((predict(modABu) - predict(lm(readout ~ as.factor(log_dose) * isSample, circ_ABl)))^2)
+  SSnonlin <- sum((predict(modABu) - predict(lm(readout ~ as.factor(log_dose) * isSample, circ_ABl)))^2, na.rm = T)
   # = RSS-SSE
   # Total SS
-  SStot <- sum((circ_ABl$readout - mean(circ_ABl$readout))^2)
+  SStot <- sum((circ_ABl$readout - mean(circ_ABl$readout, na.rm = T))^2, na.rm=T)
   # 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) /
+  F_R2_A <- sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Al)) - mean(predict(modA), na.rm = T))^2 - (predict(modA) - mean(circ_Al$readout, na.rm = T))^2, na.rm = T) /
-    (sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Al)) - circ_Al$readout)^2) / (nrow(circ_Al) - 3))
+    (sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Al)) - circ_Al$readout)^2, na.rm = T) / (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) /
+  F_R2_B <- sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Bl)) - mean(predict(modB), na.rm = T))^2 - (predict(modB) - mean(circ_Bl$readout))^2, na.rm = T) /
-    (sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2) / (nrow(circ_Bl) - 3))
+    (sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2, na.rm = T) / (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_B <- sum((predict(modB) - mean(circ_Bl$readout, na.rm = T))^2) / (sum((circ_Bl$readout - predict(modB))^2, na.rm = T) / (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_slope_A <- sum((predict(modA) - mean(circ_Al$readout, na.rm = T))^2) / (sum((circ_Al$readout - predict(modA))^2, na.rm = T) / (nrow(circ_Al) - 2))
   # F-test on regression: MSSreg/MSSE
   if (is.na(F_nonlin)) F_nonlin <- 0
   if (F_nonlin > 0) {
@@ -899,11 +945,12 @@ PlotLinPLA_FUNC <- function(circle, sigmoid, all_l2, pl_df, indS, indT) {
   }
 
   
+
   p <- ggplot(all_l2, aes(x = log_dose, y = readout, color = factor(isRef))) +
     geom_point(size = 2) +
     # labs(title=paste("linear regression model", indS,indT), color="product") +
     scale_colour_manual(labels = c("test", "reference"), values = c("#C2173F", "#4545BA")) +
-    ylim(min(all_l2$readout), max(all_l2$readout)) +
+    
     scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
     scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
     theme_bw()
@@ -937,6 +984,7 @@ PlotLinPLA_FUNC <- function(circle, sigmoid, all_l2, pl_df, indS, indT) {
     x = log_dose, y = readout, shape = factor(isRef),
     size = 5, alpha = 0.2
   ), col = c("black"), inherit.aes = FALSE) +
+    ylim(min(all_l2$readout), max(all_l2$readout)) +
     scale_shape_manual(labels = c("test", "reference"), values = c(21, 21))
   # fit intercept for test and ref and common slope
 
@@ -973,6 +1021,7 @@ PlotLinPLA_FUNC <- function(circle, sigmoid, all_l2, pl_df, indS, indT) {
       title = paste("restricted linear regression model"),
       subtitle = paste("Regression on highlighted points")
     ) +
+    ylim(min(all_l2$readout), max(all_l2$readout)) +
     theme(legend.position = "none", axis.text = element_text(size = 14))
   pr3 <- pr2 + geom_point(circle, mapping = aes(
     x = log_dose, y = readout, shape = factor(isRef),
@@ -1018,7 +1067,7 @@ pot4plFUNC <- function(ro_new, PureErrFlag) {
   all_l$readout[all_l$readout < 0] <- 0.01
   all_l$readouttrans <- log(all_l$readout)
   # browser()
-  CORdat <- cor(ro_new[, 1], ro_new[, ncol(ro_new)])
+  CORdat <- COR_FUNC(ro_new[, 1], ro_new[, ncol(ro_new)])
   if (CORdat < 0) SLOPE <- -1 else SLOPE <- 1
   #
   FITs <- Fitting_FUNC(ro_new, TransFlag = FALSE)
@@ -1112,8 +1161,8 @@ ParamCI_F <- function(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975) {
   var_log_xt <- (se_xt / xt)^2
   se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVar/(xs*xt)
 
-  lower_log_ratio <- log_xt - log_xs - qt(Conf, DFs) * se_log_ratio
+  lower_log_ratio <- log_xs - log_xt - qt(Conf, DFs) * se_log_ratio
-  upper_log_ratio <- log_xt - log_xs + qt(Conf, DFs) * se_log_ratio
+  upper_log_ratio <- log_xs - log_xt + qt(Conf, DFs) * se_log_ratio
   ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio))
   return(ci_ratio)
 }
@@ -1143,6 +1192,9 @@ ParamCI_F <- function(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975) {
 #'
 #' tests_FUNC(ro_new=dat, Lim, PureErrF)
 tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
+  
+  DatL <- list()
+  
   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)
@@ -1150,7 +1202,8 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   all_l$isSample <- isSample
   all_l$Conc <- exp(all_l$log_dose)
   all_l$readout[all_l$readout < 0] <- 0.01
-  # browser()
+  all_l <- all_l[complete.cases(all_l),]
+#browser()
   FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE)
   if (is.character(FITs)) {
     return(FITs)
@@ -1172,7 +1225,7 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   VCOVpure <- V_V * meanPureErr
   DFsPure <- FitAnova[4, 1]
 
-
+#browser()
   testPOTr <- logical()
   if (POTr_CI[1] * 100 > Lim[[9]] & POTr_CI[2] * 100 < Lim[[10]]) testPOTr <- 0 else testPOTr <- 1
 
@@ -1188,19 +1241,19 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   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((predPotU - mean(all_l$readout))^2), 5)
+  SStreat <- round(sum((predPotU - mean(all_l$readout, na.rm = T))^2, na.rm = T), 5)
-  SSregr <- round(sum((predPot - mean(all_l$readout))^2), 5)
+  SSregr <- round(sum((predPot - mean(all_l$readout, na.rm=T))^2, na.rm=T), 5)
   # non-parallelism
-  SSnonparall <- round(sum(smr$residuals^2) - sum(smu$residuals^2), 5)
+  SSnonparall <- round(sum(smr$residuals^2, na.rm=T) - sum(smu$residuals^2, na.rm=T), 5)
-  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2), 5)
+  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout, na.rm=T))^2, na.rm=T), 5)
-
+# browser()
-  RSS <- round(sum(smu$residuals^2), 5)
+  RSS <- round(sum(smu$residuals^2, na.rm=T), 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 <- sum(resid(lm(readout ~ factor(Conc) * isSample, all_l))^2, na.rm=T) # =FitAnova[4,2]
   SSE_df <- FitAnova[4, 1]
   PureMSE <- SSE / SSE_df
   RMSE_pure <- sqrt(PureMSE)
@@ -1225,12 +1278,12 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
 
   test_a <- test_b <- test_d <- test_ad <- logical()
 
-  RSS_r <- round(sum(smr$residuals^2), 5)
+  RSS_r <- round(sum(smr$residuals^2, na.rm=T), 5)
   MSE_r <- RSS_r / (nrow(all_l) - 5)
   RMSE_r <- round(sqrt(MSE_r), 6)
-  Dat$RMSE_r <- RMSE_r
+  DatL$RMSE_r <- RMSE_r
-  Dat$RMSE_pure <- RMSE_pure
+  DatL$RMSE_pure <- RMSE_pure
-  Dat$RMSE_unr <- round(RMSEunr, 6)
+  DatL$RMSE_unr <- round(RMSEunr, 6)
 
   coeffs <- smu$coefficients[, 1]
   # browser()
@@ -1242,6 +1295,7 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   lCI_laDiff <- lAs_diff - qt(0.975, smu$df[2]) * sqrt(smu$coefficients["ds", 2]^2 + smu$coefficients["dt", 2]^2)
   if (uCI_laDiff < Lim[[2]] & lCI_laDiff > Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1
 
+  #browser()
   #### EQ test on upper asymptote ratio ----
   # as <- coeffs["as"]
   # at <- coeffs["at"]
@@ -1254,11 +1308,12 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   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) - 8
-  uAsCI2 <- ParamCI_F(dt, ds, se_dt, se_ds, CoVarlog_d, DFs, Conf = 0.975)
+  uAsCI2 <- ParamCI_F(ds, dt, se_dt, se_ds, CoVarlog_d, DFs, Conf = 0.975)
   if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1
-  estUppA <- round(at / as, 5)
+  estUppA <- round(dt / ds, 5)
 
-  Dat$uAsCI <- uAsCI2
+  DatL$uAsCI <- uAsCI2
+# browser()
   
   #### EQ test on slope ratio ----
   # bs <- coeffs["bs"]
@@ -1271,11 +1326,11 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   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)
+  slopeCI2 <- ParamCI_F(bs, bt, 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)
+  estSlope <- round(abs(bt) / abs(bs), 5)
 
-  Dat$slopeRatioCI <- slopeCI2
+  DatL$slopeRatioCI <- slopeCI2
 
   #### EQ test on lower As ratio ----
 
@@ -1287,11 +1342,11 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   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)
+  lAsCI2 <- ParamCI_F(as, at, 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
+  DatL$lAsCI <- lAsCI2
 
   #### EQtest on ratio of As difference ----
   AsDiffRatio <- (dt - at) / (ds - as)
@@ -1305,11 +1360,11 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
   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(dt_at, ds_as, se_dt_at, se_ds_as, CoVar = 0, DFs, Conf = 0.975)
+  AsDiffCI2 <- ParamCI_F( ds_as,dt_at, se_dt_at, se_ds_as, CoVar = 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)
+  estDiffA <- round(dt_at /ds_as, 5)
 
-  Dat$up_lowAs <- abs(ds - as)
+  Dat$estDiffA <- estDiffA
 
   lowerCIlowerA <- lAsCI2[1]
   lowerCIupperA <- uAsCI2[1]
@@ -1337,8 +1392,8 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
     ),
     estimate = c(
       round(p_F_regr, 3), round(lAs_diff, 5),
-      estLowA, round(bs / bt, 5), estUppA, p_F_nonlin,
+      estLowA, estSlope, estUppA, p_F_nonlin,
-      round(dt_at / ds_as, 5), round(potAll2[1] * 100, 2), round(potAllU2[1] * 100, 2)
+      estDiffA, 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]]),
@@ -1375,7 +1430,10 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
 #' ANOVA4plUnresfunc(ro_new)
 #' 
 ANOVA4plUnresfunc <- function(ro_new) {
+  
   all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
+  
+  #browser()
   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)
@@ -1383,6 +1441,8 @@ ANOVA4plUnresfunc <- function(ro_new) {
   all_l$isSample <- isSample
   all_l$Conc <- exp(all_l$log_dose)
   all_l$readout[all_l$readout < 0] <- 0.01
+  all_l <- all_l[complete.cases(all_l),]
+  
   
   FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE)
   smr <- FITs[[1]]

dev/BioassayReportROUT.Rmd

@@ -0,0 +1,335 @@
+---
+output:
+  pdf_document:
+    extra_dependencies: ["float"]
+    number_sections: true
+    toc: true
+    toc_depth: 3
+header_includes:
+  -\usepackage{fancyheadr}
+  -\setlength{\headheight}{22pt}%
+  -\usepackage{lastpage}
+  -\pagestyle{fancy}
+  -\usepackage{pdflscape}
+  -\usepackage{longtable}
+  -\rhead{\includegraphics[width=.15\textwidth]{`r getwd()`/logov2.png}}
+params:
+  FileName: NA
+  author: NA
+  NoP: NA
+  Assay: NA
+  REP: NA
+  coeffs: NA
+author: "Author: `r params$author`"
+title: |
+    | ![](logov2.png){width=1in}
+    | 4PL bioassay evaluation
+subtitle: |
+  `r params$FileName`
+  
+  <left> Unique time: </left> <right> `r Sys.time()`</right>
+date: "`r paste(params$NoP, params$Assay)`"
+
+---
+
+<!-- \fancyfoot[C]{\thepage\ of \pageref{LastPage}} -->
+<!-- \newpage -->
+
+<!-- \newpage -->
+
+```{r setup, include=FALSE}
+
+knitr::opts_chunk$set(echo = TRUE)
+
+library(knitr)
+library(DT)
+library(kableExtra)
+
+REP <- params$REP
+author <- params$author
+coeffs <- params$coeffs
+
+all_l <- REP$all_l
+#ANOVAXLS <- REP$ANOVAXLS
+#XLplot4pl <- REP$XLplot4pl
+DiagnTable <- REP$DiagnTable
+UnRPLAausw <- REP$UnRPLAausw
+UnRPLBend <- REP$UnRPLBend
+PLAausw <- REP$PLAausw
+PLbend <- REP$PLBend
+pottab4plXL <- REP$pottab4plXL
+Lim <- REP$Lim
+XLdat2 <- REP$XLdat2
+PureErr <- REP$PureErr
+ro_newROUT <- REP$ro_newROUT
+ROUTplot <- REP$ROUTplot
+ANOVA_ROUT <- REP$ANOVA_ROUT
+
+CIplot <- REP$CIplot
+testsTabROUT <- REP$testsTabROUT
+relpotTestPlot <- REP$relpotTestPlot
+
+#browser()
+
+```
+
+
+# Introduction
+
+Bioassay potency estimation uses statistical methods to quantify the strength of a biological product or drug by comparing its response to that of a reference standard. Because biological responses are inherently variable, affected by assay conditions, cell systems or organisms, and measurement noise, the 4-parametric logistic regression is used to obtain reliable potency values.  
+USP<1034> recommends calculation of standard errors of ratios of the parameters using Fieller's theorem [1] or using the "delta" method (for a discussion about the "delta" method see [3]). The gradient approach using the differences on the log-scale is mathematically more stable und thus preferable compared to a ratio approach [2].
+
+# Raw data
+
+All data used for the 4PL evaluation is shown in table 1:
+
+```{r alll, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(XLdat2, format = "markdown", caption= "Uploaded data (test and reference) ", digits=3)
+
+if (!is.null(ro_newROUT)) {
+  kable(ro_newROUT, format = "markdown", caption= "Data after exclusion of suspected outliers (see [6]) ", digits=3)
+}
+
+```
+
+
+# Results
+
+## Overall result
+
+```{r Over_all, echo=FALSE, comment=NA, warning=NA, message=NA}
+
+browser()
+potFlag <- 0
+if (pottab4plXL["test_result"][[1]][1]=="failed") potFlag <- 1
+AnalysisFlag <- FALSE
+if (potFlag==1 | sum(testsTabROUT$test_results)>0) AnalysisFlag <- TRUE
+
+colFmt <- function() {
+  
+  outputFormat <- knitr::opts_knit$get("rmarkdown.pandoc.to")
+ if(AnalysisFlag) {
+    text <- paste("\\textcolor{red}{Analysis failed}",sep="")
+  } else  {
+    text <- paste("\\textcolor{black}{Analysis succeeded}",sep="")
+  }  
+  return(text)
+}
+
+
+```
+
+`r colFmt()`
+
+
+## 4pl-regression
+
+Relative potency (absolute and relative confidence limits) are shown in Table 3. `r if(PureErr) {"Pure Error is used for calculations."}` 
+`r if (!PureErr) {"RMSE of restricted model is used for confidence limit calculation."}`
+
+```{r Pot_tab4pl, echo=FALSE, comment=NA, warning=NA, message=NA}
+
+#browser()
+if (pottab4plXL["test_result"][[1]][1]==1) { cat(paste("FAILED: relative potency CL result of restricted model outside limits: ", Lim[[9]], "to" ,Lim[[10]] ))}
+if (pottab4plXL["test_result"][[1]][1]==0) { cat(paste("PASSED: relative potency CL result of restricted model within limits: ", Lim[[9]], "to" ,Lim[[10]] ))}
+kable(pottab4plXL, format = "markdown", caption= "Relative potency with absolute and relative CLs ", digits=3, row.names = F) %>%
+  kable_styling(latex_options = "hold_position")
+
+
+```
+NOTE: results of unrestricted model for Information only.
+
+
+## Plot of the data and models
+
+Plots in Figure 1 shows the restricted model.
+
+
+```{r XLplot, echo=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap="Plot of models", fig.align='left', comment=F, message=F, results='asis', fig.pos='H'}
+
+plot(ROUTplot)
+
+
+```
+
+
+## ANOVA table
+
+The ANOVA of the unconstrained model is listed in table 4. Bates and Watts [4] proposed a test on parallelism which compares the residual sum of squares of the restricted model (ResRSSE) with the residual sum of squares of the unrestricted model (UnresRSSE). If the UnresRSSE is significantly smaller than the ResRSSE, the p-value of "Non-parallelism" is smaller than 0.05 (line 4 in table 4). This test is for information only as it may be overly sensitive in case of small overall variability of the data. 
+
+```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(ANOVA_ROUT, format = "markdown", caption= "Analysis of variance", digits=3) %>%
+  kable_styling(latex_options = "hold_position")
+
+
+```
+
+## Assay suitability tests
+
+Table 5 lists the chosen suitability test results with confidence limits, where applicable. F-tests should be read with caution, if the overall variability is small, as the test gets overly sensitive.
+
+
+```{r SST_ergebn, echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE}
+
+kable(testsTabROUT, row.names = F, format = "markdown", caption="Assay suitability results", digits=4)
+
+
+```
+
+
+\footnotesize
+
+```{r Fussnote, echo=F, comment=NA}
+
+cat("*...The estimate for F-test on regression and on non-linearity is the p-value")
+cat(  "F-test on regression passes if F-value > F-crit and thus p < 0.05")
+cat(  "F-test on non-linearity passes if F-value < F-crit and thus p > 0.05")
+cat(  "Test results outcome:")
+cat("  0 ... test passed (for EQ tests: CL within limits);")
+cat("  1 ... test failed (for EQ tests: CL not within limits);")
+
+
+  
+```
+
+\normalsize
+
+
+```{r AST_Ergebn, echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE}
+
+TestsTabFlag <- FALSE
+if (sum(testsTabROUT$test_results)>0) TestsTabFlag <- TRUE
+colFmt2 <- function() {
+  
+  outputFormat <- knitr::opts_knit$get("rmarkdown.pandoc.to")
+ if(TestsTabFlag) {
+    text <- paste("\\textcolor{red}{Assay suitability tests failed}",sep="")
+  } else  {
+    text <- paste("\\textcolor{black}{Assay suitability tests succeeded}",sep="")
+  }  
+  return(text)
+}
+
+
+```
+
+
+`r colFmt2()`
+
+
+
+
+## Fitting results with curve points
+
+The results of the non-linear fitting procedure for the restricted model (5 parameters) is listed in table 5:
+
+```{r PLAausw, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(PLAausw, format = "markdown", caption= "Restricted 4PL model", digits=3, row.names = F)
+
+
+```
+
+
+
+Sebaugh et al proposed bend points for test and reference samples, that define the points with highest turning behavior. Table 6 lists these bendpoints as well as asymptote points ~ twice as far from the center as the bendpoints.
+
+```{r PLBend, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(PLbend, format = "markdown", caption= "Bendpoints and asymptote points of restricted 4PL model", digits=3)
+
+
+```
+
+The results of the non-linear fitting procedure for the unrestricted model (8 parameters) is listed in table 7:
+
+```{r UnRPLAausw, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(UnRPLAausw, format = "markdown", caption= "Unrestricted 4PL model", digits=3, row.names = F)
+
+```
+
+
+
+# Signature
+
+
+\vspace{1.5cm} 
+\noindent
+\begin{tabular}{p{6cm}p{1cm}p{6cm}}
+  \cline{1-1} \cline{3-3}
+  Date & & Signature 
+\end{tabular}
+
+
+\newpage
+
+# Appendix: Formulas
+
+## 4PL regression
+
+$$
+Y = D + \frac{A-D} {1+(\frac{C} {x})^B } + \epsilon
+$$
+
+where: x ... concentration of the analyte
+
+A: upper asymptote
+
+B: slope
+
+D: lower asymptote
+
+C ... EC50
+
+
+## log-logistic 4P regression
+
+$$
+Y = D + \frac{A-D} {1+e^{(B*(C - log(x))) }} + \epsilon
+$$
+
+
+
+## Intercept for slope at EC50
+
+$$
+I = A+\frac{D-A}{2}-B_{true}*EC50
+$$
+
+## Slope at EC50
+
+$$
+B_{true}=B*\frac{D-A}{4}
+$$
+
+## Confidence intervals
+
+In general, the confidence intervals are calculated as follows:
+$$
+ CI = \hat\theta\pm se(\hat\theta)*q^{t_{n-p}}_{1-\frac{\alpha}{2}}
+$$
+…where $\hat\theta$ is a fitted parameter or a linear combination thereof, q is the 1-alpha/2 quantile of the Student’s t-distribution with n-p degrees of freedom and se is the standard error derived from any covariance matrix.
+
+Let $\theta$ be the 4+1 parameters of the fit (a, b, d, EC50 of reference and EC50 difference). It can be shown that the least squares estimator $\hat\theta$ is normally distributed with asymptotic covariance matrix. The gradient method provides one of several ways to calculate the covariance matrix:
+
+$$
+\hat{V(\theta)}= \sigma^2(A(\hat\theta)^T*A(\hat\theta))^{-1}
+$$
+where A($\theta$) is the n x p matrix of the first partial derivatives for each parameter (i.e. gradient) realized at the fitted parameter estimates. The RMSE of the model or the pure error is used as estimate of $\sigma$. The square root of the diagonals of $\hat{V(\theta)}$ gives the standard errors and with that confidence intervals (CI) can be computed.
+
+# Literature
+
+[1] Finney, D.J.: (1978) Statistical Method in Biological Assay, London: Charles Griffin House, 3rd edition (pp. 80-82)
+
+[2] Franz, V.H.: Ratios: A short guide to confidence limits and proper use. arXiv:0710.2024v1, 10 Oct 2007
+
+[3] VerHoef, J.M.: Who invented the Delta Method? The American Statistician, 2012, 66:2, 124-127 DOI: 10.1080/00031305.2012.687494
+
+[4] Bates, D.M., Watts, D.G. (1988). Comparing models. In: Nonlinear Regression Analysis and Its Applications. New York: Wiley, pp 103-108
+
+[5] Bates, D.M., Watts, D.G. (1988) 2. In: Nonlinear Regression Analysis and Its Applications. New York: Wiley, pp 52-58
+
+[6] Motulsky, Brown Outlier testing

dev/ROUT.R

@@ -0,0 +1,267 @@
+################################################################################
+# F.Innerbichler
+# Jun 2026
+# Robust Outlier detection 
+################################################################################
+
+# Code acc. to Mutulsky and Brown hpptps://pubmed.ncbi.nlm.nih.gov/16526949/
+# except line 175 p <- 2*(1-pt(z, length(residuals)-5)) was exchanged with p <- 2*(1-pcauchy(z)) to be robust
+
+library(minpack.lm)
+library(ggplot2)
+library(rstudioapi)
+
+# Getting path for current open file
+current_path = rstudioapi::getActiveDocumentContext()$path
+setwd(dirname(current_path))
+
+outsPlot_FUN <- function(all_l,OUTs_, TS, PROC = "not specified", PROBE="", Q=0.01, par) {
+#browser()
+  #all_l_rout <- all_l[KEEP,]
+  klein <- min(log(all_l$conc))
+  gross <- max(log(all_l$conc))
+  x_seq <- seq(klein, gross, (gross-klein)/100)
+  samTRUE <- f4pl(x=x_seq, bottom=par["a"], top=par["d"],hill=par["b"], logEC50=par["cs"]-par["r"])
+  refTRUE <- f4pl(x=x_seq, bottom=par["a"], top=par["d"],hill=par["b"], logEC50=par["cs"])
+  
+  pl_T <- data.frame(cbind(x_seq, refTRUE, samTRUE))
+  p <- ggplot(all_l) +
+    geom_point(aes(x=log(conc), y=y, shape=factor(isRef)), size=2) +
+    theme_bw()
+
+  p2 <- p + geom_point(OUTs_, mapping=aes(x=log(conc), y=y, shape=factor(isRef)),
+                       size=5, color="violetred", stroke=2, inherit.aes = FALSE) +
+    scale_shape_manual(label=c("R",TS), values = c(21,25)) +
+    ggtitle(paste(PROC, PROBE, TS, "Threshold=", Q)) +
+    theme(legend.title=element_blank())
+  p3 <- p2 + geom_line(data=pl_T, aes(x=x_seq, y=refTRUE), color="blue", inherit.aes = F) +
+    geom_line(data=pl_T, aes(x=x_seq, y=samTRUE), color="red", inherit.aes = F)
+  p3
+  
+}
+
+#### 4PL model ----
+f4pl <- function(x, bottom, top, logEC50, hill) {
+  bottom+(top-bottom)/(1+exp((logEC50-x)*hill))
+}
+
+#### residuals for joint REF/TEST fit ----
+resid_4pl_joint <- function(par,x,y, is_ref) {
+  bottom <- par["bottom"]
+  top <- par["top"]
+  hill <- par["hill"]
+  le50_r <- par["logEC50_ref"]
+  le50_t <- par["logEC50_test"]
+  
+  le50 <- ifelse(is_ref, le50_r, le50_t)
+  yhat <- f4pl(x, bottom, top, le50, hill)
+  y - yhat
+}
+
+
+jac_4pl_joint <- function(par, x,y,is_ref) {
+  bottom <- par["bottom"]
+  top <- par["top"]
+  hill <- par["hill"]
+  le50_r <- par["logEC50_ref"]
+  le50_t <- par["logEC50_test"]
+  
+  le50 <- ifelse(is_ref, le50_r, le50_t)
+  
+  # useful intermediates
+  
+  texp <- exp((le50-x)*hill) # t = 10^((le50-x)*hill))
+  den <- (1-texp)
+  frac <- (top-bottom)/den
+  
+  dyhat_dbottom <- 1-1/den
+  dyhat_dtop <- 1/den
+  
+  d_invden_dle50 <- -(texp*log(2.718282)*hill)/(den^2)
+  
+  d_invden_dhill <- -(texp*log(2.718282)*(le50-x))/(den^2)
+  
+  dyhat_dle50 <- (top-bottom)* d_invden_dle50
+  dyhat_dhill <- (top-bottom)* d_invden_dhill
+  
+  J <- cbind(
+    bottom = -dyhat_dbottom,
+    topm = -dyhat_dtop,
+    hill = -dyhat_dhill,
+    logEC50_ref = ifelse(is_ref, -dyhat_dle50,0),
+    logEC50_test = ifelse(!is_ref, -dyhat_dle50,0)
+  )
+  J
+}
+
+
+robust_fit_irls <- function(par_start, x,y, is_ref,
+                            max_outer = 30, control = nls.lm.control(maxiter=200, ftol=1e-10, ptol=1e-10)) {
+  par <- par_start
+ #browser() 
+  for (k in seq_len(max_outer)) {
+    r <- resid_4pl_joint(par,x,y,is_ref)
+    
+    # Mot Brown
+    rSD <- quantile(r, 0.6827)*length(is_ref)/(length(is_ref)-5)
+    # Robust scale (MAD); fallback if MAD ~0
+    # s <- median(abs(r))/0.6745
+    
+    # Lorentzian / Cauchy weights acc. Mot/Brown
+    w <- 1/log(1+(abs(r)/rSD)^2)
+    # weighted residual function lor LM: sqrt(w)*r
+    fn_w <- function(p) {
+      rr <- resid_4pl_joint(p,x,y,is_ref)
+      sqrt(w)*rr
+      # w*rr
+    }
+    jac_w <- function(p) {
+      JJ <- jac_4pl_joint(p,x,y,is_ref)
+      JJ*sqrt(w) # row-wise scaling
+      # JJ*w
+    }
+    
+    #browser()
+    fit <- tryCatch(nls.lm(par=par, fn=fn_w, jac=jac_w, control=control),
+                    error=function(e) {
+                      paste("error at robust fit irls", k)
+                      # tau <<- tau *10
+                      #Return zero delta to avoid NaN updates
+                      return(rep(0, length(par)))
+                    })
+    
+    if (max(abs(fit$par -par)) < 1e-8) {
+      par <- fit$par
+      break
+    }
+    par <- fit$par
+  }
+  r_final <- resid_4pl_joint(par,x,y,is_ref)
+  # s_final <- median(abs(r_final))/0.6745
+  rSD <- quantile(r, 0.6827)*length(r)/(length(r)-5)
+  # Robust scale (MAD); fallback if MAD ~0
+  # s <- median(abs(r))/0.6745
+  list(par=par, residuals=r_final, scale=rSD)
+  
+}
+
+ # ROUT outlier detection via BH-FDR [1] https://cran.r-project.org/web//packages/minpack.lm/refman/minpack.lm.html)[6]
+
+rout_detect <- function(residuals, scale, Q=0.01) {
+  z <- abs(residuals/scale)
+  p <- 2*(1 - pcauchy(z))
+  # p <- 2*(1-pt(z, length(residuals)-5))
+  
+  o <- order(p)
+  p_sorted <- p[o]
+  m <- length(p)
+  # thresh <- (seq_len(m)/m) # *Q
+  keep <- p_sorted <= Q # thresh
+  (out_idx <- o[keep])
+  sort(out_idx)
+  
+}
+
+# Final OLS LM refit after removing outliers
+
+ols_refit <- function(par_start, x,y,is_ref, keep_idx,
+                      control = nls.lm.control(maxiter=400, ftol=1e-12, ptol=1e-12)) {
+  xk <- x[keep_idx]; yk <- y[keep_idx]; rk <- is_ref[keep_idx]
+  
+  fn_w <- function(p) resid_4pl_joint(p,xk,yk,rk)
+  jac_w <- function(p) jac_4pl_joint(p,xk,yk,rk)
+#browser()  
+  fit <- tryCatch(nls.lm(par=par_start, fn=fn_w, jac=jac_w, control=control),
+                  error=function(e) {
+                    paste("error at ols_refits")
+                    # tau <<- tau *10
+                    #Return zero delta to avoid NaN updates
+                    return(rep(0, length(par)))
+                  })
+  
+  list(par = fit$par, residuals = fn_w(fit$par), keep=keep_idx)
+}
+
+#. high-level ROUT 4PL analysis
+
+rout_4pl_potency <- function(df, Q=0.01, par_start=NULL, max_outer = 30) {
+  stopifnot(all(c("sample","conc","y") %in% names(df))) # "rep",
+ #browser() 
+  x <- log(df$conc)
+  y <- df$y
+  is_ref <- df$isRef == 1
+  
+  # start values if not provided
+  if (is.null(par_start)) {
+    bottom0 = min(y, na.rm = T)
+    top0 = max(y, na.rm = T)
+    if (df$y[1]>df$y[8]) hill0 <- 1 else hill0 <- -1
+    le50_r0 <- median(x[is_ref], na.rm = T)
+    le50_t0 <- median(x[!is_ref], na.rm = T)
+    par_start <- c(bottom=bottom0, top=top0, hill=hill0, logEC50_ref = le50_r0 , logEC50_test = le50_t0)
+  }
+  # Step1 : robust fit
+  robust <- robust_fit_irls(par_start, x,y,is_ref, max_outer = max_outer)
+  # Step 2: outlier detection
+  out_idx <- rout_detect(robust$residuals, robust$scale, Q=Q)
+  keep_idx <- setdiff(seq_len(nrow(df)), out_idx)
+  
+  # Step 3: OLS refit
+  ols <- ols_refit(robust$par, x,y,is_ref, keep_idx)
+  
+  # Relative potency
+  rp <- exp(ols$par["logEC50_ref"] - ols$par["logEC50_test"])*100
+  
+  list(Q=Q, outliers=out_idx,
+       kept=keep_idx,
+       par_robust = robust$par,
+       par_final = ols$par,
+       rp_percent = rp,
+       df=df)
+  
+}
+
+# sheets <- openxlsx::getSheetNames("~/plateflow/outlierFUB.xlsx")
+# all_dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile="~/plateflow/outlierFUB.xlsx")
+# names(all_dat) <- sheets
+# 
+# Plate <- 1
+# 
+# PlanteN <- sheets[Plate]
+# NoDils <- 8
+# Nreps <- 3
+# DAT <- all_dat[[Plate]]
+# colnames(DAT) <- DAT[2,]
+# REF <- DAT[3:10,1:4]
+# SAM1 <- DAT[13:20,1:4]
+# SAM2 <- DAT[23:30,1:4]
+# SAM3 <- DAT[33:40,1:4]
+# 
+# REF_ <- sapply(REF, as.numeric)
+# colnames(REF_) <- colnames(REF)
+# SAM1_ <- sapply(SAM1, as.numeric)
+# colnames(SAM1_) <- colnames(SAM1)
+# SAM2_ <- sapply(SAM2, as.numeric)
+# colnames(SAM2_) <- colnames(SAM2)
+# SAM3_ <- sapply(SAM3, as.numeric)
+# colnames(SAM3_) <- colnames(SAM3)
+# 
+# plot_df <- rbind(REF_, SAM1_, SAM2_,SAM3_)
+# 
+# plot_df_ <- cbind(as.data.frame(plot_df), sample=c(rep("R", NoDils),rep("S1R", NoDils),rep("S2", NoDils),rep("S3", NoDils)))
+# test_df <- plot_df_[1:16,]
+# test_df2 <- plot_df_[c(1:8, 17:24),]
+# test_df3 <- plot_df_[c(1:8, 25:32),]
+# 
+# all_l <- melt(test_df, id.vars = c("sample","Dose"), variable.name="replname", value.name="readout")
+# colnames(all_l) <- c("sample","conc","rep","y")
+# 
+# for (Q in c(0.01,0.015,0.02)) {
+#   res <- rout_4pl_potency(all_l, Q)
+#   OUTs_ <- all_l[res$outliers,]
+#   print(outsPlot_FUN(all_l, OUTs_, TS=all_l$sample[9], PROC="ROUT",PROBE="Sheet 1",Q,par=res$par_final))
+# }
+
+
+
+

dev/app.R

@@ -28,9 +28,10 @@ library(twopartm)
 library(car)
 library(dplyr)
 library(scales)
+library(tolerance)
 
 source("../R/Global.R")
-
+source("ROUT.R")
 
 #### ui ----
 
@@ -116,9 +117,33 @@ server <- function(input, output, session) {
           "It needs to contain 1 column with the  dilution concentrations (first or last column) and at least 2 columns of reference  and test sample readouts, respectively.",
           "The reference readout columns have to be before the test sample readout columns. The column names for reference and test are free to set, but different for all columns.",
           "The column name of the dilution concentrations can be as follows: concentration, dose, log_concentration, log_dose (first letter can be capital)",
-          "It is assumed, that the concentrations are in anti-log or in natural log mode.",
+          "If the concentrations are in logarithmized, any log base can be used.",
+          br(), br(),
+          "EXPLORE the 4pl function: visualize the meta data in the context of a 4 PL fit or a linear regression fit. ",
+          "Enter the 4 parameters of test and reference sample and see, what this means.", br(),
+           br(),
+          "OPTIMIZE the concentrations: plot all results you have and adjust the concentrations accordingly. ",
+          "Get help if you want to read pdfs in contacting us.", br(),
         ),
-        column(6, )
+        column(2,
+               style = "background: #7FAEFF88",
+               "HERE: Enter the equivalence limits for 4PL suitability tests. If you need help to set them, contact us.", 
+               numericInput("lEACratiola", "lower EAC ratio of LAs", 0.005, step = 0.001),
+               numericInput("uEACratiola", "upper EAC for ratio of LAs", 100, step = 1),
+               numericInput("lEACratioSlope", "lower EAC for ratio of slopes", 0.55, step = 0.01),
+               numericInput("uEACratioSlope", "upper EAC for ratio of slopes", 1.84, step = 0.1),
+               numericInput("lEACratioua", "lower EAC for ratio of UAs", 0.75, step = 0.1),
+               numericInput("uEACratioua", "upper EAC for ratio of UAs", 1.33, step = 0.1)
+        ),
+        column(2,
+               style = "background: #7FAEFF88",
+               numericInput("lowerPot", "lower EAC for potency", 75, step = 1),
+               numericInput("upperPot", "upper EAC for potency", 133, step = 1),
+               numericInput("lEACratioAdiff", "lower EAC of ratio of asymptote differences", 0.75, step = 0.01),
+               numericInput("uEACratioAdiff", "upper EAC of ratio of asymptote differences", 1.33, step = 0.01),
+               numericInput("lEACdiffla", "lower EAC for diff. of LA", -0.175, step = 0.001),
+               numericInput("uEACdiffla", "upper EAC for diff. of LA", 0.189, step = 0.001)
+        )
       ),
       tabPanel(
         "Documentation",
@@ -152,17 +177,17 @@ server <- function(input, output, session) {
               ),
               uiOutput(outputId = "sheetName"),
               "For data format in the EXCEL file see Data template",
-              "If no data are uploaded, the settings to the right are used for calculations.",
               tags$head(tags$style(HTML("label {font-size:80%;margin-bottom: 3px;margin-top: 3px;}"))),
               div(checkboxInput("PureErr", "Should pure error be used for calculation of CIs?", FALSE),
                 style = "font-size: 24px !important;color: #C2173F"
               ),
 
               # actionLink("selectall","SelectAll"),
-              h5("\n\n\n Author: Franz Innerbichler, InnerAnalytics")
+              #h5("\n\n\n Author: Franz Innerbichler, InnerAnalytics")
             ),
             column(
               4,
+              
               h4("Suitability tests for 4-parametric logistic regression"),
               "(potency CI test is set per default)",
               checkboxGroupInput("selectedSSTs", "Which suitability tests to be used?",
@@ -186,26 +211,15 @@ server <- function(input, output, session) {
                 ),
                 selected = c("1", "2", "3", "4", "5", "6", "7", "8")
               )
-            ),
-            column(2,
-              style = "background: #7FAEFF88",
-              numericInput("lEACratiola", "lower EAC ratio of LAs", 0.005, step = 0.001),
-              numericInput("uEACratiola", "upper EAC for ratio of LAs", 100, step = 1),
-              numericInput("lEACratioSlope", "lower EAC for ratio of slopes", 0.55, step = 0.01),
-              numericInput("uEACratioSlope", "upper EAC for ratio of slopes", 1.84, step = 0.1),
-              numericInput("lEACratioua", "lower EAC for ratio of UAs", 0.75, step = 0.1),
-              numericInput("uEACratioua", "upper EAC for ratio of UAs", 1.33, step = 0.1)
-            ),
-            column(2,
-              style = "background: #7FAEFF88",
-              numericInput("lowerPot", "lower EAC for potency", 75, step = 1),
-              numericInput("upperPot", "upper EAC for potency", 133, step = 1),
-              numericInput("lEACratioAdiff", "lower EAC of ratio of asymptote differences", 0.75, step = 0.01),
-              numericInput("uEACratioAdiff", "upper EAC of ratio of asymptote differences", 1.33, step = 0.01),
-              numericInput("lEACdiffla", "lower EAC for diff. of LA", -0.175, step = 0.001),
-              numericInput("uEACdiffla", "upper EAC for diff. of LA", 0.189, step = 0.001)
             )
           ),
+          
+          tabPanel(
+            "Uploaded data",
+            tableOutput("XLdata")
+          ),
+          
+          ###### 4pl output ----
           tabPanel(
             "4pl-Analysis",
             tags$style(HTML("pre { color: black; background-color: #FFE1FF;
@@ -253,6 +267,7 @@ server <- function(input, output, session) {
               )
             )
           ),
+          ##### linear output ----
           tabPanel(
             "linear Analysis",
             sidebarLayout(
@@ -287,7 +302,7 @@ server <- function(input, output, session) {
                     )
                   ),
                   tabPanel(
-                    "Tests and ANOVAA",
+                    "Tests and ANOVA",
                     column(
                       12,
                       h3("Tests for linear PLA:"),
@@ -312,6 +327,19 @@ server <- function(input, output, session) {
               )
             )
           ),
+          ##### ROUT output ----
+          tabPanel("robust outlier testing",
+                   downloadButton("downloadROUTReport", label = "Download ROUT report", class = "butt"),
+                   sliderInput("Qslider", "adjust Q-value in %",min=0.1, max = 10, value = 1, step=0.1),
+                   plotOutput("OutlierPlot"),
+                   tableOutput("OutlierDF"),
+                   
+          "GUIDANCE: The procedure of Motulsky & Brown allows for robust outlier testing.",
+          "Adjust the slider to mark the suspected outliers. Mostly, a Q-value of 2% is sufficient.",
+          "If the general variability of the data is high, many datapoints will be flagged, also ones that are not deemed to be outliers",
+          "An indicator for high variability is, when the Q-value needs to be increased above 5%, to flag suspected outliers.",
+          "Then, please re-consider, if the suspected outlier is not 'just' normal variability."
+          ),
           tabPanel(
             "parameter estimates",
             htmlOutput("PureErrWParEst"),
@@ -376,7 +404,7 @@ server <- function(input, output, session) {
         mainPanel(
           width = 12,
           tabsetPanel(
-            id = "tabs",
+            id = "tabs2",
             tabPanel(
               "Settings",
               h4("Settings of 4PL regression"),
@@ -552,8 +580,8 @@ server <- function(input, output, session) {
             tabPanel(
               "Report",
               h4("Settings for report"),
-              downloadButton("downloadXLReport", label = "Download PDF report", class = "butt"),
+              downloadButton("downloadXLReportMeta", label = "Download PDF report", class = "butt"),
-              tags$style(type = "text/css", "#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"),
+              tags$style(type = "text/css", "#downloadXLReportMeta {background-color: orange; color: black;font-family: COurier New}"),
             )
           )
         )
@@ -588,12 +616,15 @@ server <- function(input, output, session) {
                   fileInput("MiFile", "", accept = ".xlsx")
               ),
               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)
+              #checkboxInput("fixupper","Fix highest concentration (if unticked, the center is fixed)",FALSE),
+              sliderInput("dilmover", "Move the dilutions(+-move in log-units)", min = -3,max=3, value=0, step=0.1, round=1),
+              numericInput("TolConf","confidence", value=0.95, step=0.01),
+              numericInput("TolPop","population", value=0.9, step=0.01)
             )
           ),
           mainPanel(
             tabsetPanel(
-              id = "tabs",
+              id = "tabs3",
               tabPanel("4pl",
                 
 
@@ -620,9 +651,13 @@ server <- function(input, output, session) {
                        "Narrower dilution ranges decrease the CIs of rel. potency, and increase the CIs of upper and lower asymptote ratios, ands Hill's slope ratios",
 
               ),
-              tabPanel("Histograms",
+              tabPanel("Plots",
                 h4("Histograms of parameters"),
-                plotOutput("histCIs"),
+                plotOutput("linerangeCIs"),
+                plotOutput("ratioSlope"),
+                plotOutput("ratioLas"),
+                plotOutput("ratioUas"),
+                plotOutput("widthCIs"),
                 column(6,
                       plotOutput("histEC50REF"), 
                       plotOutput("histLasREF"),
@@ -637,7 +672,16 @@ server <- function(input, output, session) {
               ),
               tabPanel(
                 "Report",
-                h4("Settings for report"))
+                h4("Settings for report"),
+                #useShinyjs(),
+                #actionButton("btn2", "Download PDF report", icon = icon("download")),
+                downloadButton("downloadWizardData", label = "Download model and plots", class = "butt"),
+                tags$style(type = "text/css", "#downloadWizardData {background-color: #4FCBD9; color: black;font-family: Courier New}"),
+                # textInput("Author", "Author", value = ""),
+                # textInput("RepIdentifier", "Report name", value = ""),
+                # textInput("NoP", "Product name", value = ""),
+                # textInput("Assay", "Assay name", value = "")
+                )
               
             )
           ) # main panel
@@ -700,7 +744,7 @@ server <- function(input, output, session) {
     reset(id = "") # from shinyjs package
   })
 
-#### input optim XL file ----
+#### input Wizard XL file ----
   observe({
     if (!is.null(input$MiFile)) {
       MinFile <- input$MiFile
@@ -736,17 +780,19 @@ server <- function(input, output, session) {
           if (length(logI) > 0 & length(logDoseI) == 0) {
             XLdat$log_dose <- XLdat[, logI]
             XLdat2 <- XLdat[, -logI]
-            CORro <- cor(XLdat$log_dose, XLdat[, 3])
+            CORro <- COR_FUNC(XLdat$log_dose, XLdat[, 3])
           } else if (length(logI) == 0 & length(logDoseI) == 0) {
             Ind <- grep(".ilution|.ose|.onc", cn)
             XLdat$log_dose <- log(XLdat[, Ind])
-            CORro <- cor(XLdat[, Ind], XLdat[, 3])
+            CORro <- COR_FUNC(XLdat[, Ind], XLdat[, 3])
             XLdat2 <- XLdat[, -Ind]
           } else if (length(logI) > 0 & length(logDoseI) > 0) {
             XLdat2 <- XLdat
-            CORro <- cor(XLdat[, logI], XLdat[, 3])
+            CORro <- COR_FUNC(XLdat[, logI], XLdat[, 3])
           }
           Dat$EXCEL <- XLdat2
+          output$XLdata <- renderTable({ XLdat2 })
+          
           PureErrFlag <- input$PureErr
           warning_text2 <- reactive({
             ifelse(PureErrFlag, "Pure Error is selected", "")
@@ -773,9 +819,82 @@ server <- function(input, output, session) {
           # all_l$readout[all_l$readout < 0] <- 0.01
           REP$all_l <- all_l
 
-          #### XLSX eval ----
+          ##### ROUT outlier testing ----
+    #browser()
+          if(!is.null(input$Qslider)) {
+            all_lROUT <- all_l[complete.cases(all_l),]
+            colnames(all_lROUT) <- c("log_dose","sample","y","isRef","isSample","conc")
+            res <- rout_4pl_potency(all_lROUT, Q=input$Qslider/100)
+            OUTs_ <- all_lROUT[res$outliers,]
+            
+            all_l_rout <- all_lROUT[res$kept,]
+            # all_l_rout$log_dose <- log(all_l_rout$Conc)
+            # colnames(all_l_rout) <- c("log_dose","sample","y","isRef","isSample","conc")
+            if (all_l_rout$conc[1]>all_l_rout$conc[6]) {
+              if (all_l_rout$y[1]>all_l_rout$y[6]) SLOPE <- 1 else SLOPE<- -1
+            } else {
+              if (all_l_rout$y[1]>all_l_rout$y[6]) SLOPE <- -1 else SLOPE<- 1
+            }
+
+            startlist <- list(a = min(all_l_rout$y), b = SLOPE, d = max(all_l_rout$y), cs = mean(log(all_l_rout$conc)), r = 0)
+            
+            mr <- tryCatch(
+              {
+                gsl_nls(
+                  fn = y ~ a + (d - a) / (1 + exp(b * ((cs - r * isSample) - log_dose))),
+                  data = all_l_rout,
+                  start = startlist, # race=T,
+                  control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)
+                )
+              },
+              warning = function(e) {
+                mr <<- "In nlsModel singular gradient matrix"
+                
+              })
+            PAR <- summary(mr)$coefficients[,1]
+            ROUTplot <- outsPlot_FUN(all_l_rout, OUTs_, TS=all_l_rout$sample[13], PROC="ROUT",PROBE=input$sheet,Q=input$Qslider/100,par=PAR)
+        
+            
+            output$OutlierDF <- renderTable({
+             OUTs_
+            })
+            output$OutlierPlot <- renderPlot({
+              print(ROUTplot)
+            })
+            
+            
+            all_l_rout2 <- all_l_rout[,-c(4:6)]
+            ro_newROUT <- reshape(all_l_rout2, direction="wide", idvar = "log_dose", timevar="sample", varying = as.vector(unique(all_l_rout2$sample)))
+            REP$ro_newROUT <- ro_newROUT
+            
+            REP$ROUTplot <- ROUTplot
+            
+            ANOVA_ROUT <- ANOVA4plUnresfunc(ro_new = ro_newROUT)
+            REP$ANOVA_ROUT <- ANOVA_ROUT
+            
+            Limite <- list(
+              as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla),
+              as.numeric(input$lEACratiola), as.numeric(input$uEACratiola),
+              as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope),
+              as.numeric(input$lEACratioua), as.numeric(input$uEACratioua),
+              as.numeric(input$lowerPot), as.numeric(input$upperPot),
+              as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff)
+            )
+            
+            tabROUT <- tests_FUNC(ro_newROUT, Limite, PureErrFlag = PureErrFlag)
+       #browser()     
+            tabROUT[1, 6:7] <- c("-", "-")
+            
+            #tabROUT2 <- tabROUT[SelTests, ]
+            #Dat$tests_FUNC <- tabROUT
+            REP$testsTabROUT <- tabROUT
+            
+            }          
+
+  
+          ##### XLSX eval ----
           #if (CORro < 0) SLOPE <- -1 else SLOPE <- 1
-          FITs <- Fitting_FUNC(XLdat2, TransFlag = FALSE)
+          FITs <- Fitting_FUNC(XLdat2, TransFlag = FALSE, nameWS="")
 
           #### if no 4pl fit is possible ----
           if (!is.null(FITs)) {
@@ -795,7 +914,7 @@ server <- function(input, output, session) {
               })
 
               warning_textNo4PLFit <- reactive({
-                ifelse(Dat$FITsFlag, "No meaningful 4PL fit was possible. This may havea several reasons: \nA control sample was tested/\n
+                ifelse(Dat$FITsFlag, "No meaningful 4PL fit was possible. This may have a several reasons: \nA control sample was tested/\n
                        the EC50 is not catched with the dilutions/\n the assay/reader had a problem",
                   "Footnote: bendpoints (linear part) and asymptote points (point where asymptote is reached) are plotted in dashed and dotted lines. They indicate whether the linear part and asymptotes are catched with the current dilutions.
                        Black line is the true slope at EC50 of REF."
@@ -1759,17 +1878,18 @@ server <- function(input, output, session) {
       slopeTe[i, ] <- lm3Te$coefficients
     }
 
-    indS <- which(abs(slopeSt[, 2]) == max(abs(slopeSt[, 2])))
+    indS <- which(abs(slopeSt[, 2]) == max(abs(slopeSt[, 2]), na.rm=T))
-    indT <- which(abs(slopeTe[, 2]) == max(abs(slopeTe[, 2])))
+    indT <- which(abs(slopeTe[, 2]) == max(abs(slopeTe[, 2]), na.rm = T))
 
     # pl_ <- slopeSt[indS,1]+slopeSt[indS,2]*log_conc
     # pl_T <- slopeTe[indT,1]+slopeTe[indT,2]*log_conc
     # pl_df <- data.frame(lnC=log_conc, plotS=pl_, plotT=pl_T)
-
+#browser()
     all_l <- melt(data.frame(tab), 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)
+    all_l2 <- all_l2[complete.cases(all_l2),]
     all_l2S <- all_l2[all_l2$isRef == 1, ]
     all_l2T <- all_l2[all_l2$isRef == 0, ]
     all_mS <- all_l2S[order(all_l2S$log_dose, decreasing = TRUE), ]
@@ -1991,7 +2111,7 @@ server <- function(input, output, session) {
     pottab4_$`upper95%CI` <- round(as.numeric(pottab4[, 4]) * 100, 2)
     pottab4_$relative_lowerCL <- round(pottab4_[, 6] / pottab4_[, 5] * 100, 2)
     pottab4_$relative_upperCL <- round(pottab4_[, 7] / pottab4_[, 5] * 100, 2)
-
+#browser()
     if (as.numeric(pottab4_$relative_lowerCL[1]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[1]) < Lim[[10]]) {
       test_potCI <- 0
     } else {
@@ -2117,7 +2237,7 @@ server <- function(input, output, session) {
   })
 
 
-  #### Dilutions Simulator ----
+  #### Meta plots all XL ---- 
   observe({ 
     if (!is.null(Dat$Mws)) {
 
@@ -2130,59 +2250,61 @@ server <- function(input, output, session) {
     for (N_WS in 1:length(AllXL)) {
    
       datWS <- as.data.frame(AllXL[[N_WS]])
-      
+      nameWS <- names(AllXL)[N_WS]
       cn <- colnames(datWS)
       logI <- grep("log|ln", cn)
       logDoseI <- grep("log_dose", cn)
       if (length(logI) > 0 & length(logDoseI) == 0) {
         datWS$log_dose <- datWS[, logI]
         datWS2 <- datWS[, -logI]
-        CORro <- cor(datWS$log_dose, datWS[, 3])
+        CORro <- COR_FUNC(datWS$log_dose, datWS[, 3])
       } else if (length(logI) == 0 & length(logDoseI) == 0) {
         Ind <- grep(".ilution|.ose|.onc", cn)
         datWS$log_dose <- log(datWS[, Ind])
-        CORro <- cor(datWS[, Ind], datWS[, 3])
+
+        CORro <- COR_FUNC(datWS[, Ind], datWS[, 3])
         datWS2 <- datWS[, -Ind]
       } else if (length(logI) > 0 & length(logDoseI) > 0) {
         datWS2 <- datWS
-        CORro <- cor(datWS[, logI], datWS[, 3])
+        CORro <- COR_FUNC(datWS[, logI], datWS[, 3])
       }
       Dat$datWS2 <- datWS2
       
-      FITs <- Fitting_FUNC(datWS2, TransFlag = F)
+      FITs <- Fitting_FUNC(datWS2, TransFlag = F, nameWS = nameWS)
+      if (!is.character(FITs)) {
+        pot_est <- FITs[[3]]
+        potEstL[[N_WS]] <- pot_est
+        potU_est <- FITs[[4]]
+        # unrestricted
+        SU_mu <- FITs[[2]]
+        URMcoefs1 <- SU_mu$coefficients
+        URMcoefs <- t(matrix(unlist(URMcoefs1[,1])))
+        URMcoefs_ <- cbind(AllSheets[[N_WS]], URMcoefs)
+        URMcoefsL[[N_WS]] <- URMcoefs_
         
-      pot_est <- FITs[[3]]
+        SU_mr <- FITs[[1]]
-      potEstL[[N_WS]] <- pot_est
+        RMcoefs1 <- SU_mr$coefficients
-      potU_est <- FITs[[4]]
+        RMcoefs <- t(matrix(unlist(RMcoefs1[,1])))
-      # unrestricted
+        RMcoefs_ <- cbind(AllSheets[[N_WS]], RMcoefs)
-      SU_mu <- FITs[[2]]
+        RMcoefsL[[N_WS]] <- RMcoefs_
-      URMcoefs1 <- SU_mu$coefficients
-      URMcoefs <- t(matrix(unlist(URMcoefs1[,1])))
-      URMcoefs_ <- cbind(AllSheets[[N_WS]], URMcoefs)
-      URMcoefsL[[N_WS]] <- URMcoefs_
-      
-      SU_mr <- FITs[[1]]
-      RMcoefs1 <- SU_mr$coefficients
-      RMcoefs <- t(matrix(unlist(RMcoefs1[,1])))
-      RMcoefs_ <- cbind(AllSheets[[N_WS]], RMcoefs)
-      RMcoefsL[[N_WS]] <- RMcoefs_
         
         
         
-      X <- seq(min(datWS2$log_dose), max(datWS2$log_dose), 0.1)
+        X <- seq(min(datWS2$log_dose), max(datWS2$log_dose), 0.1)
-      sigRef <- URMcoefs[1,1] + (URMcoefs[1,3]-URMcoefs[1,1])/(1+exp(URMcoefs[1,2]*(URMcoefs[1,4]-X)))
+        sigRef <- URMcoefs[1,1] + (URMcoefs[1,3]-URMcoefs[1,1])/(1+exp(URMcoefs[1,2]*(URMcoefs[1,4]-X)))
-      sigTest1 <- URMcoefs[1,5] + (URMcoefs[1,7]-URMcoefs[1,5])/(1+exp(URMcoefs[1,6]*(URMcoefs[1,4] - URMcoefs[1,8]-X)))
+        sigTest1 <- URMcoefs[1,5] + (URMcoefs[1,7]-URMcoefs[1,5])/(1+exp(URMcoefs[1,6]*(URMcoefs[1,4] - URMcoefs[1,8]-X)))
-      #browser()  
+        #browser()  
-      dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, Sheet = AllSheets[[N_WS]])
+        dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, Sheet = AllSheets[[N_WS]])
-      dfPlotsigTest <- data.frame(X=X, sigTest = sigTest1, Sheet = AllSheets[[N_WS]])
+        dfPlotsigTest <- data.frame(X=X, sigTest = sigTest1, Sheet = AllSheets[[N_WS]])
-      
-      if (!exists("SIGrefDF")) SIGrefDF <- dfPlotsigRef else SIGrefDF <- rbind(SIGrefDF, dfPlotsigRef)
-      if (!exists("SIGtestDF")) SIGtestDF <- dfPlotsigTest else SIGtestDF <- rbind(SIGtestDF,dfPlotsigTest)
         
+        if (!exists("SIGrefDF")) SIGrefDF <- dfPlotsigRef else SIGrefDF <- rbind(SIGrefDF, dfPlotsigRef)
+        if (!exists("SIGtestDF")) SIGtestDF <- dfPlotsigTest else SIGtestDF <- rbind(SIGtestDF,dfPlotsigTest)
+      }
     } #for N_WS
     
-    #browser()   
+#browser()   
     URMcoefsDF <- t(matrix(unlist(URMcoefsL),nrow=9))
+    colnames(URMcoefsDF) <- c("WS name", "lowerAs REF","slope REF","upperAs REF","EC50 REF", "lowerAs TEST","slope TEST","upperAs TEST","EC50 Difference")
     EC50TEST <- as.numeric(URMcoefsDF[,5]) - as.numeric(URMcoefsDF[,9])
     # EC50TEST <- EC50TEST[!EC50TEST %in% boxplot.stats(EC50TEST)$out]
     EC50REF <- as.numeric(URMcoefsDF[,5])
@@ -2191,18 +2313,28 @@ server <- function(input, output, session) {
     # UasREF <- UasREF[!UasREF %in% boxplot.stats(UasREF)$out]
     LasREF <- as.numeric(URMcoefsDF[,2])
     # LasREF <- LasREF[!LasREF %in% boxplot.stats(LasREF)$out]
-    UasTEST <- as.numeric(URMcoefsDF[,4])
+    UasTEST <- as.numeric(URMcoefsDF[,8])
-    LasTEST <- as.numeric(URMcoefsDF[,2])
+    LasTEST <- as.numeric(URMcoefsDF[,6])
+    slopeREF <- as.numeric(URMcoefsDF[,3])
+    slopeTEST <- as.numeric(URMcoefsDF[,7])
+    
+    slopeRatio <- slopeTEST/slopeREF
+    LasRatio <- LasTEST/LasREF
+    UasRatio <- UasTEST/UasREF
+    ratioDF <- data.frame(WS_name = URMcoefsDF[,1], slopeRatio = slopeRatio, LasRatio = LasRatio, UasRatio = UasRatio)
     
     RMcoefsDF <- t(matrix(unlist(RMcoefsL),nrow=6))
-    
+    colnames(RMcoefsDF) <- c("WS name", "lower asymptote","Hill's slope","upper asymptote","log(EC50 ref)","logEC50 difference")
     Dat$URMcoefsDF <- URMcoefsDF
+    Dat$ModU <- URMcoefsDF
     Dat$RestrM <- RMcoefsDF
+    Dat$ModR <- RMcoefsDF
     
     CalcPotDF <-  t(matrix(unlist(potEstL),nrow=3))
+    colnames(CalcPotDF) <- c("rel_potency","lower_CI","upper_CI")
     Dat$CalcPot <- CalcPotDF
     # 
-    #### sigmoid plots ----
+    #### Wizard sigmoid plots ----
     
     Slope <- as.numeric(URMcoefsDF[1,3])
     if (Slope > 0) {
@@ -2211,11 +2343,15 @@ server <- function(input, output, session) {
     
     #browser()
     BoxDF <- data.frame(EC50REF = EC50REF, EC50TEST = EC50TEST, LasREF = LasREF, UasREF = UasREF)
+    UasParTolREF <- normtol.int(x = UasREF, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
+    LasParTolREF <- normtol.int(x = LasREF, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
     
     p1 <- ggplot(SIGrefDF, aes(x=X, y=sigRef, col=as.factor(Sheet))) +
       geom_line() +
       annotate("text", label="x", x=x_UA, y=UasREF, alpha=0.2) +
       annotate("text", label="o", x=x_LA, y=LasREF, alpha=0.2) +
+      geom_hline(yintercept = c(UasParTolREF[[4]], UasParTolREF[[5]]), linetype=2, col="grey") +
+      geom_hline(yintercept = c(LasParTolREF[[4]], LasParTolREF[[5]]), linetype=2, col="grey") +
       geom_vline(xintercept = EC50REF, alpha = 0.2) +
       scale_x_continuous(expand = c(0, 0)) + 
       scale_y_continuous(expand = c(0, 0)) + 
@@ -2224,78 +2360,27 @@ server <- function(input, output, session) {
       expand_limits(x = c(min(SIGrefDF$X) - 0.1 * diff(range(SIGrefDF$X)), 
                           max(SIGrefDF$X) + 0.1 * diff(range(SIGrefDF$X)))) + 
       xlab("dilutions") +
-      #ggtitle("Plot of all calculated reference fits (unrestricted model, in gray vertical lines: EC50)") +
+      ggtitle("REF sample 4PL-fits (unrestricted model, gray vertical lines: EC50)") +
       theme_bw() +
       theme(axis.text = element_text(face = "bold", size = 15),
             plot.title = element_text(size = 15, face = "bold"),
             plot.margin = unit(c(0.2, 0.2, 0.5, 0.5), "lines"))
-    # Horizontal marginal boxplot - to appear at the top of the chart
-    pBox_hor <- ggplot( BoxDF, aes(x = factor(1), y = EC50REF)) + 
-      geom_boxplot(outlier.colour = NA) + 
-      geom_jitter(position = position_jitter(width = 0.05)) + 
-      scale_y_continuous(expand = c(0, 0)) + 
-      expand_limits(y = c(min(SIGrefDF$X) - 0.1 * diff(range(SIGrefDF$X)), 
-                          max(SIGrefDF$X) + 0.1 * diff(range(SIGrefDF$X)))) + 
-      coord_flip() + 
-      theme_bw() +
-      theme(axis.text = element_blank(), 
-            axis.title = element_blank(), 
-            axis.ticks = element_blank(), 
-            plot.margin = unit(c(1, 0.2, -0.5, 0.5), "lines"))
-    
-    # Vertical marginal boxplot - to appear at the right of the chart
-    pBox_ver <- ggplot(BoxDF, aes(x = factor(1), y = UasREF)) + 
-      geom_boxplot(outlier.colour = NA) + 
-      geom_jitter(position = position_jitter(width = 0.05)) + 
-      scale_y_continuous(expand = c(0, 0)) + 
-      expand_limits(y = c(min(SIGrefDF$sigRef) - 0.1 * diff(range(SIGrefDF$sigRef)), 
-                          max(SIGrefDF$sigRef) + 0.1 * diff(range(SIGrefDF$sigRef)))) + 
-      theme_bw() +
-      theme(axis.text = element_blank(), 
-            axis.title = element_blank(), 
-            axis.ticks = element_blank(), 
-            plot.margin = unit(c(0.2, 1, 0.5, -0.5), "lines"))
-    
-    #browser()
-    gt1 <- ggplot_gtable(ggplot_build(p1))
-    gt2 <- ggplot_gtable(ggplot_build(pBox_hor))
-    gt3 <- ggplot_gtable(ggplot_build(pBox_ver))
-    
-    # Get maximum widths and heights
-    maxWidth <- unit.pmax(gt1$widths[2:3], gt2$widths[2:3])
-    maxHeight <- unit.pmax(gt1$heights[4:5], gt3$heights[4:5])
-    
-    # Set the maximums in the gtables for gt1, gt2 and gt3
-    gt1$widths[2:3] <- as.list(maxWidth)
-    gt2$widths[2:3] <- as.list(maxWidth)
-    
-    gt1$heights[4:5] <- as.list(maxHeight)
-    gt3$heights[4:5] <- as.list(maxHeight)
-    # Create a new gtable
-    gt <- gtable(widths = unit(c(7, 1), "null"), height = unit(c(1, 7), "null"))
-    
-    # Instert gt1, gt2 and gt3 into the new gtable
-    gt <- gtable_add_grob(gt, gt1, 2, 1)
-    gt <- gtable_add_grob(gt, gt2, 1, 1)
-    gt <- gtable_add_grob(gt, gt3, 2, 2)
-    
-    # grid.rect(x = 0.5, y = 0.5, height = 0.995, width = 0.995, default.units = "npc", 
-    #           gp = gpar(col = "black", fill = NA, lwd = 1))
-    # And render the plot
-    grid.newpage()
-#browser()
-    
-    output$sigPlotREF <- renderPlot({ grid.draw(gt) })
     
+    output$sigPlotREF <- renderPlot({ p1 })
     Dat$sigPlotREF <- p1
-    # 
+    
+    UasParTolTEST <- normtol.int(x = UasTEST, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
+    LasParTolTEST <- normtol.int(x = LasTEST, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
+   
     p2 <- ggplot(SIGtestDF, aes(x=X, y=sigTest, col=as.factor(Sheet))) +
       geom_line() +
       annotate("text", label="x", x=x_UA, y=UasTEST, alpha=0.2) +
+      geom_hline(yintercept = c(UasParTolTEST[[4]], UasParTolTEST[[5]]), linetype=2, col="grey") +
       annotate("text", label="o", x=x_LA, y=LasTEST, alpha=0.2) +
+      geom_hline(yintercept = c(LasParTolTEST[[4]], LasParTolTEST[[5]]), linetype=2, col="grey") +
       geom_vline(xintercept = EC50TEST, alpha = 0.2) +
       xlab("dilutions") +
-      ggtitle("Calculated test sample fits (unrestricted model, in gray vertical lines: EC50)") +
+      ggtitle("TEST sample 4PL-fits (unrestricted model, gray vertical lines: EC50)") +
       theme_bw() +
       theme(axis.text = element_text(face = "bold", size = 15),
             plot.title = element_text(size = 15, face = "bold"))
@@ -2306,26 +2391,87 @@ server <- function(input, output, session) {
   
     #### histograms right panel ----
     
-    #browser()
    
-    all_lPot <- data.frame(Cat_potency= c(rep("rel poteny",nrow(CalcPotDF)), rep("lower CI",nrow(CalcPotDF)),rep("upper CI",nrow(CalcPotDF))),
+    
+    all_lPot <- data.frame(Cat_potency= c(rep("rel_poteny",nrow(CalcPotDF)), rep("lower_CI",nrow(CalcPotDF)),rep("upper_CI",nrow(CalcPotDF))),
       Potency_and_CI = c(CalcPotDF[,1], CalcPotDF[,2],CalcPotDF[,3]))
     all_lPot[,2][all_lPot[,2] > 5] <- NA
     all_lPot[,2][all_lPot[,2] < 0.1] <- NA
     
-    P_histCI <- ggplot(all_lPot, aes(x=Potency_and_CI, fill=Cat_potency)) +
+    CalcPotDF <- as.data.frame(CalcPotDF)
+    CalcPotDF$width_CI <- CalcPotDF$upper_CI - CalcPotDF$lower_CI
+    widthCLTol <- normtol.int(x = CalcPotDF$width_CI, alpha = 1-input$TolConf, P = input$TolPop, side = 1)
+    
+    P_linerangeCI <- ggplot(CalcPotDF, aes(x=seq(1,nrow(CalcPotDF)))) + #, aes(x=Potency_and_CI,  fill=Cat_potency)
+      geom_linerange(aes(ymin=lower_CI, ymax=upper_CI), color="black") +
+      geom_point(aes(y=rel_potency), alpha = 0.1) +
+      #scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
+      ggtitle("CLs of relative potencies, standard RMSEs") +
+      # scale_x_continuous(
+      #   breaks=seq(trunc(min(all_lPot$Potency_and_CI, na.rm=T)*10)/10, max(all_lPot$Potency_and_CI, na.rm=T)*1.1, by=0.4),
+      # ) +
+      theme_bw() +
+      theme(axis.text = element_text(face="bold", size=15),
+            axis.text.x = element_text(angle=90),
+            plot.title= element_text(size=15, face="bold"))
+    #P_linerangeCI
+    output$linerangeCIs <- renderPlot({ P_linerangeCI })
+    
+    P_widthCIs <- ggplot(CalcPotDF, aes(x=width_CI, fill="blue")) +
       geom_histogram(color="#e9ecef", alpha=0.6, position = "identity") +
-      scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
+      geom_density(alpha = 0.1) +
-      ggtitle("Histogram of relative potencies, standard RMSEs") +
+      #scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
-      scale_x_continuous(
+      labs(title = "Histogram of width of CLs", subtitle = paste("with upper",input$TolConf,input$TolPop, "tolerance interval")) +
-        breaks=seq(trunc(min(all_lPot$Potency_and_CI, na.rm=T)*10)/10, max(all_lPot$Potency_and_CI, na.rm=T)*1.1, by=0.4),
+      geom_vline(xintercept = widthCLTol[[5]]) +
-      ) +
       theme_bw() +
       theme(axis.text = element_text(face="bold", size=15),
             axis.text.x = element_text(angle=90),
             plot.title= element_text(size=15, face="bold")) 
     
-    output$histCIs <- renderPlot({ P_histCI })
+    output$widthCIs <- renderPlot({ P_widthCIs })
+    
+    SlopeTol <- normtol.int(x = slopeRatio, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
+    LasTol <- normtol.int(x = LasRatio, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
+    UasTol <- normtol.int(x = UasRatio, alpha = 1-input$TolConf, P = input$TolPop, side = 2)
+    
+    P_ratioSlope <- ggplot(ratioDF, aes(x=slopeRatio, fill="turquoise")) +
+      geom_histogram(color="#e9ecef", alpha=0.6, position = "identity") +
+      geom_density(alpha = 0.1) +
+      #scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
+      labs(title = "Histogram of Hill's slope ratios ", subtitle = paste("with",input$TolConf,input$TolPop, "tolerance interval")) +
+      geom_vline(xintercept = c(SlopeTol[[4]], SlopeTol[[5]]), linetype=2, col="grey") +
+      theme_bw() +
+      theme(axis.text = element_text(face="bold", size=15),
+            axis.text.x = element_text(angle=90),
+            plot.title= element_text(size=15, face="bold"))
+    
+    output$ratioSlope <- renderPlot({ P_ratioSlope })
+    
+    P_ratioLas <- ggplot(ratioDF, aes(x=LasRatio, fill="turquoise")) +
+      geom_histogram(color="#e9ecef", alpha=0.6, position = "identity") +
+      geom_density(alpha = 0.1) +
+      #scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
+      labs(title = "Histogram of lower asymptote ratios ", subtitle = paste("with",input$TolConf,input$TolPop, "tolerance interval")) +
+      geom_vline(xintercept = c(LasTol[[4]], LasTol[[5]]), linetype=2, col="grey") +
+      theme_bw() +
+      theme(axis.text = element_text(face="bold", size=15),
+            axis.text.x = element_text(angle=90),
+            plot.title= element_text(size=15, face="bold"))
+    
+    output$ratioLas <- renderPlot({ P_ratioLas })
+    
+    P_ratioUas <- ggplot(ratioDF, aes(x=UasRatio, fill="turquoise")) +
+      geom_histogram(color="#e9ecef", alpha=0.6, position = "identity") +
+      geom_density(alpha = 0.1) +
+      #scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
+      labs(title = "Histogram of upper asymptote ratios ", subtitle = paste("with",input$TolConf,input$TolPop, "tolerance interval")) +
+      geom_vline(xintercept = c(UasTol[[4]], UasTol[[5]]), linetype=2, col="grey") +
+      theme_bw() +
+      theme(axis.text = element_text(face="bold", size=15),
+            axis.text.x = element_text(angle=90),
+            plot.title= element_text(size=15, face="bold"))
+    
+    output$ratioUas <- renderPlot({ P_ratioUas })
     
     output$histEC50REF <- renderPlot({
       hist(EC50REF, col="steelblue", border="white", main = 'Histogram of EC50REF')
@@ -2350,27 +2496,28 @@ server <- function(input, output, session) {
     Dat$histLasREF <- hist(LasREF, col="violet", border="white", main = 'Histogram of EC50REF')
     Dat$histUasREF <- hist(UasREF, col="darkturquoise", border="white", main = 'Histogram of EC50REF')
    
+    #####  Dilutions Simulator ----
     tab <- AllXL[[1]]
     dils <- tab$log_dose
-    min_y <- min(tab[, 1:2])
+    min_y <- min(tab[, 1:2], na.rm = T)
-    max_y <- max(tab[, 1:2])
+    max_y <- max(tab[, 1:2], na.rm = T)
-
+#browser()
-    if (input$fixupper) {
+    # if (input$fixupper) {
-      dils_av <- dils - max(dils)
+    #   dils_av <- dils - max(dils)
-      dils_av_ <- dils_av * (input$dilslider / 100 + 1)
+    #   dils_av_ <- dils_av * (input$dilslider / 100 + 1) + input$dilmover
-      dils2 <- round(dils_av_ + max(dils), 4)
+    #   dils2 <- round(dils_av_ + max(dils), 4)
-      dilfactors <- 1 / exp(dils2 - lag(dils2))
+    #   dilfactors <- 1 / exp(dils2 - lag(dils2))
-    } else {
+    # } else {
       if (!is.null(EC50TEST)) {
         av <- mean(EC50TEST, na.rm = TRUE)
       } else {
         av <- (min(dils) + max(dils)) / 2
       }
       dils_av <- dils - av
-      dils_avsc <- dils_av * (input$dilslider / 100 + 1)
+      dils_avsc <- dils_av * (input$dilslider / 100 + 1) + input$dilmover
       dils2 <- dils_avsc + av
       dilfactors <- 1 / exp(dils2 - lag(dils2))
-    }
+    #}
  
     
    Dat$newDils <- dils2
@@ -2423,16 +2570,20 @@ server <- function(input, output, session) {
       )
       DilsTable
     })
+##### Plot for dilution slider ----
     
    if (!is.null(p2)) {
      #p2 <- Dat$p2
      p_dil <- p2 +
-       annotate("pointrange", x = dils2, y = rep(min_y, length(dils2)), xmin = min(dils2), xmax = max(dils2)) +
+       geom_vline(xintercept = dils2, col = "red", linetype = 2, alpha=0.5) +
-       annotate("text", x = dils2, y = rep(min_y + (max_y - min_y) * 0.05, length(dils2)), label = as.character(round(dils2, 3))) +
+       annotate("pointrange", x = dils2, y = rep(min_y, length(dils2)), xmin = min(dils2), xmax = max(dils2),
+                colour = "red" ,linetype = 3, shape=24) +
+       annotate("text", x = dils2, y = rep(min_y + (max_y - min_y) * 0.05, length(dils2)), label = as.character(round(dils2, 3)),colour = "red") +
        annotate("text",
          x = dils2[-1] + (max(dils2) - min(dils2)) * 0.05,
          y = rep(min_y + (max_y - min_y) * 0.1, length(dils2[-1])),
-         label = as.character(round(dilfactors[-1], 3)))
+         label = as.character(round(dilfactors[-1], 3)),colour = "red") 
+     
        # geom_line(
        #   data = as.data.frame(pl_df), aes(x = dils2, y = SAMPLE50), color = "grey15", linetype = 2,
        #   inherit.aes = F
@@ -2441,23 +2592,6 @@ server <- function(input, output, session) {
        #   data = as.data.frame(pl_df), aes(x = dils2, y = SAMPLE200), color = "grey15", linetype = 2,
        #   inherit.aes = F
        # ) +
-       # geom_vline(xintercept = c(Xbend50, Xbend200), col = "grey15", linetype = 2) +
-       # { if (input$scenario == "scenario 6") {
-       #     annotate("pointrange",
-       #       x = optdils2, y = rep(min_y + (max_y - min_y) * 0.2, length(optdils2)),
-       #       xmin = min(optdils2), xmax = max(optdils2), color = "seagreen"
-       #     )
-       #   }
-       # } +
-       # {
-       #   if (input$scenario == "scenario 6") {
-       #     annotate("text",
-       #       x = optdils2, y = rep(min_y + (max_y - min_y) * 0.25, length(optdils2)),
-       #       label = as.character(round(optdils2, 3)), color = "seagreen"
-       #     )
-       #   }
-       # } +
-      
     
        # annotate("text",
        #   x = optdils[1], y = (max_y + min_y) * 0.5,
@@ -2468,6 +2602,7 @@ server <- function(input, output, session) {
           print(p_dil)
        })
      
+       Dat$DilPlot <- p_dil
       } # if (!is.null(p2))
     } # if !is.null Dat$Mws
 
@@ -2546,7 +2681,7 @@ server <- function(input, output, session) {
     })
   })
 
-  #### simulations ----
+  #### NOT SHOWN: simulations ----
   observe({
     observeEvent(input$goSim, {
       sd_fac_ <- as.numeric(input$sdfac)
@@ -2659,7 +2794,7 @@ server <- function(input, output, session) {
   })
 
 
-  #### simulation Histograms output ----
+  #### NOT SHOWN: simulation Histograms output ----
 
   output$plotHistuAs <- renderPlot({
     if (!is.null(Dat$resHist)) {
@@ -2790,6 +2925,124 @@ server <- function(input, output, session) {
       )
     }
   )
+    
+    #### Download ROUT report ----
+    output$downloadROUTReport <- downloadHandler(
+      filename = paste0("Report_ROUT_Evaluation.pdf"),
+      content = function(file) {
+        tpdr <- tempdir()
+        tempReport <- file.path(tpdr, "BioassayReportROUT.Rmd")
+        file.copy("BioassayReportROUT.Rmd", tempReport, overwrite = T)
+        
+        tempReportc <- file.path(tpdr, "logov2.png")
+        file.copy("logov2.png", tempReportc, overwrite = T)
+        
+        rmarkdown::render(tempReport,
+                          output_file = file,
+                          params = list(
+                            FileName = Dat$FileName,
+                            author = Dat$Author,
+                            NoP = Dat$NoP,
+                            Assay = Dat$Assay,
+                            REP = REP,
+                            coeffs = Dat$coeffs_UN
+                          ),
+                          envir = new.env(parent = globalenv())
+        )
+      }
+    )
+    
+    #### download Meta 4PL report----
+    
+    observeEvent(input$btn2, {
+      if(!Dat$FITsFlag) {
+        runjs("$('#downloadXLReportMeta')[0].click();")
+      } else {
+        showModal(modalDialog(
+          title = "No 4PL model to Download",
+          "Please select other data before download.",
+          easyClose = TRUE,
+          footer = NULL
+        ))
+      }
+    })
+    
+    output$downloadXLReportMeta <- downloadHandler(
+      filename = paste0("Report_4PLEvaluation", Dat$RepIdentifier, ".pdf"),
+      content = function(file) {
+        tpdr <- tempdir()
+        tempReport <- file.path(tpdr, "Doc_BioassayReport.Rmd")
+        file.copy("Doc_BioassayReport.Rmd", tempReport, overwrite = T)
+        
+        tempReportc <- file.path(tpdr, "logov2.png")
+        file.copy("logov2.png", tempReportc, overwrite = T)
+        
+        rmarkdown::render(tempReport,
+                          output_file = file,
+                          params = list(
+                            FileName = Dat$FileName,
+                            author = Dat$Author,
+                            NoP = Dat$NoP,
+                            Assay = Dat$Assay,
+                            REP = REP,
+                            coeffs = Dat$coeffs_UN
+                          ),
+                          envir = new.env(parent = globalenv())
+        )
+      }
+    )
+    
+    
+    #### download Wizard report ----
+    
+    output$downloadWizardData <- downloadHandler(
+      
+      filename = paste0("CompiledData", Dat$nameRep, ".zip"),
+      content = function(file) {
+        fs <- c()
+        tpdr <- tempdir()
+        filename = paste0("CompiledData", Dat$nameRep, ".zip")
+        #tempReport <- file.path(tpdr, "Doc_BioassayLinReport.Rmd")
+        #file.copy("Doc_BioassayLinReport.Rmd", tempReport, overwrite = TRUE)
+        
+        #tempReportc <- file.path(tpdr, "logov2.png")
+        #file.copy("logov2.png", tempReportc, overwrite = TRUE)
+        
+        # rmarkdown::render(tempReport,
+        #                   output_file = file,
+        #                   params = list(
+        #                     FileName = Dat$FileName,
+        #                     author = Dat$Author,
+        #                     NoP = Dat$NoP,
+        #                     Assay = Dat$Assay,
+        #                     REP = REP,
+        #                     REPlin = REPlin,
+        #                     coeffsLin = Dat$coeffs_UN
+        #                   ),
+        #                   envir = new.env(parent = globalenv())
+        # )
+     #browser()   
+        fileOutModU=paste(paste0(tpdr, sep='/', 'unrModelFits'), sep='','.csv')
+        fs=c(fs, fileOutModU)
+        ModU <- Dat$ModU
+        write.csv(ModU, fileOutModU, row.names = F)
+        
+        fileOutModR=paste(paste0(tpdr, sep='/', 'restrModelFits'), sep='','.csv')
+        fs=c(fs, fileOutModR)
+        ModR <- Dat$ModR
+        write.csv(ModR, fileOutModR, row.names = F)
+        
+        DilPlot <- Dat$DilPlot
+        fileOutDilPlot =paste(paste0(tpdr, sep='/', 'SigmoidDilutionsPlot'), sep='','.png')
+        fs=c(fs, fileOutDilPlot)
+        png(fileOutDilPlot, width=600, height=400)
+        print(DilPlot)
+        dev.off()
+    #browser()    
+        zip::zipr(zipfile=file, files=fs, include_directories = F)
+        
+      }, contentType = "application/zip"
+    )
 }
 
 shinyApp(ui, server)