# anova() on nested admFits: the ORDINARY likelihood-ratio test. # # Everything here is about the arithmetic and the refusals. There is no # reweighting to test -- that is the point of the plain form -- so what is left # is that the test statistic is the objective difference, that Df counts the # added parameters, and that the comparisons which are NOT likelihood ratios # are refused rather than reported. # A minimal stand-in for an admFit: anova() reads only `objective`, `env$method`, # `env$nNodes`, `env$nSim`, `env$strataNodes`, `env$AIC`/`env$BIC` and the # parameter names off admExtra. .lrt_fit <- function(par_names, objective, method = "adgh", nNodes = 5L, nSim = NULL, strataNodes = NULL, AIC = NA_real_, BIC = NA_real_) { e <- new.env(parent = emptyenv()) e$admExtra <- list(par_names = par_names) e$method <- method e$nNodes <- nNodes e$nSim <- nSim e$strataNodes <- strataNodes e$AIC <- AIC e$BIC <- BIC structure(list(env = e, objective = objective), class = "admFit") } test_that("the statistic is the objective difference on the added parameters", { full <- .lrt_fit(c("tcl", "tv", "b1"), 100) red <- .lrt_fit(c("tcl", "tv"), 106) a <- anova(full, red) expect_s3_class(a, "anova.admFit") # smallest model first, the way anova.lm orders its table expect_identical(a$Npar, c(2L, 3L)) expect_equal(a$dOFV[2L], 6) expect_identical(a$Df[2L], 1L) expect_equal(a$p[2L], stats::pchisq(6, 1, lower.tail = FALSE)) # the first row is the reference model and carries no test expect_true(is.na(a$dOFV[1L])) expect_true(is.na(a$p[1L])) }) test_that("Df counts EVERY added parameter, not just one", { full <- .lrt_fit(c("tcl", "tv", "b1", "b2", "b3"), 40) red <- .lrt_fit(c("tcl", "tv"), 55) a <- anova(full, red) expect_identical(a$Df[2L], 3L) expect_equal(a$p[2L], stats::pchisq(15, 3, lower.tail = FALSE)) }) test_that("a NEGATIVE dOFV is reported, not clamped to zero", { # The larger model cannot fit worse at its own optimum, so this says one of # the two did not converge. Rounding it to zero would hide that. full <- .lrt_fit(c("tcl", "tv", "b1"), 110) red <- .lrt_fit(c("tcl", "tv"), 106) a <- anova(full, red) expect_lt(a$dOFV[2L], 0) expect_true(is.na(a$p[2L])) }) test_that("a non-nested pair is refused, not approximated around", { a <- .lrt_fit(c("tcl", "tv", "b1"), 100) b <- .lrt_fit(c("tcl", "tv", "b2"), 101) expect_error(anova(a, b), "not nested") # and so is a pair with nothing to test c1 <- .lrt_fit(c("tcl", "tv"), 100) c2 <- .lrt_fit(c("tcl", "tv"), 101) expect_error(anova(c1, c2), "nothing to test") }) test_that("fits from DIFFERENT estimators are refused", { # Each scores its own approximation to the same likelihood -- FO-linearised, # quadrature, Monte Carlo -- so the difference of two objectives is not a # likelihood ratio. This returned a perfectly finite p before it was caught. full <- .lrt_fit(c("tcl", "tv", "b1"), 100, method = "adgh") red <- .lrt_fit(c("tcl", "tv"), 106, method = "adfo") expect_error(anova(full, red), "different estimators") }) test_that("fits on different node counts are refused", { # The objective moves with the grid, which is why .adghGrid refuses to change # the point count mid-fit in the first place. full <- .lrt_fit(c("tcl", "tv", "b1"), 100, nNodes = 9L) red <- .lrt_fit(c("tcl", "tv"), 106, nNodes = 5L) expect_error(anova(full, red), "node counts") }) test_that("admc/adirmc fits on different n_sim are refused", { # nNodes is NULL for these two estimators, so the grid check above is a # no-op -- n_sim is their equivalent moving part and needs its own guard. full <- .lrt_fit(c("tcl", "tv", "b1"), 100, method = "admc", nNodes = NULL, nSim = 5000L) red <- .lrt_fit(c("tcl", "tv"), 106, method = "admc", nNodes = NULL, nSim = 2000L) expect_error(anova(full, red), "n_sim") }) test_that("fits at different stratum resolutions are refused", { # The unnamed form, as a fit saved before the stamp carried covariate names. full <- .lrt_fit(c("tcl", "tv", "b1"), 100, strataNodes = 9L) red <- .lrt_fit(c("tcl", "tv"), 106, strataNodes = 5L) expect_error(anova(full, red), "stratum resolutions") # AN UNSET STAMP IS NOT A DISAGREEMENT, whether or not the other side is # named. This used to be refused here and accepted two tests down, where the # named side made the per-covariate branch fire and the intersection with an # unnamed stamp came back empty -- so the same pair was judged on whether # the covariate happened to have a name. expect_s3_class(anova(full, .lrt_fit(c("tcl", "tv"), 106)), "anova.admFit") # A fit saved BEFORE the stamp carried names still gets compared. Entering # the per-covariate branch on one named side alone, the intersection was # empty and this passed at genuinely different resolutions -- the one pair # the unnamed fallback exists to catch. expect_error(anova(full, .lrt_fit(c("tcl", "tv"), 106, strataNodes = c(WT = 5L))), "stratum resolutions") # ...and the same stamp on either side of the rename is still the same # resolution, so an integer `9` and the character `"9"` agree. expect_s3_class(anova(full, .lrt_fit(c("tcl", "tv"), 106, strataNodes = "9")), "anova.admFit") }) test_that("the resolution is compared per covariate, where the fits overlap", { # THE STAMP IS PER COVARIATE THE MODEL READS. One number per fit refused # exactly the nested pair a covariate test is made of: the null model has # dropped the term, so its sources are not cut into nodes along it and it # records nothing -- and its objective cannot depend on a resolution it # never used. Measured on a real pair, the narrowed null fit (2 studies) # and the same null refitted at the full resolution (10 studies) give dOFV # 409.162568 and 409.162615: the same test to five decimals. full <- .lrt_fit(c("tcl", "tv", "bcrcl"), 100, strataNodes = c(CRCL = 5L)) # The null model reads no covariate at all, so it stamps nothing. expect_s3_class(anova(full, .lrt_fit(c("tcl", "tv"), 106)), "anova.admFit") # Nor does a covariate only ONE of them reads put them on different scales. expect_s3_class(anova(full, .lrt_fit(c("tcl", "tv"), 106, strataNodes = c(WT = 9L))), "anova.admFit") # REFUSED where it matters. Both read CRCL and cut it differently: expect_error(anova(full, .lrt_fit(c("tcl", "tv"), 106, strataNodes = c(CRCL = 9L))), "stratum resolutions on 'CRCL'") # ...and a covariate the model READS but no source was cut on records "0", # because integrating over the whole declared distribution and scoring a set # of nodes are different data to a model that can see the difference. expect_error(anova(full, .lrt_fit(c("tcl", "tv"), 106, strataNodes = c(CRCL = "0"))), "stratum resolutions on 'CRCL'") # "0" IS NOT "1". A source cut at `strata_nodes = 1L` is pinned at the latent # median, which is a third dataset again -- and stamping both as 1 made this # pair pass. expect_error(anova(.lrt_fit(c("tcl", "tv", "bcrcl"), 100, strataNodes = c(CRCL = "1")), .lrt_fit(c("tcl", "tv"), 106, strataNodes = c(CRCL = "0"))), "stratum resolutions on 'CRCL'") # THE WHOLE MULTISET, not its maximum: two sources cut at three and nine # nodes are not the same data as both cut at nine, and max() reported both # as 9. expect_error(anova(.lrt_fit(c("tcl", "tv", "bcrcl"), 100, strataNodes = c(CRCL = "3/9")), .lrt_fit(c("tcl", "tv"), 106, strataNodes = c(CRCL = "9"))), "stratum resolutions on 'CRCL'") }) test_that("anova() needs a pair, and needs admFits", { full <- .lrt_fit(c("tcl", "tv", "b1"), 100) expect_error(anova(full), "nothing to compare") expect_error(anova(full, 1), "must be an admFit") }) test_that("printing keeps `Test` as written", { # stats::print.anova routes through printCoefmat, which calls data.matrix() # -- and that turns a CHARACTER column into its factor CODES, so "1 vs 2" # would print as 1. full <- .lrt_fit(c("tcl", "tv", "b1"), 100) red <- .lrt_fit(c("tcl", "tv"), 106) out <- utils::capture.output(print(anova(full, red))) expect_true(any(grepl("Likelihood-ratio test", out))) expect_true(any(grepl("1 vs 2", out, fixed = TRUE))) })