ROUT outlier testing added
This commit is contained in:
+267
@@ -0,0 +1,267 @@
|
||||
################################################################################
|
||||
# F.Innerbichler
|
||||
# Jun 2026
|
||||
# Robust Outlier detection
|
||||
################################################################################
|
||||
|
||||
# Code acc. to Mutulsky and Brown hpptps://pubmed.ncbi.nlm.nih.gov/16526949/
|
||||
# except line 175 p <- 2*(1-pt(z, length(residuals)-5)) was exchanged with p <- 2*(1-pcauchy(z)) to be robust
|
||||
|
||||
library(minpack.lm)
|
||||
library(ggplot2)
|
||||
library(rstudioapi)
|
||||
|
||||
# Getting path for current open file
|
||||
current_path = rstudioapi::getActiveDocumentContext()$path
|
||||
setwd(dirname(current_path))
|
||||
|
||||
outsPlot_FUN <- function(all_l,OUTs_, TS, PROC = "not specified", PROBE="", Q=0.01, par) {
|
||||
#browser()
|
||||
#all_l_rout <- all_l[KEEP,]
|
||||
klein <- min(log(all_l$conc))
|
||||
gross <- max(log(all_l$conc))
|
||||
x_seq <- seq(klein, gross, (gross-klein)/100)
|
||||
samTRUE <- f4pl(x=x_seq, bottom=par["a"], top=par["d"],hill=par["b"], logEC50=par["cs"]-par["r"])
|
||||
refTRUE <- f4pl(x=x_seq, bottom=par["a"], top=par["d"],hill=par["b"], logEC50=par["cs"])
|
||||
|
||||
pl_T <- data.frame(cbind(x_seq, refTRUE, samTRUE))
|
||||
p <- ggplot(all_l) +
|
||||
geom_point(aes(x=log(conc), y=y, shape=factor(isRef)), size=2) +
|
||||
theme_bw()
|
||||
|
||||
p2 <- p + geom_point(OUTs_, mapping=aes(x=log(conc), y=y, shape=factor(isRef)),
|
||||
size=5, color="violetred", stroke=2, inherit.aes = FALSE) +
|
||||
scale_shape_manual(label=c("R",TS), values = c(21,25)) +
|
||||
ggtitle(paste(PROC, PROBE, TS, "Threshold=", Q)) +
|
||||
theme(legend.title=element_blank())
|
||||
p3 <- p2 + geom_line(data=pl_T, aes(x=x_seq, y=refTRUE), color="blue", inherit.aes = F) +
|
||||
geom_line(data=pl_T, aes(x=x_seq, y=samTRUE), color="red", inherit.aes = F)
|
||||
p3
|
||||
|
||||
}
|
||||
|
||||
#### 4PL model ----
|
||||
f4pl <- function(x, bottom, top, logEC50, hill) {
|
||||
bottom+(top-bottom)/(1+exp((logEC50-x)*hill))
|
||||
}
|
||||
|
||||
#### residuals for joint REF/TEST fit ----
|
||||
resid_4pl_joint <- function(par,x,y, is_ref) {
|
||||
bottom <- par["bottom"]
|
||||
top <- par["top"]
|
||||
hill <- par["hill"]
|
||||
le50_r <- par["logEC50_ref"]
|
||||
le50_t <- par["logEC50_test"]
|
||||
|
||||
le50 <- ifelse(is_ref, le50_r, le50_t)
|
||||
yhat <- f4pl(x, bottom, top, le50, hill)
|
||||
y - yhat
|
||||
}
|
||||
|
||||
|
||||
jac_4pl_joint <- function(par, x,y,is_ref) {
|
||||
bottom <- par["bottom"]
|
||||
top <- par["top"]
|
||||
hill <- par["hill"]
|
||||
le50_r <- par["logEC50_ref"]
|
||||
le50_t <- par["logEC50_test"]
|
||||
|
||||
le50 <- ifelse(is_ref, le50_r, le50_t)
|
||||
|
||||
# useful intermediates
|
||||
|
||||
texp <- exp((le50-x)*hill) # t = 10^((le50-x)*hill))
|
||||
den <- (1-texp)
|
||||
frac <- (top-bottom)/den
|
||||
|
||||
dyhat_dbottom <- 1-1/den
|
||||
dyhat_dtop <- 1/den
|
||||
|
||||
d_invden_dle50 <- -(texp*log(2.718282)*hill)/(den^2)
|
||||
|
||||
d_invden_dhill <- -(texp*log(2.718282)*(le50-x))/(den^2)
|
||||
|
||||
dyhat_dle50 <- (top-bottom)* d_invden_dle50
|
||||
dyhat_dhill <- (top-bottom)* d_invden_dhill
|
||||
|
||||
J <- cbind(
|
||||
bottom = -dyhat_dbottom,
|
||||
topm = -dyhat_dtop,
|
||||
hill = -dyhat_dhill,
|
||||
logEC50_ref = ifelse(is_ref, -dyhat_dle50,0),
|
||||
logEC50_test = ifelse(!is_ref, -dyhat_dle50,0)
|
||||
)
|
||||
J
|
||||
}
|
||||
|
||||
|
||||
robust_fit_irls <- function(par_start, x,y, is_ref,
|
||||
max_outer = 30, control = nls.lm.control(maxiter=200, ftol=1e-10, ptol=1e-10)) {
|
||||
par <- par_start
|
||||
#browser()
|
||||
for (k in seq_len(max_outer)) {
|
||||
r <- resid_4pl_joint(par,x,y,is_ref)
|
||||
|
||||
# Mot Brown
|
||||
rSD <- quantile(r, 0.6827)*length(is_ref)/(length(is_ref)-5)
|
||||
# Robust scale (MAD); fallback if MAD ~0
|
||||
# s <- median(abs(r))/0.6745
|
||||
|
||||
# Lorentzian / Cauchy weights acc. Mot/Brown
|
||||
w <- 1/log(1+(abs(r)/rSD)^2)
|
||||
# weighted residual function lor LM: sqrt(w)*r
|
||||
fn_w <- function(p) {
|
||||
rr <- resid_4pl_joint(p,x,y,is_ref)
|
||||
sqrt(w)*rr
|
||||
# w*rr
|
||||
}
|
||||
jac_w <- function(p) {
|
||||
JJ <- jac_4pl_joint(p,x,y,is_ref)
|
||||
JJ*sqrt(w) # row-wise scaling
|
||||
# JJ*w
|
||||
}
|
||||
|
||||
#browser()
|
||||
fit <- tryCatch(nls.lm(par=par, fn=fn_w, jac=jac_w, control=control),
|
||||
error=function(e) {
|
||||
paste("error at robust fit irls", k)
|
||||
# tau <<- tau *10
|
||||
#Return zero delta to avoid NaN updates
|
||||
return(rep(0, length(par)))
|
||||
})
|
||||
|
||||
if (max(abs(fit$par -par)) < 1e-8) {
|
||||
par <- fit$par
|
||||
break
|
||||
}
|
||||
par <- fit$par
|
||||
}
|
||||
r_final <- resid_4pl_joint(par,x,y,is_ref)
|
||||
# s_final <- median(abs(r_final))/0.6745
|
||||
rSD <- quantile(r, 0.6827)*length(r)/(length(r)-5)
|
||||
# Robust scale (MAD); fallback if MAD ~0
|
||||
# s <- median(abs(r))/0.6745
|
||||
list(par=par, residuals=r_final, scale=rSD)
|
||||
|
||||
}
|
||||
|
||||
# ROUT outlier detection via BH-FDR [1] https://cran.r-project.org/web//packages/minpack.lm/refman/minpack.lm.html)[6]
|
||||
|
||||
rout_detect <- function(residuals, scale, Q=0.01) {
|
||||
z <- abs(residuals/scale)
|
||||
p <- 2*(1 - pcauchy(z))
|
||||
# p <- 2*(1-pt(z, length(residuals)-5))
|
||||
|
||||
o <- order(p)
|
||||
p_sorted <- p[o]
|
||||
m <- length(p)
|
||||
# thresh <- (seq_len(m)/m) # *Q
|
||||
keep <- p_sorted <= Q # thresh
|
||||
(out_idx <- o[keep])
|
||||
sort(out_idx)
|
||||
|
||||
}
|
||||
|
||||
# Final OLS LM refit after removing outliers
|
||||
|
||||
ols_refit <- function(par_start, x,y,is_ref, keep_idx,
|
||||
control = nls.lm.control(maxiter=400, ftol=1e-12, ptol=1e-12)) {
|
||||
xk <- x[keep_idx]; yk <- y[keep_idx]; rk <- is_ref[keep_idx]
|
||||
|
||||
fn_w <- function(p) resid_4pl_joint(p,xk,yk,rk)
|
||||
jac_w <- function(p) jac_4pl_joint(p,xk,yk,rk)
|
||||
#browser()
|
||||
fit <- tryCatch(nls.lm(par=par_start, fn=fn_w, jac=jac_w, control=control),
|
||||
error=function(e) {
|
||||
paste("error at ols_refits")
|
||||
# tau <<- tau *10
|
||||
#Return zero delta to avoid NaN updates
|
||||
return(rep(0, length(par)))
|
||||
})
|
||||
|
||||
list(par = fit$par, residuals = fn_w(fit$par), keep=keep_idx)
|
||||
}
|
||||
|
||||
#. high-level ROUT 4PL analysis
|
||||
|
||||
rout_4pl_potency <- function(df, Q=0.01, par_start=NULL, max_outer = 30) {
|
||||
stopifnot(all(c("sample","conc","y") %in% names(df))) # "rep",
|
||||
#browser()
|
||||
x <- log(df$conc)
|
||||
y <- df$y
|
||||
is_ref <- df$isRef == 1
|
||||
|
||||
# start values if not provided
|
||||
if (is.null(par_start)) {
|
||||
bottom0 = min(y, na.rm = T)
|
||||
top0 = max(y, na.rm = T)
|
||||
if (df$y[1]>df$y[8]) hill0 <- 1 else hill0 <- -1
|
||||
le50_r0 <- median(x[is_ref], na.rm = T)
|
||||
le50_t0 <- median(x[!is_ref], na.rm = T)
|
||||
par_start <- c(bottom=bottom0, top=top0, hill=hill0, logEC50_ref = le50_r0 , logEC50_test = le50_t0)
|
||||
}
|
||||
# Step1 : robust fit
|
||||
robust <- robust_fit_irls(par_start, x,y,is_ref, max_outer = max_outer)
|
||||
# Step 2: outlier detection
|
||||
out_idx <- rout_detect(robust$residuals, robust$scale, Q=Q)
|
||||
keep_idx <- setdiff(seq_len(nrow(df)), out_idx)
|
||||
|
||||
# Step 3: OLS refit
|
||||
ols <- ols_refit(robust$par, x,y,is_ref, keep_idx)
|
||||
|
||||
# Relative potency
|
||||
rp <- exp(ols$par["logEC50_ref"] - ols$par["logEC50_test"])*100
|
||||
|
||||
list(Q=Q, outliers=out_idx,
|
||||
kept=keep_idx,
|
||||
par_robust = robust$par,
|
||||
par_final = ols$par,
|
||||
rp_percent = rp,
|
||||
df=df)
|
||||
|
||||
}
|
||||
|
||||
# sheets <- openxlsx::getSheetNames("~/plateflow/outlierFUB.xlsx")
|
||||
# all_dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile="~/plateflow/outlierFUB.xlsx")
|
||||
# names(all_dat) <- sheets
|
||||
#
|
||||
# Plate <- 1
|
||||
#
|
||||
# PlanteN <- sheets[Plate]
|
||||
# NoDils <- 8
|
||||
# Nreps <- 3
|
||||
# DAT <- all_dat[[Plate]]
|
||||
# colnames(DAT) <- DAT[2,]
|
||||
# REF <- DAT[3:10,1:4]
|
||||
# SAM1 <- DAT[13:20,1:4]
|
||||
# SAM2 <- DAT[23:30,1:4]
|
||||
# SAM3 <- DAT[33:40,1:4]
|
||||
#
|
||||
# REF_ <- sapply(REF, as.numeric)
|
||||
# colnames(REF_) <- colnames(REF)
|
||||
# SAM1_ <- sapply(SAM1, as.numeric)
|
||||
# colnames(SAM1_) <- colnames(SAM1)
|
||||
# SAM2_ <- sapply(SAM2, as.numeric)
|
||||
# colnames(SAM2_) <- colnames(SAM2)
|
||||
# SAM3_ <- sapply(SAM3, as.numeric)
|
||||
# colnames(SAM3_) <- colnames(SAM3)
|
||||
#
|
||||
# plot_df <- rbind(REF_, SAM1_, SAM2_,SAM3_)
|
||||
#
|
||||
# plot_df_ <- cbind(as.data.frame(plot_df), sample=c(rep("R", NoDils),rep("S1R", NoDils),rep("S2", NoDils),rep("S3", NoDils)))
|
||||
# test_df <- plot_df_[1:16,]
|
||||
# test_df2 <- plot_df_[c(1:8, 17:24),]
|
||||
# test_df3 <- plot_df_[c(1:8, 25:32),]
|
||||
#
|
||||
# all_l <- melt(test_df, id.vars = c("sample","Dose"), variable.name="replname", value.name="readout")
|
||||
# colnames(all_l) <- c("sample","conc","rep","y")
|
||||
#
|
||||
# for (Q in c(0.01,0.015,0.02)) {
|
||||
# res <- rout_4pl_potency(all_l, Q)
|
||||
# OUTs_ <- all_l[res$outliers,]
|
||||
# print(outsPlot_FUN(all_l, OUTs_, TS=all_l$sample[9], PROC="ROUT",PROBE="Sheet 1",Q,par=res$par_final))
|
||||
# }
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user