# Copyright (C) 2005, 2008, 2010, 2012, 2014, 2015 Heather Turner
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 or 3 of the License
# (at your option).
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
predict.gnm <- function (object, newdata = NULL,
type = c("link", "response", "terms"), se.fit = FALSE,
dispersion = NULL, terms = NULL,
na.action = na.exclude, ...) {
type <- match.arg(type)
if (type == "terms") {
hasintercept <- attr(object$terms, "intercept") > 0L
## do not include eliminate term - cannot check estimability without
## creating full matrix, defeating point of eliminate
if (is.null(terms)) {
terms <- attr(object$terms, "term.labels")
} else {
terms <- setdiff(terms, "(eliminate)")
}
}
if (is.null(newdata)) {
pred <- switch(type, link = object$predictors,
response = object$fitted.values,
terms = {pred <- termPredictors(object)
## see 6.3.6 white book & predict.lm
if (hasintercept) {
predc <- sweep(pred, 2, colMeans(pred))
const <- sum(pred[1,]) - sum(predc[1,])
structure(predc[, terms, drop = FALSE],
constant = const)
} else structure(pred[, terms, drop = FALSE],
constant = 0)})
if (!is.null(na.act <- object$na.action)){
pred <- napredict(na.act, pred)
}
if (!inherits(pred, "matrix") && !is.null(object$table.attr))
attributes(pred) <- object$table.attr
} else {
modelTerms <- delete.response(terms(object))
if (is.null(object$eliminate)){
modelData <- model.frame(modelTerms, newdata, na.action = na.action,
xlev = object$xlevels)
} else {
## eliminate is evaluated in data/environment of formula
## => need to substitute here
modelData <-
eval(substitute(
model.frame(modelTerms, newdata,
eliminate = eliminate,
na.action = na.action, xlev = object$xlevels),
list(eliminate = object$call$eliminate)))
}
## use same contrasts as in original model
contr <- lapply(model.frame(object)[names(modelData)],
attr, "contrasts")
for (i in which(!vapply(contr, is.null, TRUE))){
modelData[[i]] <- C(modelData[[i]], contr[[i]])
}
if (length(offID <- attr(modelTerms, "offset"))){
offset <- eval(attr(modelTerms, "variables")[[offID + 1]],
newdata)
} else offset <- eval(object$call$offset, newdata)
modelTools <- gnmTools(modelTerms, modelData)
varPredictors <- modelTools$varPredictors(parameters(object))
pred <- modelTools$predictor(varPredictors, term = type == "terms")
if (type == "terms") {
rownames(pred) <- rownames(modelData)
} else names(pred) <- rownames(modelData)
if (!is.null(offset)) pred <- offset + pred
if (!is.null(object$eliminate)) {
prede <- attr(coef(object), "eliminate")
if (type != "terms") pred <- prede[modelData$`(eliminate)`] + pred
}
switch(type, response = {pred <- family(object)$linkinv(pred)},
terms = {if (hasintercept) {
predc <- sweep(pred, 2,
colMeans(termPredictors(object)))
const <- sum(pred[1,]) - sum(predc[1,])
pred <- structure(predc[, terms, drop = FALSE],
constant = const)
} else structure(pred[, terms, drop = FALSE],
constant = 0)},
link = )
if (!is.null(na.act <- attr(modelData, "na.action")))
pred <- napredict(na.act, pred)
}
if (se.fit) {
V <- vcov(object, dispersion = dispersion, with.eliminate = TRUE)
residual.scale <- as.vector(sqrt(attr(V, "dispersion")))
if (is.null(newdata)) {
X <- model.matrix(object)
elim <- object$eliminate
} else {
X <- modelTools$localDesignFunction(parameters(object),
varPredictors)
elim <- modelData$`(eliminate)`
}
covElim <- attr(V, "covElim")[elim, , drop = FALSE]
varElim <- attr(V, "varElim")[elim]
switch(type,
link = {
if (is.null(elim))
se.fit <- sqrt(diag(X %*% tcrossprod(V, X)))
else se.fit <-
sqrt(diag(X %*% tcrossprod(V, X)) +
2 * rowSums(X * covElim) + varElim)},
response = {
eta <- na.omit(c(family(object)$linkfun(pred)))
d <- family(object)$mu.eta(eta)
if (is.null(object$eliminate))
se.fit <- sqrt(diag(X %*% tcrossprod(V, X)))
else se.fit <- sqrt(diag(X %*% tcrossprod(V, X)) +
2*rowSums(X * covElim) + varElim)
se.fit <- d * se.fit},
terms = {
if (is.null(newdata)) {
assign <- split(seq(ncol(X)), attr(X, "assign"))
} else {
M <- model.matrix(object)
assign <- split(seq(ncol(X)), attr(M, "assign"))
}
if (hasintercept) {
if (is.null(newdata)) {
X <- sweep(X, 2, colMeans(X))
} else X <- sweep(X, 2, colMeans(M))
}
se.fit <- matrix(, nrow = nrow(X), ncol = length(terms))
s <- 0
adj <- hasintercept
for (i in match(terms, colnames(pred))) {
s <- s + 1
t <- assign[[i + adj]]
se.fit[, s] <-
sqrt(diag(X[, t] %*%
tcrossprod(V[t, t], X[, t])))
## check estimability of term
Xt <- X
Xt[, -t] <- 0
estimable <- checkEstimable(object, t(Xt))
is.na(se.fit)[estimable %in% c(FALSE, NA), s] <- TRUE
}
})
## check estimability of predictions
if (!is.null(newdata) && type != "terms"){
estimable <- checkEstimable(object, t(X))
is.na(se.fit)[estimable %in% c(FALSE, NA)] <- TRUE
}
if (!is.null(na.act)) {
se.fit <- napredict(na.act, se.fit)
}
if (inherits(pred, "table"))
attributes(se.fit) <- object$table.attr
else
attributes(se.fit) <- attributes(pred)
pred <- list(fit = pred, se.fit = se.fit,
residual.scale = residual.scale)
}
pred
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.