Thursday, February 11, 2016

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
}

3 comments:

  1. Hi,
    Wow I am so happy to have found your R-Code. It works like charm. I just have one question, is it possible get also the numbers of the confidence intervals that are shown in your plot.

    ReplyDelete
  2. Hi

    These are exactly the R-codes that I was looking to create the LOD plot, but I can't generate an output when I run these codes in R. Not sure if I am doing something wrong here.

    Your help is greatly appreciated.

    Ahmad

    ahmadr215@tpg.com.au

    ReplyDelete