From 9e0dae397df8c18e7333d2604cae96b2a7927420 Mon Sep 17 00:00:00 2001 From: ranke Date: Fri, 23 Jun 2006 15:33:27 +0000 Subject: - inverse.predict now has a var.s argument instead of the never tested ss argument. This is documented in the updated vignette - loq() now has w.loq and var.loq arguments, and stops with a message if neither are specified and the model has weights. - calplot doesn't stop any more for weighted regression models, but only refrains from drawing prediction bands - Added method = "din" to lod(), now that I actually have it (DIN 32645) and was able to see which approximation is used therein. - removed the demos, as the examples and tests are already partially duplicated - The vignette is more of a collection of various notes, but should certainly be helpful for the user. - Version bump to 0.1-xxx git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@16 5fad18fb-23f0-0310-ab10-e59a3bee62b4 --- R/calplot.R | 23 ++++++++++++----------- R/inverse.predict.lm.R | 35 +++++++++++++++++------------------ R/lod.R | 40 +++++++++++++++++++++++----------------- R/loq.R | 18 ++++++++++++++---- 4 files changed, 66 insertions(+), 50 deletions(-) (limited to 'R') diff --git a/R/calplot.R b/R/calplot.R index 0cc6120..6aed9c0 100644 --- a/R/calplot.R +++ b/R/calplot.R @@ -1,6 +1,6 @@ calplot <- function(object, - xlim = c("auto","auto"), ylim = c("auto","auto"), - xlab = "Concentration", ylab = "Response", alpha=0.05, + xlim = c("auto", "auto"), ylim = c("auto", "auto"), + xlab = "Concentration", ylab = "Response", alpha = 0.05, varfunc = NULL) { UseMethod("calplot") @@ -22,14 +22,6 @@ calplot.lm <- function(object, if (length(object$coef) > 2) stop("More than one independent variable in your model - not implemented") - if (length(object$weights) > 0) { - stop(paste( - "\nConfidence and prediction intervals for weighted linear models require", - "weights for the x values from which the predictions are to be generated.", - "This is not supported by the internally used predict.lm method.", - sep = "\n")) - } - if (alpha <= 0 | alpha >= 1) stop("Alpha should be between 0 and 1 (exclusive)") @@ -65,9 +57,17 @@ calplot.lm <- function(object, points(x,y, pch = 21, bg = "yellow") matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), col = c("black", "red", "red")) + if (length(object$weights) > 0) { + legend(min(x), + max(pred.lim, na.rm = TRUE), + legend = c("Fitted Line", "Confidence Bands"), + lty = c(1, 3), + lwd = 2, + col = c("black", "green4"), + horiz = FALSE, cex = 0.9, bg = "gray95") + } else { matlines(newdata[[1]], conf.lim, lty = c(1, 3, 3), col = c("black", "green4", "green4")) - legend(min(x), max(pred.lim, na.rm = TRUE), legend = c("Fitted Line", "Confidence Bands", @@ -76,4 +76,5 @@ calplot.lm <- function(object, lwd = 2, col = c("black", "green4", "red"), horiz = FALSE, cex = 0.9, bg = "gray95") + } } diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index 8352c26..927e672 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -1,21 +1,21 @@ # This is an implementation of Equation (8.28) in the Handbook of Chemometrics # and Qualimetrics, Part A, Massart et al (1997), page 200, validated with -# Example 8 on the same page +# Example 8 on the same page, extended as specified in the package vignette inverse.predict <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { UseMethod("inverse.predict") } inverse.predict.default <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { stop("Inverse prediction only implemented for univariate lm objects.") } inverse.predict.lm <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { yname = names(object$model)[[1]] xname = names(object$model)[[2]] @@ -29,11 +29,11 @@ inverse.predict.lm <- function(object, newdata, ..., w <- rep(1,length(split(object$model[[yname]],object$model[[xname]]))) } .inverse.predict(object = object, newdata = newdata, - ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname) + ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) } inverse.predict.rlm <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { yname = names(object$model)[[1]] xname = names(object$model)[[2]] @@ -43,10 +43,10 @@ inverse.predict.rlm <- function(object, newdata, ..., wx <- split(object$weights,object$model[[xname]]) w <- sapply(wx,mean) .inverse.predict(object = object, newdata = newdata, - ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname) + ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) } -.inverse.predict <- function(object, newdata, ws, alpha, ss, w, xname, yname) +.inverse.predict <- function(object, newdata, ws, alpha, var.s, w, xname, yname) { if (length(object$coef) > 2) stop("More than one independent variable in your model - not implemented") @@ -57,27 +57,26 @@ inverse.predict.rlm <- function(object, newdata, ..., ybars <- mean(newdata) m <- length(newdata) - yx <- split(object$model[[yname]],object$model[[xname]]) + yx <- split(object$model[[yname]], object$model[[xname]]) n <- length(yx) df <- n - length(object$coef) x <- as.numeric(names(yx)) ybar <- sapply(yx,mean) - yhatx <- split(object$fitted.values,object$model[[xname]]) - yhat <- sapply(yhatx,mean) - se <- sqrt(sum(w*(ybar - yhat)^2)/df) - if (ss == "auto") { - ss <- se - } else { - ss <- ss + yhatx <- split(object$fitted.values, object$model[[xname]]) + yhat <- sapply(yhatx, mean) + se <- sqrt(sum(w * (ybar - yhat)^2)/df) + + if (var.s == "auto") { + var.s <- se^2/ws } b1 <- object$coef[[xname]] ybarw <- sum(w * ybar)/sum(w) -# This is an adapted form of equation 8.28 (see package vignette) +# This is the adapted form of equation 8.28 (see package vignette) sxhats <- 1/b1 * sqrt( - ss^2/(ws * m) + + (var.s / m) + se^2 * (1/sum(w) + (ybars - ybarw)^2 * sum(w) / (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2))) diff --git a/R/lod.R b/R/lod.R index 54618c8..f5bbb7d 100644 --- a/R/lod.R +++ b/R/lod.R @@ -1,14 +1,14 @@ -lod <- function(object, ..., alpha = 0.05, beta = 0.05) +lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default") { UseMethod("lod") } -lod.default <- function(object, ..., alpha = 0.05, beta = 0.05) +lod.default <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default") { stop("lod is only implemented for univariate lm objects.") } -lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05) +lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default") { if (length(object$weights) > 0) { stop(paste( @@ -24,23 +24,29 @@ lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05) yname <- names(object$model)[[1]] newdata <- data.frame(0) names(newdata) <- xname - y0 <- predict(object, newdata, interval="prediction", - level = 1 - 2 * alpha ) + y0 <- predict(object, newdata, interval = "prediction", + level = 1 - 2 * alpha) yc <- y0[[1,"upr"]] - xc <- inverse.predict(object,yc)[["Prediction"]] - f <- function(x) - { - newdata <- data.frame(x) - names(newdata) <- xname - pi.y <- predict(object, newdata, interval = "prediction", + if (method == "din") { + y0.d <- predict(object, newdata, interval = "prediction", level = 1 - 2 * beta) - yd <- pi.y[[1,"lwr"]] - (yd - yc)^2 + deltay <- y0.d[[1, "upr"]] - y0.d[[1, "fit"]] + lod.y <- yc + deltay + lod.x <- inverse.predict(object, lod.y)$Prediction + } else { + f <- function(x) { + newdata <- data.frame(x) + names(newdata) <- xname + pi.y <- predict(object, newdata, interval = "prediction", + level = 1 - 2 * beta) + yd <- pi.y[[1,"lwr"]] + (yd - yc)^2 + } + lod.x <- optimize(f,interval=c(0,max(object$model[[xname]])))$minimum + newdata <- data.frame(x = lod.x) + names(newdata) <- xname + lod.y <- predict(object, newdata) } - lod.x <- optimize(f,interval=c(0,max(object$model[[xname]])))$minimum - newdata <- data.frame(x = lod.x) - names(newdata) <- xname - lod.y <- predict(object, newdata) lod <- list(lod.x, lod.y) names(lod) <- c(xname, yname) return(lod) diff --git a/R/loq.R b/R/loq.R index ee22d38..5776096 100644 --- a/R/loq.R +++ b/R/loq.R @@ -1,22 +1,32 @@ -loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") +loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto") { UseMethod("loq") } -loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") +loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto") { stop("loq is only implemented for univariate lm objects.") } -loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") +loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto") { + if (length(object$weights) > 0 && var.loq == "auto" && w.loq == "auto") { + stop(paste("If you are using a model from weighted regression,", + "you need to specify a reasonable approximation for the", + "weight (w.loq) or the variance (var.loq) at the", + "limit of quantification")) + } xname <- names(object$model)[[2]] yname <- names(object$model)[[1]] f <- function(x) { newdata <- data.frame(x = x) names(newdata) <- xname y <- predict(object, newdata) - p <- inverse.predict(object, rep(y, n), ws = w, alpha = alpha) + p <- inverse.predict(object, rep(y, n), ws = w.loq, + var.s = var.loq, alpha = alpha) (p[["Prediction"]] - k * p[["Confidence"]])^2 } tmp <- optimize(f,interval=c(0,max(object$model[[2]]))) -- cgit v1.2.1