
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, error = FALSE)
```

# Import of Data

We import the data and select the variables needed for the analysis.
```{r}
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 (!is.null(params$subset)) {
  calData <- subset(CalibrationData, eval(parse(text = params$subset)))
} else {
  calData <- CalibrationData
}
calData <- calData[,c(y.var, x.var)]
calData
```

# Model

We will apply the following model.
```{r}
params$formula
```

# Analysis

We now fit a generalized additive model.
```{r}
library(mgcv)
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.
```{r}
N <- 100
conc <- seq(from = min(calData[,params$concVar]), to = max(calData[,params$concVar]), len = N)
newDF <- data.frame(conc = conc)
names(newDF) <- params$concVar
resp <- as.vector(predict(fit, newdata = newDF))
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.
```{r}
library(ggplot2)
k <- params$k
modelPlot <- ggplot(calData, aes(x = .data[[params$concVar]],
                                 y = !!rlang::parse_expr(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}$.
```{r}
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.
```{r}
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.
```{r}
saveRDS(object = predFunc, file = paste0(params$filename, "_Model.rds"))
```

# Software

```{r}
sessionInfo()
```

