library(shiny) library(shinydashboard) library(shinyjs) library(shinyAce) library(shinydashboard) library(shinycssloaders) library(shinyBS) library(purrr) library(gslnls) library(tidyverse) library(ggplot2) library(reshape2) library(openxlsx) library(DT) library(ggpubr) library(gridExtra) library(drc) library(twopartm) library(car) library(dplyr) #### reactive values ---- Dat <- reactiveValues() REP <- reactiveValues() #### functions ---- # dilFUN2 <- function(cs_,dils,Faktor) { # av <- cs_ # dils_av <- dils # dils_avsc <- dils_av*Faktor # dils2 <- dils_avsc+av # dilfactors <- 1/exp(dils2-lag(dils2)) # return(dilfactors) # } #' Levenberg Marquard fit of 4 pl #' #' Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, #' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. #' #' @param ro_new The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. #' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). #' @returns Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, #' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. #' @export #' @examples #' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), #' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), #' 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)) #' TransF<- FALSE #' Dat<- list() #' te <- Fitting_FUNC(dat,TransF) #' print(te) Fitting_FUNC <- function(ro_new, TransFlag=F) { CORro <- cor(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) isSample <- rep(c(0,1),1,each=nrow(all_l)/2) all_l2 <- cbind(all_l, isRef, isSample) #browser() if (CORro<0) SLOPE <- -1 else SLOPE <- 1 if (!TransFlag) { startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0) mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), data=all_l2, start=startlist,#trace=T, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) s_mr <- tryCatch({ s_mr <- summary(mr) }, error = function(err) { s_mr <- NULL }) } else { startlist <- list(a=log(min(ro_new[,2])), b=SLOPE, d=log(max(ro_new[,2])), cs=mean(all_l$log_dose),r=0) mrT <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), data=all_l2, start=startlist, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) s_mr <- summary(mrT) } 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) tryCatch({ mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), data=all_l, start=startlistmu, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(msg){ return(0) }) Sum_u <- tryCatch({ summary(mu) }, error=function(msg){ 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) tryCatch({ muT <- gsl_nls(fn = log(readout) ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), data=all_l, start=startlistmu, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(msg){ return(0) }) Sum_u <- tryCatch({ summary(muT) }, error=function(msg){ return(0) }) } if (!TransFlag) { 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) } return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu)) } #' Plot sigmoidal curve #' #' Returns the final plots of the 4pl function as sigmoidal lines, and the single readouts as scatter, with REF in blue and TEST in red. #' Two plots are returned in a grid object: the unrestricted model and the restricted model. #' #' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. #' @param sigmoid The parameters of the 4pl fit of unrestricted model from EXCEL input. #' @param det_sig The parameters of the 4pl fit of unrestricted model from the metadata input. #' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). #' @returns A grid object either of the original scale or the natural log of the readouts. #' @export #' @examples #' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), #' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), #' 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 #' Dat<- list() #' p <- plot_f(dat,sigmoid,det_sig,TransF) #' print(p) plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig, CORdat <- cor(dat[,1],dat[,ncol(dat)]) #browser() all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) isSample <- rep(c(0,1),1,each=nrow(all_l)/2) all_l2 <- cbind(all_l, isRef, isSample) #browser() MODLS <- Fitting_FUNC(dat, TransFlag=F) s_mr <- MODLS[[1]] a <- s_mr$coefficients["a",1] b <- s_mr$coefficients["b",1] cs <- s_mr$coefficients["cs",1] d <- s_mr$coefficients["d",1] r <- s_mr$coefficients["r",1] log_dose <- unique(all_l$log_dose) seq_x <- seq(min(log_dose),max(log_dose),0.1) SAMPLE <- a+(d-a)/(1+exp(b*((cs-r)-seq_x))) REF <- a+(d-a)/(1+exp(b*((cs)-seq_x))) s_mu <- MODLS[[2]] as <- s_mu$coefficients["as",1] bs <- s_mu$coefficients["bs",1] cs <- s_mu$coefficients["cs",1] ds <- s_mu$coefficients["ds",1] at <- s_mu$coefficients["at",1] bt <- s_mu$coefficients["bt",1] ct <- s_mu$coefficients["cs",1]-s_mu$coefficients["r",1] dt <- s_mu$coefficients["dt",1] r <- s_mu$coefficients["r",1] SAMPLEtrue <- at + (dt -at)/(1+exp(bt*((cs-r-seq_x)))) REFtrue <- as + (ds -as)/(1+exp(bs*((cs-seq_x)))) #browser() pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue) all_l2$readout[all_l2$readout < 0] <- 0.01 all_l2$readouttrans <- log(all_l2$readout) slopeEC50 <- b*(d-a)/4 Intercept <- a+(d-a)/2-b*(d-a)/4*cs #browser() Xbendl3 <- cs-(1.5434/b) Xbendu3 <- cs+(1.5434/b) XbendlT <- cs-r-(1.5434/b) XbenduT <- cs-r+(1.5434/b) XasymplS <- cs-(3/b) XasympuS <- cs+(3/b) XasymplT <- cs-r-(3/b) XasympuT <- cs-r+(3/b) bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3), bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,3), asympREF_lower = round(XasymplS,3), asympREF_upper=round(XasympuS,3), asympSAMPLE_lower = round(XasymplT,3), asympSAMPLE_upper=round(XasympuT,3)) Dat$bendpoints <- bendpoints Dat$cfordils <- cs p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + geom_point(shape=factor(isRef), size=3,alpha=0.8) + labs(title = paste("restricted 4pl"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + scale_shape_manual(labels=c("test","reference")) + theme_bw() + theme(axis.text = element_text(size=14)) p2 <- p + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLE), color="#C2173F", inherit.aes = F) + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="#4545BA", inherit.aes = F) + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="#C2173F", linetype=2, alpha=0.4, inherit.aes = F) + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="#4545BA", linetype=2, alpha=0.4, inherit.aes = F) + geom_vline(xintercept=c(Xbendl3, Xbendu3), col="#4545BA",linetype=2) + geom_vline(xintercept=c(XbendlT, XbenduT), col="#C2173F",linetype=2) + geom_vline(xintercept=c( XasymplS, XasympuS), col="#4545BA55",linetype=2) + geom_vline(xintercept=c(XasymplT, XasympuT), col="#C2173F55",linetype=2) + annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) + geom_abline(slope = slopeEC50, intercept = Intercept) + theme(legend.position="none") Dat$p2 <- p2 # transformed plots p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + geom_point(shape=factor(isRef), size=3,alpha=0.8) + labs(title = paste("restricted transformed 4pl"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + theme_bw() MODLStrans <- Fitting_FUNC(dat,TransFlag = TRUE) s_mrt <- MODLStrans[[1]] a_trans <- s_mrt$coefficients["a",1] b_trans <- s_mrt$coefficients["b",1] cs_trans <- s_mrt$coefficients["cs",1] d_trans <- s_mrt$coefficients["d",1] r_trans <- s_mrt$coefficients["r",1] XbendlTrans <- cs_trans-(1.5434/b_trans) XbenduTrans <- cs_trans+(1.5434/b_trans) XbendlTransT <- cs_trans-r_trans-(1.5434/b_trans) XbenduTransT <- cs_trans-r_trans+(1.5434/b_trans) bendpointsTRANS <- c(bendREF_lower = round(XbendlTrans,3), bendREF_upper=round(XbenduTrans,3), bendSAMPLE_lower = round(XbendlTransT,3), bendSAMPLE_upper=round(XbenduTransT,3)) Dat$bendpointsTRANS <- bendpointsTRANS SAMPLEtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans-r_trans)-seq_x))) REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans)-seq_x))) pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans) p_rt2 <- p_rt + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=SAMPLEtrans), color="#C2173F", inherit.aes = F) + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="#4545BA", inherit.aes = F) + geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="#4545BA",linetype=2) + geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="#C2173F",linetype=2) + theme(legend.position = "none", axis.text=element_text(size=14)) #browser() Sum_u <- MODLS[[2]] #browser() #if (is.null(det_sig)) { ast <- Sum_u$coefficients["as",1] ate <- Sum_u$coefficients["at",1] bst <- Sum_u$coefficients["bs",1] bte <- Sum_u$coefficients["bt",1] cst <- Sum_u$coefficients["cs",1] cte <- Sum_u$coefficients["cs",1]-Sum_u$coefficients["r",1] dst <- Sum_u$coefficients["ds",1] dte <- Sum_u$coefficients["dt",1] # } else { # ast <- det_sig[5] # ate <- det_sig[6] # bst <- det_sig[1] # bte <- det_sig[2] # cst <- det_sig[7] # cte <- det_sig[8] # dst <- det_sig[3] # dte <- det_sig[4] # } REFu <- ast + (dst-ast)/(1+exp(bst*(cst-seq_x))) SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(cte-seq_x))) pl_df2 <- cbind(seq_x, SAMPLEu, REFu) #browser() pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + geom_point(size=3) + labs(title="unrestricted 4_pl-Model", color="product") + scale_color_manual(labels = c("test","reference"), values=c("#C2173F88","#4545BA88")) + theme_bw() pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu), color="#C2173F", inherit.aes = F) + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu), color="#4545BA", inherit.aes = F, show.legend = F) pu2_ <- pu2 + theme(legend.position = "none", axis.text = element_text(size=14)) putrans <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + geom_point( size=3) + labs(title="unrestricted transformed 4_pl-Model", color="product") + scale_color_manual(labels = c("test","reference"), values=c("#C2173FAA","#4545BAAA")) + theme_bw() Sum_ut <- MODLStrans[[2]] ast_t <- Sum_ut$coefficients["as",1] ate_t <- Sum_ut$coefficients["at",1] bst_t <- Sum_ut$coefficients["bs",1] bte_t <- Sum_ut$coefficients["bt",1] cst_t <- Sum_ut$coefficients["cs",1] cte_t <- Sum_ut$coefficients["cs",1]-Sum_ut$coefficients["r",1] dst_t <- Sum_ut$coefficients["ds",1] dte_t <- Sum_ut$coefficients["dt",1] REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(cst_t-seq_x))) SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(cte_t-seq_x))) pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans) pu2_t <- putrans + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=SAMPLEu_trans), color="#C2173F", inherit.aes = F) + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans), color="#4545BA", inherit.aes = F, show.legend = F) pu3_t <- pu2_t if (TransFlag) grid.arrange(p_rt2,pu3_t, nrow=1) else grid.arrange(p2,pu2_, nrow=1) } #' Simulate a well-plate readout #' #' Returns a possible well-plate result based on the metadata provided and a variability. #' Homo- and heteroscedasticity can be simulated. #' #' @param as, bs, cs, ds, at, bt, dt,r,ct, gt, gs are the parameter of the 5pl logistic regression. #' @param sd_fac The standard deviation with units in % of upper-lower asymptote. #' @param log_conc The natural log of the concentrations. #' @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. #' @export #' @examples #' as=3; bs=1; cs=-4; ds=10; at=3; bt=1; dt=10;r=0.0001;ct=cs-r;sd_fac=0.1; gt=1; gs=1; #' lnConc=c(1,0,-1,-2,-3,-4,-5,-6) #' heteroNoise <- FALSE #' noDilS <- 3 #' noD <- 8 #' Calc_DilRes(as,bs,cs,ds,at,bt,dt,r,ct,sd_fac,gt,gs,log_conc=lnConc, heteroNoise, noDilS, noD) Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct=cs-r, sd_fac=0.1, gt=1, gs=1, log_conc, heteroNoise=FALSE, noDilSeries, noDils) { yAxfac <- (ds-as) log_dose <- log_conc isRef <- rep(c(1,0),1,each=length(log_conc)*noDilSeries) isSample <- rep(c(0,1),1,each=length(log_conc)*noDilSeries) #browser() av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose))) if (heteroNoise) { # heterosc noise ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac/100))), nrow=noDils, ncol=noDilSeries*2) } else { # homosc noise ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac*yAxfac/100))), nrow=noDils, ncol=noDilSeries*2) } ro_jit <- abs(ro_jit) ro_jit2 <- cbind(ro_jit, log_dose) if (noDilSeries==3) { colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose") } else { colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose") } return(ro_jit2) } #' Calculate the potency of linear PLA #' #' The delta Method is used for calculating the potency using a restricted linear model. #' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. #' #' @param circles Data frame with the 3 consecutive concentrations and the respective readouts. #' @param Lim The limits to check if the relative CI is within the limits. #' @param PureErrFlag Boolean if Pure Error should be used. #' @returns A data-frame with potency estimate, absolute CIs, test result, relative CIs. #' @export #' @examples #' CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9, #' -3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7), #' replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3", #' "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"), #' readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2), #' isRef=c(rep(1,9), rep(0,9)), #' isSample = c(rep(0,9), rep(1,9))) #' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant #' PureErrF <- TRUE #' #' #' LinPotTab(circles=CIRC, Lim, PureErrF) LinPotTab <- function(circles, Lim, PureErrFlag) { circ_ABl <- circles circ_Al <- circ_ABl[circ_ABl$isSample ==1,] circ_Al <- circ_ABl[circ_ABl$isSample ==0,] # restr CSSI model modAB <- lm(readout ~ log_dose + isSample, circ_ABl) coeffs <- modAB$coefficients SU_modAB <- tryCatch({ SU_modAB <- summary(modAB) }, error = function(msg) { return(NA) }) # Intercept diff/slope modAB linPot <- exp(modAB$coefficients[3]/modAB$coefficients[2]) if(PureErrFlag) { FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl)) meanPureErr <- FitAnova[4,3] DFsPure <- FitAnova[4,1] VCOV <- vcov(modAB) V_V <- VCOV/SU_modAB$sigma^2 VCOVpure <- V_V*meanPureErr SEsPure <- sqrt(diag(V_V)*meanPureErr) } log_pot_delta <- deltaMethod(modAB, "isSample/log_dose") if (PureErrFlag) { V_ <- log_pot_delta$SE^2/SU_modAB$sigma^2 V_p <- V_*meanPureErr potDeltaPureSE <- sqrt(V_p) CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure)*potDeltaPureSE CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure)*potDeltaPureSE } else { CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE } #browser() ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up)) if (ExpLinPot[2]*100>Lim[[9]] & ExpLinPot[3]*100Lim[9] & ExpLinPot[3]*100>Lim[10]) test_potCI <- 0 else test_potCI <- 1 su_mod <- summary(modAB)$coefficients su_mod2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.")), su_mod) su_modU <- summary(modABu)$coefficients su_modU2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.","slope difference")), su_modU) uCI_SloDiff <- su_modU[4,1] + qt(0.975,8)*su_modU[4,2] lCI_SloDiff <- su_modU[4,1] - qt(0.975,8)*su_modU[4,2] SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiff,uCI_SloDiff) lenCirc <- nrow(circ_ABl) dfTreat <- 3 dfPrep <- 1 dfReg <- 1 dfnonP <- 1 dfRMSE <- c(lenCirc-3-1) dfTotal <- lenCirc-1 dfPureE <- lenCirc-6 dfNonLin <- dfRMSE-dfPureE RSS <- sum(resid(lm(readout ~ log_dose*isSample, circ_ABl))^2) MSE <- RSS/dfRMSE SSE <- sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2) MSpure <- SSE/dfPureE if (PureErrFlag) { FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl)) meanPureErr <- FitAnova[4,3] SU_modAB <- tryCatch({ SU_modAB <- summary(modAB) }, error = function(msg) { return(NA) }) if (length(SU_modAB)>1) s_modABcoeffs <- summary(modAB)$coefficients DFsPure <- FitAnova[4,1] VCOV <- vcov(modAB) V_V <- VCOV/SU_modAB$sigma^2 VCOVpure <- V_V*meanPureErr SEsPure <- sqrt(diag(V_V)*meanPureErr) su_mod2[,3] <- SEsPure su_mod2[,4] <- su_mod2[,2]/su_mod2[,3] su_mod2[,5] <- 2*(1-pt(abs(su_mod[,4]), FitAnova[4,1])) s_mu <- summary(modABu)$coefficients SU_modABu <- summary(modABu) VCOVu <- vcov(modABu) V_Vu <- VCOVu/SU_modABu$sigma^2 SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) su_modU2[,3] <- SEsPureU su_modU2[,4] <- su_modU2[,2]/su_modU2[,3] su_modU2[,5] <- 2*(1-pt(abs(su_modU2[,4]), FitAnova[4,1])) uCI_SloDiffP <- su_modU[4,1] + qt(0.975,8)*SEsPureU[4] lCI_SloDiffP <- su_modU[4,1] - qt(0.975,8)*SEsPureU[4] SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiffP,uCI_SloDiffP) SSRes <- SSE dfRes <- dfPureE } else { SSRes <- RSS dfRes <- dfRMSE } # treatment SStreat <- print(sum((predict(lm(readout ~ factor(log_dose)*isSample, circ_ABl))-mean(circ_ABl$readout))^2)) F_treat <- (SStreat/dfTreat)/(SSRes/dfRes) # Preparation SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl))-mean(circ_ABl$readout))^2)) F_prep <- (SSprep/dfTreat)/(SSRes/dfRes) # Regression # ANOVA tape II SS of regression SSreg <- Anova(lm(readout ~log_dose + isSample, circ_ABl))[1,1] # Non-parallelism # diff of RSS of restricted and unrestricted model SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2) F_nonpar <- SSnonpar/(sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)/(lenCirc-4)) # non-linearity SSnonlin <- sum((predict(modABu)-predict(lm(readout ~ as.factor(log_dose)*isSample, circ_ABl)))^2) # = RSS-SSE # Total SS SStot <- sum((circ_ABl$readout-mean(circ_ABl$readout))^2) # Significance of R^2 F-ratio # MSR/MSE # sample A F_R2_A <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2)/ (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - circ_Al$readout)^2)/(nrow(circ_Al)-3)) pFR2_A <- round(pf(F_R2_A,1,6),4) # sample B F_R2_B <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2)/ (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2)/(nrow(circ_Bl)-3)) pFR2_B <- round(pf(F_R2_B,1,6),4) # sign of non-lin with pure error: MSSnonlin/MSSE F_nonlin <- (SSnonlin/2)/(SSE/dfPureE) # sign of slope F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2)/(sum((circ_Bl$readout - predict(modB))^2)/(nrow(circ_Bl)-2)) F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2)/(sum((circ_Al$readout - predict(modA))^2)/(nrow(circ_Al)-2)) # F-test on regression: MSSreg/MSSE if (is.na(F_nonlin)) F_nonlin <- 0 if (F_nonlin > 0) { p_F_nonlin <- round(pf(F_nonlin,2,dfPureE, lower.tail = F),5) } else { p_F_nonlin <- "SSnonlin neg or 0"; } # significances F_regr <- (SSreg/1)/(SSRes/dfRes) p_F_regr <- round(pf(F_regr,1,dfRes, lower.tail = F),3) p_F_treat <- round(pf(F_treat,3,dfRes, lower.tail = F),3) p_F_prep <- round(pf(F_prep,1,dfRes, lower.tail = F),3) p_F_slope_A <- round(pf(F_slope_A,1,(nrow(circ_Al)-2), lower.tail = F),3) p_F_slope_B <- round(pf(F_slope_B,1,(nrow(circ_Bl)-2), lower.tail = F),3) p_F_nonp <- round(pf(F_nonpar,1,dfRes, lower.tail = F),3) p_F_LoF <- p_F_nonlin res_tab_lin <- data.frame(test = c("F-test on sign. of regression", "F_test on non-lin", "F-test on R^2 A","F_test on R^2 B", "F-test on slope A","F-test on slope B", "F-test on non-parallelism","F-test on preparation"), test_results = c(ifelse(p_F_regr<0.05,0,1),ifelse(p_F_nonlin<0.05,1,0), ifelse(pFR2_A<0.05,1,0),ifelse(pFR2_B<0.05,1,0), ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1), ifelse(p_F_nonp<0.05,1,0),ifelse(p_F_prep<0.05,0,1)), estimate = c(p_F_regr, p_F_nonlin,pFR2_A,pFR2_B,p_F_slope_A, p_F_slope_B,p_F_nonp,p_F_prep), Source = c("Treatment","Preparation","Regression","Non-parallelism", "Resid Error","Non-linearity","Pure error", "Total"), df = c(dfTreat,1,1,1,dfRMSE,2,dfPureE,lenCirc-1), SumSquares = c(round(SStreat,5),round(SSprep,5),round(SSreg,5), round(SSnonpar,5),round(RSS,5),round(SSnonlin,5), round(SSE,5),round(SStot,5)), MS = c(round(SStreat/dfTreat,5),round(SSprep,5),round(SSreg,5), round(SSnonpar,5),round(RSS/dfRMSE,5),round(SSnonlin/2,5), round(SSE/dfPureE,5),round(SStot/dfTotal,5)), "F-value" = c(round(F_treat,5), round(F_prep,5),round(F_regr,5), round(F_nonpar,5),"",round(F_nonlin,5),"",""), "p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "","")) RET <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2) return(RET) } PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) { mLin <- gsl_nls(readout ~ (intS+r)*isSample + intS*isRef + k*log_dose, data=circle, start=list(intS = 0, k=1,r=0), control = gsl_nls_control(xtol=1e-10,ftol=1e-10,gtol=1e-10)) # alternativ: modAB <- lm(readout ~ log_dose+isSample, circle) sum_mLin <- summary(mLin) log_dose <- unique(all_l2$log_dose) seq_x <- seq(min(log_dose), max(log_dose),0.1) SAMPLEtrue <- sigmoid[5] + (sigmoid[7]-sigmoid[5])/(1+exp(sigmoid[6]*((sigmoid[4]-sigmoid[8]-seq_x)))) REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[2]*((sigmoid[4]-seq_x)))) truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue) p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) + geom_point() + 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)) + theme_bw() p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="#4545BA", inherit.aes = F) + geom_line(data=pl_df,aes(x=lnC,y=plotT),color="#C2173F", inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4, inherit.aes = F) + labs(title = paste("unrestricted linear regression model",indS,indT), color="product") + theme(legend.position="none", axis.text = element_text(size=14)) p3 <- p2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef), size=5,alpha=0.2), inherit.aes = FALSE) + scale_shape_manual(labels=c("test","reference"), values=c(21,21)) # fit intercept for test and ref and common slope pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log_dose pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log_dose pl_rest <- data.frame(lnC=log_dose, plotS=pl_restS, plotT=pl_restT) pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="#4545BA", inherit.aes = F) + geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="#C2173F", inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4, inherit.aes = F) + labs(title = paste("restricted linear regression model",indS,indT), color="product") + 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), size=5,alpha=0.2), inherit.aes = FALSE) + scale_shape_manual(labels=c("test","reference"), values=c(21,21)) return(grid.arrange(p3,pr3,nrow=1)) } #' Calculates the potency of 4PL PLA of all models model #' #' The gradient method is used for calculating the potency for a restricted model, an unrestricteed model, #' and the log-transformations thereof. #' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. #' #' @param ro_new The dilution series readouts for REF and TEST + 1 column log(concentration). #' @param PureErrFlag Boolean if Pure Error should be used. #' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results #' 2) summary of unrestricted linear model 3) Slope difference +CIs #' 4) summary of restricted linear model. #' @export #' @examples #' ro_new <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), #' R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), #' T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), #' T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), #' log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) #' PureErrF <- TRUE #' #' #' pot4plFUNC(ro_new, PureErrF) pot4plFUNC <- function(ro_new, PureErrFlag) { all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) isSample <- rep(c(0,1),1,each=nrow(all_l)/2) all_l$isRef <- isRef all_l$isSample <- isSample all_l$Conc <- exp(all_l$log_dose) all_l$readout[all_l$readout < 0] <- 0.01 all_l$readouttrans <- log(all_l$readout) #browser() CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) if (CORdat<0) SLOPE <- -1 else SLOPE <- 1 # FITs <- Fitting_FUNC(ro_new, TransFlag = F) if (!PureErrFlag) { pot_est <- FITs[[3]] potU_est <- FITs[[4]] } else { FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) meanPureErr <- FitAnova[4,3] # SU_mr <- tryCatch({ # SU_mr <- summary(mr) # }, error = function(msg) { # return() # }) SU_mr <- FITs[[1]] SU_mrCoeff <- SU_mr$coefficients # if (length(SU_mr)>1) { # SU_mrCoeff <- SU_mr$coefficients # } else { SU_mrCoeff <- rep(NA,5) } #te2$cov.unscaled[1,1]*te2$sigma^2 #VCOV <- vcov(mr) #V_V <- VCOV/SU_mr$sigma^2 V_V <- SU_mr$cov.unscaled SEsPure <- sqrt(diag(V_V)*meanPureErr) pot_est <- c(exp(SU_mrCoeff['r',1]),exp(SU_mrCoeff['r',1]-qt(0.975,nrow(all_l)-5)*SEsPure['r']), exp(SU_mrCoeff['r',1]+qt(0.975,nrow(all_l)-5)*SEsPure['r'])) # unrestricted SU_mu <- FITs[[2]] s_muCoeff <- SU_mu$coefficients #SU_mu <- summary(mu) #VCOVu <- vcov(mu) V_Vu <- SU_mu$cov.unscaled SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) potU_est <- c(exp(s_muCoeff['r',1]),exp(s_muCoeff['r',1]-qt(0.975,nrow(all_l)-8)*SEsPureU['r']), + exp(s_muCoeff['r',1]+qt(0.975,nrow(all_l)-8)*SEsPureU['r'])) } # PureErrFlag FITstrans <- Fitting_FUNC(ro_new, TransFlag = TRUE) # pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic")) # potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic")) pot_est_log <- FITstrans[[3]] potU_est_log <- FITstrans[[4]] colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI") colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI") #browser() su_mr_log <- FITstrans[[1]] #summary(mr_log) Dat$RMSE_Rlog <- su_mr_log$sigma su_mu_log <- FITstrans[[2]] # summary(mu_log) Dat$RMSE_Ulog <- su_mu_log$sigma Dat$up_lowAslog <- su_mu_log$coefficients[1,1] - su_mu_log$coefficients[4,1] potALL <- rbind(round(pot_est,6), round(potU_est,6), round(pot_est_log,6), round(potU_est_log,6)) potALL2 <- cbind(c("restricted","unrestricted","ln-transformed restr","ln-transformed unrestr"), potALL) return(potALL2) } #' Calculates logarithmized CIs for ratios #' #' Ratio CIs can be unbound, and not calulatable. Therefore the ratios are logarithmized. #' The Covariance is not used, as it is 0 for comparisons between products. #' #' @param xs The estimate of the reference sample. #' @param xt The estimate of the test sample. #' @param se_xs The standard error of the estimate of the reference sample. #' @param se_xt The standard error of the estimate of the test sample. #' @param CoVar The covariance between test and reference estimate (set to 0). #' @param DFs Degrees of freedom, either from pure error term (no of readouts - no of concentrations*2), or RMSE term (no of readouts - 8 parameters). #' @param Conf The confidence of the Confidence Interval (default 95% which requires 0.975). #' @returns A data-frame with the lower and upper CI in anti-log form. #' @export #' @examples #' xs=2; xt=3.2; se_xt=0.34;se_xs=0.23; DFs=32-16 #' #' #' 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) { log_xs <- log(abs(xs)) log_xt <- log(abs(xt)) var_log_xs <- (se_xs/xs)^2 # approximate variance of log(xs) 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 ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) return(ci_ratio) } #' Calculates EQ tests and F-Tests for 4P #' #' CIs of asmptote ratios,slope ratios, asympt-difference ratio, and lower asympt-diff. calculated. #' In addition the F-test on significance of regression and non-linearity are calculated. #' #' @param ro_new Data frame with the concentrations and the respective readouts. #' @param Lim The limits to check if the relative CI is within the limits. #' @param PureErrFlag Boolean if Pure Error should be used. #' @returns A data-frame with the lower and upper CI in anti-log form. #' @export #' @examples #' #' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), #' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), #' 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)) #' Lim <- c(-1, 1, 0.005 , 2, 0.5, 2, 0.5,2,75,133, 75, 133) #' PureErrF <- FALSE #' #' tests_FUNC(ro_new, Lim,PureErrF) tests_FUNC <- function(ro_new, Lim, PureErrFlag) { all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) isSample <- rep(c(0,1),1,each=nrow(all_l)/2) all_l$isRef <- isRef all_l$isSample <- isSample all_l$Conc <- exp(all_l$log_dose) all_l$readout[all_l$readout < 0] <- 0.01 #browser() FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) POTr_CI <- FITs[[3]][2:3] potAll2 <- FITs[[3]] smr <- FITs[[1]] smu <- FITs[[2]] FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) # pure error pureSS <- FitAnova[4,2] pureSS_df <- FitAnova[4,1] meanPureErr <- FitAnova[4,3] #vcovMU <- vcov(mu) V_V <- smu$cov.unscaled SEsPure <- sqrt(diag(V_V)*meanPureErr) VCOVpure <- V_V*meanPureErr DFsPure <- FitAnova[4,1] testPOTr <- logical() if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100Lim[[9]]/100 & potAllU2[3] 0) { p_F_nonlin <- round(pf(F_nonlin, AnovaDFs[6], ERR_df, lower.tail = F),5) } else { p_F_nonlin <- "SSnonlin neg or single dilutions" } test_a <- test_b <- test_d <- test_ad <- logical() RSS_r <- round(sum(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) coeffs <- smu$coefficients[,1] #browser() ## EQ test on lower As diff as <- coeffs["as"] at <- coeffs["at"] lAs_diff <- (at-as) uCI_laDiff <- lAs_diff+qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2) 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 #### EQ test on upper asymptote ratio ---- # as <- coeffs["as"] # at <- coeffs["at"] # uAsRatio <- compParm(potU, "a","/",display=F) # uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2]) #browser() ds <- smu$coefficients["ds",1] dt <- smu$coefficients["dt",1] if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds","ds"]) else se_ds <- smu$coefficients["ds",2] if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt","dt"]) else se_dt <- smu$coefficients["dt",2] if (PureErrFlag) CoVarlog_d <- VCOVpure["dt","ds"] else CoVarlog_d <- vcovMU["dt","ds"] if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-8 uAsCI2 <- ParamCI_F(dt,ds,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 #### EQ test on slope ratio ---- # bs <- coeffs["bs"] # bt <- coeffs["bt"] # slopeRatio <- compParm(potU, "b","/",display=F) # slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2]) bs <- smu$coefficients["bs",1] bt <- smu$coefficients["bt",1] if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs","bs"]) else se_bs <- smu$coefficients["bs",2] if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt","bt"]) else se_bt <- smu$coefficients["bt",2] if (PureErrFlag) CoVarlog_b <- VCOVpure["bt","bs"] else CoVarlog_b <- vcovMU["bt","bs"] slopeCI2 <- ParamCI_F(bt,bs,se_bt, se_bs,CoVarlog_b, DFs, Conf=0.975) if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1 estUppA <- round(at/as,5) Dat$slopeRatioCI <- slopeCI2 #### EQ test on lower As ratio ---- # lAsRatio <- compParm(potU, "d","/",display=F) # slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2]) as <- smu$coefficients["as",1] at <- smu$coefficients["at",1] if (PureErrFlag) se_as <- sqrt(VCOVpure["as","as"]) else se_as <- smu$coefficients["as",2] if (PureErrFlag) se_at <- sqrt(VCOVpure["at","at"]) else se_at <- smu$coefficients["at",2] if (PureErrFlag) CoVarlog_a <- VCOVpure["at","as"] else CoVarlog_a <- vcovMU["at","as"] lAsCI2 <- ParamCI_F(at,as,se_at, se_as,CoVarlog_a, DFs, Conf=0.975) if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1 estLowA <- round(at/as,5) Dat$lAsCI <- lAsCI2 #### EQtest on ratio of As difference ---- AsDiffRatio <- (dt-at)/(ds-as) dt_at <- (dt-at) ds_as <- (ds-as) se_ds_asPure <- sqrt(VCOVpure["as","as"]+VCOVpure["ds","ds"]-2*VCOVpure["as","ds"]) se_dt_atPure <- sqrt(VCOVpure["at","at"]+VCOVpure["dt","dt"]-2*VCOVpure["at","dt"]) se_ds_asRMSE <- sqrt(vcovMU["as","as"]+vcovMU["ds","ds"]-2*vcovMU["as","ds"]) se_dt_atRMSE <- sqrt(vcovMU["at","at"]+vcovMU["dt","dt"]-2*vcovMU["at","dt"]) if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE AsDiffCI2 <- ParamCI_F(dt_at,ds_as,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) Dat$up_lowAs <- abs(ds-as) lowerCIlowerA <- lAsCI2[1]; lowerCIupperA <- uAsCI2[1]; upperCIlowerA <- lAsCI2[2]; upperCIupperA <- uAsCI2[2] test_lowA <- test_d; test_uppA <- test_a #browser() res_tab <- data.frame(test= c("F-test on sign. of regression*", "EQ test on lower asymptotes difference", "EQ test ratio of lower asymptotes", "EQ test ratio of Hill slopes", "EQ test ratio of upper asymptotes", "F-test on non-linearity*", "EQ test ratio of asymptote difference", "geom. rel. CI restr. model", "geom. rel. CI unrestr. model"), test_results = c(ifelse(p_F_regr<0.05,0,1), test_la_diff, test_lowA, test_b, test_uppA, ifelse(p_F_nonlin>1,1, ifelse(p_F_nonlin<0.05,1,0)), test_ad, testPOTr, test_c), estimate = c(round(p_F_regr, 3), round(lAs_diff, 5), estLowA, round(bs/bt,5), estUppA, p_F_nonlin, round(dt_at/ds_as, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)), lower_limit = c("-",Lim[[1]],Lim[[3]],Lim[[5]],Lim[[7]],"-",Lim[[11]],Lim[[9]],Lim[[9]]), upper_limit = c("-",Lim[[2]],Lim[[4]],Lim[[6]],Lim[[8]],"-",Lim[[12]],Lim[[10]],Lim[[10]]), lower_CI = c(RMSE_r, round(lCI_laDiff,3), round(lAsCI2[1],5), round(slopeCI2[1],5), round(uAsCI2[1],5), "-", round(AsDiffCI2[1],5), round(potAll2[2],2), round(potAllU2[2],2)), upper_CI = c(RMSE_pure, round(uCI_laDiff,3), round(lAsCI2[2],5), round(slopeCI2[2],5), round(uAsCI2[2],5), "-", round(AsDiffCI2[2],5), round(potAll2[3],2), round(potAllU2[3],2)) ) return(res_tab) } #' ANOVA unrestricted model #' #' ANOVA with unrestricted model. #' #' #' @param ro_new Data frame with the concentrations and the respective readouts. #' @returns A data-frame with the analysis of variance. #' @export #' @examples #' #' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), #' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), #' 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)) #' #' #' 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) isRef <- rep(c(1,0),1,each=all_len/2) isSample <- rep(c(0,1),1,each=all_len/2) all_l$isRef <- isRef all_l$isSample <- isSample all_l$Conc <- exp(all_l$log_dose) all_l$readout[all_l$readout < 0] <- 0.01 FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) smr <- FITs[[1]] smu <- FITs[[2]] smrPREDs <- FITs[[5]] smuPREDs <- FITs[[6]] # pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), # pmodels=data.frame(1,1,1,isSample)) # potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), # pmodels=data.frame(isSample, isSample,isSample,isSample)) SStreat <- round(sum((smuPREDs - mean(all_l$readout))^2),5) SStreat_df <- length(unique(all_l$log_dose))-1 SSregr <- round(sum((smrPREDs-mean(all_l$readout))^2),5) ## Non-parallel SSnonparallel <- round(sum(smr$residuals^2) - sum(smu$residuals^2),5) ## Preparation SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) ## Resid Err RSS <- round(sum(smu$residuals^2),5) RSS_df <- nrow(all_l)-SStreat_df-1 FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) # PureErr SSE <- round(FitAnova[4,3],5) SSE_df <- FitAnova[4,1] # Non-Linearity SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - smuPREDs)^2),4) LoF_df <- FitAnova[1,1]+FitAnova[2,1] ## Total SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5) MSE <- RSS/RSS_df noConc <- length(unique(all_l$Conc)) AnovaDFs <- c(noConc-1, 1,3,noConc-4-1, nrow(all_l)-noConc, noConc, nrow(all_l)-noConc-noConc, nrow(all_l)-1) p_SStreat <- round(pf((SStreat/AnovaDFs[1])/MSE, AnovaDFs[1],RSS_df, lower.tail = F),3) p_SSprep <- round(pf((SSprep/AnovaDFs[2])/MSE, AnovaDFs[2],RSS_df, lower.tail = F),3) p_SSregr <- round(pf((SSregr/AnovaDFs[3])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) p_SSnonp <- round(pf((SSnonparallel/AnovaDFs[4])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) p_SSLoF <- round(pf((SSnonlin/LoF_df)/(SSE/SSE_df), LoF_df,SSE_df, lower.tail = F),5) ANOVAtab <- data.frame(Source = c("Treatment","Preparation","Regression", "Non-Parallelism","Residual Error","Non-linearity", "Pure Error","Total"), DF = AnovaDFs, SumSquares = c(SStreat, SSprep,SSregr, SSnonparallel, RSS, SSnonlin,SSE, SStot), MeanSquares = c(round(SStreat/AnovaDFs[1],3), SSprep, round(SStreat/AnovaDFs[3],3),round(SSnonparallel/AnovaDFs[4],3), round(MSE,5), round(SSnonlin/LoF_df,5), round(SSE/SSE_df,5),""), "F-value" = c(round((SStreat/AnovaDFs[1])/MSE,5), round((SSprep/AnovaDFs[2])/MSE,5), round((SSregr/AnovaDFs[3])/MSE,5),round((SSnonparallel/AnovaDFs[4])/MSE,5), "",round((SSnonlin/LoF_df)/(SSE/SSE_df),5),"",""), "p_value" = c(round(p_SStreat,3), p_SSprep, round(p_SSregr,3), p_SSnonp,"",p_SSLoF,"","") ) return(ANOVAtab) } #' Calculates descriptive statistics per concentration #' #' Mean, SD, and CV% per concentration and product.. #' #' #' @param ro_new Data frame with the concentrations and the respective readouts. #' @returns A data-frame with the concentrations and metadata. #' @export #' @examples #' #' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), #' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), #' 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)) #' #' #' perConcTab(ro_new, noDilSeries=3) perConcTab <- function(ro_new, noDilSeries) { Reftab <- ro_new[,c(1:noDilSeries)] Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))] tReftab <- t(Reftab) colnames(tReftab) <- round(ro_new[,ncol(ro_new)],5) avs <- apply(tReftab,2,mean) sds <- apply(tReftab,2,sd) cv <- sds/avs*100 tReftab2 <- rbind(tReftab, avs,sds,cv) tTesttab <- t(Testtab) colnames(tTesttab) <- round(ro_new[,ncol(ro_new)],5) avs_test <- apply(tTesttab,2,mean) sds_test <- apply(tTesttab,2,sd) cv_test <- sds_test/avs_test*100 tTesttab2 <- rbind(tTesttab, avs_test,sds_test,cv_test) concTab <- rbind(tReftab2, tTesttab2) return(concTab) } #' Calculates dilution series. #' #' Recursive function to calcuate a geometric dilution series. #' #' @param x Starting concentration. #' @param Div Divisor to get next concentration. #' @param N Number of dilution step, increases. #' @param res Vector of results. #' @param noDil Final number of concentrations -1 (the initial one). #' @returns A data-frame with the concentrations and metadata. #' @export #' @examples #' #' x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7 #' #' divFUN(x,Div,N,res,noDil) divFUN <- function(x,Div,N,res,noDil) { N <- N+1 y <- x/Div res <- c(res,y) if (N==noDil) { return(res) } divFUN(y,Div,N,res,noDil) } #### ui ---- ui <- dashboardPage( dashboardHeader(title = "Plateflow"), dashboardSidebar( sidebarMenu( img(src="logo.png", width=230), menuItem("Home", tabName="home", icon=icon("home")), menuItem("Data template", tabName = "template", icon=icon("table"), menuSubItem( tags$li(a("EXCEL File", target="self",href="TestFile.xlsx"))) ), # menuItem("User Manual /Validation", tabName = "manual", icon=icon("book"), # tabName here and in dashboard body need to be identical # menuSubItem(icon = NULL, tags$li(a("Document", target="self",href="UserManual.pdf"))) # ), menuItem("EXCEL upload", tabName="Dataupload", icon=icon("magnet", lib="glyphicon")), menuItem("4PL metadata + report", tabName="fourPL", icon=icon("chart-line", lib="font-awesome")), #menuItem("XLSX diagnostics", tabName="XLdiagn", icon=icon("chart-bar", lib="font-awesome")), # menuItem("Linear regression + report", tabName="pla", icon=icon("pencil", lib="glyphicon")), menuItem("Wizard", tabName="wizard", icon=icon("chart-column", lib="font-awesome")), menuItem("Documentation", tabName="documentation", icon=icon("chart-area", lib="font-awesome")) ), tags$footer(HTML(paste(tags$strong(tags$u("InnerAnalytics")), paste(rep(" ",9), collapse=""), "Developer:", paste(rep(" ",9), collapse=""), "Host on:", paste(rep(" ",9), collapse=""), "Version:")), align="left", style= "position:fixed; bottom:0;width=100%; background: #FFC337BB; font-family: Times New Roman; font-size:100%; padding: 5px; color:#4545BA; box-sizing: border-box; z-index: 1000;")), dashboardBody( fluidPage( tabItems( tabItem(tabName = "home", htmlOutput("homePage")), tabItem(tabName = "Dataupload", uiOutput("Dataupload")), tabItem(tabName = "fourPL", uiOutput("fourPL")), #tabItem(tabName = "XLdiagn", uiOutput("XLdiagn")), tabItem(tabName = "pla", uiOutput("pla")), tabItem(tabName = "wizard", uiOutput("wizard")), tabItem(tabName = "documentation", uiOutput("docu")) ) ) ), skin="blue" ) #### server ---- server <- function(input, output, session) { ReportParS <- reactiveValues() IPReportParS <- reactiveValues() #### renderUIs ---- output$homePage <- renderUI({ navbarPage("Home", tabPanel("Introduction", tags$img(src="logo.png", class="adv_logo"), h4("Introduction to the bioassay software"), tags$mark("linear regression"), br(), column(4, tags$table(id="dose-table", 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), 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), 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) )) ), tabPanel("Documentation", h4("Introduction "), h4("Information on dilution setting"), "(for details see:", a(href="ADONIS.pdf","Whitepaper", download=NA, target="_blank"),")",tags$br(), "Bend points are calculated according to following formula:", withMathJax(" $$bp_{1/2} = \\pm\\frac{1.31696}{Hill's slope}$$")), tabPanel("Configuration", verbatimTextOutput("sessioninfo")) ) }) ##### UI XL ---- output$Dataupload <- renderUI({ navbarPage(title="Information", tabPanel(title = "Real data", tabsetPanel( tabPanel("Data input", column(3, #img(src="Screenshot.png", width=200), box(title = "Upload", status="warning",solidHeader = T, width=12, "Please upload your EXCEL file here", fileInput("iFile",'',accept=".xlsx")), 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")), column(6, h4("Suitability tests for 4-parametric logistic regression"), checkboxGroupInput("selectedSSTs", "Which suitability tests to be used?", choices= c("F-test on Regr."="1", "EQ-test on lower asymptote difference"= "2", "EQ-test on ratio of lower asymptote"= "3","EQ-test on ratio of Hill slopes"= "4", "EQ-test on ratio of upper asymptote"= "5", "F-test on non-linearity"="6", "EQ-test on ratio of asymptote differences"= "7"), selected= c("1","2","3","4","5","6","7")), h4("Suitability tests for Parallel Line Assay"), checkboxGroupInput("selectedSSTsLinear", "Which suitability tests to be used?", choices= c("F-test on Regr."="1", "F-test on non-linearity"= "2", "F-test on R^2 A"= "3","F-test on R^2 B"= "4", "F-test on slope A"= "5", "F-test on slope B"="6", "F-test on non-parallelism"= "7", "F-test on preparation"="8"), selected= c("1","2","3","4","5","6","7","8")) ) ), tabPanel("4pl-Analysis", tags$style(HTML("pre { color: black; background-color: #FFE1FF; font-family: 'Helvetica Neue', Helvetica, Arial, sans-serif;font-size: 12px;} ")), wellPanel( fluidRow( column(4, #h5("Diagnostics only shown if EXCEL is uploaded"), htmlOutput("PureErrW2"), tags$head(tags$style("#PureErrW2{color: red; font-size: 16px; font_style: italic;}")), tableOutput("Filesampl"), tableOutput("relpotTestTab"), plotOutput("relpotTestPlot", width="300px", height="150px"), # Pot CI plot h4("Akaike Information Criterion"), tableOutput("AIC"), h5("First row: restricted model; 2nd row: unrestricted model"), h5("Smaller values of AIC indicate better fit to the data"), tableOutput("VarDiagn") ), column(8, plotOutput("XLplot"), DTOutput("pottab4plXL"), plotOutput("diagnplot"), DTOutput("EQtests"), DTOutput("pottab4plTransXL"), tableOutput("ANOVAXLS") ) ))), tabPanel("linear Analysis", sidebarLayout( sidebarPanel( width=2, fluidRow( column(12, numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")), bsPopover(id="q4", title="", content="The calculated limits ...")), ) )), mainPanel( tabsetPanel(id="tabs", tabPanel("linear PLA", column(12, htmlOutput("PureErrW3"), tags$head(tags$style("#PureErrW3{color: red; font-size: 16px; font_style: italic;}")), plotOutput("plotLin"), "Delta method is used for potency CIs", DT::dataTableOutput("pottab"), h4("Unrestricted linear model (SSSI):"), tableOutput("SummaryModABu"), h4("Restricted linear model (CSSI):"), tableOutput("SummaryModAB"), h3("Tests for linear PLA):"), box(title="Suitability tests", status="primary",solidHeader = T, width=12, DTOutput("TESTSlin")), h5("The estimate is the p-value of the test"), h5("F-tests on regression, significance of slopes, and preparation need to have a p-value <0.05 to pass"), h5("All other tests pass if p-value > 0.05"), "SST CI for difference of slopes:", tableOutput("SlopeDiffCI"), h3("ANOVA for parallel line assay"), DTOutput("ANOVAlin")) ), tabPanel("Report", h4("Settings for report") )) ) ) ), tabPanel("parameter estimates", column(3,style = "background: #4FCBD922", br(), h4("Regression results restricted"), tableOutput("coeffs_r"), "Bend points restricted", tableOutput("bends_r2")), column(3,style = "background: #B5C74022", br(), h4("Regression results unrestricted"), tableOutput("coeffs_unr")), column(3,style = "background: #F9545422", h4("Regression results (ln-transformed)"), tableOutput("logcoeffs_r"), tableOutput("bends_unr2"), tableOutput("logcoeffs_unr")) ), tabPanel("Report", h4("Settings for report"), downloadButton("downloadXLReport", label="Download PDF report", class="butt"), tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"), ) ) ) ) }) ##### UI Meta ---- output$fourPL <- renderUI({ navbarPage(title="4PL+linear reg", tabPanel("Analysis and Plots", #sidebarLayout( # sidebarPanel( # width=4, # fluidRow( # ) # ), mainPanel(width=12, tabsetPanel(id="tabs", tabPanel("Settings", h4("Settings of 4PL regression"), div(checkboxInput("PureErrMeta", "Should pure error be used for calculation of CIs?", FALSE), style = "font-size: 24px !important;color: #C2173F"), h4("User help"), h5("If new dilutions are entered, please adjust EC50 to avoid calculation errors"), # numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")), # bsPopover(id="q4", title="", content="The calculated limits ...")), #h5("Diagnostics only shown if EXCEL is uploaded"), column(2,style = "background: #7FAEFF", #actionButton("StartCalc", "Click, when calculations to be started"), h4("curve settings"), numericInput("lowAsymptREF", "lower asymptote REF",10,step=1,min=0), numericInput("lowAsymptTEST", "lower asymptote TEST",10,step=1,min=0), numericInput("uppAsymptREF", "upper asymptote REF",110,step=1,min=0), numericInput("uppAsymptTEST", "upper asymptote TEST",110,step=1,min=0) ), column(2,style = "background: #7FAEFF", numericInput("slopeREF", "slope REF",1,step=0.1,min=-10), numericInput("slopeTEST", "slope TEST",1,step=0.1,min=-10), numericInput("EC50", "EC50 REF",-3.5), numericInput("potDiff", "potency difference",0) ), column(2,style = "background: #627ADD", h4("dilutions"), numericInput("CONC1", "highest concentration",0.3, min=-3.5), numericInput("CONC2", "2nd concentration",0.15), numericInput("CONC3", "3rd concentration",0.075), numericInput("CONC4", "4th concentration",0.0375), numericInput("CONC5", "5th concentration",0.01875), numericInput("CONC6", "6th concentration",0.00938) ), column(2,style = "background: #627ADD", numericInput("CONC7", "7th concentration",0.00469), numericInput("CONC8", "8thd concentration",0.00235), numericInput("CONC9", "9thd concentration",value=NA), numericInput("CONC10", "10th concentration",value=NA), numericInput("CONC11", "11th concentration",value=NA), numericInput("CONC12", "lowest concentration",NA) ), column(2,style = "background: #4FCBD9", h4("geometric dilution scheme"), numericInput("ConcStart", "starting concentration",value=NA, min=0), numericInput("dilutionFac", "dilution factor",value=NA, min=0, max=10), numericInput("NoDil", "no. of dilutions",value=NA, min=8), numericInput("NoDilSer", "no. of dil. series",value=NA), verbatimTextOutput("dilutions") ), column(2, h3("Settings"), helpText("Vary the slider to see the effect of special cause"), sliderInput("sdfac","Variability as % of lower to upper asymptote", max=10, value=3, min=0.1, step=0.1), checkboxInput("heterosked","heteroskedastic noise", FALSE), sliderInput("potencydiff","potency of test (%)", max=200, min=50, value=100, step=1), sliderInput("outlL","outlier in lower asymptote", min=0, max=1.5, value=0, step=0.1), sliderInput("outlM","outlier in mid part", min=0, max=1.5,value=0, step=0.1), sliderInput("outlU","outlier in upper asymptote", min=0, max=1.5,value=0, step=0.1) ), h4("log-dilutions from settings above"), verbatimTextOutput("logdil"), column(8, box(title = "Simulated data per log-concentration", status="warning",solidHeader = T, width=12, "incl. mean, sd and CV%", DT::dataTableOutput("ConctabMeta"))) #) ), #### 4pl fits ---- tabPanel("4pl-fit", tags$head(tags$style("#PureErrW2{color: red; font-size: 10px; font_style: italic;}")), wellPanel( fluidRow( column(10, tags$style(span(htmlOutput("PureErrW3"), style="color: red")), htmlOutput("PureErrW3"), plotOutput("plot4plMeta", width = "80%"), DTOutput("pottab4pl"), "Footnote: test performed on relative CIs.", DTOutput("EQtests4pl"), # SSTs h5("*...The estimate for F-test on regression and on non-linearity is the p-value"), h5("F-test on regression passes if F-value > F-crit and thus p < 0.05"), h5("F-test on non-linearity passes if F-value < F-crit and thus p > 0.05"), h5("Test results outcome: 0 ... test passed (for EQ tests: CI within limits); 1 ... test failed (for EQ tests CI not within limits); -1 ... calculations unbound/denominator too close to 0"), #plotOutput("CIplot, height=50%") ), column(8, "4 PL ANOVA unrestricted", box(title = "ANOVA unrestricted", status="warning",solidHeader = T, width=12, "", DT::dataTableOutput("ANOVA")), h5("Please note that for the CIs of rel. potency the RSS of the restricted model is used"), h5("RSS ... 'Residual error' in the SumSquares column"), h5("MSE ... 'Residual error' in the MeanSquaress column"), h5("SSE ... 'Pure error' in the SumSquares column"), h5("RMSE ... Square root of the 'Residual Error' in the MeanSquares column"), verbatimTextOutput("RMSE") ) )) ), tabPanel("ln-transformed y", h4("ln-transformed y-axis plots"), plotOutput("plot4plTransMeta", width = "80%"), DT::dataTableOutput("pottab4plTrans"), ), tabPanel("linear regression", h4("Evaluations from meta-data"), htmlOutput("PureErrW3"), tags$head(tags$style("#PureErrW3{color: red; font-size: 16px; font_style: italic;}")), column(12, plotOutput("plotLinMeta"), "Delta method is used for potency CIs", DTOutput("pottabMeta") ), column(5, h3("Tests for linear PLA:"), box(title="Suitability tests", status="primary",solidHeader = T,collapsible=T, width=12, DTOutput("TESTSlinMeta")), h5("The estimate is the p-value of the test"), h5("F-tests on regression, significance of slopes, and preparation need to have a p-value <0.05 to pass"), h5("All other tests pass if p-value > 0.05"), box(title="Unrestricted linear model (SSSI):", status="primary",solidHeader = T,collapsible=T, width=12, tableOutput("SummaryModABuMeta")), h4("Restricted linear model (CSSI):"), box(title="Restricted linear model (CSSI):", status="primary",solidHeader = T,collapsible=T, width=12, tableOutput("SummaryModABMeta")) ), column(6, h3("ANOVA for parallel line assay"), box(title="ANOVA for simultated data", status="primary",solidHeader = T, collapsible=T, width=12, DTOutput("ANOVAlinMeta")), " CI for difference of slopes:", tableOutput("SlopeDiffCIMeta"), ) ), tabPanel("Report", h4("Settings for report"), downloadButton("downloadXLReport", label="Download PDF report", class="butt"), tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"), ) ) ) #) ) ) }) output$pla <- renderUI({ navbarPage(title="pla", tabPanel("Analysis and Plots", ) ) }) output$wizard <- renderUI({ navbarPage(title="Dilution setting", tabPanel("Plots", sidebarLayout( sidebarPanel( width=3, fluidRow( column(6, numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")), bsPopover(id="q4", title="", content="The calculated limits ..."))) )), mainPanel( tabsetPanel(id="tabs", tabPanel("4pl", box(title="ANOVA table", status="primary",solidHeader = T, width=12, tableOutput("Anovatab")), column(4, h3("Confidence intervals"), tableOutput("CIs"), "The confidence intrval table is interaactive for changes in: variability slider (%SD), potency of test-slider, and 'Adjust the dilutions'-slider", tableOutput("optimalDils"), selectInput(inputId="scenario", label= "Select an 'optimal' scenario:", choices = c("scenario 2","scenario 3","scenario 6","steep slope"))), column(5, plotOutput("plotfordilutions"), h4("in grey: most extreme bend point lines of theoretical samples with 50% and 200% potency"), sliderInput("dilslider", "Adjust the dilutions(+-change in %)", min = -100,max=100, value=0, step=1, round=0), checkboxInput("fixupper","Fix highest concentration (if unticked, the center is fixed)",FALSE), h5("Dilution factors"), tableOutput("adjlogdil"), "Short guidance: wider dilution ranges increase the CIs of rel. potency, and decrease the CIs of upper and lower asymptote ratios, as well as Hill's slope ratios", br(), "Narrower dilution ranges decrease the CIs of rel. potency, and increase the CIs of upper and lower asymptote ratios, ands Hill's slope ratios",), column(3, h3("Bend points"), tableOutput("bps"), tableOutput("extremebps"), h4("Explanation of the plots") )), tabPanel("Report", h4("Settings for report") )) ) ))) }) v <- reactiveValues(num_dose=0, next.dose.t=0) sigmoid <- reactive({ sig <- c(input$lowAsymptREF, input$lowAsymptTEST,input$uppAsymptREF,input$uppAsymptREF, input$slopeREF,input$slopeTEST,input$EC50,input$potDiff) sig }) CONC <- reactive({ Konz_ <- c(input$CONC1,input$CONC2,input$CONC3,input$CONC4, input$CONC5,input$CONC6,input$CONC7,input$CONC8, input$CONC9,input$CONC10,input$CONC11,input$CONC12) if (any(na.omit(Konz_)==0)) Konz_[Konz_ ==0] <- 0.0000001 Konz <- na.omit(Konz_) }) Dils <- reactive({ Dilutions <- c(input$ConcStart,input$dilutionFac,input$NoDil,input$NoDilSer) }) #### input EXCEL file ---- observe({ if (!is.null(input$iFile)) { inFile <- input$iFile ext <- tools::file_ext(inFile$name) file.rename(inFile$datapath, paste(inFile$datapath, ".xlsx",sep="")) t.filelocation <- gsub('\\\\','/',paste(inFile$datapath, ext,sep=".")) sheets <- openxlsx::getSheetNames(t.filelocation) dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = t.filelocation) names(dat) <- sheets Dat$wb <- dat names(Dat$wb) <- sheets Dat$sheets <- sheets Dat$FileName <- input$iFile[["name"]] } }) output$sheetName <- renderUI({ if (!is.null(Dat$wb)) { #browser() cnSheets <- Dat$sheets cnSheets2 <- c("please choose", cnSheets) selectInput(inputId = "sheet", label="Select a sheet:",choices = cnSheets) } }) observeEvent(input$sign_out, { unlink(input$iFile$datapath) reset(id = "") # from shinyjs package }) #### process XLSX file ---- observe({ if (!is.null(input$iFile)) { if (!is.null(input$sheet)) { if (input$sheet != "please choose") { #browser() XLdat <- Dat$wb[input$sheet][[1]] if (is.null(XLdat)) XLdat <- Dat$wb[Dat$sheets[1]][[1]] cn <- colnames(XLdat) logI <- grep("log", cn) logDoseI <- grep("log_dose", cn) if (length(logI)>0 & length(logDoseI)==0) { XLdat$log_dose <- XLdat[,logI] XLdat2 <- XLdat[,-logI] CORro <- cor(XLdat$log_dose, XLdat[,3]) } else if (length(logI)==0 & length(logDoseI)==0) { Ind <- grep("dilu|dose|Dose|Conc|conc",cn) XLdat$log_dose <- log(XLdat[,Ind]) CORro <- cor(XLdat[,Ind], XLdat[,3]) XLdat2 <- XLdat[,-Ind] } else if (length(logI)>0 & length(logDoseI)>0) { XLdat2 <- XLdat CORro <- cor(XLdat[,logI], XLdat[,3]) } Dat$EXCEL <- XLdat2 PureErrFlag <- input$PureErr warning_text2 <- reactive({ ifelse(PureErrFlag, 'Pure Error is selected', '') }) output$PureErrW2 <- renderText(warning_text2()) noDilSeries <-(ncol(XLdat2)-1)/2 noDils <- nrow(XLdat2) Dat$noDilSeriesXL <- noDilSeries all_l <- melt(data.frame(XLdat2), id.vars="log_dose",variable.name = "replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(XLdat2)*noDilSeries) isSample <- rep(c(0,1),1,each=nrow(XLdat2)*noDilSeries) all_l$isRef <- isRef all_l$isSample <- isSample all_l$Conc <- exp(all_l$log_dose) # all_l$readout[all_l$readout < 0] <- 0.01 REP$all_l <- all_l #### XLSX eval ---- if (CORro<0) SLOPE <- -1 else SLOPE <- 1 FITs <- Fitting_FUNC(XLdat2, TransFlag = FALSE) Smr <- FITs[[1]] # summary(mr) Smu <- FITs[[2]] # summary(mu) coeffsMR <- Smr$coefficients[,1] coeffsMU <- Smu$coefficients[,1] Dat$coeffsMRes <- coeffsMR Dat$coeffsMUnr <- coeffsMU Dat$coeffs_UN <- coeffsMU names(coeffsMU) <- c("lowAsym REF", "slope REF","upperAsym REF","EC50 REF","lowAsym TEST","slope TEST","upperAsym TEST","r") if (!PureErrFlag) { pot_est <- FITs[[3]] potU_est <- FITs[[4]] colnames(pot_est) <- c("estimate","lowerCI","upperCI") colnames(potU_est) <- c("estimate","lowerCI","upperCI") } else { FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, all_l)) meanPureErr <- FitAnova[4,3] DFsPure <- FitAnova[4,1] #VCOV <- vcov(mr) V_V <- Smr$cov.unscaled #VCOV/Smr$sigma^2 VCOVpure <- V_V*meanPureErr SEsPure <- sqrt(diag(V_V)*meanPureErr) pot_est <- data.frame(estimate=exp(coeffsMR[5]), lowerCI = exp(coeffsMR[5]-qt(0.975,DFsPure)*SEsPure[5]), upperCI = exp(coeffsMR[5]+qt(0.975,DFsPure)*SEsPure[5])) #VCOVu <- vcov(mu) V_Vu <- Smu$cov.unscaled #VCOVpure <- V_Vu*meanPureErr SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) potU_est <- data.frame(estimate=exp(coeffsMU[7]), lowerCI = exp(coeffsMU[7]-qt(0.975,DFsPure)*SEsPureU[7]), upperCI = exp(coeffsMU[7]+qt(0.975,DFsPure)*SEsPureU[7])) } #browser() Dat$potDiffXL <- potU_est[1]*100 RMSE_unr_diagn <- Smu$sigma # sqrt(SU$resVar) RMSE_res_diagn <- Smr$sigma #sqrt(SR$resVar) up_lowDiffDiagn <- Smu$coefficients["ds",1]-Smu$coefficients["as",1] ProzSD_diagn <- RMSE_unr_diagn*100/up_lowDiffDiagn Dat$ProzSD_XL <- ProzSD_diagn observe({ pot_est3 <- data.frame(pot_est*100) MaxPl <- max(input$upperPot, pot_est3$upperCI) MinPl <- min(input$lowerPot, pot_est3$lowerCI) MaxPl_ <- MaxPl*1.2 MinPl_ <- MinPl*0.8 #browser() p_relCI <- ggplot(data=pot_est3, aes(xmin=lowerCI, xmax=upperCI, y=1)) + geom_linerange(size=4, col="darkseagreen",alpha=0.5) + geom_point(aes(x=estimate, y=1), col="grey15", shape=13, size=10) + geom_vline(xintercept = c(input$lowerPot, input$upperPot), col="indianred") + annotate("text", x=input$lowerPot-13, y=1.040, label=paste("lower EAC:", input$lowerPot), col="indianred") + annotate("text", x=input$upperPot+13, y=1.040, label=paste("upper EAC:", input$upperPot), col="indianred") + annotate("text", x=pot_est3$lowerCI-10, y=1.020, label=paste("lower CL:", round(pot_est3$lowerCI,1)), col="darkgreen") + annotate("text", x=pot_est3$upperCI+10, y=1.020, label=paste("upper CL:", round(pot_est3$upperCI,1)), col="darkgreen") + annotate("text", x=pot_est3$estimate, y=0.98, label=paste("rel. potency:", round(pot_est3$estimate,1)), col="black") + ylim(c(0.95, 1.05)) + xlim(c(MinPl_,MaxPl_)) + xlab("relative potency + confidence interval") + theme_bw() + theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) output$relpotTestPlot <- renderPlot({ p_relCI }) REP$relpotTestPlot <- p_relCI output$relpotTestTab <- renderTable({ pot_est3 }) }) ANOVAtab2 <- ANOVA4plUnresfunc(ro_new = XLdat2) output$ANOVAXLS <- renderTable({ ANOVAtab2 }) REP$ANOVAXLS <- ANOVAtab2 #browser() FITsTrans <- Fitting_FUNC(XLdat2, TransFlag = TRUE) # SUlog <- FITsTrans[[2]] SRlog <- FITsTrans[[1]] RMSE_unrlog_diagn <- SUlog$sigma RMSE_reslog_diagn <- SRlog$sigma up_lowDifflogDiagn <- SUlog$coefficients["ds",1]-SUlog$coefficients["as",1] ProzSDlog_diagn <- RMSE_unrlog_diagn * 100 / up_lowDifflogDiagn #### Diagnostic RMSE table #### DiagnTable <- data.frame(parameter = c("RMSE unrestricted", "RMSE_restr.", "Diff_upper-lowerAsymp", "%SD (unrestricted)", "RMSE log_unrestricted", "RMSE log_restr", "diff_up-lowAsymp_log", "%SD (log unrestricted)"), result = c(round(RMSE_unr_diagn, 4), round(RMSE_res_diagn, 4), round(up_lowDiffDiagn, 4), round(ProzSD_diagn, 4), round(RMSE_unrlog_diagn, 4), round(RMSE_reslog_diagn, 4), round(up_lowDifflogDiagn, 4), round(ProzSDlog_diagn, 4))) Dat$DiagnTable <- DiagnTable REP$DiagnTable <- DiagnTable logpotest <- FITsTrans[[3]] #exp(confintd(mrlog, "r", method = "asymptotic")) # compParm(logpot, "c") logpotUest <- FITsTrans[[4]] # exp(confintd(mulog, "r", method = "asymptotic")) # compParm(logpotu, "c") # Berechnung der Konfidenzintervalle (CI) # logpotCI <- c(exp(Smrlog[5,1] - qt(0.975, nrow(all_1)-5) * Smrlog[5,2]), exp(Smrlog[5,1]), exp(Smrlog[5,1] + qt(0.975, nrow(all_1)-5) * Smrlog[5,2])) colnames(logpotest) <- c("estimate", "lowerCI", "upperCI") colnames(logpotUest) <- c("estimate", "lowerCI", "upperCI") #browser() cnXL <- colnames(XLdat2) Filesample <- data.frame(Test = c("File name", "samples"), Test2=c(Dat$FileName, paste(cnXL[1], " vs ", cnXL[4]))) colnames(Filesample) <- c("", "") output$Filesampl <- renderTable({ Filesample }, rownames = F) UnRPLAausw <- data.frame(Information = c("model", "lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref", "lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Difference", "relative potency", "lower CI", "upper CI"), Results = unlist(c("UNRESTRICTED", round(coeffsMU, 3), round(potU_est*100, 3)))) # von psl_nls # "log relative potency", "log lower CI", "log upper CI", round(logpotest, 3), round(compParm(potu, "c", display = F), 3) output$coeffs_unr <- renderTable({ UnRPLAausw }) #browser() UnRPLAausw2 <- data.frame(Dat$bendpointsTRANS) if (length(UnRPLAausw2) > 0) { colnames(UnRPLAausw2) <- c("bendpoints log") UnrBendLog <- data.frame(Bendpoint = c("REF_lower","REF_upper", "TEST_lower","REF_lower"), bendpoints_logscale = UnRPLAausw2) output$bends_unr2 <- renderTable({ UnrBendLog }) } REP$UnRPLAausw <- UnRPLAausw REP$UnRPLAausw2 <- UnRPLAausw2 # browser() coeffs_R <- coeffsMR # pot$coefficients coeffs_R[5] <- coeffs_R[4] - coeffs_R[5] names(coeffs_R) <- c("lower A", "slope", "upper A", "EC50 REF", "EC50 TEST") # --- Ergebnistabelle: RESTRICTED --- PLAAusw <- data.frame( Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref", "EC50 Test", "relative potency", "lower CI", "upper CI"), Results = unlist(c("RESTRICTED", round(coeffs_R, 3), round(pot_est[1, ] * 100, 3)))) # von gs1_nls output$coeffs_r <- renderTable({ PLAAusw }) PLAAusw2 <- data.frame(Dat$bendpoints) output$bends_r2 <- renderTable({ PLAAusw2 }, digits = 3, rownames = T) REP$PLAausw <- PLAAusw REP$PLBend <- PLAAusw2 #### Koeffizienten-Extraktion ---- logcoeffs_R <- SRlog$coefficients[, 1] # logpot$coefficients names(logcoeffs_R) <- c("lower A", "Hill's slope", "upper A", "EC50 REF","EC50 DIFF") # --- Ergebnistabelle: LOG RESTRICTED --- LogPLAAusw <- data.frame( Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref", "EC50 difference", "log relative potency", "log lower CI", "log upper CI"), Results = unlist(c("LOG RESTRICTED", round(logcoeffs_R, 3), round(logpotest * 100, 3)))) # von gsl_nls output$logcoeffs_r <- renderTable({ LogPLAAusw }) REP$LogPLAausw <- LogPLAAusw logcoeffs_UNR <- SUlog$coefficients[,1] names(logcoeffs_UNR) <- c("lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref", "lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Diff" ) # --- Ergebnistabelle: LOG UNRESTRICTED --- LogUnrPLAAusw <- data.frame( Information = c("model", "lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref", "lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Diff" , "relative potency", "lower CI", "upper CI"), Results = unlist(c("LOG UNRESTRICTED", round(logcoeffs_UNR, 3), round(logpotUest * 100, 3)))) # von gsl_nls output$logcoeffs_unr <- renderTable({ LogUnrPLAAusw }) REP$LogUnrPLAausw <- LogUnrPLAAusw #browser() if (exists("Ind")) { Dat$dilution <- XLdat[,Ind] } else Dat$dilution <- exp(XLdat[,logI]) # --- Plot-Ausgabe --- output$XLplot <- renderPlot({ plot_f(XLdat2, TransFlag=F) }) REP$XLdat2 <- XLdat2 # --- Diagnose-Plots (Residualanalyse) --- output$diagnplot <- renderPlot({ op <- par(mfrow = c(2, 2), mar = c(3.2, 3.2, 2, .5), mgp = c(2, .7, 0)) PREDs <- FITs[[5]] PREDsU <- FITs[[6]] # 1. Residuals vs Fitted plot(Smr$residuals ~ PREDs, main = "Residuals restricted") abline(h = 0) qqnorm(Smr$residuals) qqline(Smr$residuals) plot(Smu$residuals ~ PREDsU, main = "Residuals unrestricted") abline(h = 0) qqnorm(Smu$residuals) qqline(Smu$residuals) par(op) # Parameter zurücksetzen }) pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(1,1,1,isSample)) potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(isSample, isSample,isSample,isSample)) output$AIC <- renderTable({ AIC <- AIC(pot, potU) }) output$VarDiagn <- renderTable({ DiagnTable }, digits=4) output$relpotplot <- renderPlot({ relpot(potU, intervall="fieller", bty="l", main="Quality of rel. potency over response") }) } # !please choose } # input$sheet } # input$iFile }) #### make geomDils reactive ---- observe({ #browser() if (is.null(input$ConcStart)) return(NULL) if (!is.na(input$ConcStart) & !is.na(input$dilutionFac) &!is.na(input$NoDil) &!is.na(input$NoDilSer)) { upR <- input$ConcStart noDil <- input$NoDil noDilSer <- input$NoDilSer Div <- input$dilutionFac res <- c() N_ <- 1 Conc <- c(upR, divFUN(upR,Div,N=N_,res,noDil)) Dat$MetaConc <- Conc } }) #### updateSlider on XLSX ---- # observe({ # if (!is.null(Dat$potDiffXL)) { # updateSliderInput(session, "potencydiff", # value=round(as.numeric(Dat$potDiffXL[[1]]),5)) # } # }) # observeEvent(input$potencydiff, { # if (!is.null(Dat$potDiffXL)) { # updateSliderInput(session, "potencydiff", # value=round(as.numeric(input$potencydiff),5)) # } # }) # observe({ # if (!is.null(Dat$ProzSD_XL)) { # updateSliderInput(session, "sdfac", # value=round(as.numeric(Dat$ProzSD_XL[[1]]),5)) # } # }) # observeEvent(input$sdfac, { # if (!is.null(Dat$ProzSD_XL)) { # updateSliderInput(session, "sdfac", # value=round(as.numeric(Dat$ProzSD_XL[[1]]),5)) # } # }) #### updaterNumeric Input ---- # observe({ # if(!is.null(Dat$coeffs_UN)) { # updateNumericInput(session, "lowAsymptREF", # value=round(as.numeric(Dat$coeffs_UN[3]),5), min=0) # updateNumericInput(session, "lowAsymptTEST", # value=round(as.numeric(Dat$coeffs_UN[4]),5), min=0) # updateNumericInput(session, "uppAsymptREF", # value=round(as.numeric(Dat$coeffs_UN[5]),5), min=0) # updateNumericInput(session, "uppAsymptTEST", # value=round(as.numeric(Dat$coeffs_UN[6]),5), min=0) # updateNumericInput(session, "slopeREF", # value=round(as.numeric(Dat$coeffs_UN[1]),5)) # updateNumericInput(session, "slopeTEST", # value=round(as.numeric(Dat$coeffs_UN[2]),5)) # updateNumericInput(session, "EC50", # value=round(as.numeric(Dat$coeffs_UN[7]),5)) # updateNumericInput(session, "potDiff", # value=round(as.numeric(Dat$coeffs_UN[7])- as.numeric(Dat$coeffs_UN[8]),5)) # } # }) # # observe({ # if(!is.null(Dat$dilution)) { # updateNumericInput(session, "CONC1", # value=as.numeric(Dat$dilution[1])) # updateNumericInput(session, "CONC2", # value=as.numeric(Dat$dilution[2])) # updateNumericInput(session, "CONC3", # value=as.numeric(Dat$dilution[3])) # updateNumericInput(session, "CONC4", # value=as.numeric(Dat$dilution[4])) # updateNumericInput(session, "CONC5", # value=as.numeric(Dat$dilution[5])) # updateNumericInput(session, "CONC6", # value=as.numeric(Dat$dilution[6])) # updateNumericInput(session, "CONC7", # value=as.numeric(Dat$dilution[7])) # updateNumericInput(session, "CONC8", # value=as.numeric(Dat$dilution[8])) # updateNumericInput(session, "CONC9", # value=as.numeric(Dat$dilution[9])) # updateNumericInput(session, "CONC10", # value=as.numeric(Dat$dilution[10])) # updateNumericInput(session, "CONC11", # value=as.numeric(Dat$dilution[11])) # updateNumericInput(session, "CONC12", # value=as.numeric(Dat$dilution[12])) # # } # }) observe({ if(!is.null(Dat$MetaConc)) { updateNumericInput(session, "CONC1", value=as.numeric(Dat$MetaConc[1])) updateNumericInput(session, "CONC2", value=as.numeric(Dat$MetaConc[2])) updateNumericInput(session, "CONC3", value=as.numeric(Dat$MetaConc[3])) updateNumericInput(session, "CONC4", value=as.numeric(Dat$MetaConc[4])) updateNumericInput(session, "CONC5", value=as.numeric(Dat$MetaConc[5])) updateNumericInput(session, "CONC6", value=as.numeric(Dat$MetaConc[6])) updateNumericInput(session, "CONC7", value=as.numeric(Dat$MetaConc[7])) updateNumericInput(session, "CONC8", value=as.numeric(Dat$MetaConc[8])) updateNumericInput(session, "CONC9", value=as.numeric(Dat$MetaConc[9])) updateNumericInput(session, "CONC10", value=as.numeric(Dat$MetaConc[10])) updateNumericInput(session, "CONC11", value=as.numeric(Dat$MetaConc[11])) updateNumericInput(session, "CONC12", value=as.numeric(Dat$MetaConc[12])) } }) #### render logDilsText ---- output$logdil <- renderText({ if (!is.null(Dat$MetaConc)) { Conc <- Dat$MetaConc } else Conc <- CONC() logdilu <-log(Conc) logdilu }) #### reactive dataset sim ---- sim <- reactive({ #browser() if(is.null(sigmoid())) return(NULL) sd_fac_ <- as.numeric(input$sdfac) r_ <- log(as.numeric(input$potencydiff)/100) as = sigmoid()[1]; bs = sigmoid()[5];cs = sigmoid()[7];ds = sigmoid()[3];at = sigmoid()[2]; bt = sigmoid()[6];r = sigmoid()[8]; ct = cs-r_; dt = sigmoid()[4]; if (!is.null(Dat$MetaConc)) Conc <- Dat$MetaConc else Conc <- CONC() log_conc <- log(Conc) av_test <- as + (ds-as)/(1+exp(bs*(cs - log_conc))) av_ref <- at + (dt-at)/(1+exp(bt*(ct - log_conc))) #browser() if (!is.na(input$NoDilSer)) { noDilSer <- input$NoDilSer } else if (!is.null(Dat$noDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3 if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(log_conc) isRef <- rep(c(1,0), 1,each=noDilSer*noDil) isSample <- rep(c(0,1), 1,each=noDilSer*noDil) #if (is.null(Dat$EXCEL)) { ro_new <- Calc_DilRes(as=as,at=at,ds=ds,dt=dt,cs=cs,ct=ct,r=r_,bt=bt,bs=bs, log_conc = log_conc, sd_fac=sd_fac_, # auslenkU=outlierU, # auslenkM=outlierM, # auslenkL=outlierL, heteroNoise = input$heterosked, noDilSeries = noDilSer, noDils = noDil) #} else ro_new <- Dat$EXCEL }) # }) ####sim2 ---- sim2 <- reactive({ tab <- sim() #if (is.null(Dat$EXCEL)) return(tab) else return(Dat$EXCEL) }) #### Plot 4pl ---- output$plot4plMeta <- renderPlot({ #browser() sigmoid <- sigmoid() det_sig=NULL plot_f(sim2(), TransFlag = F) }) #### Plot 4pl Transformed ---- output$plot4plTransMeta <- renderPlot({ #browser() sigmoid <- sigmoid() det_sig=NULL plot_f(sim2(), TransFlag = T) }) #### Testergebnisse für 4PL ---- observe({ if (is.null(sim2())) return(NULL) if (is.null(input$PureErrMeta)) return(NULL) #observeEvent(input$StartCalc,{ PureErrFlag <- input$PureErrMeta warning_text3 <- reactive({ ifelse(PureErrFlag, 'Pure error selected','') }) #browser() output$PureErrW3 <- renderText(warning_text3()) 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)) Dat$limite <- Limite #browser() tab <- tests_FUNC(sim2(), Limite, PureErrFlag = PureErrFlag) if (length(tab)>1) { tab[1,6:7] <- c("-","-") Dat$tests_FUNC <- tab REP$testsTab <- tab tab2 <- tab[1:7,] dat <- datatable(tab2,options = list( paging=TRUE, dom="t", rownames=FALSE )) %>% formatStyle("test_results", target='row', backgroundColor = styleEqual(c(-1,0,1), c("pink",'lightgreen','lightgrey'))) } else { dat <- datatable(data.frame(test_results = "Convergeance failed for the uploaded dataset")) } #browser() output$EQtests4pl <- renderDT({ dat}) }) # observe #### Testergebnisse für XLSX ---- observe({ if (is.null(Dat$EXCEL)) return(NULL) if (is.null(input$PureErr)) return(NULL) #observeEvent(input$StartCalc,{ PureErrFlag <- input$PureErr warning_text3 <- reactive({ ifelse(PureErrFlag, 'Pure error selected','') }) output$PureErrW3 <- renderText(warning_text3()) 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)) Dat$limite <- Limite #browser() SelTests <- as.numeric(input$selectedSSTs) tab <- tests_FUNC(Dat$EXCEL, Limite, PureErrFlag = PureErrFlag) tab[1,6:7] <- c("-","-") Dat$tests_FUNC <- tab REP$testsTab <- tab tab2 <- tab[SelTests,] dat <- datatable(tab2,options = list( paging=TRUE, dom="t", rownames=FALSE )) %>% formatStyle("test_results", target='row', backgroundColor = styleEqual(c(-1,0,1), c("pink",'lightgreen','lightgrey'))) output$EQtests <- renderDT({ dat }) }) # observe ####plot CIs ---- observe({ tab <- Dat$tests_FUNC if (is.null(tab)) return(NULL) tab2 <- tab[-c(1,2,3,6),] tab2[,3:7] <- apply(tab2[,3:7],2,as.numeric) tab2[4:5,3:7] <- tab2[4:5,3:7]/100 p_CIs <- ggplot(tab2,aes(x=test,y=estimate, color=test,group=test)) + geom_point() + geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=0.4) + geom_crossbar(aes(ymin=lower_limit, ymax=upper_limit), size=0.8) + coord_flip() + theme_bw() + theme(legend.position = "none",text = element_text(size=20)) output$CIplot <- renderPlot({ p_CIs}, height=200) REP$CIplot <- p_CIs }) output$simdat <- DT::renderDataTable({ tab <- sim2() if (is.character(tab)) stop(tab) tab2 <- round(tab, 5) colnames(tab2) <- c(paste("T", seq(1,(ncol(tab2)-1)/2)), paste("R", seq(1,(ncol(tab2)-1)/2)), "log_conc" ) dat <- datatable(tab2, options=list( paging=T, pageLength=20, dom="t" )) }) ##### Concentrationtab Meta ---- output$ConctabMeta <- DT::renderDataTable({ if (!is.na(Dils()[1]) & is.na(Dils()[4])) return(NULL) tab <- sim2() if (is.character(tab)) stop(tab) if (!is.na(Dils()[4])) { noDilSer <- Dils()[4] } else if (!is.null(Dat$noDilSeriesXL)) { noDilSer <- Dat$noDilSeriesXL } else { noDilSer <- 3 } Conc <- CONC() Conctab <- perConcTab(tab, noDilSeries = noDilSer) Dat$Conctab <- Conctab dat <- datatable(Conctab, options=list( paging=T, pageLength=12, dom="t" )) %>% formatStyle(0, target='row', backgroundColor = styleEqual(c("avs","sds","cv", "avs_test","sds_test","cv_test"), c('lightgrey','lightgreen','pink','lightgrey','lightgreen','pink')) ) %>% formatRound(columns=colnames(Conctab), digits=3) }) output$Conctab <- DT::renderDataTable({ if (!is.na(Dils()[1]) & is.na(Dils()[4])) return(NULL) tab <- sim2() if (is.character(tab)) stop(tab) if (!is.na(Dils()[4])) { noDilSer <- Dils()[4] } else if (!is.null(Dat$noDilSeriesXL)) { noDilSer <- Dat$noDilSeriesXL } else { noDilSer <- 3 } Conc <- CONC() Conctab <- perConcTab(tab, noDilSeries = noDilSer) Dat$Conctab <- Conctab dat <- datatable(Conctab, options=list( paging=T, pageLength=12, dom="t" )) %>% formatStyle(0, target='row', backgroundColor = styleEqual(c("avs","sds","cv", "avs_test","sds_test","cv_test"), c('lightgrey','lightgreen','pink','lightgrey','lightgreen','pink')) ) %>% formatRound(columns=colnames(Conctab), digits=3) }) #### process XL linearly, Plot output ---- output$plotLin <- renderPlot({ if (is.null(Dat$EXCEL)) return(NULL) tab <- Dat$EXCEL # tab <- sim2() if (is.character(tab)) stop(tab) #browser() log_conc <- tab$log_dose noDilSer = (ncol(tab)-1)/2 noDil <- nrow(tab) Conctab <- perConcTab(tab, noDilSer) # if (!is.na(Dils()[3])) noDil <- Dils()[3] else noDil = length(Conc) # slopeSt <- slopeTe <- matrix(NA, nrow=noDil-2,ncol=2) for (i in 1:(noDil-2)) { avs <- Conctab[noDilSer+1,] threes <- data.frame(lnC=log_conc[i:(i+2)], resp=avs[i:(i+2)]) lm3St <- lm(resp ~ lnC, data=threes) slopeSt[i,] <- lm3St$coefficients avt <- Conctab[noDilSer*2+4,] threet <- data.frame(lnC=log_conc[i:(i+2)], resp=avt[i:(i+2)]) lm3Te <- lm(resp ~ lnC, data=threet) slopeTe[i,] <- lm3Te$coefficients } indS <- which(abs(slopeSt[,2]) == max(abs(slopeSt[,2]))) indT <- which(abs(slopeTe[,2]) == max(abs(slopeTe[,2]))) 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) 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_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),] all_mT <- all_l2T[order(all_l2T$log_dose, decreasing=TRUE),] circleS <- all_mS[(indS*noDilSer-(noDilSer-1)):((indS+2)*noDilSer),] circleT <- all_mT[(indT*noDilSer-(noDilSer-1)):((indT+2)*noDilSer),] circle <- rbind(circleS,circleT) Dat$circles <- circle #browser() sigmoid <- Dat$coeffsMUnr pLin <- PlotLinPLA_FUNC(circle, sigmoid = sigmoid, all_l2, pl_df, indS, indT) pLin }) #### process metadata, Plot output ---- output$plotLinMeta <- renderPlot({ tab <- sim2() if(is.null(tab)) return(NULL) if (is.character(tab)) stop(tab) #browser() if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer = (ncol(tab)-1)/2 Conc <- CONC() log_conc <- log(Conc) Conctab <- perConcTab(tab, noDilSer) if (!is.na(Dils()[3])) noDil <- Dils()[3] else noDil = length(Conc) slopeSt <- slopeTe <- matrix(NA, nrow=noDil-2,ncol=2) for (i in 1:(noDil-2)) { avs <- Conctab[noDilSer+1,] threes <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avs[i:(i+2)]) lm3St <- lm(resp ~ lnC, data=threes) slopeSt[i,] <- lm3St$coefficients avt <- Conctab[noDilSer*2+4,] threet <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avt[i:(i+2)]) lm3Te <- lm(resp ~ lnC, data=threet) slopeTe[i,] <- lm3Te$coefficients } indS <- which(abs(slopeSt[,2]) == max(abs(slopeSt[,2]))) indT <- which(abs(slopeTe[,2]) == max(abs(slopeTe[,2]))) 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) 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_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),] all_mT <- all_l2T[order(all_l2T$log_dose, decreasing=TRUE),] circleS <- all_mS[(indS*noDilSer-(noDilSer-1)):((indS+2)*noDilSer),] circleT <- all_mT[(indT*noDilSer-(noDilSer-1)):((indT+2)*noDilSer),] circle <- rbind(circleS,circleT) Dat$circlesMeta <- circle sigmoid <- sigmoid() #browser() pLin2 <- PlotLinPLA_FUNC(circle, sigmoid = sigmoid, all_l2, pl_df,indS, indT) pLin2 }) #### linear PLA tests Metadata ---- output$TESTSlinMeta <- DT::renderDataTable({ tab <- sim2() if (is.null(tab)) return(NULL) Conc <- CONC() Limite <- Dat$limite circlesMeta <- Dat$circlesMeta PureErrFlag <- input$PureErrMeta warning_text <- reactive({ ifelse(PureErrFlag, 'Pure error is selected','') }) output$PureErrW <- renderText(warning_text()) browser() LIN <- ANOVAlintests(tab,circlesMeta,Limite,PureErrFlag=PureErrFlag) df <- LIN[[1]] su_modU <- LIN[[2]] su_mod2 <- LIN[[4]] output$SummaryModABuMeta <- renderTable({ su_modU }, digits=5) output$SummaryModABMeta <- renderTable({ su_mod2 }, digits=5) slopeDiffCI <- t(data.frame(LIN[[3]])) colnames(slopeDiffCI) <- c("slope difference","lower CI","upper CI") output$SlopeDiffCIMeta <- renderTable({ slopeDiffCI },digits=5) SelTestsL <- as.numeric(input$selectedSSTsLinear) df2 <- df Dat$ANOVAMeta <- df[,4:length(df)] dat <- datatable(df2[,1:3], options=list( paging=T, dom="t",rownames=F )) %>% formatStyle("test_results", target="row",backgroundColor = styleEqual(c(-1,0,1), c("pink","lightgreen","lightgrey"))) }) #### linear PLA tests XLinput ---- output$TESTSlin <- DT::renderDataTable({ tab <- Dat$EXCEL if (is.character(tab)) stop(tab) Conc <- exp(tab$log_dose) 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)) noDil <- nrow(tab) noDilSer <- Dat$noDilSeriesXL Conctab <- perConcTab(tab, noDilSeries = noDilSer) #browser() slopeSt <- slopeTe <- matrix(NA, nrow=noDil-2,ncol=2) for (i in 1:(noDil-2)) { avs <- Conctab[noDilSer+1,] threes <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avs[i:(i+2)]) lm3St <- lm(resp ~ lnC, data=threes) slopeSt[i,] <- lm3St$coefficients avt <- Conctab[noDilSer*2+4,] threet <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avt[i:(i+2)]) lm3Te <- lm(resp ~ lnC, data=threet) slopeTe[i,] <- lm3Te$coefficients } indS <- which(abs(slopeSt[,2]) == max(abs(slopeSt[,2]))) indT <- which(abs(slopeTe[,2]) == max(abs(slopeTe[,2]))) # 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) 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_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),] all_mT <- all_l2T[order(all_l2T$log_dose, decreasing=TRUE),] circleS <- all_mS[(indS*noDilSer-(noDilSer-1)):((indS+2)*noDilSer),] circleT <- all_mT[(indT*noDilSer-(noDilSer-1)):((indT+2)*noDilSer),] circle <- rbind(circleS,circleT) PureErrFlag <- input$PureErr warning_text <- reactive({ ifelse(PureErrFlag, 'Pure error is selected','') }) output$PureErrW3 <- renderText(warning_text()) #browser() LIN <- ANOVAlintests(tab,circle,Limite,PureErrFlag=PureErrFlag) df <- LIN[[1]] su_modU <- LIN[[2]] su_mod2 <- LIN[[4]] output$SummaryModABu <- renderTable({ su_modU }, digits=5) output$SummaryModAB <- renderTable({ su_mod2 }, digits=5) slopeDiffCI <- t(data.frame(LIN[[3]])) colnames(slopeDiffCI) <- c("slope difference","lower CI","upper CI") output$SlopeDiffCI <- renderTable({ slopeDiffCI },digits=5) SelTestsL <- as.numeric(input$selectedSSTsLinear) df2 <- df[SelTestsL,] Dat$ANOVA <- df[,4:length(df)] dat <- datatable(df2[,1:3], options=list( paging=T, dom="t",rownames=F )) %>% formatStyle("test_results", target="row",backgroundColor = styleEqual(c(-1,0,1), c("pink","lightgreen","lightgrey"))) }) #### output 4PL ANOVA tests Meta ---- output$ANOVA <- DT::renderDataTable({ sigmoid <- sigmoid() tab <- ANOVA4plUnresfunc(sim2()) # ,sigmoid dat <- datatable(tab, options=list( dom="t",rownames=F )) %>% formatStyle("p_value", target="row", backgroundColor = styleEqual(c("p_value"), c("lightgrey"))) }) #### output 4PL ANOVA tests XL ---- # not needed # output$ANOVA_XL <- DT::renderDataTable({ # tab <- Dat$EXCEL # tab <- ANOVA4plUnresfunc(sim2()) # ,sigmoid # dat <- datatable(tab, # options=list( # dom="t",rownames=F # )) %>% formatStyle("p_value", target="row", # backgroundColor = styleEqual(c("p_value"), # c("lightgrey"))) # }) #### output RMSEs ---- output$RMSE <- renderText({ paste("RMSE (unrestricted model):", Dat$RMSE_unr, "(Use it to compare against RMSE restr. model for non-parallelism)\n", "RMSE (restricted model):", Dat$RMSE_r, "\n", "Pure RMSE (unrestricted model):", Dat$RMSE_pure, "\n", "%SD (unr. model): ", Dat$RMSE_unr*100/Dat$up_lowAs, "(calculated as: RMSE/(upper-lower Asymptote)*100\n", "RMSE (log restr. model): ", Dat$RMSE_Rlog, "\n", "RMSE (log unrestr. model): ", Dat$RMSE_Ulog, "\n", "%SDlog (unr model): ", Dat$RMSE_Ulog*100/Dat$up_lowAslog ) }) output$ANOVAlin <- DT::renderDataTable({ if (is.null(Dat$ANOVA)) return(NULL) ANOVAlin <- Dat$ANOVA dat <- datatable(ANOVAlin, options=list( dom="t",rownames=F )) %>% formatStyle("p.value", target='cell', backgroundColor = styleEqual(c("p.value"), c("lightgrey"))) }) output$ANOVAlinMeta <- DT::renderDataTable({ ANOVAlin <- Dat$ANOVAMeta dat <- datatable(ANOVAlin, options=list( dom="t",rownames=F )) %>% formatStyle("p.value", target='cell', backgroundColor = styleEqual(c("p.value"), c("lightgrey"))) }) #### output Lin pot tab XL ---- output$pottab <- DT::renderDataTable({ if (is.null(Dat$circles)) return(NULL) Lim <- 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)) circles <- Dat$circles PureErrFlag <- input$PureErr pottab <- LinPotTab(circles,Lim,PureErrFlag = PureErrFlag) #browser() dat <- datatable(pottab, options=list( dom="t",rownames=F )) %>% formatStyle("test_result", target='row', backgroundColor = styleEqual(c(0,1), c("lightgrey","#F9545488"))) }) ### output pot tab Meta ---- output$pottabMeta <- DT::renderDataTable({ Lim <- 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)) circles <- Dat$circlesMeta PureErrFlag <- input$PureErrMeta pottab <- LinPotTab(circles,Lim,PureErrFlag = PureErrFlag) #browser() dat <- datatable(pottab, options=list( dom="t",rownames=F )) %>% formatStyle("test_result", target='row', backgroundColor = styleEqual(c(0,1), c("lightgrey","#F9545488"))) }) #### 4pl potency table Meta ---- observe({ #browser() if (is.null(sim2()) | is.null(Dils())) return(NULL) ro_new <- sim2() Dils_ <- Dils() if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer <- 3 PureErrFl <- input$PureErrMeta pottab4 <- pot4plFUNC(ro_new = ro_new, PureErrFlag = PureErrFl) #browser() Lim <- 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)) pottab4_ <- data.frame(pottab4) pottab4_$potency <- round(as.numeric(pottab4[,2])*100,1) pottab4_$`lower95%CI` <- round(as.numeric(pottab4[,3])*100,2) 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) if (as.numeric(pottab4_$relative_lowerCL[1]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[1]) < Lim[[10]] ) { test_potCI <- 0 } else {test_potCI <- 1 } if (as.numeric(pottab4_$relative_lowerCL[2]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[2]) < Lim[[10]] ) { test_potUCI <- 0 } else {test_potUCI <- 1 } if (as.numeric(pottab4_$relative_lowerCL[3]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[3]) < Lim[[10]] ) { test_potCI_t <- 0 } else {test_potCI_t <- 1 } if (as.numeric(pottab4_$relative_lowerCL[4]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[4]) < Lim[[10]] ) { test_potUCI_t <- 0 } else {test_potUCI_t <- 1 } pottab4_ <- cbind(pottab4_[,-(2:4)], data.frame(tests=c(test_potCI, test_potUCI,test_potCI_t,test_potUCI_t))) colnames(pottab4_) <- c("model","potency","lower95%CI","upper95%CI","relative_lower95%CI","relative_upper95%CI","test_result") output$pottab4pl <- DT::renderDataTable({ dat <- datatable(pottab4_[1:2,], options=list(digits=3, paging=T, dom="t",rownames=F )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), c("lightgreen","pink"))) }) output$pottab4plTrans <- DT::renderDataTable({ dat <- datatable(pottab4_[3:4,], options=list(digits=3, paging=T, dom="t",rownames=F )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), c("lightgreen","pink"))) }) }) #### 4pl potency table XL ---- observe({ #browser() if (is.null(Dat$EXCEL)) return(NULL) ro_new <- Dat$EXCEL noDilSer <- Dat$noDilSeriesXL PureErrFl <- input$PureErr pottab4 <- pot4plFUNC(ro_new = ro_new, PureErrFlag = PureErrFl) #browser() Lim <- 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)) pottab4_ <- data.frame(pottab4) pottab4_$potency <- round(as.numeric(pottab4[,2])*100,1) pottab4_$`lower95%CI` <- round(as.numeric(pottab4[,3])*100,2) 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) if (as.numeric(pottab4_$relative_lowerCL[1]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[1]) < Lim[[10]] ) { test_potCI <- 0 } else {test_potCI <- 1 } if (as.numeric(pottab4_$relative_lowerCL[2]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[2]) < Lim[[10]] ) { test_potUCI <- 0 } else {test_potUCI <- 1 } if (as.numeric(pottab4_$relative_lowerCL[3]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[3]) < Lim[[10]] ) { test_potCI_t <- 0 } else {test_potCI_t <- 1 } if (as.numeric(pottab4_$relative_lowerCL[4]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[4]) < Lim[[10]] ) { test_potUCI_t <- 0 } else {test_potUCI_t <- 1 } pottab4_ <- cbind(pottab4_[,-(2:4)], data.frame(tests=c(test_potCI, test_potUCI,test_potCI_t,test_potUCI_t))) colnames(pottab4_) <- c("model","potency","lower95%CI","upper95%CI","relative_lower95%CI","relative_upper95%CI","test_result") output$pottab4plXL <- DT::renderDataTable({ dat <- datatable(pottab4_[1:2,], options=list(digits=3, paging=T, dom="t",rownames=F )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), c("#B5C74055","#F9545455"))) }) output$pottab4plTransXL <- DT::renderDataTable({ dat <- datatable(pottab4_[3:4,], options=list(digits=3, paging=T, dom="t",rownames=F )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), c("#B5C74055","#F9545455"))) }) }) #### Dilutions Simulator ---- output$plotfordilutions <- renderPlot({ tab <- sim2() #browser() tab <- as.data.frame(tab) dils <- tab$log_dose min_y <- min(tab[,1:3]) max_y <- max(tab[,1:3]) 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 { if (!is.null(Dat$cfordils)) { av <- Dat$cfordils } else { av <- (min(dils) + max(dils))/2 } dils_av <- dils-av dils_avsc <- dils_av*(input$dilslider/100+1) dils2 <- dils_avsc+av dilfactors <- 1/exp(dils2-lag(dils2)) } Dat$newDils <- dils2 sigmoid <- sigmoid() #browser() BPs <- Dat$bendpoints EC50REF <- (BPs[2]+BPs[1])/2 Einh <- abs((BPs[2]-BPs[1])/5) asyml <- EC50REF-2*(EC50REF-BPs[1]) asymu <- EC50REF+2*(EC50REF-BPs[1]) det_sig <- Dat$coeffs_UN if (is.null(Dat$coeffs_UN)) { SAMPLE50 <- sigmoid[1] + (sigmoid[3] - sigmoid[1])/(1+exp(sigmoid[5]*( (sigmoid[7]+0.693147)- dils2))) SAMPLE200 <- sigmoid[1] + (sigmoid[3] - sigmoid[1])/(1+exp(sigmoid[5]*( (sigmoid[7]-0.693147)-dils2))) Xbend50l <- sigmoid[7] + 0.693147-1.5434/sigmoid[5] Xbend200l <- sigmoid[7] - 0.693147-1.5434/sigmoid[5] Xbend50u <- sigmoid[7] + 0.693147+1.5434/sigmoid[5] Xbend200u <- sigmoid[7] - 0.693147+1.5434/sigmoid[5] Xbend50 <- max(Xbend50l, Xbend50u) Xbend200 <- min(Xbend200l, Xbend200u) dummy <- plot_f(tab) } else { #browser() SAMPLE50 <- det_sig[3] + (det_sig[5] - det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]+0.693147-dils2))) SAMPLE200 <- det_sig[3] + (det_sig[5] - det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-0.693147-dils2))) Xbend50l <- det_sig[7] + 0.693147-1.5434/det_sig[1] Xbend200l <- det_sig[7] - 0.693147-1.5434/det_sig[1] Xbend50u <- det_sig[7] + 0.693147+1.5434/det_sig[1] Xbend200u <- det_sig[7] - 0.693147+1.5434/det_sig[1] Xbend50 <- max(Xbend50l, Xbend50u) Xbend200 <- min(Xbend200l, Xbend200u) dummy <- plot_f(tab) } pl_df <- cbind(dils2, SAMPLE50, SAMPLE200) #browser() # scenario2 eqSpac <- abs((BPs[1]-BPs[2])/5) optdils <- c((asyml+BPs[1])/2, BPs[1], BPs[1]+1*eqSpac, BPs[1]+2*eqSpac,BPs[1]+3*eqSpac,BPs[1]+4*eqSpac,BPs[2], (asymu+BPs[2])/2) # scenario 3 eqSpac_3 <- abs((BPs[1]-BPs[2])/3) optdils_3 <- c(BPs[1]-2*eqSpac_3, BPs[1]-eqSpac_3, BPs[1], BPs[1]+1*eqSpac_3, BPs[1]+2*eqSpac_3,BPs[2], BPs[2]+eqSpac_3, BPs[2]+2*eqSpac_3) # scenario 6 Einh2 <- abs(((BPs[2]-BPs[1])*0.7)/5) eqSpac2 <- (2*0.7/Einh)/3 optdils2 <- c((asyml+BPs[1])/2, BPs[1], EC50REF-1.5*Einh2, EC50REF-0.5*Einh2,EC50REF+0.5*Einh2,EC50REF+1.5*Einh2, BPs[2], (asymu+BPs[2])/2) # steep slope eqSpac3 <- (abs(Xbend200-Xbend50))/5 optdils3 <- c(Xbend200-eqSpac3,Xbend200, Xbend200+1*eqSpac3, Xbend200+2*eqSpac3,Xbend200+3*eqSpac3,Xbend200+4*eqSpac3,Xbend50, Xbend50+eqSpac3) output$extremebps <- renderTable({ ExtremeBPs <- c(Xbend50,Xbend200) DF2 <- data.frame(sample=c("50% sample (right)", "200% sample (left)"), Extreme_BPs=ExtremeBPs) DF2 }) optD <- data.frame(cbind(optdils, optdils_3,optdils2, optdils3)) colnames(optD) <- c("scenario2","scenario3","scenario6","steep slope") output$optimalDils <- renderTable({ optD }) output$adjlogdil <- renderTable({ adjlogdilfactors <- round(dilfactors,3) adjlogdils <- round(dils2,3) adjdils <- round(exp(dils2),3) DilsTable <- data.frame('adjusted ln(dilutions)' = adjlogdils, 'adjusted ln_dilution_factors' = adjlogdilfactors, 'adjusted dilutions' = adjdils) DilsTable }) if (!is.null(Dat$p2)) { p2 <- Dat$p2 p_dil <- p2 + annotate("pointrange",x=dils2,y=rep(min_y, length(dils2)), xmin=min(dils2), xmax=max(dils2)) + annotate("text", x=dils2,y=rep(min_y+(max_y-min_y)*0.05, length(dils2)), label=as.character(round(dils2,3))) + 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))) + geom_line(data=as.data.frame(pl_df),aes(x=dils2,y=SAMPLE50), color="grey15", linetype=2, inherit.aes = F) + geom_line(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")} + {if (input$scenario =="scenario 2") annotate("pointrange",x=optdils,y=rep(min_y+(max_y-min_y)*0.2, length(optdils)), xmin=min(optdils), xmax=max(optdils), color="seagreen")} + {if (input$scenario =="scenario 2") annotate("text",x=optdils,y=rep(min_y+(max_y-min_y)*0.25, length(optdils)), label=as.character(round(optdils,3)), color="seagreen")} + {if (input$scenario =="scenario 3") annotate("pointrange",x=optdils_3,y=rep(min_y+(max_y-min_y)*0.2, length(optdils_3)), xmin=min(optdils_3), xmax=max(optdils_3), color="seagreen")} + {if (input$scenario =="scenario 3") annotate("text",x=optdils_3,y=rep(min_y+(max_y-min_y)*0.25, length(optdils_3)), label=as.character(round(optdils_3,3)), color="seagreen")} + {if (input$scenario =="steep slope") annotate("pointrange",x=optdils3,y=rep(min_y+(max_y-min_y)*0.2, length(optdils3)), xmin=min(optdils3), xmax=max(optdils3), color="seagreen")} + {if (input$scenario =="steep slope") annotate("text",x=optdils3,y=rep(min_y+(max_y-min_y)*0.25, length(optdils3)), label=as.character(round(optdils3,3)), color="seagreen")} + annotate("text",x=optdils[1],y=(max_y+min_y)*0.5, label=paste("in green: optimal \n dilutions acc. to Whitepaper\n", input$scenario), color="seagreen", size=14/.pt,fontface="bold") } print(p_dil) }) #### Dilutions CI table ---- observe({ if (is.null(input$potencydiff)) return(NULL) output$CIs <- renderTable({ PureErrFlag <- input$PureErr if (is.null(Dat$coeffs_UN)) { # checks if an EXCEL was uploaded sigmoid <- sigmoid() det_sig=NULL ast = sigmoid()[1];bst = sigmoid()[5];cst = sigmoid()[7];dst = sigmoid()[3];ate = sigmoid()[2]; bte = sigmoid()[6];r_ = sigmoid()[8]; cte = cst-r_;dte = sigmoid()[4]; } else { sigmoid <- NULL det_sig <- Dat$coeffs_UN ast <- det_sig[3] ate <- det_sig[4] bst <- det_sig[1] bte <- det_sig[2] cst <- det_sig[7] cte <- det_sig[7] -log(input$potencydiff/100) dst <- det_sig[5] dte <- det_sig[6] r_ <- log(input$potencydiff/100) } if (!is.na(input$NoDilSer)) { noDilSer <- input$NoDilSer } else if (!is.null(Dat$NoDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3 if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(Dat$newDils) #browser() tab <- Calc_DilRes(as=ast,at=ate,ds=dst,dt=dte,cs=cst,ct=cte,r=r_,bt=bte,bs=bst, sd_fac=input$sdfac,log_conc=Dat$newDils, # auslenkU=outlierU, # auslenkM=outlierM, # auslenkL=outlierL, heteroNoise = FALSE, noDilSeries = noDilSer, noDils = noDil) 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)) CItable <- tests_FUNC(tab,Limite,PureErrFlag=PureErrFlag) CItable_ <- CItable[-c(1,2,6,8,9),-c(2,4,5)] potAll <- pot4plFUNC(tab, input$PureErr) restrPot <- potAll[1,1:4] restrPot[2:4] <- round(as.numeric(restrPot[2:4]),5) potAll_ <- rbind(CItable_, restrPot) potAll_$CIwidth <- as.numeric(potAll_[,4])-as.numeric(potAll_[,3]) potAll_[,1] <- c("ratio of lower asymptotes","ratio of slopes","ratio of upper asymptotes", "ratio of asympt. differences","restricted potency") output$bps <- renderTable({ DF <- data.frame(sample=names(Dat$bendpoints),BPs=Dat$bendpoints) DF }) return(potAll_) }) }) #### simulations ---- observe({ observeEvent(input$goSim,{ sd_fac_ <- as.numeric(input$sdfac) r_ <- log(as.numeric(input$potencydiff)/100) Conc <- Dat$MetaConc as = sigmoid()[1]; bs = sigmoid()[5];cs = sigmoid()[7];ds = sigmoid()[3];at = sigmoid()[2]; bt = sigmoid()[6];r = sigmoid()[8]; ct = cs-r_; dt = sigmoid()[4] if (!is.null(Dat$MetaConc)) { Conc <- Dat$MetaConc } else { Conc <- CONC() } log_dose <- log(Conc) yAxfac <- (ds-as) if (!is.na(input$NoDilSer)) { noDilSer <- input$NoDilSer } else if (!is.null(Dat$NoDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3 if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(Conc) isRef <- rep(c(1,0),1,each=noDilSer*noDil) isSample <- rep(c(0,1),1,each=noDilSer*noDil) N <- as.numeric(input$simN) av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose))) resHist <- matrix(NA,nrow=N, ncol=13) residualsList <- list() start.time2 <- Sys.time() withProgress(message = 'Making plot', value=0, { for (i in 1:N) { if (input$heterosked) { # heterosc noise ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac_/100))), nrow=noDil, ncol=noDilSer*2) } else { # homosc noise ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac_*yAxfac/100))), nrow=noDil, ncol=noDilSer*2) } # browser() ro_jit <- abs(ro_jit) ro_new <- cbind(ro_jit, log_dose) all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name = "replname", value.name = "readout") all_l$isRef <- isRef all_l$isSample <- isSample all_l$Conc <- exp(all_l$log_dose) pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(1,1,1,isSample)) potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE) potAll2 <- potAll[1:3] RSS <- sum(pot$predres[,2]^2) dfreed <- nrow(all_l)-5 MSE <- RSS/dfreed potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(isSample, isSample,isSample,isSample)) DF_U <- nrow(all_l)-8 uAsratio <- compParm(potU, "a",display=F) uCIuAs <- uAsratio[1]+qt(0.975,DF_U)*uAsratio[2] lCIuAs <- uAsratio[1]-qt(0.975,DF_U)*uAsratio[2] lAsratio <- compParm(potU, "d",display=F) uCIlAs <- lAsratio[1]+qt(0.975,DF_U)*lAsratio[2] lCIlAs <- lAsratio[1]-qt(0.975,DF_U)*lAsratio[2] Sloperatio <- compParm(potU, "b",display=F) uCISlo <- Sloperatio[1]+qt(0.975,DF_U)*Sloperatio[2] lCISlo <- Sloperatio[1]-qt(0.975,DF_U)*Sloperatio[2] su <- summary(potU) v <- vcov(potU)[c(5,6),c(5,6)] Vd <- vcov(potU)[c(3,4),c(3,4)] Va_d <- v+Vd A_DTEST <- su$coefficients[6,1]-su$coefficients[4,1] A_DREF <- su$coefficients[5,1]-su$coefficients[3,1] if (abs(at/(sqrt(Va_d[2,2]/3))) > qt(0.95,2)) { try(Fie_ad <- round(FiellerRatio(A_DREF,A_DTEST, Va_d),5)) } if (!exists("Fie_ad")) Fie_ad <- NA resHist[i,] <- c(potAll2, sqrt(MSE),Sloperatio[1],lCISlo, uCISlo, uAsratio[1], lCIuAs, uCIuAs, Fie_ad[1],Fie_ad[2],Fie_ad[3]) colnames(resHist) <- c("pot4pl","lCI4pl","uCI4pl","RMSE","estSlope_ratio", "lCISlope_ratio","uCISlope_ratio","estuAs_ratio", "lCIuAs_ratio","uCIuAs_ratio","estAsyDiff_ratio", "lCIAsyDiff_ratio", "uCIAsyDiff_ratio") incProgress(1/N, detail=paste("Doing simulations",i)) } # withProgress }) end.time2 <- Sys.time() Dat$resHist <- resHist }) }) #### simulation Histograms output ---- output$plotHistuAs <- renderPlot({ if (!is.null(Dat$resHist)) { resHist <- Dat$resHist #browser() resHistuAs <- as.data.frame(resHist[,8:10]) resHistuAs_l <- melt(data.frame(resHistuAs), variable.name="ratio_CIs", value.name = "readout") #browser() lowquant_uAs <- quantile(resHistuAs[,2], probs=as.numeric(input$lowQuant)/100) upquant_uAs <- quantile(resHistuAs[,3], probs=as.numeric(input$uppQuant)/100) p_uAs <- ggplot(resHistuAs_l) + geom_histogram(aes(readout, fill=ratio_CIs),alpha=0.5,position="identity") + labs(title = paste("upper asymptote ratio EACs:", round(lowquant_uAs,3), " to ", round(upquant_uAs,3))) + geom_vline(xintercept = c(lowquant_uAs, upquant_uAs), color="black", linetype="dashed", linewidth=1) + geom_vline(xintercept = c(input$lEACratioua , input$uEACratioua), color="red", linetype="dashed", linewidth=1) + theme_bw() # asympt diff ratio resHistAsDiff <- as.data.frame(resHist[,11:13]) resHistAsDiff_l <- melt(data.frame(resHistAsDiff), variable.name="ratio_CIs", value.name = "readout") lowquant_AsDiff <- quantile(resHistAsDiff[,2], probs=as.numeric(input$lowQuant)/100) upquant_AsDiff <- quantile(resHistAsDiff[,3], probs=as.numeric(input$uppQuant)/100) p_AsDiff <- ggplot(resHistAsDiff_l, aes(readout, fill=ratio_CIs)) + geom_histogram(alpha=0.5,position="identity") + labs(title = paste("asymptote diff. ratio EACs:", round(lowquant_AsDiff,3), " to ", round(upquant_AsDiff,3))) + geom_vline(xintercept = c(lowquant_AsDiff, upquant_AsDiff), color="black", linetype="dashed", linewidth=1) + geom_vline(xintercept = c(input$lEACratioAdiff , input$uEACratioAdiff), color="red", linetype="dashed", linewidth=1) + theme_bw() # Slope ratio resHistSlo <- as.data.frame(resHist[,5:7]) resHistSlo_l <- melt(data.frame(resHistSlo), variable.name="ratio_CIs", value.name = "readout") lowquant_Slo <- quantile(resHistSlo[,2], probs=as.numeric(input$lowQuant)/100) upquant_Slo <- quantile(resHistSlo[,3], probs=as.numeric(input$uppQuant)/100) p_Slo <- ggplot(resHistSlo_l, aes(readout, fill=ratio_CIs)) + geom_histogram(alpha=0.5,position="identity") + labs(title = paste("Slope ratio EACs:", round(lowquant_Slo,3), " to ", round(upquant_Slo,3))) + geom_vline(xintercept = c(lowquant_Slo, upquant_Slo), color="black", linetype="dashed", linewidth=1) + geom_vline(xintercept = c(input$lEACratioSlope , input$uEACratioSlope), color="red", linetype="dashed", linewidth=1) + theme_bw() # poency ratio resHistPot <- as.data.frame(resHist[,1:3]) resHistPot_l <- melt(data.frame(resHistPot), variable.name="ratio_CIs", value.name = "readout") lowquant_Pot <- quantile(resHistPot[,2], probs=as.numeric(input$lowQuant)/100) upquant_Pot <- quantile(resHistPot[,3], probs=as.numeric(input$uppQuant)/100) #browser() p_Pot <- ggplot(resHistPot_l, aes(readout, fill=ratio_CIs)) + geom_histogram(alpha=0.5,position="identity") + labs(title = paste("Poency ratio EACs:", round(lowquant_Pot,3), " to ", round(upquant_Pot,3))) + geom_vline(xintercept = c(lowquant_Pot, upquant_Pot), color="black", linetype="dashed", linewidth=1) + geom_vline(xintercept = c(input$lowerPot/100, input$upperPot/100), color="red", linetype="dashed", linewidth=1) + theme_bw() grid.arrange(p_Slo, p_AsDiff, p_uAs, p_Pot, nrow=1) } }) #### download XL report---- output$downloadXLReport <- downloadHandler( filename= paste0("Report_4PLEvaluation", Dat$FileName,".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,"logo.png") file.copy("logo.png", tempReportc, overwrite = T) rmarkdown::render(tempReport, output_file = file, params = list(FileName = Dat$FileName, author = Dat$author, REP = REP, coeffs = Dat$coeffs_UN), envir = new.env(parent = globalenv())) } ) } shinyApp(ui, server)