aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/calplot.R52
-rw-r--r--R/chemCal.R95
-rw-r--r--R/inverse.predict.lm.R95
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)
}

Contact - Imprint