[r-cran-estimability] 01/01: Imported Upstream version 1.1
Jonathon Love
jon at thon.cc
Fri Aug 28 10:58:16 UTC 2015
This is an automated email from the git hooks/post-receive script.
jonathon-guest pushed a commit to branch master
in repository r-cran-estimability.
commit 8da926d7c34fb73817778b9bc3ebdc4340e8cce5
Author: Jonathon Love <jon at thon.cc>
Date: Fri Aug 28 06:52:27 2015 -0400
Imported Upstream version 1.1
---
DESCRIPTION | 16 ++++++
MD5 | 8 +++
NAMESPACE | 17 ++++++
R/epredict.lm.R | 97 +++++++++++++++++++++++++++++++++++
R/estimability.R | 63 +++++++++++++++++++++++
inst/NEWS | 15 ++++++
man/epredict.lm.Rd | 122 ++++++++++++++++++++++++++++++++++++++++++++
man/estimability-package.Rd | 32 ++++++++++++
man/nonest.basis.Rd | 76 +++++++++++++++++++++++++++
9 files changed, 446 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..df70f59
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,16 @@
+Package: estimability
+Type: Package
+Title: Estimability Tools for Linear Models
+Version: 1.1
+Date: 2015-02-11
+Author: Russell V. Lenth
+Maintainer: Russell V. Lenth <russell-lenth at uiowa.edu>
+Depends: stats
+Description: Provides tools for determining estimability of linear functions of regression coefficients,
+ and 'epredict' methods that handle non-estimable cases correctly.
+LazyData: Yes
+License: GPL-2
+Packaged: 2015-02-11 23:07:28 UTC; rlenth
+NeedsCompilation: no
+Repository: CRAN
+Date/Publication: 2015-02-12 07:39:33
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..6534e1b
--- /dev/null
+++ b/MD5
@@ -0,0 +1,8 @@
+1aeae29d661228c4b1b2d4be9d53cdbd *DESCRIPTION
+df62720bde52c7b274298c8425d6f365 *NAMESPACE
+ef35de60ea07685e778cc287e71b5289 *R/epredict.lm.R
+68cf3128cfe31ab716169b84ef7f6066 *R/estimability.R
+14f2193f4249591397d66e700503bc95 *inst/NEWS
+894059c01d85207fb5cbfe01e399a45c *man/epredict.lm.Rd
+34dfc7c7a3beb3e4ab147120840da44b *man/estimability-package.Rd
+3a4b30946a2cea7c393f1c58c18859a0 *man/nonest.basis.Rd
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..1d7e492
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,17 @@
+# Exports from estimability package
+
+export(all.estble)
+export(epredict)
+export(eupdate)
+export(is.estble)
+export(nonest.basis)
+
+S3method(epredict, lm)
+S3method(epredict, mlm)
+S3method(epredict, glm)
+
+S3method(eupdate, lm)
+
+S3method(nonest.basis, qr)
+S3method(nonest.basis, matrix)
+S3method(nonest.basis, lm)
diff --git a/R/epredict.lm.R b/R/epredict.lm.R
new file mode 100644
index 0000000..61841aa
--- /dev/null
+++ b/R/epredict.lm.R
@@ -0,0 +1,97 @@
+# Patch for predict.lm, predict.glm, predict.mlm
+#
+# If newdata present, and fit is rank-deficient,
+# we check estimability and replace any non-est cases with NA
+# Use options(estimability.quiet = TRUE) to suppress message
+# Use options(estimability.suppress = TRUE) to override this patch
+
+# Main workhorse -- call with stats-library predict function
+.patch.predict = function(object, newdata, type,
+ nonest.tol = 1e-8, nbasis = object$nonest, ...) {
+
+ if(missing(newdata))
+ predict(object = object, type = type, ...)
+
+ else {
+ type = match.arg(type, c("response", "link", "terms", "matrix", "estimability"))
+ if (all(!is.na(object$coefficients)) && (type != "matrix"))
+ if (type == "estimability")
+ return (rep(TRUE, nrow(newdata)))
+ else
+ return (predict(object = object, newdata = newdata, type = type, ...))
+
+ if(is.null(nbasis)) {
+ if (!is.null(qr <- object$qr))
+ nbasis = nonest.basis(qr)
+ else
+ nbasis = nonest.basis(model.matrix(object))
+ }
+ trms = delete.response(terms(object))
+ m = model.frame(trms, newdata, na.action = na.pass, xlev = object$xlevels)
+ X = model.matrix(trms, m, contrasts.arg = object$contrasts)
+ nonest = !is.estble(X, nbasis, nonest.tol)
+
+ if (type == "estimability")
+ return (!nonest)
+ else if (type == "matrix") {
+ attr(X, "estble") = !nonest
+ return (X)
+ }
+ # (else) we have a type anticipated by stats::predict
+
+ w.handler <- function(w){ # suppress the incorrect warning
+ if (!is.na(pmatch("prediction from a rank-deficient", w$message)))
+ invokeRestart("muffleWarning")
+ }
+ result = withCallingHandlers(
+ predict(object = object, newdata = newdata, type = type, ...),
+ warning = w.handler)
+
+ if (any(nonest)) {
+ if (is.matrix(result))
+ result[nonest, ] = NA
+ else if (is.list(result)) {
+ result$fit[nonest] = NA
+ result$se.fit[nonest] = NA
+ }
+ else
+ result[nonest] = NA
+ if(getOption("estimability.quiet", FALSE))
+ message("Non-estimable cases are replaced by 'NA'")
+ }
+
+ result
+ }
+}
+
+
+# Generic for epredict
+epredict = function(object, ...)
+ UseMethod("epredict")
+
+epredict.lm = function(object, newdata, ...,
+ type = c("response", "terms", "matrix", "estimability"),
+ nonest.tol = 1e-8, nbasis = object$nonest)
+ .patch.predict(object, newdata, type[1], nonest.tol, nbasis, ...)
+
+epredict.glm = function(object, newdata, ...,
+ type = c("link", "response", "terms", "matrix", "estimability"),
+ nonest.tol = 1e-8, nbasis = object$nonest)
+ .patch.predict(object, newdata, type[1], nonest.tol, nbasis, ...)
+
+epredict.mlm = function(object, newdata, ...,
+ type = c("response", "matrix", "estimability"),
+ nonest.tol = 1e-8, nbasis = object$nonest)
+ .patch.predict(object, newdata, type[1], nonest.tol, nbasis, ...)
+
+
+# Generic for eupdate -- adds nonest basis to object
+eupdate = function(object, ...)
+ UseMethod("eupdate")
+
+eupdate.lm = function(object, ...) {
+ if (length(list(...)) > 0)
+ object = do.call("update", list(object = object, ...))
+ object$nonest = nonest.basis(object)
+ object
+}
diff --git a/R/estimability.R b/R/estimability.R
new file mode 100644
index 0000000..54d2aba
--- /dev/null
+++ b/R/estimability.R
@@ -0,0 +1,63 @@
+# Obtain an orthonormal basis for nonestimable functions
+
+# Generic
+nonest.basis = function(x, ...)
+ UseMethod("nonest.basis")
+
+
+# Main method -- for class "qr"
+nonest.basis.qr = function(x, ...) {
+ rank = x$rank
+ tR = t(qr.R(x))
+ p = nrow(tR)
+ if (rank == p)
+ return (all.estble)
+
+ # null space of X is same as null space of R in QR decomp
+ if (ncol(tR) < p) # add columns if not square
+ tR = cbind(tR, matrix(0, nrow=p, ncol=p-ncol(tR)))
+
+ # last few rows are zero -- add a diagonal of 1s
+ extras = rank + seq_len(p - rank)
+ tR[extras, extras] = diag(1, p - rank)
+
+ # nbasis is last p - rank cols of Q in QR decomp of tR
+ nbasis = qr.Q(qr(tR))[ , extras, drop = FALSE]
+
+ # permute the rows via pivot
+ nbasis[x$pivot, ] = nbasis
+
+ nbasis
+}
+
+
+nonest.basis.matrix = function(x, ...)
+ nonest.basis(qr(x), ...)
+
+
+nonest.basis.lm = function(x, ...) {
+ if (is.null(x$qr))
+ x = update(x, method = "qr", qr = TRUE)
+ nonest.basis(x$qr)
+}
+
+
+
+# utility to check estimability of x'beta, given nonest.basis
+is.estble = function(x, nbasis, tol = 1e-8) {
+ if (is.matrix(x))
+ return(apply(x, 1, is.estble, nbasis, tol))
+ if(is.na(nbasis[1]))
+ TRUE
+ else {
+ chk = as.numeric(crossprod(nbasis, x))
+ ssqx = sum(x*x) # BEFORE subsetting x
+ # If x really small, don't scale chk'chk
+ if (ssqx < tol) ssqx = 1
+ sum(chk*chk) < tol * ssqx
+ }
+}
+
+
+# nonestimability basis that makes everything estimable
+all.estble = matrix(NA)
diff --git a/inst/NEWS b/inst/NEWS
new file mode 100644
index 0000000..4686e5c
--- /dev/null
+++ b/inst/NEWS
@@ -0,0 +1,15 @@
+Update history for **estimability***
+
+1.1
+ Design improvements to aid in potential scope and usability:
+ * Made 'nonest.basis' a generic, with provided methods for
+ "qr", "matrix", and "lm"
+ * Added 'eupdate' generic and 'lm' method for updating a
+ model object and including its nonestimability basis as
+ part of the object
+ Added 'type = "matrix"' and 'type = "estimability"'
+ options for 'epredict'
+
+
+1.0-2
+ Initial version on CRAN
diff --git a/man/epredict.lm.Rd b/man/epredict.lm.Rd
new file mode 100644
index 0000000..cbc4134
--- /dev/null
+++ b/man/epredict.lm.Rd
@@ -0,0 +1,122 @@
+\name{epredict}
+\alias{epredict}
+\alias{epredict.lm}
+\alias{epredict.glm}
+\alias{epredict.mlm}
+\alias{eupdate}
+\alias{eupdate.lm}
+
+\title{
+Estimability Enhancements for \code{lm} and Relatives
+}
+\description{
+These functions call the corresponding S3 \code{predict} methods in the \pkg{stats} package, but with a check for estimability of new predictions, and with appropriate actions for non-estimable cases.
+}
+\usage{
+\S3method{epredict}{lm}(object, newdata, ...,
+ type = c("response", "terms", "matrix", "estimability"),
+ nonest.tol = 1e-8, nbasis = object$nonest)
+\S3method{epredict}{glm}(object, newdata, ...,
+ type = c("link", "response", "terms", "matrix", "estimability"),
+ nonest.tol = 1e-8, nbasis = object$nonest)
+\S3method{epredict}{mlm}(object, newdata, ...,
+ type = c("response", "matrix", "estimability"),
+ nonest.tol = 1e-8, nbasis = object$nonest)
+
+eupdate(object, ...)
+}
+\arguments{
+ \item{object}{An object ingeriting from \code{lm}}
+ \item{newdata}{A \code{data.frame} containing predictor combinations for new predictions}
+ \item{\dots}{Arguments passed to \code{\link{predict}} or \code{\link{update}}}
+ \item{nonest.tol}{Tolerance used by \code{\link{is.estble}} to check estimability of new predictions}
+ \item{type}{Character string specifying the desired result. See Details.}
+ \item{nbasis}{Basis for the null space, e.g., a result of a call to \code{\link{nonest.basis}}. If \code{nbasis} is \code{NULL}, a basis is constructed from \code{object}.}
+}
+\details{
+If \code{newdata} is missing or \code{object} is not rank-deficient, this method passes its arguments directly to the same method in the \pkg{stats} library. In rank-deficient cases with \code{newdata} provided, each row of \code{newdata} is tested for estimability against the null basis provided in \code{nbasis}. Any non-estimable cases found are replaced with \code{NA}s.
+
+The \code{type} argument is passed to \code{\link[stats]{predict}} when it is one of \code{"response"}, \code{"link"}, or \code{"terms"}. With \code{newdata} present and \code{type = "matrix"}, the model matrix for \code{newdata} is returned, with an attribute \code{"estble"} that is a logical vector of length \samp{nrow(newdata)} indicating whether each row is estimable. With \code{type = "estimability"}, just the logical vector is returned.
+
+If you anticipate making several \code{epredict} calls with new data, it improves efficiency to either obtain the null basis and provide it in the call, or add it to \code{object} with the name \code{"nonest"} (perhaps via a call to \code{eupdate}).
+
+\code{eupdate} is an S3 generic function with a method provided for \code{"lm"} objects. It updates the object according to any arguments in \code{...}, then obtains the updated object's nonestimable basis and returns it in \code{object$nonest}.
+}
+
+\value{
+The same as the result of a call to the \code{predict} method in the \pkg{stats} package, except rows or elements corresponding to non-estimable predictor combinations are set to \code{NA}. The value for \code{type} is \code{"matrix"} or \code{"estimability"} is explained under details.}
+
+\author{
+Russell V. Lenth <russell-lenth at uiowa.edu>
+}
+
+\note{
+The usual rank-deficiency warning from \code{stats::predict} is suppressed; but when non-estimable cases are found, a message is displayed explaining that these results were replaced by \code{NA}. If you wish that message suppressed, use \samp{options(estimability.quiet = TRUE)}.
+}
+
+
+\seealso{
+ \code{\link[stats]{predict.lm}} in the \pkg{stats} package;
+ \code{\link{nonest.basis}}.
+}
+\examples{
+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), ])
+
+\dontrun{
+### 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 of \dontrun
+
+}
+
+\keyword{ models }
+\keyword{ regression }
diff --git a/man/estimability-package.Rd b/man/estimability-package.Rd
new file mode 100644
index 0000000..1d3b504
--- /dev/null
+++ b/man/estimability-package.Rd
@@ -0,0 +1,32 @@
+\name{estimability-package}
+\alias{estimability-package}
+\alias{estimability}
+\docType{package}
+\title{
+Estimability Tools for Linear Models
+}
+\description{
+Provides tools for determining estimability of linear functions of regression coefficients,
+and alternative \code{epredict} methods for \code{lm}, \code{glm}, and \code{mlm} objects that handle non-estimable cases correctly.
+}
+\details{
+\tabular{ll}{
+Package: \tab estimability\cr
+Type: \tab Package\cr
+Details: \tab See DESCRIPTION file\cr
+}
+When a linear model is not of full rank, the regression coefficients are not uniquely estimable. However, the predicted values are unique, as are other linear combinations where the coefficients lie in the row space of the data matrix. Thus, estimability of a linear function of regression coefficients can be determined by testing whether the coefficients lie in this row space -- or equivalently, are orthogonal to the corresponding null space.
+This package provides functions \code{\link{nonest.basis}} and \code{\link{is.estble}} to facilitate such an estimability test. Package developers may find these useful for incorporating in their \code{predict} methods when new predictor settings are involved.
+
+The package also provides \code{\link{epredict}} methods -- alternatives to the \code{\link{predict}} methods in the \pkg{stats} package for \code{"lm"}, \code{"glm"}, and \code{"mlm"} objects. When the \code{newdata} argument is specified, estimability of each new prediction is checked and any non-estimable cases are replaced by \code{NA}.
+}
+\author{
+Russell V. Lenth <russell-lenth at uiowa.edu>
+}
+\references{
+Monahan, John F. (2008) \emph{A Primer on Linear Models}, CRC Press. (Chapter 3)
+}
+
+\keyword{ package }
+\keyword{ models }
+\keyword{ regression }
diff --git a/man/nonest.basis.Rd b/man/nonest.basis.Rd
new file mode 100644
index 0000000..9299d6f
--- /dev/null
+++ b/man/nonest.basis.Rd
@@ -0,0 +1,76 @@
+\name{nonest.basis}
+\alias{nonest.basis}
+\alias{nonest.basis.qr}
+\alias{nonest.basis.matrix}
+\alias{nonest.basis.lm}
+\alias{all.estble}
+\alias{is.estble}
+
+\title{Estimability Tools}
+\description{
+This documents the functions needed to test estimability of linear functions of regression coefficients.
+}
+\usage{
+nonest.basis(x, ...)
+\S3method{nonest.basis}{qr}(x, ...)
+\S3method{nonest.basis}{matrix}(x, ...)
+\S3method{nonest.basis}{lm}(x, ...)
+
+all.estble
+
+is.estble(x, nbasis, tol = 1e-8)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{x}{For \code{nonest.basis}, an object of a class in \samp{methods("nonest.basis")}. Or, in \code{is.estble}, a numeric vector or matrix for assessing estimability of \samp{sum(x * beta)}, where \code{beta} is the vector of regression coefficients.}
+ \item{nbasis}{Matrix whose columns span the null space of the model matrix. Such a matrix is returned by \code{nonest.basis}.}
+ \item{tol}{Numeric tolerance for assessing nonestimability. For nonzero \eqn{x}, estimability of \eqn{\beta'x} is assessed by whether or not \eqn{||N'x||^2 < \tau ||x'x||^2}, where \eqn{N} and \eqn{\tau} denote \code{nbasis} and \code{tol}, respectively.}
+ \item{\dots}{Additional arguments, currently ignored.}
+}
+
+\details{
+Consider a linear model \eqn{y = X\beta + E}. If \eqn{X} is not of full rank, it is not possible to estimate \eqn{\beta} uniquely. However, \eqn{X\beta} \emph{is} uniquely estimable, and so is \eqn{a'X\beta} for any conformable vector \eqn{a}. Since \eqn{a'X} comprises a linear combination of the rows of \eqn{X}, it follows that we can estimate any linear function where the coefficients lie in the row space of \eqn{X}. Equivalently, we can check to ensure that the coefficients are orthog [...]
+
+The constant \code{all.estble} is simply a 1 x 1 matrix of \code{NA}. This specifies a trivial non-estimability basis, and using it as \code{nbasis} will cause everything to test as estimable.
+}
+
+
+\value{
+When \eqn{X} is not full-rank, the methods for \code{nonest.basis} return a basis for the null space of \eqn{X}. The number of rows is equal to the number of regression coefficients (\emph{including} any \code{NA}s); and the number of columns is equal to the rank deficiency of the model matrix. The columns are orthonormal. If the model is full-rank, then \code{nonest.basis} returns \code{all.estble}. The \code{matrix} method uses \eqn{X} itself, the \code{qr} method uses the \eqn{QR} dec [...]
+
+The function \code{is.estble} returns a logical value (or vector, if \code{x} is a matrix) that is \code{TRUE} if the function is estimable and \code{FALSE} if not.
+}
+
+\references{
+Monahan, John F. (2008) \emph{A Primer on Linear Models}, CRC Press. (Chapter 3)
+}
+
+\author{
+Russell V. Lenth <russell-lenth at uiowa.edu>
+}
+
+
+\examples{
+require(estimability)
+
+X <- cbind(rep(1,5), 1:5, 5:1, 2:6)
+( nb <- nonest.basis(X) )
+
+# Test estimability of some linear functions for this X matrix
+lfs <- rbind(c(1,4,2,5), c(2,3,9,5), c(1,2,2,1), c(0,1,-1,1))
+is.estble(lfs, nb)
+
+# Illustration on 'lm' objects:
+warp.lm1 <- lm(breaks ~ wool * tension, data = warpbreaks,
+ subset = -(26:38),
+ contrasts = list(wool = "contr.treatment", tension = "contr.treatment"))
+zapsmall(nonest.basis(warp.lm1))
+
+warp.lm2 <- update(warp.lm1,
+ contrasts = list(wool = "contr.sum", tension = "contr.helmert"))
+zapsmall(nonest.basis(warp.lm2))
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ models }
+\keyword{ regression }
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/r-cran-estimability.git
More information about the debian-science-commits
mailing list