aboutsummaryrefslogtreecommitdiff
path: root/R/inverse.predict.lm.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/inverse.predict.lm.R')
-rw-r--r--R/inverse.predict.lm.R97
1 files changed, 97 insertions, 0 deletions
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)
+}

Contact - Imprint