diff options
Diffstat (limited to 'R/inverse.predict.lm.R')
| -rw-r--r-- | R/inverse.predict.lm.R | 35 |
1 files changed, 17 insertions, 18 deletions
diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index 8352c26..927e672 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -1,21 +1,21 @@ # 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 +# Example 8 on the same page, extended as specified in the package vignette inverse.predict <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { UseMethod("inverse.predict") } inverse.predict.default <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, ss = "auto") + 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, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { yname = names(object$model)[[1]] xname = names(object$model)[[2]] @@ -29,11 +29,11 @@ inverse.predict.lm <- function(object, newdata, ..., w <- rep(1,length(split(object$model[[yname]],object$model[[xname]]))) } .inverse.predict(object = object, newdata = newdata, - ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname) + 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, ss = "auto") + ws = "auto", alpha = 0.05, var.s = "auto") { yname = names(object$model)[[1]] xname = names(object$model)[[2]] @@ -43,10 +43,10 @@ inverse.predict.rlm <- function(object, newdata, ..., wx <- split(object$weights,object$model[[xname]]) w <- sapply(wx,mean) .inverse.predict(object = object, newdata = newdata, - ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname) + ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) } -.inverse.predict <- function(object, newdata, ws, alpha, ss, w, xname, 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") @@ -57,27 +57,26 @@ inverse.predict.rlm <- function(object, newdata, ..., ybars <- mean(newdata) m <- length(newdata) - yx <- split(object$model[[yname]],object$model[[xname]]) + 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 (ss == "auto") { - ss <- se - } else { - ss <- ss + 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 an adapted form of equation 8.28 (see package vignette) +# This is the adapted form of equation 8.28 (see package vignette) sxhats <- 1/b1 * sqrt( - ss^2/(ws * m) + + (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))) |