From 1b3af203a2afdb0ca8f82b638523809ed98d657e Mon Sep 17 00:00:00 2001 From: Simon Innerbichler Date: Tue, 19 May 2026 11:21:39 +0200 Subject: [PATCH] formatted Global.R --- R/Global.R | 1840 +++++++++++++++++++++++++++++----------------------- 1 file changed, 1015 insertions(+), 825 deletions(-) diff --git a/R/Global.R b/R/Global.R index 7d74982..06fad47 100644 --- a/R/Global.R +++ b/R/Global.R @@ -14,95 +14,130 @@ #' 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) +#' 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) +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 + # 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" - }) + 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 - }) + 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)) + 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) }) + 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) }) + 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")) + 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")) + pot_est <- exp(confintd(mrT, "r", method = "asymptotic")) + potU_est <- exp(confintd(muT, "r", method = "asymptotic")) PRED <- predict(mrT) PREDu <- predict(muT) } @@ -113,70 +148,85 @@ Fitting_FUNC <- function(ro_new, TransFlag=F) { #' 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. +#' 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)) +#' 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) +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() - + # 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 + 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")) + + # 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)) - + # 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. +#' 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. @@ -185,151 +235,169 @@ plotSingularity <- function(dat) { #sigmoid,det_sig, #' @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) +#' 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) +#' 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) +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) + # 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] - + 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))) - + 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() + + 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)) + 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")) + + # 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) + + # 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") + 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")) + + 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) + + 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)) + 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))) - + 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() + 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] + # 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] @@ -340,57 +408,65 @@ plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig, # 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))) + 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")) + + # 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 <- 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(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))) + + 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) + + 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) + 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. +#' 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. @@ -401,37 +477,47 @@ plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig, #' @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) +#' 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) +#' 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))) + 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) + 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 <- 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") + 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") + colnames(ro_jit2) <- c("R_dil1", "R_dil2", "T_dil1", "T_dil2", "log_dose") } return(ro_jit2) } @@ -439,7 +525,7 @@ Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct #' 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. +#' 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. @@ -447,66 +533,75 @@ Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct #' @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 +#' 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(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,] + 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) - }) + 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] + 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) + 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 + 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 + 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 + 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() + # 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]] & relLinpotCI[3] < Lim[[10]]) test_potCI <- 0 else test_potCI <- 1 + + pottab <- cbind( + round(linPot * 100, 3), round(ExpLinPot[2] * 100, 3), round(ExpLinPot[3] * 100, 3), + round(test_potCI, 3), round(relLinpotCI[2], 3), round(relLinpotCI[3], 3) + ) + colnames(pottab) <- c("Potency", "lower 95%CI", "upper 95%CI", "test_result", "lowerRel95%CI", "upperRel95%CI") return(pottab) } @@ -514,198 +609,220 @@ LinPotTab <- function(circles, Lim, PureErrFlag) { #' Calculate the ANOVA of linear PLA and restricted and unrestr. model summaries #' #' 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. +#' 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 1) data.frame with F-tests ANOVA and test results +#' @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 -#' dat <- 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)) -#' CIRC <- data.frame(log_dose=c( -2.590267, -2.590267, -3.2834 , -3.2834, -3.97656, -3.97656, -2.590267, -2.590267,-3.2834, -3.2834, -3.97656, -3.97656), -#' replname= c("R_dil1","R_dil2","R_dil1", "R_dil2","R_dil1","R_dil2", "T_dil1","T_dil2","T_dil1", "T_dil2","T_dil1","T_dil2"), -#' readout = c(31993, 31925, 49336, 49003 , 68332, 68034 , 33915, 34007, 52131, 52179 , 70617, 69962), -#' isRef=c(rep(1,6), rep(0,6)), -#' isSample = c(rep(0,6), rep(1,6))) -#' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant +#' dat <- 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) +#' ) +#' CIRC <- data.frame( +#' log_dose = c(-2.590267, -2.590267, -3.2834, -3.2834, -3.97656, -3.97656, -2.590267, -2.590267, -3.2834, -3.2834, -3.97656, -3.97656), +#' replname = c("R_dil1", "R_dil2", "R_dil1", "R_dil2", "R_dil1", "R_dil2", "T_dil1", "T_dil2", "T_dil1", "T_dil2", "T_dil1", "T_dil2"), +#' readout = c(31993, 31925, 49336, 49003, 68332, 68034, 33915, 34007, 52131, 52179, 70617, 69962), +#' isRef = c(rep(1, 6), rep(0, 6)), +#' isSample = c(rep(0, 6), rep(1, 6)) +#' ) +#' Lim <- c(rep(0, 8), 70, 130) # only Lim 9 and 10 relevant #' PureErrF <- TRUE -#' -#' +#' +#' #' ANOVAlintests(ro_new, circles, Lim, PureErrF) - - ANOVAlintests <- function(ro_new, circles, 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 <- 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_lA <- all_l[all_l$isSample == 1,] # TEST - all_lB <- all_l[all_l$isSample == 0,] # REF - #browser() + all_lA <- all_l[all_l$isSample == 1, ] # TEST + all_lB <- all_l[all_l$isSample == 0, ] # REF + # browser() circ_ABl <- circles - circ_Al <- circ_ABl[circ_ABl$isSample ==1,] - circ_Bl <- circ_ABl[circ_ABl$isSample ==0,] + circ_Al <- circ_ABl[circ_ABl$isSample == 1, ] + circ_Bl <- circ_ABl[circ_ABl$isSample == 0, ] # restr CSSI model modAB <- lm(readout ~ log_dose + isSample, circ_ABl) # unrestr with interact term SSSI - modABu <- lm(readout ~ log_dose + isSample + log_dose*isSample, circ_ABl) + modABu <- lm(readout ~ log_dose + isSample + log_dose * isSample, circ_ABl) modA <- lm(readout ~ log_dose + isSample, circ_Al) modB <- lm(readout ~ log_dose + isSample, circ_Bl) - + log_pot_delta <- deltaMethod(modAB, "isSample/log_dose") - 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 + 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 ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up)) - if (ExpLinPot[2]*100>Lim[9] & ExpLinPot[3]*100>Lim[10]) test_potCI <- 0 else test_potCI <- 1 - + if (ExpLinPot[2] * 100 > 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_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) - + 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 - + 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] + 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])) - + 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) - + 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) + 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) + 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] + 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)) - + 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) + 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) + 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) + 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) + 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) - + 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_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"; } - + 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) + 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) + + 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) } @@ -713,8 +830,8 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) { #' 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. +#' 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. @@ -727,133 +844,169 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) { #' @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) +#' 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)) +#' 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) + 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) + 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)) + + + + 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)) + 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)) -} + 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. +#' 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 +#' @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)) +#' 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 <- 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 - # + # 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] + 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) { @@ -864,42 +1017,46 @@ pot4plFUNC <- function(ro_new, PureErrFlag) { # 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 + # 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'])) + 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) + # 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'])) + 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) + 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) + 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) } @@ -918,22 +1075,23 @@ pot4plFUNC <- function(ro_new, PureErrFlag) { #' @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) { +#' 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 + 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 + + 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) } @@ -950,202 +1108,225 @@ ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) { #' @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) +#' +#' 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(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 <- 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() + # browser() FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE) - if (is.character(FITs)) return(FITs) # if singularity - + 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)) + 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) + 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] - - - + 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]*100 Lim[[9]] & POTr_CI[2] * 100 < Lim[[10]]) testPOTr <- 0 else testPOTr <- 1 + potAllU2 <- FITs[[4]] test_c <- logical() - if((potAllU2[2] >Lim[[9]]/100 & potAllU2[3] Lim[[9]] / 100 & potAllU2[3] < Lim[[10]] / 100)) test_c <- 0 else test_c <- 1 predPot <- FITs[[5]] predPotU <- FITs[[6]] - - vcovMU <- V_V*smu$sigma^2 - + + vcovMU <- V_V * smu$sigma^2 + #### ANOVA ---- noConc <- length(unique(all_l$Conc)) nofitted <- noConc - AnovaDFs <- c(nofitted-1,1,3,nofitted-4-1,nrow(all_l)-nofitted, nofitted,nrow(all_l)-2*nofitted,nrow(all_l)-1) - SStreat <- round(sum((predPotU-mean(all_l$readout))^2),5) - SSregr <- round(sum((predPot-mean(all_l$readout))^2),5) + AnovaDFs <- c(nofitted - 1, 1, 3, nofitted - 4 - 1, nrow(all_l) - nofitted, nofitted, nrow(all_l) - 2 * nofitted, nrow(all_l) - 1) + SStreat <- round(sum((predPotU - mean(all_l$readout))^2), 5) + SSregr <- round(sum((predPot - mean(all_l$readout))^2), 5) # non-parallelism - SSnonparall <- round(sum(smr$residuals^2)-sum(smu$residuals^2),5) - SSprep <- round(sum((predict(lm(readout ~ isSample, all_l))-mean(all_l$readout))^2),5) - - RSS <- round(sum(smu$residuals^2),5) + SSnonparall <- round(sum(smr$residuals^2) - sum(smu$residuals^2), 5) + SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2), 5) + + RSS <- round(sum(smu$residuals^2), 5) RSS_df <- AnovaDFs[5] - MSEunr <- RSS/RSS_df - RMSEunr <- sqrt(RSS/RSS_df) + MSEunr <- RSS / RSS_df + RMSEunr <- sqrt(RSS / RSS_df) # Pure Err - FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l)) - SSE <- sum(resid(lm(readout ~ factor(Conc)*isSample, all_l))^2) # =FitAnova[4,2] - SSE_df <- FitAnova[4,1] - PureMSE <- SSE/SSE_df + FitAnova <- anova(lm(readout ~ factor(Conc) * isSample, all_l)) + SSE <- sum(resid(lm(readout ~ factor(Conc) * isSample, all_l))^2) # =FitAnova[4,2] + SSE_df <- FitAnova[4, 1] + PureMSE <- SSE / SSE_df RMSE_pure <- sqrt(PureMSE) ## non-lin = LoF - if (PureErrFlag) { ERR <- PureMSE; ERR_df <- SSE_df } else { ERR <- MSEunr; ERR_df <- RSS_df } - SSnonlin <- sum((predict(lm(readout ~ factor(Conc)*isSample, all_l))-predPotU)^2) - LoF_df <- FitAnova[1,1]+FitAnova[2,1] - F_regr <- (SSregr/AnovaDFs[3])/ERR - p_F_regr <- round(pf(F_regr, AnovaDFs[3], ERR_df, lower.tail = F),5) - if (ncol(ro_new)<4) F_nonlin <- 0 else F_nonlin <- (SSnonlin/AnovaDFs[6])/ERR + if (PureErrFlag) { + ERR <- PureMSE + ERR_df <- SSE_df + } else { + ERR <- MSEunr + ERR_df <- RSS_df + } + SSnonlin <- sum((predict(lm(readout ~ factor(Conc) * isSample, all_l)) - predPotU)^2) + LoF_df <- FitAnova[1, 1] + FitAnova[2, 1] + F_regr <- (SSregr / AnovaDFs[3]) / ERR + p_F_regr <- round(pf(F_regr, AnovaDFs[3], ERR_df, lower.tail = F), 5) + if (ncol(ro_new) < 4) F_nonlin <- 0 else F_nonlin <- (SSnonlin / AnovaDFs[6]) / ERR if (F_nonlin > 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" } - + 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) + + 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() + 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) + 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) + # 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) - + 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) + + 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) - + 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) + + 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) - + 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"]) + 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) + + 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)) + 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) } @@ -1159,26 +1340,26 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) { #' @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)) -#' -#' +#' +#' 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_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) + 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]] @@ -1188,48 +1369,57 @@ ANOVA4plUnresfunc <- function(ro_new) { # 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) + 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) + 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) + 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)) + 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] + 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] + 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 + 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,"","") + 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) } @@ -1242,38 +1432,36 @@ ANOVA4plUnresfunc <- function(ro_new) { #' @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) - - - +#' +#' 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))] + 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) - + 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) + 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. @@ -1288,18 +1476,20 @@ perConcTab <- function(ro_new, noDilSeries) { #' @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) +#' +#' 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) } -