Start of the integration of nonlinear calibration models

git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@19 5fad18fb-23f0-0310-ab10-e59a3bee62b4

5 files changed, 120 insertions, 3 deletions

diff --git a/R/calfunctions.R b/R/calfunctions.R
new file mode 100644
index 0000000..6ce29f7
--- /dev/null
+++ b/R/calfunctions.R
@@ -0,0 +1,2 @@
+powfunc <- function(x,a,b) a * x^b
+ipowfunc <- function(y,a,b) (y/a)^1/b
diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index 927e672..d57275c 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -11,7 +11,7 @@ inverse.predict <- function(object, newdata, ...,
 inverse.predict.default <- function(object, newdata, ...,
   ws = "auto", alpha = 0.05, var.s = "auto")
 {
-  stop("Inverse prediction only implemented for univariate lm objects.")
+  stop("Inverse prediction only implemented for univariate lm and nls objects.")
 }
 
 inverse.predict.lm <- function(object, newdata, ...,
@@ -46,6 +46,24 @@ inverse.predict.rlm <- function(object, newdata, ...,
     ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname)
 }
 
+inverse.predict.nls <- function(object, newdata, ...,
+  ws = "auto", alpha = 0.05, var.s = "auto")
+{
+  yname = names(object$model)[[1]]
+  xname = names(object$model)[[2]]
+  if (ws == "auto") {
+    ws <- ifelse(length(object$weights) > 0, mean(object$weights), 1)
+  }
+  if (length(object$weights) > 0) {
+    wx <- split(object$weights,object$model[[xname]])
+    w <- sapply(wx,mean)
+  } else {
+    w <- rep(1,length(split(object$model[[yname]],object$model[[xname]])))
+  }
+  if (length(object$coef) > 2)
+    stop("More than one independent variable in your model - not implemented")
+}
+
 .inverse.predict <- function(object, newdata, ws, alpha, var.s, w, xname, yname)
 {
   if (length(object$coef) > 2)
diff --git a/inst/doc/chemCal.Rnw b/inst/doc/chemCal.Rnw
index 8cdc97c..956c664 100644
--- a/inst/doc/chemCal.Rnw
+++ b/inst/doc/chemCal.Rnw
@@ -29,8 +29,8 @@ on this subject. However, I did not encounter any proof or explanation of the
 formula cited below yet, so I can't really confirm that Massart's method is correct.
 
 When calibrating an analytical method, the first task is to generate a suitable
-model. If we want to use the \texttt{chemCal} functions, we will have to restrict
-ourselves to univariate, possibly weighted, linear regression so far.
+model. If we want to use the \texttt{chemCal} functions, we will have to
+restrict ourselves to univariate, possibly weighted, linear regression so far.
 
 Once such a model has been created, the calibration can be graphically
 shown by using the \texttt{calplot} function:
@@ -119,6 +119,38 @@ where I interpret $\frac{{s_s}^2}{w_s}$ as an estimator of the variance at locat
 $\hat{x_s}$, which can be replaced by a user-specified value using the argument
 \texttt{var.s} of the function \texttt{inverse.predict}.
 
+\section*{Fitting and using a variance function}
+
+In the R package \texttt{nlme} variance functions are used for weighted 
+regressions. But since the \texttt{predict.nlme} method does not calculate
+prediction intervals, this is not useful for the \texttt{calplot} function.
+
+Two approaches could be used for fitting variance functions, one based on
+residuals from an unweighted fit, and one based on just the variances
+of the different samples along the x axis. If we used the residuals for
+fitting, a bias of the model in a certain area would result in a higher
+variance, so it seems preferable to choose the second approach. Of course,
+a prerequisite is to have sufficient repetitions for each sample in any
+case.
+
+Let's take the above example and estimate a variance function
+
+<<>>=
+massart97ex3
+massart97ex3$x <- factor(massart97ex3$x)
+summary <- summaryBy(y~x, data = massart97ex3,FUN=c(mean,sd,var))
+summary$x <- as.numeric(as.vector((summary$x)))
+plot(summary$x, summary$y.var,
+	xlim=c(0,50), 
+	ylim=c(0,max(summary$y.var)))
+varModel <- lm(y.var ~ I(x^2) + x, data=summary)
+varCurve <- predict(varModel, newdata=data.frame(x=0:5000/100))
+lines(0:5000/100,varCurve)
+
+
+
+
+
 \begin{thebibliography}{1}
 \bibitem{massart97}
 Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J.,
diff --git a/man/ipowfunc.Rd b/man/ipowfunc.Rd
new file mode 100644
index 0000000..c0f946b
--- /dev/null
+++ b/man/ipowfunc.Rd
@@ -0,0 +1,33 @@
+\name{ipowfunc}
+\alias{ipowfunc}
+\title{Power function}
+\description{
+  Inverse of the arithmetic power function \code{\link{powfunc}} used for
+  modelling univariate nonlinear calibration data. }
+\usage{
+  powfunc(x,a,b)
+}
+\arguments{
+  \item{x}{
+    Independent variable}
+  \item{a}{
+    Coefficient}
+  \item{b}{
+    Exponent}
+}
+\value{
+  The result of evaluating the function
+    \deqn{f(x) = \frac{y}{a}^\frac{1}{b}}{f(x) = y/a^1/b}
+  which is the inverse of the function defined by \code{\link{powfunc}}
+} 
+\author{
+  Johannes Ranke 
+  \email{jranke@uni-bremen.de} 
+  \url{http://www.uft.uni-bremen.de/chemie/ranke}
+}
+\seealso{
+  The original function \code{\link{powfunc}}.
+}
+\keyword{models}
+\keyword{regression}
+\keyword{nonlinear}
diff --git a/man/powfunc.Rd b/man/powfunc.Rd
new file mode 100644
index 0000000..73fe3b0
--- /dev/null
+++ b/man/powfunc.Rd
@@ -0,0 +1,32 @@
+\name{powfunc}
+\alias{powfunc}
+\title{Power function}
+\description{
+  Arithmetic power function for modelling univariate nonlinear calibration data.
+}
+\usage{
+  powfunc(x,a,b)
+}
+\arguments{
+  \item{x}{
+    Independent variable}
+  \item{a}{
+    Coefficient}
+  \item{b}{
+    Exponent}
+}
+\value{
+  The result of evaluating the function
+    \deqn{f(x) = a x^b}{f(x) = a * x^b}
+} 
+\author{
+  Johannes Ranke 
+  \email{jranke@uni-bremen.de} 
+  \url{http://www.uft.uni-bremen.de/chemie/ranke}
+}
+\seealso{
+  The inverse of this function \code{\link{ipowfunc}}.
+}
+\keyword{models}
+\keyword{regression}
+\keyword{nonlinear}