################################################################################ # All functions for plateflow ################################################################################ #' 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 <- tryCatch({gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))), data=all_l2, start=startlist,#race=T, control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6)) }, warning = function(e) { mr <<- "In nlsModel singular gradient matrix" }) # Stop if singular gradient matrix if (is.character(mr)) return(mr) 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 when the fitting has a singularity #' #' Returns the scatter plot of the data, with REF in blue and TEST in red. #' One plots are returned. #' #' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration. #' @returns A grid object with 2 linearity plots, restricted and unrestricted model. #' @export #' @examples #' data.frame(R_dil1 = c(10.0651024695491, 10.9844983291817, 10.7635586089293, 10.4597656321327, 10.3898668457823, 10.8171761349909, #' 10.319758021908, 10.1304854046653), #' R_dil2 = c(10.9649145494504, 10.0202868589385, 10.8424145955735, 10.9311360356894, 10.3284659026404, #' 10.6890147558796, 10.3014450252305, 10.9594838595181), #' R_dil3 = c(10.4630510824383, 10.4566715089363, 10.2350765290036, 10.3300581874798, 10.9648088137065, #' 10.286893755805, 10.4856643841389, 10.5275521552307), #' T_dil1 = c(12.732175566336, 12.7756403995095, 12.1672539684741, 12.7060603907892, 12.8000685682832, #' 12.8800092157515, 12.7160581291873, 12.6996878912416), #' T_dil2 = c(12.3923194313831, 12.0943488144175, 12.7955302154828, 12.4825917078735, 12.6856540203788, #' 12.7348548498556, 12.9222470610476, 12.1186618671252), #' T_dil3 = c(12.7899182255274, 12.9722600411128, 12.7078445380891, 12.4913523531941, 12.1718281909609, #' 12.5313873615133, 12.952802332772, 12.5960321394342), #' log_dose = c(0, -1.09861228866811, -2.19722457733622, -3.29583686600433, -4.39444915467244, #' -5.49306144334055, -6.59167373200866, -7.69028602067677)) #' #' #' p <- plotSingularity(dat) #' print(p) #' plotSingularity <- function(dat) { #sigmoid,det_sig, CORdat <- cor(dat[,1],dat[,ncol(dat)]) #browser() all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout") isRef <- rep(c(1,0),1,each=nrow(all_l)/2) isSample <- rep(c(0,1),1,each=nrow(all_l)/2) all_l2 <- cbind(all_l, isRef, isSample) #browser() log_dose <- unique(all_l$log_dose) seq_x <- seq(min(log_dose),max(log_dose),0.1) #browser() #all_l2$readout[all_l2$readout < 0] <- 0.01 all_l2$readouttrans <- log(all_l2$readout) #browser() pSing <- 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("No 4pl fit possible"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + scale_shape_manual(labels=c("test","reference")) + scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) + #theme_bw() + theme(axis.text = element_text(size=14)) return(pSing) } #' 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 model"), color="product") + scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) + scale_shape_manual(labels=c("test","reference")) + scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) + #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="#4545BABB",linetype=3) + geom_vline(xintercept=c(XasymplT, XasympuT), col="#C2173FBB",linetype=3) + 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 4pl 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) if (!is.null(sigmoid)) { 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) } else { SAMPLEtrue <- NULL REFtrue <- NULL truePL_df <- NULL } 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)) + scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) + 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) + {if (!is.null(truePL_df)) geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, inherit.aes = F) } + {if (!is.null(truePL_df)) 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), col=c("black"), 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) + {if (!is.null(truePL_df)) geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4, inherit.aes = F) } + {if (!is.null(truePL_df)) 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("black"), 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) if (is.character(FITs)) return(FITs) # if singularity 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) }