diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/calfunctions.R | 2 | ||||
| -rw-r--r-- | R/calplot.R | 80 | ||||
| -rw-r--r-- | R/inverse.predict.lm.R | 115 | ||||
| -rw-r--r-- | R/lod.R | 53 | ||||
| -rw-r--r-- | R/loq.R | 40 |
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) -} |