aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/calfunctions.R2
-rw-r--r--R/calplot.R80
-rw-r--r--R/inverse.predict.lm.R115
-rw-r--r--R/lod.R53
-rw-r--r--R/loq.R40
5 files changed, 0 insertions, 290 deletions
diff --git a/R/calfunctions.R b/R/calfunctions.R
deleted file mode 100644
index 6ce29f7..0000000
--- a/R/calfunctions.R
+++ /dev/null
@@ -1,2 +0,0 @@
-powfunc <- function(x,a,b) a * x^b
-ipowfunc <- function(y,a,b) (y/a)^1/b
diff --git a/R/calplot.R b/R/calplot.R
deleted file mode 100644
index 6aed9c0..0000000
--- a/R/calplot.R
+++ /dev/null
@@ -1,80 +0,0 @@
-calplot <- function(object,
- xlim = c("auto", "auto"), ylim = c("auto", "auto"),
- xlab = "Concentration", ylab = "Response", alpha = 0.05,
- varfunc = NULL)
-{
- UseMethod("calplot")
-}
-
-calplot.default <- function(object,
- xlim = c("auto","auto"), ylim = c("auto","auto"),
- xlab = "Concentration", ylab = "Response",
- alpha=0.05, varfunc = NULL)
-{
- stop("Calibration plots only implemented for univariate lm objects.")
-}
-
-calplot.lm <- function(object,
- xlim = c("auto","auto"), ylim = c("auto","auto"),
- xlab = "Concentration", ylab = "Response", alpha=0.05,
- varfunc = NULL)
-{
- if (length(object$coef) > 2)
- stop("More than one independent variable in your model - not implemented")
-
- if (alpha <= 0 | alpha >= 1)
- stop("Alpha should be between 0 and 1 (exclusive)")
-
- m <- object
- level <- 1 - alpha
- y <- m$model[[1]]
- x <- m$model[[2]]
- if (xlim[1] == "auto") xlim[1] <- 0
- if (xlim[2] == "auto") xlim[2] <- max(x)
- xlim <- as.numeric(xlim)
- newdata <- list(
- x = seq(from = xlim[[1]], to = xlim[[2]], length=250))
- names(newdata) <- names(m$model)[[2]]
- if (is.null(varfunc)) {
- varfunc <- if (length(m$weights)) {
- function(variable) mean(m$weights)
- } else function(variable) rep(1,250)
- }
- pred.lim <- predict(m, newdata, interval = "prediction",
- level=level, weights.newdata = varfunc(m))
- conf.lim <- predict(m, newdata, interval = "confidence",
- level=level)
- yrange.auto <- range(c(0,pred.lim))
- if (ylim[1] == "auto") ylim[1] <- yrange.auto[1]
- if (ylim[2] == "auto") ylim[2] <- yrange.auto[2]
- plot(1,
- type = "n",
- xlab = xlab,
- ylab = ylab,
- xlim = as.numeric(xlim),
- ylim = as.numeric(ylim)
- )
- 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",
- "Prediction Bands"),
- lty = c(1, 3, 4),
- 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
deleted file mode 100644
index d57275c..0000000
--- a/R/inverse.predict.lm.R
+++ /dev/null
@@ -1,115 +0,0 @@
-# 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, extended as specified in the package vignette
-
-inverse.predict <- function(object, newdata, ...,
- ws = "auto", alpha = 0.05, var.s = "auto")
-{
- UseMethod("inverse.predict")
-}
-
-inverse.predict.default <- function(object, newdata, ...,
- ws = "auto", alpha = 0.05, var.s = "auto")
-{
- stop("Inverse prediction only implemented for univariate lm and nls objects.")
-}
-
-inverse.predict.lm <- 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]])))
- }
- .inverse.predict(object = object, newdata = newdata,
- 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, var.s = "auto")
-{
- yname = names(object$model)[[1]]
- xname = names(object$model)[[2]]
- if (ws == "auto") {
- ws <- mean(object$w)
- }
- wx <- split(object$weights,object$model[[xname]])
- w <- sapply(wx,mean)
- .inverse.predict(object = object, newdata = 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)
- stop("More than one independent variable in your model - not implemented")
-
- if (alpha <= 0 | alpha >= 1)
- stop("Alpha should be between 0 and 1 (exclusive)")
-
- ybars <- mean(newdata)
- m <- length(newdata)
-
- 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 (var.s == "auto") {
- var.s <- se^2/ws
- }
-
- b1 <- object$coef[[xname]]
-
- ybarw <- sum(w * ybar)/sum(w)
-
-# This is the adapted form of equation 8.28 (see package vignette)
- sxhats <- 1/b1 * sqrt(
- (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)))
- )
-
- if (names(object$coef)[1] == "(Intercept)") {
- b0 <- object$coef[["(Intercept)"]]
- } else {
- b0 <- 0
- }
-
- xs <- (ybars - b0) / b1
- t <- qt(1-0.5*alpha, n - 2)
- conf <- t * sxhats
- result <- list("Prediction"=xs,"Standard Error"=sxhats,
- "Confidence"=conf, "Confidence Limits"=c(xs - conf, xs + conf))
- return(result)
-}
diff --git a/R/lod.R b/R/lod.R
deleted file mode 100644
index f5bbb7d..0000000
--- a/R/lod.R
+++ /dev/null
@@ -1,53 +0,0 @@
-lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default")
-{
- UseMethod("lod")
-}
-
-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, method = "default")
-{
- if (length(object$weights) > 0) {
- stop(paste(
- "\nThe detemination of a lod from calibration models obtained by",
- "weighted linear regression requires confidence intervals for",
- "predicted y values taking into account 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"
- ))
- }
- xname <- names(object$model)[[2]]
- yname <- names(object$model)[[1]]
- newdata <- data.frame(0)
- names(newdata) <- xname
- y0 <- predict(object, newdata, interval = "prediction",
- level = 1 - 2 * alpha)
- yc <- y0[[1,"upr"]]
- if (method == "din") {
- y0.d <- predict(object, newdata, interval = "prediction",
- level = 1 - 2 * beta)
- 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 <- list(lod.x, lod.y)
- names(lod) <- c(xname, yname)
- return(lod)
-}
diff --git a/R/loq.R b/R/loq.R
deleted file mode 100644
index 5776096..0000000
--- a/R/loq.R
+++ /dev/null
@@ -1,40 +0,0 @@
-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.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.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.loq,
- var.s = var.loq, alpha = alpha)
- (p[["Prediction"]] - k * p[["Confidence"]])^2
- }
- tmp <- optimize(f,interval=c(0,max(object$model[[2]])))
- loq.x <- tmp$minimum
- newdata <- data.frame(x = loq.x)
- names(newdata) <- xname
- loq.y <- predict(object, newdata)
- loq <- list(loq.x, loq.y)
- names(loq) <- c(xname, yname)
- return(loq)
-}

Contact - Imprint