diff --git a/Doc_BioassayLinReport.Rmd b/Doc_BioassayLinReport.Rmd new file mode 100644 index 0000000..6f058cf --- /dev/null +++ b/Doc_BioassayLinReport.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 + REPLin: NA + coeffsLin: 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)`" + +--- + + + + + + +```{r setup, include=FALSE} + +knitr::opts_chunk$set(echo = TRUE) + +library(knitr) +library(DT) + +REP <- params$REP +REPLin <- params$REPLin +author <- params$author +coeffsLin <- params$coeffsLin + +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 + + + + + + + + + + + + + + + +## 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/Global.R b/Global.R new file mode 100644 index 0000000..3118936 --- /dev/null +++ b/Global.R @@ -0,0 +1,1221 @@ + + +Dat <- reactiveValues() +REP <- reactiveValues() + +#' Levenberg Marquard fit of 4 pl +#' +#' Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, +#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. +#' +#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. +#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). +#' @returns Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, +#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. +#' @export +#' @examples +#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' TransF<- FALSE +#' Dat<- list() +#' te <- Fitting_FUNC(dat,TransF) +#' print(te) + + +Fitting_FUNC <- function(ro_new, TransFlag=F) { + CORro <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) + #browser() + all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") + isRef <- rep(c(1,0),1,each=nrow(all_l)/2) + isSample <- rep(c(0,1),1,each=nrow(all_l)/2) + all_l2 <- cbind(all_l, isRef, isSample) + #browser() + if (CORro<0) SLOPE <- -1 else SLOPE <- 1 + if (!TransFlag) { + startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0) + + mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), + data=all_l2, + start=startlist,#trace=T, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + s_mr <- tryCatch({ + s_mr <- summary(mr) + }, + error = function(err) { + s_mr <- NULL + }) + } else { + startlist <- list(a=log(min(ro_new[,2])), b=SLOPE, d=log(max(ro_new[,2])), cs=mean(all_l$log_dose),r=0) + mrT <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), + data=all_l2, + start=startlist, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + s_mr <- summary(mrT) + } + + if (!TransFlag) { + startlistmu <- list(as=min(ro_new[,2]), bs=SLOPE, ds=max(ro_new[,2]), cs=mean(all_l$log_dose), + at=min(ro_new[,2]), bt=SLOPE, dt=max(ro_new[,2]), r=0) + tryCatch({ + mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ + (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), + data=all_l, + start=startlistmu, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + }, + error = function(msg){ + return(0) }) + + Sum_u <- tryCatch({ summary(mu) }, + error=function(msg){ + return(0) }) + } else { + startlistmu <- list(as=log(min(ro_new[,2])), bs=SLOPE, ds=log(max(ro_new[,2])), cs=mean(all_l$log_dose), + at=log(min(ro_new[,2])), bt=SLOPE, dt=log(max(ro_new[,2])), r=0) + tryCatch({ + muT <- gsl_nls(fn = log(readout) ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ + (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), + data=all_l, + start=startlistmu, + control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) + }, + error = function(msg){ + return(0) }) + + Sum_u <- tryCatch({ summary(muT) }, + error=function(msg){ + return(0) }) + } + if (!TransFlag) { + pot_est <- exp(confintd(mr, "r", method="asymptotic")) + potU_est <- exp(confintd(mu, "r", method="asymptotic")) + PRED <- predict(mr) + PREDu <- predict(mu) + } else { + pot_est <- exp(confintd(mrT, "r", method="asymptotic")) + potU_est <- exp(confintd(muT, "r", method="asymptotic")) + PRED <- predict(mrT) + PREDu <- predict(muT) + } + return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu)) +} + + +#' Plot sigmoidal curve +#' +#' Returns the final plots of the 4pl function as sigmoidal lines, and the single readouts as scatter, with REF in blue and TEST in red. +#' Two plots are returned in a grid object: the unrestricted model and the restricted model. +#' +#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. +#' @param sigmoid The parameters of the 4pl fit of unrestricted model from EXCEL input. +#' @param det_sig The parameters of the 4pl fit of unrestricted model from the metadata input. +#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). +#' @returns A grid object either of the original scale or the natural log of the readouts. +#' @export +#' @examples +#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' sigmoid <- c(0.7163324, 0.5636804, 10.6156340 , 9.9784160, -0.7504673, -0.7108692, -3.5788141, -0.6662962) +#' det_sig <- FALSE +#' TransF<- FALSE +#' Dat<- list() +#' p <- plot_f(dat,sigmoid,det_sig,TransF) +#' print(p) + + +plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig, + CORdat <- cor(dat[,1],dat[,ncol(dat)]) + #browser() + all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") + isRef <- rep(c(1,0),1,each=nrow(all_l)/2) + isSample <- rep(c(0,1),1,each=nrow(all_l)/2) + all_l2 <- cbind(all_l, isRef, isSample) + #browser() + MODLS <- Fitting_FUNC(dat, TransFlag=F) + s_mr <- MODLS[[1]] + a <- s_mr$coefficients["a",1] + b <- s_mr$coefficients["b",1] + cs <- s_mr$coefficients["cs",1] + d <- s_mr$coefficients["d",1] + r <- s_mr$coefficients["r",1] + + log_dose <- unique(all_l$log_dose) + seq_x <- seq(min(log_dose),max(log_dose),0.1) + SAMPLE <- a+(d-a)/(1+exp(b*((cs-r)-seq_x))) + REF <- a+(d-a)/(1+exp(b*((cs)-seq_x))) + + s_mu <- MODLS[[2]] + + as <- s_mu$coefficients["as",1] + bs <- s_mu$coefficients["bs",1] + cs <- s_mu$coefficients["cs",1] + ds <- s_mu$coefficients["ds",1] + at <- s_mu$coefficients["at",1] + bt <- s_mu$coefficients["bt",1] + ct <- s_mu$coefficients["cs",1]-s_mu$coefficients["r",1] + dt <- s_mu$coefficients["dt",1] + r <- s_mu$coefficients["r",1] + + + SAMPLEtrue <- at + (dt -at)/(1+exp(bt*((cs-r-seq_x)))) + REFtrue <- as + (ds -as)/(1+exp(bs*((cs-seq_x)))) + + #browser() + pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue) + all_l2$readout[all_l2$readout < 0] <- 0.01 + all_l2$readouttrans <- log(all_l2$readout) + slopeEC50 <- b*(d-a)/4 + Intercept <- a+(d-a)/2-b*(d-a)/4*cs + #browser() + Xbendl3 <- cs-(1.5434/b) + Xbendu3 <- cs+(1.5434/b) + XbendlT <- cs-r-(1.5434/b) + XbenduT <- cs-r+(1.5434/b) + XasymplS <- cs-(3/b) + XasympuS <- cs+(3/b) + XasymplT <- cs-r-(3/b) + XasympuT <- cs-r+(3/b) + bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3), + bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,3), + asympREF_lower = round(XasymplS,3), asympREF_upper=round(XasympuS,3), + asympSAMPLE_lower = round(XasymplT,3), asympSAMPLE_upper=round(XasympuT,3)) + Dat$bendpoints <- bendpoints + Dat$cfordils <- cs + #browser() + p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + + geom_point(shape=factor(isRef), size=3,alpha=0.8) + + labs(title = paste("restricted 4pl"), + color="product") + + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + + scale_shape_manual(labels=c("test","reference")) + + theme_bw() + + theme(axis.text = element_text(size=14)) + + p2 <- p + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLE), color="#C2173F", + inherit.aes = F) + + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="#4545BA", + inherit.aes = F) + + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="#C2173F", linetype=2, alpha=0.4, + inherit.aes = F) + + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="#4545BA", linetype=2, alpha=0.4, + inherit.aes = F) + + geom_vline(xintercept=c(Xbendl3, Xbendu3), col="#4545BA",linetype=2) + + geom_vline(xintercept=c(XbendlT, XbenduT), col="#C2173F",linetype=2) + + geom_vline(xintercept=c( XasymplS, XasympuS), col="#4545BA55",linetype=2) + + geom_vline(xintercept=c(XasymplT, XasympuT), col="#C2173F55",linetype=2) + + annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) + + geom_abline(slope = slopeEC50, intercept = Intercept) + + theme(legend.position="none") + Dat$p2 <- p2 + + # transformed plots + p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + + geom_point(shape=factor(isRef), size=3,alpha=0.8) + + labs(title = paste("restricted transformed 4pl"), color="product") + + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + + theme_bw() + + MODLStrans <- Fitting_FUNC(dat,TransFlag = TRUE) + s_mrt <- MODLStrans[[1]] + a_trans <- s_mrt$coefficients["a",1] + b_trans <- s_mrt$coefficients["b",1] + cs_trans <- s_mrt$coefficients["cs",1] + d_trans <- s_mrt$coefficients["d",1] + r_trans <- s_mrt$coefficients["r",1] + + XbendlTrans <- cs_trans-(1.5434/b_trans) + XbenduTrans <- cs_trans+(1.5434/b_trans) + XbendlTransT <- cs_trans-r_trans-(1.5434/b_trans) + XbenduTransT <- cs_trans-r_trans+(1.5434/b_trans) + bendpointsTRANS <- c(bendREF_lower = round(XbendlTrans,3), bendREF_upper=round(XbenduTrans,3), + bendSAMPLE_lower = round(XbendlTransT,3), bendSAMPLE_upper=round(XbenduTransT,3)) + Dat$bendpointsTRANS <- bendpointsTRANS + SAMPLEtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans-r_trans)-seq_x))) + REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans)-seq_x))) + + pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans) + p_rt2 <- p_rt + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=SAMPLEtrans), color="#C2173F", + inherit.aes = F) + + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="#4545BA", + inherit.aes = F) + + geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="#4545BA",linetype=2) + + geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="#C2173F",linetype=2) + + theme(legend.position = "none", axis.text=element_text(size=14)) + + #browser() + Sum_u <- MODLS[[2]] + #browser() + #if (is.null(det_sig)) { + ast <- Sum_u$coefficients["as",1] + ate <- Sum_u$coefficients["at",1] + bst <- Sum_u$coefficients["bs",1] + bte <- Sum_u$coefficients["bt",1] + cst <- Sum_u$coefficients["cs",1] + cte <- Sum_u$coefficients["cs",1]-Sum_u$coefficients["r",1] + dst <- Sum_u$coefficients["ds",1] + dte <- Sum_u$coefficients["dt",1] + # } else { + # ast <- det_sig[5] + # ate <- det_sig[6] + # bst <- det_sig[1] + # bte <- det_sig[2] + # cst <- det_sig[7] + # cte <- det_sig[8] + # dst <- det_sig[3] + # dte <- det_sig[4] + # } + REFu <- ast + (dst-ast)/(1+exp(bst*(cst-seq_x))) + SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(cte-seq_x))) + pl_df2 <- cbind(seq_x, SAMPLEu, REFu) + #browser() + pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + + geom_point(size=3) + + labs(title="unrestricted 4_pl-Model", color="product") + + scale_color_manual(labels = c("test","reference"), values=c("#C2173F88","#4545BA88")) + + theme_bw() + pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu), + color="#C2173F", inherit.aes = F) + + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu), + color="#4545BA", inherit.aes = F, + show.legend = F) + pu2_ <- pu2 + + theme(legend.position = "none", axis.text = element_text(size=14)) + + putrans <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + + geom_point( size=3) + + labs(title="unrestricted transformed 4_pl-Model", color="product") + + scale_color_manual(labels = c("test","reference"), values=c("#C2173FAA","#4545BAAA")) + + theme_bw() + + Sum_ut <- MODLStrans[[2]] + + ast_t <- Sum_ut$coefficients["as",1] + ate_t <- Sum_ut$coefficients["at",1] + bst_t <- Sum_ut$coefficients["bs",1] + bte_t <- Sum_ut$coefficients["bt",1] + cst_t <- Sum_ut$coefficients["cs",1] + cte_t <- Sum_ut$coefficients["cs",1]-Sum_ut$coefficients["r",1] + dst_t <- Sum_ut$coefficients["ds",1] + dte_t <- Sum_ut$coefficients["dt",1] + + REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(cst_t-seq_x))) + SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(cte_t-seq_x))) + pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans) + + pu2_t <- putrans + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=SAMPLEu_trans), + color="#C2173F", inherit.aes = F) + + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans), + color="#4545BA", inherit.aes = F, + show.legend = F) + pu3_t <- pu2_t + if (TransFlag) grid.arrange(p_rt2,pu3_t, nrow=1) else grid.arrange(p2,pu2_, nrow=1) +} + +#' Simulate a well-plate readout +#' +#' Returns a possible well-plate result based on the metadata provided and a variability. +#' Homo- and heteroscedasticity can be simulated. +#' +#' @param as, bs, cs, ds, at, bt, dt,r,ct, gt, gs are the parameter of the 5pl logistic regression. +#' @param sd_fac The standard deviation with units in % of upper-lower asymptote. +#' @param log_conc The natural log of the concentrations. +#' @param heteroNoise Boolean if heteroscedstic noise should be added. +#' @param noDilSeries A number >1 indicating how many dilution series per product should be simulated. +#' @param noDils A number >7 indicating how many dilutions steps per product should be simulated. +#' @returns A data-frame with readouts and natural log of concentrations. +#' @export +#' @examples +#' as=3; bs=1; cs=-4; ds=10; at=3; bt=1; dt=10;r=0.0001;ct=cs-r;sd_fac=0.1; gt=1; gs=1; +#' lnConc=c(1,0,-1,-2,-3,-4,-5,-6) +#' heteroNoise <- FALSE +#' noDilS <- 3 +#' noD <- 8 +#' Calc_DilRes(as,bs,cs,ds,at,bt,dt,r,ct,sd_fac,gt,gs,log_conc=lnConc, heteroNoise, noDilS, noD) + +Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct=cs-r, + sd_fac=0.1, gt=1, gs=1, log_conc, + heteroNoise=FALSE, noDilSeries, noDils) { + yAxfac <- (ds-as) + log_dose <- log_conc + isRef <- rep(c(1,0),1,each=length(log_conc)*noDilSeries) + isSample <- rep(c(0,1),1,each=length(log_conc)*noDilSeries) + #browser() + av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ + (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose))) + if (heteroNoise) { + # heterosc noise + ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac/100))), nrow=noDils, ncol=noDilSeries*2) + } else { + # homosc noise + ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac*yAxfac/100))), nrow=noDils, ncol=noDilSeries*2) + } + ro_jit <- abs(ro_jit) + + ro_jit2 <- cbind(ro_jit, log_dose) + if (noDilSeries==3) { + colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose") + } else { + colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose") + } + return(ro_jit2) +} + +#' Calculate the potency of linear PLA +#' +#' The delta Method is used for calculating the potency using a restricted linear model. +#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. +#' +#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts. +#' @param Lim The limits to check if the relative CI is within the limits. +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with potency estimate, absolute CIs, test result, relative CIs. +#' @export +#' @examples +#' CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9, +#' -3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7), +#' replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3", +#' "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"), +#' readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2), +#' isRef=c(rep(1,9), rep(0,9)), +#' isSample = c(rep(0,9), rep(1,9))) +#' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant +#' PureErrF <- TRUE +#' +#' +#' LinPotTab(circles=CIRC, Lim, PureErrF) + + +LinPotTab <- function(circles, Lim, PureErrFlag) { + circ_ABl <- circles + circ_Al <- circ_ABl[circ_ABl$isSample ==1,] + circ_Al <- circ_ABl[circ_ABl$isSample ==0,] + # restr CSSI model + modAB <- lm(readout ~ log_dose + isSample, circ_ABl) + coeffs <- modAB$coefficients + SU_modAB <- tryCatch({ + SU_modAB <- summary(modAB) + }, error = function(msg) { + return(NA) + }) + # Intercept diff/slope modAB + linPot <- exp(modAB$coefficients[3]/modAB$coefficients[2]) + + if(PureErrFlag) { + FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl)) + meanPureErr <- FitAnova[4,3] + DFsPure <- FitAnova[4,1] + VCOV <- vcov(modAB) + V_V <- VCOV/SU_modAB$sigma^2 + VCOVpure <- V_V*meanPureErr + SEsPure <- sqrt(diag(V_V)*meanPureErr) + } + + log_pot_delta <- deltaMethod(modAB, "isSample/log_dose") + if (PureErrFlag) { + V_ <- log_pot_delta$SE^2/SU_modAB$sigma^2 + V_p <- V_*meanPureErr + potDeltaPureSE <- sqrt(V_p) + CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure)*potDeltaPureSE + CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure)*potDeltaPureSE + } else { + CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE + CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE + } + #browser() + ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up)) + + # Rel pot CI + relLinpotCI <- ExpLinPot/linPot*100 + if (relLinpotCI[2]>Lim[[9]] & relLinpotCI[3]Lim[9] & ExpLinPot[3]*100>Lim[10]) test_potCI <- 0 else test_potCI <- 1 + + su_mod <- summary(modAB)$coefficients + su_mod2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.")), su_mod) + su_modU <- summary(modABu)$coefficients + su_modU2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.","slope difference")), su_modU) + + uCI_SloDiff <- su_modU[4,1] + qt(0.975,8)*su_modU[4,2] + lCI_SloDiff <- su_modU[4,1] - qt(0.975,8)*su_modU[4,2] + SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiff,uCI_SloDiff) + + lenCirc <- nrow(circ_ABl) + dfTreat <- 3 + dfPrep <- 1 + dfReg <- 1 + dfnonP <- 1 + dfRMSE <- c(lenCirc-3-1) + dfTotal <- lenCirc-1 + dfPureE <- lenCirc-6 + dfNonLin <- dfRMSE-dfPureE + + RSS <- sum(resid(lm(readout ~ log_dose*isSample, circ_ABl))^2) + MSE <- RSS/dfRMSE + SSE <- sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2) + MSpure <- SSE/dfPureE + + if (PureErrFlag) { + FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl)) + meanPureErr <- FitAnova[4,3] + SU_modAB <- tryCatch({ + SU_modAB <- summary(modAB) + }, error = function(msg) { + return(NA) + }) + if (length(SU_modAB)>1) s_modABcoeffs <- summary(modAB)$coefficients + + DFsPure <- FitAnova[4,1] + VCOV <- vcov(modAB) + V_V <- VCOV/SU_modAB$sigma^2 + VCOVpure <- V_V*meanPureErr + SEsPure <- sqrt(diag(V_V)*meanPureErr) + su_mod2[,3] <- SEsPure + su_mod2[,4] <- su_mod2[,2]/su_mod2[,3] + su_mod2[,5] <- 2*(1-pt(abs(su_mod[,4]), FitAnova[4,1])) + + s_mu <- summary(modABu)$coefficients + SU_modABu <- summary(modABu) + VCOVu <- vcov(modABu) + V_Vu <- VCOVu/SU_modABu$sigma^2 + SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) + + su_modU2[,3] <- SEsPureU + su_modU2[,4] <- su_modU2[,2]/su_modU2[,3] + su_modU2[,5] <- 2*(1-pt(abs(su_modU2[,4]), FitAnova[4,1])) + + uCI_SloDiffP <- su_modU[4,1] + qt(0.975,8)*SEsPureU[4] + lCI_SloDiffP <- su_modU[4,1] - qt(0.975,8)*SEsPureU[4] + SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiffP,uCI_SloDiffP) + + SSRes <- SSE + dfRes <- dfPureE + + } else { + SSRes <- RSS + dfRes <- dfRMSE + } + + # treatment + SStreat <- print(sum((predict(lm(readout ~ factor(log_dose)*isSample, circ_ABl))-mean(circ_ABl$readout))^2)) + F_treat <- (SStreat/dfTreat)/(SSRes/dfRes) + # Preparation + SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl))-mean(circ_ABl$readout))^2)) + F_prep <- (SSprep/dfTreat)/(SSRes/dfRes) + # Regression + # ANOVA tape II SS of regression + SSreg <- Anova(lm(readout ~log_dose + isSample, circ_ABl))[1,1] + # Non-parallelism + # diff of RSS of restricted and unrestricted model + SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2) + F_nonpar <- SSnonpar/(sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)/(lenCirc-4)) + + # non-linearity + SSnonlin <- sum((predict(modABu)-predict(lm(readout ~ as.factor(log_dose)*isSample, circ_ABl)))^2) + # = RSS-SSE + # Total SS + SStot <- sum((circ_ABl$readout-mean(circ_ABl$readout))^2) + # Significance of R^2 F-ratio + # MSR/MSE + # sample A + F_R2_A <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2)/ + (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - circ_Al$readout)^2)/(nrow(circ_Al)-3)) + pFR2_A <- round(pf(F_R2_A,1,6),4) + # sample B + F_R2_B <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2)/ + (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2)/(nrow(circ_Bl)-3)) + pFR2_B <- round(pf(F_R2_B,1,6),4) + # sign of non-lin with pure error: MSSnonlin/MSSE + F_nonlin <- (SSnonlin/2)/(SSE/dfPureE) + + # sign of slope + F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2)/(sum((circ_Bl$readout - predict(modB))^2)/(nrow(circ_Bl)-2)) + F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2)/(sum((circ_Al$readout - predict(modA))^2)/(nrow(circ_Al)-2)) + # F-test on regression: MSSreg/MSSE + if (is.na(F_nonlin)) F_nonlin <- 0 + if (F_nonlin > 0) { + p_F_nonlin <- round(pf(F_nonlin,2,dfPureE, lower.tail = F),5) + } else { p_F_nonlin <- "SSnonlin neg or 0"; } + + # significances + F_regr <- (SSreg/1)/(SSRes/dfRes) + p_F_regr <- round(pf(F_regr,1,dfRes, lower.tail = F),3) + p_F_treat <- round(pf(F_treat,3,dfRes, lower.tail = F),3) + p_F_prep <- round(pf(F_prep,1,dfRes, lower.tail = F),3) + p_F_slope_A <- round(pf(F_slope_A,1,(nrow(circ_Al)-2), lower.tail = F),3) + p_F_slope_B <- round(pf(F_slope_B,1,(nrow(circ_Bl)-2), lower.tail = F),3) + p_F_nonp <- round(pf(F_nonpar,1,dfRes, lower.tail = F),3) + p_F_LoF <- p_F_nonlin + + res_tab_lin <- data.frame(test = c("F-test on sign. of regression", "F_test on non-lin", + "F-test on R^2 A","F_test on R^2 B", + "F-test on slope A","F-test on slope B", + "F-test on non-parallelism","F-test on preparation"), + test_results = c(ifelse(p_F_regr<0.05,0,1),ifelse(p_F_nonlin<0.05,1,0), + ifelse(pFR2_A<0.05,1,0),ifelse(pFR2_B<0.05,1,0), + ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1), + ifelse(p_F_nonp<0.05,1,0),ifelse(p_F_prep<0.05,0,1)), + estimate = c(p_F_regr, p_F_nonlin,pFR2_A,pFR2_B,p_F_slope_A, + p_F_slope_B,p_F_nonp,p_F_prep), + Source = c("Treatment","Preparation","Regression","Non-parallelism", + "Resid Error","Non-linearity","Pure error", "Total"), + df = c(dfTreat,1,1,1,dfRMSE,2,dfPureE,lenCirc-1), + SumSquares = c(round(SStreat,5),round(SSprep,5),round(SSreg,5), + round(SSnonpar,5),round(RSS,5),round(SSnonlin,5), + round(SSE,5),round(SStot,5)), + MS = c(round(SStreat/dfTreat,5),round(SSprep,5),round(SSreg,5), + round(SSnonpar,5),round(RSS/dfRMSE,5),round(SSnonlin/2,5), + round(SSE/dfPureE,5),round(SStot/dfTotal,5)), + "F-value" = c(round(F_treat,5), round(F_prep,5),round(F_regr,5), + round(F_nonpar,5),"",round(F_nonlin,5),"",""), + "p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "","")) + RET <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2) + return(RET) +} + + +#' Plots linear regression on the scatterplot. +#' +#' Restricted and unrestricted model are plotted. The number of dilution steps used for the regression +#' are printed in the title. Reference is blue, red is the sample. The points used for regression +#' are highlighted with a circle. +#' +#' @param circle Data frame with the 3 consecutive concentrations and the respective readouts. +#' @param sigmoid Parameters of the unrestricted fit of the full dataset for drawing the curves. +#' @param all_l2 Long format data.frame of the readout data incl log(concentration) and isRef (0s and 1s). +#' @param pl_df Data.frame for plotting the regression lines. +#' @param indS Index of dilution, where regression starts for reference sample. +#' @param indT Index of dilution, where regression starts for test sample. +#' @returns A grid object of 2 plots of unrestricted and restricted linear PLA models. +#' @export +#' @examples +#' circle <- read.csv("~/plateflow/circle.csv") +#' all_l2 <- read.csv("~/plateflow/all_l2.csv") +#' sigmoid <- c(10.0, 10.0, 110.0, 110.0, 1.0, 1.0, -3.5, 0.0) +#' indS <- 3 +#' indT <- 3 +#' pl_df <- data.frame(lnC=c(-1.203973,-1.897120 ,-2.590267,-3.283414,-3.976562,-4.669176,-5.362323,-6.053340), +#' plotS = c(113.772511,97.668371,81.564231,65.460091,49.355952,33.264200,17.160060,1.105405), +#' plotT = c(114.213375,97.588663,80.963951,64.339239,47.714527,31.102604,14.477892,-2.095735)) +#' +#' +#' PlotLinPLA_FUNC(circle, sigmoid, all_l2, pl_df, indS,indT) + + +PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) { + #browser() + mLin <- gsl_nls(readout ~ (intS+r)*isSample + intS*isRef + k*log_dose, + data=circle, + start=list(intS = 0, k=1,r=0), + control = gsl_nls_control(xtol=1e-10,ftol=1e-10,gtol=1e-10)) + # alternativ: modAB <- lm(readout ~ log_dose+isSample, circle) + sum_mLin <- summary(mLin) + + log_dose <- unique(all_l2$log_dose) + seq_x <- seq(min(log_dose), max(log_dose),0.1) + SAMPLEtrue <- sigmoid[5] + (sigmoid[7]-sigmoid[5])/(1+exp(sigmoid[6]*((sigmoid[4]-sigmoid[8]-seq_x)))) + REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[2]*((sigmoid[4]-seq_x)))) + + truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue) + + p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) + + geom_point(size=2) + + #labs(title=paste("linear regression model", indS,indT), color="product") + + scale_colour_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) + + ylim(min(all_l2$readout),max(all_l2$readout)) + + theme_bw() + p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="#4545BA", + inherit.aes = F) + + geom_line(data=pl_df,aes(x=lnC,y=plotT),color="#C2173F", + inherit.aes = F) + + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, + inherit.aes = F) + + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4, + inherit.aes = F) + + labs(title = paste("unrestricted PLA model"), subtitle = paste("Regression starts for reference sample:",indS, "for test sample:",indT)) + + theme(legend.position="none", axis.text = element_text(size=14)) + p3 <- p2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef), + size=5,alpha=0.2), inherit.aes = FALSE) + + scale_shape_manual(labels=c("test","reference"), values=c(21,21)) + # fit intercept for test and ref and common slope + + + pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log_dose + pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log_dose + pl_rest <- data.frame(lnC=log_dose, plotS=pl_restS, plotT=pl_restT) + + pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="#4545BA", + inherit.aes = F) + + geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="#C2173F", + inherit.aes = F) + + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, + inherit.aes = F) + + geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4, + inherit.aes = F) + + labs(title = paste("restricted linear regression model"), + subtitle = paste("Regression on highlighted points")) + + 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), col=c("#1FE8B2"), stroke=2, inherit.aes = FALSE) + + scale_shape_manual(labels=c("test","reference"), values=c(21,21)) + return(grid.arrange(p3,pr3,nrow=1)) +} + + + +#' Calculates the potency of 4PL PLA of all models model +#' +#' The gradient method is used for calculating the potency for a restricted model, an unrestricteed model, +#' and the log-transformations thereof. +#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. +#' +#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(concentration). +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results +#' 2) summary of unrestricted linear model 3) Slope difference +CIs +#' 4) summary of restricted linear model. +#' @export +#' @examples +#' ro_new <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), +#' R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), +#' T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), +#' T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), +#' log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) +#' PureErrF <- TRUE +#' +#' +#' pot4plFUNC(ro_new, PureErrF) + + +pot4plFUNC <- function(ro_new, PureErrFlag) { + all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") + isRef <- rep(c(1,0),1,each=nrow(all_l)/2) + isSample <- rep(c(0,1),1,each=nrow(all_l)/2) + all_l$isRef <- isRef + all_l$isSample <- isSample + all_l$Conc <- exp(all_l$log_dose) + all_l$readout[all_l$readout < 0] <- 0.01 + all_l$readouttrans <- log(all_l$readout) + #browser() + CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) + if (CORdat<0) SLOPE <- -1 else SLOPE <- 1 + # + FITs <- Fitting_FUNC(ro_new, TransFlag = F) + if (!PureErrFlag) { + pot_est <- FITs[[3]] + potU_est <- FITs[[4]] + } else { + FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) + meanPureErr <- FitAnova[4,3] + DFsPure <- FitAnova[4,1] + # SU_mr <- tryCatch({ + # SU_mr <- summary(mr) + # }, error = function(msg) { + # return() + # }) + SU_mr <- FITs[[1]] + SU_mrCoeff <- SU_mr$coefficients + # if (length(SU_mr)>1) { + # SU_mrCoeff <- SU_mr$coefficients + # } else { SU_mrCoeff <- rep(NA,5) } + #te2$cov.unscaled[1,1]*te2$sigma^2 + #VCOV <- vcov(mr) + #V_V <- VCOV/SU_mr$sigma^2 + V_V <- SU_mr$cov.unscaled + SEsPure <- sqrt(diag(V_V)*meanPureErr) + pot_est <- c(exp(SU_mrCoeff['r',1]),exp(SU_mrCoeff['r',1]-qt(0.975,DFsPure)*SEsPure['r']), + exp(SU_mrCoeff['r',1]+qt(0.975,DFsPure)*SEsPure['r'])) + # unrestricted + SU_mu <- FITs[[2]] + s_muCoeff <- SU_mu$coefficients + #SU_mu <- summary(mu) + #VCOVu <- vcov(mu) + V_Vu <- SU_mu$cov.unscaled + SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) + potU_est <- c(exp(s_muCoeff['r',1]),exp(s_muCoeff['r',1]-qt(0.975,DFsPure)*SEsPureU['r']), + + exp(s_muCoeff['r',1]+qt(0.975,DFsPure)*SEsPureU['r'])) + } # PureErrFlag + + FITstrans <- Fitting_FUNC(ro_new, TransFlag = TRUE) + + + # pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic")) + # potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic")) + pot_est_log <- FITstrans[[3]] + potU_est_log <- FITstrans[[4]] + colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI") + colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI") + #browser() + su_mr_log <- FITstrans[[1]] #summary(mr_log) + Dat$RMSE_Rlog <- su_mr_log$sigma + su_mu_log <- FITstrans[[2]] # summary(mu_log) + Dat$RMSE_Ulog <- su_mu_log$sigma + Dat$up_lowAslog <- su_mu_log$coefficients[1,1] - su_mu_log$coefficients[4,1] + potALL <- rbind(round(pot_est,6), round(potU_est,6), round(pot_est_log,6), round(potU_est_log,6)) + + potALL2 <- cbind(c("restricted","unrestricted","ln-transformed restr","ln-transformed unrestr"), potALL) + return(potALL2) +} + +#' Calculates logarithmized CIs for ratios +#' +#' Ratio CIs can be unbound, and not calulatable. Therefore the ratios are logarithmized. +#' The Covariance is not used, as it is 0 for comparisons between products. +#' +#' @param xs The estimate of the reference sample. +#' @param xt The estimate of the test sample. +#' @param se_xs The standard error of the estimate of the reference sample. +#' @param se_xt The standard error of the estimate of the test sample. +#' @param CoVar The covariance between test and reference estimate (set to 0). +#' @param DFs Degrees of freedom, either from pure error term (no of readouts - no of concentrations*2), or RMSE term (no of readouts - 8 parameters). +#' @param Conf The confidence of the Confidence Interval (default 95% which requires 0.975). +#' @returns A data-frame with the lower and upper CI in anti-log form. +#' @export +#' @examples +#' xs=2; xt=3.2; se_xt=0.34;se_xs=0.23; DFs=32-16 +#' +#' +#' ParamCI_F(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) + + + +ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) { + log_xs <- log(abs(xs)) + log_xt <- log(abs(xt)) + var_log_xs <- (se_xs/xs)^2 # approximate variance of log(xs) + var_log_xt <- (se_xt/xt)^2 + se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVar/(xs*xt) + + lower_log_ratio <- log_xt-log_xs - qt(Conf,DFs)*se_log_ratio + upper_log_ratio <- log_xt-log_xs + qt(Conf,DFs)*se_log_ratio + ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) + return(ci_ratio) +} + + +#' Calculates EQ tests and F-Tests for 4P +#' +#' CIs of asmptote ratios,slope ratios, asympt-difference ratio, and lower asympt-diff. calculated. +#' In addition the F-test on significance of regression and non-linearity are calculated. +#' +#' @param ro_new Data frame with the concentrations and the respective readouts. +#' @param Lim The limits to check if the relative CI is within the limits. +#' @param PureErrFlag Boolean if Pure Error should be used. +#' @returns A data-frame with the lower and upper CI in anti-log form. +#' @export +#' @examples +#' +#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' Lim <- c(-1, 1, 0.005 , 2, 0.5, 2, 0.5,2,75,133, 75, 133) +#' PureErrF <- FALSE +#' +#' tests_FUNC(ro_new, Lim,PureErrF) + + +tests_FUNC <- function(ro_new, Lim, PureErrFlag) { + + all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") + isRef <- rep(c(1,0),1,each=nrow(all_l)/2) + isSample <- rep(c(0,1),1,each=nrow(all_l)/2) + all_l$isRef <- isRef + all_l$isSample <- isSample + all_l$Conc <- exp(all_l$log_dose) + all_l$readout[all_l$readout < 0] <- 0.01 + #browser() + FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) + POTr_CI <- FITs[[3]][2:3] + potAll2 <- FITs[[3]] + + smr <- FITs[[1]] + smu <- FITs[[2]] + FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) + # pure error + pureSS <- FitAnova[4,2] + pureSS_df <- FitAnova[4,1] + meanPureErr <- FitAnova[4,3] + #vcovMU <- vcov(mu) + V_V <- smu$cov.unscaled + SEsPure <- sqrt(diag(V_V)*meanPureErr) + VCOVpure <- V_V*meanPureErr + DFsPure <- FitAnova[4,1] + + + + testPOTr <- logical() + if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100Lim[[9]]/100 & potAllU2[3] 0) { + p_F_nonlin <- round(pf(F_nonlin, AnovaDFs[6], ERR_df, lower.tail = F),5) + } else { p_F_nonlin <- "SSnonlin neg or single dilutions" } + + test_a <- test_b <- test_d <- test_ad <- logical() + + RSS_r <- round(sum(smr$residuals^2),5) + MSE_r <- RSS_r/(nrow(all_l)-5) + RMSE_r <- round(sqrt(MSE_r),6) + Dat$RMSE_r <- RMSE_r + Dat$RMSE_pure <- RMSE_pure + Dat$RMSE_unr <- round(RMSEunr,6) + + coeffs <- smu$coefficients[,1] + #browser() + ## EQ test on lower As diff + as <- coeffs["as"] + at <- coeffs["at"] + lAs_diff <- (at-as) + uCI_laDiff <- lAs_diff+qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2) + lCI_laDiff <- lAs_diff-qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2) + if (uCI_laDiff < Lim[[2]] & lCI_laDiff > Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1 + + #### EQ test on upper asymptote ratio ---- + # as <- coeffs["as"] + # at <- coeffs["at"] + # uAsRatio <- compParm(potU, "a","/",display=F) + # uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2]) + #browser() + ds <- smu$coefficients["ds",1] + dt <- smu$coefficients["dt",1] + if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds","ds"]) else se_ds <- smu$coefficients["ds",2] + if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt","dt"]) else se_dt <- smu$coefficients["dt",2] + if (PureErrFlag) CoVarlog_d <- VCOVpure["dt","ds"] else CoVarlog_d <- vcovMU["dt","ds"] + if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-8 + uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.975) + if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1 + estUppA <- round(at/as,5) + + Dat$uAsCI <- uAsCI2 + + #### EQ test on slope ratio ---- + # bs <- coeffs["bs"] + # bt <- coeffs["bt"] + # slopeRatio <- compParm(potU, "b","/",display=F) + # slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2]) + + bs <- smu$coefficients["bs",1] + bt <- smu$coefficients["bt",1] + if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs","bs"]) else se_bs <- smu$coefficients["bs",2] + if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt","bt"]) else se_bt <- smu$coefficients["bt",2] + if (PureErrFlag) CoVarlog_b <- VCOVpure["bt","bs"] else CoVarlog_b <- vcovMU["bt","bs"] + slopeCI2 <- ParamCI_F(bt,bs,se_bt, se_bs,CoVarlog_b, DFs, Conf=0.975) + if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1 + estUppA <- round(at/as,5) + + Dat$slopeRatioCI <- slopeCI2 + + #### EQ test on lower As ratio ---- + + # lAsRatio <- compParm(potU, "d","/",display=F) + # slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2]) + + as <- smu$coefficients["as",1] + at <- smu$coefficients["at",1] + if (PureErrFlag) se_as <- sqrt(VCOVpure["as","as"]) else se_as <- smu$coefficients["as",2] + if (PureErrFlag) se_at <- sqrt(VCOVpure["at","at"]) else se_at <- smu$coefficients["at",2] + if (PureErrFlag) CoVarlog_a <- VCOVpure["at","as"] else CoVarlog_a <- vcovMU["at","as"] + lAsCI2 <- ParamCI_F(at,as,se_at, se_as,CoVarlog_a, DFs, Conf=0.975) + if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1 + estLowA <- round(at/as,5) + + Dat$lAsCI <- lAsCI2 + + #### EQtest on ratio of As difference ---- + AsDiffRatio <- (dt-at)/(ds-as) + + dt_at <- (dt-at) + ds_as <- (ds-as) + se_ds_asPure <- sqrt(VCOVpure["as","as"]+VCOVpure["ds","ds"]-2*VCOVpure["as","ds"]) + se_dt_atPure <- sqrt(VCOVpure["at","at"]+VCOVpure["dt","dt"]-2*VCOVpure["at","dt"]) + se_ds_asRMSE <- sqrt(vcovMU["as","as"]+vcovMU["ds","ds"]-2*vcovMU["as","ds"]) + se_dt_atRMSE <- sqrt(vcovMU["at","at"]+vcovMU["dt","dt"]-2*vcovMU["at","dt"]) + if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE + if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE + + AsDiffCI2 <- ParamCI_F(dt_at,ds_as,se_dt_at, se_ds_as,CoVar=0, DFs, Conf=0.975) + if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1 + estLowA <- round(at/as,5) + + Dat$up_lowAs <- abs(ds-as) + + lowerCIlowerA <- lAsCI2[1]; lowerCIupperA <- uAsCI2[1]; upperCIlowerA <- lAsCI2[2]; upperCIupperA <- uAsCI2[2] + test_lowA <- test_d; test_uppA <- test_a + #browser() + res_tab <- data.frame(test= c("F-test on sign. of regression*", + "EQ test on lower asymptotes difference", + "EQ test ratio of lower asymptotes", + "EQ test ratio of Hill slopes", + "EQ test ratio of upper asymptotes", + "F-test on non-linearity*", + "EQ test ratio of asymptote difference", + "geom. rel. CI restr. model", + "geom. rel. CI unrestr. model"), + test_results = c(ifelse(p_F_regr<0.05,0,1), test_la_diff, test_lowA, test_b, test_uppA, + ifelse(p_F_nonlin>1,1, ifelse(p_F_nonlin<0.05,1,0)), test_ad, + testPOTr, test_c), + estimate = c(round(p_F_regr, 3), round(lAs_diff, 5), + estLowA, round(bs/bt,5), estUppA, p_F_nonlin, + round(dt_at/ds_as, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)), + lower_limit = c("-",Lim[[1]],Lim[[3]],Lim[[5]],Lim[[7]],"-",Lim[[11]],Lim[[9]],Lim[[9]]), + upper_limit = c("-",Lim[[2]],Lim[[4]],Lim[[6]],Lim[[8]],"-",Lim[[12]],Lim[[10]],Lim[[10]]), + lower_CI = c(RMSE_r, round(lCI_laDiff,3), round(lAsCI2[1],5), round(slopeCI2[1],5), + round(uAsCI2[1],5), "-", round(AsDiffCI2[1],5), round(potAll2[2],2), round(potAllU2[2],2)), + upper_CI = c(RMSE_pure, round(uCI_laDiff,3), round(lAsCI2[2],5), round(slopeCI2[2],5), + round(uAsCI2[2],5), "-", round(AsDiffCI2[2],5), round(potAll2[3],2), round(potAllU2[3],2)) + ) + return(res_tab) +} + +#' ANOVA unrestricted model +#' +#' ANOVA with unrestricted model. +#' +#' +#' @param ro_new Data frame with the concentrations and the respective readouts. +#' @returns A data-frame with the analysis of variance. +#' @export +#' @examples +#' +#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' +#' +#' ANOVA4plUnresfunc(ro_new) + + +ANOVA4plUnresfunc <- function(ro_new) { + all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") + all_len <- nrow(all_l) + isRef <- rep(c(1,0),1,each=all_len/2) + isSample <- rep(c(0,1),1,each=all_len/2) + all_l$isRef <- isRef + all_l$isSample <- isSample + all_l$Conc <- exp(all_l$log_dose) + all_l$readout[all_l$readout < 0] <- 0.01 + + FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) + smr <- FITs[[1]] + smu <- FITs[[2]] + smrPREDs <- FITs[[5]] + smuPREDs <- FITs[[6]] + # pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), + # pmodels=data.frame(1,1,1,isSample)) + # potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), + # pmodels=data.frame(isSample, isSample,isSample,isSample)) + SStreat <- round(sum((smuPREDs - mean(all_l$readout))^2),5) + SStreat_df <- length(unique(all_l$log_dose))-1 + SSregr <- round(sum((smrPREDs-mean(all_l$readout))^2),5) + ## Non-parallel + SSnonparallel <- round(sum(smr$residuals^2) - sum(smu$residuals^2),5) + ## Preparation + SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) + ## Resid Err + RSS <- round(sum(smu$residuals^2),5) + RSS_df <- nrow(all_l)-SStreat_df-1 + FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) + # PureErr + SSE <- round(FitAnova[4,3],5) + SSE_df <- FitAnova[4,1] + # Non-Linearity + SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - smuPREDs)^2),4) + LoF_df <- FitAnova[1,1]+FitAnova[2,1] + ## Total + SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5) + MSE <- RSS/RSS_df + noConc <- length(unique(all_l$Conc)) + AnovaDFs <- c(noConc-1, 1,3,noConc-4-1, nrow(all_l)-noConc, noConc, nrow(all_l)-noConc-noConc, nrow(all_l)-1) + p_SStreat <- round(pf((SStreat/AnovaDFs[1])/MSE, AnovaDFs[1],RSS_df, lower.tail = F),3) + p_SSprep <- round(pf((SSprep/AnovaDFs[2])/MSE, AnovaDFs[2],RSS_df, lower.tail = F),3) + p_SSregr <- round(pf((SSregr/AnovaDFs[3])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) + p_SSnonp <- round(pf((SSnonparallel/AnovaDFs[4])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) + p_SSLoF <- round(pf((SSnonlin/LoF_df)/(SSE/SSE_df), LoF_df,SSE_df, lower.tail = F),5) + + ANOVAtab <- data.frame(Source = c("Treatment","Preparation","Regression", + "Non-Parallelism","Residual Error","Non-linearity", + "Pure Error","Total"), + DF = AnovaDFs, + SumSquares = c(SStreat, SSprep,SSregr, SSnonparallel, + RSS, SSnonlin,SSE, SStot), + MeanSquares = c(round(SStreat/AnovaDFs[1],3), SSprep, round(SStreat/AnovaDFs[3],3),round(SSnonparallel/AnovaDFs[4],3), + round(MSE,5), round(SSnonlin/LoF_df,5), round(SSE/SSE_df,5),""), + "F-value" = c(round((SStreat/AnovaDFs[1])/MSE,5), round((SSprep/AnovaDFs[2])/MSE,5), + round((SSregr/AnovaDFs[3])/MSE,5),round((SSnonparallel/AnovaDFs[4])/MSE,5), + "",round((SSnonlin/LoF_df)/(SSE/SSE_df),5),"",""), + "p_value" = c(round(p_SStreat,3), p_SSprep, round(p_SSregr,3), p_SSnonp,"",p_SSLoF,"","") + ) + + return(ANOVAtab) +} + +#' Calculates descriptive statistics per concentration +#' +#' Mean, SD, and CV% per concentration and product.. +#' +#' +#' @param ro_new Data frame with the concentrations and the respective readouts. +#' @returns A data-frame with the concentrations and metadata. +#' @export +#' @examples +#' +#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), +#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), +#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), +#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) +#' +#' +#' perConcTab(ro_new, noDilSeries=3) + + + +perConcTab <- function(ro_new, noDilSeries) { + Reftab <- ro_new[,c(1:noDilSeries)] + Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))] + tReftab <- t(Reftab) + colnames(tReftab) <- round(ro_new[,ncol(ro_new)],5) + + avs <- apply(tReftab,2,mean) + sds <- apply(tReftab,2,sd) + cv <- sds/avs*100 + tReftab2 <- rbind(tReftab, avs,sds,cv) + + tTesttab <- t(Testtab) + colnames(tTesttab) <- round(ro_new[,ncol(ro_new)],5) + + avs_test <- apply(tTesttab,2,mean) + sds_test <- apply(tTesttab,2,sd) + cv_test <- sds_test/avs_test*100 + tTesttab2 <- rbind(tTesttab, avs_test,sds_test,cv_test) + concTab <- rbind(tReftab2, tTesttab2) + return(concTab) + +} + +#' Calculates dilution series. +#' +#' Recursive function to calcuate a geometric dilution series. +#' +#' @param x Starting concentration. +#' @param Div Divisor to get next concentration. +#' @param N Number of dilution step, increases. +#' @param res Vector of results. +#' @param noDil Final number of concentrations -1 (the initial one). +#' @returns A data-frame with the concentrations and metadata. +#' @export +#' @examples +#' +#' x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7 +#' +#' divFUN(x,Div,N,res,noDil) + + + +divFUN <- function(x,Div,N,res,noDil) { + N <- N+1 + y <- x/Div + res <- c(res,y) + if (N==noDil) { return(res) } + divFUN(y,Div,N,res,noDil) +} + diff --git a/RoxygenTestFile.Rmd b/RoxygenTestFile.Rmd new file mode 100644 index 0000000..72a1580 --- /dev/null +++ b/RoxygenTestFile.Rmd @@ -0,0 +1,51 @@ +--- +title: "RoxygenFileTest" +author: "F Innerbichler" +date: "2026-05-13" +output: pdf_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +library(shiny) +source("Global.R") + +library(testthat) + +``` + +```{r LinPotTab, echo=F} + + +library(car) +CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9, +-3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7), +replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3", + "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"), + readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2), + isRef=c(rep(1,9), rep(0,9)), + isSample = c(rep(0,9), rep(1,9))) +Lim <- c(rep(0,8), 70,130) +PureErrF <- TRUE + + +TestTab <- as.matrix(LinPotTab(circles=CIRC, Lim, PureErrF)) +# SolTab is the Solution table +SolTab <- matrix(c(104.959, 87.994, 125.196, 0, 83.836, 119.281),nrow=1) +colnames(SolTab) <- c("Potency", "lower 95%CI", "upper 95%CI", "test_result", "lowerRel95%CI", "upperRel95%CI") +#all.equal(TestTab, SolTab) + +expect_equal(TestTab, SolTab, check.attributes = FALSE) +# no output means "all_equal" + +``` + + + + + +```{r pressure, echo=FALSE} + +``` + +Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot. diff --git a/TestFileRoxygen.R b/TestFileRoxygen.R index 48d6737..fd73aec 100644 --- a/TestFileRoxygen.R +++ b/TestFileRoxygen.R @@ -202,3 +202,19 @@ x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7 divFUN(x,Div,N,res,noDil) # [1] 0.3333333333 0.1111111111 0.0370370370 0.0123456790 0.0041152263 0.0013717421 0.0004572474 + + +circle <- read.csv("~/plateflow/circle.csv") +all_l2 <- read.csv("~/plateflow/all_l2.csv") +sigmoid <- c(10.0, 10.0, 110.0, 110.0, 1.0, 1.0, -3.5, 0.0) +indS <- 3 +indT <- 3 +pl_df <- data.frame(lnC=c(-1.203973,-1.897120 ,-2.590267,-3.283414,-3.976562,-4.669176,-5.362323,-6.053340), +plotS = c(113.772511,97.668371,81.564231,65.460091,49.355952,33.264200,17.160060,1.105405), +plotT = c(114.213375,97.588663,80.963951,64.339239,47.714527,31.102604,14.477892,-2.095735)) + + +PLOT <- PlotLinPLA_FUNC(circle, sigmoid, all_l2, pl_df, indS,indT) +PlotLinTest.png + + diff --git a/app.R b/app.R index 723cd91..6c8e123 100644 --- a/app.R +++ b/app.R @@ -24,1202 +24,7 @@ library(dplyr) Dat <- reactiveValues() REP <- reactiveValues() -#### functions ---- - - - -# dilFUN2 <- function(cs_,dils,Faktor) { -# av <- cs_ -# dils_av <- dils -# dils_avsc <- dils_av*Faktor -# dils2 <- dils_avsc+av -# dilfactors <- 1/exp(dils2-lag(dils2)) -# return(dilfactors) -# } - -#' Levenberg Marquard fit of 4 pl -#' -#' Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, -#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. -#' -#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. -#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). -#' @returns Returns list of following: summary of restricted and unrestricted model fits of the 4pl function, -#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof. -#' @export -#' @examples -#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), -#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), -#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), -#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) -#' TransF<- FALSE -#' Dat<- list() -#' te <- Fitting_FUNC(dat,TransF) -#' print(te) - - -Fitting_FUNC <- function(ro_new, TransFlag=F) { - CORro <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) - #browser() - all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") - isRef <- rep(c(1,0),1,each=nrow(all_l)/2) - isSample <- rep(c(0,1),1,each=nrow(all_l)/2) - all_l2 <- cbind(all_l, isRef, isSample) -#browser() - if (CORro<0) SLOPE <- -1 else SLOPE <- 1 - if (!TransFlag) { - startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0) - - mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), - data=all_l2, - start=startlist,#trace=T, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - s_mr <- tryCatch({ - s_mr <- summary(mr) - }, - error = function(err) { - s_mr <- NULL - }) - } else { - startlist <- list(a=log(min(ro_new[,2])), b=SLOPE, d=log(max(ro_new[,2])), cs=mean(all_l$log_dose),r=0) - mrT <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), - data=all_l2, - start=startlist, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - s_mr <- summary(mrT) - } - - if (!TransFlag) { - startlistmu <- list(as=min(ro_new[,2]), bs=SLOPE, ds=max(ro_new[,2]), cs=mean(all_l$log_dose), - at=min(ro_new[,2]), bt=SLOPE, dt=max(ro_new[,2]), r=0) - tryCatch({ - mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(msg){ - return(0) }) - - Sum_u <- tryCatch({ summary(mu) }, - error=function(msg){ - return(0) }) - } else { - startlistmu <- list(as=log(min(ro_new[,2])), bs=SLOPE, ds=log(max(ro_new[,2])), cs=mean(all_l$log_dose), - at=log(min(ro_new[,2])), bt=SLOPE, dt=log(max(ro_new[,2])), r=0) - tryCatch({ - muT <- gsl_nls(fn = log(readout) ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))), - data=all_l, - start=startlistmu, - control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) - }, - error = function(msg){ - return(0) }) - - Sum_u <- tryCatch({ summary(muT) }, - error=function(msg){ - return(0) }) - } - if (!TransFlag) { - pot_est <- exp(confintd(mr, "r", method="asymptotic")) - potU_est <- exp(confintd(mu, "r", method="asymptotic")) - PRED <- predict(mr) - PREDu <- predict(mu) - } else { - pot_est <- exp(confintd(mrT, "r", method="asymptotic")) - potU_est <- exp(confintd(muT, "r", method="asymptotic")) - PRED <- predict(mrT) - PREDu <- predict(muT) - } - return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu)) -} - - -#' Plot sigmoidal curve -#' -#' Returns the final plots of the 4pl function as sigmoidal lines, and the single readouts as scatter, with REF in blue and TEST in red. -#' Two plots are returned in a grid object: the unrestricted model and the restricted model. -#' -#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. -#' @param sigmoid The parameters of the 4pl fit of unrestricted model from EXCEL input. -#' @param det_sig The parameters of the 4pl fit of unrestricted model from the metadata input. -#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE). -#' @returns A grid object either of the original scale or the natural log of the readouts. -#' @export -#' @examples -#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), -#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), -#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), -#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) -#' sigmoid <- c(0.7163324, 0.5636804, 10.6156340 , 9.9784160, -0.7504673, -0.7108692, -3.5788141, -0.6662962) -#' det_sig <- FALSE -#' TransF<- FALSE -#' Dat<- list() -#' p <- plot_f(dat,sigmoid,det_sig,TransF) -#' print(p) - - -plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig, - CORdat <- cor(dat[,1],dat[,ncol(dat)]) - #browser() - all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") - isRef <- rep(c(1,0),1,each=nrow(all_l)/2) - isSample <- rep(c(0,1),1,each=nrow(all_l)/2) - all_l2 <- cbind(all_l, isRef, isSample) - #browser() - MODLS <- Fitting_FUNC(dat, TransFlag=F) - s_mr <- MODLS[[1]] - a <- s_mr$coefficients["a",1] - b <- s_mr$coefficients["b",1] - cs <- s_mr$coefficients["cs",1] - d <- s_mr$coefficients["d",1] - r <- s_mr$coefficients["r",1] - - log_dose <- unique(all_l$log_dose) - seq_x <- seq(min(log_dose),max(log_dose),0.1) - SAMPLE <- a+(d-a)/(1+exp(b*((cs-r)-seq_x))) - REF <- a+(d-a)/(1+exp(b*((cs)-seq_x))) - - s_mu <- MODLS[[2]] - - as <- s_mu$coefficients["as",1] - bs <- s_mu$coefficients["bs",1] - cs <- s_mu$coefficients["cs",1] - ds <- s_mu$coefficients["ds",1] - at <- s_mu$coefficients["at",1] - bt <- s_mu$coefficients["bt",1] - ct <- s_mu$coefficients["cs",1]-s_mu$coefficients["r",1] - dt <- s_mu$coefficients["dt",1] - r <- s_mu$coefficients["r",1] - - - SAMPLEtrue <- at + (dt -at)/(1+exp(bt*((cs-r-seq_x)))) - REFtrue <- as + (ds -as)/(1+exp(bs*((cs-seq_x)))) - - #browser() - pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue) - all_l2$readout[all_l2$readout < 0] <- 0.01 - all_l2$readouttrans <- log(all_l2$readout) - slopeEC50 <- b*(d-a)/4 - Intercept <- a+(d-a)/2-b*(d-a)/4*cs -#browser() - Xbendl3 <- cs-(1.5434/b) - Xbendu3 <- cs+(1.5434/b) - XbendlT <- cs-r-(1.5434/b) - XbenduT <- cs-r+(1.5434/b) - XasymplS <- cs-(3/b) - XasympuS <- cs+(3/b) - XasymplT <- cs-r-(3/b) - XasympuT <- cs-r+(3/b) - bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3), - bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,3), - asympREF_lower = round(XasymplS,3), asympREF_upper=round(XasympuS,3), - asympSAMPLE_lower = round(XasymplT,3), asympSAMPLE_upper=round(XasympuT,3)) - Dat$bendpoints <- bendpoints - Dat$cfordils <- cs - - p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + - geom_point(shape=factor(isRef), size=3,alpha=0.8) + - labs(title = paste("restricted 4pl"), - color="product") + - scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + - scale_shape_manual(labels=c("test","reference")) + - theme_bw() + - theme(axis.text = element_text(size=14)) - - p2 <- p + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLE), color="#C2173F", - inherit.aes = F) + - geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="#4545BA", - inherit.aes = F) + - geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="#C2173F", linetype=2, alpha=0.4, - inherit.aes = F) + - geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="#4545BA", linetype=2, alpha=0.4, - inherit.aes = F) + - geom_vline(xintercept=c(Xbendl3, Xbendu3), col="#4545BA",linetype=2) + - geom_vline(xintercept=c(XbendlT, XbenduT), col="#C2173F",linetype=2) + - geom_vline(xintercept=c( XasymplS, XasympuS), col="#4545BA55",linetype=2) + - geom_vline(xintercept=c(XasymplT, XasympuT), col="#C2173F55",linetype=2) + - annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) + - geom_abline(slope = slopeEC50, intercept = Intercept) + - theme(legend.position="none") - Dat$p2 <- p2 - - # transformed plots - p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + - geom_point(shape=factor(isRef), size=3,alpha=0.8) + - labs(title = paste("restricted transformed 4pl"), color="product") + - scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + - theme_bw() - - MODLStrans <- Fitting_FUNC(dat,TransFlag = TRUE) - s_mrt <- MODLStrans[[1]] - a_trans <- s_mrt$coefficients["a",1] - b_trans <- s_mrt$coefficients["b",1] - cs_trans <- s_mrt$coefficients["cs",1] - d_trans <- s_mrt$coefficients["d",1] - r_trans <- s_mrt$coefficients["r",1] - - XbendlTrans <- cs_trans-(1.5434/b_trans) - XbenduTrans <- cs_trans+(1.5434/b_trans) - XbendlTransT <- cs_trans-r_trans-(1.5434/b_trans) - XbenduTransT <- cs_trans-r_trans+(1.5434/b_trans) - bendpointsTRANS <- c(bendREF_lower = round(XbendlTrans,3), bendREF_upper=round(XbenduTrans,3), - bendSAMPLE_lower = round(XbendlTransT,3), bendSAMPLE_upper=round(XbenduTransT,3)) - Dat$bendpointsTRANS <- bendpointsTRANS - SAMPLEtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans-r_trans)-seq_x))) - REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans)-seq_x))) - - pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans) - p_rt2 <- p_rt + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=SAMPLEtrans), color="#C2173F", - inherit.aes = F) + - geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="#4545BA", - inherit.aes = F) + - geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="#4545BA",linetype=2) + - geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="#C2173F",linetype=2) + - theme(legend.position = "none", axis.text=element_text(size=14)) - - #browser() - Sum_u <- MODLS[[2]] - #browser() - #if (is.null(det_sig)) { - ast <- Sum_u$coefficients["as",1] - ate <- Sum_u$coefficients["at",1] - bst <- Sum_u$coefficients["bs",1] - bte <- Sum_u$coefficients["bt",1] - cst <- Sum_u$coefficients["cs",1] - cte <- Sum_u$coefficients["cs",1]-Sum_u$coefficients["r",1] - dst <- Sum_u$coefficients["ds",1] - dte <- Sum_u$coefficients["dt",1] - # } else { - # ast <- det_sig[5] - # ate <- det_sig[6] - # bst <- det_sig[1] - # bte <- det_sig[2] - # cst <- det_sig[7] - # cte <- det_sig[8] - # dst <- det_sig[3] - # dte <- det_sig[4] - # } - REFu <- ast + (dst-ast)/(1+exp(bst*(cst-seq_x))) - SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(cte-seq_x))) - pl_df2 <- cbind(seq_x, SAMPLEu, REFu) - #browser() - pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) + - geom_point(size=3) + - labs(title="unrestricted 4_pl-Model", color="product") + - scale_color_manual(labels = c("test","reference"), values=c("#C2173F88","#4545BA88")) + - theme_bw() - pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu), - color="#C2173F", inherit.aes = F) + - geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu), - color="#4545BA", inherit.aes = F, - show.legend = F) - pu2_ <- pu2 + - theme(legend.position = "none", axis.text = element_text(size=14)) - - putrans <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) + - geom_point( size=3) + - labs(title="unrestricted transformed 4_pl-Model", color="product") + - scale_color_manual(labels = c("test","reference"), values=c("#C2173FAA","#4545BAAA")) + - theme_bw() - - Sum_ut <- MODLStrans[[2]] - - ast_t <- Sum_ut$coefficients["as",1] - ate_t <- Sum_ut$coefficients["at",1] - bst_t <- Sum_ut$coefficients["bs",1] - bte_t <- Sum_ut$coefficients["bt",1] - cst_t <- Sum_ut$coefficients["cs",1] - cte_t <- Sum_ut$coefficients["cs",1]-Sum_ut$coefficients["r",1] - dst_t <- Sum_ut$coefficients["ds",1] - dte_t <- Sum_ut$coefficients["dt",1] - - REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(cst_t-seq_x))) - SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(cte_t-seq_x))) - pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans) - - pu2_t <- putrans + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=SAMPLEu_trans), - color="#C2173F", inherit.aes = F) + - geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans), - color="#4545BA", inherit.aes = F, - show.legend = F) - pu3_t <- pu2_t - if (TransFlag) grid.arrange(p_rt2,pu3_t, nrow=1) else grid.arrange(p2,pu2_, nrow=1) -} - -#' Simulate a well-plate readout -#' -#' Returns a possible well-plate result based on the metadata provided and a variability. -#' Homo- and heteroscedasticity can be simulated. -#' -#' @param as, bs, cs, ds, at, bt, dt,r,ct, gt, gs are the parameter of the 5pl logistic regression. -#' @param sd_fac The standard deviation with units in % of upper-lower asymptote. -#' @param log_conc The natural log of the concentrations. -#' @param heteroNoise Boolean if heteroscedstic noise should be added. -#' @param noDilSeries A number >1 indicating how many dilution series per product should be simulated. -#' @param noDils A number >7 indicating how many dilutions steps per product should be simulated. -#' @returns A data-frame with readouts and natural log of concentrations. -#' @export -#' @examples -#' as=3; bs=1; cs=-4; ds=10; at=3; bt=1; dt=10;r=0.0001;ct=cs-r;sd_fac=0.1; gt=1; gs=1; -#' lnConc=c(1,0,-1,-2,-3,-4,-5,-6) -#' heteroNoise <- FALSE -#' noDilS <- 3 -#' noD <- 8 -#' Calc_DilRes(as,bs,cs,ds,at,bt,dt,r,ct,sd_fac,gt,gs,log_conc=lnConc, heteroNoise, noDilS, noD) - -Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct=cs-r, - sd_fac=0.1, gt=1, gs=1, log_conc, - heteroNoise=FALSE, noDilSeries, noDils) { - yAxfac <- (ds-as) - log_dose <- log_conc - isRef <- rep(c(1,0),1,each=length(log_conc)*noDilSeries) - isSample <- rep(c(0,1),1,each=length(log_conc)*noDilSeries) - #browser() - av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/ - (1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose))) - if (heteroNoise) { - # heterosc noise - ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac/100))), nrow=noDils, ncol=noDilSeries*2) - } else { - # homosc noise - ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac*yAxfac/100))), nrow=noDils, ncol=noDilSeries*2) - } - ro_jit <- abs(ro_jit) - - ro_jit2 <- cbind(ro_jit, log_dose) - if (noDilSeries==3) { - colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose") - } else { - colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose") - } - return(ro_jit2) -} - -#' Calculate the potency of linear PLA -#' -#' The delta Method is used for calculating the potency using a restricted linear model. -#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. -#' -#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts. -#' @param Lim The limits to check if the relative CI is within the limits. -#' @param PureErrFlag Boolean if Pure Error should be used. -#' @returns A data-frame with potency estimate, absolute CIs, test result, relative CIs. -#' @export -#' @examples -#' CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9, -#' -3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7), -#' replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3", -#' "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"), -#' readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2), -#' isRef=c(rep(1,9), rep(0,9)), -#' isSample = c(rep(0,9), rep(1,9))) -#' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant -#' PureErrF <- TRUE -#' -#' -#' LinPotTab(circles=CIRC, Lim, PureErrF) - - -LinPotTab <- function(circles, Lim, PureErrFlag) { - circ_ABl <- circles - circ_Al <- circ_ABl[circ_ABl$isSample ==1,] - circ_Al <- circ_ABl[circ_ABl$isSample ==0,] - # restr CSSI model - modAB <- lm(readout ~ log_dose + isSample, circ_ABl) - coeffs <- modAB$coefficients - SU_modAB <- tryCatch({ - SU_modAB <- summary(modAB) - }, error = function(msg) { - return(NA) - }) - # Intercept diff/slope modAB - linPot <- exp(modAB$coefficients[3]/modAB$coefficients[2]) - - if(PureErrFlag) { - FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl)) - meanPureErr <- FitAnova[4,3] - DFsPure <- FitAnova[4,1] - VCOV <- vcov(modAB) - V_V <- VCOV/SU_modAB$sigma^2 - VCOVpure <- V_V*meanPureErr - SEsPure <- sqrt(diag(V_V)*meanPureErr) - } - - log_pot_delta <- deltaMethod(modAB, "isSample/log_dose") - if (PureErrFlag) { - V_ <- log_pot_delta$SE^2/SU_modAB$sigma^2 - V_p <- V_*meanPureErr - potDeltaPureSE <- sqrt(V_p) - CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure)*potDeltaPureSE - CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure)*potDeltaPureSE - } else { - CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE - CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE - } - #browser() - ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up)) - if (ExpLinPot[2]*100>Lim[[9]] & ExpLinPot[3]*100Lim[9] & ExpLinPot[3]*100>Lim[10]) test_potCI <- 0 else test_potCI <- 1 - - su_mod <- summary(modAB)$coefficients - su_mod2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.")), su_mod) - su_modU <- summary(modABu)$coefficients - su_modU2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.","slope difference")), su_modU) - - uCI_SloDiff <- su_modU[4,1] + qt(0.975,8)*su_modU[4,2] - lCI_SloDiff <- su_modU[4,1] - qt(0.975,8)*su_modU[4,2] - SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiff,uCI_SloDiff) - - lenCirc <- nrow(circ_ABl) - dfTreat <- 3 - dfPrep <- 1 - dfReg <- 1 - dfnonP <- 1 - dfRMSE <- c(lenCirc-3-1) - dfTotal <- lenCirc-1 - dfPureE <- lenCirc-6 - dfNonLin <- dfRMSE-dfPureE - - RSS <- sum(resid(lm(readout ~ log_dose*isSample, circ_ABl))^2) - MSE <- RSS/dfRMSE - SSE <- sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2) - MSpure <- SSE/dfPureE - - if (PureErrFlag) { - FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl)) - meanPureErr <- FitAnova[4,3] - SU_modAB <- tryCatch({ - SU_modAB <- summary(modAB) - }, error = function(msg) { - return(NA) - }) - if (length(SU_modAB)>1) s_modABcoeffs <- summary(modAB)$coefficients - - DFsPure <- FitAnova[4,1] - VCOV <- vcov(modAB) - V_V <- VCOV/SU_modAB$sigma^2 - VCOVpure <- V_V*meanPureErr - SEsPure <- sqrt(diag(V_V)*meanPureErr) - su_mod2[,3] <- SEsPure - su_mod2[,4] <- su_mod2[,2]/su_mod2[,3] - su_mod2[,5] <- 2*(1-pt(abs(su_mod[,4]), FitAnova[4,1])) - - s_mu <- summary(modABu)$coefficients - SU_modABu <- summary(modABu) - VCOVu <- vcov(modABu) - V_Vu <- VCOVu/SU_modABu$sigma^2 - SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) - - su_modU2[,3] <- SEsPureU - su_modU2[,4] <- su_modU2[,2]/su_modU2[,3] - su_modU2[,5] <- 2*(1-pt(abs(su_modU2[,4]), FitAnova[4,1])) - - uCI_SloDiffP <- su_modU[4,1] + qt(0.975,8)*SEsPureU[4] - lCI_SloDiffP <- su_modU[4,1] - qt(0.975,8)*SEsPureU[4] - SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiffP,uCI_SloDiffP) - - SSRes <- SSE - dfRes <- dfPureE - - } else { - SSRes <- RSS - dfRes <- dfRMSE - } - - # treatment - SStreat <- print(sum((predict(lm(readout ~ factor(log_dose)*isSample, circ_ABl))-mean(circ_ABl$readout))^2)) - F_treat <- (SStreat/dfTreat)/(SSRes/dfRes) - # Preparation - SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl))-mean(circ_ABl$readout))^2)) - F_prep <- (SSprep/dfTreat)/(SSRes/dfRes) - # Regression - # ANOVA tape II SS of regression - SSreg <- Anova(lm(readout ~log_dose + isSample, circ_ABl))[1,1] - # Non-parallelism - # diff of RSS of restricted and unrestricted model - SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2) - F_nonpar <- SSnonpar/(sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)/(lenCirc-4)) - - # non-linearity - SSnonlin <- sum((predict(modABu)-predict(lm(readout ~ as.factor(log_dose)*isSample, circ_ABl)))^2) - # = RSS-SSE - # Total SS - SStot <- sum((circ_ABl$readout-mean(circ_ABl$readout))^2) - # Significance of R^2 F-ratio - # MSR/MSE - # sample A - F_R2_A <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2)/ - (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - circ_Al$readout)^2)/(nrow(circ_Al)-3)) - pFR2_A <- round(pf(F_R2_A,1,6),4) - # sample B - F_R2_B <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2)/ - (sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2)/(nrow(circ_Bl)-3)) - pFR2_B <- round(pf(F_R2_B,1,6),4) - # sign of non-lin with pure error: MSSnonlin/MSSE - F_nonlin <- (SSnonlin/2)/(SSE/dfPureE) - - # sign of slope - F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2)/(sum((circ_Bl$readout - predict(modB))^2)/(nrow(circ_Bl)-2)) - F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2)/(sum((circ_Al$readout - predict(modA))^2)/(nrow(circ_Al)-2)) - # F-test on regression: MSSreg/MSSE - if (is.na(F_nonlin)) F_nonlin <- 0 - if (F_nonlin > 0) { - p_F_nonlin <- round(pf(F_nonlin,2,dfPureE, lower.tail = F),5) - } else { p_F_nonlin <- "SSnonlin neg or 0"; } - - # significances - F_regr <- (SSreg/1)/(SSRes/dfRes) - p_F_regr <- round(pf(F_regr,1,dfRes, lower.tail = F),3) - p_F_treat <- round(pf(F_treat,3,dfRes, lower.tail = F),3) - p_F_prep <- round(pf(F_prep,1,dfRes, lower.tail = F),3) - p_F_slope_A <- round(pf(F_slope_A,1,(nrow(circ_Al)-2), lower.tail = F),3) - p_F_slope_B <- round(pf(F_slope_B,1,(nrow(circ_Bl)-2), lower.tail = F),3) - p_F_nonp <- round(pf(F_nonpar,1,dfRes, lower.tail = F),3) - p_F_LoF <- p_F_nonlin - - res_tab_lin <- data.frame(test = c("F-test on sign. of regression", "F_test on non-lin", - "F-test on R^2 A","F_test on R^2 B", - "F-test on slope A","F-test on slope B", - "F-test on non-parallelism","F-test on preparation"), - test_results = c(ifelse(p_F_regr<0.05,0,1),ifelse(p_F_nonlin<0.05,1,0), - ifelse(pFR2_A<0.05,1,0),ifelse(pFR2_B<0.05,1,0), - ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1), - ifelse(p_F_nonp<0.05,1,0),ifelse(p_F_prep<0.05,0,1)), - estimate = c(p_F_regr, p_F_nonlin,pFR2_A,pFR2_B,p_F_slope_A, - p_F_slope_B,p_F_nonp,p_F_prep), - Source = c("Treatment","Preparation","Regression","Non-parallelism", - "Resid Error","Non-linearity","Pure error", "Total"), - df = c(dfTreat,1,1,1,dfRMSE,2,dfPureE,lenCirc-1), - SumSquares = c(round(SStreat,5),round(SSprep,5),round(SSreg,5), - round(SSnonpar,5),round(RSS,5),round(SSnonlin,5), - round(SSE,5),round(SStot,5)), - MS = c(round(SStreat/dfTreat,5),round(SSprep,5),round(SSreg,5), - round(SSnonpar,5),round(RSS/dfRMSE,5),round(SSnonlin/2,5), - round(SSE/dfPureE,5),round(SStot/dfTotal,5)), - "F-value" = c(round(F_treat,5), round(F_prep,5),round(F_regr,5), - round(F_nonpar,5),"",round(F_nonlin,5),"",""), - "p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "","")) - RET <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2) - return(RET) -} - -PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) { - - mLin <- gsl_nls(readout ~ (intS+r)*isSample + intS*isRef + k*log_dose, - data=circle, - start=list(intS = 0, k=1,r=0), - control = gsl_nls_control(xtol=1e-10,ftol=1e-10,gtol=1e-10)) - # alternativ: modAB <- lm(readout ~ log_dose+isSample, circle) - sum_mLin <- summary(mLin) - - log_dose <- unique(all_l2$log_dose) - seq_x <- seq(min(log_dose), max(log_dose),0.1) - SAMPLEtrue <- sigmoid[5] + (sigmoid[7]-sigmoid[5])/(1+exp(sigmoid[6]*((sigmoid[4]-sigmoid[8]-seq_x)))) - REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[2]*((sigmoid[4]-seq_x)))) - - truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue) - - p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) + - geom_point() + - labs(title=paste("linear regression model", indS,indT), color="product") + - scale_colour_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) + - ylim(min(all_l2$readout),max(all_l2$readout)) + - theme_bw() - p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="#4545BA", - inherit.aes = F) + - geom_line(data=pl_df,aes(x=lnC,y=plotT),color="#C2173F", - inherit.aes = F) + - geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, - inherit.aes = F) + - geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4, - inherit.aes = F) + - labs(title = paste("unrestricted linear regression model",indS,indT), color="product") + - theme(legend.position="none", axis.text = element_text(size=14)) - p3 <- p2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef), - size=5,alpha=0.2), inherit.aes = FALSE) + - scale_shape_manual(labels=c("test","reference"), values=c(21,21)) - # fit intercept for test and ref and common slope - - - pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log_dose - pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log_dose - pl_rest <- data.frame(lnC=log_dose, plotS=pl_restS, plotT=pl_restT) - - pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="#4545BA", - inherit.aes = F) + - geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="#C2173F", - inherit.aes = F) + - geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, - inherit.aes = F) + - geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4, - inherit.aes = F) + - labs(title = paste("restricted linear regression model",indS,indT), color="product") + - theme(legend.position="none", axis.text = element_text(size=14)) - pr3 <- pr2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef), - size=5,alpha=0.2), inherit.aes = FALSE) + - scale_shape_manual(labels=c("test","reference"), values=c(21,21)) - return(grid.arrange(p3,pr3,nrow=1)) -} - - - -#' Calculates the potency of 4PL PLA of all models model -#' -#' The gradient method is used for calculating the potency for a restricted model, an unrestricteed model, -#' and the log-transformations thereof. -#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used. -#' -#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(concentration). -#' @param PureErrFlag Boolean if Pure Error should be used. -#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results -#' 2) summary of unrestricted linear model 3) Slope difference +CIs -#' 4) summary of restricted linear model. -#' @export -#' @examples -#' ro_new <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229), -#' R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608), -#' T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791), -#' T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655), -#' log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)) -#' PureErrF <- TRUE -#' -#' -#' pot4plFUNC(ro_new, PureErrF) - - -pot4plFUNC <- function(ro_new, PureErrFlag) { - all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") - isRef <- rep(c(1,0),1,each=nrow(all_l)/2) - isSample <- rep(c(0,1),1,each=nrow(all_l)/2) - all_l$isRef <- isRef - all_l$isSample <- isSample - all_l$Conc <- exp(all_l$log_dose) - all_l$readout[all_l$readout < 0] <- 0.01 - all_l$readouttrans <- log(all_l$readout) -#browser() - CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)]) - if (CORdat<0) SLOPE <- -1 else SLOPE <- 1 - # - FITs <- Fitting_FUNC(ro_new, TransFlag = F) - if (!PureErrFlag) { - pot_est <- FITs[[3]] - potU_est <- FITs[[4]] - } else { - FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) - meanPureErr <- FitAnova[4,3] - # SU_mr <- tryCatch({ - # SU_mr <- summary(mr) - # }, error = function(msg) { - # return() - # }) - SU_mr <- FITs[[1]] - SU_mrCoeff <- SU_mr$coefficients - # if (length(SU_mr)>1) { - # SU_mrCoeff <- SU_mr$coefficients - # } else { SU_mrCoeff <- rep(NA,5) } - #te2$cov.unscaled[1,1]*te2$sigma^2 - #VCOV <- vcov(mr) - #V_V <- VCOV/SU_mr$sigma^2 - V_V <- SU_mr$cov.unscaled - SEsPure <- sqrt(diag(V_V)*meanPureErr) - pot_est <- c(exp(SU_mrCoeff['r',1]),exp(SU_mrCoeff['r',1]-qt(0.975,nrow(all_l)-5)*SEsPure['r']), - exp(SU_mrCoeff['r',1]+qt(0.975,nrow(all_l)-5)*SEsPure['r'])) - # unrestricted - SU_mu <- FITs[[2]] - s_muCoeff <- SU_mu$coefficients - #SU_mu <- summary(mu) - #VCOVu <- vcov(mu) - V_Vu <- SU_mu$cov.unscaled - SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) - potU_est <- c(exp(s_muCoeff['r',1]),exp(s_muCoeff['r',1]-qt(0.975,nrow(all_l)-8)*SEsPureU['r']), - + exp(s_muCoeff['r',1]+qt(0.975,nrow(all_l)-8)*SEsPureU['r'])) - } # PureErrFlag - - FITstrans <- Fitting_FUNC(ro_new, TransFlag = TRUE) - - - # pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic")) - # potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic")) - pot_est_log <- FITstrans[[3]] - potU_est_log <- FITstrans[[4]] - colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI") - colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI") - #browser() - su_mr_log <- FITstrans[[1]] #summary(mr_log) - Dat$RMSE_Rlog <- su_mr_log$sigma - su_mu_log <- FITstrans[[2]] # summary(mu_log) - Dat$RMSE_Ulog <- su_mu_log$sigma - Dat$up_lowAslog <- su_mu_log$coefficients[1,1] - su_mu_log$coefficients[4,1] - potALL <- rbind(round(pot_est,6), round(potU_est,6), round(pot_est_log,6), round(potU_est_log,6)) - - potALL2 <- cbind(c("restricted","unrestricted","ln-transformed restr","ln-transformed unrestr"), potALL) - return(potALL2) -} - -#' Calculates logarithmized CIs for ratios -#' -#' Ratio CIs can be unbound, and not calulatable. Therefore the ratios are logarithmized. -#' The Covariance is not used, as it is 0 for comparisons between products. -#' -#' @param xs The estimate of the reference sample. -#' @param xt The estimate of the test sample. -#' @param se_xs The standard error of the estimate of the reference sample. -#' @param se_xt The standard error of the estimate of the test sample. -#' @param CoVar The covariance between test and reference estimate (set to 0). -#' @param DFs Degrees of freedom, either from pure error term (no of readouts - no of concentrations*2), or RMSE term (no of readouts - 8 parameters). -#' @param Conf The confidence of the Confidence Interval (default 95% which requires 0.975). -#' @returns A data-frame with the lower and upper CI in anti-log form. -#' @export -#' @examples -#' xs=2; xt=3.2; se_xt=0.34;se_xs=0.23; DFs=32-16 -#' -#' -#' ParamCI_F(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) - - - -ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) { - log_xs <- log(abs(xs)) - log_xt <- log(abs(xt)) - var_log_xs <- (se_xs/xs)^2 # approximate variance of log(xs) - var_log_xt <- (se_xt/xt)^2 - se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVar/(xs*xt) - - lower_log_ratio <- log_xt-log_xs - qt(Conf,DFs)*se_log_ratio - upper_log_ratio <- log_xt-log_xs + qt(Conf,DFs)*se_log_ratio - ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) - return(ci_ratio) -} - - -#' Calculates EQ tests and F-Tests for 4P -#' -#' CIs of asmptote ratios,slope ratios, asympt-difference ratio, and lower asympt-diff. calculated. -#' In addition the F-test on significance of regression and non-linearity are calculated. -#' -#' @param ro_new Data frame with the concentrations and the respective readouts. -#' @param Lim The limits to check if the relative CI is within the limits. -#' @param PureErrFlag Boolean if Pure Error should be used. -#' @returns A data-frame with the lower and upper CI in anti-log form. -#' @export -#' @examples -#' -#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), -#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), -#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), -#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) -#' Lim <- c(-1, 1, 0.005 , 2, 0.5, 2, 0.5,2,75,133, 75, 133) -#' PureErrF <- FALSE -#' -#' tests_FUNC(ro_new, Lim,PureErrF) - - -tests_FUNC <- function(ro_new, Lim, PureErrFlag) { - - all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") - isRef <- rep(c(1,0),1,each=nrow(all_l)/2) - isSample <- rep(c(0,1),1,each=nrow(all_l)/2) - all_l$isRef <- isRef - all_l$isSample <- isSample - all_l$Conc <- exp(all_l$log_dose) - all_l$readout[all_l$readout < 0] <- 0.01 - #browser() - FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) - POTr_CI <- FITs[[3]][2:3] - potAll2 <- FITs[[3]] - - smr <- FITs[[1]] - smu <- FITs[[2]] - FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) - # pure error - pureSS <- FitAnova[4,2] - pureSS_df <- FitAnova[4,1] - meanPureErr <- FitAnova[4,3] - #vcovMU <- vcov(mu) - V_V <- smu$cov.unscaled - SEsPure <- sqrt(diag(V_V)*meanPureErr) - VCOVpure <- V_V*meanPureErr - DFsPure <- FitAnova[4,1] - - - - testPOTr <- logical() - if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100Lim[[9]]/100 & potAllU2[3] 0) { - p_F_nonlin <- round(pf(F_nonlin, AnovaDFs[6], ERR_df, lower.tail = F),5) - } else { p_F_nonlin <- "SSnonlin neg or single dilutions" } - - test_a <- test_b <- test_d <- test_ad <- logical() - - RSS_r <- round(sum(smr$residuals^2),5) - MSE_r <- RSS_r/(nrow(all_l)-5) - RMSE_r <- round(sqrt(MSE_r),6) - Dat$RMSE_r <- RMSE_r - Dat$RMSE_pure <- RMSE_pure - Dat$RMSE_unr <- round(RMSEunr,6) - - coeffs <- smu$coefficients[,1] - #browser() - ## EQ test on lower As diff - as <- coeffs["as"] - at <- coeffs["at"] - lAs_diff <- (at-as) - uCI_laDiff <- lAs_diff+qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2) - lCI_laDiff <- lAs_diff-qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2) - if (uCI_laDiff < Lim[[2]] & lCI_laDiff > Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1 - - #### EQ test on upper asymptote ratio ---- - # as <- coeffs["as"] - # at <- coeffs["at"] - # uAsRatio <- compParm(potU, "a","/",display=F) - # uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2]) - #browser() - ds <- smu$coefficients["ds",1] - dt <- smu$coefficients["dt",1] - if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds","ds"]) else se_ds <- smu$coefficients["ds",2] - if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt","dt"]) else se_dt <- smu$coefficients["dt",2] - if (PureErrFlag) CoVarlog_d <- VCOVpure["dt","ds"] else CoVarlog_d <- vcovMU["dt","ds"] - if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-8 - uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.975) - if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1 - estUppA <- round(at/as,5) - - Dat$uAsCI <- uAsCI2 - - #### EQ test on slope ratio ---- - # bs <- coeffs["bs"] - # bt <- coeffs["bt"] - # slopeRatio <- compParm(potU, "b","/",display=F) - # slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2]) - - bs <- smu$coefficients["bs",1] - bt <- smu$coefficients["bt",1] - if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs","bs"]) else se_bs <- smu$coefficients["bs",2] - if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt","bt"]) else se_bt <- smu$coefficients["bt",2] - if (PureErrFlag) CoVarlog_b <- VCOVpure["bt","bs"] else CoVarlog_b <- vcovMU["bt","bs"] - slopeCI2 <- ParamCI_F(bt,bs,se_bt, se_bs,CoVarlog_b, DFs, Conf=0.975) - if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1 - estUppA <- round(at/as,5) - - Dat$slopeRatioCI <- slopeCI2 - - #### EQ test on lower As ratio ---- - - # lAsRatio <- compParm(potU, "d","/",display=F) - # slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2]) - - as <- smu$coefficients["as",1] - at <- smu$coefficients["at",1] - if (PureErrFlag) se_as <- sqrt(VCOVpure["as","as"]) else se_as <- smu$coefficients["as",2] - if (PureErrFlag) se_at <- sqrt(VCOVpure["at","at"]) else se_at <- smu$coefficients["at",2] - if (PureErrFlag) CoVarlog_a <- VCOVpure["at","as"] else CoVarlog_a <- vcovMU["at","as"] - lAsCI2 <- ParamCI_F(at,as,se_at, se_as,CoVarlog_a, DFs, Conf=0.975) - if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1 - estLowA <- round(at/as,5) - - Dat$lAsCI <- lAsCI2 - - #### EQtest on ratio of As difference ---- - AsDiffRatio <- (dt-at)/(ds-as) - - dt_at <- (dt-at) - ds_as <- (ds-as) - se_ds_asPure <- sqrt(VCOVpure["as","as"]+VCOVpure["ds","ds"]-2*VCOVpure["as","ds"]) - se_dt_atPure <- sqrt(VCOVpure["at","at"]+VCOVpure["dt","dt"]-2*VCOVpure["at","dt"]) - se_ds_asRMSE <- sqrt(vcovMU["as","as"]+vcovMU["ds","ds"]-2*vcovMU["as","ds"]) - se_dt_atRMSE <- sqrt(vcovMU["at","at"]+vcovMU["dt","dt"]-2*vcovMU["at","dt"]) - if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE - if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE - - AsDiffCI2 <- ParamCI_F(dt_at,ds_as,se_dt_at, se_ds_as,CoVar=0, DFs, Conf=0.975) - if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1 - estLowA <- round(at/as,5) - - Dat$up_lowAs <- abs(ds-as) - - lowerCIlowerA <- lAsCI2[1]; lowerCIupperA <- uAsCI2[1]; upperCIlowerA <- lAsCI2[2]; upperCIupperA <- uAsCI2[2] - test_lowA <- test_d; test_uppA <- test_a - #browser() - res_tab <- data.frame(test= c("F-test on sign. of regression*", - "EQ test on lower asymptotes difference", - "EQ test ratio of lower asymptotes", - "EQ test ratio of Hill slopes", - "EQ test ratio of upper asymptotes", - "F-test on non-linearity*", - "EQ test ratio of asymptote difference", - "geom. rel. CI restr. model", - "geom. rel. CI unrestr. model"), - test_results = c(ifelse(p_F_regr<0.05,0,1), test_la_diff, test_lowA, test_b, test_uppA, - ifelse(p_F_nonlin>1,1, ifelse(p_F_nonlin<0.05,1,0)), test_ad, - testPOTr, test_c), - estimate = c(round(p_F_regr, 3), round(lAs_diff, 5), - estLowA, round(bs/bt,5), estUppA, p_F_nonlin, - round(dt_at/ds_as, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)), - lower_limit = c("-",Lim[[1]],Lim[[3]],Lim[[5]],Lim[[7]],"-",Lim[[11]],Lim[[9]],Lim[[9]]), - upper_limit = c("-",Lim[[2]],Lim[[4]],Lim[[6]],Lim[[8]],"-",Lim[[12]],Lim[[10]],Lim[[10]]), - lower_CI = c(RMSE_r, round(lCI_laDiff,3), round(lAsCI2[1],5), round(slopeCI2[1],5), - round(uAsCI2[1],5), "-", round(AsDiffCI2[1],5), round(potAll2[2],2), round(potAllU2[2],2)), - upper_CI = c(RMSE_pure, round(uCI_laDiff,3), round(lAsCI2[2],5), round(slopeCI2[2],5), - round(uAsCI2[2],5), "-", round(AsDiffCI2[2],5), round(potAll2[3],2), round(potAllU2[3],2)) - ) - return(res_tab) -} - -#' ANOVA unrestricted model -#' -#' ANOVA with unrestricted model. -#' -#' -#' @param ro_new Data frame with the concentrations and the respective readouts. -#' @returns A data-frame with the analysis of variance. -#' @export -#' @examples -#' -#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), -#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), -#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), -#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) -#' -#' -#' ANOVA4plUnresfunc(ro_new) - - -ANOVA4plUnresfunc <- function(ro_new) { - all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout") - all_len <- nrow(all_l) - isRef <- rep(c(1,0),1,each=all_len/2) - isSample <- rep(c(0,1),1,each=all_len/2) - all_l$isRef <- isRef - all_l$isSample <- isSample - all_l$Conc <- exp(all_l$log_dose) - all_l$readout[all_l$readout < 0] <- 0.01 - - FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) - smr <- FITs[[1]] - smu <- FITs[[2]] - smrPREDs <- FITs[[5]] - smuPREDs <- FITs[[6]] - # pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - # pmodels=data.frame(1,1,1,isSample)) - # potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")), - # pmodels=data.frame(isSample, isSample,isSample,isSample)) - SStreat <- round(sum((smuPREDs - mean(all_l$readout))^2),5) - SStreat_df <- length(unique(all_l$log_dose))-1 - SSregr <- round(sum((smrPREDs-mean(all_l$readout))^2),5) - ## Non-parallel - SSnonparallel <- round(sum(smr$residuals^2) - sum(smu$residuals^2),5) - ## Preparation - SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5) - ## Resid Err - RSS <- round(sum(smu$residuals^2),5) - RSS_df <- nrow(all_l)-SStreat_df-1 - FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) - # PureErr - SSE <- round(FitAnova[4,3],5) - SSE_df <- FitAnova[4,1] - # Non-Linearity - SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - smuPREDs)^2),4) - LoF_df <- FitAnova[1,1]+FitAnova[2,1] - ## Total - SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5) - MSE <- RSS/RSS_df - noConc <- length(unique(all_l$Conc)) - AnovaDFs <- c(noConc-1, 1,3,noConc-4-1, nrow(all_l)-noConc, noConc, nrow(all_l)-noConc-noConc, nrow(all_l)-1) - p_SStreat <- round(pf((SStreat/AnovaDFs[1])/MSE, AnovaDFs[1],RSS_df, lower.tail = F),3) - p_SSprep <- round(pf((SSprep/AnovaDFs[2])/MSE, AnovaDFs[2],RSS_df, lower.tail = F),3) - p_SSregr <- round(pf((SSregr/AnovaDFs[3])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) - p_SSnonp <- round(pf((SSnonparallel/AnovaDFs[4])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3) - p_SSLoF <- round(pf((SSnonlin/LoF_df)/(SSE/SSE_df), LoF_df,SSE_df, lower.tail = F),5) - - ANOVAtab <- data.frame(Source = c("Treatment","Preparation","Regression", - "Non-Parallelism","Residual Error","Non-linearity", - "Pure Error","Total"), - DF = AnovaDFs, - SumSquares = c(SStreat, SSprep,SSregr, SSnonparallel, - RSS, SSnonlin,SSE, SStot), - MeanSquares = c(round(SStreat/AnovaDFs[1],3), SSprep, round(SStreat/AnovaDFs[3],3),round(SSnonparallel/AnovaDFs[4],3), - round(MSE,5), round(SSnonlin/LoF_df,5), round(SSE/SSE_df,5),""), - "F-value" = c(round((SStreat/AnovaDFs[1])/MSE,5), round((SSprep/AnovaDFs[2])/MSE,5), - round((SSregr/AnovaDFs[3])/MSE,5),round((SSnonparallel/AnovaDFs[4])/MSE,5), - "",round((SSnonlin/LoF_df)/(SSE/SSE_df),5),"",""), - "p_value" = c(round(p_SStreat,3), p_SSprep, round(p_SSregr,3), p_SSnonp,"",p_SSLoF,"","") - ) - - return(ANOVAtab) -} - -#' Calculates descriptive statistics per concentration -#' -#' Mean, SD, and CV% per concentration and product.. -#' -#' -#' @param ro_new Data frame with the concentrations and the respective readouts. -#' @returns A data-frame with the concentrations and metadata. -#' @export -#' @examples -#' -#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497), -#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665), -#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591), -#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764)) -#' -#' -#' perConcTab(ro_new, noDilSeries=3) - - - -perConcTab <- function(ro_new, noDilSeries) { - Reftab <- ro_new[,c(1:noDilSeries)] - Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))] - tReftab <- t(Reftab) - colnames(tReftab) <- round(ro_new[,ncol(ro_new)],5) - - avs <- apply(tReftab,2,mean) - sds <- apply(tReftab,2,sd) - cv <- sds/avs*100 - tReftab2 <- rbind(tReftab, avs,sds,cv) - - tTesttab <- t(Testtab) - colnames(tTesttab) <- round(ro_new[,ncol(ro_new)],5) - - avs_test <- apply(tTesttab,2,mean) - sds_test <- apply(tTesttab,2,sd) - cv_test <- sds_test/avs_test*100 - tTesttab2 <- rbind(tTesttab, avs_test,sds_test,cv_test) - concTab <- rbind(tReftab2, tTesttab2) - return(concTab) - -} - -#' Calculates dilution series. -#' -#' Recursive function to calcuate a geometric dilution series. -#' -#' @param x Starting concentration. -#' @param Div Divisor to get next concentration. -#' @param N Number of dilution step, increases. -#' @param res Vector of results. -#' @param noDil Final number of concentrations -1 (the initial one). -#' @returns A data-frame with the concentrations and metadata. -#' @export -#' @examples -#' -#' x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7 -#' -#' divFUN(x,Div,N,res,noDil) - - - -divFUN <- function(x,Div,N,res,noDil) { - N <- N+1 - y <- x/Div - res <- c(res,y) - if (N==noDil) { return(res) } - divFUN(y,Div,N,res,noDil) -} - +source("Global.R") #### ui ---- @@ -1396,8 +201,8 @@ server <- function(input, output, session) { tabsetPanel(id="tabs", tabPanel("linear PLA", column(12, - htmlOutput("PureErrW3"), - tags$head(tags$style("#PureErrW3{color: red; + htmlOutput("PureErrWLinXL"), + tags$head(tags$style("#PureErrWLinXL{color: red; font-size: 16px; font_style: italic;}")), plotOutput("plotLin"), @@ -1421,12 +226,17 @@ server <- function(input, output, session) { DTOutput("ANOVAlin")) ), tabPanel("Report", - h4("Settings for report") + h4("Settings for report"), + )) ) ) ), tabPanel("parameter estimates", + htmlOutput("PureErrWParEst"), + tags$head(tags$style("#PureErrWParEst{color: red; + font-size: 16px; + font_style: italic;}")), column(3,style = "background: #4FCBD922", br(), h4("Regression results restricted"), @@ -1446,8 +256,10 @@ server <- function(input, output, session) { ), tabPanel("Report", h4("Settings for report"), - downloadButton("downloadXLReport", label="Download PDF report", class="butt"), - tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"), + downloadButton("downloadXLReport", label="Download 4PL PDF report", class="butt"), + tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: Courier New}"), + downloadButton("downloadXLReportLin", label="Download linear PLA PDF report", class="butt"), + tags$style(type="text/css","#downloadXLReportLin {background-color: #4FCBD9; color: black;font-family: Courier New}"), ) ) ) @@ -1478,6 +290,16 @@ server <- function(input, output, session) { # numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")), # bsPopover(id="q4", title="", content="The calculated limits ...")), #h5("Diagnostics only shown if EXCEL is uploaded"), + column(2, + h3("Settings"), + helpText("Vary the slider to see the effect of special cause"), + sliderInput("sdfac","Variability as % of lower to upper asymptote", max=10, value=3, min=0.1, step=0.1), + checkboxInput("heterosked","heteroskedastic noise", FALSE), + sliderInput("potencydiff","potency of test (%)", max=200, min=50, value=100, step=1), + # sliderInput("outlL","outlier in lower asymptote", min=0, max=1.5, value=0, step=0.1), + # sliderInput("outlM","outlier in mid part", min=0, max=1.5,value=0, step=0.1), + # sliderInput("outlU","outlier in upper asymptote", min=0, max=1.5,value=0, step=0.1) + ), column(2,style = "background: #7FAEFF", #actionButton("StartCalc", "Click, when calculations to be started"), h4("curve settings"), @@ -1521,35 +343,26 @@ server <- function(input, output, session) { verbatimTextOutput("dilutions") ), - column(2, - h3("Settings"), - helpText("Vary the slider to see the effect of special cause"), - sliderInput("sdfac","Variability as % of lower to upper asymptote", max=10, value=3, min=0.1, step=0.1), - checkboxInput("heterosked","heteroskedastic noise", FALSE), - sliderInput("potencydiff","potency of test (%)", max=200, min=50, value=100, step=1), - sliderInput("outlL","outlier in lower asymptote", min=0, max=1.5, value=0, step=0.1), - sliderInput("outlM","outlier in mid part", min=0, max=1.5,value=0, step=0.1), - sliderInput("outlU","outlier in upper asymptote", min=0, max=1.5,value=0, step=0.1) - ), + h4("log-dilutions from settings above"), - verbatimTextOutput("logdil"), + column(8, box(title = "Simulated data per log-concentration", status="warning",solidHeader = T, width=12, "incl. mean, sd and CV%", - DT::dataTableOutput("ConctabMeta"))) + DT::dataTableOutput("ConctabMeta")), + verbatimTextOutput("logdil")) #) ), #### 4pl fits ---- tabPanel("4pl-fit", - tags$head(tags$style("#PureErrW2{color: red; - font-size: 10px; - font_style: italic;}")), + wellPanel( fluidRow( column(10, - tags$style(span(htmlOutput("PureErrW3"), style="color: red")), - htmlOutput("PureErrW3"), - + htmlOutput("PureErrW4plMeta"), + tags$head(tags$style("#PureErrW4plMeta{color: red; + font-size: 16px; + font_style: italic;}")), plotOutput("plot4plMeta", width = "80%"), DTOutput("pottab4pl"), "Footnote: test performed on relative CIs.", @@ -1578,13 +391,18 @@ server <- function(input, output, session) { ), tabPanel("ln-transformed y", + htmlOutput("PureErrWLogMeta"), + tags$head(tags$style("#PureErrWLogMeta{color: red; + font-size: 16px; + font_style: italic;}")), h4("ln-transformed y-axis plots"), plotOutput("plot4plTransMeta", width = "80%"), DT::dataTableOutput("pottab4plTrans"), ), tabPanel("linear regression", - h4("Evaluations from meta-data"), + htmlOutput("PureErrW3"), + h4("Evaluations from meta-data"), tags$head(tags$style("#PureErrW3{color: red; font-size: 16px; font_style: italic;}")), @@ -1689,7 +507,7 @@ server <- 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, + sig <- c(input$lowAsymptREF, input$lowAsymptTEST,input$uppAsymptREF,input$uppAsymptTEST, input$slopeREF,input$slopeTEST,input$EC50,input$potDiff) sig }) @@ -1769,6 +587,11 @@ server <- function(input, output, session) { ifelse(PureErrFlag, 'Pure Error is selected', '') }) output$PureErrW2 <- renderText(warning_text2()) + + warning_textParEst <- reactive({ + ifelse(PureErrFlag, 'Pure Error is selected', '') + }) + output$PureErrWParEst <- renderText(warning_textParEst()) noDilSeries <-(ncol(XLdat2)-1)/2 noDils <- nrow(XLdat2) @@ -1808,16 +631,17 @@ server <- function(input, output, session) { DFsPure <- FitAnova[4,1] #VCOV <- vcov(mr) V_V <- Smr$cov.unscaled #VCOV/Smr$sigma^2 - VCOVpure <- V_V*meanPureErr + #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])) + #browser() #VCOVu <- vcov(mu) V_Vu <- Smu$cov.unscaled #VCOVpure <- V_Vu*meanPureErr SEsPureU <- sqrt(diag(V_Vu)*meanPureErr) - potU_est <- data.frame(estimate=exp(coeffsMU[7]), lowerCI = exp(coeffsMU[7]-qt(0.975,DFsPure)*SEsPureU[7]), - upperCI = exp(coeffsMU[7]+qt(0.975,DFsPure)*SEsPureU[7])) + potU_est <- data.frame(estimate=exp(coeffsMU[8]), lowerCI = exp(coeffsMU[8]-qt(0.975,DFsPure)*SEsPureU[8]), + upperCI = exp(coeffsMU[8]+qt(0.975,DFsPure)*SEsPureU[8])) } @@ -1862,6 +686,7 @@ server <- function(input, output, session) { }) + ANOVAtab2 <- ANOVA4plUnresfunc(ro_new = XLdat2) output$ANOVAXLS <- renderTable({ ANOVAtab2 }) @@ -2218,16 +1043,28 @@ server <- function(input, output, session) { }) - #### Plot 4pl ---- + #### Plot 4pl Meta ---- output$plot4plMeta <- renderPlot({ - #browser() + PureErrFlag <- input$PureErrMeta + warning_text3 <- reactive({ + ifelse(PureErrFlag, 'Pure error selected','') + }) + + output$PureErrW4plMeta <- renderText(warning_text3()) + sigmoid <- sigmoid() det_sig=NULL plot_f(sim2(), TransFlag = F) }) - #### Plot 4pl Transformed ---- + #### Plot 4pl Meta Transformed ---- output$plot4plTransMeta <- renderPlot({ + PureErrFlag <- input$PureErrMeta + warning_text3 <- reactive({ + ifelse(PureErrFlag, 'Pure error selected','') + }) + + output$PureErrWLogMeta <- renderText(warning_text3()) #browser() sigmoid <- sigmoid() det_sig=NULL @@ -2235,7 +1072,7 @@ server <- function(input, output, session) { }) - #### Testergebnisse für 4PL ---- + #### Testergebnisse für 4PL Meta ---- observe({ if (is.null(sim2())) return(NULL) if (is.null(input$PureErrMeta)) return(NULL) @@ -2320,7 +1157,7 @@ server <- function(input, output, session) { }) # observe - ####plot CIs ---- + #### plot CIs XL---- observe({ tab <- Dat$tests_FUNC if (is.null(tab)) return(NULL) @@ -2342,7 +1179,7 @@ server <- function(input, output, session) { REP$CIplot <- p_CIs }) - + #### simulated data tab Meta ---- output$simdat <- DT::renderDataTable({ tab <- sim2() if (is.character(tab)) stop(tab) @@ -2529,7 +1366,7 @@ server <- function(input, output, session) { ifelse(PureErrFlag, 'Pure error is selected','') }) output$PureErrW <- renderText(warning_text()) - browser() + #browser() LIN <- ANOVAlintests(tab,circlesMeta,Limite,PureErrFlag=PureErrFlag) df <- LIN[[1]] su_modU <- LIN[[2]] @@ -2700,13 +1537,19 @@ server <- function(input, output, session) { circles <- Dat$circles PureErrFlag <- input$PureErr + + warning_text2 <- reactive({ + ifelse(PureErrFlag, 'Pure Error is selected', '') + }) + output$PureErrWLinXL <- renderText(warning_text2()) + pottab <- LinPotTab(circles,Lim,PureErrFlag = PureErrFlag) #browser() dat <- datatable(pottab, options=list( dom="t",rownames=F )) %>% formatStyle("test_result", target='row', - backgroundColor = styleEqual(c(0,1), c("lightgrey","#F9545488"))) + backgroundColor = styleEqual(c(0,1), c("#B5C74055","#F9545488"))) }) @@ -2718,7 +1561,7 @@ server <- function(input, output, session) { as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), as.numeric(input$lowerPot), as.numeric(input$upperPot)) - + #browser() circles <- Dat$circlesMeta PureErrFlag <- input$PureErrMeta pottab <- LinPotTab(circles,Lim,PureErrFlag = PureErrFlag) @@ -2727,7 +1570,7 @@ server <- function(input, output, session) { options=list( dom="t",rownames=F )) %>% formatStyle("test_result", target='row', - backgroundColor = styleEqual(c(0,1), c("lightgrey","#F9545488"))) + backgroundColor = styleEqual(c(0,1), c("#B5C74055","#F9545488"))) }) #### 4pl potency table Meta ---- @@ -3211,7 +2054,7 @@ server <- function(input, output, session) { } }) - #### download XL report---- + #### download XL 4PL report---- output$downloadXLReport <- downloadHandler( filename= paste0("Report_4PLEvaluation", Dat$FileName,".pdf"), @@ -3234,7 +2077,28 @@ server <- function(input, output, session) { } ) + #### download XL Lin report---- + output$downloadXLReportLin <- downloadHandler( + filename= paste0("Report_linPLA",Dat$nameRep ,".pdf"), + + content = function(file) { + tpdr <- tempdir() + tempReport <- file.path(tpdr,"Doc_BioassayLinReport.Rmd") + file.copy("Doc_BioassayLinReport.Rmd", tempReport, overwrite = T) + + tempReportc <- file.path(tpdr,"logo.png") + file.copy("logo.png", tempReportc, overwrite = T) + + rmarkdown::render(tempReport, output_file = file, + params = list(FileName = Dat$FileName, + author = Dat$author, + REP = REP, + coeffs = Dat$coeffs_UN), + envir = new.env(parent = globalenv())) + + } + ) } shinyApp(ui, server)