diff --git a/DESCRIPTION b/DESCRIPTION index e27fd3d751..898e583598 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,6 +20,7 @@ Authors@R: c( person("Chendi", "Liao", role = "aut"), person("Jennifer", "Li", role = "aut"), person("David", "Munoz Tord", , "david.munoztord@mailbox.org", role = "ctb"), + person("Wojciech", "Wojciak", , "wojciech.wojciak@gmail.com", role = "ctb"), person("F. Hoffmann-La Roche AG", role = c("cph", "fnd")) ) Description: Table, Listings, and Graphs (TLG) library for common outputs diff --git a/NEWS.md b/NEWS.md index c8a5d64dba..5516a00a6e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,13 @@ # tern 0.9.12 * Fixing CRAN release issue. +* Fixed `prop_diff_cmh()` giving wrong Sato and Miettinen-Nurminen results when + some strata have only one group, and failing with a single stratum. (#1535) +* `prop_diff_cmh()` returns `weights`, `n1` and `n2` for all strata now, with `NA` + weights for empty strata. It returns `NA` instead of `0` when no stratum has + both groups. (#1535) +* `prop_cmh()` returns `NA` instead of `1` for the Sato p-value when the response + has no variation. (#1535) # tern 0.9.11 diff --git a/R/prop_diff.R b/R/prop_diff.R index fc16759b55..8fc01cda37 100644 --- a/R/prop_diff.R +++ b/R/prop_diff.R @@ -760,6 +760,7 @@ NULL #' @param correct (`flag`)\cr whether to include the continuity correction. For further #' information, see [stats::prop.test()]. #' +#' @order 1 #' @examples #' # Wald confidence interval #' set.seed(2) @@ -804,6 +805,7 @@ prop_diff_wald <- function(rsp, #' @describeIn h_prop_diff Anderson-Hauck confidence interval \insertCite{HauckAnderson1986}{tern}. #' +#' @order 2 #' @examples #' # Anderson-Hauck confidence interval #' ## "Mid" case: 3/4 respond in group A, 1/2 respond in group B. @@ -840,8 +842,10 @@ prop_diff_ha <- function(rsp, } #' @describeIn h_prop_diff Newcombe confidence interval. It is based on -#' the Wilson score confidence interval for a single binomial proportion \insertCite{Newcombe1998}{tern}. +#' the Wilson score confidence interval for a single binomial proportion +#' \insertCite{Newcombe1998}{tern}. #' +#' @order 3 #' @examples #' # Newcombe confidence interval #' @@ -884,91 +888,19 @@ prop_diff_nc <- function(rsp, ) } -#' @describeIn h_prop_diff Helper function to calculate the CMH weighted -#' difference in proportions. -#' -#' @param tbl (`array`)\cr 3-dimensional array with dimensions corresponding to -#' group, response, and strata. The second dimension (response) should have names -#' "TRUE" and "FALSE". -#' -#' @keywords internal -h_diff_cmh <- function(tbl) { - checkmate::assert_array(tbl) - checkmate::assert_integer(c(ncol(tbl), nrow(tbl)), lower = 2, upper = 2) - checkmate::assert_integer(length(dim(tbl)), lower = 3, upper = 3) - checkmate::assert_true(identical(dimnames(tbl)[[2]], c("TRUE", "FALSE"))) - - n1 <- colSums(tbl[1, 1:2, ]) - n2 <- colSums(tbl[2, 1:2, ]) - x1 <- tbl[1, 1, ] - p1 <- x1 / n1 - x2 <- tbl[2, 1, ] - p2 <- x2 / n2 - # CMH weights - use_stratum <- (n1 > 0) & (n2 > 0) - n1 <- n1[use_stratum] - n2 <- n2[use_stratum] - p1 <- p1[use_stratum] - p2 <- p2[use_stratum] - wt <- (n1 * n2 / (n1 + n2)) - wt_normalized <- wt / sum(wt) - est1 <- sum(wt_normalized * p1) - est2 <- sum(wt_normalized * p2) - diff_est <- est2 - est1 - - list( - x1 = x1, - n1 = n1, - p1 = p1, - x2 = x2, - n2 = n2, - p2 = p2, - wt = wt, - wt_normalized = wt_normalized, - est1 = est1, - est2 = est2, - diff_est = diff_est - ) -} - -#' @describeIn h_prop_diff Helper function to calculate the standard error for the -#' CMH weighted difference in proportions. -#' -#' @param cmh_results (`list`)\cr output of [h_diff_cmh()]. -#' -#' @keywords internal -h_diff_cmh_se <- function(cmh_results, diff_se = c("standard", "sato")) { - checkmate::assert_list(cmh_results, types = "numeric", any.missing = FALSE, names = "unique") - diff_se <- match.arg(diff_se) - l <- cmh_results # For easier readability of the formulas below. - if (diff_se == "standard") { - terms1 <- l$p1 * (1 - l$p1) / l$n1 - terms2 <- l$p2 * (1 - l$p2) / l$n2 - sqrt(sum((terms1 + terms2) * l$wt_normalized^2)) - } else { - # Sato variance estimator. - p_terms_num <- l$n2^2 * l$x1 - - l$n1^2 * l$x2 + - l$n1 * l$n2 * (l$n1 - l$n2) / 2 - p_terms <- p_terms_num / (l$n1 + l$n2)^2 - q_terms_num <- l$x1 * (l$n2 - l$x2) + l$x2 * (l$n1 - l$x1) - q_terms <- q_terms_num / (2 * (l$n1 + l$n2)) - num <- l$diff_est * sum(p_terms) + sum(q_terms) - denom <- sum(l$wt)^2 - sqrt(num / denom) - } -} - -#' @describeIn h_prop_diff Calculates the weighted difference. This is defined as the difference in -#' response rates between the experimental treatment group and the control treatment group, adjusted -#' for stratification factors by applying Cochran-Mantel-Haenszel (CMH) weights. For the CMH chi-squared +#' @describeIn h_prop_diff Calculates the weighted difference. This is defined +#' as the difference in response rates between the experimental treatment +#' group and the control treatment group, adjusted for stratification factors +#' by applying Cochran-Mantel-Haenszel (CMH) weights. For the CMH chi-squared #' test, use [stats::mantelhaen.test()]. #' -#' @param strata (`factor`)\cr variable with one level per stratum and same length as `rsp`. -#' @param diff_se (`string`)\cr method to estimate the standard error for the difference, either -#' `standard`, `sato` \insertCite{Sato1989}{tern} or +#' @param strata (`factor`)\cr variable with one level per stratum and same +#' length as `rsp`. +#' @param diff_se (`string`)\cr method to estimate the standard error for the +#' difference, either `standard`, `sato` \insertCite{Sato1989}{tern} or #' `miettinen_nurminen` \insertCite{MiettinenNurminen1985}{tern}. #' +#' @order 4 #' @examples #' # Cochran-Mantel-Haenszel confidence interval #' @@ -997,142 +929,51 @@ prop_diff_cmh <- function(rsp, strata, conf_level = 0.95, diff_se = c("standard", "sato", "miettinen_nurminen")) { + diff_se <- match.arg(diff_se) + grp <- as_factor_keep_attributes(grp) strata <- as_factor_keep_attributes(strata) - diff_se <- match.arg(diff_se) check_diff_prop_ci( rsp = rsp, grp = grp, conf_level = conf_level, strata = strata ) - if (any(tapply(rsp, strata, length, default = 0) < 5)) { - warning("Less than 5 observations in some strata.") - } - - # first dimension: CONTROL, TX + # 1st dimension: CONTROL, TX # 2nd dimension: TRUE, FALSE # 3rd dimension: levels of strata - # Note: rsp needs to be a factor to handle edge case of - # no FALSE (or TRUE) rsp records. - t_tbl <- table( - grp, - factor(rsp, levels = c("TRUE", "FALSE")), - strata - ) - cmh <- h_diff_cmh(t_tbl) + # Note: rsp needs to be a factor to handle edge case of no FALSE (or TRUE). + tbl <- table(grp, factor(rsp, levels = c("TRUE", "FALSE")), strata) - estimate <- c(cmh$est1, cmh$est2) - names(estimate) <- levels(grp) - se1 <- sqrt(sum(cmh$wt_normalized^2 * cmh$p1 * (1 - cmh$p1) / cmh$n1)) - se2 <- sqrt(sum(cmh$wt_normalized^2 * cmh$p2 * (1 - cmh$p2) / cmh$n2)) - z <- stats::qnorm((1 + conf_level) / 2) - err1 <- z * se1 - err2 <- z * se2 - ci1 <- c((cmh$est1 - err1), (cmh$est1 + err1)) - ci2 <- c((cmh$est2 - err2), (cmh$est2 + err2)) - estimate_ci <- list(ci1, ci2) - names(estimate_ci) <- levels(grp) + if (any(marginSums(tbl, margin = 3L) < 5L)) { + warning("Less than 5 observations in some strata.") + } + + prop <- h_prop_cmh(tbl, conf_level = conf_level) + prop_diff_est <- unname(prop$est2 - prop$est1) if (diff_se %in% c("standard", "sato")) { - se_diff <- h_diff_cmh_se(cmh, diff_se = diff_se) - diff_ci <- c(cmh$diff_est - z * se_diff, cmh$diff_est + z * se_diff) - } else { - # Miettinen and Nurminen method is used. - z_stat_fun <- function(delta) { - var_est <- h_miettinen_nurminen_var_est( - n1 = cmh$n1, n2 = cmh$n2, - x1 = cmh$x1, x2 = cmh$x2, - diff_par = delta - )$var_est - num <- sum(cmh$wt * (cmh$p2 - cmh$p1 - delta)) - denom <- sqrt(sum(cmh$wt^2 * var_est)) - num / denom + prop_diff_var <- if (diff_se == "standard") { + unname(prop$var1 + prop$var2) + } else { # "sato" + h_cmh_sato_var(prop) } - # Find upper and lower confidence limits by root finding such that - # z_stat_fun(limit) = +/- z quantile: - root_lower <- function(delta) z_stat_fun(delta) - z - root_upper <- function(delta) z_stat_fun(delta) + z - diff_ci <- c( - stats::uniroot(root_lower, interval = c(-0.99, cmh$diff_est))$root, - stats::uniroot(root_upper, interval = c(cmh$diff_est, 0.99))$root - ) - # Calculate the standard error separately. - var_est <- h_miettinen_nurminen_var_est( - n1 = cmh$n1, n2 = cmh$n2, - x1 = cmh$x1, x2 = cmh$x2, - diff_par = cmh$diff_est - )$var_est - se_diff <- sqrt(sum(cmh$wt_normalized^2 * var_est)) + prop_diff_se <- sqrt(prop_diff_var) + z <- stats::qnorm((1 + conf_level) / 2) + prop_diff_ci <- prop_diff_est + c(-1, 1) * z * prop_diff_se + } else { # "miettinen_nurminen" + mn <- h_miettinen_nurminen_stratified_ci(prop, conf_level = conf_level) + prop_diff_se <- mn$se + prop_diff_ci <- mn$ci } list( - prop = estimate, - prop_ci = estimate_ci, - diff = cmh$diff_est, - diff_ci = diff_ci, - se_diff = se_diff, - weights = cmh$wt_normalized, - n1 = cmh$n1, - n2 = cmh$n2 - ) -} - -#' Variance Estimates in Strata following Miettinen and Nurminen -#' -#' The variable names in this function follow the notation in the original -#' paper by \insertCite{MiettinenNurminen1985;textual}{tern}, cf. Appendix 1. -#' -#' @param n1 (`numeric`)\cr sample sizes in group 1. -#' @param n2 (`numeric`)\cr sample sizes in group 2. -#' @param x1 (`numeric`)\cr number of responders in group 1. -#' @param x2 (`numeric`)\cr number of responders in group 2. -#' @param diff_par (`numeric`)\cr assumed difference in true proportions -#' (group 2 minus group 1). -#' @return A named `list` with elements: -#' -#' - `p1_hat`: estimated proportion in group 1 -#' - `p2_hat`: estimated proportion in group 2 -#' - `var_est`: variance estimate of the difference in proportions -#' -#' @keywords internal -#' @references -#' \insertAllCited{} -h_miettinen_nurminen_var_est <- function(n1, n2, x1, x2, diff_par) { - # nolint start - # Translate to the notation in the paper. - S0 <- n1 - S1 <- n2 - c0 <- x1 - c1 <- x2 - RD <- diff_par - - # Further definitions. - S <- S0 + S1 - c <- c0 + c1 - - # Coefficients of the third-degree polynomial. - L3 <- S - L2 <- (S1 + 2 * S0) * RD - S - c - L1 <- (S0 * RD - S - 2 * c0) * RD + c - L0 <- c0 * RD * (1 - RD) - # nolint end - - # Solution for group 1 proportion. - q <- L2^3 / (3 * L3)^3 - L1 * L2 / (6 * L3^2) + L0 / (2 * L3) - p <- sign(q) * sqrt(L2^2 / (3 * L3)^2 - L1 / (3 * L3)) - a <- (1 / 3) * (base::pi + acos(q / p^3)) - p1_hat <- 2 * p * cos(a) - L2 / (3 * L3) - - # Estimated group 2 proportion. - p2_hat <- p1_hat + RD - - # Variance estimate. - var_est <- (p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2) * - S / (S - 1) - - list( - p1_hat = p1_hat, - p2_hat = p2_hat, - var_est = var_est + prop = prop$est_both_groups, + prop_ci = prop$ci_both_groups, + diff = prop_diff_est, + diff_ci = prop_diff_ci, + se_diff = prop_diff_se, + weights = prop$w_normalized, + n1 = prop$n1, + n2 = prop$n2 ) } @@ -1149,6 +990,7 @@ h_miettinen_nurminen_var_est <- function(n1, n2, x1, x2, diff_par) { #' the Cochran-Mantel-Haenszel method, while `"wilson_h"` uses the heuristic #' weights proposed by [prop_strat_wilson()]. #' +#' @order 5 #' @examples #' # Stratified Newcombe confidence interval #' @@ -1249,6 +1091,535 @@ prop_diff_strat_nc <- function(rsp, ) } + +#' @describeIn h_prop_diff Unconditional exact confidence interval for the difference in +#' proportions by inverting one-sided tail tests over a nuisance parameter. This is +#' the "tail method" described by Santner and Snell \insertCite{SantnerSnell1980}{tern}. +#' +#' @order 6 +#' @examples +#' # Unconditional exact confidence interval +#' n11 <- 40 +#' n21 <- 5 +#' n1 <- 78 +#' n2 <- 17 +#' rsp <- c(rep(TRUE, n21), rep(FALSE, n2 - n21), rep(TRUE, n11), rep(FALSE, n1 - n11)) +#' grp <- factor(c(rep("B", n2), rep("A", n1)), levels = c("B", "A")) +#' +#' prop_diff_uncond_exact(rsp = rsp, grp = grp, conf_level = 0.95) +#' +#' @export +prop_diff_uncond_exact <- function(rsp, + grp, + conf_level = 0.95) { + grp <- as_factor_keep_attributes(grp) + check_diff_prop_ci(rsp = rsp, grp = grp, conf_level = conf_level) + + alpha <- 1 - conf_level + cutoff <- alpha / 2 + + tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) + + # Step 0: Calculate the observed difference in proportions + # and the observed test statistic value. + n2 <- sum(tbl[1, ]) + n1 <- sum(tbl[2, ]) + + if (n1 == 0 || n2 == 0) { + return(list( + diff = NaN, + diff_ci = c(NaN, NaN) + )) + } + + n21_obs <- tbl[1, 1] + n11_obs <- tbl[2, 1] + diff_est <- n11_obs / n1 - n21_obs / n2 + + # Step 1: Enumerate all tables in A with fixed row margins + # n1 and n2. + if (n1 * n2 > 1e5) { + warning("uncond_exact_diff: Large sample sizes may lead to long computation time.") + } + tables <- expand.grid( + n11 = 0:n1, + n21 = 0:n2 + ) + + # Step 2: Compute T(a) = n11 / n1 - n21 / n2 for each table a in A. + t_values <- tables$n11 / n1 - tables$n21 / n2 + t0 <- diff_est + + # Step 3: For each hypothesized difference d*, compute the worst-case + # tail probabilities P_U(d*) and P_L(d*) by maximizing over the nuisance + # parameter p2. + p_upper <- function(d_star) { + # Step 4a: Compute worst-case one-sided tail probability: + # P_U(d*) = sup_p2 sum_{T(a) >= t0} f(...) + h_worst_case_tail_probability( + d_star = d_star, + n1 = n1, + n2 = n2, + t_values = t_values, + t0 = t0, + tables = tables, + tail = "upper" + ) + } + p_lower <- function(d_star) { + # Step 4b: Compute worst-case one-sided tail probability: + # P_L(d*) = sup_p2 sum_{T(a) <= t0} f(...) + h_worst_case_tail_probability( + d_star = d_star, + n1 = n1, + n2 = n2, + t_values = t_values, + t0 = t0, + tables = tables, + tail = "lower" + ) + } + + # Step 5: Invert one-sided tests to obtain the two-sided + # 100 * (1 - alpha)% CI for d = p1 - p2. + # For monotone one-sided p-value functions, use uniroot to solve + # P_U(d) = alpha/2 and P_L(d) = alpha/2 directly. + diff_ci <- c( + h_find_ci_bound_uniroot(p_upper, cutoff = cutoff, direction = "increasing"), + h_find_ci_bound_uniroot(p_lower, cutoff = cutoff, direction = "decreasing") + ) + + list( + diff = diff_est, + diff_ci = diff_ci + ) +} + +#' Helper function to calculate the CMH-weighted proportions and their +#' confidence intervals. +#' +#' @description `r lifecycle::badge("stable")` +#' +#' @param tbl (`array`)\cr +#' A three-dimensional contingency table containing counts for each +#' combination of group, response, and stratum, in that order. +#' The first two dimensions must each have exactly two levels, and the second +#' dimension (response) must have names `"TRUE"` and `"FALSE"`. +#' At least one stratum must be present. Strata with all cell counts equal to +#' zero are allowed. +#' All cell values must be finite, non-missing integer counts. +#' @param conf_level (`number(1)`)\cr +#' Confidence level for the confidence intervals. +#' +#' @return A named list containing the CMH-weighted proportion estimates, +#' confidence intervals, and intermediate quantities. +#' The stratum-specific quantities `x1`, `n1`, `p1`, `x2`, `n2`, `p2`, `w`, +#' and `w_normalized` are vectors with a length equal to the number of strata +#' in `tbl` and retain the same stratum order. +#' Some of these quantities may be `NA` for strata where they are not defined. +#' In particular, `p1` or `p2` is `NA` when the corresponding group has no +#' observations in that stratum. +#' +#' `est1` and `est2` are the overall CMH-weighted proportion estimates for the +#' two groups, respectively. `est_both_groups` contains these two estimates in +#' group order. +#' `ci_both_groups` contains the corresponding confidence intervals in the same +#' group order. +#' +#' If no stratum contains observations in both groups, the CMH weights +#' cannot be normalized and the overall estimates and confidence intervals +#' are `NA` +#' +#' @seealso [prop_diff_cmh()] +#' @keywords internal +#' +h_prop_cmh <- function(tbl, conf_level = 0.95) { + checkmate::assert_array(tbl, mode = "integerish", any.missing = FALSE, d = 3L) + checkmate::assert_true(nrow(tbl) == 2L) + checkmate::assert_true(ncol(tbl) == 2L) + checkmate::assert_true(dim(tbl)[3L] > 0L) + checkmate::assert_true(identical(dimnames(tbl)[[2]], c("TRUE", "FALSE"))) + checkmate::assert_true(all(tbl >= 0)) + checkmate::assert_true(all(is.finite(tbl))) + assert_proportion_value(conf_level) + + strata_names <- dimnames(tbl)[[3L]] # Can be NULL. + + x1 <- setNames(tbl[1L, "TRUE", ], strata_names) + x2 <- setNames(tbl[2L, "TRUE", ], strata_names) + n1 <- apply(tbl[1L, , , drop = FALSE], MARGIN = 3L, sum) + n2 <- apply(tbl[2L, , , drop = FALSE], MARGIN = 3L, sum) + p1 <- ifelse(n1 > 0, x1 / n1, NA_real_) + p2 <- ifelse(n2 > 0, x2 / n2, NA_real_) + + # CMH weights. + w <- ifelse(n1 + n2 > 0, (n1 * n2) / (n1 + n2), NA_real_) + w_sum <- sum(w, na.rm = TRUE) + + if (w_sum > 0) { + # In addition to ensuring a non-zero denominator, w_sum > 0 ensures that + # for at least one stratum h, w[h] is non-NA and > 0, and therefore, + # n1[h] > 0 and n2[h] > 0. + # Consequently, p1[h], p2[h], and w_normalized[h] are all non-NA. + # Thus, all four sums below contain at least one non-NA element and cannot + # yield an unjustified 0. + # This is important to note because sum(numeric(0)) returns 0. + w_normalized <- w / w_sum + est1 <- sum(w_normalized * p1, na.rm = TRUE) + est2 <- sum(w_normalized * p2, na.rm = TRUE) + + var1 <- sum(w_normalized^2 * p1 * (1 - p1) / n1, na.rm = TRUE) + var2 <- sum(w_normalized^2 * p2 * (1 - p2) / n2, na.rm = TRUE) + z <- stats::qnorm((1 + conf_level) / 2) + ci1 <- est1 + c(-1, 1) * z * sqrt(var1) + ci2 <- est2 + c(-1, 1) * z * sqrt(var2) + } else { + w_normalized <- setNames(rep(NA_real_, dim(tbl)[3L]), strata_names) + est1 <- est2 <- var1 <- var2 <- NA_real_ + ci1 <- ci2 <- c(NA_real_, NA_real_) + } + + group_names <- dimnames(tbl)[[1L]] # Can be NULL. + + list( + x1 = x1, n1 = n1, p1 = p1, # Quantities for group 1. + x2 = x2, n2 = n2, p2 = p2, # Quantities for group 2. + w = w, + w_normalized = w_normalized, + est1 = setNames(est1, group_names[1L]), + est2 = setNames(est2, group_names[2L]), + est_both_groups = setNames(c(est1, est2), group_names), + var1 = setNames(var1, group_names[1L]), + var2 = setNames(var2, group_names[2L]), + ci_both_groups = setNames(list(ci1, ci2), group_names) + ) +} + +#' Sato Variance Estimate for the CMH-weighted Difference in Proportions +#' +#' @description `r lifecycle::badge("stable")` +#' +#' Calculates the Sato variance estimate for the difference between two +#' Cochran-Mantel-Haenszel (CMH)-weighted proportions. The estimate is used to +#' obtain the standard error and confidence interval for the stratified +#' difference in response proportions. +#' +#' The calculation follows the variance estimator proposed by +#' \insertCite{Sato1989;textual}{tern}. The required stratum-specific counts, +#' sample sizes, CMH weights, and overall CMH-weighted proportion estimates are +#' supplied in the `prop` object returned by [h_prop_cmh()]. +#' +#' @details +#' `h_cmh_sato_var()` takes a `prop` list as returned by [h_prop_cmh()]. +#' The `prop` object must contain vectors `est1`, `est2`, `x1`, `x2`, +#' `n1`, `n2`, and `w`, which provide the overall CMH-weighted estimates +#' and the stratum-specific quantities required for the Sato variance +#' calculation. +#' +#' @param prop (`list`)\cr +#' A named list returned by [h_prop_cmh()]. It must contain the following +#' atomic vectors: +#' \describe{ +#' \item{`est1`}{CMH-weighted estimated proportion for group 1. +#' May be `NA_real_` when a CMH-weighted estimate cannot be calculated.} +#' \item{`est2`}{CMH-weighted estimated proportion for group 2. +#' May be `NA_real_` when a CMH-weighted estimate cannot be calculated.} +#' \item{`x1`}{Number of responders in group 1 for each stratum.} +#' \item{`x2`}{Number of responders in group 2 for each stratum.} +#' \item{`n1`}{Number of observations in group 1 for each stratum.} +#' \item{`n2`}{Number of observations in group 2 for each stratum.} +#' \item{`w`}{Unnormalized CMH weights for each stratum.} +#' } +#' +#' The vectors `x1`, `x2`, `n1`, `n2`, and `w` must be of the same length. +#' +#' The unnormalized CMH weights for stratum \eqn{i} given by +#' \deqn{ +#' \frac{n_{1i} n_{2i}}{n_{1i} + n_{2i}}, +#' } +#' for \eqn{n_{1i} + n_{2i} > 0}, where \eqn{n_{1i}} is the total number of +#' observations in group \eqn{1} in stratum \eqn{i}, and \eqn{n_{2i}} is the +#' total number of observations in group \eqn{2} in stratum \eqn{i}. +#' +#' Missing weights in `w` are allowed and are ignored when calculating their +#' sum. +#' +#' @return A `numeric(1)` containing the Sato estimate of the variance of +#' the CMH-weighted difference in proportions. Returns `NA_real_` when +#' the variance cannot be estimated because there are no usable strata +#' or the sum of the supplied CMH weights is zero. +#' +#' @seealso [prop_diff_cmh()], [h_prop_cmh()] +#' +#' @references +#' \insertAllCited{} +#' +#' @keywords internal +#' +h_cmh_sato_var <- function(prop) { + checkmate::assert_list(prop, min.len = 7L, names = "named") + checkmate::assert_subset(c("est1", "est2", "x1", "x2", "n1", "n2", "w"), names(prop)) + checkmate::assert_number(prop$est1, lower = -1, upper = 1, na.ok = TRUE, finite = TRUE) + checkmate::assert_number(prop$est2, lower = -1, upper = 1, na.ok = TRUE, finite = TRUE) + checkmate::assert_integerish(prop$x1, min.len = 1L, lower = 0, any.missing = FALSE) + checkmate::assert_integerish(prop$x2, len = length(prop$x1), lower = 0, any.missing = FALSE) + checkmate::assert_integerish(prop$n1, len = length(prop$x1), lower = 0, any.missing = FALSE) + checkmate::assert_integerish(prop$n2, len = length(prop$x1), lower = 0, any.missing = FALSE) + checkmate::assert_true(all(prop$x1 <= prop$n1)) + checkmate::assert_true(all(prop$x2 <= prop$n2)) + checkmate::assert_numeric(prop$w, len = length(prop$x1), lower = 0, finite = TRUE) + + # For easier readability of the formulas below. + est1 <- prop$est1 + est2 <- prop$est2 + x1 <- prop$x1 + x2 <- prop$x2 + n1 <- prop$n1 + n2 <- prop$n2 + w_unnormalized <- prop$w + + n <- n1 + n2 + + p_numerator <- n2^2 * x1 - n1^2 * x2 + n1 * n2 * (n1 - n2) / 2 + p <- ifelse(n > 0, p_numerator / n^2, NA_real_) + + q_numerator <- x1 * (n2 - x2) + x2 * (n1 - x1) + q <- ifelse(n > 0, q_numerator / (2 * n), NA_real_) + + w_sum <- sum(w_unnormalized, na.rm = TRUE) + if (any(!is.na(p)) && w_sum > 0) { # Note: any(!is.na(p)) == TRUE <=> any(!is.na(q)) == TRUE. + num <- (est2 - est1) * sum(p, na.rm = TRUE) + sum(q, na.rm = TRUE) + unname(num) / w_sum^2 + } else { + NA_real_ + } +} + +#' Variance Estimate Following Miettinen and Nurminen +#' +#' @description `r lifecycle::badge("stable")` +#' +#' Calculates the \insertCite{MiettinenNurminen1985;textual}{tern} variance +#' estimate for the difference between two proportions. The estimate is based on +#' the constrained maximum likelihood estimates of the two proportions under the +#' specified risk difference and is used to obtain the standard error for the +#' Miettinen-Nurminen confidence interval. +#' +#' @details +#' The risk difference is defined as `est2` - `est1`. For each stratum, the +#' function calculates the constrained maximum likelihood estimate for the +#' proportion in group 1 and obtains the corresponding estimate for group 2 +#' by adding the risk difference. The variance is then calculated from these +#' estimates using the Miettinen-Nurminen variance formula. +#' +#' The variance is returned as `NA_real_` for strata where the variance cannot +#' be calculated. +#' +#' The variable names in this function follow the notation in the original +#' paper by \insertCite{MiettinenNurminen1985;textual}{tern}, cf. Appendix 1. +#' +#' @param est1 (`numeric(1)`) \cr +#' Estimated proportion for group 1. Used together with `est2` to define the +#' risk difference. +#' May be `NA_real_` when an estimate cannot be calculated. +#' @param est2 (`numeric(1)`) \cr +#' Estimated proportion for group 2. Used together with `est1` to define the +#' risk difference. +#' May be `NA_real_` when an estimate cannot be calculated. +#' @param x1 (`numeric`) \cr +#' Number of responders in group 1 for each stratum. +#' Must have length at least 1. +#' @param x2 (`numeric`) \cr +#' Number of responders in group 2 for each stratum. +#' Must have the same length as `x1`. +#' @param n1 (`numeric`) \cr +#' Number of observations in group 1 for each stratum. +#' Must have the same length as `x1`. +#' @param n2 (`numeric`) \cr +#' Number of observations in group 2 for each stratum. +#' Must have the same length as `x1`. +#' +#' @return A named `list` with elements: +#' +#' - `p1_est`: constrained maximum likelihood estimate of the proportion in +#' group 1 for each stratum. +#' - `p2_est`: constrained maximum likelihood estimate of the proportion in +#' group 2 for each stratum. +#' - `var_est`: Miettinen-Nurminen variance estimate for each stratum. +#' +#' @seealso [prop_diff_cmh()], [h_prop_cmh()], [h_miettinen_nurminen_stratified_ci()] +#' @references +#' \insertAllCited{} +#' +#' @keywords internal +#' +h_miettinen_nurminen_var <- function(est1, est2, x1, x2, n1, n2) { + checkmate::assert_number(est1, lower = -1, upper = 1, na.ok = TRUE, finite = TRUE) + checkmate::assert_number(est2, lower = -1, upper = 1, na.ok = TRUE, finite = TRUE) + checkmate::assert_integerish(x1, min.len = 1L, lower = 0, any.missing = FALSE) + checkmate::assert_integerish(x2, len = length(x1), lower = 0, any.missing = FALSE) + checkmate::assert_integerish(n1, len = length(x1), lower = 0, any.missing = FALSE) + checkmate::assert_integerish(n2, len = length(x1), lower = 0, any.missing = FALSE) + checkmate::assert_true(all(x1 <= n1)) + checkmate::assert_true(all(x2 <= n2)) + + # nolint start + # Translate to the notation in the paper. + S0 <- n1 + S1 <- n2 + c0 <- x1 + c1 <- x2 + RD <- est2 - est1 + + # Further definitions. + S <- S0 + S1 + c <- c0 + c1 + + # Coefficients of the third-degree polynomial. + L3 <- S + L2 <- (S1 + 2 * S0) * RD - S - c + L1 <- (S0 * RD - S - 2 * c0) * RD + c + L0 <- c0 * RD * (1 - RD) + # nolint end + + # Solution for group 1 proportion. + q <- L2^3 / (3 * L3)^3 - L1 * L2 / (6 * L3^2) + L0 / (2 * L3) + p <- sign(q) * sqrt(L2^2 / (3 * L3)^2 - L1 / (3 * L3)) + a <- (1 / 3) * (base::pi + acos(q / p^3)) + p1_mle <- 2 * p * cos(a) - L2 / (3 * L3) + + # Estimated group 2 proportion. + p2_mle <- p1_mle + RD + + # Variance estimate. + var_est <- ifelse( + n1 > 0 & n2 > 0 & S > 1, + (p1_mle * (1 - p1_mle) / n1 + p2_mle * (1 - p2_mle) / n2) * S / (S - 1), + NA_real_ + ) + + list( + p1_est = p1_mle, + p2_est = p2_mle, + var_est = var_est + ) +} + +#' Stratified Miettinen-Nurminen Confidence Interval +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' Calculates the stratified Miettinen-Nurminen confidence interval and +#' standard error for the difference in proportions. The method uses +#' constrained maximum likelihood estimates within each stratum and combines +#' the stratum-specific variance estimates using the normalized CMH weights. +#' +#' @details +#' The difference in proportions is defined as the proportion in group 2 +#' minus the proportion in group 1. +#' +#' @param prop (`list`)\cr +#' A named list returned by [h_prop_cmh()]. It must contain the following +#' atomic vectors: +#' \describe{ +#' \item{`est1`}{CMH-weighted estimated proportion for group 1. +#' May be `NA_real_` when a CMH-weighted estimate cannot be calculated.} +#' \item{`est2`}{CMH-weighted estimated proportion for group 2. +#' May be `NA_real_` when a CMH-weighted estimate cannot be calculated.} +#' \item{`x1`}{Number of responders in group 1 for each stratum.} +#' \item{`x2`}{Number of responders in group 2 for each stratum.} +#' \item{`n1`}{Number of observations in group 1 for each stratum.} +#' \item{`n2`}{Number of observations in group 2 for each stratum.} +#' \item{`p1`}{Observed response proportion in group 1 for each stratum.} +#' \item{`p2`}{Observed response proportion in group 2 for each stratum.} +#' \item{`w`}{Unnormalized CMH weight for each stratum.} +#' \item{`w_normalized`}{Normalized CMH weight for each stratum.} +#' } +#' +#' @param conf_level (`number(1)`) \cr +#' Confidence level for the confidence interval. +#' +#' @return A named list containing: +#' \describe{ +#' \item{`ci`}{(`numeric(2)`) Lower and upper confidence limits for the +#' stratified difference in proportions.} +#' \item{`se`}{(`numeric(1)`) Standard error of the stratified difference +#' in proportions.} +#' } +#' +#' @seealso [prop_diff_cmh()], [h_prop_cmh()], [h_miettinen_nurminen_var()] +#' +#' @keywords internal +h_miettinen_nurminen_stratified_ci <- function(prop, conf_level = 0.95) { + checkmate::assert_list(prop, min.len = 10L, names = "named") + checkmate::assert_subset( + c("est1", "est2", "x1", "x2", "n1", "n2", "p1", "p2", "w", "w_normalized"), + names(prop) + ) + checkmate::assert_number(prop$est1, lower = -1, upper = 1, na.ok = TRUE, finite = TRUE) + checkmate::assert_number(prop$est2, lower = -1, upper = 1, na.ok = TRUE, finite = TRUE) + checkmate::assert_numeric(prop$p1, min.len = 1L, lower = 0, upper = 1, finite = TRUE) + checkmate::assert_numeric(prop$p2, len = length(prop$p1), lower = 0, upper = 1, finite = TRUE) + checkmate::assert_numeric(prop$w, len = length(prop$p1), lower = 0, finite = TRUE) + checkmate::assert_numeric(prop$w_normalized, len = length(prop$p1), lower = 0, finite = TRUE) + assert_proportion_value(conf_level) + + p <- prop # For easier readability of the formulas below. + + if (is.na(p$est1) || is.na(p$est2)) { + return(list(ci = c(NA_real_, NA_real_), se = NA_real_)) + } + + # Calculate the standard error. + var_est <- h_miettinen_nurminen_var( + est1 = p$est1, est2 = p$est2, + x1 = p$x1, x2 = p$x2, + n1 = p$n1, n2 = p$n2 + )$var_est + + w_var <- p$w_normalized^2 * var_est + se <- if (any(!is.na(w_var))) { + sqrt(sum(w_var, na.rm = TRUE)) + } else { + NA_real_ + } + + # Calculate the confidence interval. + + # Stratified Miettinen-Nurminen score function. + score_fun <- function(delta) { + var_est <- h_miettinen_nurminen_var( + est1 = 0, est2 = delta, + x1 = p$x1, x2 = p$x2, + n1 = p$n1, n2 = p$n2 + )$var_est + + # Ensure that both the numerator and denominator are computed + # using the same set of strata. + non_na <- !is.na(p$w) & !is.na(p$p1) & !is.na(p$p2) & !is.na(var_est) + + denom <- sqrt(sum(p$w[non_na]^2 * var_est[non_na])) + if (any(non_na) && denom > 0) { + sum(p$w[non_na] * (p$p2[non_na] - p$p1[non_na] - delta)) / denom + } else { + NA_real_ + } + } + + # Confidence interval consists of all values of delta for which + # score_fun(delta) falls in the two-sided acceptance region, + # {delta: -z <= score_fun(delta) <= z}, where z = z_{1 - alpha/2}. + z <- stats::qnorm((1 + conf_level) / 2) + root_lower <- function(delta) score_fun(delta) - z + root_upper <- function(delta) score_fun(delta) + z + ci <- c( + uniroot_catch_na(root_lower, interval = c(-0.99, p$est2 - p$est1)), + uniroot_catch_na(root_upper, interval = c(p$est2 - p$est1, 0.99)) + ) + + list(ci = ci, se = se) +} + #' Worst case tail probability for unconditional exact CI calculation #' #' This function is an internal helper for [prop_diff_uncond_exact()]. @@ -1379,105 +1750,3 @@ h_find_ci_bound_uniroot <- function(p_value_function, maxiter = maxiter )$root } - -#' @describeIn h_prop_diff Unconditional exact confidence interval for the difference in -#' proportions by inverting one-sided tail tests over a nuisance parameter. This is -#' the "tail method" described by Santner and Snell \insertCite{SantnerSnell1980}{tern}. -#' -#' @examples -#' # Unconditional exact confidence interval -#' n11 <- 40 -#' n21 <- 5 -#' n1 <- 78 -#' n2 <- 17 -#' rsp <- c(rep(TRUE, n21), rep(FALSE, n2 - n21), rep(TRUE, n11), rep(FALSE, n1 - n11)) -#' grp <- factor(c(rep("B", n2), rep("A", n1)), levels = c("B", "A")) -#' -#' prop_diff_uncond_exact(rsp = rsp, grp = grp, conf_level = 0.95) -#' -#' @export -prop_diff_uncond_exact <- function(rsp, - grp, - conf_level = 0.95) { - grp <- as_factor_keep_attributes(grp) - check_diff_prop_ci(rsp = rsp, grp = grp, conf_level = conf_level) - - alpha <- 1 - conf_level - cutoff <- alpha / 2 - - tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) - - # Step 0: Calculate the observed difference in proportions - # and the observed test statistic value. - n2 <- sum(tbl[1, ]) - n1 <- sum(tbl[2, ]) - - if (n1 == 0 || n2 == 0) { - return(list( - diff = NaN, - diff_ci = c(NaN, NaN) - )) - } - - n21_obs <- tbl[1, 1] - n11_obs <- tbl[2, 1] - diff_est <- n11_obs / n1 - n21_obs / n2 - - # Step 1: Enumerate all tables in A with fixed row margins - # n1 and n2. - if (n1 * n2 > 1e5) { - warning("uncond_exact_diff: Large sample sizes may lead to long computation time.") - } - tables <- expand.grid( - n11 = 0:n1, - n21 = 0:n2 - ) - - # Step 2: Compute T(a) = n11 / n1 - n21 / n2 for each table a in A. - t_values <- tables$n11 / n1 - tables$n21 / n2 - t0 <- diff_est - - # Step 3: For each hypothesized difference d*, compute the worst-case - # tail probabilities P_U(d*) and P_L(d*) by maximizing over the nuisance - # parameter p2. - p_upper <- function(d_star) { - # Step 4a: Compute worst-case one-sided tail probability: - # P_U(d*) = sup_p2 sum_{T(a) >= t0} f(...) - h_worst_case_tail_probability( - d_star = d_star, - n1 = n1, - n2 = n2, - t_values = t_values, - t0 = t0, - tables = tables, - tail = "upper" - ) - } - p_lower <- function(d_star) { - # Step 4b: Compute worst-case one-sided tail probability: - # P_L(d*) = sup_p2 sum_{T(a) <= t0} f(...) - h_worst_case_tail_probability( - d_star = d_star, - n1 = n1, - n2 = n2, - t_values = t_values, - t0 = t0, - tables = tables, - tail = "lower" - ) - } - - # Step 5: Invert one-sided tests to obtain the two-sided - # 100 * (1 - alpha)% CI for d = p1 - p2. - # For monotone one-sided p-value functions, use uniroot to solve - # P_U(d) = alpha/2 and P_L(d) = alpha/2 directly. - diff_ci <- c( - h_find_ci_bound_uniroot(p_upper, cutoff = cutoff, direction = "increasing"), - h_find_ci_bound_uniroot(p_lower, cutoff = cutoff, direction = "decreasing") - ) - - list( - diff = diff_est, - diff_ci = diff_ci - ) -} diff --git a/R/prop_diff_test.R b/R/prop_diff_test.R index 25ef568734..c39c25e9b7 100644 --- a/R/prop_diff_test.R +++ b/R/prop_diff_test.R @@ -412,9 +412,9 @@ prop_cmh <- function(ary, sqrt(unname(mh_res$statistic)) * stat_sign } else { # Use the Sato variance estimator. - cmh <- h_diff_cmh(ary) - cmh_se <- h_diff_cmh_se(cmh, diff_se = "sato") - cmh$diff_est / cmh_se + prop <- h_prop_cmh(ary) + prop_diff_var <- h_cmh_sato_var(prop) + unname(prop$est2 - prop$est1) / sqrt(prop_diff_var) } if (transform == "wilson_hilferty") { diff --git a/R/utils.R b/R/utils.R index 70e5f9ad45..f6a89ed01d 100644 --- a/R/utils.R +++ b/R/utils.R @@ -539,3 +539,30 @@ get_complete_cases <- function(df, quiet = FALSE, additional_message = ".") { df } } + +#' Find a root while returning NA for missing function values +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' A wrapper around [stats::uniroot()] that returns `NA_real_` when `f` is `NA` +#' at either end of `interval`. +#' +#' @param f (`function`)\cr function for which the root is sought. +#' @param interval (`numeric(2)`)\cr end points of the interval to be searched. +#' @param ... further arguments passed to [stats::uniroot()]. +#' +#' @return +#' A numeric scalar containing the found root, or `NA_real_` if `f` returns +#' `NA` at either endpoint of `interval`. +#' +#' @seealso [stats::uniroot()] +#' @keywords internal +uniroot_catch_na <- function(f, interval, ...) { + # Checked here, as the uniroot() error message is translated. + f_lower <- f(min(interval)) + f_upper <- f(max(interval)) + if (is.na(f_lower) || is.na(f_upper)) { + return(NA_real_) + } + stats::uniroot(f, interval = interval, f.lower = f_lower, f.upper = f_upper, ...)$root +} diff --git a/man/h_cmh_sato_var.Rd b/man/h_cmh_sato_var.Rd new file mode 100644 index 0000000000..25fc58fbc3 --- /dev/null +++ b/man/h_cmh_sato_var.Rd @@ -0,0 +1,70 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prop_diff.R +\name{h_cmh_sato_var} +\alias{h_cmh_sato_var} +\title{Sato Variance Estimate for the CMH-weighted Difference in Proportions} +\usage{ +h_cmh_sato_var(prop) +} +\arguments{ +\item{prop}{(\code{list})\cr +A named list returned by \code{\link[=h_prop_cmh]{h_prop_cmh()}}. It must contain the following +atomic vectors: +\describe{ +\item{\code{est1}}{CMH-weighted estimated proportion for group 1. +May be \code{NA_real_} when a CMH-weighted estimate cannot be calculated.} +\item{\code{est2}}{CMH-weighted estimated proportion for group 2. +May be \code{NA_real_} when a CMH-weighted estimate cannot be calculated.} +\item{\code{x1}}{Number of responders in group 1 for each stratum.} +\item{\code{x2}}{Number of responders in group 2 for each stratum.} +\item{\code{n1}}{Number of observations in group 1 for each stratum.} +\item{\code{n2}}{Number of observations in group 2 for each stratum.} +\item{\code{w}}{Unnormalized CMH weights for each stratum.} +} + +The vectors \code{x1}, \code{x2}, \code{n1}, \code{n2}, and \code{w} must be of the same length. + +The unnormalized CMH weights for stratum \eqn{i} given by +\deqn{ + \frac{n_{1i} n_{2i}}{n_{1i} + n_{2i}}, + } +for \eqn{n_{1i} + n_{2i} > 0}, where \eqn{n_{1i}} is the total number of +observations in group \eqn{1} in stratum \eqn{i}, and \eqn{n_{2i}} is the +total number of observations in group \eqn{2} in stratum \eqn{i}. + +Missing weights in \code{w} are allowed and are ignored when calculating their +sum.} +} +\value{ +A \code{numeric(1)} containing the Sato estimate of the variance of +the CMH-weighted difference in proportions. Returns \code{NA_real_} when +the variance cannot be estimated because there are no usable strata +or the sum of the supplied CMH weights is zero. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} + +Calculates the Sato variance estimate for the difference between two +Cochran-Mantel-Haenszel (CMH)-weighted proportions. The estimate is used to +obtain the standard error and confidence interval for the stratified +difference in response proportions. + +The calculation follows the variance estimator proposed by +\insertCite{Sato1989;textual}{tern}. The required stratum-specific counts, +sample sizes, CMH weights, and overall CMH-weighted proportion estimates are +supplied in the \code{prop} object returned by \code{\link[=h_prop_cmh]{h_prop_cmh()}}. +} +\details{ +\code{h_cmh_sato_var()} takes a \code{prop} list as returned by \code{\link[=h_prop_cmh]{h_prop_cmh()}}. +The \code{prop} object must contain vectors \code{est1}, \code{est2}, \code{x1}, \code{x2}, +\code{n1}, \code{n2}, and \code{w}, which provide the overall CMH-weighted estimates +and the stratum-specific quantities required for the Sato variance +calculation. +} +\references{ +\insertAllCited{} +} +\seealso{ +\code{\link[=prop_diff_cmh]{prop_diff_cmh()}}, \code{\link[=h_prop_cmh]{h_prop_cmh()}} +} +\keyword{internal} diff --git a/man/h_miettinen_nurminen_stratified_ci.Rd b/man/h_miettinen_nurminen_stratified_ci.Rd new file mode 100644 index 0000000000..e5a5bf5690 --- /dev/null +++ b/man/h_miettinen_nurminen_stratified_ci.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prop_diff.R +\name{h_miettinen_nurminen_stratified_ci} +\alias{h_miettinen_nurminen_stratified_ci} +\title{Stratified Miettinen-Nurminen Confidence Interval} +\usage{ +h_miettinen_nurminen_stratified_ci(prop, conf_level = 0.95) +} +\arguments{ +\item{prop}{(\code{list})\cr +A named list returned by \code{\link[=h_prop_cmh]{h_prop_cmh()}}. It must contain the following +atomic vectors: +\describe{ +\item{\code{est1}}{CMH-weighted estimated proportion for group 1. +May be \code{NA_real_} when a CMH-weighted estimate cannot be calculated.} +\item{\code{est2}}{CMH-weighted estimated proportion for group 2. +May be \code{NA_real_} when a CMH-weighted estimate cannot be calculated.} +\item{\code{x1}}{Number of responders in group 1 for each stratum.} +\item{\code{x2}}{Number of responders in group 2 for each stratum.} +\item{\code{n1}}{Number of observations in group 1 for each stratum.} +\item{\code{n2}}{Number of observations in group 2 for each stratum.} +\item{\code{p1}}{Observed response proportion in group 1 for each stratum.} +\item{\code{p2}}{Observed response proportion in group 2 for each stratum.} +\item{\code{w}}{Unnormalized CMH weight for each stratum.} +\item{\code{w_normalized}}{Normalized CMH weight for each stratum.} +}} + +\item{conf_level}{(\code{number(1)}) \cr +Confidence level for the confidence interval.} +} +\value{ +A named list containing: +\describe{ +\item{\code{ci}}{(\code{numeric(2)}) Lower and upper confidence limits for the +stratified difference in proportions.} +\item{\code{se}}{(\code{numeric(1)}) Standard error of the stratified difference +in proportions.} +} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +Calculates the stratified Miettinen-Nurminen confidence interval and +standard error for the difference in proportions. The method uses +constrained maximum likelihood estimates within each stratum and combines +the stratum-specific variance estimates using the normalized CMH weights. +} +\details{ +The difference in proportions is defined as the proportion in group 2 +minus the proportion in group 1. +} +\seealso{ +\code{\link[=prop_diff_cmh]{prop_diff_cmh()}}, \code{\link[=h_prop_cmh]{h_prop_cmh()}}, \code{\link[=h_miettinen_nurminen_var]{h_miettinen_nurminen_var()}} +} +\keyword{internal} diff --git a/man/h_miettinen_nurminen_var.Rd b/man/h_miettinen_nurminen_var.Rd new file mode 100644 index 0000000000..d3c2896d62 --- /dev/null +++ b/man/h_miettinen_nurminen_var.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prop_diff.R +\name{h_miettinen_nurminen_var} +\alias{h_miettinen_nurminen_var} +\title{Variance Estimate Following Miettinen and Nurminen} +\usage{ +h_miettinen_nurminen_var(est1, est2, x1, x2, n1, n2) +} +\arguments{ +\item{est1}{(\code{numeric(1)}) \cr +Estimated proportion for group 1. Used together with \code{est2} to define the +risk difference. +May be \code{NA_real_} when an estimate cannot be calculated.} + +\item{est2}{(\code{numeric(1)}) \cr +Estimated proportion for group 2. Used together with \code{est1} to define the +risk difference. +May be \code{NA_real_} when an estimate cannot be calculated.} + +\item{x1}{(\code{numeric}) \cr +Number of responders in group 1 for each stratum. +Must have length at least 1.} + +\item{x2}{(\code{numeric}) \cr +Number of responders in group 2 for each stratum. +Must have the same length as \code{x1}.} + +\item{n1}{(\code{numeric}) \cr +Number of observations in group 1 for each stratum. +Must have the same length as \code{x1}.} + +\item{n2}{(\code{numeric}) \cr +Number of observations in group 2 for each stratum. +Must have the same length as \code{x1}.} +} +\value{ +A named \code{list} with elements: +\itemize{ +\item \code{p1_est}: constrained maximum likelihood estimate of the proportion in +group 1 for each stratum. +\item \code{p2_est}: constrained maximum likelihood estimate of the proportion in +group 2 for each stratum. +\item \code{var_est}: Miettinen-Nurminen variance estimate for each stratum. +} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} + +Calculates the \insertCite{MiettinenNurminen1985;textual}{tern} variance +estimate for the difference between two proportions. The estimate is based on +the constrained maximum likelihood estimates of the two proportions under the +specified risk difference and is used to obtain the standard error for the +Miettinen-Nurminen confidence interval. +} +\details{ +The risk difference is defined as \code{est2} - \code{est1}. For each stratum, the +function calculates the constrained maximum likelihood estimate for the +proportion in group 1 and obtains the corresponding estimate for group 2 +by adding the risk difference. The variance is then calculated from these +estimates using the Miettinen-Nurminen variance formula. + +The variance is returned as \code{NA_real_} for strata where the variance cannot +be calculated. + +The variable names in this function follow the notation in the original +paper by \insertCite{MiettinenNurminen1985;textual}{tern}, cf. Appendix 1. +} +\references{ +\insertAllCited{} +} +\seealso{ +\code{\link[=prop_diff_cmh]{prop_diff_cmh()}}, \code{\link[=h_prop_cmh]{h_prop_cmh()}}, \code{\link[=h_miettinen_nurminen_stratified_ci]{h_miettinen_nurminen_stratified_ci()}} +} +\keyword{internal} diff --git a/man/h_miettinen_nurminen_var_est.Rd b/man/h_miettinen_nurminen_var_est.Rd deleted file mode 100644 index 447d31b70b..0000000000 --- a/man/h_miettinen_nurminen_var_est.Rd +++ /dev/null @@ -1,36 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prop_diff.R -\name{h_miettinen_nurminen_var_est} -\alias{h_miettinen_nurminen_var_est} -\title{Variance Estimates in Strata following Miettinen and Nurminen} -\usage{ -h_miettinen_nurminen_var_est(n1, n2, x1, x2, diff_par) -} -\arguments{ -\item{n1}{(\code{numeric})\cr sample sizes in group 1.} - -\item{n2}{(\code{numeric})\cr sample sizes in group 2.} - -\item{x1}{(\code{numeric})\cr number of responders in group 1.} - -\item{x2}{(\code{numeric})\cr number of responders in group 2.} - -\item{diff_par}{(\code{numeric})\cr assumed difference in true proportions -(group 2 minus group 1).} -} -\value{ -A named \code{list} with elements: -\itemize{ -\item \code{p1_hat}: estimated proportion in group 1 -\item \code{p2_hat}: estimated proportion in group 2 -\item \code{var_est}: variance estimate of the difference in proportions -} -} -\description{ -The variable names in this function follow the notation in the original -paper by \insertCite{MiettinenNurminen1985;textual}{tern}, cf. Appendix 1. -} -\references{ -\insertAllCited{} -} -\keyword{internal} diff --git a/man/h_prop_cmh.Rd b/man/h_prop_cmh.Rd new file mode 100644 index 0000000000..a5e2e15cb9 --- /dev/null +++ b/man/h_prop_cmh.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prop_diff.R +\name{h_prop_cmh} +\alias{h_prop_cmh} +\title{Helper function to calculate the CMH-weighted proportions and their +confidence intervals.} +\usage{ +h_prop_cmh(tbl, conf_level = 0.95) +} +\arguments{ +\item{tbl}{(\code{array})\cr +A three-dimensional contingency table containing counts for each +combination of group, response, and stratum, in that order. +The first two dimensions must each have exactly two levels, and the second +dimension (response) must have names \code{"TRUE"} and \code{"FALSE"}. +At least one stratum must be present. Strata with all cell counts equal to +zero are allowed. +All cell values must be finite, non-missing integer counts.} + +\item{conf_level}{(\code{number(1)})\cr +Confidence level for the confidence intervals.} +} +\value{ +A named list containing the CMH-weighted proportion estimates, +confidence intervals, and intermediate quantities. +The stratum-specific quantities \code{x1}, \code{n1}, \code{p1}, \code{x2}, \code{n2}, \code{p2}, \code{w}, +and \code{w_normalized} are vectors with a length equal to the number of strata +in \code{tbl} and retain the same stratum order. +Some of these quantities may be \code{NA} for strata where they are not defined. +In particular, \code{p1} or \code{p2} is \code{NA} when the corresponding group has no +observations in that stratum. + +\code{est1} and \code{est2} are the overall CMH-weighted proportion estimates for the +two groups, respectively. \code{est_both_groups} contains these two estimates in +group order. +\code{ci_both_groups} contains the corresponding confidence intervals in the same +group order. + +If no stratum contains observations in both groups, the CMH weights +cannot be normalized and the overall estimates and confidence intervals +are \code{NA} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +} +\seealso{ +\code{\link[=prop_diff_cmh]{prop_diff_cmh()}} +} +\keyword{internal} diff --git a/man/h_prop_diff.Rd b/man/h_prop_diff.Rd index 8305fce438..92d476387b 100644 --- a/man/h_prop_diff.Rd +++ b/man/h_prop_diff.Rd @@ -1,15 +1,13 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/prop_diff.R -\name{h_prop_diff} -\alias{h_prop_diff} +\name{prop_diff_wald} \alias{prop_diff_wald} \alias{prop_diff_ha} \alias{prop_diff_nc} -\alias{h_diff_cmh} -\alias{h_diff_cmh_se} \alias{prop_diff_cmh} \alias{prop_diff_strat_nc} \alias{prop_diff_uncond_exact} +\alias{h_prop_diff} \title{Helper functions to calculate proportion difference} \usage{ prop_diff_wald(rsp, grp, conf_level = 0.95, correct = FALSE) @@ -18,10 +16,6 @@ prop_diff_ha(rsp, grp, conf_level) prop_diff_nc(rsp, grp, conf_level, correct = FALSE) -h_diff_cmh(tbl) - -h_diff_cmh_se(cmh_results, diff_se = c("standard", "sato")) - prop_diff_cmh( rsp, grp, @@ -52,18 +46,12 @@ prop_diff_uncond_exact(rsp, grp, conf_level = 0.95) \item{correct}{(\code{flag})\cr whether to include the continuity correction. For further information, see \code{\link[stats:prop.test]{stats::prop.test()}}.} -\item{tbl}{(\code{array})\cr 3-dimensional array with dimensions corresponding to -group, response, and strata. The second dimension (response) should have names -"TRUE" and "FALSE".} - -\item{cmh_results}{(\code{list})\cr output of \code{\link[=h_diff_cmh]{h_diff_cmh()}}.} +\item{strata}{(\code{factor})\cr variable with one level per stratum and same length as \code{rsp}.} -\item{diff_se}{(\code{string})\cr method to estimate the standard error for the difference, either -\code{standard}, \code{sato} \insertCite{Sato1989}{tern} or +\item{diff_se}{(\code{string})\cr method to estimate the standard error for the +difference, either \code{standard}, \code{sato} \insertCite{Sato1989}{tern} or \code{miettinen_nurminen} \insertCite{MiettinenNurminen1985}{tern}.} -\item{strata}{(\code{factor})\cr variable with one level per stratum and same length as \code{rsp}.} - \item{weights_method}{(\code{string})\cr method used to estimate the weights for stratified Newcombe method. Must be either \code{"cmh"} or \code{"wilson_h"}. \code{"cmh"} uses weights derived from @@ -87,17 +75,13 @@ interval. \item \code{prop_diff_ha()}: Anderson-Hauck confidence interval \insertCite{HauckAnderson1986}{tern}. \item \code{prop_diff_nc()}: Newcombe confidence interval. It is based on -the Wilson score confidence interval for a single binomial proportion \insertCite{Newcombe1998}{tern}. - -\item \code{h_diff_cmh()}: Helper function to calculate the CMH weighted -difference in proportions. - -\item \code{h_diff_cmh_se()}: Helper function to calculate the standard error for the -CMH weighted difference in proportions. +the Wilson score confidence interval for a single binomial proportion +\insertCite{Newcombe1998}{tern}. -\item \code{prop_diff_cmh()}: Calculates the weighted difference. This is defined as the difference in -response rates between the experimental treatment group and the control treatment group, adjusted -for stratification factors by applying Cochran-Mantel-Haenszel (CMH) weights. For the CMH chi-squared +\item \code{prop_diff_cmh()}: Calculates the weighted difference. This is defined +as the difference in response rates between the experimental treatment +group and the control treatment group, adjusted for stratification factors +by applying Cochran-Mantel-Haenszel (CMH) weights. For the CMH chi-squared test, use \code{\link[stats:mantelhaen.test]{stats::mantelhaen.test()}}. \item \code{prop_diff_strat_nc()}: Calculates the stratified Newcombe confidence interval and difference in response @@ -205,4 +189,3 @@ prop_diff_uncond_exact(rsp = rsp, grp = grp, conf_level = 0.95) \seealso{ \code{\link[=prop_diff]{prop_diff()}} for implementation of these helper functions. } -\keyword{internal} diff --git a/man/tern-package.Rd b/man/tern-package.Rd index 6d1ba5f2ef..41599e6ec5 100644 --- a/man/tern-package.Rd +++ b/man/tern-package.Rd @@ -39,6 +39,7 @@ Authors: Other contributors: \itemize{ \item David Munoz Tord \email{david.munoztord@mailbox.org} [contributor] + \item Wojciech Wojciak \email{wojciech.wojciak@gmail.com} [contributor] \item F. Hoffmann-La Roche AG [copyright holder, funder] } diff --git a/man/uniroot_catch_na.Rd b/man/uniroot_catch_na.Rd new file mode 100644 index 0000000000..999f9d2a93 --- /dev/null +++ b/man/uniroot_catch_na.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{uniroot_catch_na} +\alias{uniroot_catch_na} +\title{Find a root while returning NA for missing function values} +\usage{ +uniroot_catch_na(f, interval, ...) +} +\arguments{ +\item{f}{(\code{function})\cr function for which the root is sought.} + +\item{interval}{(\code{numeric(2)})\cr end points of the interval to be searched.} + +\item{...}{further arguments passed to \code{\link[stats:uniroot]{stats::uniroot()}}.} +} +\value{ +A numeric scalar containing the found root, or \code{NA_real_} if \code{f} returns +\code{NA} at either endpoint of \code{interval}. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +A wrapper around \code{\link[stats:uniroot]{stats::uniroot()}} that returns \code{NA_real_} when \code{f} is \code{NA} +at either end of \code{interval}. +} +\seealso{ +\code{\link[stats:uniroot]{stats::uniroot()}} +} +\keyword{internal} diff --git a/tests/testthat/_snaps/prop_diff.md b/tests/testthat/_snaps/prop_diff.md index 4a75a6ce3d..db555e2467 100644 --- a/tests/testthat/_snaps/prop_diff.md +++ b/tests/testthat/_snaps/prop_diff.md @@ -1,4 +1,4 @@ -# `prop_diff_ha` (proportion difference by Anderson-Hauck) +# prop_diff_ha (proportion difference by Anderson-Hauck) Code res @@ -22,7 +22,7 @@ [1] -0.8451161 0.8451161 -# `prop_diff_nc` (proportion difference by Newcombe) +# prop_diff_nc (proportion difference by Newcombe) Code res @@ -46,7 +46,7 @@ [1] -0.361619 0.361619 -# `prop_diff_wald` (proportion difference by Wald's test: with correction) +# prop_diff_wald (proportion difference by Wald's test: with correction) Code res @@ -82,7 +82,7 @@ [1] -0.375 0.375 -# `prop_diff_wald` (proportion difference by Wald's test: without correction) +# prop_diff_wald (proportion difference by Wald's test: without correction) Code res @@ -118,7 +118,7 @@ [1] 0 0 -# `prop_diff_cmh` (proportion difference by CMH) +# prop_diff_cmh (proportion difference by CMH) Code res @@ -157,7 +157,7 @@ 8 9 4 9 6 6 -# `prop_diff_cmh` with Sato variance estimator for difference +# prop_diff_cmh with Sato variance estimator for difference Code res @@ -196,17 +196,6 @@ 8 9 4 9 6 6 -# h_miettinen_nurminen_var_est works as expected - - list(p1_hat = 0.342213591803752, p2_hat = 0.442213591803752, - var_est = 0.0405774934104561) - ---- - - list(p1_hat = c(0.342213591803752, 0.265846883932378), p2_hat = c(0.442213591803752, - 0.365846883932378), var_est = c(0.0405774934104561, 0.0301587022300622 - )) - # prop_diff_cmh works correctly when some strata don't have both groups Code @@ -234,16 +223,16 @@ [1] 0.09489839 $weights - b.x a.y b.y a.z b.z - 0.2408257 0.1297378 0.2408257 0.1997279 0.1888829 + a.x b.x a.y b.y a.z b.z + 0.0000000 0.2408257 0.1297378 0.2408257 0.1997279 0.1888829 $n1 - b.x a.y b.y a.z b.z - 11 8 11 13 11 + a.x b.x a.y b.y a.z b.z + 12 11 8 11 13 11 $n2 - b.x a.y b.y a.z b.z - 9 4 9 6 6 + a.x b.x a.y b.y a.z b.z + 0 9 4 9 6 6 # prop_diff_cmh works correctly when strata combinations are empty @@ -273,16 +262,16 @@ [1] 0.09489839 $weights - b.x a.y b.y a.z b.z - 0.2408257 0.1297378 0.2408257 0.1997279 0.1888829 + a.x b.x a.y b.y a.z b.z + NA 0.2408257 0.1297378 0.2408257 0.1997279 0.1888829 $n1 - b.x a.y b.y a.z b.z - 11 8 11 13 11 + a.x b.x a.y b.y a.z b.z + 0 11 8 11 13 11 $n2 - b.x a.y b.y a.z b.z - 9 4 9 6 6 + a.x b.x a.y b.y a.z b.z + 0 9 4 9 6 6 # prop_diff_strat_nc output matches equivalent SAS function output @@ -293,7 +282,1196 @@ value lower upper 0.25390590 0.03467969 0.44544132 -# `estimate_proportion_diff` is compatible with `rtables` +# h_prop_cmh works as expected with non-sparse tables + + Code + res + Output + $x1 + S1 S2 S3 + 12 15 9 + + $n1 + S1 S2 S3 + 30 35 23 + + $p1 + S1 S2 S3 + 0.4000000 0.4285714 0.3913043 + + $x2 + S1 S2 S3 + 8 10 6 + + $n2 + S1 S2 S3 + 30 35 27 + + $p2 + S1 S2 S3 + 0.2666667 0.2857143 0.2222222 + + $w + S1 S2 S3 + 15.00 17.50 12.42 + + $w_normalized + S1 S2 S3 + 0.3339270 0.3895815 0.2764915 + + $est1 + ref + 0.4087266 + + $est2 + Not-ref + 0.2617988 + + $est_both_groups + ref Not-ref + 0.4087266 0.2617988 + + $var1 + ref + 0.002745713 + + $var2 + Not-ref + 0.002101216 + + $ci_both_groups + $ci_both_groups$ref + [1] 0.3060254 0.5114279 + + $ci_both_groups$`Not-ref` + [1] 0.1719559 0.3516416 + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 1 42 0 + + $n1 + S1 S2 S3 + 8 50 95 + + $p1 + S1 S2 S3 + 0.125 0.840 0.000 + + $x2 + S1 S2 S3 + 0 3 2 + + $n2 + S1 S2 S3 + 13 70 13 + + $p2 + S1 S2 S3 + 0.00000000 0.04285714 0.15384615 + + $w + S1 S2 S3 + 4.952381 29.166667 11.435185 + + $w_normalized + S1 S2 S3 + 0.1087140 0.6402625 0.2510235 + + $est1 + ref + 0.5514097 + + $est2 + Not-ref + 0.06605883 + + $est_both_groups + ref Not-ref + 0.55140974 0.06605883 + + $var1 + ref + 0.001263492 + + $var2 + Not-ref + 0.0008712136 + + $ci_both_groups + $ci_both_groups$ref + [1] 0.4817416 0.6210779 + + $ci_both_groups$`Not-ref` + [1] 0.00820789 0.12390977 + + + +--- + + Code + res + Output + $x1 + S1 + 1 + + $n1 + S1 + 8 + + $p1 + S1 + 0.125 + + $x2 + S1 + 0 + + $n2 + S1 + 13 + + $p2 + S1 + 0 + + $w + S1 + 4.952381 + + $w_normalized + S1 + 1 + + $est1 + ref + 0.125 + + $est2 + Not-ref + 0 + + $est_both_groups + ref Not-ref + 0.125 0.000 + + $var1 + ref + 0.01367188 + + $var2 + Not-ref + 0 + + $ci_both_groups + $ci_both_groups$ref + [1] -0.1041723 0.3541723 + + $ci_both_groups$`Not-ref` + [1] 0 0 + + + +# h_prop_cmh handles empty and sparse contingency tables + + Code + res + Output + $x1 + S1 S2 S3 + 0 0 0 + + $n1 + S1 S2 S3 + 0 0 0 + + $p1 + S1 S2 S3 + NA NA NA + + $x2 + S1 S2 S3 + 0 0 0 + + $n2 + S1 S2 S3 + 0 0 0 + + $p2 + S1 S2 S3 + NA NA NA + + $w + S1 S2 S3 + NA NA NA + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 1 0 0 + + $n1 + S1 S2 S3 + 1 0 0 + + $p1 + S1 S2 S3 + 1 NA NA + + $x2 + S1 S2 S3 + 0 0 0 + + $n2 + S1 S2 S3 + 0 0 0 + + $p2 + S1 S2 S3 + NA NA NA + + $w + S1 S2 S3 + 0 NA NA + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 0 1 0 + + $n1 + S1 S2 S3 + 0 1 0 + + $p1 + S1 S2 S3 + NA 1 NA + + $x2 + S1 S2 S3 + 0 0 0 + + $n2 + S1 S2 S3 + 0 0 0 + + $p2 + S1 S2 S3 + NA NA NA + + $w + S1 S2 S3 + NA 0 NA + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 0 0 0 + + $n1 + S1 S2 S3 + 3 0 0 + + $p1 + S1 S2 S3 + 0 NA NA + + $x2 + S1 S2 S3 + 0 0 0 + + $n2 + S1 S2 S3 + 0 0 0 + + $p2 + S1 S2 S3 + NA NA NA + + $w + S1 S2 S3 + 0 NA NA + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 0 0 0 + + $n1 + S1 S2 S3 + 0 0 1 + + $p1 + S1 S2 S3 + NA NA 0 + + $x2 + S1 S2 S3 + 0 0 0 + + $n2 + S1 S2 S3 + 0 0 0 + + $p2 + S1 S2 S3 + NA NA NA + + $w + S1 S2 S3 + NA NA 0 + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 4 0 0 + + $n1 + S1 S2 S3 + 4 0 0 + + $p1 + S1 S2 S3 + 1 NA NA + + $x2 + S1 S2 S3 + 0 7 0 + + $n2 + S1 S2 S3 + 0 7 0 + + $p2 + S1 S2 S3 + NA 1 NA + + $w + S1 S2 S3 + 0 0 NA + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 0 0 0 + + $n1 + S1 S2 S3 + 9 0 0 + + $p1 + S1 S2 S3 + 0 NA NA + + $x2 + S1 S2 S3 + 0 0 0 + + $n2 + S1 S2 S3 + 0 12 0 + + $p2 + S1 S2 S3 + NA 0 NA + + $w + S1 S2 S3 + 0 0 NA + + $w_normalized + S1 S2 S3 + NA NA NA + + $est1 + ref + NA + + $est2 + Not-ref + NA + + $est_both_groups + ref Not-ref + NA NA + + $var1 + ref + NA + + $var2 + Not-ref + NA + + $ci_both_groups + $ci_both_groups$ref + [1] NA NA + + $ci_both_groups$`Not-ref` + [1] NA NA + + + +--- + + Code + res + Output + $x1 + S1 S2 S3 + 4 0 10 + + $n1 + S1 S2 S3 + 4 0 22 + + $p1 + S1 S2 S3 + 1.0000000 NA 0.4545455 + + $x2 + S1 S2 S3 + 0 0 40 + + $n2 + S1 S2 S3 + 0 0 83 + + $p2 + S1 S2 S3 + NA NA 0.4819277 + + $w + S1 S2 S3 + 0.00000 NA 17.39048 + + $w_normalized + S1 S2 S3 + 0 NA 1 + + $est1 + ref + 0.4545455 + + $est2 + Not-ref + 0.4819277 + + $est_both_groups + ref Not-ref + 0.4545455 0.4819277 + + $var1 + ref + 0.01126972 + + $var2 + Not-ref + 0.003008113 + + $ci_both_groups + $ci_both_groups$ref + [1] 0.2464777 0.6626132 + + $ci_both_groups$`Not-ref` + [1] 0.3744310 0.5894244 + + + +# h_prop_cmh respects a custom confidence level + + Code + res + Output + $x1 + S1 S2 S3 + 12 15 9 + + $n1 + S1 S2 S3 + 30 35 23 + + $p1 + S1 S2 S3 + 0.4000000 0.4285714 0.3913043 + + $x2 + S1 S2 S3 + 8 10 6 + + $n2 + S1 S2 S3 + 30 35 27 + + $p2 + S1 S2 S3 + 0.2666667 0.2857143 0.2222222 + + $w + S1 S2 S3 + 15.00 17.50 12.42 + + $w_normalized + S1 S2 S3 + 0.3339270 0.3895815 0.2764915 + + $est1 + ref + 0.4087266 + + $est2 + Not-ref + 0.2617988 + + $est_both_groups + ref Not-ref + 0.4087266 0.2617988 + + $var1 + ref + 0.002745713 + + $var2 + Not-ref + 0.002101216 + + $ci_both_groups + $ci_both_groups$ref + [1] 0.3415739 0.4758794 + + $ci_both_groups$`Not-ref` + [1] 0.2030537 0.3205438 + + + +# h_cmh_sato_var works as expected with non-sparse tables + + 0.00484709089617089 + +--- + + 0.00301389352651087 + +--- + + 0.013671875 + +# h_cmh_sato_var empty and sparse contingency tables + + NA_real_ + +--- + + NA_real_ + +--- + + NA_real_ + +--- + + NA_real_ + +--- + + NA_real_ + +--- + + NA_real_ + +--- + + NA_real_ + +--- + + 0.0142778351745434 + +# h_miettinen_nurminen_var works as expected with non-sparse tables + + Code + res1 + Output + $p1_est + S1 S2 S3 + 0.4075605 0.4307970 0.3778246 + + $p2_est + S1 S2 S3 + 0.2606327 0.2838692 0.2308967 + + $var_est + S1 S2 S3 + 0.01471723 0.01299995 0.01714055 + + +--- + + Code + res2 + Output + $p1_est + S1 S2 S3 + 0.4853509 0.6010255 0.4953824 + + $p2_est + S1 S2 S3 + 0.00000000 0.11567462 0.01003146 + + $var_est + S1 S2 S3 + 0.032784334 0.006309801 0.003426996 + + +--- + + Code + res3 + Output + $p1_est + S1 + 0.125 + + $p2_est + S1 + 1.110223e-16 + + $var_est + S1 + 0.01435547 + + +# h_miettinen_nurminen_var empty and sparse contingency tables + + Code + res1 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res2 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res3 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res4 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res5 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res6 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res7 + Output + $p1_est + S1 S2 S3 + NA NA NA + + $p2_est + S1 S2 S3 + NA NA NA + + $var_est + S1 S2 S3 + NA NA NA + + +--- + + Code + res8 + Output + $p1_est + S1 S2 S3 + 0.9726177 NaN 0.4545455 + + $p2_est + S1 S2 S3 + 1.0000000 NaN 0.4819277 + + $var_est + S1 S2 S3 + NA NA 0.01441512 + + +# h_miettinen_nurminen_var works as expected + + list(p1_est = 0.342213591803752, p2_est = 0.442213591803752, + var_est = 0.0405774934104561) + +--- + + list(p1_est = c(0.342213591803752, 0.265846883932378), p2_est = c(0.442213591803752, + 0.365846883932378), var_est = c(0.0405774934104561, 0.0301587022300622 + )) + +# h_miettinen_nurminen_stratified_ci works as expected with non-sparse tables + + Code + res1 + Output + $ci + [1] -0.281568259 -0.008104148 + + $se + [1] 0.07017465 + + +--- + + Code + res2 + Output + $ci + [1] -0.5899410 -0.3703714 + + $se + [1] 0.05648034 + + +--- + + Code + res3 + Output + $ci + [1] -0.4797396 0.1311335 + + $se + [1] 0.1198143 + + +# h_miettinen_nurminen_stratified_ci empty and sparse contingency tables + + Code + res1 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res2 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res3 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res4 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res5 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res6 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res7 + Output + $ci + [1] NA NA + + $se + [1] NA + + +--- + + Code + res8 + Output + $ci + [1] -0.2015304 0.2461176 + + $se + [1] 0.120063 + + +# h_miettinen_nurminen_stratified_ci respects a custom confidence level + + Code + res + Output + $ci + [1] -0.2357582 -0.0563385 + + $se + [1] 0.07017465 + + +# estimate_proportion_diff is compatible with rtables Code res @@ -303,7 +1481,7 @@ Difference in Response rate (%) 25.0 90% CI (Anderson-Hauck) (-92.0, 100.0) -# `estimate_proportion_diff` and cmh is compatible with `rtables` +# estimate_proportion_diff and cmh is compatible with rtables Code res @@ -452,7 +1630,7 @@ [1] "Difference in Response rate (%) and 95% CI (Wald, without correction)" -# `estimate_proportion_diff` with diff_est_ci builds single-row table +# estimate_proportion_diff with diff_est_ci builds single-row table Code res diff --git a/tests/testthat/_snaps/test_proportion_diff.md b/tests/testthat/_snaps/test_proportion_diff.md index e0002a3519..d7280b7b4b 100644 --- a/tests/testthat/_snaps/test_proportion_diff.md +++ b/tests/testthat/_snaps/test_proportion_diff.md @@ -292,8 +292,8 @@ Code res Output - B A - ———————————————————————————————————————————————————————————————————————————————————— - Variable Label - p-value (Cochran-Mantel-Haenszel Test with Sato Variance Estimator) 1.0000 + B A + ———————————————————————————————————————————————————————————————————————————————— + Variable Label + p-value (Cochran-Mantel-Haenszel Test with Sato Variance Estimator) NA diff --git a/tests/testthat/helper-prop_diff.R b/tests/testthat/helper-prop_diff.R new file mode 100644 index 0000000000..77dd58320d --- /dev/null +++ b/tests/testthat/helper-prop_diff.R @@ -0,0 +1,102 @@ +h_get_prop_data <- function(sparse = FALSE) { + checkmate::assert_flag(sparse) + + dimnames <- list(grp = c("ref", "Not-ref"), rsp = c("TRUE", "FALSE"), strata = c("S1", "S2", "S3")) + + tables <- if (!sparse) { + list( + tbl1 = array( + c( + 12, 8, 18, 22, # S1 + 15, 10, 20, 25, # S2 + 9, 6, 14, 21 # S3 + ), + dim = c(2L, 2L, 3L), + dimnames = dimnames + ), + tbl2 = array( + c( + 1, 0, 7, 13, # S1 + 42, 3, 8, 67, # S2 + 0, 2, 95, 11 # S3 + ), + dim = c(2L, 2L, 3L), + dimnames = dimnames + ), + tbl3 = array( + c(1, 0, 7, 13), # S1 + dim = c(2L, 2L, 1L), + dimnames = list(grp = c("ref", "Not-ref"), rsp = c("TRUE", "FALSE"), strata = "S1") + ) + ) + } else { + list( + tbl1 = array( + c( + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl2 = array( + c( + 1, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl3 = array( + c( + 0, 0, 0, 0, + 1, 0, 0, 0, + 0, 0, 0, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl4 = array( + c( + 0, 0, 3, 0, + 0, 0, 0, 0, + 0, 0, 0, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl5 = array( + c( + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 1, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl6 = array( + c( + 4, 0, 0, 0, + 0, 7, 0, 0, + 0, 0, 0, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl7 = array( + c( + 0, 0, 9, 0, + 0, 0, 0, 12, + 0, 0, 0, 0 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ), + tbl8 = array( + c( + 4, 0, 0, 0, + 0, 0, 0, 0, + 10, 40, 12, 43 + ), + dim = c(2L, 2L, 3L), dimnames = dimnames + ) + ) + } + + tables +} diff --git a/tests/testthat/test-prop_diff.R b/tests/testthat/test-prop_diff.R index afaff53951..dbf949e1e4 100644 --- a/tests/testthat/test-prop_diff.R +++ b/tests/testthat/test-prop_diff.R @@ -1,4 +1,4 @@ -testthat::test_that("`prop_diff_ha` (proportion difference by Anderson-Hauck)", { +test_that("prop_diff_ha (proportion difference by Anderson-Hauck)", { # "Mid" case: 3/4 respond in group A, 1/2 respond in group B. rsp <- c(TRUE, FALSE, FALSE, TRUE, TRUE, TRUE) grp <- factor(c("A", "B", "A", "B", "A", "A"), levels = c("B", "A")) @@ -17,7 +17,7 @@ testthat::test_that("`prop_diff_ha` (proportion difference by Anderson-Hauck)", testthat::expect_snapshot(res) }) -testthat::test_that("`prop_diff_nc` (proportion difference by Newcombe)", { +testthat::test_that("prop_diff_nc (proportion difference by Newcombe)", { # "Mid" case: 3/4 respond in group A, 1/2 respond in group B. rsp <- c(TRUE, FALSE, FALSE, TRUE, TRUE, TRUE) grp <- factor(c("A", "B", "A", "B", "A", "A"), levels = c("B", "A")) @@ -38,7 +38,7 @@ testthat::test_that("`prop_diff_nc` (proportion difference by Newcombe)", { testthat::expect_snapshot(res) }) -testthat::test_that("`prop_diff_wald` (proportion difference by Wald's test: with correction)", { +testthat::test_that("prop_diff_wald (proportion difference by Wald's test: with correction)", { # "Mid" case: 3/4 respond in group A, 1/2 respond in group B. rsp <- c(TRUE, FALSE, FALSE, TRUE, TRUE, TRUE) grp <- factor(c("A", "B", "A", "B", "A", "A"), levels = c("B", "A")) @@ -68,7 +68,7 @@ testthat::test_that("`prop_diff_wald` (proportion difference by Wald's test: wit testthat::expect_snapshot(res) }) -testthat::test_that("`prop_diff_wald` (proportion difference by Wald's test: without correction)", { +testthat::test_that("prop_diff_wald (proportion difference by Wald's test: without correction)", { # "Mid" case: 3/4 respond in group A, 1/2 respond in group B. rsp <- c(TRUE, FALSE, FALSE, TRUE, TRUE, TRUE) grp <- factor(c("A", "B", "A", "B", "A", "A"), levels = c("B", "A")) @@ -100,7 +100,7 @@ testthat::test_that("`prop_diff_wald` (proportion difference by Wald's test: wit testthat::expect_snapshot(res) }) -testthat::test_that("`prop_diff_cmh` (proportion difference by CMH)", { +testthat::test_that("prop_diff_cmh (proportion difference by CMH)", { set.seed(2, kind = "Mersenne-Twister") rsp <- sample(c(TRUE, FALSE), 100, TRUE) grp <- sample(c("Placebo", "Treatment"), 100, TRUE) @@ -124,7 +124,7 @@ testthat::test_that("`prop_diff_cmh` (proportion difference by CMH)", { )) }) -testthat::test_that("`prop_diff_cmh` with Sato variance estimator for difference", { +testthat::test_that("prop_diff_cmh with Sato variance estimator for difference", { set.seed(2, kind = "Mersenne-Twister") rsp <- sample(c(TRUE, FALSE), 100, TRUE) grp <- sample(c("Placebo", "Treatment"), 100, TRUE) @@ -148,20 +148,6 @@ testthat::test_that("`prop_diff_cmh` with Sato variance estimator for difference testthat::expect_snapshot(res) }) -testthat::test_that("h_miettinen_nurminen_var_est works as expected", { - result <- h_miettinen_nurminen_var_est( - n1 = 10, n2 = 15, - x1 = 4, x2 = 6, diff_par = 0.1 - ) - expect_snapshot_value(result, style = "deparse", tolerance = 1e-4) - - result2 <- h_miettinen_nurminen_var_est( - n1 = c(10, 12), n2 = c(15, 18), - x1 = c(4, 2), x2 = c(6, 8), diff_par = 0.1 - ) - expect_snapshot_value(result2, style = "deparse", tolerance = 1e-4) -}) - testthat::test_that("prop_diff_cmh works correctly with Miettinen-Nurminen variance estimator", { # Example from the Lu (2008) paper, described in Melikov and Mosier (2025), # https://pharmasug.org/proceedings/2025/SA/PharmaSUG-2025-SA-198.pdf @@ -260,7 +246,20 @@ testthat::test_that("prop_diff_cmh works correctly when strata combinations are testthat::expect_snapshot(res) }) -testthat::test_that("`prop_strat_nc` (proportion difference by stratified Newcombe) with cmh weights", { +testthat::test_that("prop_diff_cmh ignores strata with only one group", { + rsp <- rep(rep(c(TRUE, FALSE), 5), c(6, 4, 3, 7, 5, 5, 8, 2, 4, 2)) + grp <- factor(rep(c("a", "b", "a", "b", "a"), c(10, 10, 10, 10, 6))) + strata <- factor(rep(c("s1", "s1", "s2", "s2", "s3"), c(10, 10, 10, 10, 6))) + keep <- strata != "s3" # s3 has group "a" only + + for (diff_se in c("standard", "sato", "miettinen_nurminen")) { + res <- prop_diff_cmh(rsp, grp, strata, diff_se = diff_se) + res_keep <- prop_diff_cmh(rsp[keep], grp[keep], droplevels(strata[keep]), diff_se = diff_se) + testthat::expect_equal(res[c("diff", "diff_ci", "se_diff")], res_keep[c("diff", "diff_ci", "se_diff")]) + } +}) + +testthat::test_that("prop_strat_nc (proportion difference by stratified Newcombe) with cmh weights", { set.seed(1) rsp <- c( sample(c(TRUE, FALSE), size = 40, prob = c(3 / 4, 1 / 4), replace = TRUE), @@ -286,7 +285,7 @@ testthat::test_that("`prop_strat_nc` (proportion difference by stratified Newcom expect_equal(as.numeric(results$diff_ci), c(0.0347, 0.4454), tolerance = 1e-3) }) -testthat::test_that("`prop_strat_nc` (proportion difference by stratified Newcombe) with wilson_h weights", { +testthat::test_that("prop_strat_nc (proportion difference by stratified Newcombe) with wilson_h weights", { set.seed(1) rsp <- c( sample(c(TRUE, FALSE), size = 40, prob = c(3 / 4, 1 / 4), replace = TRUE), @@ -336,7 +335,7 @@ testthat::test_that("prop_diff_strat_nc output matches equivalent SAS function o }) -testthat::test_that("`prop_diff_uncond_exact` matches reference values and works with edge cases", { +testthat::test_that("prop_diff_uncond_exact matches reference values and works with edge cases", { mk_data <- function(n11, n21, n1, n2) { rsp <- c(rep(TRUE, n21), rep(FALSE, n2 - n21), rep(TRUE, n11), rep(FALSE, n1 - n11)) grp <- factor(c(rep("B", n2), rep("A", n1)), levels = c("B", "A")) @@ -392,6 +391,170 @@ testthat::test_that("`prop_diff_uncond_exact` matches reference values and works ) }) +test_that("h_prop_cmh works as expected with non-sparse tables", { + tables <- h_get_prop_data() + + for (tbl in tables) { + res <- h_prop_cmh(tbl) + expect_snapshot(res) + } +}) + +test_that("h_prop_cmh handles empty and sparse contingency tables", { + tables <- h_get_prop_data(sparse = TRUE) + + for (tbl in tables) { + res <- h_prop_cmh(tbl) + expect_snapshot(res) + } +}) + +test_that("h_prop_cmh respects a custom confidence level", { + tables <- h_get_prop_data() + + res <- h_prop_cmh(tables$tbl1, conf_level = 0.8) + expect_snapshot(res) +}) + +test_that("h_cmh_sato_var works as expected with non-sparse tables", { + tables <- h_get_prop_data() + + res1 <- h_cmh_sato_var(h_prop_cmh(tables$tbl1)) + res2 <- h_cmh_sato_var(h_prop_cmh(tables$tbl2)) + res3 <- h_cmh_sato_var(h_prop_cmh(tables$tbl3)) + + expect_snapshot_value(res1, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res2, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res3, style = "deparse", tolerance = 1e-6) +}) + +test_that("h_cmh_sato_var empty and sparse contingency tables", { + tables <- h_get_prop_data(sparse = TRUE) + + res1 <- h_cmh_sato_var(h_prop_cmh(tables$tbl1)) + res2 <- h_cmh_sato_var(h_prop_cmh(tables$tbl2)) + res3 <- h_cmh_sato_var(h_prop_cmh(tables$tbl3)) + res4 <- h_cmh_sato_var(h_prop_cmh(tables$tbl4)) + res5 <- h_cmh_sato_var(h_prop_cmh(tables$tbl5)) + res6 <- h_cmh_sato_var(h_prop_cmh(tables$tbl6)) + res7 <- h_cmh_sato_var(h_prop_cmh(tables$tbl7)) + res8 <- h_cmh_sato_var(h_prop_cmh(tables$tbl8)) + + expect_snapshot_value(res1, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res2, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res3, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res4, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res5, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res6, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res7, style = "deparse", tolerance = 1e-6) + expect_snapshot_value(res8, style = "deparse", tolerance = 1e-6) +}) + +test_that("h_miettinen_nurminen_var works as expected with non-sparse tables", { + tables <- h_get_prop_data() + t1 <- h_prop_cmh(tables$tbl1) + t2 <- h_prop_cmh(tables$tbl2) + t3 <- h_prop_cmh(tables$tbl3) + + res1 <- h_miettinen_nurminen_var(t1$est1, t1$est2, t1$x1, t1$x2, t1$n1, t1$n2) + res2 <- h_miettinen_nurminen_var(t2$est1, t2$est2, t2$x1, t2$x2, t2$n1, t2$n2) + res3 <- h_miettinen_nurminen_var(t3$est1, t3$est2, t3$x1, t3$x2, t3$n1, t3$n2) + + expect_snapshot(res1) + expect_snapshot(res2) + expect_snapshot(res3) +}) + +test_that("h_miettinen_nurminen_var empty and sparse contingency tables", { + tables <- h_get_prop_data(sparse = TRUE) + t1 <- h_prop_cmh(tables$tbl1) + t2 <- h_prop_cmh(tables$tbl2) + t3 <- h_prop_cmh(tables$tbl3) + t4 <- h_prop_cmh(tables$tbl4) + t5 <- h_prop_cmh(tables$tbl5) + t6 <- h_prop_cmh(tables$tbl6) + t7 <- h_prop_cmh(tables$tbl7) + t8 <- h_prop_cmh(tables$tbl8) + + res1 <- h_miettinen_nurminen_var(t1$est1, t1$est2, t1$x1, t1$x2, t1$n1, t1$n2) + res2 <- h_miettinen_nurminen_var(t2$est1, t2$est2, t2$x1, t2$x2, t2$n1, t2$n2) + res3 <- h_miettinen_nurminen_var(t3$est1, t3$est2, t3$x1, t3$x2, t3$n1, t3$n2) + res4 <- h_miettinen_nurminen_var(t4$est1, t4$est2, t4$x1, t4$x2, t4$n1, t4$n2) + res5 <- h_miettinen_nurminen_var(t5$est1, t5$est2, t5$x1, t5$x2, t5$n1, t5$n2) + res6 <- h_miettinen_nurminen_var(t6$est1, t6$est2, t6$x1, t6$x2, t6$n1, t6$n2) + res7 <- h_miettinen_nurminen_var(t7$est1, t7$est2, t7$x1, t7$x2, t7$n1, t7$n2) + res8 <- h_miettinen_nurminen_var(t8$est1, t8$est2, t8$x1, t8$x2, t8$n1, t8$n2) + + expect_snapshot(res1) + expect_snapshot(res2) + expect_snapshot(res3) + expect_snapshot(res4) + expect_snapshot(res5) + expect_snapshot(res6) + expect_snapshot(res7) + expect_snapshot(res8) +}) + +testthat::test_that("h_miettinen_nurminen_var works as expected", { + result <- h_miettinen_nurminen_var( + est2 = 0.2, est1 = 0.1, + x1 = 4, x2 = 6, + n1 = 10, n2 = 15 + ) + expect_snapshot_value(result, style = "deparse", tolerance = 1e-4) + + result2 <- h_miettinen_nurminen_var( + est2 = 0.2, est1 = 0.1, + x1 = c(4, 2), x2 = c(6, 8), + n1 = c(10, 12), n2 = c(15, 18) + ) + expect_snapshot_value(result2, style = "deparse", tolerance = 1e-4) +}) + +test_that("h_miettinen_nurminen_stratified_ci works as expected with non-sparse tables", { + tables <- h_get_prop_data() + + res1 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl1)) + res2 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl2)) + res3 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl3)) + + expect_snapshot(res1) + expect_snapshot(res2) + expect_snapshot(res3) +}) + +test_that("h_miettinen_nurminen_stratified_ci empty and sparse contingency tables", { + tables <- h_get_prop_data(sparse = TRUE) + + res1 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl1)) + res2 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl2)) + res3 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl3)) + res4 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl4)) + res5 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl5)) + res6 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl6)) + res7 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl7)) + res8 <- h_miettinen_nurminen_stratified_ci(h_prop_cmh(tables$tbl8)) + + expect_snapshot(res1) + expect_snapshot(res2) + expect_snapshot(res3) + expect_snapshot(res4) + expect_snapshot(res5) + expect_snapshot(res6) + expect_snapshot(res7) + expect_snapshot(res8) +}) + +test_that("h_miettinen_nurminen_stratified_ci respects a custom confidence level", { + tables <- h_get_prop_data() + + res <- h_miettinen_nurminen_stratified_ci( + h_prop_cmh(tables$tbl1), + conf_level = 0.8 + ) + expect_snapshot(res) +}) + testthat::test_that("h_worst_case_tail_probability returns valid tail probabilities", { n1 <- 2 n2 <- 2 @@ -508,7 +671,7 @@ test_that("d_proportion_diff returns correct descriptions", { ) }) -testthat::test_that("`estimate_proportion_diff` is compatible with `rtables`", { +testthat::test_that("estimate_proportion_diff is compatible with rtables", { # "Mid" case: 3/4 respond in group A, 1/2 respond in group B. dta <- data.frame( rsp = c(TRUE, FALSE, FALSE, TRUE, TRUE, TRUE), @@ -529,7 +692,7 @@ testthat::test_that("`estimate_proportion_diff` is compatible with `rtables`", { testthat::expect_snapshot(res) }) -testthat::test_that("`estimate_proportion_diff` and cmh is compatible with `rtables`", { +testthat::test_that("estimate_proportion_diff and cmh is compatible with rtables", { set.seed(1) nex <- 100 # Number of test rows dta <- data.frame( @@ -555,7 +718,7 @@ testthat::test_that("`estimate_proportion_diff` and cmh is compatible with `rtab testthat::expect_snapshot(res) }) -testthat::test_that("`estimate_proportion_diff` and strat_newcombe is compatible with `rtables`", { +testthat::test_that("estimate_proportion_diff and strat_newcombe is compatible with rtables", { set.seed(1) rsp <- c( sample(c(TRUE, FALSE), size = 40, prob = c(3 / 4, 1 / 4), replace = TRUE), @@ -759,7 +922,7 @@ test_that("s_proportion_diff errors when stratified method is chosen without str ) }) -test_that("s_proportion_diff errors when strata are provided with the non-stratified method `uncond_exact_diff`", { +test_that("s_proportion_diff errors when strata are provided with the non-stratified method uncond_exact_diff", { dta <- data.frame( rsp = sample(c("Y", "N"), 10, TRUE), grp = factor(rep(c("A", "B"), each = 5)), @@ -850,7 +1013,7 @@ testthat::test_that("s_proportion_diff ref column returns empty diff_est_ci", { testthat::expect_snapshot(res) }) -testthat::test_that("`estimate_proportion_diff` with diff_est_ci builds single-row table", { +testthat::test_that("estimate_proportion_diff with diff_est_ci builds single-row table", { set.seed(42, kind = "Mersenne-Twister") dta <- data.frame( rsp = sample(c(TRUE, FALSE), 100, TRUE), diff --git a/tests/testthat/test-test_proportion_diff.R b/tests/testthat/test-test_proportion_diff.R index 0b7eabe6ea..42bb924734 100644 --- a/tests/testthat/test-test_proportion_diff.R +++ b/tests/testthat/test-test_proportion_diff.R @@ -1,4 +1,4 @@ -testthat::test_that("prop_chisq returns right result", { +test_that("prop_chisq returns right result", { set.seed(1, kind = "Mersenne-Twister") rsp <- c( sample(c(TRUE, FALSE), size = 20, prob = c(3 / 4, 1 / 4), replace = TRUE), diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index b380dd9155..e2b4c1924c 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -167,9 +167,7 @@ testthat::test_that("n_available works as expected", { testthat::expect_snapshot(res) }) -################## -## range_noinf -################## +# range_noinf ---- # INTEGER no zero-len data, no NAs, no Inf @@ -777,3 +775,14 @@ testthat::test_that( testthat::expect_equal(mod, mod2) } ) + +testthat::test_that("uniroot_catch_na returns the root on success", { + res <- uniroot_catch_na(function(x) x^2 - 2, interval = c(0, 2)) + testthat::expect_equal(res, 1.414213, tolerance = 1e-6) +}) + +testthat::test_that("uniroot_catch_na returns NA when function evaluates to NA", { + f <- function(x) if (x > 0.5) NA_real_ else x - 0.25 + res <- uniroot_catch_na(f, interval = c(0, 1)) + testthat::expect_identical(res, NA_real_) +})