knitr::opts_chunk$set(echo = TRUE)
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
We will apply the following model.
params$formula
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
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
We save the results.
save(fit, LOB, LOD, LOQ, file = paste0(params$filename, "_Results.RData"))
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"))
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.