bugfix: 4PL fit made robust
This commit is contained in:
+83
-34
@@ -23,6 +23,25 @@ library(car)
|
||||
library(dplyr)
|
||||
library(scales)
|
||||
|
||||
#' Estimate correlations
|
||||
#'
|
||||
#' returns the correlation of 2 vectors
|
||||
#'
|
||||
#' @param vec1 The 1st vector.
|
||||
#' @param vec2 The 2nd vector.
|
||||
#' @returns A float as correlatioin estimate
|
||||
#' @export
|
||||
#' @examples
|
||||
#' suppressMessages(source("../../dev/setup.R"))
|
||||
#' vector1 <- c(1,2,3,4,5)
|
||||
#' vector2 <- c(5.1,4.3,NA,1.9,1.2)
|
||||
#' te <- COR_FUNC(vector1,vector2)
|
||||
#' print(te)
|
||||
COR_FUNC <- function(vec1, vec2) {
|
||||
df <- data.frame(v1 = vec1, v2 = vec2)
|
||||
df2 <- df[complete.cases(df),]
|
||||
return(cor(df2[,1],df2[,2]))
|
||||
}
|
||||
|
||||
#' Levenberg Marquard fit of 4 pl
|
||||
#'
|
||||
@@ -46,8 +65,9 @@ library(scales)
|
||||
#' Dat <- list()
|
||||
#' te <- Fitting_FUNC(dat, TransF)
|
||||
#' print(te)
|
||||
Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
CORro <- cor(ro_new[, 1], ro_new[, ncol(ro_new)])
|
||||
Fitting_FUNC <- function(ro_new, TransFlag = FALSE, nameWS ="") {
|
||||
#browser()
|
||||
CORro <- COR_FUNC(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)
|
||||
@@ -69,6 +89,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
},
|
||||
warning = function(e) {
|
||||
mr <<- "In nlsModel singular gradient matrix"
|
||||
|
||||
}
|
||||
)
|
||||
# Stop if singular gradient matrix
|
||||
@@ -82,6 +103,12 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
},
|
||||
error = function(err) {
|
||||
s_mr <- NULL
|
||||
showModal(modalDialog(
|
||||
title = " fit",
|
||||
paste("fit not possible: EC50 outside concentration range for dataset", nameWS),
|
||||
easyClose = TRUE,
|
||||
footer = NULL
|
||||
))
|
||||
}
|
||||
)
|
||||
} else {
|
||||
@@ -98,7 +125,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
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
|
||||
at = min(ro_new[, 4]), bt = SLOPE, dt = max(ro_new[, 4]), r = 0
|
||||
)
|
||||
tryCatch(
|
||||
{
|
||||
@@ -120,13 +147,19 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
summary(mu)
|
||||
},
|
||||
error = function(msg) {
|
||||
showModal(modalDialog(
|
||||
title = "4PL fit",
|
||||
paste("fit not possible: EC50 outside concentration range for dataset", nameWS),
|
||||
easyClose = TRUE,
|
||||
footer = NULL
|
||||
))
|
||||
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
|
||||
at = log(min(ro_new[, 4])), bt = SLOPE, dt = log(max(ro_new[, 4])), r = 0
|
||||
)
|
||||
tryCatch(
|
||||
{
|
||||
@@ -152,16 +185,27 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
}
|
||||
)
|
||||
}
|
||||
|
||||
#browser()
|
||||
if (!TransFlag) {
|
||||
pot_est <- exp(confintd(mr, "r", method = "asymptotic"))
|
||||
potU_est <- exp(confintd(mu, "r", method = "asymptotic"))
|
||||
PRED <- predict(mr)
|
||||
PREDu <- predict(mu)
|
||||
#browser()
|
||||
if (length(s_mr) ==1 | length(Sum_u) ==1) {
|
||||
return("failed")
|
||||
} else {
|
||||
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)
|
||||
if (length(s_mr) ==1 | length(Sum_u) ==1) {
|
||||
return("failed")
|
||||
}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))
|
||||
}
|
||||
@@ -212,7 +256,7 @@ Fitting_FUNC <- function(ro_new, TransFlag = FALSE) {
|
||||
#' p <- plotSingularity(dat)
|
||||
#' print(p)
|
||||
plotSingularity <- function(dat) { # sigmoid,det_sig,
|
||||
CORdat <- cor(dat[, 1], dat[, ncol(dat)])
|
||||
CORdat <- COR_FUNC(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)
|
||||
@@ -272,7 +316,7 @@ plotSingularity <- function(dat) { # sigmoid,det_sig,
|
||||
#' p <- plot_f(dat, TransFlag)
|
||||
#' print(p)
|
||||
plot_f <- function(dat, TransFlag = FALSE) { # sigmoid,det_sig,
|
||||
CORdat <- cor(dat[, 1], dat[, ncol(dat)])
|
||||
CORdat <- COR_FUNC(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)
|
||||
@@ -1018,7 +1062,7 @@ pot4plFUNC <- function(ro_new, PureErrFlag) {
|
||||
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)])
|
||||
CORdat <- COR_FUNC(ro_new[, 1], ro_new[, ncol(ro_new)])
|
||||
if (CORdat < 0) SLOPE <- -1 else SLOPE <- 1
|
||||
#
|
||||
FITs <- Fitting_FUNC(ro_new, TransFlag = FALSE)
|
||||
@@ -1112,8 +1156,8 @@ ParamCI_F <- function(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975) {
|
||||
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_xs - log_xt - qt(Conf, DFs) * se_log_ratio
|
||||
upper_log_ratio <- log_xs - log_xt + qt(Conf, DFs) * se_log_ratio
|
||||
ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio))
|
||||
return(ci_ratio)
|
||||
}
|
||||
@@ -1143,6 +1187,9 @@ ParamCI_F <- function(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975) {
|
||||
#'
|
||||
#' tests_FUNC(ro_new=dat, Lim, PureErrF)
|
||||
tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
|
||||
DatL <- list()
|
||||
|
||||
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)
|
||||
@@ -1228,9 +1275,9 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
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)
|
||||
DatL$RMSE_r <- RMSE_r
|
||||
DatL$RMSE_pure <- RMSE_pure
|
||||
DatL$RMSE_unr <- round(RMSEunr, 6)
|
||||
|
||||
coeffs <- smu$coefficients[, 1]
|
||||
# browser()
|
||||
@@ -1242,6 +1289,7 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
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
|
||||
|
||||
#browser()
|
||||
#### EQ test on upper asymptote ratio ----
|
||||
# as <- coeffs["as"]
|
||||
# at <- coeffs["at"]
|
||||
@@ -1254,12 +1302,13 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
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)
|
||||
uAsCI2 <- ParamCI_F(ds, dt, 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
|
||||
estUppA <- round(dt / ds, 5)
|
||||
|
||||
DatL$uAsCI <- uAsCI2
|
||||
# browser()
|
||||
|
||||
#### EQ test on slope ratio ----
|
||||
# bs <- coeffs["bs"]
|
||||
# bt <- coeffs["bt"]
|
||||
@@ -1271,11 +1320,11 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
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)
|
||||
slopeCI2 <- ParamCI_F(bs, bt, 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)
|
||||
estSlope <- round(abs(bt) / abs(bs), 5)
|
||||
|
||||
Dat$slopeRatioCI <- slopeCI2
|
||||
DatL$slopeRatioCI <- slopeCI2
|
||||
|
||||
#### EQ test on lower As ratio ----
|
||||
|
||||
@@ -1287,11 +1336,11 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
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)
|
||||
lAsCI2 <- ParamCI_F(as, at, 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
|
||||
DatL$lAsCI <- lAsCI2
|
||||
|
||||
#### EQtest on ratio of As difference ----
|
||||
AsDiffRatio <- (dt - at) / (ds - as)
|
||||
@@ -1305,11 +1354,11 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
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( ds_as,dt_at, 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)
|
||||
estDiffA <- round(dt_at /ds_as, 5)
|
||||
|
||||
Dat$up_lowAs <- abs(ds - as)
|
||||
Dat$estDiffA <- estDiffA
|
||||
|
||||
lowerCIlowerA <- lAsCI2[1]
|
||||
lowerCIupperA <- uAsCI2[1]
|
||||
@@ -1337,8 +1386,8 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
|
||||
),
|
||||
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)
|
||||
estLowA, estSlope, estUppA, p_F_nonlin,
|
||||
estDiffA, 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]]),
|
||||
|
||||
Reference in New Issue
Block a user