lexibel wenn fehlende Werte

R/Global.R

@@ -711,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, ]
@@ -798,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) {
@@ -1201,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)
@@ -1239,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)
@@ -1276,7 +1278,7 @@ 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)
   DatL$RMSE_r <- RMSE_r
@@ -1428,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)
@@ -1436,7 +1441,9 @@ 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]]
   smu <- FITs[[2]]

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/app.R

@@ -218,6 +218,8 @@ server <- function(input, output, session) {
             "Uploaded data",
             tableOutput("XLdata")
           ),
+          
+          ###### 4pl output ----
           tabPanel(
             "4pl-Analysis",
             tags$style(HTML("pre { color: black; background-color: #FFE1FF;
@@ -265,6 +267,7 @@ server <- function(input, output, session) {
               )
             )
           ),
+          ##### linear output ----
           tabPanel(
             "linear Analysis",
             sidebarLayout(
@@ -324,10 +327,13 @@ server <- function(input, output, session) {
               )
             )
           ),
+          ##### ROUT output ----
           tabPanel("robust outlier testing",
-                   sliderInput("Qslider", "adjust Q-value in %",min=0.1, max = 20, value = 1, step=0.1),
+                   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",
@@ -816,7 +822,7 @@ server <- function(input, output, session) {
           ##### ROUT outlier testing ----
     #browser()
           if(!is.null(input$Qslider)) {
-            all_lROUT <- all_l
+            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,]
@@ -847,12 +853,42 @@ server <- function(input, output, session) {
               })
             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
+            
             }          
 
   
@@ -1842,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), ]
@@ -2889,6 +2926,31 @@ 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----