Probit Model for Limit of Detection Estimation per CLSI EP17A2E Guidance with Application to Quantitative Molecular Measurement Procedures

Following is the R function:

lod.model <- function(dat)
{
    ## PURPOSE: Make a Probit Model Fitting and Ploting for LOD Data
    ## ----------------------------------------------------------------------
    ## ARGUMENT:
    ##   dat: a data frame with three columns "conc", "success", "failure"
    ##        "conc" is the Target Concentration Levels.
    ##        "success" is the Number of Positive Results.
    ##        "failure" is the Number of Negative Results. 
    ## ----------------------------------------------------------------------
    ## RETURN: Plot, GLM model object, and 95% LOD Estimate. 
    ## ----------------------------------------------------------------------
    ## AUTHOR: Feiming Chen,  Date: 10 Feb 2016, 14:28

    ## Per "CLSI EP17A2E" Guidance, apply default recommendation:
    ##   1. Remove Zero Predictor Values.
    ##   2. Appply "log10" transformation on predictors.
    ##   3. Use Probit Model
    ##   4. Report 95% LOD (concentration level)

    ## remove 0 levels in "conc" since will do "log10" transform on predictor. 
    dat <- subset(dat, conc > 0)

    r <- glm(cbind(success, failure) ~ log10(conc), data=dat, family=binomial(link="probit"))

    grid <- seq(0.001, max(1, dat$conc), length.out=200)
    preds <- predict(r, newdata=data.frame(conc=grid), type="link", se.fit=TRUE)

    se <- 1.96 * preds$se.fit
    upr <- unlist(r$family$linkinv(preds$fit + se))
    lwr <- unlist(r$family$linkinv(preds$fit - se))
    fit <- unlist(r$family$linkinv(preds$fit))

    ce <- coef(r)
    lod95 <- 10^((r$family$linkfun(0.95) - ce[1])/ce[2]) # 95% LOD Estimate
    
    ## Plot
    with(dat, plot(conc, success/(success+failure),ylab="Probability", xlab="Concentration", log="x", pch=19, ylim=c(0, 1)))
    title(main=paste("Probit Model, 95% LOD =", round(lod95, 3)))

    lines(grid, fit)
    lines(grid, lwr, lty=2)
    lines(grid, upr, lty=2)
    abline(h=0.95, lty=2, col="red")
    abline(v=lod95, lty=2, col="red")
    axis(2, 0.95, labels = "0.95")

    invisible(list(r=r, lod95=lod95))
}
if (F) {                                # Unit Test 
    dat <- data.frame(conc=c(5.4e5, 1.1e5, 5e4, 1e4, 5e3, 1e3), 
                      success=c(6, 6, 6, 3, 0, 0), 
                      failure=c(0, 0, 0, 3, 6, 6))
    lod.model(dat)

    ## CLSI Example: EP17A2E Table C1 & C2, Lot 1, LOD=0.077 CFU/mL. 
    dat <- data.frame(conc=c(0, 0.025, 0.05, 0.15, 0.3, 0.5, 0.006, 0.014),
                      success=c(0, 23, 29, 32, 32, 32, 11, 15), 
                      failure=c(22, 9, 3, 0, 0, 0, 19, 15))
    r <- lod.model(dat)                 # Verify 95% LOD = 0.077
}