diff options
| -rw-r--r-- | R/calfunctions.R | 2 | ||||
| -rw-r--r-- | R/inverse.predict.lm.R | 20 | ||||
| -rw-r--r-- | inst/doc/chemCal.Rnw | 36 | ||||
| -rw-r--r-- | man/ipowfunc.Rd | 33 | ||||
| -rw-r--r-- | man/powfunc.Rd | 32 |
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} |