# R code to accompany FTIR calibration model video
library(investr)
library(dplyr)
library(tidyr)
library(openxlsx)

# Read data from Excel file (download and save in "data" subfolder)
obs <- read.xlsx("data/APP3-A.xlsx", sheet=1)
standard <- read.xlsx("data/APP3-A.xlsx", sheet=2)

#### Inverse regression and calibration analysis
temp <- obs %>% pivot_longer(cols = c("A", "B", "C", "D", "E"), 
                             names_to="speciment", values_to="carbon") %>% 
  group_by(Day, speciment) %>% 
  summarize(carbon_avg=mean(carbon)) %>% 
  left_join(standard, by="speciment")

dat <- data.frame(y=temp$carbon_avg, x=temp$standard)

#### Initially fit a straight line calibration curve which shows some deficiencies
mod <- lm(y ~ x, data = dat)
summary(mod)
res <- invest(mod, y0 = 1, mean.response = TRUE)
plotFit(mod, interval = "prediction", shade = TRUE, col.pred = "darkolivegreen1",
        extend.range = TRUE)
abline(h = 1, v = c(res$lower, res$estimate, res$upper), lty = 2)


# Find the lower detection limit, p.28, the minimum amount of carbon can be detected
invest(mod, y0 = 0, mean.response = TRUE, lower=-0.01)

#### Alternatively fit a polynomial calibration curve, to accomodate lower D measurements
# Use invest() rather than calibrate() which relies on analytic solutions, not always work
mod <- lm(y ~ poly(x, 3), data = dat)
summary(mod)
res <- invest(mod, y0 = 1, mean.response = TRUE)
res
plotFit(mod, interval = "prediction", shade = TRUE, col.pred = "darkolivegreen1",
        extend.range = TRUE, xlab="碳样品标准量 (ppma)", ylab="仪器读取量 (ppma)")
abline(h = 1, v = c(res$lower, res$estimate, res$upper), lty = 2)


# Here is the direct way to find x based of the fitted polynomial y=1, ?not exactly as above
uniroot(function(x) predict(mod, newdata = list("x" = x))-1,
        lower = min(dat$x), upper = max(dat$x))$root

# Zoom-in plot
plotFit(mod, interval = "prediction", shade = TRUE, col.pred = "darkolivegreen1",
        extend.range = TRUE, xlim=c(1.1, 1.17), ylim=c(0.9, 1.1),
        xlab="碳样品标准量 (ppma)", ylab="仪器读取量 (ppma)")
abline(h = 1, v = c(res$lower, res$estimate, res$upper), lty = 2)

# Find the lower detection limit, p.28, the minimum amount of carbon can be detected
invest(mod, y0 = 0, mean.response = TRUE, lower=-0.01)