aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--R/calfunctions.R2
-rw-r--r--R/inverse.predict.lm.R20
-rw-r--r--inst/doc/chemCal.Rnw36
-rw-r--r--man/ipowfunc.Rd33
-rw-r--r--man/powfunc.Rd32
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}

Contact - Imprint