predictCI: Calculate a predicted value matrix (fit, lwr, upr) for a regression, either lm or glm, on either link or response scale.

Description

This adapts code from predict.glm and predict.lm. I eliminated type = "terms" from consideration.

Usage

predictCI(object, newdata = NULL, type = c("response", "link"),
  interval = c("none", "confidence", "prediction"), dispersion = NULL,
  scale = NULL, na.action = na.pass, level = 0.95, ...)

Arguments

object

Regression object, class must include glm or lm.

newdata

Data frame including focal values for predictors

type

One of c("response", "link"), defaults to former.

interval

One of c("none", "confidence", "prediction"). "prediction" is defined only for lm objects, not for glm.

dispersion

Will be estimated if not provided. The variance coefficient of the glm, same as scale squared. Dispersion is allowed as an argument in predict.glm.

scale

The square root of dispersion. In an lm, this is the RMSE, called sigma in summary.lm.

na.action

What to do with missing values

level

Optional. Default = 0.95. Specify whatever confidence level one desires.

...

Other arguments to be passed to predict

Value

c(fit, lwr, upr), and possibly more.

Details

R's predict.glm does not have an interval argument. There are about 50 methods to calculate CIs for predicted values of GLMs, that's a major worry. This function takes the simplest route, calculating the (fit, lwr, upr) in the linear predictor scale, and then if type= "response", those 3 columns are put through linkinv(). This is the same method that SAS manuals suggest they use, same as Ben Bolker suggests in r-help (2010). I'd rather use one of the fancy tools like Edgeworth expansion, but that R code is not available (but is promised).

Use predict.lm with se.fit = TRUE to calculate fit and se.fit. Then calculate lwr and upr as fit +/- tval * se.fit. If model is lm, the model df.residual will be used to get tval. If glm, this is a normal approximation, so we thugishly assert tval = 1.98.

There's some confusing term translation. I wish R lm and glm would be brought into line. For lm, residual.scale = sigma. For glm, residual.scale = sqrt(dispersion)