probit的LOD计算

原理

具体的原理见下图

使用probit 进行拟合，得到预测某个位置的具体LOD数值

代码

suppressMessages(library(ggplot2))
suppressMessages(library(reshape))
suppressMessages(library(ggpubr))
suppressMessages(library(DescTools))



setwd("C:\\Users\\Desktop")
infile <- "Raw.txt"
outfile <- "Raw.pdf"
max_concentration <- 5
test <- "two.sided"
times <- 100
ci.level <- 0.95
method <- "wald"
lod_res <- "PanelM.T7.Raw.res.txt"
marker <- "Concentration (%)"
title <- "SNV DowmSample"

df <- read.table(infile, header= T)
if(any(c('percent', 'concentration') != colnames(df))){
  stop("Check first line of your input file.\n", call.=FALSE)
}

############################################### fit your curve #####################################
myprobit1 = glm(percent ~ concentration , df ,family = binomial(link ="probit"))

fit=predict(myprobit1 , data.frame(concentration  = seq(0 , max_concentration ,length.out = 10000)), type = "response")
data_s1 = data.frame(concentration = seq(0 , max_concentration ,length.out = 10000) , fit=fit)

data_s1$up = BinomCI(data_s1$fit * times , times , sides = test , conf.level = ci.level, method = method)[,3]
data_s1$low = BinomCI(data_s1$fit*times , times,sides = test,conf.level = ci.level, method = method)[,2]

######################################## calculate lod 90/95 #######################################
lod_95_low = round(BinomCI(0.95*times , times , sides = test , conf.level = ci.level, method = method)[1,2],3)
lod_95_up = round(BinomCI(0.95*times , times , sides = test , conf.level = ci.level, method = method)[1,3],3)

lod_95 = round(min(data_s1[data_s1$fit > 0.95,'concentration']),4)
vaf_95_low = round(min(data_s1[data_s1$up > 0.95,'concentration']),4)
vaf_95_up = round(min(data_s1[data_s1$low > 0.95,'concentration']),4)
lod_90_low = round(BinomCI(0.90*times,times,sides = test,conf.level = 0.9, method = method)[1,2],3)
lod_90_up = round(BinomCI(0.90*times,times,sides = test,conf.level = 0.9, method = method)[1,3],3)

lod_90 = round(min(data_s1[data_s1$fit > 0.9,'concentration']),4)
vaf_90_low = round(min(data_s1[data_s1$up > 0.9,'concentration']),4)
vaf_90_up = round(min(data_s1[data_s1$low > 0.9,'concentration']),4)

data_out2 = data.frame('type' = marker, 'lod90' = lod_90, 'lwr90' = vaf_90_low,
                       'upr90' = vaf_90_up, 'cilwr90' = lod_90_low, 'ciupr90' = lod_90_up,
                       'lod95' = lod_95, 'lwr95' = vaf_95_low, 'upr95' = vaf_95_up,
                       'cilwr95' = lod_95_low, 'ciupr95' = lod_95_up)
write.table(data_out2 , lod_res , sep = '\t' , quote = F,row.names = F)

########################################################## plot ####################################
mytheme = theme(plot.title = element_text(face = "bold",size = 18,color = "brown",hjust = 0.44445),
                axis.title = element_text(face = "bold.italic",size = 12,color = "brown"),
                axis.text = element_text(face = "bold",size = 10,color = "darkblue"),
                legend.position = 'none')

p1 <- ggplot()+ geom_line(data = data_s1 , aes(x = concentration,y = fit)) +
  geom_ribbon(data = data_s1 , aes(x = concentration , ymin = low, ymax = up), alpha = 0.2) +
  geom_label(aes(x = lod_95,y = 0.8), label = paste0(' LOD = ' , format(lod_95,digits = 3)) , color = 'forestgreen' , hjust = 0) +
  geom_point(data = df , aes(x = concentration,y = percent) , size = 1,color = 'forestgreen') +
  geom_hline(yintercept = 0.95, color = 'forestgreen' , linetype = 2) +
  geom_vline(xintercept = lod_95 ,color = 'forestgreen' , linetype = 2) +
  labs(x = marker , y = 'Detection probability',title = title) +
  scale_x_continuous(limits = c(0, max_concentration )) +
  theme_bw() +
  mytheme + scale_color_manual(values = c('forestgreen') , aesthetics = c("colour", "fill")) + guides(color = F)
p1
ggsave(filename = outfile , p1,height = 4, width = 4)