MuMIn (version 1.3.6)

predict.averaging: Predict Method for the Averaged Model

Description

Model-averaged predictions with optional standard errors.

Usage

## S3 method for class 'averaging':
predict(object, newdata, se.fit = NULL, interval = NULL,
    type=c("link", "response"), ...)

Arguments

object

An object returned by model.avg.

newdata

An optional data frame in which to look for variables with which to predict. If omitted, the fitted values are used.

se.fit

logical, indicates if standard errors should be returned. This has any effect only if the predict methods for each of the component models support it.

interval

Currently not used.

type

Predictions on response scale are only possible if all component models use the same family. See predict.glm

...

arguments to be passed to respective predict method (e.g. level for lme model).

Value

• An object of class averaging with following elements:
• summary
  a data.frame with deviance, AICc, Delta and weights for the component models.
• coefficients, variance
  matrices of component models' coefficients and their variances
• variable.codes
  names of the variables with numerical codes used in the summary
• avg.model
  the averaged model summary, (data.frame containing averaged coefficients, unconditional standard error, adjusted SE, and confidence intervals)
• importance
  the relative importance of variables
• beta
  (logical) were standardized coefficients used?
• term.names
  character vector giving names of all terms in the model
• residuals
  the residuals (response minus fitted values).
• x, formula
  the model matrix and formula analogical to those that would be used in a single model.
• method
  how the missing terms were handled ("NA" or "0").
• call
  the matched call.

encoding

utf-8

Details

predict.averaging supports method = "NA" only for linear, fixed effect models. In other cases (e.g. nonlinear or mixed models), prediction is obtained using brute force, i.e. by calling predict on each component model and weighted averaging the results, which is equivalent to assuming that missing coefficients equal zero (method = "0").

Besides predict and coef, other generic methods such as formula, residuals and vcov are supported. logLik method returns a list of logLik objects for the component models.

See Also

model.avg See par.avg for details of model-averaged parameter calculation.

Examples

require(graphics)

# Example from Burnham and Anderson (2002), page 100:
data(Cement)
lm1 <- lm(y ~ ., data = Cement)
dd <- dredge(lm1)
top.models <- get.models(dd, subset=cumsum(weight) <= .95)
avgm <- model.avg(top.models)

# helper function
nseq <- function(x, len=length(x)) seq(min(x, na.rm=TRUE),
    max(x, na.rm=TRUE),    length=len)

# New predictors: X1 along the range of original data, other 
# variables held constant at their means
newdata <- as.data.frame(lapply(lapply(Cement[1:5], mean), rep, 25))
newdata$X1 <- nseq(Cement$X1, nrow(newdata))

# Predictions from each of the models in a set:
pred <- sapply(top.models, predict, newdata=newdata)
# Add predictions from the models averaged using two methods:
pred <- cbind(pred,
    averaged.0=predict(avgm, newdata),
    averaged.NA=predict(update(avgm, method="NA"), newdata))

matplot(x=newdata$X1, y=pred, type="l", lwd=c(rep(1,ncol(pred)-2), 2, 2),
    xlab="X1", ylab="y")

legend("topleft",
    legend=c(lapply(top.models, formula),
        paste("Averaged model (method=", c("0", "NA"), ")", sep="")),
    col=1:6, lty=1:5, lwd=c(rep(1,ncol(pred)-2), 2, 2), cex = .75)


# Example with gam models (based on "example(gam)")
require(mgcv)
dat <- gamSim(1, n = 500, dist="poisson", scale=0.1)

gam1 <- gam(y ~ s(x0) + s(x1) + s(x2) +  s(x3) + (x1 + x2 + x3)^2,
    family = poisson, data = dat, method = "REML")

cat(dQuote(getAllTerms(gam1)), "\n")

# include only models with either smooth OR linear term (but not both)
# for each predictor variable:
dd <- dredge(gam1, subset=xor(`s(x1)`, x1) & xor(`s(x2)`, x2) & 
    xor(`s(x3)`, x3))
# ...this may take a while.

subset(dd, cumsum(weight) < .95)

top.models <- get.models(dd, cumsum(weight) <= .95)

newdata <- as.data.frame(lapply(lapply(dat, mean), rep, 50))
newdata$x1 <- nseq(dat$x1, nrow(newdata))
pred <- cbind(
    sapply(top.models, predict, newdata=newdata),
    averaged=predict(model.avg(top.models), newdata))

matplot(x=newdata$x1, y=pred, type="l", xlab="x1", ylab="y"
    lwd=c(rep(1, ncol(pred) - 2), 2, 2))