lmpredint <- function(object, newdata, k, level = 0.95, alternative="two.sided", quantile=NULL, absError=0.001, interval=c(0, 100)){
tt <- terms(object)
if (!inherits(object, "lm")) warning("calling predict.lm(<fake-lm-object>) ...")
if (missing(newdata) || is.null(newdata)) {
stop("newdata is missing!")
} else {
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.pass, xlev = object$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
offset <- rep(0, nrow(X))
if (!is.null(off.num <- attr(tt, "offset"))) for (i in off.num) offset <- offset + eval(attr(tt, "variables")[[i + 1]], newdata)
if (!is.null(object$call$offset)) offset <- offset + eval(object$call$offset, newdata)
mmDone <- FALSE
}
if (k > nrow(newdata) | k < 1) stop("k needs to be between 0 < k <= nrow(newdata)")
n <- length(object$residuals)
p <- object$rank
p1 <- seq_len(p)
piv <- if (p) eval(parse(text="stats:::qr.lm(object)$pivot[p1]")) # ::: not allowed anymore
if (p < ncol(X) && !(missing(newdata) || is.null(newdata))) warning("prediction from a rank-deficient fit may be misleading")
beta <- object$coefficients
predictor <- drop(X[, piv, drop = FALSE] %*% beta[piv])
if (!is.null(offset)) predictor <- predictor + offset
w <- eval(parse(text="stats:::weights.default(object)")) # ::: not allowed anymore
if (!is.null(w)) {
weights <- w
warning("assuming prediction variance inversely proportional to weights used for fitting\n")
}
weights <- 1
if (!is.null(object$weights)) warning("Assuming constant prediction variance even though model fit is weighted\n")
if (inherits(weights, "formula")) {
if (length(weights) != 2L) stop("'weights' as formula should be one-sided")
d <- newdata
weights <- eval(weights[[2L]], d, environment(weights))
}
r <- object$residuals
w <- object$weights
rss <- sum(if (is.null(w)) r^2 else r^2 * w)
df <- object$df.residual
res.var <- rss/df
if (p > 0) {
XRinv <- X[, piv] %*% qr.solve(qr.R(eval(parse(text="stats:::qr.lm(object)")))[p1, p1])
ip <- drop(XRinv^2 %*% rep(res.var, p))
} else ip <- rep(0, n)
if (is.null(quantile)){
if (alternative == "two.sided"){
tfrac <- -PIcritval(k,m=nrow(newdata),n,level, df, absError=absError, interval=interval)
} else {
tfrac <- -PIonesided(k,m=nrow(newdata),n,level, df, absError=absError, interval=interval)
}
} else {
tfrac <- -quantile
}
pred.var = res.var/weights
hwid <- tfrac * sqrt(ip + pred.var)
predictor <- data.frame(cbind(predictor, predictor + hwid %o% c(1, -1)))
colnames(predictor) <- c("fit", "lower", "upper")
if (alternative == "less") predictor[,2] <- -Inf
if (alternative == "greater") predictor[,3] <- Inf
new(Class="PIlm",
prediction = predictor,
newdata = newdata,
model = object,
k = k,
level = level,
alternative = alternative)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.