1306 lines
57 KiB
R
1306 lines
57 KiB
R
################################################################################
|
|
# All functions for plateflow
|
|
################################################################################
|
|
|
|
|
|
#' Levenberg Marquard fit of 4 pl
|
|
#'
|
|
#' Returns list of following: summary of restricted and unrestricted model fits of the 4pl function,
|
|
#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof.
|
|
#'
|
|
#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration.
|
|
#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE).
|
|
#' @returns Returns list of following: summary of restricted and unrestricted model fits of the 4pl function,
|
|
#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof.
|
|
#' @export
|
|
#' @examples
|
|
#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
|
|
#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
|
|
#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
|
|
#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764))
|
|
#' TransF<- FALSE
|
|
#' Dat<- list()
|
|
#' te <- Fitting_FUNC(dat,TransF)
|
|
#' print(te)
|
|
|
|
|
|
Fitting_FUNC <- function(ro_new, TransFlag=F) {
|
|
CORro <- cor(ro_new[,1],ro_new[,ncol(ro_new)])
|
|
#browser()
|
|
all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
|
|
isRef <- rep(c(1,0),1,each=nrow(all_l)/2)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l2 <- cbind(all_l, isRef, isSample)
|
|
#browser()
|
|
if (CORro<0) SLOPE <- -1 else SLOPE <- 1
|
|
if (!TransFlag) {
|
|
startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0)
|
|
|
|
mr <- tryCatch({gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))),
|
|
data=all_l2,
|
|
start=startlist,#race=T,
|
|
control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
|
|
},
|
|
warning = function(e) {
|
|
mr <<- "In nlsModel singular gradient matrix"
|
|
})
|
|
# Stop if singular gradient matrix
|
|
if (is.character(mr)) return(mr)
|
|
|
|
s_mr <- tryCatch({
|
|
s_mr <- summary(mr)
|
|
},
|
|
error = function(err) {
|
|
s_mr <- NULL
|
|
})
|
|
} else {
|
|
startlist <- list(a=log(min(ro_new[,2])), b=SLOPE, d=log(max(ro_new[,2])), cs=mean(all_l$log_dose),r=0)
|
|
mrT <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))),
|
|
data=all_l2,
|
|
start=startlist,
|
|
control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
|
|
s_mr <- summary(mrT)
|
|
}
|
|
|
|
if (!TransFlag) {
|
|
startlistmu <- list(as=min(ro_new[,2]), bs=SLOPE, ds=max(ro_new[,2]), cs=mean(all_l$log_dose),
|
|
at=min(ro_new[,2]), bt=SLOPE, dt=max(ro_new[,2]), r=0)
|
|
tryCatch({
|
|
mu <- gsl_nls(fn = readout ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
|
|
(1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))),
|
|
data=all_l,
|
|
start=startlistmu,
|
|
control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
|
|
},
|
|
error = function(msg){
|
|
return(0) })
|
|
|
|
Sum_u <- tryCatch({ summary(mu) },
|
|
error=function(msg){
|
|
return(0) })
|
|
} else {
|
|
startlistmu <- list(as=log(min(ro_new[,2])), bs=SLOPE, ds=log(max(ro_new[,2])), cs=mean(all_l$log_dose),
|
|
at=log(min(ro_new[,2])), bt=SLOPE, dt=log(max(ro_new[,2])), r=0)
|
|
tryCatch({
|
|
muT <- gsl_nls(fn = log(readout) ~ as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
|
|
(1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(cs-r*isSample-log_dose))),
|
|
data=all_l,
|
|
start=startlistmu,
|
|
control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
|
|
},
|
|
error = function(msg){
|
|
return(0) })
|
|
|
|
Sum_u <- tryCatch({ summary(muT) },
|
|
error=function(msg){
|
|
return(0) })
|
|
}
|
|
if (!TransFlag) {
|
|
pot_est <- exp(confintd(mr, "r", method="asymptotic"))
|
|
potU_est <- exp(confintd(mu, "r", method="asymptotic"))
|
|
PRED <- predict(mr)
|
|
PREDu <- predict(mu)
|
|
} else {
|
|
pot_est <- exp(confintd(mrT, "r", method="asymptotic"))
|
|
potU_est <- exp(confintd(muT, "r", method="asymptotic"))
|
|
PRED <- predict(mrT)
|
|
PREDu <- predict(muT)
|
|
}
|
|
return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu))
|
|
}
|
|
|
|
|
|
#' Plot when the fitting has a singularity
|
|
#'
|
|
#' Returns the scatter plot of the data, with REF in blue and TEST in red.
|
|
#' One plots are returned.
|
|
#'
|
|
#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration.
|
|
#' @returns A grid object with 2 linearity plots, restricted and unrestricted model.
|
|
#' @export
|
|
#' @examples
|
|
#' data.frame(R_dil1 = c(10.0651024695491, 10.9844983291817, 10.7635586089293, 10.4597656321327, 10.3898668457823, 10.8171761349909,
|
|
#' 10.319758021908, 10.1304854046653),
|
|
#' R_dil2 = c(10.9649145494504, 10.0202868589385, 10.8424145955735, 10.9311360356894, 10.3284659026404,
|
|
#' 10.6890147558796, 10.3014450252305, 10.9594838595181),
|
|
#' R_dil3 = c(10.4630510824383, 10.4566715089363, 10.2350765290036, 10.3300581874798, 10.9648088137065,
|
|
#' 10.286893755805, 10.4856643841389, 10.5275521552307),
|
|
#' T_dil1 = c(12.732175566336, 12.7756403995095, 12.1672539684741, 12.7060603907892, 12.8000685682832,
|
|
#' 12.8800092157515, 12.7160581291873, 12.6996878912416),
|
|
#' T_dil2 = c(12.3923194313831, 12.0943488144175, 12.7955302154828, 12.4825917078735, 12.6856540203788,
|
|
#' 12.7348548498556, 12.9222470610476, 12.1186618671252),
|
|
#' T_dil3 = c(12.7899182255274, 12.9722600411128, 12.7078445380891, 12.4913523531941, 12.1718281909609,
|
|
#' 12.5313873615133, 12.952802332772, 12.5960321394342),
|
|
#' log_dose = c(0, -1.09861228866811, -2.19722457733622, -3.29583686600433, -4.39444915467244,
|
|
#' -5.49306144334055, -6.59167373200866, -7.69028602067677))
|
|
#'
|
|
#'
|
|
#' p <- plotSingularity(dat)
|
|
#' print(p)
|
|
#'
|
|
plotSingularity <- function(dat) { #sigmoid,det_sig,
|
|
CORdat <- cor(dat[,1],dat[,ncol(dat)])
|
|
#browser()
|
|
all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout")
|
|
isRef <- rep(c(1,0),1,each=nrow(all_l)/2)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l2 <- cbind(all_l, isRef, isSample)
|
|
#browser()
|
|
|
|
log_dose <- unique(all_l$log_dose)
|
|
seq_x <- seq(min(log_dose),max(log_dose),0.1)
|
|
|
|
|
|
#browser()
|
|
#all_l2$readout[all_l2$readout < 0] <- 0.01
|
|
all_l2$readouttrans <- log(all_l2$readout)
|
|
|
|
#browser()
|
|
pSing <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
|
|
geom_point(shape=factor(isRef), size=3,alpha=0.8) +
|
|
labs(title = paste("No 4pl fit possible"),
|
|
color="product") +
|
|
scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) +
|
|
scale_shape_manual(labels=c("test","reference")) +
|
|
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
|
|
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
|
|
#theme_bw() +
|
|
theme(axis.text = element_text(size=14))
|
|
|
|
return(pSing)
|
|
}
|
|
|
|
|
|
|
|
|
|
#' Plot sigmoidal curve
|
|
#'
|
|
#' Returns the final plots of the 4pl function as sigmoidal lines, and the single readouts as scatter, with REF in blue and TEST in red.
|
|
#' Two plots are returned in a grid object: the unrestricted model and the restricted model.
|
|
#'
|
|
#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration.
|
|
#' @param sigmoid The parameters of the 4pl fit of unrestricted model from EXCEL input.
|
|
#' @param det_sig The parameters of the 4pl fit of unrestricted model from the metadata input.
|
|
#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE).
|
|
#' @returns A grid object either of the original scale or the natural log of the readouts.
|
|
#' @export
|
|
#' @examples
|
|
#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
|
|
#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
|
|
#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
|
|
#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764))
|
|
#' sigmoid <- c(0.7163324, 0.5636804, 10.6156340 , 9.9784160, -0.7504673, -0.7108692, -3.5788141, -0.6662962)
|
|
#' det_sig <- FALSE
|
|
#' TransF<- FALSE
|
|
#' Dat<- list()
|
|
#' p <- plot_f(dat,sigmoid,det_sig,TransF)
|
|
#' print(p)
|
|
|
|
|
|
plot_f <- function(dat, TransFlag=F) { #sigmoid,det_sig,
|
|
CORdat <- cor(dat[,1],dat[,ncol(dat)])
|
|
#browser()
|
|
all_l <- melt(data.frame(dat), id.vars="log_dose", variable.name="replname", value.name = "readout")
|
|
isRef <- rep(c(1,0),1,each=nrow(all_l)/2)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l2 <- cbind(all_l, isRef, isSample)
|
|
#browser()
|
|
MODLS <- Fitting_FUNC(dat, TransFlag=F)
|
|
s_mr <- MODLS[[1]]
|
|
a <- s_mr$coefficients["a",1]
|
|
b <- s_mr$coefficients["b",1]
|
|
cs <- s_mr$coefficients["cs",1]
|
|
d <- s_mr$coefficients["d",1]
|
|
r <- s_mr$coefficients["r",1]
|
|
|
|
log_dose <- unique(all_l$log_dose)
|
|
seq_x <- seq(min(log_dose),max(log_dose),0.1)
|
|
SAMPLE <- a+(d-a)/(1+exp(b*((cs-r)-seq_x)))
|
|
REF <- a+(d-a)/(1+exp(b*((cs)-seq_x)))
|
|
|
|
s_mu <- MODLS[[2]]
|
|
|
|
as <- s_mu$coefficients["as",1]
|
|
bs <- s_mu$coefficients["bs",1]
|
|
cs <- s_mu$coefficients["cs",1]
|
|
ds <- s_mu$coefficients["ds",1]
|
|
at <- s_mu$coefficients["at",1]
|
|
bt <- s_mu$coefficients["bt",1]
|
|
ct <- s_mu$coefficients["cs",1]-s_mu$coefficients["r",1]
|
|
dt <- s_mu$coefficients["dt",1]
|
|
r <- s_mu$coefficients["r",1]
|
|
|
|
|
|
SAMPLEtrue <- at + (dt -at)/(1+exp(bt*((cs-r-seq_x))))
|
|
REFtrue <- as + (ds -as)/(1+exp(bs*((cs-seq_x))))
|
|
|
|
#browser()
|
|
pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue)
|
|
all_l2$readout[all_l2$readout < 0] <- 0.01
|
|
all_l2$readouttrans <- log(all_l2$readout)
|
|
slopeEC50 <- b*(d-a)/4
|
|
Intercept <- a+(d-a)/2-b*(d-a)/4*cs
|
|
#browser()
|
|
Xbendl3 <- cs-(1.5434/b)
|
|
Xbendu3 <- cs+(1.5434/b)
|
|
XbendlT <- cs-r-(1.5434/b)
|
|
XbenduT <- cs-r+(1.5434/b)
|
|
XasymplS <- cs-(3/b)
|
|
XasympuS <- cs+(3/b)
|
|
XasymplT <- cs-r-(3/b)
|
|
XasympuT <- cs-r+(3/b)
|
|
bendpoints <- c(bendREF_lower = round(Xbendl3,3), bendREF_upper=round(Xbendu3,3),
|
|
bendSAMPLE_lower = round(XbendlT,3), bendSAMPLE_upper=round(XbenduT,3),
|
|
asympREF_lower = round(XasymplS,3), asympREF_upper=round(XasympuS,3),
|
|
asympSAMPLE_lower = round(XasymplT,3), asympSAMPLE_upper=round(XasympuT,3))
|
|
Dat$bendpoints <- bendpoints
|
|
Dat$cfordils <- cs
|
|
#browser()
|
|
p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
|
|
geom_point(shape=factor(isRef), size=3,alpha=0.8) +
|
|
labs(title = paste("restricted 4pl model"),
|
|
color="product") +
|
|
scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) +
|
|
scale_shape_manual(labels=c("test","reference")) +
|
|
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
|
|
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
|
|
#theme_bw() +
|
|
theme(axis.text = element_text(size=14))
|
|
|
|
p2 <- p + geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLE), color="#C2173F",
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REF), color="#4545BA",
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=SAMPLEtrue), color="#C2173F", linetype=2, alpha=0.4,
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df), aes(x=seq_x, y=REFtrue), color="#4545BA", linetype=2, alpha=0.4,
|
|
inherit.aes = F) +
|
|
geom_vline(xintercept=c(Xbendl3, Xbendu3), col="#4545BA",linetype=2) +
|
|
geom_vline(xintercept=c(XbendlT, XbenduT), col="#C2173F",linetype=2) +
|
|
geom_vline(xintercept=c( XasymplS, XasympuS), col="#4545BABB",linetype=3) +
|
|
geom_vline(xintercept=c(XasymplT, XasympuT), col="#C2173FBB",linetype=3) +
|
|
annotate("text", x=cs, y=a+(d-a)/2, label="0", size=5) +
|
|
geom_abline(slope = slopeEC50, intercept = Intercept) +
|
|
theme(legend.position="none")
|
|
Dat$p2 <- p2
|
|
|
|
# transformed plots
|
|
p_rt <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) +
|
|
geom_point(shape=factor(isRef), size=3,alpha=0.8) +
|
|
labs(title = paste("restricted transformed 4pl"), color="product") +
|
|
scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) +
|
|
theme_bw()
|
|
|
|
MODLStrans <- Fitting_FUNC(dat,TransFlag = TRUE)
|
|
s_mrt <- MODLStrans[[1]]
|
|
a_trans <- s_mrt$coefficients["a",1]
|
|
b_trans <- s_mrt$coefficients["b",1]
|
|
cs_trans <- s_mrt$coefficients["cs",1]
|
|
d_trans <- s_mrt$coefficients["d",1]
|
|
r_trans <- s_mrt$coefficients["r",1]
|
|
|
|
XbendlTrans <- cs_trans-(1.5434/b_trans)
|
|
XbenduTrans <- cs_trans+(1.5434/b_trans)
|
|
XbendlTransT <- cs_trans-r_trans-(1.5434/b_trans)
|
|
XbenduTransT <- cs_trans-r_trans+(1.5434/b_trans)
|
|
bendpointsTRANS <- c(bendREF_lower = round(XbendlTrans,3), bendREF_upper=round(XbenduTrans,3),
|
|
bendSAMPLE_lower = round(XbendlTransT,3), bendSAMPLE_upper=round(XbenduTransT,3))
|
|
Dat$bendpointsTRANS <- bendpointsTRANS
|
|
SAMPLEtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans-r_trans)-seq_x)))
|
|
REFtrans <- a_trans+(d_trans-a_trans)/(1+exp(b_trans*((cs_trans)-seq_x)))
|
|
|
|
pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans)
|
|
p_rt2 <- p_rt + geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=SAMPLEtrans), color="#C2173F",
|
|
inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df_trans), aes(x=seq_x, y=REFtrans), color="#4545BA",
|
|
inherit.aes = F) +
|
|
geom_vline(xintercept=c(XbendlTrans, XbenduTrans), col="#4545BA",linetype=2) +
|
|
geom_vline(xintercept=c(XbendlTransT, XbenduTransT), col="#C2173F",linetype=2) +
|
|
theme(legend.position = "none", axis.text=element_text(size=14))
|
|
|
|
#browser()
|
|
Sum_u <- MODLS[[2]]
|
|
#browser()
|
|
#if (is.null(det_sig)) {
|
|
ast <- Sum_u$coefficients["as",1]
|
|
ate <- Sum_u$coefficients["at",1]
|
|
bst <- Sum_u$coefficients["bs",1]
|
|
bte <- Sum_u$coefficients["bt",1]
|
|
cst <- Sum_u$coefficients["cs",1]
|
|
cte <- Sum_u$coefficients["cs",1]-Sum_u$coefficients["r",1]
|
|
dst <- Sum_u$coefficients["ds",1]
|
|
dte <- Sum_u$coefficients["dt",1]
|
|
# } else {
|
|
# ast <- det_sig[5]
|
|
# ate <- det_sig[6]
|
|
# bst <- det_sig[1]
|
|
# bte <- det_sig[2]
|
|
# cst <- det_sig[7]
|
|
# cte <- det_sig[8]
|
|
# dst <- det_sig[3]
|
|
# dte <- det_sig[4]
|
|
# }
|
|
REFu <- ast + (dst-ast)/(1+exp(bst*(cst-seq_x)))
|
|
SAMPLEu <- ate + (dte-ate)/(1+exp(bte*(cte-seq_x)))
|
|
pl_df2 <- cbind(seq_x, SAMPLEu, REFu)
|
|
#browser()
|
|
pu <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
|
|
geom_point(size=3) +
|
|
labs(title="unrestricted 4pl model", color="product") +
|
|
scale_color_manual(labels = c("test","reference"), values=c("#C2173F88","#4545BA88")) +
|
|
theme_bw()
|
|
pu2 <- pu + geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=SAMPLEu),
|
|
color="#C2173F", inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df2), aes(x=seq_x, y=REFu),
|
|
color="#4545BA", inherit.aes = F,
|
|
show.legend = F)
|
|
pu2_ <- pu2 +
|
|
theme(legend.position = "none", axis.text = element_text(size=14))
|
|
|
|
putrans <- ggplot(all_l2, aes(x=log_dose, y=readouttrans, color=factor(isRef))) +
|
|
geom_point( size=3) +
|
|
labs(title="unrestricted transformed 4_pl-Model", color="product") +
|
|
scale_color_manual(labels = c("test","reference"), values=c("#C2173FAA","#4545BAAA")) +
|
|
theme_bw()
|
|
|
|
Sum_ut <- MODLStrans[[2]]
|
|
|
|
ast_t <- Sum_ut$coefficients["as",1]
|
|
ate_t <- Sum_ut$coefficients["at",1]
|
|
bst_t <- Sum_ut$coefficients["bs",1]
|
|
bte_t <- Sum_ut$coefficients["bt",1]
|
|
cst_t <- Sum_ut$coefficients["cs",1]
|
|
cte_t <- Sum_ut$coefficients["cs",1]-Sum_ut$coefficients["r",1]
|
|
dst_t <- Sum_ut$coefficients["ds",1]
|
|
dte_t <- Sum_ut$coefficients["dt",1]
|
|
|
|
REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(cst_t-seq_x)))
|
|
SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(cte_t-seq_x)))
|
|
pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans)
|
|
|
|
pu2_t <- putrans + geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=SAMPLEu_trans),
|
|
color="#C2173F", inherit.aes = F) +
|
|
geom_line(data=as.data.frame(pl_df2u_t), aes(x=seq_x, y=REFu_trans),
|
|
color="#4545BA", inherit.aes = F,
|
|
show.legend = F)
|
|
pu3_t <- pu2_t
|
|
if (TransFlag) grid.arrange(p_rt2,pu3_t, nrow=1) else grid.arrange(p2,pu2_, nrow=1)
|
|
}
|
|
|
|
#' Simulate a well-plate readout
|
|
#'
|
|
#' Returns a possible well-plate result based on the metadata provided and a variability.
|
|
#' Homo- and heteroscedasticity can be simulated.
|
|
#'
|
|
#' @param as, bs, cs, ds, at, bt, dt,r,ct, gt, gs are the parameter of the 5pl logistic regression.
|
|
#' @param sd_fac The standard deviation with units in % of upper-lower asymptote.
|
|
#' @param log_conc The natural log of the concentrations.
|
|
#' @param heteroNoise Boolean if heteroscedstic noise should be added.
|
|
#' @param noDilSeries A number >1 indicating how many dilution series per product should be simulated.
|
|
#' @param noDils A number >7 indicating how many dilutions steps per product should be simulated.
|
|
#' @returns A data-frame with readouts and natural log of concentrations.
|
|
#' @export
|
|
#' @examples
|
|
#' as=3; bs=1; cs=-4; ds=10; at=3; bt=1; dt=10;r=0.0001;ct=cs-r;sd_fac=0.1; gt=1; gs=1;
|
|
#' lnConc=c(1,0,-1,-2,-3,-4,-5,-6)
|
|
#' heteroNoise <- FALSE
|
|
#' noDilS <- 3
|
|
#' noD <- 8
|
|
#' Calc_DilRes(as,bs,cs,ds,at,bt,dt,r,ct,sd_fac,gt,gs,log_conc=lnConc, heteroNoise, noDilS, noD)
|
|
|
|
Calc_DilRes <- function(as=3, bs=1, cs=-4, ds=10, at=3, bt=1, dt=10,r=0.0001,ct=cs-r,
|
|
sd_fac=0.1, gt=1, gs=1, log_conc,
|
|
heteroNoise=FALSE, noDilSeries, noDils) {
|
|
yAxfac <- (ds-as)
|
|
log_dose <- log_conc
|
|
isRef <- rep(c(1,0),1,each=length(log_conc)*noDilSeries)
|
|
isSample <- rep(c(0,1),1,each=length(log_conc)*noDilSeries)
|
|
#browser()
|
|
av <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
|
|
(1+isRef*exp(bs*(cs - log_dose)) + isSample*exp(bt*(ct-log_dose)))
|
|
if (heteroNoise) {
|
|
# heterosc noise
|
|
ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac/100))), nrow=noDils, ncol=noDilSeries*2)
|
|
} else {
|
|
# homosc noise
|
|
ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac*yAxfac/100))), nrow=noDils, ncol=noDilSeries*2)
|
|
}
|
|
ro_jit <- abs(ro_jit)
|
|
|
|
ro_jit2 <- cbind(ro_jit, log_dose)
|
|
if (noDilSeries==3) {
|
|
colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose")
|
|
} else {
|
|
colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose")
|
|
}
|
|
return(ro_jit2)
|
|
}
|
|
|
|
#' Calculate the potency of linear PLA
|
|
#'
|
|
#' The delta Method is used for calculating the potency using a restricted linear model.
|
|
#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used.
|
|
#'
|
|
#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts.
|
|
#' @param Lim The limits to check if the relative CI is within the limits.
|
|
#' @param PureErrFlag Boolean if Pure Error should be used.
|
|
#' @returns A data-frame with potency estimate, absolute CIs, test result, relative CIs.
|
|
#' @export
|
|
#' @examples
|
|
#' CIRC <- data.frame(log_dose = c(-2.5,-2.5,-2.5, -3.2,-3.2,-3.2,-3.9,-3.9,-3.9,
|
|
#' -3.2,-3.2,-3.2,-3.9,-3.9,-3.9,-4.7,-4.7,-4.7),
|
|
#' replname= c("R_dil1","R_dil1","R_dil1", "R_dil2","R_dil2","R_dil2", "R_dil3","R_dil3","R_dil3",
|
|
#' "T_dil1","T_dil1","T_dil1", "T_dil2","T_dil2","T_dil2", "T_dil3","T_dil3","T_dil3"),
|
|
#' readout = c(72.1,75.8,76.04,59.8,61,62.7,43.6,45,41.5,53.5,62.2,65.9,48.3,43.8,43.14,28.17,29.2,31.2),
|
|
#' isRef=c(rep(1,9), rep(0,9)),
|
|
#' isSample = c(rep(0,9), rep(1,9)))
|
|
#' Lim <- c(rep(0,8),70,130) # only Lim 9 and 10 relevant
|
|
#' PureErrF <- TRUE
|
|
#'
|
|
#'
|
|
#' LinPotTab(circles=CIRC, Lim, PureErrF)
|
|
|
|
|
|
LinPotTab <- function(circles, Lim, PureErrFlag) {
|
|
circ_ABl <- circles
|
|
circ_Al <- circ_ABl[circ_ABl$isSample ==1,]
|
|
circ_Al <- circ_ABl[circ_ABl$isSample ==0,]
|
|
# restr CSSI model
|
|
modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
|
|
coeffs <- modAB$coefficients
|
|
SU_modAB <- tryCatch({
|
|
SU_modAB <- summary(modAB)
|
|
}, error = function(msg) {
|
|
return(NA)
|
|
})
|
|
# Intercept diff/slope modAB
|
|
linPot <- exp(modAB$coefficients[3]/modAB$coefficients[2])
|
|
|
|
if(PureErrFlag) {
|
|
FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
|
|
meanPureErr <- FitAnova[4,3]
|
|
DFsPure <- FitAnova[4,1]
|
|
VCOV <- vcov(modAB)
|
|
V_V <- VCOV/SU_modAB$sigma^2
|
|
VCOVpure <- V_V*meanPureErr
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
}
|
|
|
|
log_pot_delta <- deltaMethod(modAB, "isSample/log_dose")
|
|
if (PureErrFlag) {
|
|
V_ <- log_pot_delta$SE^2/SU_modAB$sigma^2
|
|
V_p <- V_*meanPureErr
|
|
potDeltaPureSE <- sqrt(V_p)
|
|
CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure)*potDeltaPureSE
|
|
CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure)*potDeltaPureSE
|
|
} else {
|
|
CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB))*log_pot_delta$SE
|
|
CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB))*log_pot_delta$SE
|
|
}
|
|
#browser()
|
|
ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))
|
|
|
|
# Rel pot CI
|
|
relLinpotCI <- ExpLinPot/linPot*100
|
|
if (relLinpotCI[2]>Lim[[9]] & relLinpotCI[3]<Lim[[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)
|
|
}
|
|
|
|
|
|
#' 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.
|
|
#'
|
|
#' @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
|
|
#' 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
|
|
#' 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$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()
|
|
circ_ABl <- circles
|
|
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)
|
|
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
|
|
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
|
|
|
|
su_mod <- summary(modAB)$coefficients
|
|
su_mod2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.")), su_mod)
|
|
su_modU <- summary(modABu)$coefficients
|
|
su_modU2 <- cbind(data.frame(parameter = c("intercept REF","slope REF","intercepts diff.","slope difference")), su_modU)
|
|
|
|
uCI_SloDiff <- su_modU[4,1] + qt(0.975,8)*su_modU[4,2]
|
|
lCI_SloDiff <- su_modU[4,1] - qt(0.975,8)*su_modU[4,2]
|
|
SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiff,uCI_SloDiff)
|
|
|
|
lenCirc <- nrow(circ_ABl)
|
|
dfTreat <- 3
|
|
dfPrep <- 1
|
|
dfReg <- 1
|
|
dfnonP <- 1
|
|
dfRMSE <- c(lenCirc-3-1)
|
|
dfTotal <- lenCirc-1
|
|
dfPureE <- lenCirc-6
|
|
dfNonLin <- dfRMSE-dfPureE
|
|
|
|
RSS <- sum(resid(lm(readout ~ log_dose*isSample, circ_ABl))^2)
|
|
MSE <- RSS/dfRMSE
|
|
SSE <- sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)
|
|
MSpure <- SSE/dfPureE
|
|
|
|
if (PureErrFlag) {
|
|
FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, circ_ABl))
|
|
meanPureErr <- FitAnova[4,3]
|
|
SU_modAB <- tryCatch({
|
|
SU_modAB <- summary(modAB)
|
|
}, error = function(msg) {
|
|
return(NA)
|
|
})
|
|
if (length(SU_modAB)>1) s_modABcoeffs <- summary(modAB)$coefficients
|
|
|
|
DFsPure <- FitAnova[4,1]
|
|
VCOV <- vcov(modAB)
|
|
V_V <- VCOV/SU_modAB$sigma^2
|
|
VCOVpure <- V_V*meanPureErr
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
su_mod2[,3] <- SEsPure
|
|
su_mod2[,4] <- su_mod2[,2]/su_mod2[,3]
|
|
su_mod2[,5] <- 2*(1-pt(abs(su_mod[,4]), FitAnova[4,1]))
|
|
|
|
s_mu <- summary(modABu)$coefficients
|
|
SU_modABu <- summary(modABu)
|
|
VCOVu <- vcov(modABu)
|
|
V_Vu <- VCOVu/SU_modABu$sigma^2
|
|
SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
|
|
|
|
su_modU2[,3] <- SEsPureU
|
|
su_modU2[,4] <- su_modU2[,2]/su_modU2[,3]
|
|
su_modU2[,5] <- 2*(1-pt(abs(su_modU2[,4]), FitAnova[4,1]))
|
|
|
|
uCI_SloDiffP <- su_modU[4,1] + qt(0.975,8)*SEsPureU[4]
|
|
lCI_SloDiffP <- su_modU[4,1] - qt(0.975,8)*SEsPureU[4]
|
|
SlopeDiffCI <- c(su_modU[4,1], lCI_SloDiffP,uCI_SloDiffP)
|
|
|
|
SSRes <- SSE
|
|
dfRes <- dfPureE
|
|
|
|
} else {
|
|
SSRes <- RSS
|
|
dfRes <- dfRMSE
|
|
}
|
|
|
|
# treatment
|
|
SStreat <- print(sum((predict(lm(readout ~ factor(log_dose)*isSample, circ_ABl))-mean(circ_ABl$readout))^2))
|
|
F_treat <- (SStreat/dfTreat)/(SSRes/dfRes)
|
|
# Preparation
|
|
SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl))-mean(circ_ABl$readout))^2))
|
|
F_prep <- (SSprep/dfTreat)/(SSRes/dfRes)
|
|
# Regression
|
|
# ANOVA tape II SS of regression
|
|
SSreg <- Anova(lm(readout ~log_dose + isSample, circ_ABl))[1,1]
|
|
# Non-parallelism
|
|
# diff of RSS of restricted and unrestricted model
|
|
SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2)
|
|
F_nonpar <- SSnonpar/(sum(resid(lm(readout ~ factor(log_dose)*isSample, circ_ABl))^2)/(lenCirc-4))
|
|
|
|
# non-linearity
|
|
SSnonlin <- sum((predict(modABu)-predict(lm(readout ~ as.factor(log_dose)*isSample, circ_ABl)))^2)
|
|
# = RSS-SSE
|
|
# Total SS
|
|
SStot <- sum((circ_ABl$readout-mean(circ_ABl$readout))^2)
|
|
# Significance of R^2 F-ratio
|
|
# MSR/MSE
|
|
# sample A
|
|
F_R2_A <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2)/
|
|
(sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Al)) - circ_Al$readout)^2)/(nrow(circ_Al)-3))
|
|
pFR2_A <- round(pf(F_R2_A,1,6),4)
|
|
# sample B
|
|
F_R2_B <- sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2)/
|
|
(sum((predict(lm(readout ~ log_dose+ I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2)/(nrow(circ_Bl)-3))
|
|
pFR2_B <- round(pf(F_R2_B,1,6),4)
|
|
# sign of non-lin with pure error: MSSnonlin/MSSE
|
|
F_nonlin <- (SSnonlin/2)/(SSE/dfPureE)
|
|
|
|
# sign of slope
|
|
F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2)/(sum((circ_Bl$readout - predict(modB))^2)/(nrow(circ_Bl)-2))
|
|
F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2)/(sum((circ_Al$readout - predict(modA))^2)/(nrow(circ_Al)-2))
|
|
# F-test on regression: MSSreg/MSSE
|
|
if (is.na(F_nonlin)) F_nonlin <- 0
|
|
if (F_nonlin > 0) {
|
|
p_F_nonlin <- round(pf(F_nonlin,2,dfPureE, lower.tail = F),5)
|
|
} else { p_F_nonlin <- "SSnonlin neg or 0"; }
|
|
|
|
# significances
|
|
F_regr <- (SSreg/1)/(SSRes/dfRes)
|
|
p_F_regr <- round(pf(F_regr,1,dfRes, lower.tail = F),3)
|
|
p_F_treat <- round(pf(F_treat,3,dfRes, lower.tail = F),3)
|
|
p_F_prep <- round(pf(F_prep,1,dfRes, lower.tail = F),3)
|
|
p_F_slope_A <- round(pf(F_slope_A,1,(nrow(circ_Al)-2), lower.tail = F),3)
|
|
p_F_slope_B <- round(pf(F_slope_B,1,(nrow(circ_Bl)-2), lower.tail = F),3)
|
|
p_F_nonp <- round(pf(F_nonpar,1,dfRes, lower.tail = F),3)
|
|
p_F_LoF <- p_F_nonlin
|
|
|
|
res_tab_lin <- data.frame(test = c("F-test on sign. of regression", "F_test on non-lin",
|
|
"F-test on R^2 A","F_test on R^2 B",
|
|
"F-test on slope A","F-test on slope B",
|
|
"F-test on non-parallelism","F-test on preparation"),
|
|
test_results = c(ifelse(p_F_regr<0.05,0,1),ifelse(p_F_nonlin<0.05,1,0),
|
|
ifelse(pFR2_A<0.05,1,0),ifelse(pFR2_B<0.05,1,0),
|
|
ifelse(p_F_slope_A<0.05,0,1),ifelse(p_F_slope_B<0.05,0,1),
|
|
ifelse(p_F_nonp<0.05,1,0),ifelse(p_F_prep<0.05,0,1)),
|
|
estimate = c(p_F_regr, p_F_nonlin,pFR2_A,pFR2_B,p_F_slope_A,
|
|
p_F_slope_B,p_F_nonp,p_F_prep),
|
|
Source = c("Treatment","Preparation","Regression","Non-parallelism",
|
|
"Resid Error","Non-linearity","Pure error", "Total"),
|
|
df = c(dfTreat,1,1,1,dfRMSE,2,dfPureE,lenCirc-1),
|
|
SumSquares = c(round(SStreat,5),round(SSprep,5),round(SSreg,5),
|
|
round(SSnonpar,5),round(RSS,5),round(SSnonlin,5),
|
|
round(SSE,5),round(SStot,5)),
|
|
MS = c(round(SStreat/dfTreat,5),round(SSprep,5),round(SSreg,5),
|
|
round(SSnonpar,5),round(RSS/dfRMSE,5),round(SSnonlin/2,5),
|
|
round(SSE/dfPureE,5),round(SStot/dfTotal,5)),
|
|
"F-value" = c(round(F_treat,5), round(F_prep,5),round(F_regr,5),
|
|
round(F_nonpar,5),"",round(F_nonlin,5),"",""),
|
|
"p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "",""))
|
|
RET <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2)
|
|
return(RET)
|
|
}
|
|
|
|
|
|
#' Plots linear regression on the scatterplot.
|
|
#'
|
|
#' Restricted and unrestricted model are plotted. The number of dilution steps used for the regression
|
|
#' are printed in the title. Reference is blue, red is the sample. The points used for regression
|
|
#' are highlighted with a circle.
|
|
#'
|
|
#' @param circle Data frame with the 3 consecutive concentrations and the respective readouts.
|
|
#' @param sigmoid Parameters of the unrestricted fit of the full dataset for drawing the curves.
|
|
#' @param all_l2 Long format data.frame of the readout data incl log(concentration) and isRef (0s and 1s).
|
|
#' @param pl_df Data.frame for plotting the regression lines.
|
|
#' @param indS Index of dilution, where regression starts for reference sample.
|
|
#' @param indT Index of dilution, where regression starts for test sample.
|
|
#' @returns A grid object of 2 plots of unrestricted and restricted linear PLA models.
|
|
#' @export
|
|
#' @examples
|
|
#' circle <- read.csv("~/plateflow/circle.csv")
|
|
#' all_l2 <- read.csv("~/plateflow/all_l2.csv")
|
|
#' sigmoid <- c(10.0, 10.0, 110.0, 110.0, 1.0, 1.0, -3.5, 0.0)
|
|
#' indS <- 3
|
|
#' indT <- 3
|
|
#' pl_df <- data.frame(lnC=c(-1.203973,-1.897120 ,-2.590267,-3.283414,-3.976562,-4.669176,-5.362323,-6.053340),
|
|
#' plotS = c(113.772511,97.668371,81.564231,65.460091,49.355952,33.264200,17.160060,1.105405),
|
|
#' plotT = c(114.213375,97.588663,80.963951,64.339239,47.714527,31.102604,14.477892,-2.095735))
|
|
#'
|
|
#'
|
|
#' PlotLinPLA_FUNC(circle, sigmoid, all_l2, pl_df, indS,indT)
|
|
|
|
|
|
PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) {
|
|
#browser()
|
|
mLin <- gsl_nls(readout ~ (intS+r)*isSample + intS*isRef + k*log_dose,
|
|
data=circle,
|
|
start=list(intS = 0, k=1,r=0),
|
|
control = gsl_nls_control(xtol=1e-10,ftol=1e-10,gtol=1e-10))
|
|
# alternativ: modAB <- lm(readout ~ log_dose+isSample, circle)
|
|
sum_mLin <- summary(mLin)
|
|
|
|
log_dose <- unique(all_l2$log_dose)
|
|
seq_x <- seq(min(log_dose), max(log_dose),0.1)
|
|
if (!is.null(sigmoid)) {
|
|
SAMPLEtrue <- sigmoid[5] + (sigmoid[7]-sigmoid[5])/(1+exp(sigmoid[6]*((sigmoid[4]-sigmoid[8]-seq_x))))
|
|
REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[2]*((sigmoid[4]-seq_x))))
|
|
truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue)
|
|
} else {
|
|
SAMPLEtrue <- NULL
|
|
REFtrue <- NULL
|
|
truePL_df <- NULL
|
|
}
|
|
|
|
|
|
|
|
|
|
p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) +
|
|
geom_point(size=2) +
|
|
#labs(title=paste("linear regression model", indS,indT), color="product") +
|
|
scale_colour_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) +
|
|
ylim(min(all_l2$readout),max(all_l2$readout)) +
|
|
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
|
|
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
|
|
theme_bw()
|
|
p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="#4545BA",
|
|
inherit.aes = F) +
|
|
geom_line(data=pl_df,aes(x=lnC,y=plotT),color="#C2173F",
|
|
inherit.aes = F) +
|
|
{if (!is.null(truePL_df)) geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4,
|
|
inherit.aes = F) } +
|
|
{if (!is.null(truePL_df)) geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4,
|
|
inherit.aes = F)} +
|
|
labs(title = paste("unrestricted PLA model"), subtitle = paste("Regression starts for reference sample:",indS, "for test sample:",indT)) +
|
|
theme(legend.position="none", axis.text = element_text(size=14))
|
|
p3 <- p2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef),
|
|
size=5,alpha=0.2), col=c("black"), inherit.aes = FALSE) +
|
|
scale_shape_manual(labels=c("test","reference"), values=c(21,21))
|
|
# fit intercept for test and ref and common slope
|
|
|
|
|
|
pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log_dose
|
|
pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log_dose
|
|
pl_rest <- data.frame(lnC=log_dose, plotS=pl_restS, plotT=pl_restT)
|
|
|
|
pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="#4545BA",
|
|
inherit.aes = F) +
|
|
geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="#C2173F",
|
|
inherit.aes = F) +
|
|
{if (!is.null(truePL_df)) geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4,
|
|
inherit.aes = F) } +
|
|
{if (!is.null(truePL_df)) geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4,
|
|
inherit.aes = F) } +
|
|
labs(title = paste("restricted linear regression model"),
|
|
subtitle = paste("Regression on highlighted points")) +
|
|
theme(legend.position="none", axis.text = element_text(size=14))
|
|
pr3 <- pr2 + geom_point(circle, mapping=aes(x=log_dose, y=readout, shape=factor(isRef),
|
|
size=5, alpha=0.2), col=c("black"), inherit.aes = FALSE) +
|
|
scale_shape_manual(labels=c("test","reference"), values=c(21,21))
|
|
return(grid.arrange(p3,pr3,nrow=1))
|
|
}
|
|
|
|
|
|
|
|
#' Calculates the potency of 4PL PLA of all models model
|
|
#'
|
|
#' The gradient method is used for calculating the potency for a restricted model, an unrestricteed model,
|
|
#' and the log-transformations thereof.
|
|
#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used.
|
|
#'
|
|
#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(concentration).
|
|
#' @param PureErrFlag Boolean if Pure Error should be used.
|
|
#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results
|
|
#' 2) summary of unrestricted linear model 3) Slope difference +CIs
|
|
#' 4) summary of restricted linear model.
|
|
#' @export
|
|
#' @examples
|
|
#' ro_new <- data.frame(R_dil1 =c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229),
|
|
#' R_dil2=c( 10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608),
|
|
#' T_dil1=c( 10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791),
|
|
#' T_dil2=c( 11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655),
|
|
#' log_dose=c( -1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334))
|
|
#' PureErrF <- TRUE
|
|
#'
|
|
#'
|
|
#' pot4plFUNC(ro_new, PureErrF)
|
|
|
|
|
|
pot4plFUNC <- function(ro_new, PureErrFlag) {
|
|
all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
|
|
isRef <- rep(c(1,0),1,each=nrow(all_l)/2)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
all_l$readout[all_l$readout < 0] <- 0.01
|
|
all_l$readouttrans <- log(all_l$readout)
|
|
#browser()
|
|
CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)])
|
|
if (CORdat<0) SLOPE <- -1 else SLOPE <- 1
|
|
#
|
|
FITs <- Fitting_FUNC(ro_new, TransFlag = F)
|
|
if (!PureErrFlag) {
|
|
pot_est <- FITs[[3]]
|
|
potU_est <- FITs[[4]]
|
|
} else {
|
|
FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
|
|
meanPureErr <- FitAnova[4,3]
|
|
DFsPure <- FitAnova[4,1]
|
|
# SU_mr <- tryCatch({
|
|
# SU_mr <- summary(mr)
|
|
# }, error = function(msg) {
|
|
# return()
|
|
# })
|
|
SU_mr <- FITs[[1]]
|
|
SU_mrCoeff <- SU_mr$coefficients
|
|
# if (length(SU_mr)>1) {
|
|
# SU_mrCoeff <- SU_mr$coefficients
|
|
# } else { SU_mrCoeff <- rep(NA,5) }
|
|
#te2$cov.unscaled[1,1]*te2$sigma^2
|
|
#VCOV <- vcov(mr)
|
|
#V_V <- VCOV/SU_mr$sigma^2
|
|
V_V <- SU_mr$cov.unscaled
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
pot_est <- c(exp(SU_mrCoeff['r',1]),exp(SU_mrCoeff['r',1]-qt(0.975,DFsPure)*SEsPure['r']),
|
|
exp(SU_mrCoeff['r',1]+qt(0.975,DFsPure)*SEsPure['r']))
|
|
# unrestricted
|
|
SU_mu <- FITs[[2]]
|
|
s_muCoeff <- SU_mu$coefficients
|
|
#SU_mu <- summary(mu)
|
|
#VCOVu <- vcov(mu)
|
|
V_Vu <- SU_mu$cov.unscaled
|
|
SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
|
|
potU_est <- c(exp(s_muCoeff['r',1]),exp(s_muCoeff['r',1]-qt(0.975,DFsPure)*SEsPureU['r']),
|
|
+ exp(s_muCoeff['r',1]+qt(0.975,DFsPure)*SEsPureU['r']))
|
|
} # PureErrFlag
|
|
|
|
FITstrans <- Fitting_FUNC(ro_new, TransFlag = TRUE)
|
|
|
|
|
|
# pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic"))
|
|
# potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic"))
|
|
pot_est_log <- FITstrans[[3]]
|
|
potU_est_log <- FITstrans[[4]]
|
|
colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI")
|
|
colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI")
|
|
#browser()
|
|
su_mr_log <- FITstrans[[1]] #summary(mr_log)
|
|
Dat$RMSE_Rlog <- su_mr_log$sigma
|
|
su_mu_log <- FITstrans[[2]] # summary(mu_log)
|
|
Dat$RMSE_Ulog <- su_mu_log$sigma
|
|
Dat$up_lowAslog <- su_mu_log$coefficients[1,1] - su_mu_log$coefficients[4,1]
|
|
potALL <- rbind(round(pot_est,6), round(potU_est,6), round(pot_est_log,6), round(potU_est_log,6))
|
|
|
|
potALL2 <- cbind(c("restricted","unrestricted","ln-transformed restr","ln-transformed unrestr"), potALL)
|
|
return(potALL2)
|
|
}
|
|
|
|
#' Calculates logarithmized CIs for ratios
|
|
#'
|
|
#' Ratio CIs can be unbound, and not calulatable. Therefore the ratios are logarithmized.
|
|
#' The Covariance is not used, as it is 0 for comparisons between products.
|
|
#'
|
|
#' @param xs The estimate of the reference sample.
|
|
#' @param xt The estimate of the test sample.
|
|
#' @param se_xs The standard error of the estimate of the reference sample.
|
|
#' @param se_xt The standard error of the estimate of the test sample.
|
|
#' @param CoVar The covariance between test and reference estimate (set to 0).
|
|
#' @param DFs Degrees of freedom, either from pure error term (no of readouts - no of concentrations*2), or RMSE term (no of readouts - 8 parameters).
|
|
#' @param Conf The confidence of the Confidence Interval (default 95% which requires 0.975).
|
|
#' @returns A data-frame with the lower and upper CI in anti-log form.
|
|
#' @export
|
|
#' @examples
|
|
#' xs=2; xt=3.2; se_xt=0.34;se_xs=0.23; DFs=32-16
|
|
#'
|
|
#'
|
|
#' ParamCI_F(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975)
|
|
|
|
|
|
|
|
ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVar,DFs, Conf=0.975) {
|
|
log_xs <- log(abs(xs))
|
|
log_xt <- log(abs(xt))
|
|
var_log_xs <- (se_xs/xs)^2 # approximate variance of log(xs)
|
|
var_log_xt <- (se_xt/xt)^2
|
|
se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVar/(xs*xt)
|
|
|
|
lower_log_ratio <- log_xt-log_xs - qt(Conf,DFs)*se_log_ratio
|
|
upper_log_ratio <- log_xt-log_xs + qt(Conf,DFs)*se_log_ratio
|
|
ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio))
|
|
return(ci_ratio)
|
|
}
|
|
|
|
|
|
#' Calculates EQ tests and F-Tests for 4P
|
|
#'
|
|
#' CIs of asmptote ratios,slope ratios, asympt-difference ratio, and lower asympt-diff. calculated.
|
|
#' In addition the F-test on significance of regression and non-linearity are calculated.
|
|
#'
|
|
#' @param ro_new Data frame with the concentrations and the respective readouts.
|
|
#' @param Lim The limits to check if the relative CI is within the limits.
|
|
#' @param PureErrFlag Boolean if Pure Error should be used.
|
|
#' @returns A data-frame with the lower and upper CI in anti-log form.
|
|
#' @export
|
|
#' @examples
|
|
#'
|
|
#' dat <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
|
|
#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
|
|
#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
|
|
#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764))
|
|
#' Lim <- c(-1, 1, 0.005 , 2, 0.5, 2, 0.5,2,75,133, 75, 133)
|
|
#' PureErrF <- FALSE
|
|
#'
|
|
#' tests_FUNC(ro_new, Lim,PureErrF)
|
|
|
|
|
|
tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
|
|
|
all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
|
|
isRef <- rep(c(1,0),1,each=nrow(all_l)/2)
|
|
isSample <- rep(c(0,1),1,each=nrow(all_l)/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
all_l$readout[all_l$readout < 0] <- 0.01
|
|
#browser()
|
|
FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE)
|
|
if (is.character(FITs)) return(FITs) # if singularity
|
|
|
|
POTr_CI <- FITs[[3]][2:3]
|
|
potAll2 <- FITs[[3]]
|
|
|
|
smr <- FITs[[1]]
|
|
smu <- FITs[[2]]
|
|
FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
|
|
# pure error
|
|
pureSS <- FitAnova[4,2]
|
|
pureSS_df <- FitAnova[4,1]
|
|
meanPureErr <- FitAnova[4,3]
|
|
#vcovMU <- vcov(mu)
|
|
V_V <- smu$cov.unscaled
|
|
SEsPure <- sqrt(diag(V_V)*meanPureErr)
|
|
VCOVpure <- V_V*meanPureErr
|
|
DFsPure <- FitAnova[4,1]
|
|
|
|
|
|
|
|
testPOTr <- logical()
|
|
if (POTr_CI[1]*100>Lim[[9]] & POTr_CI[2]*100<Lim[[10]] ) testPOTr <- 0 else testPOTr <- 1
|
|
|
|
potAllU2 <- FITs[[4]]
|
|
test_c <- logical()
|
|
if((potAllU2[2] >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
|
|
|
|
#### 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)
|
|
# 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)
|
|
RSS_df <- AnovaDFs[5]
|
|
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
|
|
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 (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" }
|
|
|
|
test_a <- test_b <- test_d <- test_ad <- logical()
|
|
|
|
RSS_r <- round(sum(smr$residuals^2),5)
|
|
MSE_r <- RSS_r/(nrow(all_l)-5)
|
|
RMSE_r <- round(sqrt(MSE_r),6)
|
|
Dat$RMSE_r <- RMSE_r
|
|
Dat$RMSE_pure <- RMSE_pure
|
|
Dat$RMSE_unr <- round(RMSEunr,6)
|
|
|
|
coeffs <- smu$coefficients[,1]
|
|
#browser()
|
|
## EQ test on lower As diff
|
|
as <- coeffs["as"]
|
|
at <- coeffs["at"]
|
|
lAs_diff <- (at-as)
|
|
uCI_laDiff <- lAs_diff+qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2)
|
|
lCI_laDiff <- lAs_diff-qt(0.975,smu$df[2])*sqrt(smu$coefficients["ds",2]^2+smu$coefficients["dt",2]^2)
|
|
if (uCI_laDiff < Lim[[2]] & lCI_laDiff > Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1
|
|
|
|
#### EQ test on upper asymptote ratio ----
|
|
# as <- coeffs["as"]
|
|
# at <- coeffs["at"]
|
|
# uAsRatio <- compParm(potU, "a","/",display=F)
|
|
# uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2])
|
|
#browser()
|
|
ds <- smu$coefficients["ds",1]
|
|
dt <- smu$coefficients["dt",1]
|
|
if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds","ds"]) else se_ds <- smu$coefficients["ds",2]
|
|
if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt","dt"]) else se_dt <- smu$coefficients["dt",2]
|
|
if (PureErrFlag) CoVarlog_d <- VCOVpure["dt","ds"] else CoVarlog_d <- vcovMU["dt","ds"]
|
|
if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l)-8
|
|
uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.975)
|
|
if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1
|
|
estUppA <- round(at/as,5)
|
|
|
|
Dat$uAsCI <- uAsCI2
|
|
|
|
#### EQ test on slope ratio ----
|
|
# bs <- coeffs["bs"]
|
|
# bt <- coeffs["bt"]
|
|
# slopeRatio <- compParm(potU, "b","/",display=F)
|
|
# slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2])
|
|
|
|
bs <- smu$coefficients["bs",1]
|
|
bt <- smu$coefficients["bt",1]
|
|
if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs","bs"]) else se_bs <- smu$coefficients["bs",2]
|
|
if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt","bt"]) else se_bt <- smu$coefficients["bt",2]
|
|
if (PureErrFlag) CoVarlog_b <- VCOVpure["bt","bs"] else CoVarlog_b <- vcovMU["bt","bs"]
|
|
slopeCI2 <- ParamCI_F(bt,bs,se_bt, se_bs,CoVarlog_b, DFs, Conf=0.975)
|
|
if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1
|
|
estUppA <- round(at/as,5)
|
|
|
|
Dat$slopeRatioCI <- slopeCI2
|
|
|
|
#### EQ test on lower As ratio ----
|
|
|
|
# lAsRatio <- compParm(potU, "d","/",display=F)
|
|
# slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2])
|
|
|
|
as <- smu$coefficients["as",1]
|
|
at <- smu$coefficients["at",1]
|
|
if (PureErrFlag) se_as <- sqrt(VCOVpure["as","as"]) else se_as <- smu$coefficients["as",2]
|
|
if (PureErrFlag) se_at <- sqrt(VCOVpure["at","at"]) else se_at <- smu$coefficients["at",2]
|
|
if (PureErrFlag) CoVarlog_a <- VCOVpure["at","as"] else CoVarlog_a <- vcovMU["at","as"]
|
|
lAsCI2 <- ParamCI_F(at,as,se_at, se_as,CoVarlog_a, DFs, Conf=0.975)
|
|
if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1
|
|
estLowA <- round(at/as,5)
|
|
|
|
Dat$lAsCI <- lAsCI2
|
|
|
|
#### EQtest on ratio of As difference ----
|
|
AsDiffRatio <- (dt-at)/(ds-as)
|
|
|
|
dt_at <- (dt-at)
|
|
ds_as <- (ds-as)
|
|
se_ds_asPure <- sqrt(VCOVpure["as","as"]+VCOVpure["ds","ds"]-2*VCOVpure["as","ds"])
|
|
se_dt_atPure <- sqrt(VCOVpure["at","at"]+VCOVpure["dt","dt"]-2*VCOVpure["at","dt"])
|
|
se_ds_asRMSE <- sqrt(vcovMU["as","as"]+vcovMU["ds","ds"]-2*vcovMU["as","ds"])
|
|
se_dt_atRMSE <- sqrt(vcovMU["at","at"]+vcovMU["dt","dt"]-2*vcovMU["at","dt"])
|
|
if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE
|
|
if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE
|
|
|
|
AsDiffCI2 <- ParamCI_F(dt_at,ds_as,se_dt_at, se_ds_as,CoVar=0, DFs, Conf=0.975)
|
|
if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1
|
|
estLowA <- round(at/as,5)
|
|
|
|
Dat$up_lowAs <- abs(ds-as)
|
|
|
|
lowerCIlowerA <- lAsCI2[1]; lowerCIupperA <- uAsCI2[1]; upperCIlowerA <- lAsCI2[2]; upperCIupperA <- uAsCI2[2]
|
|
test_lowA <- test_d; test_uppA <- test_a
|
|
#browser()
|
|
res_tab <- data.frame(test= c("F-test on sign. of regression*",
|
|
"EQ test on lower asymptotes difference",
|
|
"EQ test ratio of lower asymptotes",
|
|
"EQ test ratio of Hill slopes",
|
|
"EQ test ratio of upper asymptotes",
|
|
"F-test on non-linearity*",
|
|
"EQ test ratio of asymptote difference",
|
|
"geom. rel. CI restr. model",
|
|
"geom. rel. CI unrestr. model"),
|
|
test_results = c(ifelse(p_F_regr<0.05,0,1), test_la_diff, test_lowA, test_b, test_uppA,
|
|
ifelse(p_F_nonlin>1,1, ifelse(p_F_nonlin<0.05,1,0)), test_ad,
|
|
testPOTr, test_c),
|
|
estimate = c(round(p_F_regr, 3), round(lAs_diff, 5),
|
|
estLowA, round(bs/bt,5), estUppA, p_F_nonlin,
|
|
round(dt_at/ds_as, 5), round(potAll2[1]*100,2),round(potAllU2[1]*100,2)),
|
|
lower_limit = c("-",Lim[[1]],Lim[[3]],Lim[[5]],Lim[[7]],"-",Lim[[11]],Lim[[9]],Lim[[9]]),
|
|
upper_limit = c("-",Lim[[2]],Lim[[4]],Lim[[6]],Lim[[8]],"-",Lim[[12]],Lim[[10]],Lim[[10]]),
|
|
lower_CI = c(RMSE_r, round(lCI_laDiff,3), round(lAsCI2[1],5), round(slopeCI2[1],5),
|
|
round(uAsCI2[1],5), "-", round(AsDiffCI2[1],5), round(potAll2[2],2), round(potAllU2[2],2)),
|
|
upper_CI = c(RMSE_pure, round(uCI_laDiff,3), round(lAsCI2[2],5), round(slopeCI2[2],5),
|
|
round(uAsCI2[2],5), "-", round(AsDiffCI2[2],5), round(potAll2[3],2), round(potAllU2[3],2))
|
|
)
|
|
return(res_tab)
|
|
}
|
|
|
|
#' ANOVA unrestricted model
|
|
#'
|
|
#' ANOVA with unrestricted model.
|
|
#'
|
|
#'
|
|
#' @param ro_new Data frame with the concentrations and the respective readouts.
|
|
#' @returns A data-frame with the analysis of variance.
|
|
#' @export
|
|
#' @examples
|
|
#'
|
|
#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
|
|
#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
|
|
#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
|
|
#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764))
|
|
#'
|
|
#'
|
|
#' ANOVA4plUnresfunc(ro_new)
|
|
|
|
|
|
ANOVA4plUnresfunc <- function(ro_new) {
|
|
all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name="replname", value.name = "readout")
|
|
all_len <- nrow(all_l)
|
|
isRef <- rep(c(1,0),1,each=all_len/2)
|
|
isSample <- rep(c(0,1),1,each=all_len/2)
|
|
all_l$isRef <- isRef
|
|
all_l$isSample <- isSample
|
|
all_l$Conc <- exp(all_l$log_dose)
|
|
all_l$readout[all_l$readout < 0] <- 0.01
|
|
|
|
FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE)
|
|
smr <- FITs[[1]]
|
|
smu <- FITs[[2]]
|
|
smrPREDs <- FITs[[5]]
|
|
smuPREDs <- FITs[[6]]
|
|
# pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
|
|
# pmodels=data.frame(1,1,1,isSample))
|
|
# potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
|
|
# pmodels=data.frame(isSample, isSample,isSample,isSample))
|
|
SStreat <- round(sum((smuPREDs - mean(all_l$readout))^2),5)
|
|
SStreat_df <- length(unique(all_l$log_dose))-1
|
|
SSregr <- round(sum((smrPREDs-mean(all_l$readout))^2),5)
|
|
## Non-parallel
|
|
SSnonparallel <- round(sum(smr$residuals^2) - sum(smu$residuals^2),5)
|
|
## Preparation
|
|
SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5)
|
|
## Resid Err
|
|
RSS <- round(sum(smu$residuals^2),5)
|
|
RSS_df <- nrow(all_l)-SStreat_df-1
|
|
FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
|
|
# PureErr
|
|
SSE <- round(FitAnova[4,3],5)
|
|
SSE_df <- FitAnova[4,1]
|
|
# Non-Linearity
|
|
SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - smuPREDs)^2),4)
|
|
LoF_df <- FitAnova[1,1]+FitAnova[2,1]
|
|
## Total
|
|
SStot <- round(sum((all_l$readout -mean(all_l$readout))^2),5)
|
|
MSE <- RSS/RSS_df
|
|
noConc <- length(unique(all_l$Conc))
|
|
AnovaDFs <- c(noConc-1, 1,3,noConc-4-1, nrow(all_l)-noConc, noConc, nrow(all_l)-noConc-noConc, nrow(all_l)-1)
|
|
p_SStreat <- round(pf((SStreat/AnovaDFs[1])/MSE, AnovaDFs[1],RSS_df, lower.tail = F),3)
|
|
p_SSprep <- round(pf((SSprep/AnovaDFs[2])/MSE, AnovaDFs[2],RSS_df, lower.tail = F),3)
|
|
p_SSregr <- round(pf((SSregr/AnovaDFs[3])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3)
|
|
p_SSnonp <- round(pf((SSnonparallel/AnovaDFs[4])/MSE, AnovaDFs[3],RSS_df, lower.tail = F),3)
|
|
p_SSLoF <- round(pf((SSnonlin/LoF_df)/(SSE/SSE_df), LoF_df,SSE_df, lower.tail = F),5)
|
|
|
|
ANOVAtab <- data.frame(Source = c("Treatment","Preparation","Regression",
|
|
"Non-Parallelism","Residual Error","Non-linearity",
|
|
"Pure Error","Total"),
|
|
DF = AnovaDFs,
|
|
SumSquares = c(SStreat, SSprep,SSregr, SSnonparallel,
|
|
RSS, SSnonlin,SSE, SStot),
|
|
MeanSquares = c(round(SStreat/AnovaDFs[1],3), SSprep, round(SStreat/AnovaDFs[3],3),round(SSnonparallel/AnovaDFs[4],3),
|
|
round(MSE,5), round(SSnonlin/LoF_df,5), round(SSE/SSE_df,5),""),
|
|
"F-value" = c(round((SStreat/AnovaDFs[1])/MSE,5), round((SSprep/AnovaDFs[2])/MSE,5),
|
|
round((SSregr/AnovaDFs[3])/MSE,5),round((SSnonparallel/AnovaDFs[4])/MSE,5),
|
|
"",round((SSnonlin/LoF_df)/(SSE/SSE_df),5),"",""),
|
|
"p_value" = c(round(p_SStreat,3), p_SSprep, round(p_SSregr,3), p_SSnonp,"",p_SSLoF,"","")
|
|
)
|
|
|
|
return(ANOVAtab)
|
|
}
|
|
|
|
#' Calculates descriptive statistics per concentration
|
|
#'
|
|
#' Mean, SD, and CV% per concentration and product..
|
|
#'
|
|
#'
|
|
#' @param ro_new Data frame with the concentrations and the respective readouts.
|
|
#' @returns A data-frame with the concentrations and metadata.
|
|
#' @export
|
|
#' @examples
|
|
#'
|
|
#' ro_new <- data.frame(REF1=c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2=c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
|
|
#' REF3=c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
|
|
#' TEST1=c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2=c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
|
|
#' TEST3=c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose=c(5.01,3.401,2.708,2.015,1.32176,0.62861,-0.0645385,-1.6739764))
|
|
#'
|
|
#'
|
|
#' perConcTab(ro_new, noDilSeries=3)
|
|
|
|
|
|
|
|
perConcTab <- function(ro_new, noDilSeries) {
|
|
Reftab <- ro_new[,c(1:noDilSeries)]
|
|
Testtab <- ro_new[,c((noDilSeries+1):(2*noDilSeries))]
|
|
tReftab <- t(Reftab)
|
|
colnames(tReftab) <- round(ro_new[,ncol(ro_new)],5)
|
|
|
|
avs <- apply(tReftab,2,mean)
|
|
sds <- apply(tReftab,2,sd)
|
|
cv <- sds/avs*100
|
|
tReftab2 <- rbind(tReftab, avs,sds,cv)
|
|
|
|
tTesttab <- t(Testtab)
|
|
colnames(tTesttab) <- round(ro_new[,ncol(ro_new)],5)
|
|
|
|
avs_test <- apply(tTesttab,2,mean)
|
|
sds_test <- apply(tTesttab,2,sd)
|
|
cv_test <- sds_test/avs_test*100
|
|
tTesttab2 <- rbind(tTesttab, avs_test,sds_test,cv_test)
|
|
concTab <- rbind(tReftab2, tTesttab2)
|
|
return(concTab)
|
|
|
|
}
|
|
|
|
#' Calculates dilution series.
|
|
#'
|
|
#' Recursive function to calcuate a geometric dilution series.
|
|
#'
|
|
#' @param x Starting concentration.
|
|
#' @param Div Divisor to get next concentration.
|
|
#' @param N Number of dilution step, increases.
|
|
#' @param res Vector of results.
|
|
#' @param noDil Final number of concentrations -1 (the initial one).
|
|
#' @returns A data-frame with the concentrations and metadata.
|
|
#' @export
|
|
#' @examples
|
|
#'
|
|
#' x <- 1; Div <- 3;N <- 0; res <- c(); noDil <- 7
|
|
#'
|
|
#' divFUN(x,Div,N,res,noDil)
|
|
|
|
|
|
|
|
divFUN <- function(x,Div,N,res,noDil) {
|
|
N <- N+1
|
|
y <- x/Div
|
|
res <- c(res,y)
|
|
if (N==noDil) { return(res) }
|
|
divFUN(y,Div,N,res,noDil)
|
|
}
|
|
|