diff --git a/NAMESPACE b/NAMESPACE index 0afd85824e..f7cd7b31b5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -248,6 +248,7 @@ export(labels_use_control) export(level_order) export(logistic_regression_cols) export(logistic_summary_by_flag) +export(mantel_fleiss_crit) export(month2day) export(or_clogit) export(or_glm) diff --git a/NEWS.md b/NEWS.md index 14a87ad62f..b702fc5747 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,9 @@ # tern 0.9.11.9000 ### Enhancements +* Added `mantel_fleiss_crit()` to check the Mantel-Fleiss criterion + for stratified 2 x 2 contingency tables, together with a vignette + demonstrating its use. (#1512) * Added a `method_only` argument to `d_proportion()`, `d_proportion_diff()`, and `d_test_proportion_diff()` to allow returning method labels without an additional description. (#1525) @@ -10,6 +13,10 @@ stored in a single column. (#1499) * Added the `exclude_rows` argument to `g_forest()` to allow excluding selected rows from the forest plot before plotting. (#1498) + +# tern 0.9.11 + +### Enhancements * Added `factor_level_method` argument to `df_explicit_na()` to control factor level ordering when converting character or logical columns. Supported methods: `"sort_auto"` (default, locale-aware, preserves original behavior), `"sort_radix"` (byte-order / ASCII sort), and diff --git a/R/prop_diff_test.R b/R/prop_diff_test.R index 5b744e9963..b798f072c2 100644 --- a/R/prop_diff_test.R +++ b/R/prop_diff_test.R @@ -500,3 +500,160 @@ prop_fisher <- function(tbl, alternative = c("two.sided", "less", "greater")) { tbl <- tbl[, c("TRUE", "FALSE")] stats::fisher.test(tbl, alternative = alternative)$p.value } + +#' Check the Mantel-Fleiss Criterion +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' Checks the Mantel-Fleiss criterion for stratified 2 x 2 contingency tables. +#' +#' @details +#' The Mantel-Fleiss statistic is calculated as +#' +#' \deqn{ +#' MF = \min\left( +#' [\sum_h m_{11h} - \sum_h {(n_{11h})}_L],\ +#' [\sum_h {(n_{11h})}_U - \sum_h m_{11h}] +#' \right), +#' } +#' +#' where \eqn{h} indexes the non-empty strata. For each stratum \eqn{h}, the +#' expected frequency of cell \eqn{(1, 1)} in table \eqn{h}, under the +#' hypothesis of no association between group and response, is +#' +#' \deqn{ +#' m_{11h} = \frac{n_{1.h} n_{.1h}}{n_h}. +#' } +#' +#' The lower and upper bounds for \eqn{n_{11h}}, given the marginal totals, +#' are: +#' +#' \deqn{ +#' {(n_{11h})}_L = \max(0, n_{1.h} - n_{.2h}), +#' } +#' \deqn{ +#' {(n_{11h})}_U = \min(n_{.1h}, n_{1.h}). +#' } +#' +#' The Mantel-Fleiss criterion is satisfied when \eqn{MF \ge} `threshold`. +#' By default, `threshold = 5`, corresponding to the criterion described +#' by Mantel and Fleiss (1980). +#' +#' Strata with all cell counts equal to zero are excluded from the +#' calculation. If all strata contain zero observations, there are no +#' non-empty strata over which to calculate the Mantel-Fleiss statistic, and +#' the statistic is therefore undefined. In this case, the function returns +#' `NA`. +#' +#' @param tbl (`array`)\cr +#' A three-dimensional contingency table containing the counts for each +#' combination of group, response, and stratum. The first two dimensions +#' must correspond to the two variables defining the 2 x 2 contingency +#' table (group and response), in either order. The third dimension must +#' correspond to the strata. The first two dimensions must each have exactly +#' two levels. All cell values must be finite, non-missing integer counts. +#' @param include_value (`logical(1)`)\cr +#' Whether to include the calculated Mantel-Fleiss statistic as an attribute +#' of the result. +#' @param threshold (`numeric(1)`)\cr +#' The minimum Mantel-Fleiss statistic required for the criterion to be +#' considered satisfied. +#' +#' @return A logical value indicating whether the Mantel-Fleiss criterion +#' is satisfied. If `include_value = TRUE`, the result also contains a +#' value attribute with the calculated Mantel-Fleiss statistic. If there +#' are no non-empty strata, the result is `NA` and the value attribute is +#' `NA_real_`. +#' +#' @examples +#' set.seed(123) +#' n <- 40 +#' +#' grp <- factor(sample(c("Active", "Control"), n, replace = TRUE)) +#' rsp <- sample(c(TRUE, FALSE), n, replace = TRUE) +#' strata1 <- factor(sample(c("A", "B"), n, replace = TRUE)) +#' strata2 <- factor(sample(c("x", "y"), n, replace = TRUE)) +#' strata <- interaction(strata1, strata2) +#' +#' tbl <- table(grp, rsp, strata) +#' tbl +#' +#' is_mf_satisfied <- mantel_fleiss_crit(tbl) +#' is_mf_satisfied +#' mantel_fleiss_crit(tbl, include_value = TRUE) +#' +#' # Examples of use. +#' +#' if (is_mf_satisfied) { +#' print("CMH") +#' prop_diff_cmh(rsp, grp, strata)$prop +#' } else { +#' print("Exact") +#' prop_diff_uncond_exact(rsp, grp)$prop +#' } +#' +#' if (is_mf_satisfied) { +#' print("CMH") +#' prop_cmh(tbl) +#' } else { +#' print("Exact") +#' prop_fisher(table(grp, rsp)) +#' } +#' +#' @references +#' Mantel, N., and Fleiss, J. L. (1980). +#' Minimum Expected Cell Size Requirements for the Mantel-Haenszel +#' One-Degree-of-Freedom Chi-Square Test and a Related Rapid Procedure. +#' \emph{American Journal of Epidemiology}, 112(1), 129--134. +#' +#' @export +mantel_fleiss_crit <- function(tbl, include_value = FALSE, threshold = 5L) { + checkmate::assert_array(tbl, mode = "integerish", any.missing = FALSE, d = 3L) + checkmate::assert_true(all(tbl >= 0L)) + checkmate::assert_true(all(is.finite(tbl))) + checkmate::assert_true(nrow(tbl) == 2L) + checkmate::assert_true(ncol(tbl) == 2L) + checkmate::assert_flag(include_value) + checkmate::assert_number(threshold) + + # Drop strata with no observations. + tbl <- tbl[, , apply(tbl, 3L, sum) > 0, drop = FALSE] + + # Add marginal totals over the group and response dimensions, + # retaining the stratum dimension. + tbl_mrgn <- stats::addmargins(tbl, margin = 1:2) + + # If there are no non-empty strata, the Mantel-Fleiss criterion is undefined + # because there are no strata over which to calculate it. + if (dim(tbl)[3L] == 0L) { + is_satisfied <- NA + if (include_value) { + attr(is_satisfied, "value") <- NA_real_ + } + return(is_satisfied) + } + + n_1dot <- tbl_mrgn[1L, "Sum", ] + n_dot1 <- tbl_mrgn["Sum", 1L, ] + n_dot2 <- tbl_mrgn["Sum", 2L, ] + n <- tbl_mrgn["Sum", "Sum", ] + + # Expected value of n_11 under the hypothesis of no association + # between group and response (within a given stratum). + m_11 <- (n_1dot * n_dot1) / n + # Lower and upper bounds for n_11 given the marginal totals (within a given stratum). + n_11_lwr <- pmax(0L, n_1dot - n_dot2) + n_11_upr <- pmin(n_dot1, n_1dot) + + mf_value <- min( + sum(m_11) - sum(n_11_lwr), + sum(n_11_upr) - sum(m_11) + ) + + is_satisfied <- mf_value >= threshold + if (include_value) { + attr(is_satisfied, "value") <- mf_value + } + + is_satisfied +} diff --git a/_pkgdown.yml b/_pkgdown.yml index d6e0e13c6f..06ce089982 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -113,6 +113,7 @@ reference: - -h_xticks - -prop_diff - check_diff_prop_ci + - mantel_fleiss_crit - title: rtables Helper Functions desc: These functions help to work with the `rtables` package and may be diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 957ec82ad2..11e0d09a00 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -115,6 +115,16 @@ @article{MantelHaenszel1959 year = {1959} } +@article{MantelFleiss1980, + title = {Minimum Expected Cell Size Requirements for the Mantel-Haenszel One-Degree-of-Freedom Chi-Square Test and a Related Rapid Procedure}, + author = {Mantel, N. and Fleiss, J. L.}, + journal = {American Journal of Epidemiology}, + volume = {112}, + number = {1}, + pages = {129--134}, + year = {1980} +} + @article{Sato1989, title = {On the variance estimator for the Mantel-Haenszel Risk Difference}, author = {Sato, Tosiya and Greenland, Sander and Robins, James M.}, diff --git a/inst/WORDLIST b/inst/WORDLIST index fd4ab057f3..fb9ae96457 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -14,7 +14,9 @@ CDISC CMH CQ Clopper +Cochran Coull +Fleiss Haenszel Hauck Hilferty @@ -22,6 +24,7 @@ Hoffmann Jeffreys Kaplan Kenward +MF MMRM MedDRA Miettinen diff --git a/man/mantel_fleiss_crit.Rd b/man/mantel_fleiss_crit.Rd new file mode 100644 index 0000000000..ccb7d5cbf2 --- /dev/null +++ b/man/mantel_fleiss_crit.Rd @@ -0,0 +1,117 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prop_diff_test.R +\name{mantel_fleiss_crit} +\alias{mantel_fleiss_crit} +\title{Check the Mantel-Fleiss Criterion} +\usage{ +mantel_fleiss_crit(tbl, include_value = FALSE, threshold = 5L) +} +\arguments{ +\item{tbl}{(\code{array})\cr +A three-dimensional contingency table containing the counts for each +combination of group, response, and stratum. The first two dimensions +must correspond to the two variables defining the 2 x 2 contingency +table (group and response), in either order. The third dimension must +correspond to the strata. The first two dimensions must each have exactly +two levels. All cell values must be finite, non-missing integer counts.} + +\item{include_value}{(\code{logical(1)})\cr +Whether to include the calculated Mantel-Fleiss statistic as an attribute +of the result.} + +\item{threshold}{(\code{numeric(1)})\cr +The minimum Mantel-Fleiss statistic required for the criterion to be +considered satisfied.} +} +\value{ +A logical value indicating whether the Mantel-Fleiss criterion +is satisfied. If \code{include_value = TRUE}, the result also contains a +value attribute with the calculated Mantel-Fleiss statistic. If there +are no non-empty strata, the result is \code{NA} and the value attribute is +\code{NA_real_}. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +Checks the Mantel-Fleiss criterion for stratified 2 x 2 contingency tables. +} +\details{ +The Mantel-Fleiss statistic is calculated as + +\deqn{ +MF = \min\left( +[\sum_h m_{11h} - \sum_h {(n_{11h})}_L],\ +[\sum_h {(n_{11h})}_U - \sum_h m_{11h}] +\right), +} + +where \eqn{h} indexes the non-empty strata. For each stratum \eqn{h}, the +expected frequency of cell \eqn{(1, 1)} in table \eqn{h}, under the +hypothesis of no association between group and response, is + +\deqn{ +m_{11h} = \frac{n_{1.h} n_{.1h}}{n_h}. +} + +The lower and upper bounds for \eqn{n_{11h}}, given the marginal totals, +are: + +\deqn{ +{(n_{11h})}_L = \max(0, n_{1.h} - n_{.2h}), +} +\deqn{ +{(n_{11h})}_U = \min(n_{.1h}, n_{1.h}). +} + +The Mantel-Fleiss criterion is satisfied when \eqn{MF \ge} \code{threshold}. +By default, \code{threshold = 5}, corresponding to the criterion described +by Mantel and Fleiss (1980). + +Strata with all cell counts equal to zero are excluded from the +calculation. If all strata contain zero observations, there are no +non-empty strata over which to calculate the Mantel-Fleiss statistic, and +the statistic is therefore undefined. In this case, the function returns +\code{NA}. +} +\examples{ +set.seed(123) +n <- 40 + +grp <- factor(sample(c("Active", "Control"), n, replace = TRUE)) +rsp <- sample(c(TRUE, FALSE), n, replace = TRUE) +strata1 <- factor(sample(c("A", "B"), n, replace = TRUE)) +strata2 <- factor(sample(c("x", "y"), n, replace = TRUE)) +strata <- interaction(strata1, strata2) + +tbl <- table(grp, rsp, strata) +tbl + +is_mf_satisfied <- mantel_fleiss_crit(tbl) +is_mf_satisfied +mantel_fleiss_crit(tbl, include_value = TRUE) + +# Examples of use. + +if (is_mf_satisfied) { + print("CMH") + prop_diff_cmh(rsp, grp, strata)$prop +} else { + print("Exact") + prop_diff_uncond_exact(rsp, grp)$prop +} + +if (is_mf_satisfied) { + print("CMH") + prop_cmh(tbl) +} else { + print("Exact") + prop_fisher(table(grp, rsp)) +} + +} +\references{ +Mantel, N., and Fleiss, J. L. (1980). +Minimum Expected Cell Size Requirements for the Mantel-Haenszel +One-Degree-of-Freedom Chi-Square Test and a Related Rapid Procedure. +\emph{American Journal of Epidemiology}, 112(1), 129--134. +} diff --git a/tests/testthat/test-mantel_fleiss_crit.R b/tests/testthat/test-mantel_fleiss_crit.R new file mode 100644 index 0000000000..93b07df21b --- /dev/null +++ b/tests/testthat/test-mantel_fleiss_crit.R @@ -0,0 +1,189 @@ +test_that("mantel_fleiss_crit() works with multiple observations and strata", { + tbl <- array( + c(9L, 8L, 6L, 9L, 6L, 5L, 8L, 5L, 5L, 5L, 5L, 4L, 11L, 5L, 7L, 2L), + dim = c(2L, 2L, 4L) + ) + + expect_silent( + result <- mantel_fleiss_crit(tbl) + ) + expect_silent( + result_val <- mantel_fleiss_crit(tbl, TRUE) + ) + + expect_identical(result, TRUE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_equal(attributes(result_val), list(value = 20.16857), tolerance = 1e-6) +}) + +test_that("mantel_fleiss_crit() works with small stratified data and dimnames", { + tbl <- array( + c(2L, 2L, 1L, 2L, 0L, 1L, 2L, 1L, 1L, 1L, 3L, 1L, 1L, 0L, 1L, 1L), + dim = c(2L, 2L, 4L), + dimnames = list(grp = c("Gr1", "Gr2"), rsp = c("T", "F"), strata = LETTERS[1:4]) + ) + + expect_silent( + result <- mantel_fleiss_crit(tbl) + ) + expect_silent( + result_val <- mantel_fleiss_crit(tbl, TRUE) + ) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_equal(attributes(result_val), list(value = 2.785714), tolerance = 1e-6) +}) + +test_that("mantel_fleiss_crit() works with 1 stratum", { + tbl <- array(c(4L, 4L, 7L, 5L), dim = c(2L, 2L, 1L)) + + expect_silent( + result <- mantel_fleiss_crit(tbl) + ) + expect_silent( + result_val <- mantel_fleiss_crit(tbl, TRUE) + ) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_equal(attributes(result_val), list(value = 3.6), tolerance = 1e-6) +}) + +test_that("mantel_fleiss_crit() ignores unobserved strata levels", { + tbl <- array( + c(1L, 4L, 3L, 3L, 5L, 0L, 7L, 4L, 0L, 0L, 0L, 0L, 3L, 1L, 0L, 6L), + dim = c(2L, 2L, 4L) + ) + + expect_silent( + result <- mantel_fleiss_crit(tbl) + ) + expect_silent( + result_val <- mantel_fleiss_crit(tbl, TRUE) + ) + + expect_identical(result, TRUE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_equal(attributes(result_val), list(value = 5.231818), tolerance = 1e-6) +}) + +test_that("mantel_fleiss_crit() returns NA when all cell counts equal zero", { + tbl <- array(rep(0L, 16L), dim = c(2L, 2L, 4L)) + + expect_silent( + result <- mantel_fleiss_crit(tbl) + ) + expect_silent( + result_val <- mantel_fleiss_crit(tbl, TRUE) + ) + + expect_identical(result, NA) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_identical(attributes(result_val), list(value = NA_real_)) +}) + +test_that("mantel_fleiss_crit() handles a stratum with observations in one cell only", { + tbl <- array(c(1L, 1L, 1L, 1L, 0L, 4L, 0L, 0L), dim = c(2L, 2L, 2L)) + + result <- mantel_fleiss_crit(tbl) + result_val <- mantel_fleiss_crit(tbl, TRUE) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_equal(attributes(result_val), list(value = 1), tolerance = 1e-6) +}) + +test_that("mantel_fleiss_crit() includes the MF = 5 boundary", { + tbl <- array(c(5L, 5L, 10L, 10L), dim = c(2L, 2L, 1L)) + + result <- mantel_fleiss_crit(tbl, include_value = TRUE) + + expect_identical(result, TRUE, ignore_attr = TRUE) + expect_identical(attributes(result), list(value = 5)) +}) + +test_that("mantel_fleiss_crit() handles data with no non-responses", { + tbl <- array(c(2L, 4L, 0L, 0L), dim = c(2L, 2L, 2L)) + + result <- mantel_fleiss_crit(tbl) + result_val <- mantel_fleiss_crit(tbl, TRUE) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_identical(attributes(result_val), list(value = 0)) +}) + +test_that("mantel_fleiss_crit() handles data with no responses", { + tbl <- array(c(0L, 0L, 4L, 0L), dim = c(2L, 2L, 2L)) + + result <- mantel_fleiss_crit(tbl) + result_val <- mantel_fleiss_crit(tbl, TRUE) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_identical(attributes(result_val), list(value = 0)) +}) + +test_that("mantel_fleiss_crit() works with observations from one group only (Gr1)", { + tbl <- array(c(46L, 0L, 4L, 0L), dim = c(2L, 2L, 2L)) + + result <- mantel_fleiss_crit(tbl) + result_val <- mantel_fleiss_crit(tbl, TRUE) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_identical(attributes(result_val), list(value = 0)) +}) + +test_that("mantel_fleiss_crit() works with observations from one group only (Gr2)", { + tbl <- array(c(0L, 30L, 0L, 3L), dim = c(2L, 2L, 2L)) + + result <- mantel_fleiss_crit(tbl) + result_val <- mantel_fleiss_crit(tbl, TRUE) + + expect_identical(result, FALSE) + expect_identical(result_val, result, ignore_attr = TRUE) + expect_identical(attributes(result_val), list(value = 0)) +}) + +test_that("mantel_fleiss_crit() respects the threshold", { + tbl <- array(c(5L, 5L, 10L, 10L), dim = c(2L, 2L, 1L)) + + expect_identical(mantel_fleiss_crit(tbl, threshold = 5), TRUE) + expect_identical(mantel_fleiss_crit(tbl, threshold = 6), FALSE) + + result <- mantel_fleiss_crit(tbl, include_value = TRUE, threshold = 6) + expect_identical(result, FALSE, ignore_attr = TRUE) + expect_identical(attributes(result), list(value = 5)) +}) + +test_that("mantel_fleiss_crit() validates inputs", { + # tbl + expect_error(mantel_fleiss_crit(matrix(1L, nrow = 2, ncol = 2))) + expect_error(mantel_fleiss_crit(array(1L, dim = c(2L, 2L, 2L, 1L)))) + expect_error(mantel_fleiss_crit(array(1L, dim = c(3L, 2L, 2L)))) + expect_error(mantel_fleiss_crit(array(1L, dim = c(2L, 3L, 2L)))) + + # Missing / invalid values. + dim3d <- c(2L, 2L, 2L) + expect_error(mantel_fleiss_crit(array(NA_integer_, dim = dim3d))) + expect_error(mantel_fleiss_crit(array("1", dim = dim3d))) + expect_error(mantel_fleiss_crit(array(NA_real_, dim = dim3d))) + expect_error(mantel_fleiss_crit(array(NaN, dim = dim3d))) + expect_error(mantel_fleiss_crit(array(-1, dim = dim3d))) + expect_error(mantel_fleiss_crit(array(-1L, dim = dim3d))) + expect_error(mantel_fleiss_crit(array(Inf, dim = dim3d))) + + # include_value + tbl <- array(1L, dim = c(2L, 2L, 3L)) + expect_error(mantel_fleiss_crit(tbl, include_value = c(TRUE, FALSE))) + expect_error(mantel_fleiss_crit(tbl, include_value = 1L)) + expect_error(mantel_fleiss_crit(tbl, include_value = 1)) + + # threshold + expect_error(mantel_fleiss_crit(tbl, threshold = c(5, 10))) + expect_error(mantel_fleiss_crit(tbl, threshold = "5")) + expect_error(mantel_fleiss_crit(tbl, threshold = NA_real_)) + expect_error(mantel_fleiss_crit(tbl, threshold = NaN)) +}) diff --git a/vignettes/mantel_fleiss_criterion.Rmd b/vignettes/mantel_fleiss_criterion.Rmd new file mode 100644 index 0000000000..a374f5e1b7 --- /dev/null +++ b/vignettes/mantel_fleiss_criterion.Rmd @@ -0,0 +1,188 @@ +--- +title: "The Mantel-Fleiss Criterion" +date: "2026-09-09" +output: + rmarkdown::html_document: + theme: "spacelab" + highlight: "kate" + toc: true + toc_float: true +vignette: > + %\VignetteIndexEntry{The Mantel-Fleiss Criterion} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +bibliography: ../inst/REFERENCES.bib +editor_options: + markdown: + wrap: 72 +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(tern) +``` + +## Introduction + +When comparing a binary response between two groups while adjusting for a +stratification variable, the Cochran-Mantel-Haenszel (CMH) test is a common +choice. Like other large-sample procedures, the CMH test relies on an asymptotic +(chi-square) approximation, which can be unreliable when the stratified +$2 \times 2$ tables are sparse. In those situations an exact method is +preferable. + +The **Mantel-Fleiss criterion** [@MantelFleiss1980] is a simple, quick check +that tells you whether the sample is large enough for the asymptotic CMH +approximation to be trustworthy. The `mantel_fleiss_crit()` function in `tern` +evaluates this criterion for a stratified $2 \times 2$ contingency table and +returns whether it is satisfied. You can use the result to decide, in a data +driven way, whether to run the CMH test or fall back to an exact procedure. + +`mantel_fleiss_crit()` is a standalone utility: it does not perform any test +itself. It is meant to be used alongside the proportion functions in `tern` +(such as `prop_diff_cmh()`, `prop_cmh()`, `prop_diff_uncond_exact()`, and +`prop_fisher()`) when writing custom analysis functions. + +## The criterion + +Consider a stratified $2 \times 2$ table where $h$ indexes the strata. Within +stratum $h$, write the cell and margin counts as + +$$ +\begin{array}{c|cc|c} +& \text{Response} & \text{No response} & \text{Row total} \\ +\hline +\text{Group 1} & n_{11h} & n_{12h} & n_{1 \cdot h} \\ +\text{Group 2} & n_{21h} & n_{22h} & n_{2 \cdot h} \\ +\hline +\text{Column total} & n_{\cdot 1 h} & n_{\cdot 2 h} & n_{h} +\end{array} +$$ + +Under the hypothesis of no association between group and response, the expected +count in cell $(1, 1)$ of stratum $h$ is + +$$ +m_{11h} = \frac{n_{1 \cdot h}\, n_{\,\cdot 1 h}}{n_{h}} . +$$ + +Given the fixed margins, the observed count $n_{11h}$ can range between the +bounds + +$$ +(n_{11h})_L = \max(0,\ n_{1 \cdot h} - n_{\, \cdot 2 h}), \qquad +(n_{11h})_U = \min(n_{\, \cdot 1 h},\ n_{1 \cdot h}) . +$$ + +The Mantel-Fleiss statistic aggregates these quantities across the non-empty +strata: + +$$ +MF = \min \left( +\left[ \sum_h m_{11h} - \sum_h (n_{11h})_L \right],\ +\left[ \sum_h (n_{11h})_U - \sum_h m_{11h} \right] +\right) . +$$ + +The criterion is considered **satisfied** when $MF \ge$ `threshold`. The +default `threshold = 5` corresponds to the rule proposed by @MantelFleiss1980: +when $MF \ge 5$, the asymptotic CMH approximation is generally adequate. + +## Basic usage + +`mantel_fleiss_crit()` expects a three-dimensional contingency table (an +`array`) whose first two dimensions are the group and response (each with two +levels, in either order) and whose third dimension is the stratum. + +```{r} +set.seed(123) +n <- 80 + +grp <- factor(sample(c("Active", "Control"), n, replace = TRUE)) +rsp <- sample(c(TRUE, FALSE), n, replace = TRUE) +strata1 <- factor(sample(c("A", "B"), n, replace = TRUE)) +strata2 <- factor(sample(c("x", "y"), n, replace = TRUE)) +strata <- interaction(strata1, strata2) + +tbl <- table(grp, rsp, strata) +tbl +``` + +Passing the table to `mantel_fleiss_crit()` returns a single logical value: + +```{r} +mantel_fleiss_crit(tbl) +``` + +To see the underlying value of the `MF` statistic, set `include_value = TRUE`. +The Mantel-Fleiss value is then attached to the result as a `"value"` attribute: + +```{r} +mantel_fleiss_crit(tbl, include_value = TRUE) +``` + +The `threshold` argument controls how large the statistic must be for the +criterion to hold. Raising it makes the criterion more conservative: + +```{r} +mantel_fleiss_crit(tbl, threshold = 15, include_value = TRUE) +``` + +## Choosing a test based on the criterion + +The typical use case is to branch between an asymptotic and an exact method +depending on whether the criterion is satisfied. The example below estimates +the stratified difference in proportions with the CMH method when the criterion +holds, and with the unconditional exact method otherwise: + +```{r} +is_mf_satisfied <- mantel_fleiss_crit(tbl) + +if (is_mf_satisfied) { + # Large enough sample: use the asymptotic CMH estimate. + prop_diff_cmh(rsp, grp, strata)$diff +} else { + # Sparse data: fall back to the exact (unstratified) method. + prop_diff_uncond_exact(rsp, grp)$diff +} +``` + +The same idea can be used to select a test statistic. Here the CMH test is used +when the criterion holds, and Fisher's exact test on the collapsed table +otherwise: + +```{r} +if (is_mf_satisfied) { + prop_cmh(tbl) +} else { + prop_fisher(table(grp, rsp)) +} +``` + +## Empty strata + +Strata that contain no observations carry no information and are dropped before +the statistic is computed. If *every* stratum is empty there is nothing to +compute, so the criterion is undefined and `mantel_fleiss_crit()` returns `NA` +(with an `NA` value attribute when `include_value = TRUE`): + +```{r} +empty_tbl <- table( + factor(character(0), levels = c("Active", "Control")), + factor(logical(0), levels = c("TRUE", "FALSE")), + factor(character(0), levels = "A") +) + +mantel_fleiss_crit(empty_tbl, include_value = TRUE) +``` + +When branching on the result, remember to handle this `NA` case explicitly if +your data can produce fully empty tables. + +## References