\documentclass[english, a4paper, 10pt, headings=small, DIV11]{scrartcl}
\SweaveInput{vignettes.defs}
\SweaveOpts{prefix.string=flu-fig-}
\hypersetup{pdftitle={hyperSpec vignette flu},
 pdfauthor={C. Beleites},
 pdfsubject={Vignette on hyperSpec's flu data set},
 pdfkeywords={hyperSpec, flu, calibration}}

% \VignetteIndexEntry{flu: Vignette on the fluorescence data set, how to write an import function, and example of a linear calibration}
% \VignetteKeywords{hyperSpec, flu, calibration}
% \VignettePackage{hyperSpec}

\begin{document}
\title{Calibration of Quinine Fluorescence Emission\\
Vignette for the Data Set \Robject{flu} of the R package \Rpackage{hyperSpec} }
\maketitle

\warnbox{Reproducing the Examples in this Vignette}{
The spectra files are shipped with \phy. This allows reproduction of the whole
vignette (source code and spectra files are in the package's documentation directory and its \texttt{rawdata} subdirectory).

For reproducing the examples in a live session, the full file names of the 
spectra can be found with the command:\\
\texttt{\small\textsl{> Sys.glob (paste (.libPaths (),
"hyperSpec/doc/rawdata/flu*.txt", sep = "/"))}}. 
Note that loading the package, and some definitions e.g. of the color palettes are executed in \texttt{vignettes.def}.
}

\tableofcontents

This vignette gives an example how to
\begin{itemize}
\item write an import function for a spectrometer manufacturer's proprietary ASCII files,
\item add further data columns to the spectra, and
\item set up a linear calibration (inverse least squares).
\end{itemize}

The data set \Robject{flu} in \Rpackage{hyperSpec} consists of 6 fluorescence emission spectra of quinine solutions.
They were acquired during an student practicum and were kindly provided by M. Kammer.

The concentrations of the solutions range from 0.05 mg/l to 0.30 mg/l. Spectra were acquired with a
Perkin Elmer LS50-B fluorescence spectrometer at 350 nm excitation.

\section{Writing an Import Function}
The raw spectra are in Perkin Elmer's ASCII file format, one spectrum per file. The files are completely ASCII text, with the actual spectra starting at line 55.

The function to import these files,  \Rfunction{scan.txt.PerkinElmer}, is discussed in the ``FileIO'' vignette, please refer to that document for details.

It  needs to be \Rfunction{source}d before use:
<<read.txt.PE>>=
source ("scan.txt.PerkinElmer.R")
flu <- scan.txt.PerkinElmer ("rawdata/flu?.txt", skip = 54)
@
Now the spectra are in a \Rclass{hyperSpec} object and can be examined e.g. by
<<rawspc>>=
flu
@
<<rawfig, fig = TRUE, width = 6, height = 3>>=
plot (flu)
@
\section{Adding further Data Columns}
The calibration model needs the quinine concentrations for the spectra.
This information can be stored together with the spectra, and also gets an appropriate label:
<<newdata>>=
flu$c <- seq (from = 0.05, to = 0.30, by = 0.05)
labels (flu, "c") <- "c / (mg / l)"
flu
save (flu, file = 'flu.rda')
@

Now the \Rclass{hyperSpec}  object \Robject{flu} contains two data columns, holding the actual spectra and the
respective concentrations. The dollar operator returns such a data column:
<<newc>>=
flu$c
@

\section{Dropping data columns}
\label{sec:dropp-data-columns}
\Rfunction{scan.txt.PerkinElmer} added a column with the file names that we don't need.
It is therefore deleted:
<<delcol>>=
flu$file <- NULL
@

\section{Linear Calibration}
As R is developed for the purpose of statistical analysis, tools for a least squares calibration
model are readily available.

The original spectra range from \Sexpr{min(wl(flu))} to \Sexpr{max(wl(flu))} nm.
However, the intensities at 450~nm are perfect for a univariate calibration. Plotting them over the
concentration is done by:

\setkeys{Gin}{width = .33\textwidth}
<<calplot1, fig=TRUE, width = 6, height = 6>>=
plotc (flu[,,450])
@

The square bracket operator extracts parts of a \Rclass{hyperSpec} object.
The first coordinate defines which spectra are to be used, the second which data columns, and the
third gives the spectral range.

We discard all the wavelengths but 450 nm:

<<cutspc>>=
flu <- flu [,,450]
labels (flu, "spc") <- expression (I ["450 nm"] / a.u.)
@

The plot could be enhanced by annotating the ordinate with the emission wavelength.
Also the axes should start at the origin, so that it is easier to see whether the calibration
function will go through the origin:
<<calplot2, fig=TRUE, width = 6, height = 6>>=
plotc (flu, xlim = range (0, flu$c), ylim = range (0, flu$spc))
@


The actual calibration is a linear model, which can be fitted by the R function \Rfunction{lm}.
\Rfunction{lm} needs a \emph{formula} that specifies which data columns are dependent and independent
variables. 

The normal calibration plot gives the emission intensity as a function of the concentration, and
the calibration function thus models $I = f (c)$, i. e.  $I = m c + b$ for a linear calibration.
This is then solved for $c$ when the calibration is used.

However, R's linear model is a quite strict in predicting: a model set up as $I = f (c)$ will predict
the intensity as a function of the concentration but not the other way round. Thus we set up an
inverse calibration model\footnote{As we can safely assume that the error on the concentrations is
far larger than the error on the instrument signal, it is actually the correct type of model from
the least squares fit point of view.}:  $c = f (I)$. The corresponding \R formula is \verb+c ~ I+, or in our case
\verb+c ~ spc+, as the intensities are stored in the data column \Rcode{\$spc}:

In addition, \Rfunction{lm} (like most R model building functions) expects the data to be a \Rclass{data.frame}.

There are three abbreviations that help to get the parts  of the \Rclass{hyperSpec} object that are frequently  needed:
\begin{description}
\setlength{\labelwidth}{2.5cm}
\setlength{\itemindent}{1.7cm}
 \item[\Rcode{flu[[]]}] returns the spectra matrix. It takes the same indices as \Rcode{[]}.
 \item[\Rcode{flu\$.}] returns the data as a \Rclass{data.frame}
 \item[\Rcode{flu\$..}] returns a  \Rclass{data.frame} that has all data columns but the spectra
\end{description}
<<abbrev>>=
flu[[]]
flu$.
flu$..
@

Putting this together, the calibration model is calculated:
<<cal>>=
calibration <- lm (c ~ spc, data = flu$.)
@

The \Rfunction{summary} gives a good overview of our model:
<<summarymodel>>=
summary (calibration)
@

In order to get predictions for new measurements, a new \Rclass{data.frame} with the same independent
variables (in columns with the same names) as in the calibration data are needed. Then the function
\Rfunction{predict} can be used. It can also calculate the prediction interval.
If we observe e.g. an intensity of 125 or 400 units, the corresponding concentrations and their 99 \% prediction intervals are:

<<pred>>=
I <- c (125, 400)
conc <- predict (calibration, newdata = list (spc = as.matrix(I)), interval = "prediction", 
                 level = .99)
conc
@

Finally, we can draw the calibration function and its 99 \% confidence interval
(also via \Rfunction{predict}) together with the prediction example. In order to draw the confidence interval into the calibration graph, we can either use a customized panel function:
<<calplot3, fig=TRUE, pdf = TRUE, eps = FALSE, width = 6, height = 6>>=

int <- list (spc = as.matrix(seq (min (flu), max(flu), length.out = 25)))
ci <- predict (calibration, newdata = int, interval = "confidence", level = 0.99)

panel.ci <-  function (x, y, ...,
                       intensity, ci.lwr, ci.upr, ci.col = "#606060") {
   panel.xyplot (x, y, ...)
   panel.lmline (x, y,...)
   panel.lines (ci.lwr, intensity, col = ci.col)
   panel.lines (ci.upr, intensity, col = ci.col)
}

plotc (flu, panel = panel.ci,
       intensity = int$spc, ci.lwr = ci [, 2], ci.upr = ci [, 3])
## # extrapolate to lower intensities
## int <- list (spc = as.matrix(0 : min (flu)))
## ci <- predict (calibration, newdata = int, interval = "confidence", level = 0.99)
## matlines (ci, int$spc, col = c ("red","#606060","#606060"), lty = 3)

# our example
#lines (conc[-1], rep(I, 2), col = "blue")
#points (conc[1], I, col = "blue", pch = 4, cex = 0.5)
@

Or, we can add the respective data to the \chy object. The meaning of the data can be saved in a new extra data column that acts as grouping variable for the plot.

First, the spectral range of \Robject{flu} is cut to contain the fluorescence emission at 450\,nm only, and the new column is introduced for the original data:
<<calplot4.1>>=
flu$type <- "data points"
@ 
Next, the calculated confidence intervals are appended:
<<calplot4.2>>=
tmp <- new ("hyperSpec", spc = as.matrix(seq (min (flu), max(flu), length.out = 25)),
                         wavelength = 450)
ci <-  predict (calibration, newdata = tmp$., interval = "confidence", level = 0.99)
tmp <- tmp [rep (seq (tmp, index = TRUE), 3)]
tmp$c <- as.numeric (ci)
tmp$type <- rep (colnames (ci), each = 25)

flu <- rbind (flu, tmp)
@

Finally, the resulting object is plotted. Our prediction example is handled by another customized panel function:
<<calplot4, fig=TRUE, pdf = TRUE, eps = FALSE, width = 6, height = 6>>=
panel.predict <- function (x, y, ..., 
                 intensity, ci, pred.col = "red", pred.pch = 19, pred.cex = 1) {
   panel.xyplot (x, y, ...)
   mapply (function (i, lwr, upr, ...) {
                 panel.lines (c (lwr, upr), rep (i, 2), ...)
              }, 
           intensity, ci [, 2], ci [, 3], MoreArgs = list (col = pred.col))
   panel.xyplot (ci [, 1], intensity, col = pred.col, pch = pred.pch, cex = pred.cex, type = "p")
}


plotc (flu, groups = type, type = c("l", "p"),
       col = c ("black", "black", "#606060", "#606060"), 
       pch = c (19, NA, NA, NA), cex = 0.5, 
       lty = c (0, 1, 1, 1),
       panel = panel.predict,
       intensity = I,
       ci = conc,
       pred.cex = 0.5)
@ 

\section*{Session Info}
\begin{footnotesize}
<<echo = FALSE>>=
<<cleanup>>
@ 
\end{footnotesize}
\end{document}


