knitr::opts_chunk$set(echo = TRUE)
We import the data and select the variables needed for the analysis.
load(paste0(PATH.OUT, "/CalibrationData.RData")) tmp <- strsplit(FORMULA, "~")[[1]] y.var <- names(unlist(sapply(colnames(CalibrationData), grep, x = tmp[1]))) x.vars <- 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.vars)] calData
We will apply the following model.
as.formula(FORMULA)
We now run a linear regression analysis including stepwise model selection by the Akaike information criterion (AIC).
In the fist step, we compute the full model.
library(MASS) full.fit <- lm(as.formula(FORMULA), data = calData) summary(full.fit)
We perform a stepwise model selection starting with the full model.
AIC.fit <- stepAIC(full.fit, trace = 0) AIC.fit AIC.fit.sum <- summary(AIC.fit) AIC.fit.sum
We plot the given concentrations against the fitted values.
library(ggplot2) R2 <- round(AIC.fit.sum$r.squared, 3) adj.R2 <- round(AIC.fit.sum$adj.r.squared, 3) DF <- data.frame(Observed = AIC.fit$model[,1], Fitted = fitted(AIC.fit)) ggplot(DF, aes(x = Observed, y = Fitted)) + geom_point() + geom_abline(slope = 1, intercept = 0) + xlab("Observed values") + ylab("Fitted values") + annotate("text", x=min(DF$Observed), y = max(DF$Fitted), label = substitute(paste(R^2, " = ", R2, ", adj. ", R^2, " = ", adj.R2), list(R2 = R2, adj.R2 = adj.R2)), vjust=1, hjust=0, size = 5)
We compute limit of blank (LOB), limit of detection (LOD) and limit of quantification (LOQ) using the linear regression fit. We get the LOB as the upper bound of the one-sided 95\% confidence interval of the intercept. Based on the LOB and the standard error of the residuals, we compute LOD and finally LOQ based on LOD.
LOB <- confint(AIC.fit, level = 0.90)[1,2] names(LOB) <- "LOB" LOB LOD <- LOB + 1.645*summary(AIC.fit)$sigma names(LOD) <- "LOD" LOD LOQ <- 3.3*LOD names(LOQ) <- "LOQ" LOQ
We save the results.
save(AIC.fit, LOB, LOD, LOQ, file = paste0(PATH.OUT, "/CalibrationResults.RData"))
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.