diff --git a/NEWS.md b/NEWS.md index 2ac125f9..412bdfb6 100644 --- a/NEWS.md +++ b/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` diff --git a/R/diagtest.R b/R/diagtest.R index a4b9e4d1..5c9dde14 100644 --- a/R/diagtest.R +++ b/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" diff --git a/R/estimate.default.R b/R/estimate.default.R index ce9780b9..98993c8b 100644 --- a/R/estimate.default.R +++ b/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" diff --git a/R/estimate_calculus.R b/R/estimate_calculus.R index 47416459..7c5d7963 100644 --- a/R/estimate_calculus.R +++ b/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)) } diff --git a/man/estimate.default.Rd b/man/estimate.default.Rd index 91c019c1..9a62a268 100644 --- a/man/estimate.default.Rd +++ b/man/estimate.default.Rd @@ -62,9 +62,11 @@ expression or variable name)} \item{IC}{if TRUE (default) the influence function decompositions are also returned (extract with \code{IC} method)} -\item{type}{type of small-sample correction} +\item{type}{type of small-sample correction (deprecated; see +\link{summary.estimate})} -\item{var.adj}{variance adjustment parameter for small-sample correction} +\item{var.adj}{variance adjustment parameter for small-sample correction +(deprecated; see \link{summary.estimate})} \item{keep}{(optional) index of parameters to keep from final result} @@ -75,9 +77,11 @@ arguments} \item{ignore.case}{Ignore case-sensitiveness in regular expression} -\item{contrast}{(optional) Contrast matrix for final Wald test} +\item{contrast}{(optional) Contrast matrix for final Wald test +(deprecated; see \link{summary.estimate})} -\item{null}{(optional) null hypothesis to test} +\item{null}{(optional) null hypothesis to test (deprecated; see +\link{summary.estimate})} \item{vcov}{(optional) covariance matrix of parameter estimates or logical. If TRUE, then \link[stats:vcov]{stats::vcov} is used to obtain the covariance matrix from the @@ -125,16 +129,15 @@ estimate(g, g0) 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 diff --git a/man/summary.estimate.Rd b/man/summary.estimate.Rd new file mode 100644 index 00000000..b6f9710c --- /dev/null +++ b/man/summary.estimate.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/estimate.default.R +\name{summary.estimate} +\alias{summary.estimate} +\title{Summary of estimate objects} +\usage{ +\method{summary}{estimate}(object, contrast, null, type, var.adj, ...) +} +\arguments{ +\item{object}{an \code{estimate} object.} + +\item{contrast}{(optional) contrast matrix for the final Wald test.} + +\item{null}{(optional) null hypothesis to test.} + +\item{type}{type of small-sample correction. Requires the estimate +to have been computed with \code{IC=TRUE} (the default).} + +\item{var.adj}{variance adjustment parameter for small-sample +correction. Requires the estimate to have been computed with +\code{IC=TRUE} (the default).} + +\item{...}{additional arguments passed to \code{\link[=estimate]{estimate()}}.} +} +\description{ +Computes hypothesis tests, contrasts, and small-sample corrections for +an \link{estimate} object. The arguments \code{null}, \code{contrast}, \code{type}, and +\code{var.adj} were previously available on \code{\link[=estimate.default]{estimate.default()}} and have +been moved here. +} +\seealso{ +\code{\link[=estimate.default]{estimate.default()}} +} diff --git a/tests/testthat/test-estimate-deprecation.R b/tests/testthat/test-estimate-deprecation.R new file mode 100644 index 00000000..721fcad3 --- /dev/null +++ b/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) +}) diff --git a/tests/testthat/test-estimate-summary-equivalence.R b/tests/testthat/test-estimate-summary-equivalence.R new file mode 100644 index 00000000..ad6f11b2 --- /dev/null +++ b/tests/testthat/test-estimate-summary-equivalence.R @@ -0,0 +1,91 @@ +## Equivalence tests: estimate(g, X=...) should produce the same coefmat/vcov +## as summary(estimate(g), X=...) during the soft-deprecation window. + +context("estimate / summary equivalence for migrated args") + +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()) +} + +## Helper: get coefmat from either deprecated or new path +deprecated_estimate <- function(...) { + suppressWarnings(estimate(...)) +} + +test_that("null= gives identical coefmat via estimate() and summary()", { + g <- make_glm() + e1 <- deprecated_estimate(g, null = 1) + e2 <- summary(estimate(g), null = 1) + expect_equal(e1$coefmat, e2$coefmat) + expect_equal(unname(e1$vcov), unname(e2$vcov)) +}) + +test_that("contrast= gives identical coefmat via estimate() and summary()", { + g <- make_glm() + B <- rbind(c(1, 0, 0), c(0, 1, -1)) + e1 <- deprecated_estimate(g, contrast = B) + e2 <- summary(estimate(g), contrast = B) + expect_equal(e1$coefmat, e2$coefmat) + expect_equal(unname(e1$vcov), unname(e2$vcov)) +}) + +test_that("null + contrast gives identical results via estimate() and summary()", { + g <- make_glm() + B <- rbind(c(1, 0, 0), c(0, 1, -1)) + e1 <- deprecated_estimate(g, contrast = B, null = c(0, 0)) + e2 <- summary(estimate(g), contrast = B, null = c(0, 0)) + expect_equal(e1$coefmat, e2$coefmat) + expect_equal(unname(e1$vcov), unname(e2$vcov)) +}) + +test_that("type='df' gives identical coefmat via estimate() and summary()", { + g <- make_glm() + e1 <- deprecated_estimate(g, type = "df") + e2 <- summary(estimate(g), type = "df") + expect_equal(e1$coefmat, e2$coefmat) + expect_equal(unname(e1$vcov), unname(e2$vcov)) +}) + +test_that("type='df1' gives identical coefmat via estimate() and summary()", { + g <- make_glm() + e1 <- deprecated_estimate(g, type = "df1") + e2 <- summary(estimate(g), type = "df1") + expect_equal(e1$coefmat, e2$coefmat) +}) + +test_that("type='hc3' + var.adj gives identical results via estimate() and summary()", { + g <- make_glm() + e1 <- deprecated_estimate(g, type = "hc3", var.adj = 0.5) + e2 <- summary(estimate(g), type = "hc3", var.adj = 0.5) + expect_equal(e1$coefmat, e2$coefmat) + expect_equal(unname(e1$vcov), unname(e2$vcov)) +}) + +test_that("type + null gives identical results via estimate() and summary()", { + g <- make_glm() + e1 <- deprecated_estimate(g, type = "df", null = 0.2) + e2 <- summary(estimate(g), type = "df", null = 0.2) + expect_equal(e1$coefmat, e2$coefmat) +}) + +test_that("df + null gives identical results via estimate() and summary()", { + g <- make_glm() + ## Per Q11: df is stored on the estimate object and reused by summary(). + e1 <- deprecated_estimate(g, df = 10, null = 1) + e2 <- summary(estimate(g, df = 10), null = 1) + expect_equal(e1$coefmat, e2$coefmat) +}) + +test_that("df is stored on the estimate object", { + g <- make_glm() + e <- estimate(g, df = 7) + expect_equal(e$df, 7) + e0 <- estimate(g) + expect_null(e0$df) +}) diff --git a/tests/testthat/test-estimate_default.R b/tests/testthat/test-estimate_default.R index f82b4450..d6467d69 100644 --- a/tests/testthat/test-estimate_default.R +++ b/tests/testthat/test-estimate_default.R @@ -82,8 +82,8 @@ test_that("summary.estimate compared with estimate", { q <- compute_wald(B, coef(a3d), vcov(a3d), null) df <- attr(q, "df") e1 <- estimate(a3d, f=B) - e1b <- estimate(a3d, f=B, null=null) - e1a <- estimate(e1, f=diag(3), null=null) + e1b <- suppressWarnings(estimate(a3d, f=B, null=null)) + e1a <- suppressWarnings(estimate(e1, f=diag(3), null=null)) expect_equal(unname(df), unname(e1a$compare$parameter)) expect_equal(unname(df), unname(e1b$compare$parameter)) expect_true(abs(pchisq(q, df=df, lower.tail=FALSE) - e1a$compare$p.value)<1e-16) @@ -92,8 +92,8 @@ test_that("summary.estimate compared with estimate", { expect_true(abs(q - e1b$compare$statistic) < 1e-9) # function spec. e2 <- estimate(a3d, function(p) c(p[1]-p[2], p[2]-p[3], p[1]-p[3])) - e2b <- estimate(a3d, function(p) c(p[1]-p[2], p[2]-p[3], p[1]-p[3]), null=null) - e2a <- estimate(e2, f=diag(3), null=null) + e2b <- suppressWarnings(estimate(a3d, function(p) c(p[1]-p[2], p[2]-p[3], p[1]-p[3]), null=null)) + e2a <- suppressWarnings(estimate(e2, f=diag(3), null=null)) expect_equal(unname(df), unname(e2a$compare$parameter)) expect_equal(unname(df), unname(e2b$compare$parameter)) expect_true(abs(pchisq(q, df=df, lower.tail=FALSE) - e2a$compare$p.value)<1e-16) @@ -113,8 +113,8 @@ test_that("summary.estimate compared with estimate", { q <- compute_wald(B, coef(a3d), vcov(a3d), null) df <- attr(q, "df") e1 <- estimate(a3d, f=B) - e1b <- estimate(a3d, f=B, null=null) - e1a <- estimate(e1, f=diag(4), null=null) + e1b <- suppressWarnings(estimate(a3d, f=B, null=null)) + e1a <- suppressWarnings(estimate(e1, f=diag(4), null=null)) expect_equal(unname(df), unname(e1a$compare$parameter)) expect_equal(unname(df), unname(e1b$compare$parameter)) expect_true(abs(pchisq(q, df=df, lower.tail=FALSE) - e1a$compare$p.value)<1e-16) @@ -131,7 +131,7 @@ test_that("summary.estimate compared with estimate", { test_that("1D: Identity contrast, null = 0", { B <- matrix(1, nrow = 1, ncol = 1) null <- 0 - e <- estimate(a1, f = B, null = null) + e <- suppressWarnings(estimate(a1, f = B, null = null)) p <- coef(a1) S <- vcov(a1) expected <- compute_wald(B, p, S, null) @@ -141,7 +141,7 @@ test_that("1D: Identity contrast, null = 0", { test_that("1D: Identity contrast, null = 1", { B <- matrix(1, nrow = 1, ncol = 1) null <- 1 - e <- estimate(a1, f = B, null = null) + e <- suppressWarnings(estimate(a1, f = B, null = null)) p <- coef(a1) S <- vcov(a1) expected <- compute_wald(B, p, S, null) @@ -151,7 +151,7 @@ test_that("1D: Identity contrast, null = 1", { test_that("1D: Scalar contrast (scaling), null = 0", { B <- matrix(2, nrow = 1, ncol = 1) null <- 0 - e <- estimate(a1, f = B, null = null) + e <- suppressWarnings(estimate(a1, f = B, null = null)) p <- coef(a1) S <- vcov(a1) expected <- compute_wald(B, p, S, null) @@ -161,7 +161,7 @@ test_that("1D: Scalar contrast (scaling), null = 0", { test_that("2D: Identity contrast, null = c(0, 0)", { B <- diag(2) null <- c(0, 0) - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -171,7 +171,7 @@ test_that("2D: Identity contrast, null = c(0, 0)", { test_that("2D: Identity contrast, null equals true coefs", { B <- diag(2) null <- c(1, 2) - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -184,7 +184,7 @@ test_that("2D: Identity contrast, null equals true coefs", { test_that("2D: Identity contrast, null = c(3, 4)", { B <- diag(2) null <- c(3, 4) - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -194,7 +194,7 @@ test_that("2D: Identity contrast, null = c(3, 4)", { test_that("2D: Difference contrast c(1,-1), null = 0", { B <- matrix(c(1, -1), nrow = 1) null <- 0 - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -204,7 +204,7 @@ test_that("2D: Difference contrast c(1,-1), null = 0", { test_that("2D: Difference contrast c(1,-1), null = -1", { B <- matrix(c(1, -1), nrow = 1) null <- -1 # Testing H0: a1 - a2 = -1 (true, since 1 - 2 = -1) - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -216,7 +216,7 @@ test_that("2D: Difference contrast c(1,-1), null = -1", { test_that("2D: Sum contrast c(1,1), null = 3", { B <- matrix(c(1, 1), nrow = 1) null <- 3 # True sum = 1 + 2 = 3 - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -228,7 +228,7 @@ test_that("2D: Sum contrast c(1,1), null = 3", { test_that("2D: Scaling contrast, null = c(2, 6)", { B <- 2 * diag(2) null <- c(2, 6) # True: 2*c(1,2) = c(2,4), so null != true - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) p <- coef(a) S <- vcov(a) expected <- compute_wald(B, p, S, null) @@ -238,7 +238,7 @@ test_that("2D: Scaling contrast, null = c(2, 6)", { test_that("3D: Identity contrast, null = true coefs", { B <- diag(3) null <- c(1, 2, 3) - e <- estimate(a3d, f = B, null = null) + e <- suppressWarnings(estimate(a3d, f = B, null = null)) p <- coef(a3d) S <- vcov(a3d) expected <- compute_wald(B, p, S, null) @@ -249,7 +249,7 @@ test_that("3D: Identity contrast, null = true coefs", { test_that("3D: Identity contrast, null = c(0, 0, 0)", { B <- diag(3) null <- c(0, 0, 0) - e <- estimate(a3d, f = B, null = null) + e <- suppressWarnings(estimate(a3d, f = B, null = null)) p <- coef(a3d) S <- vcov(a3d) expected <- compute_wald(B, p, S, null) @@ -263,7 +263,7 @@ test_that("3D: Pairwise differences, null = c(-1, -2)", { 0, 1, -1 ), nrow = 2, byrow = TRUE) null <- c(-1, -1) - e <- estimate(a3d, f = B, null = null) + e <- suppressWarnings(estimate(a3d, f = B, null = null)) p <- coef(a3d) S <- vcov(a3d) expected <- compute_wald(B, p, S, null) @@ -275,7 +275,7 @@ test_that("3D: Single-row contrast, null = 6", { # H0: a1 + a2 + a3 = 6 (true: 1+2+3=6) B <- matrix(c(1, 1, 1), nrow = 1) null <- 6 - e <- estimate(a3d, f = B, null = null) + e <- suppressWarnings(estimate(a3d, f = B, null = null)) p <- coef(a3d) S <- vcov(a3d) expected <- compute_wald(B, p, S, null) @@ -286,7 +286,7 @@ test_that("3D: Single-row contrast, null = 6", { test_that("4D: Identity contrast, null equals true coefs", { B <- diag(4) null <- c(1, 2, 3, 4) - e <- estimate(a4d, f = B, null = null) + e <- suppressWarnings(estimate(a4d, f = B, null = null)) p <- coef(a4d) S <- vcov(a4d) expected <- compute_wald(B, p, S, null) @@ -301,7 +301,7 @@ test_that("4D: Non-square contrast (2x4), null = c(0, 0)", { 0, 0, 1, -1 ), nrow = 2, byrow = TRUE) null <- c(0, 0) - e <- estimate(a4d, f = B, null = null) + e <- suppressWarnings(estimate(a4d, f = B, null = null)) p <- coef(a4d) S <- vcov(a4d) expected <- compute_wald(B, p, S, null) @@ -315,7 +315,7 @@ test_that("4D: Non-square contrast, null equals true values", { 0, 0, 1, -1 ), nrow = 2, byrow = TRUE) null <- c(-1, -1) - e <- estimate(a4d, f = B, null = null) + e <- suppressWarnings(estimate(a4d, f = B, null = null)) p <- coef(a4d) S <- vcov(a4d) expected <- compute_wald(B, p, S, null) @@ -326,11 +326,11 @@ test_that("4D: Non-square contrast, null equals true values", { test_that("Degrees of freedom equals nrow(B)", { B <- diag(2) null <- c(0, 0) - e <- estimate(a, f = B, null = null) + e <- suppressWarnings(estimate(a, f = B, null = null)) expect_identical(unname(e$compare$parameter), nrow(B)) B <- matrix(c(1, -1, 0, 0, 1, -1), nrow = 2, byrow = TRUE) null <- c(0, 0) - e <- estimate(a3d, f = B, null = null) + e <- suppressWarnings(estimate(a3d, f = B, null = null)) expect_equal(unname(e$compare$parameter), nrow(B)) }) @@ -423,7 +423,7 @@ test_that("estimate.default parsedesign dispatch", { expect_equal(e$coefmat[, 1], 2.0) # Character contrasts combine with null hypothesis vector. - e <- estimate(a3d, "a1", "a1", null = c(0, 1)) + e <- suppressWarnings(estimate(a3d, "a1", "a1", null = c(0, 1))) expect_true(e$coefmat[1, "P-value"] != e$coefmat[2, "P-value"]) }) diff --git a/tests/testthat/test-summary-requires-IC.R b/tests/testthat/test-summary-requires-IC.R new file mode 100644 index 00000000..76ebf304 --- /dev/null +++ b/tests/testthat/test-summary-requires-IC.R @@ -0,0 +1,41 @@ +## summary.estimate() requires IC=TRUE on the input estimate when applying +## small-sample corrections (type / var.adj). Per Q5. + +context("summary.estimate IC requirement") + +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("summary(..., type=) errors when input estimate has IC=FALSE", { + g <- make_glm() + e <- estimate(g, IC = FALSE) + expect_error( + summary(e, type = "df"), + "IC=TRUE", + fixed = TRUE + ) +}) + +test_that("summary(..., var.adj=) errors when input estimate has IC=FALSE", { + g <- make_glm() + e <- estimate(g, IC = FALSE) + expect_error( + summary(e, type = "hc3", var.adj = 0.5), + "IC=TRUE", + fixed = TRUE + ) +}) + +test_that("summary(..., null=) works when input estimate has IC=FALSE", { + ## null/contrast do NOT require IC, only type/var.adj do. + g <- make_glm() + e <- estimate(g, IC = FALSE) + expect_silent(summary(e, null = 1)) +}) diff --git a/vignettes/influencefunction.Rmd b/vignettes/influencefunction.Rmd index 75dfd81e..c4126a55 100644 --- a/vignettes/influencefunction.Rmd +++ b/vignettes/influencefunction.Rmd @@ -254,7 +254,7 @@ coef(e) ## ## Asymptotic (robust) variance estimate vcov(e) ## Matrix with estimates and confidence limits -estimate(e, null=0, level = 0.99) |> parameter() +estimate(e, level = 0.99) |> parameter() ## Influence curve IC(e) |> head() ## Join estimates @@ -445,7 +445,7 @@ As we have access to the joint asymptotic distribution we can for example test for whether the odds-ratio is the same for the two responses: ```{r hypo1} -estimate(e, cbind(0,1,0,-1), null=0) +summary(estimate(e), contrast=cbind(0,1,0,-1), null=0) ``` More details an be found in the Section on [hypothesis testing](#linear-contrasts-and-hypothesis-testing). @@ -855,7 +855,7 @@ estimate(gg, B) The $\bm{b}_0$ vector (default assumed to be zero) can be specified via the `null` argument ```{r estcontrast1} -estimate(gg, B, null=1) +summary(estimate(gg, B), null=1) ``` For testing multiple hypotheses we use that @@ -873,7 +873,7 @@ estimate(gg, B) Such linear statistics can also be specified directly as expressions of the parameter names ```{r} -estimate(gg, a + a.1, 2*a - a.2, a, null=c(2,1,1)) +summary(estimate(gg, a + a.1, 2*a - a.2, a), null=c(2,1,1)) ``` We refer to the function `lava::contr` and `lava::parsedesigns` for defining contrast matrices. @@ -885,7 +885,7 @@ comparisons of different exposure estimates: ```{r pairwise.diff} pairwise.diff(3) -estimate(gg, pairwise.diff(3), null=c(1,1,1), use=c(2,4,6)) +summary(estimate(gg, pairwise.diff(3), use=c(2,4,6)), null=c(1,1,1)) ``` When conducting multiple tests each at a nominal-level of $\alpha$ the overall @@ -913,7 +913,7 @@ in [@RitzPipper_multcomp] the joint distribution of $Z_{1},\ldots,Z_{p}$ can be estimated from the IFs. This is implemented in the `alpha_zmax` method ```{r pcorrect, eval=has_mets } -gg0 <- estimate(gg, use="^a", regex=TRUE, null=rep(.8, 3)) +gg0 <- suppressWarnings(estimate(gg, use="^a", regex=TRUE, null=rep(.8, 3))) alpha_zmax(gg0) ```