Estimability Enhancements for lm and Relatives
These functions call the corresponding S3 predict methods in the stats package, but with a check for estimability of new predictions, and with appropriate actions for non-estimable cases.
## S3 method for class 'lm'
epredict(object, newdata, ...,
type = c("response", "terms", "matrix", "estimability"),
nonest.tol = 1e-8, nbasis = object$nonest)
## S3 method for class 'glm'
epredict(object, newdata, ...,
type = c("link", "response", "terms", "matrix", "estimability"),
nonest.tol = 1e-8, nbasis = object$nonest)
## S3 method for class 'mlm'
epredict(object, newdata, ...,
type = c("response", "matrix", "estimability"),
nonest.tol = 1e-8, nbasis = object$nonest)
eupdate(object, ...)object |
An object ingeriting from |
newdata |
A |
... |
|
nonest.tol |
Tolerance used by |
type |
Character string specifying the desired result. See Details. |
nbasis |
Basis for the null space, e.g., a result of a call to |
If newdata is missing or object is not rank-deficient, this method passes its arguments directly to the same method in the stats library. In rank-deficient cases with newdata provided, each row of newdata is tested for estimability against the null basis provided in nbasis. Any non-estimable cases found are replaced with NAs.
The type argument is passed to predict when it is one of "response", "link", or "terms". With newdata present and type = "matrix", the model matrix for newdata is returned, with an attribute "estble" that is a logical vector of length nrow(newdata) indicating whether each row is estimable. With type = "estimability", just the logical vector is returned.
If you anticipate making several epredict calls with new data, it improves efficiency to either obtain the null basis and provide it in the call, or add it to object with the name "nonest" (perhaps via a call to eupdate).
eupdate is an S3 generic function with a method provided for "lm" objects. It updates the object according to any arguments in ..., then obtains the updated object's nonestimable basis and returns it in object$nonest.
The same as the result of a call to the predict method in the stats package, except rows or elements corresponding to non-estimable predictor combinations are set to NA. The value for type is "matrix" or "estimability" is explained under details.
The usual rank-deficiency warning from stats::predict is suppressed; but when non-estimable cases are found, a message is displayed explaining that these results were replaced by NA. If you wish that message suppressed, use options(estimability.quiet = TRUE).
Russell V. Lenth <russell-lenth@uiowa.edu>
predict.lm in the stats package;
nonest.basis.
require("estimability")
# Fake data where x3 and x4 depend on x1, x2, and intercept
x1 <- -4:4
x2 <- c(-2,1,-1,2,0,2,-1,1,-2)
x3 <- 3*x1 - 2*x2
x4 <- x2 - x1 + 4
y <- 1 + x1 + x2 + x3 + x4 + c(-.5,.5,.5,-.5,0,.5,-.5,-.5,.5)
# Different orderings of predictors produce different solutions
mod1234 <- lm(y ~ x1 + x2 + x3 + x4)
mod4321 <- eupdate(lm(y ~ x4 + x3 + x2 + x1))
# (Estimability checking with mod4321 will be more efficient because
# it will not need to recreate the basis)
mod4321$nonest
# test data:
testset <- data.frame(
x1 = c(3, 6, 6, 0, 0, 1),
x2 = c(1, 2, 2, 0, 0, 2),
x3 = c(7, 14, 14, 0, 0, 3),
x4 = c(2, 4, 0, 4, 0, 4))
# Look at predictions when we don't check estimability
suppressWarnings( # Disable the warning from stats::predict.lm
rbind(p1234 = predict(mod1234, newdata = testset),
p4321 = predict(mod4321, newdata = testset)))
# Compare with results when we do check:
rbind(p1234 = epredict(mod1234, newdata = testset),
p4321 = epredict(mod4321, newdata = testset))
# Note that estimable cases have the same predictions
# change mod1234 and include nonest basis
mod134 <- eupdate(mod1234, . ~ . - x2, subset = -c(3, 7))
mod134$nonest
# When row spaces are the same, bases are interchangeable
# so long as you account for the ordering of parameters:
epredict(mod4321, newdata = testset, type = "estimability",
nbasis = nonest.basis(mod1234)[c(1,5:2), ])
## Not run:
### Additional illustration
example(nonest.basis) ## creates model objects warp.lm1 and warp.lm2
# The two models have different contrast specs. But the empty cell
# is correctly identified in both:
fac.cmb <- expand.grid(wool = c("A", "B"), tension = c("L", "M", "H"))
cbind(fac.cmb,
pred1 = epredict(warp.lm1, newdata = fac.cmb),
pred2 = epredict(warp.lm2, newdata = fac.cmb))
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.