diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/calplot.R | 80 | ||||
| -rw-r--r-- | R/inverse.predict.lm.R | 97 | ||||
| -rw-r--r-- | R/lod.R | 55 | ||||
| -rw-r--r-- | R/loq.R | 41 |
4 files changed, 273 insertions, 0 deletions
diff --git a/R/calplot.R b/R/calplot.R new file mode 100644 index 0000000..6aed9c0 --- /dev/null +++ b/R/calplot.R @@ -0,0 +1,80 @@ +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 new file mode 100644 index 0000000..927e672 --- /dev/null +++ b/R/inverse.predict.lm.R @@ -0,0 +1,97 @@ +# 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 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 <- 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) +} @@ -0,0 +1,55 @@ +lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") +{ + UseMethod("lod") +} + +lod.default <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") +{ + stop("lod is only implemented for univariate lm objects.") +} + +lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "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]] + xvalues <- 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 + } + if (tol == "default") tol = min(xvalues[xvalues !=0]) / 1000 + lod.x <- optimize(f, interval = c(0, max(xvalues) * 10), tol = tol)$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) +} @@ -0,0 +1,41 @@ +loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto", tol = "default") +{ + UseMethod("loq") +} + +loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto", tol = "default") +{ + 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", tol = "default") +{ + 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]] + xvalues <- 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 + } + if (tol == "default") tol = min(xvalues[xvalues !=0]) / 1000 + loq.x <- optimize(f, interval = c(0, max(xvalues) * 10), tol = tol)$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) +} |