Catch singularities when fitting 4pl, SCRUM jobs listed

This commit is contained in:
2026-05-17 11:13:02 +02:00
parent ec13d95387
commit 9ff1a360d4
5 changed files with 136 additions and 103 deletions
+64 -13
View File
@@ -37,10 +37,17 @@ Fitting_FUNC <- function(ro_new, TransFlag=F) {
if (!TransFlag) {
startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0)
mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))),
mr <- tryCatch({gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))),
data=all_l2,
start=startlist,#trace=T,
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)
},
@@ -103,6 +110,40 @@ Fitting_FUNC <- function(ro_new, TransFlag=F) {
return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu))
}
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
#'
@@ -680,10 +721,18 @@ PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) {
log_dose <- unique(all_l2$log_dose)
seq_x <- seq(min(log_dose), max(log_dose),0.1)
SAMPLEtrue <- sigmoid[5] + (sigmoid[7]-sigmoid[5])/(1+exp(sigmoid[6]*((sigmoid[4]-sigmoid[8]-seq_x))))
REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[2]*((sigmoid[4]-seq_x))))
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
}
truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue)
p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) +
geom_point(size=2) +
@@ -697,10 +746,10 @@ PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) {
inherit.aes = F) +
geom_line(data=pl_df,aes(x=lnC,y=plotT),color="#C2173F",
inherit.aes = F) +
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4,
inherit.aes = F) +
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4,
inherit.aes = F) +
{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),
@@ -717,10 +766,10 @@ PlotLinPLA_FUNC <-function(circle, sigmoid, all_l2, pl_df, indS, indT) {
inherit.aes = F) +
geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="#C2173F",
inherit.aes = F) +
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="#C2173F", linetype=2,alpha=0.4,
inherit.aes = F) +
geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="#4545BA", linetype=2,alpha=0.4,
inherit.aes = F) +
{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))
@@ -895,6 +944,8 @@ tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
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]]