commit d02cea0698773ae3e58305a8acd0c6eec82772b4 Author: Franz Innerbichler Date: Tue Apr 14 21:46:10 2026 +0200 initial commit diff --git a/Doc_BioassayReport.Rmd b/Doc_BioassayReport.Rmd new file mode 100644 index 0000000..769170f --- /dev/null +++ b/Doc_BioassayReport.Rmd @@ -0,0 +1,433 @@ +--- +output: + pdf_document: + extra_dependencies: ["float"] + number_sections: true + toc: true + toc_depth: 3 +header_includes: + -\usepackage{fancyheadr} + -\setlength{\headheight}{22pt}% + -\usepackage{lastpage} + -\pagestyle{fancy} + -\usepackage{pdflscape} + -\usepackage{longtable} + -\rhead{\includegraphics[width=.15\textwidth]{`r getwd()`/logo.png}} +params: + FileName: NA + newTitle: NA + author: NA + REP: NA + coeffs: NA +author: "Author: `r params$author`" +title: | + | ![](logo.png){width=1in} + | 4PL bioassay evaluation +subtitle: | + `r params$FileName` + + Unique time: `r Sys.time()` +date: "`r paste(params$Subway, params$Version)`" + +--- + +\fancyfoot[C]{\thepage\ of \pageref{LastPage}} +\newpage + +\newpage + +```{r setup, include=FALSE} + +knitr::opts_chunk$set(echo = TRUE) + +library(knitr) +library(DT) + +REP <- params$REP +author <- params$author +coeffs <- params$coeffs + +all_l <- REP$all_l +ANOVAXLS <- REP$ANOVAXLS +DiagnTable <- REP$DiagnTable +UnRPLAausw <- REP$UnRPLAausw +UnRPLBend <- REP$UnRPLBend +PLAausw <- REP$PLAausw +PLBend <- REP$PLBend +LogPLAausw <- REP$LogPLAausw +LogUnrPLAausw <- REP$LogUnrPLAausw + +XLdat2 <- REP$XLdat2 + +CIplot <- REP$CIplot +testsTab <- REP$testsTab +relpotTestPlot <- REP$relpotTestPlot + + + +``` + + +# Introduction + +Bioassay potency estimation uses statistical methods to quantify the strength of a biological product or drug by comparing its response to that of a reference standard. Because biological responses are inherently variable, affected by assay conditions, cell systems or organisms, and measurement noise, the 4-parametric logistic regression is used to obtain reliable potency values. The variance for confidence interval calculation is coming from the regression procedure itself and is an excellent predictor for the variability of any future potency determinations. +USP<1034> recommends calculation of standard errors of ratios of the parameters using Fieller's theorem [Finney D.J. 1978] or using the "delta" method (for a discussion about the "delta" method see [Ver Hoef 2012]). However, the presented gradient approach using the differences on the log-scale is methematically more stable und thus preferable compared to any ratio approach ([Franz, V.H. 2007]). + +# Results + +All data used for the 4PL evaluation is shown in table 1: + +```{r alll, echo=FALSE, warning=FALSE, results='asis'} + +kable(all_l, format = "markdown", caption= "Uploaded data (test and reference) in long format", digits=3) + +``` + +The following 4 plots show all 4 models: restricted and unrestricted, and log transformed, respectively. + +You can also embed plots, for example: + +```{r XLplot, echo=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap="Plot of models", fig.align='left'} + +plot_f <- function(dat, sigmoid,det_sig) { + CORdat <- cor(dat[,1],dat[,ncol(dat)]) + + 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) + + if(is.null(det_sig)) { + if (CORdat<0) { + startlist <- list(a=sigmoid[3], b=-sigmoid[5],cs=sigmoid[7], + d=sigmoid[1],r=sigmoid[8]) + } else { + startlist <- list(a=sigmoid[3],b=sigmoid[5],cs=sigmoid[7], + d=sigmoid[1],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]) + } + #browser() + tryCatch({ + mr <- gsl_nls(fn = readout ~ 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)) + }, + error = function(err) { + err$message + }) + s_mr <- summary(mr) + 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[4] + (sigmoid[2] -sigmoid[4])/(1+exp(sigmoid[6]*(seq_x-(sigmoid[7]-sigmoid[8])))) + REFtrue <- sigmoid[3] + (sigmoid[1] -sigmoid[3])/(1+exp(sigmoid[5]*(seq_x-(sigmoid[7])))) + } else { + SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(-det_sig[2]*(seq_x-(det_sig[8])))) + REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(-det_sig[1]*(seq_x-(det_sig[7])))) + } + + 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)) + + 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") + + + # 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)) + + 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) +} + +plot_f(XLdat2, sigmoid=NULL, det_sig=coeffs) + + +``` + +The ANOVA of the unconstrained model is listed in table 2: + +```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'} + +kable(ANOVAXLS, format = "markdown", caption= "ANOVA table unrestricted", digits=3) + + +``` + + +```{r SST_ergebn, fig.align='center', fig.pos='htb!', echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE} + +kable(testsTab[1:7,], row.names = F, format = "markdown", caption="SST results") + + +``` + + +*...The estimate for F-test on regression and on non-linearity is the p-value + F-test on regression passes if F-value > F-crit and thus p < 0.05 + F-test on non-linearity passes if F-value < F-crit and thus p > 0.05 + 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 + + +```{r, label= 'CIplot', echo=FALSE, warning=FALSE, fig.width=100, fig.cap='Selected SSt confidence intervals with entered limits', fig.align='center'} + +png("CIplot.png") +print(CIplot) +dev.off() + + +``` + + +![](CIplot.png){width=60%} + + +## Fitting results of the 4 models with bend points + +The results of the non-linear fitting procedure for the restricted model (5 parameters) is listed in table 4: + +```{r PLAausw, echo=FALSE, warning=FALSE, results='asis'} + +kable(PLAausw, format = "markdown", caption= "Restricted 4PL evaluation", digits=3, row.names = F) + + +``` + + +A depiction of the CI and corresponding limits of relative potency is shown here: + +```{r, label='relpotPlot', echo=FALSE, warning=FALSE, fig.height=2, fig.width=3.5, fig.cap="Rel potency with CIs and limits", fig.align='left', results='asis'} + +print(relpotTestPlot) + + +``` + + +The bend points for test and reference sample are in table 5: + +```{r PLBend, echo=FALSE, warning=FALSE, results='asis'} + +kable(PLBend, format = "markdown", caption= "Bendpoints (Sebaugh) of restricted 4PL", digits=3) + + +``` + +The results of the non-linear fitting procedure for the unrestricted model (8 parameters) is listed in table 6: + +```{r UnRPLAausw, echo=FALSE, warning=FALSE, results='asis'} + +kable(UnRPLAausw, format = "markdown", caption= "UNrestricted 4PL evaluation", digits=3, row.names = F) + + +``` + + +```{r UnRPLBend, echo=FALSE, warning=FALSE, results='asis'} + +kable(UnRPLBend, format = "markdown", caption= "Bend points of 4PL unrestricted", digits=3, row.names = F) + + +``` + + + +```{r LogPLAausw, echo=FALSE, warning=FALSE, results='asis'} + +kable(LogPLAausw, format = "markdown", caption= "Restricted 4PL evaluation with log-transformed response", digits=3) + + +``` + + + +```{r LogUnRPLAausw, echo=FALSE, warning=FALSE, results='asis'} + +kable(LogUnrPLAausw, format = "markdown", caption= "Unrestricted 4PL evaluation with log-transformed response", digits=3) + + +``` + + +# Appendix: Formulas + +## 4PL regression + +$$ +Y = D + \frac{A-D} {1+(\frac{C} {x})^B } + \epsilon +$$ + + +## log-logistic 4P regression + +$$ +Y = D + \frac{A-D} {1+e^{(B*(C - log(x))) }} + \epsilon +$$ + +where: x ... concentration of the analyte + +A: upper asymptote + +B: slope + +D: lower asymptote + +C ... EC50 + +# Literature + +Finney, D.J.: (1978) Statistical Method in Biological Assay, London: Charles Griffin House, 3rd edition (pp. 80-82) + +Franz, V.H.: Ratios: A short guide to confidence limits and proper use. arXiv:0710.2024v1, 10 Oct 2007 + +VerHoef, J.M.: Who invented the Delta Method? The American Statistician, 2012, 66:2, 124-127 DOI: 10.1080/00031305.2012.687494 + + + + diff --git a/TestF4PL.xlsx b/TestF4PL.xlsx new file mode 100644 index 0000000..d349bfb Binary files /dev/null and b/TestF4PL.xlsx differ diff --git a/TestFnegSl.xlsx b/TestFnegSl.xlsx new file mode 100644 index 0000000..f8d78b5 Binary files /dev/null and b/TestFnegSl.xlsx differ diff --git a/server.R b/server.R new file mode 100644 index 0000000..33920e8 --- /dev/null +++ b/server.R @@ -0,0 +1,2306 @@ + + +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, circ_ABl)) + 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) + + outputcoeffs_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$diagplot <- 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", input$iFile,".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 + + + + + + + diff --git a/ui.R b/ui.R new file mode 100644 index 0000000..de2b31b --- /dev/null +++ b/ui.R @@ -0,0 +1,276 @@ +library(shiny) +library(shinyjs) +library(shinyAce) +library(bslib) + +version <- "1.0.0" + +function(req) { + fluidPage( + useShinyjs(), + titlePanel("Dilution Wizard"), + tags$head(tags$style((".shiny-notification {position: fixed; top: 20%;left: 50%; margin-top: -100px; + margin-left: -250px; color: blue;font-size: 20px; font-style: italic;}")), + tabsetPanel(id="toptab", + tabPanel("Data entry", + wellPanel( + fluidRow( + tags$head(tags$style(HTML("label {font-size:80%;margin-bottom: 3px;margin-top: 3px;}"))), + column(2, + tags$image(src="logo.png", class="adv_logo"), + h4("upload EXCEL"), + fileInput(inputId = "iFile", label = "", accept="ms-excel"), + uiOutput(outputId = "sheetName"), + div(checkboxInput("PureErr", "Should pure error be used for calculation of CIs?", FALSE), + style = "font-size: 24px !important;color: red"), + verbatimTextOutput("stats"), + h5("\n\n\n Author: Franz Innerbichler, InnerAnalytics")), + div(id="parameter", + column(1,style = "background: lightgrey", + 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(1,style = "background: lightgrey", + 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(1,style = "background: lightgreen", + 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(1,style = "background: lightgreen", + + numericInput("CONC7", "7th concentration",0.00469), + numericInput("CONC8", "8thd concentration",0.00235), + numericInput("CONC9", "9thd concentration",NA), + numericInput("CONC10", "10th concentration",NA), + numericInput("CONC11", "11th concentration",NA), + + numericInput("CONC12", "lowest concentration",NA) + ), + column(1,style = "background: lightblue", + h4("geometric dilution scheme"), + numericInput("ConcStart", "starting concentration",NA), + numericInput("dilutionFac", "dilution factor",NA), + numericInput("NoDil", "no. of dilutions",NA), + numericInput("NoDilSer", "no. of dil. series",NA), + verbatimTextOutput("dilutions") + ), + column(4, + h4("log-dilution scheme"), + verbatimTextOutput("logdil"), + h4("User help"), + h5("If new dilutions are entered, please adjust EC50 to avoid calculation errors")) + ) + ) + )), + + #### XLSX diagnostics ---- + tabPanel ("XLSX diagnostics", + tags$style(HTML("pre { color: black; background-color: #FFE1FF; + font-family: 'Helvetica Neue', Helvetica, Arial, sans-serif;font-size: 12px;} ")), + wellPanel( + fluidRow( + column(2, + 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"), + downloadButton("downloadXLReport", label="Download PDF report", class="butt"), + tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"), + plotOutput("relPotTestPlot", width="300px", height="150px"), + 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(4, + plotOutput("XLplot"), + column(6, + br(), + "Regression results restricted", + tableOutput("coeffs_r"), + "Bend points restricted", + tableOutput("coeffs_r2")), + column(6, + br(), + "Regression results unrestricted", + tableOutput("coeffs_unr"))), + column(4, + plotOutput("diagnplot"), + column(6, + tableOutput("logcoeffs_r"), + tableOutput("coeffs_unr2")), + column(6, + tableOutput("logcoeffs_unr"))), + column(2, + tableOutput("ANOVAXLS"), + DT::renderDataTable("EQtests")) + ) + )), + #### 4pl fits ---- + tabPanel("4pl-fit", + wellPanel( + fluidRow( + column(6, + tags$style(span(htmlOutput("PureErrW3"), style="color: red")), + htmlOutput("PureErrW3"), + tags$head(tags$style("#PureErrW2{color: red; + font-size: 16px; + font_style: italic;}")), + plotOutput("plot"), + DT::dataTableOutput("pottab4pl"), + "Footnote: test performed on absolute CIs."), + column(6, + DT::dataTableOutput("secondary"), + 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%")) + )) + ), + #### Linear PLA ---- + tabPanel("Linear PLA", + wellPanel( + fluidRow( + column(6, + htmlOutput("PureErrW3"), + tags$head(tags$style("#PureErrW2{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")), + column(3, + h3("Tests for linear PLA):"), + DT::dataTableOutput("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")), + column(3, + h3("ANOVA for parallel line assay"), + DT::dataTableOutput("ANOVAlin")) + ) + )), + #### Dilution simulator ---- + tabPanel("Dilution simulator", + wellPanel( + fluidRow( + 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") + ) + ) + )), + #### EQ limits simulation ---- + tabPanel("EQ-limit simulatioin", + wellPanel( + fluidRow( + column(1, style="background: lightblue", + numericInput("simN","number of datasets to simulate", 100, step=10, min=10, max=1000), + numericInput("lowQuant", "lower quantile", 1, min = 0.01,max=20), + numericInput("uppQuant", "upper quantile", 99, min = 801,max=99.99), + actionButton('goSim', "Start simulation",class="btn-success") + ), + column(4, + h4("estimation of EACs for slope ratio, dynamic range-ratio, upper asymptote ratio, and potency"), + plotOutput("plotHistuAs", width = "1400px", height = "400px")) + ) + )), + #### Documentation ---- + tabPanel("Documentation", + 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}$$"))), + + #### setting lower panel ---- + wellPanel( + fluidRow( + 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=1), + sliderInput("outlM","outlier in mid part", min=0, max=1.5,value=0, step=1), + sliderInput("outlU","outlier in upper asymptote", min=0, max=1.5,value=0, step=1) + ), + column(1, + 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) + )), + column(4, + "4 PL ANOVA unrestricted", + DT::renderDataTable("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 MSs column"), + h5("SSE ... 'Pure error' in the SumSquares column"), + h5("RMSE ... Square root of the 'Residual Error' in the MS column"), + verbatimTextOutput("RMSE") + ), + column(5, + DT::dataTableOutput("Conctab")) + ) + ) + + + ) + ) +} \ No newline at end of file