feat(estimate.default): move hypothesis testing into summary.estimate method

NEWS.md

@@ -1,5 +1,10 @@
 # lava 1.9.1
   Development version
+  - `estimate.default`: the `null`, `contrast`, `type`, and `var.adj`
+    arguments are soft-deprecated. Use
+    `summary(estimate(...), null=, contrast=, type=, var.adj=)` instead.
+    The default-path Wald p-value (H0: beta = 0) continues to be reported by
+    `estimate()`.
   - `estimate.default`: removed `score.deriv` and `folds` arguments.
   - `estimate.default`: removed `R` and `null.sim` arguments
   - Safe evaluation of rank in `wald_test`

R/diagtest.R

@@ -156,7 +156,7 @@ diagtest <- function(table,positive=2,exact=FALSE,p0=NA,confint=c("logit","arcsi
                 }
             }
         }
-        res <- estimate(M,calc_diag,print=prfun,null=c(rep(p0,7),0),back.transform=btransform,...)
+        res <- suppressWarnings(estimate(M,calc_diag,print=prfun,null=c(rep(p0,7),0),back.transform=btransform,...))
     }
     CI <- confint[1]
     if (exact) CI <- "exact"

R/estimate.default.R

@@ -20,15 +20,19 @@ estimate <- function(x, ...) UseMethod("estimate")
 ##' @param level level of confidence limits
 ##' @param IC if TRUE (default) the influence function decompositions are also
 ##'   returned (extract with \code{IC} method)
-##' @param type type of small-sample correction
+##' @param type type of small-sample correction (deprecated; see
+##'   [summary.estimate])
 ##' @param var.adj variance adjustment parameter for small-sample correction
+##'   (deprecated; see [summary.estimate])
 ##' @param keep (optional) index of parameters to keep from final result
 ##' @param use (optional) index of parameters to use in calculations
 ##' @param regex If TRUE use regular expression (perl compatible) for keep, use
 ##'   arguments
 ##' @param ignore.case Ignore case-sensitiveness in regular expression
 ##' @param contrast (optional) Contrast matrix for final Wald test
-##' @param null (optional) null hypothesis to test
+##'   (deprecated; see [summary.estimate])
+##' @param null (optional) null hypothesis to test (deprecated; see
+##'   [summary.estimate])
 ##' @param vcov (optional) covariance matrix of parameter estimates or logical.
 ##' If TRUE, then [stats::vcov] is used to obtain the covariance matrix from the
 ##' provided model object `x`
@@ -66,16 +70,15 @@ estimate <- function(x, ...) UseMethod("estimate")
 ##' estimate(g)
 ##'
 ##' ## Testing contrasts
-##' estimate(g, null=0)
 ##' estimate(g, rbind(c(1,1,0), c(1,0,2)))
-##' estimate(g, rbind(c(1,1,0), c(1,0,2)), null=c(1,2))
+##' summary(estimate(g, rbind(c(1,1,0), c(1,0,2))), null=c(1,2))
 ##' estimate(g, 2:3) ## same as cbind(0,1,-1)
 ##' estimate(g, as.list(2:3)) ## same as rbind(c(0,1,0),c(0,0,1))
 ##' ## Alternative syntax
 ##' estimate(g, "z", "z"-"x", 2*"z"-3*"x")
 ##' estimate(g, "?")  ## Wildcards
 ##' estimate(g, "*Int*", "z")
-##' estimate(g, "1", "2"-"3", null = c(0,1))
+##' summary(estimate(g, "1", "2"-"3"), null = c(0,1))
 ##' estimate(g, 2, 3)
 ##'
 ##' ## Usual (non-robust) confidence intervals
@@ -219,6 +222,16 @@ estimate.default <- function(x=NULL, f=NULL, ..., data, id,
       "Use the 'vcov' argument to compute model-based SEs."
     )
   }
+  if (!missing(null) || !missing(contrast) ||
+      !missing(type) || !missing(var.adj)) {
+    .Deprecated(
+      msg = paste0(
+        "The 'null', 'contrast', 'type', and 'var.adj' arguments of ",
+        "estimate.default() are deprecated. Use ",
+        "summary(estimate(...), null=, contrast=, type=, var.adj=) instead."
+      )
+    )
+  }
   if (!missing(use)) {
     p0 <- c(
       "f", "contrast", "only.coef",
@@ -376,6 +389,8 @@ estimate.default <- function(x=NULL, f=NULL, ..., data, id,
   }
   if (inherits(x, "lm") && family(x)$family == "gaussian"
       && is.null(df) && !missing(vcov)) {
+    # defaults to t-distribution when calculating p-values with model-based
+    # SEs
     df <- x$df.residual
   }
   if (!is.null(ic_theta) && (missing(vcov) || is.null(vcov))) {
@@ -568,6 +583,8 @@ estimate.default <- function(x=NULL, f=NULL, ..., data, id,
     if (!missing(null) && missing(contrast))
       beta0 <- beta0-null
     if (!is.null(df)) {
+      # the idea should be to assign df as an attribute to the returned
+      # object, and re-use it inside summary.estimate
       za <- qt(1-alpha/2, df=df)
       pval <- 2*pt(abs(res[, 1]/res[, 2]), df=df, lower.tail=FALSE)
     } else {
@@ -596,7 +613,7 @@ estimate.default <- function(x=NULL, f=NULL, ..., data, id,
     names(coefs) <- rownames(res)
   }
   res <- structure(list(coef=coefs, coefmat=res, vcov=V,
-                        IC=NULL, print=print, id=idstack),
+                        IC=NULL, print=print, id=idstack, df=df),
                    class="estimate")
   if (IC) ## && is.null(back.transform))
     res$IC <- ic_theta
@@ -636,7 +653,7 @@ estimate.default <- function(x=NULL, f=NULL, ..., data, id,
       rownames(res$coefmat) <- names(res$coef)
     }
   }
-
+  
   if (!is.null(back.transform)) {
     res$coefmat[, c(1, 3, 4)] <- do.call(back.transform,
                                          list(res$coefmat[, c(1, 3, 4)]))
@@ -821,20 +838,78 @@ with_unique_warnings <- function(expr) {
   })
 }
 
+##' Summary of estimate objects
+##'
+##' Computes hypothesis tests, contrasts, and small-sample corrections for
+##' an [estimate] object. The arguments `null`, `contrast`, `type`, and
+##' `var.adj` were previously available on [estimate.default()] and have
+##' been moved here.
+##'
+##' @param object an `estimate` object.
+##' @param contrast (optional) contrast matrix for the final Wald test.
+##' @param null (optional) null hypothesis to test.
+##' @param type type of small-sample correction. Requires the estimate
+##'   to have been computed with `IC=TRUE` (the default).
+##' @param var.adj variance adjustment parameter for small-sample
+##'   correction. Requires the estimate to have been computed with
+##'   `IC=TRUE` (the default).
+##' @param ... additional arguments passed to [estimate()].
+##' @seealso [estimate.default()]
 ##' @export
 summary.estimate <- function(object,
                              contrast,
+                             null,
+                             type,
+                             var.adj,
                              ...) {
   with_unique_warnings({
     p <- coef(object)
-    if (missing(contrast)) contrast <- diag(1,nrow=length(p))#as.list(seq_along(p))
-    test <- estimate(coef=p,
-                     vcov=vcov(object),
-                     f = FALSE,
-                     contrast = contrast,
-                     ...)
+    df <- object$df
+    correction <- !missing(type) || !missing(var.adj)
+    user_contrast <- !missing(contrast)
+    user_null <- !missing(null)
+    if (missing(contrast)) contrast <- diag(1, nrow=length(p))
+    args <- list(
+      coef = p,
+      contrast = contrast,
+      df = df
+    )
+    ## When the user supplies any of null/contrast/type/var.adj we want
+    ## the recall to run the full Wald-test override (contrast.transform
+    ## must be TRUE). In the no-op case we preserve the historical
+    ## behavior of passing f=FALSE so the H0:beta=0 p-values in coefmat
+    ## are kept verbatim.
+    if (!user_contrast && !user_null && !correction) {
+      args$f <- FALSE
+    }
+    if (user_null) args$null <- null
+    if (correction) {
+      if (is.null(object$IC)) {
+        stop(
+          "Small-sample corrections in summary() require an estimate ",
+          "computed with IC=TRUE."
+        )
+      }
+      args$IC <- object$IC
+      args$id <- object$id
+      if (!missing(type)) args$type <- type
+      if (!missing(var.adj)) args$var.adj <- var.adj
+    } else {
+      args$vcov <- vcov(object)
+    }
+    args <- c(args, list(...))
+    test <- suppressWarnings(do.call(estimate, args))
     class(test) <- "NULL"
     test$compare <- test$compare
+    ## Preserve dimnames on vcov: the recall path through type/var.adj
+    ## corrections builds V from the influence function and may drop
+    ## dimnames. Restore them from the original coef names so the
+    ## summary output is consistent with the deprecated estimate() path.
+    if (!is.null(test$vcov) && is.null(dimnames(test$vcov)) &&
+        !is.null(names(test$coef)) &&
+        nrow(test$vcov) == length(test$coef)) {
+      dimnames(test$vcov) <- list(names(test$coef), names(test$coef))
+    }
     object <- test[c("coef", "coefmat", "vcov", "call",
                      "ncluster", "model.index", "compare")]
     class(object) <- "summary.estimate"

R/estimate_calculus.R

@@ -584,5 +584,5 @@ operator_grad <- function(x, y, x_const, y_const, dx, dy) {
   if (!(is.numeric(e1) || is.numeric(e2))) stop("numeric comperator needed")
   null <- if (is.numeric(e1)) e1 else e2
   e <- if (is.numeric(e1)) e2 else e1
-  estimate(e, null=null)
+  suppressWarnings(estimate(e, null=null))
 }

tests/testthat/test-estimate-deprecation.R

@@ -0,0 +1,65 @@
+## Tests for soft deprecation of null/contrast/type/var.adj in estimate.default
+
+context("estimate.default deprecation warnings")
+
+make_glm <- function() {
+  set.seed(1)
+  d <- data.frame(
+    y = rbinom(100, 1, 0.5),
+    x = rnorm(100),
+    z = rnorm(100)
+  )
+  glm(y ~ x + z, data = d, family = binomial())
+}
+
+test_that("estimate(..., null=) emits deprecation warning", {
+  g <- make_glm()
+  expect_warning(estimate(g, null = 0), class = "deprecatedWarning")
+})
+
+test_that("estimate(..., null=0) warns even when null=0 (Q7)", {
+  g <- make_glm()
+  expect_warning(estimate(g, null = 0), class = "deprecatedWarning")
+})
+
+test_that("estimate(..., contrast=) emits deprecation warning", {
+  g <- make_glm()
+  B <- rbind(c(1, 0, 0), c(0, 1, 0))
+  expect_warning(estimate(g, contrast = B), class = "deprecatedWarning")
+})
+
+test_that("estimate(..., type=) emits deprecation warning", {
+  g <- make_glm()
+  expect_warning(estimate(g, type = "df"), class = "deprecatedWarning")
+})
+
+test_that("estimate(..., var.adj=) emits deprecation warning", {
+  g <- make_glm()
+  expect_warning(
+    estimate(g, type = "hc3", var.adj = 0.5),
+    class = "deprecatedWarning"
+  )
+})
+
+test_that("multiple deprecated args produce exactly one deprecation warning", {
+  g <- make_glm()
+  ws <- capture_warnings(estimate(g, null = 1, type = "df"))
+  dep_msgs <- grep("deprecated", ws, value = TRUE)
+  expect_equal(length(dep_msgs), 1L)
+})
+
+test_that("internally-derived contrast from f does NOT warn (Q8)", {
+  g <- make_glm()
+  ws <- capture_warnings(estimate(g, "x"))
+  expect_equal(length(grep("deprecated", ws)), 0L)
+  ws <- capture_warnings(estimate(g, 2:3))
+  expect_equal(length(grep("deprecated", ws)), 0L)
+  ws <- capture_warnings(estimate(g, rbind(c(1, 0, 0), c(0, 1, 0))))
+  expect_equal(length(grep("deprecated", ws)), 0L)
+})
+
+test_that("plain estimate() emits no deprecation warning (Q4)", {
+  g <- make_glm()
+  ws <- capture_warnings(estimate(g))
+  expect_equal(length(grep("deprecated", ws)), 0L)
+})