Dilution simulator updated

2026-06-05 11:32:15 +02:00
parent f5575f8429
commit abf02ef4ec
3 changed files with 97 additions and 40 deletions
@@ -177,7 +177,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
 #' @export
 #' @examples
 #' suppressMessages(source("../../dev/setup.R"))
-#' data.frame(
+#' dat <- data.frame(
 #'   R_dil1 = c(
 #'     10.0651024695491, 10.9844983291817, 10.7635586089293, 10.4597656321327, 10.3898668457823, 10.8171761349909,
 #'     10.319758021908, 10.1304854046653
@@ -265,11 +265,11 @@ plotSingularity <- function(dat) { # sigmoid,det_sig,
 #'   TEST1 = c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2 = c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
 #'   TEST3 = c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose = c(5.01, 3.401, 2.708, 2.015, 1.32176, 0.62861, -0.0645385, -1.6739764)
 #' )
-#' sigmoid <- c(0.7163324, 0.5636804, 10.6156340, 9.9784160, -0.7504673, -0.7108692, -3.5788141, -0.6662962)
-#' det_sig <- FALSE
-#' TransF <- FALSE
+#' 
+#' 
+#' TransFlag <- FALSE
 #' Dat <- list()
-#' p <- plot_f(dat, sigmoid, det_sig, TransF)
+#' p <- plot_f(dat, TransFlag)
 #' print(p)
 plot_f <- function(dat, TransFlag = FALSE) { # sigmoid,det_sig,
  CORdat <- cor(dat[, 1], dat[, ncol(dat)])
@@ -490,7 +490,7 @@ plot_f <- function(dat, TransFlag = FALSE) { # sigmoid,det_sig,
 #' @param heteroNoise Boolean if heteroscedstic noise should be added.
 #' @param noDilSeries A number >1 indicating how many dilution series per product should be simulated.
 #' @param noDils A number >7 indicating how many dilutions steps per product should be simulated.
-#' @returns A data-frame with readouts and natural log of concentrations.
+#' @returns A data-frame with readouts and natural log of concentrations. avs: means per concentration; sds: standard deviation per concentration;cvs: coefficients of variation
 #' @export
 #' @examples
 #' suppressMessages(source("../../dev/setup.R"))
@@ -656,7 +656,8 @@ LinPotTab <- function(circles, Lim, PureErrFlag) {
 #' PureErrF <- TRUE
 #'
 #'
-#' ANOVAlintests(ro_new, circles, Lim, PureErrF)
+#' ANOVAlintests(dat, circles, Lim, PureErrF)
+#' 
 ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
@@ -873,6 +874,7 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
 #' )
 #'
 #' PlotLinPLA_FUNC(circle, sigmoid, all_l2, pl_df, indS, indT)
+#' 
 PlotLinPLA_FUNC <- function(circle, sigmoid, all_l2, pl_df, indS, indT) {
  # browser()
  mLin <- gsl_nls(readout ~ (intS + r) * isSample + intS * isRef + k * log_dose,
@@ -1100,7 +1102,7 @@ pot4plFUNC <- function(ro_new, PureErrFlag) {
 #' se_xt <- 0.34
 #' se_xs <- 0.23
 #' DFs <- 32 - 16
-#'
+#' CoVar <- 0
 #'
 #' ParamCI_F(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975)
 ParamCI_F <- function(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975) {
@@ -1371,6 +1373,7 @@ 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")
  all_len <- nrow(all_l)
@@ -579,36 +579,32 @@ server <- function(input, output, session) {
          sidebarPanel(
            width = 3,
            fluidRow(
-              column(
-                6,
+              
                box(
                  title = "Upload multiple worksheets", status = "warning", solidHeader = TRUE, width = 12, "Please upload your EXCEL file here",
                  fileInput("MiFile", "", accept = ".xlsx")
              )
-
-              )
            )
          ),
          mainPanel(
            tabsetPanel(
              id = "tabs",
-              tabPanel(
-                "4pl",
+              tabPanel("4pl",
                box(
                  title = "ANOVA table", status = "primary", solidHeader = TRUE, width = 12,
                  tableOutput("Anovatab")
                ),

-                column(8,
+                column(12,
                  # h3("Confidence intervals"),
                  # tableOutput("CIs"),
-                  # "The confidence interval table is interaactive for changes in: variability slider (%SD), potency of test-slider,
+                  # "The confidence interval table is interactive for changes in: variability slider (%SD), potency of test-slider,
                  #                      and 'Adjust the dilutions'-slider",
                  # tableOutput("optimalDils"),
                  plotOutput("sigPlotREF"),
-                  selectInput(inputId = "scenario", label = "Select an 'optimal' scenario:", choices = c("scenario 2", "scenario 3", "scenario 6", "steep slope"))
-                ),
-                column(5,
+                  plotOutput("sigPlotTEST"),
+                  #selectInput(inputId = "scenario", label = "Select an 'optimal' scenario:", choices = c("scenario 2", "scenario 3", "scenario 6", "steep slope"))
+                
                   plotOutput("plotfordilutions"),
                  # h4("in grey: most extreme bend point lines of theoretical samples with 50% and 200% potency"),
                  # sliderInput("dilslider", "Adjust the dilutions(+-change in %)", min = -100, max = 100, value = 0, step = 1, round = 0),
@@ -627,6 +623,14 @@ server <- function(input, output, session) {
                  h4("Explanation of the plots")
                )
              ),
+              tabPanel("Histograms",
+                h4("Histograms of parameters"),
+                plotOutput("histCIs"),
+                plotOutput("histEC50REF"), 
+                plotOutput("histLasREF"),
+                plotOutput("histUasREF"),
+                
+              ),
              tabPanel(
                "Report",
                h4("Settings for report")
@@ -766,7 +770,7 @@ server <- function(input, output, session) {
          REP$all_l <- all_l

          #### XLSX eval ----
-          if (CORro < 0) SLOPE <- -1 else SLOPE <- 1
+          #if (CORro < 0) SLOPE <- -1 else SLOPE <- 1
          FITs <- Fitting_FUNC(XLdat2, TransFlag = FALSE)

          #### if no 4pl fit is possible ----
@@ -2117,7 +2121,7 @@ server <- function(input, output, session) {
    AllXL <- Dat$Mws
    AllSheets <- Dat$Msheets
    
-    URMcoefsL <- list()
+    URMcoefsL <- RMcoefsL <- potEstL <- list()
    
    for (N_WS in 1:length(AllXL)) {
   
@@ -2144,6 +2148,7 @@ server <- function(input, output, session) {
    FITs <- Fitting_FUNC(datWS2, TransFlag = F)
 
    pot_est <- FITs[[3]]
+    potEstL[[N_WS]] <- pot_est
    potU_est <- FITs[[4]]
    # unrestricted
    SU_mu <- FITs[[2]]
@@ -2152,6 +2157,14 @@ server <- function(input, output, session) {
    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)
    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)))
@@ -2174,10 +2187,16 @@ 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]
-    # 
-    # Dat$URMcoefsDF <- URMcoefsDF
-    # Dat$RestrM <- RestrM
-    # Dat$CalcPot <- CalcPot
+    UasTEST <- as.numeric(URMcoefsDF[,4])
+    LasTEST <- as.numeric(URMcoefsDF[,2])
+    
+    RMcoefsDF <- t(matrix(unlist(RMcoefsL),nrow=6))
+    
+    Dat$URMcoefsDF <- URMcoefsDF
+    Dat$RestrM <- RMcoefsDF
+    
+    CalcPotDF <-  t(matrix(unlist(potEstL),nrow=3))
+    Dat$CalcPot <- CalcPotDF
    # 
    #### sigmoid plots ----
    
@@ -2199,23 +2218,58 @@ server <- function(input, output, session) {

    output$sigPlotREF <- renderPlot({ p1 })

-    # PLOTS$sigPlotREF <- p1
+    Dat$sigPlotREF <- p1
    # 
-    # p2 <- ggplot(SIGtestDF, aes(x_X, y=sigTest, col=as.factor(Prod))) +
-    #   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_vline(xintercept = EC50TEST, alpha = 0.2) +
-    #   xlab("dilutions") +
-    #   ggtitle("Plot of all calculated reference fits (unrestricted model, in gray vertical lines: EC50)") +
-    #   theme_bw() +
-    #   theme(axis.text = element_text(face = "bold", size = 15),
-    #         plot.title = element_text(size = 15, face = "bold"))
-    # 
-    # output$sigPlotREF <- renderPlot({ p2 })
-    # 
-    # PLOTS$sigPlotTEST <- p2
+    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) +
+      annotate("text", label="o", x=x_LA, y=LasTEST, alpha=0.2) +
+      geom_vline(xintercept = EC50TEST, alpha = 0.2) +
+      xlab("dilutions") +
+      ggtitle("Calculated test sample fits (unrestricted model, in gray vertical lines: EC50)") +
+      theme_bw() +
+      theme(axis.text = element_text(face = "bold", size = 15),
+            plot.title = element_text(size = 15, face = "bold"))

+    output$sigPlotTEST <- renderPlot({ p2 })
+
+    Dat$sigPlotTEST <- p2
+    
+    #### 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))),
+      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)) +
+      geom_histogram(color="#e9ecef", alpha=0.6, position = "identity") +
+      scale_fill_manual(values=c("darkgreen","darkblue","salmon2","tomato3")) +
+      ggtitle("Histogram 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"))
+    
+    output$histCIs <- renderPlot({ P_histCI })
+    
+    output$histEC50REF <- renderPlot({
+      hist(EC50REF, col="steelblue", border="white", main = 'Histogram of EC50REF')
+    })
+    output$histLasREF <- renderPlot({
+      hist(LasREF, col="violet", border="white",main = 'Histogram of lower asymptotes REF')
+    })
+    output$histUasREF <- renderPlot({
+      hist(UasREF, col="darkturquoise", border="white",main = 'Histogram of upper asymptotes REF')
+    })
+    Dat$histEC50REF <- hist(EC50REF, col="steelblue", border="white", main = 'Histogram of EC50REF')
+    Dat$histLasREF <- hist(LasREF, col="violet", border="white", main = 'Histogram of EC50REF')
+    Dat$histUasREF <- hist(UasREF, col="darkturquoise", border="white", main = 'Histogram of EC50REF')
    
    # dils <- tab$log_dose
    # min_y <- min(tab[, 1:3])