diff --git a/.DS_Store b/.DS_Store index 4fe77ec..0c391a1 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/MultiL2.numbers b/MultiL2.numbers new file mode 100644 index 0000000..e219dbc Binary files /dev/null and b/MultiL2.numbers differ diff --git a/MultiL_del.xlsx b/MultiL_del.xlsx new file mode 100644 index 0000000..a1385e7 Binary files /dev/null and b/MultiL_del.xlsx differ diff --git a/R/Global.R b/R/Global.R index 9b6128a..4bf9614 100644 --- a/R/Global.R +++ b/R/Global.R @@ -23,6 +23,25 @@ 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),] + return(cor(df2[,1],df2[,2])) +} #' Levenberg Marquard fit of 4 pl #' @@ -46,8 +65,9 @@ library(scales) #' Dat <- list() #' te <- Fitting_FUNC(dat, TransF) #' print(te) -Fitting_FUNC <- function(ro_new, TransFlag = FALSE) { - CORro <- cor(ro_new[, 1], ro_new[, ncol(ro_new)]) +Fitting_FUNC <- function(ro_new, TransFlag = FALSE, nameWS ="") { + #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 +89,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) { }, warning = function(e) { mr <<- "In nlsModel singular gradient matrix" + } ) # Stop if singular gradient matrix @@ -82,6 +103,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 +125,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 +147,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 +185,27 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) { } ) } + +#browser() if (!TransFlag) { - pot_est <- exp(confintd(mr, "r", method = "asymptotic")) - potU_est <- exp(confintd(mu, "r", method = "asymptotic")) - PRED <- predict(mr) - PREDu <- predict(mu) +#browser() + if (length(s_mr) ==1 | length(Sum_u) ==1) { + return("failed") + } 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")) - potU_est <- exp(confintd(muT, "r", method = "asymptotic")) - PRED <- predict(mrT) - PREDu <- predict(muT) + if (length(s_mr) ==1 | length(Sum_u) ==1) { + return("failed") + }else { + 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 +256,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 +316,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) @@ -1018,7 +1062,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 +1156,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 - upper_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_xs - log_xt + qt(Conf, DFs) * se_log_ratio ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) return(ci_ratio) } @@ -1143,6 +1187,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) @@ -1228,9 +1275,9 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { RSS_r <- round(sum(smr$residuals^2), 5) MSE_r <- RSS_r / (nrow(all_l) - 5) RMSE_r <- round(sqrt(MSE_r), 6) - Dat$RMSE_r <- RMSE_r - Dat$RMSE_pure <- RMSE_pure - Dat$RMSE_unr <- round(RMSEunr, 6) + DatL$RMSE_r <- RMSE_r + DatL$RMSE_pure <- RMSE_pure + DatL$RMSE_unr <- round(RMSEunr, 6) coeffs <- smu$coefficients[, 1] # browser() @@ -1242,6 +1289,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,12 +1302,13 @@ 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) - - Dat$uAsCI <- uAsCI2 + estUppA <- round(dt / ds, 5) + DatL$uAsCI <- uAsCI2 +# browser() + #### EQ test on slope ratio ---- # bs <- coeffs["bs"] # bt <- coeffs["bt"] @@ -1271,11 +1320,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 +1336,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 +1354,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 +1386,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, - round(dt_at / ds_as, 5), round(potAll2[1] * 100, 2), round(potAllU2[1] * 100, 2) + estLowA, estSlope, estUppA, p_F_nonlin, + 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]]), diff --git a/TestF4PLGHZ2.xlsx b/T4PLGHZ200_1MAK_4P.xlsx similarity index 100% rename from TestF4PLGHZ2.xlsx rename to T4PLGHZ200_1MAK_4P.xlsx diff --git a/TestF4PLGHZ1.xlsx b/T4PLGHZ50_1MAK_4P.xlsx similarity index 100% rename from TestF4PLGHZ1.xlsx rename to T4PLGHZ50_1MAK_4P.xlsx diff --git a/Tests3Plates.numbers b/Tests3Plates.numbers index 75c69d9..adf4ae2 100644 Binary files a/Tests3Plates.numbers and b/Tests3Plates.numbers differ diff --git a/Tests4Plates_fail.numbers b/Tests4Plates_fail.numbers new file mode 100644 index 0000000..c8d07f5 Binary files /dev/null and b/Tests4Plates_fail.numbers differ diff --git a/Tests4Plates_fail.xlsx b/Tests4Plates_fail.xlsx new file mode 100644 index 0000000..eae7fb0 Binary files /dev/null and b/Tests4Plates_fail.xlsx differ diff --git a/dev/app.R b/dev/app.R index 3d99621..2ad4cb1 100644 --- a/dev/app.R +++ b/dev/app.R @@ -28,6 +28,7 @@ library(twopartm) library(car) library(dplyr) library(scales) +library(tolerance) source("../R/Global.R") @@ -116,7 +117,7 @@ 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.", ), column(6, ) ), @@ -159,10 +160,11 @@ server <- function(input, output, session) { ), # 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?", @@ -206,6 +208,11 @@ server <- function(input, output, session) { numericInput("uEACdiffla", "upper EAC for diff. of LA", 0.189, step = 0.001) ) ), + + tabPanel( + "Uploaded data", + tableOutput("XLdata") + ), tabPanel( "4pl-Analysis", tags$style(HTML("pre { color: black; background-color: #FFE1FF; @@ -376,7 +383,7 @@ server <- function(input, output, session) { mainPanel( width = 12, tabsetPanel( - id = "tabs", + id = "tabs2", tabPanel( "Settings", h4("Settings of 4PL regression"), @@ -588,12 +595,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", @@ -736,17 +746,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", "") @@ -775,7 +787,7 @@ server <- function(input, output, session) { #### 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 +807,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." @@ -2130,55 +2142,56 @@ 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) - - 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_ - - 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))) - #browser() - dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, 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) - + 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_ + + 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))) + #browser() + dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, 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) + } } #for N_WS #browser() @@ -2211,11 +2224,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)) + @@ -2229,70 +2246,19 @@ server <- function(input, output, session) { 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)") + @@ -2354,23 +2320,23 @@ server <- function(input, output, session) { dils <- tab$log_dose min_y <- min(tab[, 1:2]) max_y <- max(tab[, 1:2]) - - if (input$fixupper) { - dils_av <- dils - max(dils) - dils_av_ <- dils_av * (input$dilslider / 100 + 1) - dils2 <- round(dils_av_ + max(dils), 4) - dilfactors <- 1 / exp(dils2 - lag(dils2)) - } else { +#browser() + # if (input$fixupper) { + # dils_av <- dils - max(dils) + # dils_av_ <- dils_av * (input$dilslider / 100 + 1) + input$dilmover + # dils2 <- round(dils_av_ + max(dils), 4) + # dilfactors <- 1 / exp(dils2 - lag(dils2)) + # } 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 diff --git a/tests/EQ tests GHZ1_2.pages b/tests/EQ tests GHZ1_2.pages index df5dcbe..c6e9d2e 100644 Binary files a/tests/EQ tests GHZ1_2.pages and b/tests/EQ tests GHZ1_2.pages differ