library(shiny) library(shinyjs) library(shinyAce) 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) library(tinytex) Dat <- reactiveValues() REP <- reactiveValues() #### Functions ---- dilFUN2 <- function(cs_,dils,Faktor) { av <- cs_ dils_av <- dils_av dils_avsc <- dils_av*Faktor dils2 <- dils_avsc+av dilfactors <- 1/exp(dils2-lag(dils2)) return(dilfactors) } plot_f <- function(dat, 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() if(is.null(det_sig)) { if (CORdat<0) { startlist <- list(a=sigmoid[1], b=-sigmoid[5],cs=sigmoid[7], d=sigmoid[3],r=sigmoid[8]) } else { startlist <- list(a=sigmoid[1],b=sigmoid[5],cs=sigmoid[7], d=sigmoid[3],r=sigmoid[8]) } } else { startlist <- list(a=det_sig[5], b=det_sig[1],cs=det_sig[7], d=det_sig[3],r=det_sig[7] - det_sig[8]) } mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(log_dose-(cs-r*isSample)))), data=all_l, start=startlist, 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 }) a <- s_mr$coefficients[1,1] b <- s_mr$coefficients[2,1] cs <- s_mr$coefficients[3,1] d <- s_mr$coefficients[4,1] r <- s_mr$coefficients[5,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*(seq_x-(cs-r)))) REF <- a+(d-a)/(1+exp(b*(seq_x-(cs)))) if (is.null(det_sig)) { SAMPLEtrue <- sigmoid[2] + (sigmoid[4] -sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-sigmoid[8]-seq_x)))) REFtrue <- sigmoid[1] + (sigmoid[3] -sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x)))) } else { SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(det_sig[2]*(det_sig[8]-seq_x))) REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-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*(a-d)/4 Xbendl3 <- cs-(1.31696/b) Xbendu3 <- cs+(1.31696/b) XbendlT <- cs-r-(1.31696/b) XbenduT <- cs-r+(1.31696/b) bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3), bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,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), alpha=0.8) + labs(title = paste("restricted 4pl; bendp:", round(Xbendl3,3),round(Xbendu3,3),round(XbendlT,3),round(XbenduT,3)), color="product") + scale_color_manual(labels=c("test","reference"), values=c("red","blue")) + 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="red", inherit.aes = F) + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="blue", inherit.aes = F) + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="red", linetype=2, alpha=0.4, inherit.aes = F) + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="blue", linetype=2, alpha=0.4, inherit.aes = F) + geom_vline(xintercept=c(Xbendl3, Xbendu3), col="blue",linetype=2) + geom_vline(xintercept=c(XbendlT, XbenduT), col="red",linetype=2) + annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) + 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), alpha=0.8) + labs(title = paste("restricted transformed 4pl"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("red","blue")) + theme_bw() mrt <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(log_dose-(cs-r*isSample)))), data=all_l2, start=startlist, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) s_mrt <- summary(mrt) a_trans <- s_mrt$coefficients[1,1] b_trans <- s_mrt$coefficients[2,1] cs_trans <- s_mrt$coefficients[3,1] d_trans <- s_mrt$coefficients[4,1] r_trans <- s_mrt$coefficients[5,1] XbendlTrans <- cs_trans-(1.31696/b_trans) XbenduTrans <- cs_trans+(1.31696/b_trans) XbendlTransT <- cs_trans-r_trans-(1.31696/b_trans) XbenduTransT <- cs_trans-r_trans+(1.31696/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*(seq_x-(cs_trans-r_trans)))) REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*(seq_x-(cs_trans)))) 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="red", inherit.aes = F) + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="blue", inherit.aes = F) + geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="blue",linetype=2) + geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="red",linetype=2) + theme(legend.position = "none", axis.text=element_text(size=14)) if (is.null(det_sig)) { unrestr <- drm(readout ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(), pmodels=data.frame(isSample, isSample,isSample,isSample)) Sum_u <- summary(unrestr) ast <- Sum_u$coefficients[3,1] ate <- Sum_u$coefficients[4,1] bst <- Sum_u$coefficients[1,1] bte <- Sum_u$coefficients[2,1] cst <- log(Sum_u$coefficients[7,1]) cte <- log(Sum_u$coefficients[8,1]) dst <- Sum_u$coefficients[5,1] dte <- Sum_u$coefficients[6,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*(seq_x-cst))) SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(seq_x-cte))) pl_df2 <- cbind(seq_x, SAMPLEu, REFu) #browser() pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + geom_point() + labs(title="unrestricted 4_pl-Model", color="product") + scale_color_manual(labels = c("test","reference"), values=c("red","blue")) + theme_bw() pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu), color="red", inherit.aes = F) + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu), color="blue", 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() + labs(title="unrestricted transformed 4_pl-Model", color="product") + scale_color_manual(labels = c("test","reference"), values=c("red","blue")) + theme_bw() unrestr_trans <- drm(readouttrans ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(), pmodels=data.frame(isSample, isSample,isSample,isSample)) Sum_ut <- summary(unrestr_trans) ast_t <- Sum_ut$coefficients[3,1] ate_t <- Sum_ut$coefficients[4,1] bst_t <- Sum_ut$coefficients[1,1] bte_t <- Sum_ut$coefficients[2,1] cst_t <- log(Sum_ut$coefficients[7,1]) cte_t <- log(Sum_ut$coefficients[8,1]) dst_t <- Sum_ut$coefficients[5,1] dte_t <- Sum_ut$coefficients[6,1] REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(seq_x-cst_t))) SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(seq_x-cte_t))) 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="red", inherit.aes = F) + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans), color="blue", inherit.aes = F, show.legend = F) pu3_t <- pu2_t grid.arrange(p2,p_rt2,pu2_,pu3_t, nrow=2) } 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) } 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,0,1), ifelse(pFR2_A<0.05,0,1),ifelse(pFR2_B<0.05,0,1), ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1), ifelse(p_F_nonp<0.05,0,1),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) } 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) CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) if (CORdat<0) SLOPE <- -1 else SLOPE <- 1 startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0) tryCatch({ mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), data=all_l, start=startlist, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { err$message }) startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose), at=max(ro_new[,2]), bt=SLOPE, dt=min(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(err) { err$message }) if (!PureErrFlag) { pot_est <- exp(confintd(mr, "r", method="asymptotic")) potU_est <- exp(confintd(mu, "r", method="asymptotic")) } 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() }) #browser() if (length(SU_mr)>1) { s_mr <- SU_mr$coefficients } else { SU_mr <- rep(NA,5) } VCOV <- vcov(mr) V_V <- VCOV/SU_mr$sigma^2 SEsPure <- sqrt(diag(V_V)*meanPureErr) pot_est <- c(exp(s_mr[5,1]),exp(s_mr[5,1]-qt(0.975,nrow(all_l)-5)*SEsPure[5]), exp(s_mr[5,1]+qt(0.975,nrow(all_l)-5)*SEsPure[5])) # unrestricted s_mu <- summary(mu)$coefficients SU_mu <- summary(mu) VCOVu <- vcov(mu) V_Vu <- VCOVu/SU_mu$sigma^2 SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) potU_est <- c(exp(s_mu[7,1]),exp(s_mu[7,1]-qt(0.975,nrow(all_l)-8)*SEsPureU[7]), exp(s_mu[7,1]+qt(0.975,nrow(all_l)-8)*SEsPureU[7])) } # PureErrFlag startlistmr_log <- list(a=max(all_l$readouttrans), b=SLOPE, d=min(all_l$readouttrans), cs=mean(all_l$log_dose),r=0) tryCatch({ mr_log <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), data=all_l, start=startlistmr_log, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { err$message }) startlistmu_log <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose), at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0) tryCatch({ mu_log <- gsl_nls(fn = readouttrans ~ 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_log, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { err$message }) pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic")) potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic")) colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI") colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI") su_mr_log <- summary(mr_log) Dat$RMSE_Rlog <- su_mr_log$sigma su_mu_log <- 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(pot_est, potU_est, pot_est_log, potU_est_log) potALL2 <- cbind(c("restricted","unrestricted","transformed restr","untransf restr"), potALL) return(potALL2) } ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVarlog,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(bs) var_log_xt <- (se_xt/xt)^2 se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVarlog) 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) } 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 pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(1,1,1,isSample)) compParm(pot, "c",display=T) ED50 <- ED(pot,c(50), interval="delta") PotEst <- ED50[1,1]/ED50[2,1] potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE) potAll2 <- potAll[1:3] CORro <- cor(ro_new[,1], ro_new[,ncol(ro_new)]) if (CORro<0) SLOPE <- -1 else SLOPE <- 1 startlist <- list(a=max(ro_new[,2]), b=SLOPE, d=min(ro_new[,2]), cs=mean(all_l$log_dose),r=0) tryCatch({ mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), data=all_l, start=startlist, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { err$message }) startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose), at=max(ro_new[,2]), bt=SLOPE, dt=min(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(err) { err$message }) smu <- tryCatch({ summary(mu) }, error=function(err) { return(0) }) POTr_CI <- potAll2[2:3] 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 <- vcovMU/smu$sigma^2 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(pot$predres[,2]^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) #browser() ## EQ test on lower As diff ds <- coeffs["ds"] dt <- coeffs["dt"] lAs_diff <- (dt-ds) uCI_laDiff <- lAs_diff+qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,2]^2) lCI_laDiff <- lAs_diff-qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,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)-noConc uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.9975) 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 <- slopeCI #### 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 <- (at-dt)/(as-ds) at_dt <- (at-dt) as_ds <- (as-ds) 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(at_dt,as_ds,se_dt_at, se_ds_as,CoVarlog=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 <- (as-ds) 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 significance of lower asymptotes", "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(at_dt/as_ds, 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) } ANOVA4plUnresfunc <- function(ro_new, sigmoid) { 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 pot <- drm(readout ~ Conc, isSample, data=all_l2, 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((potU$predres[,1] - mean(all_l$readout))^2),5) SStreat_df <- length(unique(all_l$log_dose))-1 SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5) ## Non-parallel SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5) ## Preparation SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) ## Resid Err RSS <- round(sum(potU$predres[,2]^2),5) RSS_df <- nrow(all_l)-SStreat_df-1 FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) # PureErr SSE <- FitAnova[4,3] SSE_df <- FitAnova[4,1] # Non-Linearity SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^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) } 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) } 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) } #_______________________________________________________________________________ #### Main function ---- #_______________________________________________________________________________ function(input,output, session) { 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 = cnSheets2) } }) 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") { 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 ec50est <- (max(all_l$log_dose)+min(all_l$log_dose))/2 startlist <- list(a=min(all_l$readout), b=SLOPE, d=max(all_l$readout), cs=ec50est,r=0) tryCatch({ mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), data=all_l, start=startlist, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { err$message }) startlistmu <- list(as=min(all_l$readout), bs=SLOPE, ds=max(all_l$readout), cs=ec50est, at=min(all_l$readout), bt=SLOPE, dt=max(all_l$readout), 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(err) { err$message }) Smr <- summary(mr) Smu <- summary(mu) coeffsMR <- Smr$coefficients[,1] coeffsMU <- Smu$coefficients[,1] Dat$coeffsMRes <- coeffsMR names(coeffsMU) <- c("lowAsym REF", "slope REF","upperAsym REF","EC50 REF","lowAsym TEST","slope TEST","upperAsym TEST","r") if (!PureErrFlag) { pot_est <- exp(confintd(mr, "r", method="asymptotic")) potU_est <- exp(confintd(mu, "r", method="asymptotic")) 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 <- 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 <- VCOVu/Smu$sigma^2 #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])) } 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)) SR <- summary(pot) SU <- summary(potU) coeffs_UN <- potU$coefficients coeffs_UN[1] <- ifelse(xor(CORro>0, coeffs_UN[1]>0), -coeffs_UN[1],coeffs_UN[1]) coeffs_UN[2] <- ifelse(xor(CORro>0, coeffs_UN[2]>0), -coeffs_UN[2],coeffs_UN[2]) coeffs_UN[7:8] <- log(coeffs_UN[7:8]) POTU <- EDcomp(potU, percVec = c(50,50), interval="delta",display=F) Dat$potDiffXL <- POTU[1]*100 RMSE_unr_diagn <- sqrt(SU$resVar) RMSE_res_diagn <- sqrt(SR$resVar) up_lowDiffDiagn <- SU$coefficients[5,1] - SU$coefficients[3,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 <- max(input$lowerPot, pot_est3$lowerCI) MaxPl_ <- MaxPl*1.1 MinPl_ <- MinPl*0.9 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$lupperPot), col="indianred") + annotate("text", x=pot_est3$lowerPot-10, y=1.020, label=paste("lower CL:", round(pot_est3$lowerCI,1)), col="darkgreen") + annotate("text", x=pot_est3$upperPot+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 }) SStreat <- round(sum((potU$predres[,1] - mean(all_l$readout))^2),5) SStreat_df <- length(unique(all_l$log_dose))-1 SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5) ## Non-parallel SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5) ## Preparation SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) ## Resid Err RSS <- round(sum(potU$predres[,2]^2),5) RSS_df <- nrow(all_l)-SStreat_df-1 FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) # PureErr SSE <- FitAnova[4,3] SSE_df <- FitAnova[4,1] # Non-Linearity SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^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) ANOVAtab2 <- data.frame(Source = c("Treatment","Preparation","Regression", "Non-Parallelism","Residual Error","Non-linearity", "Pure Error","Total"), DF = round(AnovaDFs,0), 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,"","") ) output$ANOVAXLS <- renderTable({ ANOVAtab2 }) REP$ANOVAXLS <- ANOVAtab2 #browser() startlistlog <- list(a=min(log(all_l$readout)), b=SLOPE, d=max(log(all_l$readout)), cs=ec50est,r=0) tryCatch({ mrlog <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))), data=all_l, start=startlistlog, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { print("Error in mrlog gsl_nls") }) startlistmulog <- list(as=min(log(all_l$readout)), bs=SLOPE, ds=max(log(all_l$readout)), cs=ec50est, at=min(log(all_l$readout)), bt=SLOPE, dt=max(log(all_l$readout)), r=0) tryCatch({ mulog <- 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=startlistmulog, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, error = function(err) { print("Error in murlog gsl_nls") }) logpot <- drm(log(readout) ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(1,1,1,isSample)) logpotU <- drm(log(readout) ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), pmodels=data.frame(isSample, isSample,isSample,isSample)) Smrlog <- summary(mrlog)$coefficients Smulog <- summary(mulog)$coefficients SUlog <- summary(logpotU) SRlog <- summary(logpot) RMSE_unrlog_diagn <- sqrt(SUlog$resVar) RMSE_reslog_diagn <- sqrt(SRlog$resVar) up_lowDifflogDiagn <- SUlog$coefficients[5, 1] - SUlog$coefficients[3, 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))) DatSDiagnTable <- DiagnTable REPSDiagnTable <- DiagnTable logpotest <- exp(confintd(mrlog, "r", method = "asymptotic")) # compParm(logpot, "c") logpotuest <- 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") cnXL <- colnames(XLdat2) Filesample <- data.frame(c("File name", "samples"), c(Dat$Filename, paste(cnXL[1], " vs ", cnXL[4], "\n"))) colnames(Filesample) <- c("", "") outputIFilesampl <- renderTable({ Filesample }) 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, 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") outputcoeffs_unr2 <- renderTable({ UnRRPLAausw2 }) 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") # coeffs_R[4] <- log(coeffs_R[4]) # coeffs_R[5] <- log(coeffs_R[5]) # --- Ergebnistabelle: RESTRICTED (Eingeschränktes Modell) --- 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$coeffs_r2 <- renderTable({ PLAAusw2 }, digits = 3, rownames = T) REP$PLAausw <- PLAAusw REP$PLBend <- PLAAusw2 # --- Koeffizienten-Extraktion --- logcoeffs_R <- Smrlog[, 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 gs1_nls output$logcoeffs_r <- renderTable({ LogPLAAusw }) REP$LogPLAausw <- LogPLAAusw logcoeffs_UNR <- Smulog[,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(logpotest * 100, 3)))) # von gs1_nls output$logcoeffs_unr <- renderTable({ LogUnrPLAAusw }) REP$LogUnrPLAausw <- LogUnrPLAAusw #browser() Dat$coeffs_UN <- coeffs_UN if (exists("Ind")) { Dat$dilution <- XLdat[,Ind] } else Dat$dilution <- exp(XLdat[,logI]) # --- Plot-Ausgabe --- output$XLplot <- renderPlot({ plot_f(XLdat2, sigmoid = NULL, det_sig = coeffs_UN) }) REP$XLdat2 <- XLdat2 # --- Diagnose-Plots (Residualanalyse) --- output$diagnplot <- renderPlot({ op <- par(mfrow = c(1, 2), mar = c(3.2, 3.2, 2, .5), mgp = c(2, .7, 0)) # 1. Residuals vs Fitted plot(residuals(pot) ~ fitted(pot), main = "Residuals restricted") abline(h = 0) qqnorm(residuals(pot)) qqline(residuals(pot)) plot(residuals(potU) ~ fitted(potU), main = "Residuals unrestricted") abline(h = 0) qqnorm(residuals(potU)) qqline(residuals(potU)) par(op) # Parameter zurücksetzen }) 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 } # inout$iFile }) #### make geomDils reactive ---- observe({ #browser() 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, "sdfacf", 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({ 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$plot <- renderPlot({ #browser() sigmoid <- sigmoid() det_sig=NULL plot_f(sim2(),sigmoid, det_sig) }) #### Testergebnisse ---- observe({ #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 tab <- tests_FUNC(sim2(), Limite, PureErrFlag = PureErrFlag) 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'))) #browser() output$secondary <- DT::renderDataTable({ dat}) output$EQtests <- DT::renderDataTable( dat, options = list(lengthChange=FALSE) ) }) # 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" )) }) 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) }) #### linear Plot output ---- output$plotLin <- renderPlot({ tab <- sim2() if (is.character(tab)) stop(tab) if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer = (ncol(tab)-1)/2 Conc <- CONC() Conctab <- Dat$Conctab 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 sigmoid <- sigmoid() log_dose <- unique(all_l$log_dose) seq_x <- seq(min(log_dose), max(log_dose),0.1) SAMPLEtrue <- sigmoid[2] + (sigmoid[4]-sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-log(input$potencydiff/100)-seq_x)))) REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-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("red","blue")) + ylim(min(all_l2$readout),max(all_l2$readout)) + theme_bw() p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="blue", inherit.aes = F) + geom_line(data=pl_df,aes(x=lnC,y=plotT),color="red", inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4, inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", 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)) 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) pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log(Conc) pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log(Conc) pl_rest <- data.frame(lnC=log(Conc), plotS=pl_restS, plotT=pl_restT) pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="blue", inherit.aes = F) + geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="red", inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4, inherit.aes = F) + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", 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)) grid.arrange(p3,pr3,nrow=1) }) #### linear PLA tests ---- output$TESTSlin <- DT::renderDataTable({ tab <- sim2() if (is.character(tab)) stop(tab) Conc <- CONC() 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)) circles <- Dat$circles PureErrFlag <- input$PureErr warning_text <- reactive({ ifelse(PureErrFlag, 'Pure error is selected','') }) output$PureErrW <- renderText(warning_text()) LIN <- ANOVAlintests(tab,circles,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) #browser() Dat$ANOVA <- df[,4:length(df)] dat <- datatable(df[,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 --- output$ANOVA <- DT::renderDataTable({ sigmoid <- sigmoid() tab <- ANOVA4plUnresfunc(sim2(),sigmid) dat <- datatable(df[,1:3], 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, "(~ entered % upper-lower asymptote)\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({ 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 pot tab ---- output$pottab <- 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$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"))) }) #### 4pl potency table ---- output$pottab4pl <- DT::renderDataTable({ ro_new <- sim2() Dils_ <- Dils() if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer <- 3 PureErrFl <- input$PureErr pottab4 <- pot4plFUNC(ro_new = ro_new, PureErrFlag = PureErrFl) 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)) if (as.numeric(pottab4[1,3])*100>Lim[[9]] & as.numeric(pottab4[1,4])*100Lim[[9]] & as.numeric(pottab4[2,4])*100Lim[[9]] & as.numeric(pottab4[3,4])*100Lim[[9]] & as.numeric(pottab4[4,4])*100% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1), c("lightgreen","pink"))) }) #### Dilutions Simulator ---- output$plotfordilutions <- renderPlot({ tab <- sim() 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.31696/sigmoid[5] Xbend200l <- sigmoid[7] - 0.693147-1.31696/sigmoid[5] Xbend50u <- sigmoid[7] + 0.693147+1.31696/sigmoid[5] Xbend200u <- sigmoid[7] - 0.693147+1.31696/sigmoid[5] Xbend50 <- max(Xbend50l, Xbend50u) Xbend200 <- min(Xbend200l, Xbend200u) dummy <- plot_f(tab,sigmoid,det_sig=NULL) } 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.31696/det_sig[1] Xbend200l <- det_sig[7] - 0.693147-1.31696/det_sig[1] Xbend50u <- det_sig[7] + 0.693147+1.31696/det_sig[1] Xbend200u <- det_sig[7] - 0.693147+1.31696/det_sig[1] Xbend50 <- max(Xbend50l, Xbend50u) Xbend200 <- min(Xbend200l, Xbend200u) dummy <- plot_f(tab,sigmoid=NULL,det_sig=det_sig) } 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({ 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())) } ) # observe({ # pfad <- getwd() # Z <- read.csv('./Zaehlfile.csv') # Z[,3] <- sapply(Z[,3], as.character) # nrowsZ <- nrow(Z) # # Z[nrowsZ+1,1] <- Z[nrowsZ,1]+1 # Z[nrowsZ+1,1] <- Sys.Date() # Z[nrowsZ+1,1] <- NA # write.csv(Z, file= './Zaahlfile.csv', row.names = F) # }) # observe({ # Z <- read.csv('./Zaehlfile.csv') # Dat$Zaehl <- read.csv('./Zaehlfile.csv') # nrowsZ <- nrow(Z) # message2 <- paste("No. of calls:", Z[nrowsZ,1]) # output$stats <- renderPrint(cat(message2)) # }) # output$downloadZaehl <- downloadHandler( # filename = function() { # paste("Zaehlfile", Sys.Date(),".csv", sep="") # # }, # content = function(file) { write.csv(Dat$Zaehl, file = file, row.names = F)}, # contentType = "text/csv" # ) } # end of session function