diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/calplot.R | 11 | ||||
| -rw-r--r-- | R/inverse.predict.lm.R | 22 | ||||
| -rw-r--r-- | R/lod.R | 23 | ||||
| -rw-r--r-- | R/loq.R | 2 |
4 files changed, 35 insertions, 23 deletions
diff --git a/R/calplot.R b/R/calplot.R index feb9727..2deed5a 100644 --- a/R/calplot.R +++ b/R/calplot.R @@ -21,9 +21,10 @@ calplot.lm <- function(object, xlim = "auto", ylim = "auto", m <- object level <- 1 - alpha - x <- m$model$x - y <- m$model$y - newdata <- data.frame(x = seq(0,max(x),length=250)) + y <- m$model[[1]] + x <- m$model[[2]] + newdata <- list(x = seq(0,max(x),length=250)) + names(newdata) <- names(m$model)[[2]] 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)) @@ -36,9 +37,9 @@ calplot.lm <- function(object, xlim = "auto", ylim = "auto", ylim = ylim ) points(x,y, pch = 21, bg = "yellow") - matlines(newdata$x, pred.lim, lty = c(1, 4, 4), + matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), col = c("black", "red", "red")) - matlines(newdata$x, conf.lim, lty = c(1, 3, 3), + matlines(newdata[[1]], conf.lim, lty = c(1, 3, 3), col = c("black", "green4", "green4")) legend(min(x), diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index 2ead9fb..e5f014c 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -17,32 +17,36 @@ inverse.predict.default <- function(object, newdata, ..., inverse.predict.lm <- function(object, newdata, ..., ws = "auto", alpha = 0.05, ss = "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$x) + wx <- split(object$weights,object$model[[xname]]) w <- sapply(wx,mean) } else { - w <- rep(1,length(split(object$model$y,object$model$x))) + 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) + ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname) } inverse.predict.rlm <- function(object, newdata, ..., ws = "auto", alpha = 0.05, ss = "auto") { + yname = names(object$model)[[1]] + xname = names(object$model)[[2]] if (ws == "auto") { ws <- mean(object$w) } - wx <- split(object$weights,object$model$x) + 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) + ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname) } -.inverse.predict <- function(object, newdata, ws, alpha, ss, w) +.inverse.predict <- function(object, newdata, ws, alpha, ss, w, xname, yname) { if (length(object$coef) > 2) stop("More than one independent variable in your model - not implemented") @@ -53,12 +57,12 @@ inverse.predict.rlm <- function(object, newdata, ..., ybars <- mean(newdata) m <- length(newdata) - yx <- split(object$model$y,object$model$x) + yx <- split(object$model[[yname]],object$model[[xname]]) n <- length(yx) df <- n - length(objects$coef) x <- as.numeric(names(yx)) ybar <- sapply(yx,mean) - yhatx <- split(object$fitted.values,object$model$x) + yhatx <- split(object$fitted.values,object$model[[xname]]) yhat <- sapply(yhatx,mean) se <- sqrt(sum(w*(ybar - yhat)^2)/df) if (ss == "auto") { @@ -67,7 +71,7 @@ inverse.predict.rlm <- function(object, newdata, ..., ss <- ss } - b1 <- object$coef[["x"]] + b1 <- object$coef[[xname]] ybarw <- sum(w * ybar)/sum(w) @@ -10,18 +10,25 @@ lod.default <- function(object, ..., alpha = 0.05, beta = 0.05) lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05) { - y0 <- predict(object, data.frame(x = 0), interval="prediction", - level = 1 - alpha) + xname <- names(object$model)[[2]] + newdata <- data.frame(0) + names(newdata) <- xname + y0 <- predict(object, newdata, interval="prediction", + level = 1 - 2 * alpha ) yc <- y0[[1,"upr"]] xc <- inverse.predict(object,yc)[["Prediction"]] f <- function(x) { - # Here I need the variance of y values as a function of x or - # y values - # Strangely, setting the confidence level to 0.5 does not result - # in a zero confidence or prediction interval + 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 } - lod.x <- optimize(f,interval=c(0,max(object$model$x)))$minimum - lod.y <- predict(object, data.frame(x = lod.x)) + lod.x <- optimize(f,interval=c(0,max(object$model[[xname]])))$minimum + newdata <- data.frame(x = lod.x) + names(lod.x) <- xname + lod.y <- predict(object, data.frame(lod.x)) return(list(x = lod.x, y = lod.y)) } @@ -5,7 +5,7 @@ loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") { - stop("lod is only implemented for univariate lm objects.") + stop("loq is only implemented for univariate lm objects.") } loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") |