knitr::opts_chunk$set(echo = TRUE)

Import of Data

We import the data and select the variables needed for the analysis.

load(paste0(params$filename, "_Data.RData"))
tmp <- strsplit(params$formula, " ~ ")[[1]]

y.var <- names(unlist(sapply(colnames(CalibrationData), 
                             grep, x = tmp[1])))
x.var <- names(unlist(sapply(colnames(CalibrationData), 
                              grep, x = tmp[2])))
if(SUBSET != ""){
  calData <- eval(call("subset", x =  CalibrationData,
                       subset = parse(text = SUBSET)))
}else{
  calData <- CalibrationData
}
calData <- calData[,c(y.var, x.var)]
calData

Model

We will apply the following model.

params$formula

Analysis

We now fit a generalized additive model.

fit <- gam(as.formula(params$formula), data = calData, method = "REML")
summary(fit)

We compute the predicted concentration values and determine the inverse of the fitted model by spline interpolation.

N <- 100
ab <- coef(fit)
names(ab) <- NULL
conc <- seq(from = 0, to = max(calData[,tmp[2]]), len = N)
newDF <- data.frame(conc = conc)
names(newDF) <- tmp[2]
resp <- predict(fit, newdata = newDF)
names(resp) <- NULL
predFunc <- function(newdata){}
body(predFunc) <- substitute({ Y <- with(newdata, eval(y))
                      fun <- splinefun(x = resp, y = conc, method = "monoH.FC")
                      fun(Y) }, list(y = parse(text = tmp[1]),
                                     resp = resp,
                                     conc = conc))
y0 <- predict(fit, newdata = newDF[1,,drop = FALSE], se = TRUE)
fun <- splinefun(x = resp, y = conc, method = "monoH.FC")

We plot the given concentrations against the fitted values.

library(ggplot2)
modelPlot <- ggplot(calData, aes_string(x = tmp[2], y = tmp[1])) +
                geom_point() + geom_smooth(method = "gam", formula = y ~ s(x, k = k))
modelPlot

Computation of LOB, LOD and LOQ

We compute limit of blank (LOB), limit of detection (LOD) and limit of quantification (LOQ) by inverting the regression fit. We get the LOB by inverting the upper bound of the one-sided 95\% confidence interval at concentration $0$. In case of LOD, the upper-bound of the 99.95\% confidence interval at concentration $0$ is inverted. LOQ is determined as $3\times\textrm{LOD}$.

if(resp[1] < resp[length(resp)]){
  LOB <- fun(resp[1] + 1.645*y0$se.fit)
  LOD <- fun(resp[1] + 3.3*y0$se.fit)
}else{
  LOB <- fun(resp[1] - 1.645*y0$se.fit)
  LOD <- fun(resp[1] - 3.3*y0$se.fit)
}
names(LOB) <- "LOB"
names(LOD) <- "LOD"
LOQ <- 3*LOD
names(LOQ) <- "LOQ"
LOB
LOD
LOQ

Save Results

We save the results.

save(fit, LOB, LOD, LOQ, file = paste0(params$filename, "_Results.RData"))

Save Model

We save the inverse of the fitted model to be able to apply it for predicting concentrations.

saveRDS(object = predFunc, file = paste0(params$filename, "_Model.rds"))

Software

sessionInfo()


Try the LFApp package in your browser

Any scripts or data that you put into this service are public.

LFApp documentation built on Nov. 6, 2023, 5:08 p.m.