diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/calplot.R | 52 | ||||
| -rw-r--r-- | R/chemCal.R | 95 | ||||
| -rw-r--r-- | R/inverse.predict.lm.R | 95 |
3 files changed, 124 insertions, 118 deletions
diff --git a/R/calplot.R b/R/calplot.R new file mode 100644 index 0000000..cea1149 --- /dev/null +++ b/R/calplot.R @@ -0,0 +1,52 @@ +calplot <- function(object, xlim = "auto", ylim = "auto", + xlab = "Concentration", ylab = "Response", alpha=0.05) +{ + UseMethod("calplot") +} + +calplot.default <- function(object, xlim = "auto", ylim = "auto", + xlab = "Concentration", ylab = "Response", alpha=0.05) +{ + stop("Calibration plots only implemented for univariate lm objects.") +} + +calplot.lm <- function(object, xlim = "auto", ylim = "auto", + xlab = "Concentration", ylab = "Response", alpha=0.05) +{ + 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 + x <- m$model$x + y <- m$model$y + newdata <- data.frame(x = seq(0,max(x),length=250)) + pred.lim <- predict(m, newdata, interval = "prediction",level=level) + conf.lim <- predict(m, newdata, interval = "confidence",level=level) + if (xlim == "auto") xlim = c(0,max(x)) + if (ylim == "auto") ylim = range(c(pred.lim,y)) + plot(1, + type = "n", + xlab = xlab, + ylab = ylab, + xlim = xlim, + ylim = ylim + ) + points(x,y, pch = 21, bg = "yellow") + matlines(newdata$x, pred.lim, lty = c(1, 4, 4), + col = c("black", "red", "red")) + matlines(newdata$x, 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/chemCal.R b/R/chemCal.R deleted file mode 100644 index fab7db4..0000000 --- a/R/chemCal.R +++ /dev/null @@ -1,95 +0,0 @@ -calm <- function(data) -{ - y <- data[[2]] - x <- data[[1]] - m <- lm(y ~ x) - s <- summary(m) - if (s$coefficients[1,4] > 0.05) - { - m <- lm(y ~ x - 1) - s <- summary(m) - m$intercept <- FALSE - } else { - m$intercept <- TRUE - } - class(m) <- "calm" - m$yname <- names(data)[[1]] - m$xname <- names(data)[[2]] - return(m) -} -predict.calm <- predict.lm -print.calm <- print.lm -summary.calm <- summary.lm -plot.calm <- function(x,..., - xunit="",yunit="",measurand="", - level=0.95) -{ - m <- x - x <- m$model$x - y <- m$model$y - newdata <- data.frame(x = seq(0,max(x),length=250)) - pred.lim <- predict(m, newdata, interval = "prediction",level=level) - conf.lim <- predict(m, newdata, interval = "confidence",level=level) - if (xunit!="") { - xlab <- paste("Concentration in ",xunit) - } else xlab <- m$xname - if (yunit=="") yunit <- m$yname - if (measurand!="") { - main <- paste("Calibration for",measurand) - } else main <- "Calibration" - plot(1, - xlab = xlab, - ylab = yunit, - type = "n", - main = main, - xlim = c(0,max(x)), - ylim = range(pred.lim) - ) - points(x,y, pch = 21, bg = "yellow") - matlines(newdata$x, pred.lim, lty = c(1, 4, 4), - col = c("black", "red", "red")) - matlines(newdata$x, 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") -} -predictx <- function(m,yobs,level=0.95) -{ - s <- summary(m) - xi <- m$model$x - yi <- m$model$y - n <- length(yi) - p <- length(yobs) - if (p > 1) { - varyobs <- var(yobs) - } else { - varyobs <- 0 - } - if (!m$intercept) { - b1 <- summary(m)$coef["x","Estimate"] - varb1 <- summary(m)$coef["x","Std. Error"]^2 - xpred <- mean(yobs)/b1 - varxpred <- (varyobs + xpred^2 * varb1) / b1^2 - sdxpred <- sqrt(varxpred) - } else - { - b0 <- summary(m)$coef["(Intercept)","Estimate"] - b1 <- summary(m)$coef["x","Estimate"] - S <- summary(m)$sigma - xpred <- (mean(yobs) - b0)/b1 - sumxxbar <- sum((xi - mean(xi))^2) - yybar <- (mean(yobs) - mean(yi))^2 - sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5 - } - t <- qt((1 + level)/2,n - 2) - confxpred <- t * sdxpred - - result <- c(estimate=xpred,sdxpred=sdxpred,confxpred=confxpred) -} diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index f438d99..f1921e4 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -1,40 +1,89 @@ # This is an implementation of Equation (8.28) in the Handbook of Chemometrics -# and Qualimetrics, Part A, Massart et al, page 200, validated with Example 8 -# on the same page +# and Qualimetrics, Part A, Massart et al (1997), page 200, validated with +# Example 8 on the same page -inverse.predict <- function(object, newdata, alpha=0.05) +inverse.predict <- function(object, newdata, + ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), + alpha=0.05, ss = "auto") { UseMethod("inverse.predict") } -inverse.predict.default <- function(object, newdata, alpha=0.05) +inverse.predict.default <- function(object, newdata, + ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), + alpha=0.05, ss = "auto") { stop("Inverse prediction only implemented for univariate lm objects.") } -inverse.predict.lm <- function(object, newdata, alpha=0.05) +inverse.predict.lm <- function(object, newdata, + ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), + alpha=0.05, ss = "auto") { - if (is.list(newdata)) { - if (!is.null(newdata$y)) - newdata <- newdata$y - else - stop("Newdata list should contain element newdata$y") - } + 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) - if (is.matrix(newdata)) { - Y <- newdata - Ybar <- apply(Y,1,mean) - nrepl <- ncol(Y) + yx <- split(object$model$y,object$model$x) + 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) } - else { - Y <- as.vector(newdata) - Ybar <- Y - nrepl <- 1 + yhatx <- split(object$fitted.values,object$model$x) + yhat <- sapply(yhatx,mean) + se <- sqrt(sum(w*(ybar - yhat)^2)/(n-2)) + if (ss == "auto") { + ss <- se + } else { + ss <- ss } - if (length(object$coef) > 2) - stop("Inverse prediction not yet implemented for more than one independent variable") + b1 <- object$coef[["x"]] - if (alpha <= 0 | alpha >= 1) - stop("Alpha should be between 0 and 1 (exclusive)") + 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 + sxhats <- 1/b1 * sqrt( + ss^2/(ws * 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) } |