From 6d118690c0cae02fc5cd4b28c1a67eecde4d9f60 Mon Sep 17 00:00:00 2001 From: ranke Date: Thu, 11 May 2006 15:53:07 +0000 Subject: - The vignette is in a publisheable state - In addition to the Massart examples, the sample data from dintest (DIN 32645) has been tested - inverse.predict and calplot now also work on glm objects git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@7 5fad18fb-23f0-0310-ab10-e59a3bee62b4 --- R/inverse.predict.lm.R | 43 ++++++++++++++++++++++--------------------- 1 file changed, 22 insertions(+), 21 deletions(-) (limited to 'R') diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index f1921e4..f09dc3e 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -19,6 +19,27 @@ inverse.predict.default <- function(object, newdata, inverse.predict.lm <- function(object, newdata, ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), alpha=0.05, ss = "auto") +{ + if (length(object$weights) > 0) { + wx <- split(object$model$y,object$model$x) + w <- sapply(wx,mean) + } else { + w <- rep(1,length(split(object$model$y,object$model$x))) + } + .inverse.predict(object, newdata, ws, alpha, ss, w) +} + +inverse.predict.rlm <- function(object, newdata, + ws = mean(object$w), alpha=0.05, ss = "auto") +{ + wx <- split(object$weights,object$model$x) + w <- sapply(wx,mean) + .inverse.predict(object, newdata, ws, alpha, ss, w) +} + +.inverse.predict <- function(object, newdata, + ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), + alpha=0.05, ss = "auto", w) { if (length(object$coef) > 2) stop("More than one independent variable in your model - not implemented") @@ -33,12 +54,6 @@ inverse.predict.lm <- function(object, newdata, n <- length(yx) x <- as.numeric(names(yx)) ybar <- sapply(yx,mean) - if (length(object$weights) > 0) { - wx <- split(object$weights,object$model$x) - w <- sapply(wx,mean) - } else { - w <- rep(1,n) - } yhatx <- split(object$fitted.values,object$model$x) yhat <- sapply(yhatx,mean) se <- sqrt(sum(w*(ybar - yhat)^2)/(n-2)) @@ -52,21 +67,7 @@ inverse.predict.lm <- function(object, newdata, ybarw <- sum(w * ybar)/sum(w) -# The commented out code for sxhats is equation 8.28 without changes. It has -# been replaced by the code below, in order to be able to take into account a -# precision in the sample measurements that differs from the precision in the -# calibration samples. - -# sxhats <- se/b1 * sqrt( -# 1/(ws * m) + -# 1/sum(w) + -# (ybars - ybarw)^2 * sum(w) / -# (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2)) -# ) - -# This is equation 8.28, but with the possibility to take into account a -# different precision measurement of the sample and standard solutions -# in analogy to equation 8.26 +# This is an adapted form of equation 8.28 (see package vignette) sxhats <- 1/b1 * sqrt( ss^2/(ws * m) + se^2 * (1/sum(w) + -- cgit v1.2.1